2.1. Shape of the domain and choice of the invariants

We choose the domain of the calculation to be a flux tube, a slender tube around a magnetic field line that is everywhere parallel to the magnetic field. The slenderness of the domain allows us to replace relevant functions by their linear expansions in the coordinates perpendicular to the magnetic field, which can be expressed as

(2.1) \begin{equation} \boldsymbol{B} = \boldsymbol{\nabla }\psi \times \boldsymbol{\nabla }\alpha , \end{equation}

where $\psi$ is the toroidal flux enclosed by a magnetic surface divided by $2 \pi$ , and $\alpha = \theta - \iota \varphi$ is the Clebsch angular coordinate labelling the different field lines on such surfaces with $\theta$ being a straight-field-line poloidal angle and $\iota$ being the rotational transform (related to the safety factor as $q = 1/\iota$ ). This representation suggests that $\psi$ and $\alpha$ be used to parametrise the directions perpendicular to the central magnetic field line $\boldsymbol{B}_0$ , and one can furthermore employ the arc length along this field line, $\ell$ , as a third coordinate. Away from the central field line of the flux tube, this coordinate is taken to be constant on planes perpendicular to $\boldsymbol{B}_0$ , i.e. $\boldsymbol{B}_0 \times \boldsymbol{\nabla }\ell =\boldsymbol{0}$ . All metric components associated with the coordinates are then defined and depend, to leading order in the width of the flux tube, only on the coordinate $\ell$ . Other functions may readily be expressed by their linear expansions in the perpendicular coordinates $\psi$ and $\alpha$ , where the expansion coefficients only depend on the parallel coordinate $\ell$ , and we have thus parametrised $\boldsymbol{r}=(\psi ,\alpha ,\ell )$ .

There are two main branches of the ITG instability: the slab mode and the curvature-driven one. In the following analysis, the available energy of the ions is calculated using an ordering that singles out the latter instability branch by taking the parallel ion dynamics to be slower than the perpendicular dynamics, i.e. $k_\| v_T \ll \boldsymbol{k}_\perp \boldsymbol{\cdot }\boldsymbol{v}_D$ . Here, $k_\|$ denotes the parallel wave number, $v_T$ the thermal ion speed, $\boldsymbol{k}_\perp$ the wave number perpendicular to the magnetic field and $\boldsymbol{v}_D$ the magnetic drift velocity (Biglari, Diamond & Rosenbluth Reference Biglari, Diamond and Rosenbluth1989; Plunk et al. Reference Plunk, Helander, Xanthopoulos and Connor2014). As such, we take the parallel velocity $v_\|$ of each ion to be invariant, as the perpendicular drifts dominate the dynamics. Furthermore, we assume that the ion magnetic moment $\mu = m v_\perp ^2/2B$ is invariant, where $m$ is the ion mass, $v_\perp$ is the perpendicular velocity and $B$ is the magnetic field strength. Finally, we take the parallel coordinate $\ell$ as being invariant as well, implying that parallel movement of the ions is much slower than typical turbulent time scales. In other words, Gardner re-stacking is only allowed on planes perpendicular to $\boldsymbol{B}_0$ ,Footnote ² since the parallel ion motion is unimportant. It is worthwhile to briefly consider when such an ordering is valid. Define the ‘transit frequency’ of a particle travelling at thermal speed $v_{T}$ to be $\omega _{\textrm {tran}} \sim v_{T}/\textit{qR}$ , i.e. how often it traverses a ‘connection length’ $\textit{qR}$ per unit time, where $R$ is the major radius of the toroidal device. Further taking the characteristic frequency of the turbulence to be $\omega _{\text{turb}} \sim v_T/L_T$ , where $L_T \sim T/|\boldsymbol{\nabla }T|$ , requiring $\omega _{\text{tran}} \ll \omega _{\text{turb}}$ necessitates

(2.2) \begin{equation} L_{T} \ll \textit{qR}. \end{equation}

In other words, if the gradients in temperature are sufficiently large, the parallel dynamics may safely be neglected.

We recall that the cross-section of the domain at constant $\ell$ cannot be chosen freely if it is assumed that both the initial distribution function and the ground state can be replaced by their linear expansions, where the full argument is given by Mackenbach et al. (Reference Mackenbach, Proll, Wakelkamp and Helander2023b ), of which we give a brief overview. Since the phase-space Jacobian is

(2.3) \begin{equation} \begin{aligned} \mathrm{d} \boldsymbol{r}\, \mathrm{d} \boldsymbol{v} = \frac {\pi }{m} \; \mathrm{d} \psi\, \mathrm{d} \alpha\, \mathrm{d}\ell\, \mathrm{d} \mu \,\mathrm{d} v_{\|}, \end{aligned} \end{equation}

the following quantity is constant in time $t$ as a consequence of the Casimir invariants (since $\mu$ , $v_\|$ and $\ell$ are invariant):

(2.4) \begin{equation} \mathcal{A} = \int _{\mathcal{D}} \mathrm{d} \psi\, \mathrm{d} \alpha \; \varTheta \left [ f(t,\psi ,\alpha ,\ell ,\mu ,v_{\|}) - \phi \right ], \end{equation}

where $\varTheta [x]$ is the Heaviside function, $f$ the plasma distribution function, $\phi$ any real scalar constant and $\mathcal{D}$ denotes the integration domain. If we furthermore assume that any distribution function may be written as $f \approx f_0 + f_\psi (\psi -\psi _0) + f_\alpha (\alpha -\alpha _0)$ (where subscripts to $f$ denote partial derivatives, i.e. $\partial _x f \equiv f_x$ , while $\psi _0$ and $\alpha _0$ denote corresponding values at the centre of the domain), the condition $\partial _t \mathcal{A}(t,\phi ,\mu ,v_{\|},\ell ) = 0$ cannot be satisfied in all domain shapes $\mathcal{D}$ , where a subscript zero refers to the fact that the function is evaluated at the centre of the domain. A geometric argument can be used to show that the only sufficiently smooth domain $\mathcal{D}$ that conserves $\mathcal{A}$ is elliptical in the $(\psi ,\alpha )$ -plane. In terms of the rescaled coordinates $x = (\psi - \psi _0) / \Delta \psi$ and $y = (\alpha - \alpha _0)/\Delta \alpha$ , the domain satisfies the equation

(2.5) \begin{equation} \mathcal{D} = \left \{ (x,y) \colon \left ( \frac {x'}{a} \right )^2 + \left ( \frac {y'}{b} \right )^2 \leqslant 1 \right \}, \end{equation}

where $x' = x \cos \vartheta - y \sin \vartheta$ and $y' = x \sin \vartheta + y \cos \vartheta$ for some rotation-angle $\vartheta$ (which is different from the poloidal angle $\theta$ ). If one violates this condition by choosing e.g. a square cross-section, still approximating the ground state with its linear expansion, unphysical results arise such as negative available energies under certain conditions. We highlight that $\Delta \psi$ , $\Delta \alpha$ and the angle $\vartheta$ can, in general, depend on the invariants.

This domain shape corresponds to the most general case and the available energy may be calculated in this tilted ellipse. However, for simplicity, we specialise to the case of a circular domain in $(x,y)$ -space, thus taking $a=b=1$ .Footnote ³ For any function $h(x,y)$ that can be approximated by a linear expansion, $h \approx h_x x + h_y y$ on the domain $\mathcal{D}$ , we note the following integral required for (1.4):

(2.6) \begin{equation} \begin{aligned} \int _{\mathcal{D}} \mathrm{d} x \,\mathrm{d} y \; \delta \left [ h(x,y) \right ] = \frac {2}{\sqrt {h_x^2 + h_y^2}}. \end{aligned} \end{equation}

2.2. Available energy

With the domain shape $\mathcal{D}$ given in (2.5) and setting $a=b=1$ , we now turn to the calculation of the available energy associated with an initially Maxwellian distribution function in local thermodynamic equilibrium (where we assume the temperature and density are flux functions, i.e. $n = n(\psi )$ and $T=T(\psi )$ ),

(2.7) \begin{equation} f_M(\psi ) = n(\psi ) \left [ \frac {m}{2 \pi T(\psi )} \right ]^{3/2} e^{-(\mu B + mv_{\|}^2/2)/T(\psi )}. \end{equation}

It is straightforward to evaluate the perpendicular derivatives of this plasma distribution function. For economy of notation, we write

(2.8a) \begin{align} \omega _*^T(\epsilon ) &\equiv T_0 \left [ \frac {\mathrm{d} \ln n}{\mathrm{d} \psi } + \frac {\mathrm{d} \ln T}{\mathrm{d} \psi } \left ( \frac {\epsilon }{T} - \frac {3}{2} \right )\right ], \end{align}

(2.8b) \begin{align} \omega _\alpha (\mu ,\ell ) &\equiv \mu B_0 \left (\frac {\partial \ln B}{\partial \psi } \right )_{\alpha ,\ell }, \end{align}

(2.8c) \begin{align} \omega _\psi (\mu ,\ell ) &\equiv -\mu B_0 \left ( \frac {\partial \ln B}{\partial \alpha } \right )_{\psi ,\ell }, \end{align}

where the particle energy is $\epsilon = \mu B + m v_{\|}^2/2$ , quantities with a subscript zero are evaluated at the centre of the elliptical domain $(\psi _0,\alpha _0)$ and understand that the derivatives are evaluated at the centre of the flux tube too. We note that the perpendicular derivatives of the magnetic field are related to the particle drifts at the centre of the flux tube, as may be seen by expressing the quantity $\boldsymbol{v}_{\boldsymbol{\nabla }}= \boldsymbol{b} \times \boldsymbol{\nabla }\ln B$ in $(\psi ,\alpha ,\ell )$ coordinates,

(2.9) \begin{align} \boldsymbol{v}_{\boldsymbol{\nabla }}= \frac {1}{B^2} \left ( \frac {\partial B}{\partial \psi } \boldsymbol{B} \times \boldsymbol{\nabla }\psi + \frac {\partial B}{\partial \alpha } \boldsymbol{B} \times \boldsymbol{\nabla }\alpha \right ), \end{align}

where we have employed $\boldsymbol{B} \times \boldsymbol{\nabla }\ell = \boldsymbol{0}$ . Taking inner products with $\boldsymbol{\nabla }\psi$ and $\boldsymbol{\nabla }\alpha$ , one finds the relations

(2.10) \begin{align} \left (\frac {\partial B}{\partial \psi }\right )_{\alpha ,\ell } = \boldsymbol{v}_{\boldsymbol{\nabla }}\boldsymbol{\cdot }\boldsymbol{\nabla }\alpha , \quad \left (\frac {\partial B}{\partial \alpha }\right )_{\psi ,\ell } = - \boldsymbol{v}_{\boldsymbol{\nabla }}\boldsymbol{\cdot }\boldsymbol{\nabla }\psi . \end{align}

Now, to solve the ground-state equation, we require perpendicular derivatives of both the particle energy and the initial distribution function, and definitions of (2.10) allow us to succinctly express them as

(2.11a) \begin{align} \left (\frac {\partial [\epsilon + \kappa ]}{\partial \psi } \right )_{\mu ,v_{\|},\alpha ,\ell } &= \omega _\alpha + \kappa _\alpha ,\end{align}

(2.11b) \begin{align} \left (\frac {\partial [\epsilon +\kappa ]}{\partial \alpha } \right )_{\mu ,v_{\|},\psi ,\ell } &= - \omega _\psi + \kappa _\psi , \end{align}

(2.11c) \begin{align} \left (\frac {\partial f_M}{\partial \psi } \right )_{\mu ,v_{\|},\alpha ,\ell } &= \frac {f_{M,0}}{T_0} \left ( \omega _*^T - \omega _\alpha \right ), \end{align}

(2.11d) \begin{align} \left (\frac {\partial f_M}{\partial \alpha } \right )_{\mu ,v_{\|},\psi ,\ell } &= \frac {f_{M,0}}{T_0} \omega _\psi , \end{align}

where $\partial _\psi \kappa |_{\alpha ,\ell } = \kappa _\alpha$ and $\partial _\alpha \kappa |_{\psi ,\ell } = \kappa _\psi$ are, for the moment, unknown. Without loss of generality, we set $\psi _0 = \alpha _0 = 0$ henceforth, for ease of notation. To calculate the numerator in the ground state (1.4) (derived by equating the Casimirs), we need the following integral as a function of the particle energy $w$ :

(2.12) \begin{equation} \begin{aligned} \frac { \pi }{m} \int _{\mathcal{D}} \mathrm{d}\psi \,\mathrm{d} \alpha \; \delta \left [ w - \epsilon - \kappa \right ] = \frac {\pi \Delta \psi \Delta \alpha }{m} \frac {2}{ \sqrt { (\Delta \psi )^2 (\omega _\alpha + \kappa _\alpha )^2 +(\Delta \alpha )^2 (\omega _\psi - \kappa _\psi )^2} }, \end{aligned} \end{equation}

where we have employed (2.6), and $w - \epsilon (\psi _0,\alpha _0) - \kappa (\psi _0,\alpha _0)$ vanishes since the ground state and initial distribution function are identical to leading order in slenderness of the flux tube. We also require the denominator, i.e. the ‘waterbag distribution’ (Ewart, Nastac & Schekochihin Reference Ewart, Nastac and Schekochihin2023),

(2.13) \begin{equation} \begin{aligned} \frac {\pi }{m} \int _{\mathcal{D}} \mathrm{d}\psi \,\mathrm{d}\alpha \; \delta \left [ f_M - F \right ] = \frac {\pi \Delta \psi \Delta \alpha T_0}{m f_{M,0}} \frac {2}{\sqrt { (\Delta \psi )^2 (\omega _*^T - \omega _\alpha )^2 +(\Delta \alpha )^2 (\omega _\psi )^2}}, \end{aligned} \end{equation}

where, to leading-order, $f_{M,0}=F$ . Substituting these results in (1.4), we find that the ground state becomes

(2.14) \begin{equation} F'[\epsilon + \kappa ;\mu ,v_\|,\ell ] = - \frac {f_{M,0}}{T_0} G, \end{equation}

where a prime denotes differentiation with respect to the first argument, and $G$ is defined by

(2.15) \begin{equation} G[\kappa _\psi (\ell ),\kappa _\alpha (\ell )] \equiv \sqrt {\frac {(\Delta \psi )^2 (\omega _\alpha - \omega _*^T )^2 +(\Delta \alpha )^2 (\omega _\psi )^2}{ (\Delta \psi )^2 (\omega _\alpha + \kappa _\alpha )^2 +(\Delta \alpha )^2 (\omega _\psi - \kappa _\psi )^2}}. \end{equation}

The functions $\kappa _\psi$ and $\kappa _\alpha$ are determined by the condition that the density remains fixed,

(2.16) \begin{equation} \int \mathrm{d} \mu \,\mathrm{d} v_\| \; (F[\epsilon +\kappa ;\mu ,v_\|,\ell ]-f_M) = 0, \end{equation}

and we wish to manipulate this equation to find explicit expressions for these functions. Taking partial derivatives of (2.16) with respect to $\psi$ and $\alpha$ , and invoking (2.11) and (2.14), gives two coupled nonlinear algebraic equations for $\kappa _\alpha$ and $\kappa _\psi$ ,

(2.17a) \begin{align} \int \mathrm{d} \mu \,\mathrm{d} v_\| \; \frac {f_{M,0}}{T_0} \left ( -\kappa _\alpha -\omega _\alpha \right )G(\kappa _\psi ,\kappa _\alpha ) &= \int \mathrm{d} \mu \,\mathrm{d} v_\| \; \frac {f_{M,0}}{T_0}\left ( \omega _*^T - \omega _\alpha \right ) ,\end{align}

(2.17b) \begin{align} \int \mathrm{d} \mu\, \mathrm{d} v_\| \; \frac {f_{M,0}}{T_0} \left (\omega _\psi - \kappa _\psi \right ) G(\kappa _\psi ,\kappa _\alpha ) &= \int \mathrm{d} \mu\, \mathrm{d} v_\| \;\frac {f_{M,0}}{T_0} \omega _\psi ,\end{align}

and we prove uniqueness of the solution pair $(\kappa _\psi ,\kappa _\alpha )$ in Appendix A.

With the ground-state given in (2.14), we can now evaluate the available energy,

(2.18) \begin{equation} A = \frac {\pi }{m} \int _{\mathcal{D} \times \mathcal{V}} \mathrm{d} \psi\, \mathrm{d} \alpha \,\mathrm{d} \ell\, \mathrm{d} \mu \,\mathrm{d} v_{\|} \; (f_M - F) \epsilon , \end{equation}

where $\mathcal{V}: \: (\mu \geqslant 0,\, v_\| \in \mathbb{R},\, \ell \in \mathbb{R})$ . As in the case of the TÆ, we note that $\int (f_M-F)$ $\mathrm{d}\psi \,\mathrm{d}\alpha = 0$ , because the number of particles with given values of $(\mu ,v_{\|},\ell )$ is conserved. The integral can thus be evaluated to leading order as

(2.19) \begin{equation} \begin{aligned} A =& \frac {\pi }{m}\int _{\mathcal{D} \times \mathcal{V}} \mathrm{d} \psi\, \mathrm{d} \alpha \,\mathrm{d} \ell \,\mathrm{d} \mu\, \mathrm{d} v_{\|} \; \left [ \frac {\partial (\epsilon +\kappa )}{\partial \psi } \frac {\partial \left ( f_M - F \right )}{\partial \psi } \psi ^2 + \frac {\partial (\epsilon +\kappa )}{\partial \alpha } \frac {\partial \left ( f_M - F \right )}{\partial \alpha } \alpha ^2 \right ] \\ &- \frac {\pi }{m}\int _{\mathcal{D} \times \mathcal{V}} \mathrm{d} \psi\, \mathrm{d} \alpha\, \mathrm{d} \ell\, \mathrm{d} \mu\, \mathrm{d} v_{\|} \; \left [ \frac {\partial \kappa }{\partial \psi } \frac {\partial \left ( f_M - F \right )}{\partial \psi } \psi ^2 + \frac {\partial \kappa }{\partial \alpha } \frac {\partial \left ( f_M - F \right )}{\partial \alpha } \alpha ^2 \right ], \end{aligned} \end{equation}

having simply added and subtracted the same term, and understanding that the partial derivatives are performed in $(\psi ,\alpha ,\ell ,\mu ,v_\|)$ -coordinates with remaining coordinates held fixed. The second term on the right-hand side of the equation vanishes as a consequence of (2.16), and writing out the derivatives in full, we find

(2.20) \begin{equation} \begin{aligned} A =-&\frac {\pi }{m} \int _{\mathcal{D} \times \mathcal{V}} \mathrm{d} \psi \,\mathrm{d} \alpha\, \mathrm{d} \ell \,\mathrm{d} \mu \,\mathrm{d} v_{\|} \; \frac {f_{M,0}}{T_0} \left [(\omega _\alpha + \kappa _\alpha ) (\omega _\alpha - \omega _*^T) \psi ^2 + (\omega _\psi - \kappa _\psi ) \omega _\psi \alpha ^2 \right ] \\ & +\frac {\pi }{m} \int _{\mathcal{D} \times \mathcal{V}} \mathrm{d} \psi\, \mathrm{d} \alpha \,\mathrm{d} \ell\, \mathrm{d} \mu\, \mathrm{d} v_{\|} \; \frac {f_{M,0}}{T_0} G \left [ \left (\omega _\alpha + \kappa _\alpha \right )^2 \psi ^2 + (\omega _\psi - \kappa _\psi )^2 \alpha ^2 \right ] . \end{aligned} \end{equation}

The integrals over the elliptical cross-section in the $(\psi ,\alpha )$ -plane can readily be evaluated, $\int _{\mathcal{D}} \psi ^2 \: \mathrm{d} \psi\, \mathrm{d} \alpha = \pi (\Delta \psi )^3 (\Delta \alpha )/4$ and $\int _{\mathcal{D}} \alpha ^2 \: \mathrm{d} \psi\, \mathrm{d} \alpha = \pi (\Delta \psi ) (\Delta \alpha )^3/4$ , so that the available energy becomes (using (2.14))

(2.21) \begin{equation} \begin{aligned} A = & \frac {\pi ^2 \Delta \psi \Delta \alpha }{4m} \int \mathrm{d} \ell\, \mathrm{d} \mu\, \mathrm{d} v_{\|} \;\frac {f_{M,0}}{T_0} \\ &\times\Bigg [ \sqrt {(\Delta \psi )^2 (\omega _\alpha + \kappa _\alpha )^2 +(\Delta \alpha )^2 (\omega _\psi - \kappa _\psi )^2} \sqrt {(\Delta \psi )^2 (\omega _\alpha - \omega _*^T )^2 +(\Delta \alpha )^2 (\omega _\psi )^2} \\ &-(\Delta \psi )^2(\omega _\alpha + \kappa _\alpha ) (\omega _\alpha - \omega _*^T) - (\Delta \alpha )^2 (\omega _\psi - \kappa _\psi ) \omega _\psi \Bigg ],\\[-18pt] \end{aligned} \end{equation}

dropping the domain of integration $\mathcal{V}$ from here on. As expected, the expression for the available energy is always positive definite, as may be verified by the following argument. Note that the function in square brackets in (2.21) is of the form

(2.22) \begin{equation} \sqrt {a^2 + b^2}\sqrt {c^2 + d^2} - ac - bd, \end{equation}

where $a = \Delta \psi (\omega _\alpha + \kappa _\alpha )$ , $b= \Delta \alpha (\omega _\psi - \kappa _\psi )$ , $c =\Delta \psi (\omega _\alpha - \omega _*^T)$ and $d= \Delta \alpha (\omega _\psi )$ . To prove positive definiteness, we must show $\sqrt {(a^2 + b^2)(c^2 + d^2)} \geqslant ac + bd$ , which follows readily as squaring both sides results in $(ad)^2 + (bc)^2 \geqslant 2 (ad) (bc)$ , simplifying to $(ad - bc)^2 \geqslant 0$ , as required. We have now fully specified the available energy: solving the nonlinear algebraic (2.17) for each $\ell$ gives the functions $\kappa _\alpha (\ell )$ and $\kappa _\psi (\ell )$ , which are required to calculate the available energy given in (2.21).

2.3. Conditions for stability

It is worthwhile to consider under what conditions the available energy vanishes (apart from the trivial case of vanishing gradients of density and temperature, or $\boldsymbol{v}_{\boldsymbol{\nabla }}\boldsymbol{\cdot }\boldsymbol{\nabla }\psi = \boldsymbol{v}_{\boldsymbol{\nabla }}\boldsymbol{\cdot }\boldsymbol{\nabla }\alpha = 0$ ), and we shall refer to this as stability. In the notation of (2.22), this occurs if and only if $ad = bc$ , and hence,

(2.23) \begin{equation} \kappa _\psi (\ell )\left [ \omega _*^T(\epsilon ) - \omega _\alpha (\mu ,\ell ) \right ] = \omega _\psi (\mu ,\ell ) \left [\omega _*^T(\epsilon ) + \kappa _\alpha (\ell )\right ]. \end{equation}

Taking the derivative of (2.23) with respect to $v_\|$ at fixed $(\mu ,\ell )$ gives

(2.24) \begin{equation} [\kappa _\psi (\ell )-\omega _\psi (\mu ,\ell )] \partial _\epsilon \omega _{*}^{T}(\epsilon )=0 \end{equation}

(having divided out $\partial \epsilon /\partial v_\|$ ), which can only be satisfied under two conditions.

1. (i) If the temperature gradient vanishes, $ \partial _\epsilon \omega _*^T = \mathrm{d} \ln T / \mathrm{d} \psi = 0$ according to (2.8). Note that there may still be a non-zero density gradient.

2. (ii) If the magnetic field is isodynamic (see Helander Reference Helander2014) so that $\omega _{\psi } = 0$ everywhere and consequently $\kappa _\psi = 0$ , consistently with (2.17b ).

Let us now investigate the first condition.

2.3.1. Pure density gradient

Under the assumption of no temperature gradient, the derivative of (2.23) with respect to $(\mu ,\ell )$ simplifies to

(2.25) \begin{equation} \kappa _\psi \omega _\alpha + (\kappa _\alpha + \omega _*^T) \omega _\psi = 0, \end{equation}

giving a direct relation between $\kappa _\psi$ and $\kappa _\alpha$ . Substituting this back into (2.23), we find

(2.26) \begin{equation} \frac {\omega _*^T \omega _\psi }{\omega _\alpha }(\kappa _\alpha + \omega _*^T) = 0. \end{equation}

Assuming the drifts do not vanish, we thus have

(2.27) \begin{align} \kappa _\alpha = - \omega _*^T, \quad \; \kappa _\psi = 0. \end{align}

It remains to be checked if the derived conditions are compatible with the constraint on the density. Substituting (2.28) into (2.15), we find that $G=1$ , and (2.17) is indeed satisfied. Hence, if there is no ITG present, the available energy always vanishes, as may be verified by filling in the results of (2.27) into (2.21). This is in line with Helander (Reference Helander2020), who found that the available energy always vanishes in a Maxwellian plasma with no temperature gradient and a fixed, slightly varying density profile (with no invariants aside from the Casimirs). Further constraining the ground state by means of additional invariants can only lower the available energy, which then must remain zero. As we shall see in the case of an isodynamic field, simply satisfying (2.23) is not sufficient for having zero available energy, and more care is required.

2.3.2. Isodynamic magnetic field

In the case where $\omega _\psi$ and $\kappa _\psi$ both vanish, the available energy given in (2.21) becomes

(2.28) \begin{equation} \begin{aligned} A = & \frac {\pi ^2 \Delta \psi \Delta \alpha }{2m} \int \mathrm{d} \ell\, \mathrm{d} \mu\, \mathrm{d} v_{\|} \; \frac {f_{M,0}}{T_0} (\Delta \psi )^2 \mathcal{R} \left [(\kappa _\alpha +\omega _\alpha ) (\omega _*^T - \omega _\alpha ) \right ], \end{aligned} \end{equation}

where the ramp function is defined as $\mathcal{R}(x) = (x + |x|)/2$ . For stability, we thus require that $\kappa _\alpha + \omega _\alpha$ and $\omega _*^T-\omega _\alpha$ have opposite signs for all $(\ell ,\mu ,v_\|)$ , where it is important to remember that $\kappa _\psi$ and $\kappa _\alpha$ only depend on $\ell$ while $\omega _\alpha$ is independent of $v_\|$ . We note that

(2.29) \begin{equation} \omega _*^T - \omega _\alpha = T_0\frac {\mathrm{d} \ln n}{\mathrm{d} \psi } \left [ 1 -\frac {3}{2}\eta + \frac {m v_\|^2}{2 T_0} \eta + \frac {\mu B_0}{T_0} \left ( \eta - \eta _B \right ) \right ], \end{equation}

where we have defined $\eta = \mathrm{d} \ln T/\mathrm{d} \ln n$ and $\eta _B= (\partial \ln B / \partial \psi )_{\alpha ,\ell }/(\mathrm{d} \ln n / \mathrm{d} \psi )$ . If $\eta \gt 2/3$ , the constant term in the square brackets of (2.29) is negative and we must require $mv_\|^2\eta /2T_0\lt 0$ to keep the sign fixed (at fixed $\mu$ and $\ell$ ), resulting in incompatible conditions. In contrast, if however $\eta \leqslant 2/3$ , we need $mv_\|^2\eta /2T_0 \geqslant 0$ , giving one condition for stability:

(2.30) \begin{equation} 0 \leqslant \eta \leqslant \frac {2}{3}. \end{equation}

This is sufficient to maintain the same sign of (2.29) for $\mu = 0$ and any $v_\|$ , and we now seek conditions that ensure the sign stays fixed for all $\mu$ . Suppose momentarily that $\omega _*^T - \omega _\alpha$ changes sign for some $\mu$ (so that $\mu B_0 (\eta - \eta _B)/T_0\lt 0$ ). To keep the product $(\kappa _\alpha +\omega _\alpha )(\omega _\alpha - \omega _*^T)$ of the same sign, $\kappa _\alpha + \omega _\alpha$ must change sign too, implying that both factors must share the same zero-line in phase space. The zero-line of the first factor is given by $\mu B_0/T_0 = - \kappa _\alpha / T_0 \partial _\psi \ln B$ , whereas the second factor gives

(2.31) \begin{equation} \frac {\mu B_0}{T_0} = \frac {1 - ({3}/{2}) \eta + ({mv_\|^2}/{2T_0}) \eta }{\eta _B - \eta }. \end{equation}

Equating the two, we must have that

(2.32) \begin{equation} \frac {\kappa _\alpha }{T_0} =\frac {1 - ({3}/{2})\eta + ({m v_\|^2}/{2T_0}) \eta }{\eta - \eta _B} \partial _\psi \ln B, \end{equation}

and we arrive at a contradiction since $\kappa _\alpha$ is a function of the arc-length alone. Thus, if $\mu B_0 (\eta - \eta _B)/T_0 \lt 0$ , there is non-zero available energy present.

Assuming now instead $\mu B(\eta - \eta _B)/T_0 \geqslant 0$ , so that maintaining the sign of (2.29) requires

(2.33) \begin{equation} \eta \geqslant \eta _B. \end{equation}

Since $\eta _B = \eta _B(\ell )$ , ensuring the highest value of $\eta _B(\ell )$ satisfies this equality is necessary for stability,

(2.34) \begin{equation} \eta \geqslant \max _{\ell } [\eta _B(\ell )], \end{equation}

where $\max [h(x)]$ returns the highest value the function $h(x)$ takes over the domain. Let us investigate what this inequality implies in typical conditions, relating it to the particle drifts via (2.10). For standard density and temperature profiles that decrease as functions of minor radius, i.e. $\partial _\psi \ln T \lt 0$ and $\partial _\psi \ln n \lt 0$ , the inequality simplifies to

(2.35) \begin{equation} \min \left [ \frac {\boldsymbol{v}_{\boldsymbol{\nabla }}\boldsymbol{\cdot }\boldsymbol{\nabla }\alpha }{B_0} \right ] \geqslant -|\partial _\psi \ln T|. \end{equation}

Equation (2.35) has a surprising corollary: if it is not satisfied, one can increase the magnitude of the ITG to satisfy it, and, if one is stable under the conditions derived, this would imply decreasing the available energy. This condition for stability is furthermore automatically met for a magnetic field with good curvature only, which we define as $\partial _\psi B |_{\alpha ,\ell } = \boldsymbol{v}_{\boldsymbol{\nabla }}\boldsymbol{\cdot }\boldsymbol{\nabla }\alpha \gt 0$ everywhere.

We have now derived conditions that keep the sign of (2.29) fixed and we must now verify whether the sign of $\kappa _\alpha + \omega _\alpha$ is always opposite to that of $\omega _*^T - \omega _\alpha$ , so that the ramp function in (2.28) is always greater than or equal to zero. To this end, we investigate the equation for $\kappa _\alpha$ , (2.17a ), which in the isodynamic limit reduces to

(2.36) \begin{equation} \int \mathrm{d} \mu\, \mathrm{d} v_\| \; \frac {f_{M,0}}{T_0} \sigma (\kappa _\alpha + \omega _\alpha ) |\omega _*^T - \omega _\alpha | = - \int \mathrm{d} \mu\, \mathrm{d} v_\| \; \frac {f_{M,0}}{T_0}\left ( \omega _*^T - \omega _\alpha \right ) , \end{equation}

where $\sigma (x) = x/|x|$ returns the sign of $x$ . Under the necessary conditions for stability derived so far, we may set

(2.37) \begin{equation} \omega _*^T - \omega _\alpha = \sigma (\omega _*^T - \omega _\alpha ) |\omega _*^T - \omega _\alpha |, \end{equation}

where $\sigma (\omega _*^T - \omega _\alpha )$ is a constant. Consequently, (2.36) becomes

(2.38) \begin{equation} \int \mathrm{d} \mu \,\mathrm{d} v_\| \; \frac {f_{M,0}}{T_0} \sigma (\kappa _\alpha + \omega _\alpha ) |\omega _*^T - \omega _\alpha | = - \sigma (\omega _*^T - \omega _\alpha ) \int \mathrm{d} \mu \,\mathrm{d} v_\| \; \frac {f_{M,0}}{T_0}|\omega _*^T - \omega _\alpha |, \end{equation}

and we have that

(2.39) \begin{equation} \sigma (\kappa _\alpha + \omega _\alpha ) = - \sigma (\omega _*^T - \omega _\alpha ), \end{equation}

as required for stability.

Figure 1.

Stability diagram in an isodynamic field where we have assumed $\partial _\psi \ln B \lt 0$ (typically the case locally on the outboard mid-plane of an up-down symmetric tokamak, i.e. the bad-curvature side). The dashed line denotes a pure density gradient, which is always stable. The region between the dashed, dotted and full line satisfies $0 \leqslant \eta \leqslant 2/3$ , but not $\eta \geqslant \eta _B$ , and thus has non-zero available energy.

Summarising, given an isodynamic magnetic field, we have found that necessary and sufficient conditions for stability are $0 \leqslant \eta \leqslant 2/3$ , and furthermore $\eta \geqslant \eta _B$ .Footnote ⁴ A clarifying plot may be found in figure 1, where the stability boundaries are plotted as a function of the gradient for $\partial _\psi \ln B \gt 0$ . We note close correspondence to results found by Costello & Plunk (Reference Costello and Plunk2025b ) (see also Plunk Reference Plunk2015), where one important difference pertains to the curvature drift. This is taken into account in the calculation by Costello & Plunk (Reference Costello and Plunk2025b ), but does not enter in the available energy. This is a simple consequence of the fact that all variation in $\epsilon = \mu B + mv_\|^2/2$ at fixed $\mu$ and $v_\|$ is due to $B$ , manifesting as gradient drifts.

2.3.3. Effect of the constraint on density

It is helpful to investigate what effect the assumed constancy of the density gradient has on the stability criteria found. This can be easily investigated, as setting $\kappa _\alpha = \kappa _\psi = 0$ is equivalent to the unconstrained case, where the density profile is such that it minimises the thermal energy of the initial Maxwellian distribution function. Equation (2.23) becomes

(2.40) \begin{equation} \omega _\psi (\mu ,\ell ) \omega _*^T(\epsilon ) = 0, \end{equation}

which can only be be satisfied in isodynamic fields, disregarding the trivial case of having no gradients in density and temperature. In isodynamic fields, the available energy can be found by setting $\kappa _\alpha = 0$ in (2.28), resulting in

(2.41) \begin{equation} \begin{aligned} A = & \frac {\pi ^2 \Delta \psi \Delta \alpha }{2m} \int \mathrm{d} \ell \,\mathrm{d} \mu\, \mathrm{d} v_{\|} \; \frac {f_{M,0}}{T_0} (\Delta \psi )^2 \omega _\alpha ^2 \mathcal{R} \left (\frac {\omega _*^T}{\omega _\alpha } - 1 \right ) . \end{aligned} \end{equation}

To have zero available energy, the argument of the ramp function must be negative everywhere, and we thus require

(2.42) \begin{equation} \frac {\omega _*^T}{\omega _\alpha } - 1 = \frac {T_0}{\eta _B \mu B_0} \left [ 1 -\frac {3}{2}\eta + \frac {m v_\|^2}{2 T_0} \eta + \frac {\mu B_0}{T_0} \left ( \eta - \eta _B \right ) \right ] \leqslant 0. \end{equation}

As already derived, for the factor in brackets to have the same sign (namely positive) everywhere, we require $0 \leqslant \eta \leqslant 2/3$ and $\eta _B \leqslant \eta$ . However, to ensure the argument is negative, we require $\eta _B \leqslant 0$ , which is a more stringent condition. In terms of a device with radially decreasing density profile, for stability, we thus require an isodynamic field with

(2.43) \begin{equation} \partial _\psi \ln B \geqslant 0 \end{equation}

everywhere to have vanishing available energy; hence, the device has good curvature everywhere. Conversely, the bad-curvature region of an isodynamic magnetic field will always have energy available for non-zero gradients, if the density is unconstrained. Physically, this may be understood by considering a particle that is moved radially to flatten the density gradient. When the density gradient is parallel to the magnetic field strength gradient, flattening the density gradient is accomplished by moving particles from regions of high to low field strength. The kinetic energy released by moving a particle from high to low field strength at constant $\mu$ and $v_\|$ ,

(2.44) \begin{align} \Delta \epsilon = \mu \left ( B_{\text{high}} - B_{\text{low}} \right ), \end{align}

is positive, implying that this process feeds the instability. Such a process is not allowed if the density gradient is kept fixed, resulting in broader regions of stability.

2.4. Limit of steep density and temperature gradients

Having derived conditions for zero available energy, we now focus on the limiting case where the gradients in density and temperature are very large. We thus set $\omega _*^T \rightarrow \omega _*^T/\varepsilon$ , and expand the relevant equations under the assumption that the positive expansion parameter $\varepsilon \ll 1$ . However, it is not clear a priori what scaling $\kappa _\alpha$ and $\kappa _\psi$ should follow. As such, let us expand these scalars generally as

(2.45) \begin{align} \kappa _{\alpha } \approx \frac {\kappa _{\alpha ,-1}}{\varepsilon } + \kappa _{\alpha ,0} + \ldots, \qquad \kappa _{\psi } \approx \frac {\kappa _{\psi ,-1}}{\varepsilon } + \kappa _{\psi ,0} + \ldots , \end{align}

and we now solve equations with different assumptions on the expansions of $\kappa _\psi$ and $\kappa _\alpha$ .

2.4.1. Weak $\kappa$ -scaling

Let us first assume a weak scaling of both $\kappa _\psi$ and $\kappa _\alpha$ with $\varepsilon$ , where we set $\kappa _{\alpha ,-1}=\kappa _{\psi ,-1}=0$ so that (2.17) is, to leading order, given as

(2.46a) \begin{align} \int \mathrm{d} \mu\, \mathrm{d} v_\| \; \frac {f_{M,0}}{T_0 \varepsilon } \left [ \omega _*^T + \frac {|\omega _*^T| (\Delta \psi ) (\kappa _{\alpha ,0} + \omega _{\alpha }) }{\sqrt {(\Delta \psi )^2 (\kappa _{\alpha ,0} + \omega _\alpha )^2 + (\Delta \alpha )^2 (\omega _\psi - \kappa _{\psi ,0})^2}} \right ] &\approx 0,\end{align}

(2.46b) \begin{align} \int \mathrm{d} \mu\, \mathrm{d} v_\| \; \frac {f_{M,0}}{T_0 \varepsilon } \frac {|\omega _*^T| (\Delta \alpha ) \left (\omega _\psi - \kappa _{\psi ,0} \right )}{\sqrt {(\Delta \psi )^2 (\kappa _{\alpha ,0} + \omega _\alpha )^2 + (\Delta \alpha )^2 (\omega _\psi - \kappa _{\psi ,0})^2}} &\approx 0.\end{align}

Similarly expanding the available energy from (2.21) gives

(2.47) \begin{equation} \begin{aligned} A \approx \frac {\pi ^2 \Delta \psi \Delta \alpha }{4m \varepsilon } \int \mathrm{d} \ell\, \mathrm{d} \mu\, \mathrm{d} v_{\|} \; & \frac {f_{M,0}}{T_0} \Big [ (\Delta \psi )^2 \omega _*^T\left ( \omega _\alpha + \kappa _{\alpha ,0} \right ) \\ & +\Delta \psi |\omega _*^T| \sqrt {(\Delta \psi )^2 (\omega _\alpha + \kappa _{\alpha ,0})^2 + (\Delta \alpha )^2 (\omega _\psi - \kappa _{\psi ,0})^2} \Big ] .\\[-12pt]& \end{aligned} \end{equation}

Equations (2.46) and (2.47) can be solved in general numerically. However, we note that the current scaling is unable to describe the solution if $0 \leqslant \eta \leqslant 2/3$ , and a stronger scaling of $\kappa$ with $\varepsilon$ gives the correct solution, which may be seen as follows. When $0 \leqslant \eta \leqslant 2/3$ , the quantity $\omega _*^T$ is of the same sign everywhere and one may set $\omega _*^T = \sigma (\omega _*^T) |\omega _*^T|$ with $\sigma (\omega _*^T) = \pm 1$ being constant, and (2.46a ) is given by

(2.48) \begin{align} & \int \mathrm{d} \mu\, \mathrm{d} v_\| \; \frac {f_{M,0}}{T_0} \frac {|\omega _*^T| \Delta \psi (\kappa _{\alpha ,0}+ \omega _\alpha )}{\sqrt {(\Delta \psi )^2 (\kappa _{\alpha ,0}+ \omega _\alpha )^2 + (\Delta \alpha )^2 (\omega _\psi - \kappa _{\psi ,0})^2}}\nonumber \\ &\quad = -\sigma (\omega _*^T) \int \mathrm{d} \mu \,\mathrm{d} v_\|\; \frac {f_{M,0}}{T_0}|\omega _*^T| \end{align}

We see that the fraction on the left-hand side is bounded in magnitude from above by

(2.49) \begin{align} \frac {|\Delta \psi (\kappa _{\alpha ,0}+ \omega _\alpha )|}{\sqrt {(\Delta \psi )^2 (\kappa _{\alpha ,0}+ \omega _\alpha )^2 + (\Delta \alpha )^2 (\omega _\psi - \kappa _{\psi ,0})^2}} \lt 1. \end{align}

To solve (2.48), the fraction must approach $\pm 1$ , for which we require $\kappa _{\alpha ,0} \rightarrow \pm \infty$ . This invalidates the weak ordering presented, and one requires a stronger scaling of $\kappa _{\alpha }$ with $\varepsilon$ .

2.4.2. Strong $\kappa$ -scaling

Let us therefore take $\kappa _{\psi ,-1} \neq 0$ and $\kappa _{\alpha ,-1} \neq 0$ , so that both variables have leading-order dependence that scales as $1/\varepsilon$ . Expanding (2.17b ) gives rise to

(2.50) \begin{equation} \frac {\kappa _{\psi ,-1}}{\varepsilon } \int \mathrm{d} \mu\, \mathrm{d} v_{\|} \; \frac {f_{M,0}}{T_0} \frac {\Delta \psi |\omega _*^T|}{\sqrt {(\Delta \psi )^2\kappa _{\alpha ,-1}^2 + (\Delta \alpha )^2\kappa _{\psi ,-1}^2}} \approx 0, \end{equation}

and we see that the integrand is always positive definite. As such, the only way to satisfy (2.50) is by setting $\kappa _{\psi ,-1} = 0$ , assuming that $\kappa _{\alpha ,-1} \neq 0$ . Turning our attention to (2.17a ), its leading-order expansion is given as

(2.51) \begin{equation} \frac {1}{\varepsilon } \int \mathrm{d} \mu\, \mathrm{d} v_{\|} \; \frac {f_{M,0}}{T_0} \left [ \sigma (\kappa _{\alpha ,-1}) |\omega _*^T| + \omega _*^T \right ] \approx 0 \implies \sigma (\kappa _{\alpha ,-1}) \approx - \frac {\int \mathrm{d} \mu \,\mathrm{d} v_{\|} \; f_{M,0} \omega _*^T }{\int \mathrm{d} \mu \,\mathrm{d} v_{\|} \; f_{M,0} |\omega _*^T| }, \end{equation}

which can only be satisfied if $\omega _*^T$ is of the same sign everywhere, necessitating that the ratio of gradients obeys $0\leqslant \eta \leqslant 2/3$ , in line with what was found in the weak $\kappa$ -scaling.

Expanding the available energy accordingly gives

(2.52) \begin{equation} \begin{aligned} A = \frac {\pi ^2 \Delta \psi \Delta \alpha }{2 m \varepsilon ^2} & \int \mathrm{d} \ell\, \mathrm{d} \mu\, \mathrm{d} v_{\|} \; \frac {f_{M,0}}{T_0} (\Delta \psi )^2 \Bigg \{ \mathcal{R} \left ( \kappa _{\alpha ,-1}\omega _*^T \right ) + \\ &\frac { \left [ \omega _*^T (\kappa _{\alpha ,0} + \omega _\alpha ) - \kappa _{\alpha ,-1} \omega _\alpha \right ] \left [ 1 + \sigma (\kappa _{\alpha ,-1} \omega _*^T) \right ]}{2 } \varepsilon + \mathcal{O}(\varepsilon ^2) \Bigg \} , \end{aligned} \end{equation}

where we note that, since $\omega _*^T$ is of the same sign everywhere when $0\leqslant \eta \leqslant 2/3$ , $\kappa _{\alpha ,-1}$ must have the opposite sign, as in (2.51). The first two terms in the curly brackets of (2.52) thus vanish, and the available energy becomes of order $\mathcal{O}(1)$ . We thus conclude that $0 \leqslant \eta \leqslant 2/3$ makes the available energy approach a constant value in the limit of large gradients in density and temperature, and this result holds for any magnetic geometry. It is worthwhile to note that, for linear stability of the curvature-driven ITG mode, it is sufficient to have $0 \leqslant \eta \leqslant 2/3$ (Biglari et al. Reference Biglari, Diamond and Rosenbluth1989), and this criterion thus has nonlinear significance too.

2.5. Perpendicular and parallel length scales

Having reached an understanding of the behaviour of the available energy in various limiting conditions, we proceed to discuss its possible relation to gyrokinetic turbulence, which requires one to specify the perpendicular length scales of the flux tube, $\Delta \psi$ and $\Delta \alpha$ (Mackenbach et al. Reference Mackenbach, Proll and Helander2022). The choice strongly affects the available energy, as the substitution $(\Delta \psi ,\Delta \alpha ) \rightarrow C (\Delta \psi ,\Delta \alpha )$ changes the available energy as $A \rightarrow C^4 A$ (see 2.21). As in TÆ, we take this length scale to be the correlation length, the length scale over which gradients in density and temperature are flattened (in turn liberating thermal energy). In situations where micro-turbulence is dominant, the correlation length is of the order of the gyroradius $\rho _{\text{gyro}} = m v_T / |q|B_{0}$ , where $q$ is the charge of the particle, $v_T = \sqrt {2T_0/m}$ is the thermal velocity and $B_0$ is a reference magnetic field strength, typically taken to be $B_0 = 2\psi _a/a_{\text{minor}}^2$ with $\psi _a$ being the toroidal flux over $2\pi$ at the last-closed flux surface, and $a_{\text{minor}}$ being the minor radius of the device. As in TÆ, we estimate these length scales by defining

(2.53) \begin{align} \rho (\psi ) = \sqrt {\psi /\psi _a}, \qquad s = \rho (\psi _0) \alpha , \end{align}

where $\rho$ is the typical normalised radial coordinate and $s$ is a newly introduced binormal coordinate (reminding ourselves of the fact that $\alpha = \theta - \iota \varphi$ ). Consequently, setting $\Delta \rho = \Delta s = \rho _{\text{gyro}}/a_{\text{minor}} = \rho _*$ implies that the correlation length is a (normalised) gyroradius, which may be related to $\Delta \psi$ and $\Delta \alpha$ in (2.21).

The available energy depends on a parallel length scale too, in non-isodynamic magnetic fields that have non-zero magnetic shear. This can readily be seen by investigating the gradient of $\alpha$ ,

(2.54) \begin{equation} \boldsymbol{\nabla }\alpha = \boldsymbol{\nabla }(\theta - \iota _0 \varphi ) - \varphi \iota _0' \boldsymbol{\nabla }\psi , \end{equation}

where we have defined $\iota _0 = \iota (\psi _0)$ and $\iota _0'= \iota '(\psi _0)$ . Hence, $\omega _\alpha$ (related to the gradient drift $\boldsymbol{v}_{\boldsymbol{\nabla }}\boldsymbol{\cdot }\boldsymbol{\nabla }\alpha$ via 2.10) grows without bound as $\varphi \to \pm \infty$ if $\iota _0' \neq 0$ and $\boldsymbol{v}_{\boldsymbol{\nabla }}\boldsymbol{\cdot }\boldsymbol{\nabla }\psi \neq 0$ . Consequently, the available energy given in (2.21) will depend on the considered range of $\ell$ , which we denote as $\Delta \ell$ , and one needs to choose an appropriate range a priori.Footnote ⁵ If the device is a tokamak, a natural parallel length scale is the connection length, i.e. the total length required to make one full poloidal turn,

(2.55) \begin{equation} \Delta \ell = \int _{-\pi }^{\pi } \frac {B}{\boldsymbol{B} \boldsymbol{\cdot }\boldsymbol{\nabla }\theta }\, \mathrm{d} \theta . \end{equation}

In a stellarator, the appropriate choice of length scale is less evident, but one is helped by the fact that the magnetic shear is typically much smaller than in tokamaks. As a consequence, in stellarators, the choice is less important and we shall simply take it to be a fixed parallel length, i.e. $\Delta \ell = \mathrm{constant}$ .

3.1. Parallel and phase-space structure of the available energy

Here, we first investigate the local properties of the available energy, i.e. how it depends on the local drifts at each $\ell$ . This may be achieved by considering a vanishingly small parallel extent $\Delta \ell$ . The available energy then becomes a function of four parameters

(3.3) \begin{equation} \widehat {A} = \widehat {A}(\hat {\omega }_n,\hat {\omega }_T,\hat {\omega }_\psi ,\hat {\omega }_\alpha ), \end{equation}

and we investigate its dependence here. The result is plotted in figure 2, where we have furthermore included a plot of the solutions to (2.17), $\hat {\kappa }_\alpha$ and $\hat {\kappa }_\psi$ .

Let us take note of two general symmetries here. First, note invariance of the available energy and the density constraint equations under $(\kappa _\psi ,\omega _\psi ) \rightarrow -(\kappa _\psi ,\omega _\psi )$ . Thus, the available energy depends only on the magnitude of the radial drift and not its sign. In terms of figure 2, this manifests as a mirror symmetry about the $\hat {\omega }_\psi = 0$ axis. Furthermore, the substitution $(\omega _\alpha ,\kappa _\alpha ,\omega _*^T) \rightarrow -(\omega _\alpha ,\kappa _\alpha ,\omega _*^T)$ similarly leaves the density constraint equations and the available energy unchanged. This tells us that the available energy depends only on the relative sign of $\omega _*^T$ and $\omega _\alpha$ . Indeed, if one flips the sign of $\omega _*^T$ , one finds that figure 2 is mirrored with respect to the $\hat {\omega }_\alpha = 0$ axis.

Figure 2.

Available energy $\widehat {A}$ (with white contour lines added so different regions can be seen clearly), and the solutions to the density constraint equations $\hat {\kappa }_\alpha$ and $\hat {\kappa }_\psi$ , as a function of the radial ( $\hat {\omega }_{\alpha }$ ) and binormal derivative ( $\hat {\omega }_\psi$ ) of the magnetic field. Here, $\hat {\omega }_T = 10$ and $\hat {\omega }_n = 2$ . Note that since $\hat {\omega }_T\sim - \partial _\rho \ln T$ (and similarly for $\hat {\omega }_n$ ) and furthermore, $\omega _\alpha \sim \partial _\rho \ln B$ . Positive/negative $\omega _\alpha$ means the magnetic field gradient is anti-/co-aligned with the density and temperature gradient.

Now, investigating the available energy in detail, it is clear that the most available energy is present for negative $\hat {\omega }_\alpha$ , i.e. if the density and temperature gradient are anti-aligned with the magnetic field gradient as in the bad-curvature region. The good-curvature region, conversely, has only very little available energy present. It is furthermore interesting to notice that if the magnitude of the binormal drift becomes sufficiently large in the bad-curvature region, the available energy decreases. Physically, this may be interpreted as follows: if the particles are drifting too fast, they are unable to couple to and resonate with the drift-wave, resulting in low available energy. For $\omega _\psi = 0$ , we see that there exists a ‘sweet-spot’ where the drifting particles are maximally resonant with the drift wave, freeing up the most available energy, at $\hat {\omega }_\alpha \approx -1.3$ (though this value of course depends on the gradients in density and temperature). The dependence of the available energy on $\hat {\omega }_\psi$ is less intricate: increasing the radial drift simply increases the total available energy.

Having investigated the dependence of the available energy at fixed $\ell$ , we now investigate its dependence on $\hat {v}_\|$ and $\hat {v}_\perp$ , i.e. we investigate $\int \mathrm{d} \psi\, \mathrm{d} \alpha\, \mathrm{d} \ell \; \epsilon (f_M - F) / E_T \rho _*^2$ for $\Delta \ell \rightarrow 0$ , which becomes (see (B8))

(3.4) \begin{equation} \begin{aligned} \widehat {\mathcal{A}}(\hat {v}_\perp , \hat {v}_\|) =& \frac {e^{ - \hat {v}_\|^2 - \hat {v}_\perp ^2} \hat {v}_\perp }{3 \sqrt {\pi }} \\ &\times \Bigg [ \sqrt { (\hat {v}_\perp ^2 \hat {\omega }_\alpha + \hat {\kappa }_\alpha )^2 + (\hat {v}_\perp ^2\hat {\omega }_\psi - \hat {\kappa }_\psi )^2} \sqrt {(\hat {v}_\perp ^2\hat {\omega }_\alpha - \hat {\omega }_*^T )^2 + (\hat {v}_\perp ^2\hat {\omega }_\psi )^2} \\ &-(\hat {v}_\perp ^2\hat {\omega }_\alpha + \hat {\kappa }_\alpha ) (\hat {v}_\perp ^2\hat {\omega }_\alpha - \hat {\omega }_*^T) - (\hat {v}_\perp ^2\hat {\omega }_\psi - \hat {\kappa }_\psi ) (\hat {v}_\perp ^2 \hat {\omega }_\psi ) \Bigg ] , \end{aligned} \end{equation}

so that $ \widehat {A} = \int \mathrm{d} v_\perp \,\mathrm{d} v_\| \; \widehat {\mathcal{A}}(\hat {v}_\perp , \hat {v}_\|)$ . This allows us to probe the structure of the ground state, its (an)isotropic features and where it differs the most from the initial Maxwellian. A plot of the phase-space structure is given in figure 3, where $\eta \gt 2/3$ . It is clear that the overall phase-space structure is rich, showing a region of much available energy near the origin and in a ‘bean’ at $\hat {v}_\perp \approx 2$ . The structure is furthermore fairly isotropic in the bad-curvature region with small drifts compared with the density and temperature gradients (negative and small $\hat {\omega }_\alpha$ , as in the top-right plot), and gains more anisotropic features in the good-curvature region (positive $\hat {\omega }_\alpha$ ). Note that there is nonetheless energy available, even in the good-curvature region. The anisotropic features partly answer how any electrostatic instability can exist in the good-curvature region, just as Ivanov et al. (Reference Ivanov, Luhadiya, Adkins and Schekochihin2025) argued that such instabilities cannot exist under certain assumptions of isotropy. The phase-space structure indicates that isotropy is broken, opening up new routes for an instability to exist.

All in all, we see that the available energy exhibits rich behaviour in terms of its dependencies on magnetic geometry and its velocity-space coordinates. The full calculation of the available energy is a weighted integral of available energies at different $\ell$ , as $\hat {\omega }_\alpha (\ell )$ and $\hat {\omega }_\psi (\ell )$ describe a trajectory in the space of figure 2 (e.g. the trajectory of a zero-shear large-aspect-ratio circular tokamak would be a circle centred at the origin, $\hat {\omega }_\alpha \sim -\cos \theta$ and $\omega _\psi \sim \sin \theta$ with $\theta =0$ corresponding to the outboard mid-plane). As we have shown here, regions of bad curvature will contribute the most to such a trajectory in $(\hat {\omega }_\alpha ,\hat {\omega }_\psi )$ -space.

3.2. Comparison to nonlinear local gyrokinetics

In this section, we compare the available energy with results from local, gyrokinetic, ITG-driven turbulence simulations. A database of such simulations in ‘random’ stellarator geometries with adiabatic electrons has recently been constructed by Landreman et al. (Reference Landreman, Choi, Alves, Balaprakash, Churchill, Conlin and Roberg-Clark2025), where the results from ${\gt} 2\, {\times}\, 10^5$ simulations are stored. The magnetic geometry of the stellarators in that article are sampled in three distinct manners, which we briefly review (note that one can sample multiple flux tubes from a single geometry).

1. (i) One sampling method imposes that the stellarators are heliotron-like (i.e. the cross-section is a rotating ellipse), where the number of field periods ( $N_{\text{fp}}$ ), the aspect ratio, the elongation, the axis torsion and the plasma beta are sampled randomly. Roughly, $27\,\% / 25\,\%$ of the stellarators/flux tubes are sampled in this manner. We note that these geometries generally do not confine trapped particles.

2. (ii) Another method simply samples equilibria from the quasr database (Giuliani Reference Giuliani2024), in which (approximately) quasi-symmetric geometries are stored (i.e. where $B$ has an invariant direction on a flux surface in Boozer-coordinates, resulting in confined trapped particles). This constitutes some $19\,\%/24\,\%$ of the stellarators/flux tubes.

3. (iii) The final method samples random boundary shapes by generating its Fourier modes from probability distributions that have been fit to a dataset of known stellarator shapes. Some $54 \,\% / 51\,\%$ of the stellarators/flux tubes are sampled in this manner.

The combined set of stellarators is diverse, and includes aspect ratios ranging from 2.9 to 10, and volume-averaged betas from 0 % to 5 %. The number of field periods, $N_{\text{fp}}$ , varies between $2$ and $8$ , which results in varying magnetic-well lengths. We have expanded this database with nonlinear simulations performed in ‘random’ tokamak geometries. These simulations, like those by Landreman et al. (Reference Landreman, Choi, Alves, Balaprakash, Churchill, Conlin and Roberg-Clark2025), were performed with the gx-code (Mandell et al. Reference Mandell, Dorland, Abel, Gaur, Kim, Martin and Qian2024), assuming adiabatic electrons and accounting for an ion-temperature and density gradient. More details, such as the construction of the probability density in parameter space where tokamaks are sampled using the analytical Miller equilibrium (Miller et al. Reference Miller, Chu, Greene, Lin-Liu and Waltz1998) and a numerical convergence study, are provided in Appendix D.

Both simulation sets contain a fixed-gradient subset with $\partial _\rho \ln T= -3.0$ and $\partial _\rho \ln n= -0.9$ , where $\rho = r/a_{\text{minor}}$ with $r$ being the minor radial coordinate labelling the flux surfaces and $a_{\text{minor}}$ is the minor radius of the device. The complement of the fixed-gradient subset has gradients that are chosen randomly, according to a sampling procedure. These two subsets are constructed so that one may systematically assess both the effect of the geometry and the effect of the gradients in density and temperature separately. To aid in the latter, some simulations are performed in the same geometry twice: once with the density and temperature gradients fixed, and once with randomly sampled values.

Finally, it should be noted that a crude argument on the basis of energy conservation relates the available energy to the energy flux (as also argued in TÆ). Consider the invariant energy in the long wavelength, electrostatic limit of gyrokinetics (Dubin et al. Reference Dubin, Krommes, Oberman and Lee1983; Brizard & Hahm Reference Brizard and Hahm2007),

(3.5) \begin{equation} \begin{aligned} E_{\text{tot}} \approx & \int \mathrm{d} \boldsymbol{x} \; \epsilon (\delta f) + \frac {m_i n_0}{2 } \int \mathrm{d} \boldsymbol{r}\; \left | \boldsymbol{v}_{\boldsymbol{E}} \right |^2, \end{aligned} \end{equation}

where $\delta f$ is the fluctuating part of the ion distribution function and $\boldsymbol{v}_{\boldsymbol{E}}$ is the $\boldsymbol{E} \times \boldsymbol{B}$ drift velocity. The available energy bounds by how much the first term of (3.5) can decrease, and it thus provides an upper bound on the kinetic ‘sloshing’ energy (the second term). Thus, one finds $\boldsymbol{v}_{\boldsymbol{E}}^2 \sim A$ , and since the energy flux goes as $Q \sim \epsilon (\delta f) |\boldsymbol{v}_{\boldsymbol{E}}|$ (to be made precise in § 3.2.1), we have $Q \sim A^{3/2}$ . Results shall often be compared with this expected scaling.

3.2.1. Fixed-gradient subset

We start our analysis by focussing on the fixed-gradient subset. The result of this analysis is plotted in figure 4, where the time-averaged nonlinear radial energy flux is plotted against the available energy. The former is defined as

(3.6) \begin{equation} Q = \frac {1}{t_{\text{sat}}} \int _{t_{\text{sat}}} \mathrm{d} t \; \left \langle \frac {1}{L_\rho L_s} \int \mathrm{d} \rho\, \mathrm{d}s\, \mathrm{d} \boldsymbol{v} \; \epsilon \left ( \delta f \right ) \boldsymbol{v_E} \boldsymbol{\cdot }\boldsymbol{\nabla }\rho \right \rangle , \end{equation}

where the radial and binormal coordinates are defined in (2.53), and $t_{\text{sat}}$ is a time window in which the turbulence is saturated. Furthermore, $\delta f$ is the fluctuating part of the distribution function, the electric-field drift $\boldsymbol{v}_{\boldsymbol{E}}$ is the $\boldsymbol{E} \times \boldsymbol{B}/B^2$ , and $L_\rho$ and $L_s$ measure the perpendicular dimensions of the flux tube so that the integral over the perpendicular coordinates is averaged. Finally, the angular brackets define a flux-surface average $\langle \ldots \rangle = \int \mathrm{d} \ell \; \ldots B^{-1}/\int \mathrm{d} \ell \; B^{-1}$ .

Figure 4.

Scatter plot comparing the available energy and the nonlinear, time-averaged, radial energy flux (from nonlinear gyrokinetic simulations) for the fixed-gradient subset (i.e. $\partial _\rho \ln T = -3.0$ and $\partial _\rho \ln n = -0.9$ ). The different outer panels correspond to different number of field-periods, and all the field-periods are plotted together in the centre plot. The black dotted line shows the expected power-law $Q \propto \widehat {A}^{3/2}$ , and the colour of the scatter indicates whether the geometry comes from the quasr-database (blue) or not (red). Furthermore $N_{\textrm {fp}}=0$ corresponds to the tokamak case.

Figure 5.

Left: scatter plot of the available energy against the nonlinear radial energy flux for tokamaks alone, where the points are coloured according to their connection length. Right: scatter plot of the (logarithm of the) ratio $Q/A^{3/2}$ and the connection length. Both show results for the fixed-gradient subset.

It is evident that, though there is certainly correlation if $N_{\text{fp}}$ is small, for high $N_{\text{fp}}$ , there is very little correlation in the bulk of the dataset, and the quasr equilibria showcase an even smaller amount of correlation (quantitative values will be provided in § 3.2.3). A likely reason for this mismatch is the dependence of the importance of parallel dynamics (which we have ignored) on the periodicity, as a typical magnetic well size scales as $L_{\text{well}} \sim R/N_{\text{fp}}$ . Approximating the parallel wavenumber as $k_{\|} \sim 1/L_{\text{well}} \sim N_{\text{fp}}/R$ , one violates $k_{\|} v_{T} \ll \boldsymbol{k}_{\perp } \boldsymbol{\cdot }\boldsymbol{v}_{D}$ to an increasing degree when raising $N_{\text{fp}}$ , invalidating a central assumption of the available energy derived. To assess the validity of this hypothesis, we furthermore include a plot of the tokamak dataset with fixed gradients (see figure 5), where the parallel length scale can readily be defined via the connection length $L_{\text{con}} \equiv \int \mathrm{d} \theta \; B/\boldsymbol{B} \boldsymbol{\cdot }\boldsymbol{\nabla }\theta$ , the total length required to make one poloidal turn. It may be seen that the simulations with a very short connection length (and thus high $k_\| \sim 1/L_{\text{con}}$ ) deviate strongly from the predicted power law. Conversely, simulations that are characterised by long connection lengths stay closer to the predicted power law, although certainly with deviations. This corroborates the idea that the parallel dynamics, i.e. non-negligible ion motion over a connection length, affects the turbulence and violates our key assumption concerning the available energy. Using the right-most plot of figure 5, we see that the decrease in $\log _{10}Q/\widehat {A}^{3/2}$ with decreasing logarithmic connection length is roughly linear, until $\log _{10}( L_{\text{con}}/a_{\text{minor}}) \approx 0.7$ is approached from above, after which the heat-flux drops strongly. This gives a quantitative estimate of when parallel dynamics is non-negligible (namely at $L_{\text{con}}/a_{\text{minor}} \approx 5$ ) and for shorter connection lengths, the available energy is likely no longer of practical use (for this combination of gradients in tokamaks).

Figure 6.

Distribution of $\log _{10} \widehat {A}$ for ‘stable’ ( $Q\leqslant 0.1$ , blue) and ‘unstable’ ( $Q\gt 0.1$ , orange) simulation data. The $y$ -axis denotes the probability density so that the total area under the blue and orange distributions evaluate to unity.

Figure 7.

Scatter plot of the available energy against the nonlinear, time-averaged, radial energy flux for the varying-gradient subset. The outer panels correspond to different number of field-periods, and all the field-periods are plotted together in the centre plot. The black dotted line shows the expected power-law $Q \propto \widehat {A}^{3/2}$ , and the colour of the scatter indicates whether the geometry comes from the quasr-database (blue) or not (red).

Figures 4 and 5 do not show the full extent of the dataset, as stable simulations are characterised by exponential decay of the heat flux, and thus depend critically on the decay rate, simulation time and initial condition, rendering the precise value of $Q$ meaningless. To investigate whether available energy may be useful in predicting the nonlinear stability of a simulation, we plot histograms of $\widehat {A}$ for $Q \leqslant 0.1$ (‘stable’ cases) and $Q\gt 0.1$ (‘unstable’ cases). The result is shown in figure 6. One sees that for all $N_{\text{fp}}$ , except tokamaks ( $N_{\text{fp}}=0$ ), sufficiently small available energy implies that the simulation is stable. This bolsters some hope for stellarator equilibrium optimisation purposes, as optimising for sufficiently low available energy will then likely imply that the device is nonlinearly stable. Upon closer scrutiny, one reason for the mismatch in tokamaks becomes apparent. As may be seen in figure 5, simulations with very short connection lengths tend to be stable and the available energy is a poor predictor in this regime. Furthermore, the shortest connection lengths $L_{\text{con}} \sim q R$ are characterised by small major radius and thus strong curvature drive as $\omega _\alpha \sim 1/R$ . The points with the strongest curvature drive have large available energy (as in a strongly driven limit $A \sim \omega _\alpha$ ), but are in reality stabilised by the short connection length, thus explaining the mismatch. Indeed, repeating the same analysis on the tokamak set but only allowing for data with $3\lt R/a_{\text{minor}}\lt 4$ , this mismatch is reduced somewhat, and the stable and unstable distributions end up looking similar, making available energy at best a poor predictor of nonlinear stability in tokamaks.

A similar phenomenon seems to occur for the stable simulations with high available energy, present for $N_{\text{fp}} = 8$ and $N_{\text{fp}} = 9$ . These data are characterised by having a magnetic field with a high mirror ratio (leading to stronger mirror forces) and a short dominant parallel wavelength (increasing $k_\|$ ), both enhancing the parallel dynamics. The available energy strongly overestimates the turbulent transport of these points since due to the large variation in magnetic field, the $\boldsymbol{\nabla }B$ -drifts are significantly enhanced, increasing the curvature drive. The parallel dynamics, however, seemingly stabilises turbulence resulting in the mismatch.

3.2.2. Varying-gradient subset

We next consider the data with varying gradients, excluding the fixed-gradient subset discussed in the previous section. The result is shown in figure 7. It is evident that the scatter is much smaller in this case, highlighting that available energy is able to capture the dependence of energy flux on the gradients in density and temperature. As before, the scatter for the tokamak subset is worse than the stellarators. To further assess how the gradient dependency is captured, we include a scatter of the data against these gradients. The result is shown in figure 8, where one sees expected trends: increasing the strength of the temperature gradient tends to increase the available energy and nonlinear energy flux. However, if $\eta$ is sufficiently close to or smaller than $2/3$ , both the nonlinear energy flux and available energy are significantly reduced. It is noteworthy that the critical $\eta$ below which the energy flux is stable is approximately $\eta \approx 1$ , which could be due to other instabilities (e.g. the slab branch) dominating or the Dimits-shift stabilising the turbulence. The current available energy model does not capture such effects, and instead shows a significant reduction if $\eta \lt 2/3$ as expected from the asymptotic theory discussed in § 2.4.

Figure 8.

Top row: the energy flux (left) and available energy (right) scattered against the density ( $\hat {\omega }_n$ ) and temperature ( $\hat {\omega }_T$ ) gradients. Bottom row: the energy flux and available energy scatter against $\arctan (\hat {\omega }_T/\hat {\omega }_n)$ . The $\eta =2/3$ line is added in all plots as a dashed red line, and an $\eta =1$ line is added as a dashed grey line. Data of the varying-gradient subset.

Table 1.

Various quantitative measures of correlation, for datasets of stellarators and tokamaks, split between the fixed ( $\mathfrak{F}$ ) and random ( $\mathfrak{R}$ ) gradient subsets. In the first column, the correlation measure and analysis type, regression (regr.) or classification (clas.), is stated. The ‘Landreman’ column has best scoring values of Landreman et al. (Reference Landreman, Choi, Alves, Balaprakash, Churchill, Conlin and Roberg-Clark2025). The ‘Stellarators’ column analyses the predictive capabilities of $\widehat {A}$ for stellarators alone. The ‘Tokamaks’ column does the same for tokamaks alone. For the XGBoost analyses, each value is the mean score on held-out data with fivefold cross-validation.

3.2.3. Quantitative comparison

Some quantitative measures of correlation are computed as a final step allowing for direct comparison to results found by Landreman et al. (Reference Landreman, Choi, Alves, Balaprakash, Churchill, Conlin and Roberg-Clark2025), who constructed a set of scalar ‘features’ that correlate especially well with the energy flux. Three measures of correlation are investigated.

1. (i) The Spearman coefficient, which lies between $-1$ and $1$ , where a higher absolute value indicates tighter correlation and a perfect monotonic relationship between two variables returns a magnitude of $1$ , see Spearman (Reference Spearman1904). For analyses involving the Spearman coefficient, whenever $Q\lt 0.1$ , it is set to zero as the exact value of such ‘stable’ simulations is meaningless, though we note that leaving $Q$ unchanged does not alter the results significantly. The Spearman coefficient then measures the correlation between $\widehat {A}$ and $Q$ (or equivalently between $\ln \widehat {A}$ and $\ln Q$ ).

2. (ii) For regression analyses, the coefficient of determination $R^2 \in [0,1]$ is used where a higher value indicates better correlation. The regression is performed between the true value $\ln Q$ and its predicted value $\ln Q_p$ . To predict the energy flux, the machine learning gradient-boosted decision-tree package XGBoost is used (Chen & Guestrin Reference Chen and Guestrin2016), which attempts to find a function $X$ that maps inputs $\boldsymbol{I}$ to the output $\ln Q$ , resulting in the prediction $X(\boldsymbol{I}) = \ln Q_p$ . For the fixed-gradient subset the input has only one component, namely $\boldsymbol{I}=\{ \widehat {A} \}$ , whereas the varying gradient subset uses $\boldsymbol{I} = \{\widehat {A},\, \hat {\omega }_T,\, \hat {\omega }_n \}$ .

3. (iii) Finally, the log-loss (also known as cross-entropy) is used for scoring classification predictions of (in)stability, where a lower value indicates more accurate predictions. The classification task attempts to predict whether a simulation is ‘stable’ ( $ Q \gt 0.1$ ) or ‘unstable’ ( $Q \leqslant 0.1$ ). This prediction is again performed using XGBoost, where the fixed-gradient subset uses only $\boldsymbol{I}=\{ \widehat {A} \}$ as an input and the random-gradient subset uses $\boldsymbol{I} = \{\widehat {A},\, \hat {\omega }_T,\, \hat {\omega }_n \}$ , resulting in the mapping $X(\boldsymbol{I})=\text{stable/unstable}$ . Accuracy scores () for the classification analyses are included as well, but these can be misleading as most of the data are unstable so that simply predicting each datum to be unstable can lead to high accuracy.

Finally, the tokamak and stellarator subsets are analysed separately due to their differing behaviours. All these correlation measures are compared with their best-scoring counterpart of Landreman et al. (Reference Landreman, Choi, Alves, Balaprakash, Churchill, Conlin and Roberg-Clark2025), who found that the nonlinear energy flux correlates especially well with

(3.7) \begin{equation} f_Q = \mathrm{mean} \left ( \left [ \varTheta\! \left ( \boldsymbol{B} \times \boldsymbol{\kappa } \boldsymbol{\cdot }\boldsymbol{\nabla }s \right ) + 0.2 \right ] \frac {|\boldsymbol{\nabla }\rho |^3}{B} \right ), \end{equation}

where $\boldsymbol{\kappa } = \boldsymbol{b} \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{b}$ is the curvature vector and the operator $\mathrm{mean}$ takes the average of the input array. As such, the correlation analysis with Spearman coefficient and the regression analysis is performed as mentioned, where $f_Q$ replaces $\widehat {A}$ . To predict stability instead, the feature

(3.8) \begin{equation} f_{\text{stab }} = \mathrm{mean} \left ( \left [ \varTheta\! \left ( \boldsymbol{B} \times \boldsymbol{\nabla }B \boldsymbol{\cdot }\boldsymbol{\nabla }s \right ) + 0.4 \right ] \frac {|\boldsymbol{\nabla }\rho |}{\sqrt {B}} \right ) \end{equation}

is used. The analysis for stability is then performed as mentioned, but with $f_{\text{stab}}$ replacing $\widehat {A}$ .

The results are presented in table 1, where four columns are presented. The first column states the correlation measure and whether it was performed on the fixed ( $\mathfrak{F}$ ) or random/varying ( $\mathfrak{R}$ ) gradient subset. It is furthermore stated if it is a regression (regr.) or classification (clas.) analysis. Best-scoring values for different correlation measures presented by Landreman et al. (Reference Landreman, Choi, Alves, Balaprakash, Churchill, Conlin and Roberg-Clark2025) are included in the second column, (‘Landreman’). The third column (‘Stellarators’) analyses the predictive capabilities of available energy for stellarators alone, where we stress again that the fixed gradient dataset only uses $\widehat {A}$ to predict $Q$ , and the varying-gradient database uses the density and temperature gradients as well to predict the energy flux. The fourth column (‘Tokamaks’) does the same for tokamaks, excluding stellarators. The correlation coefficients exhibit the behaviour discussed in previous sections: the available energy correlates relatively poorly with the energy flux in tokamaks in general. It correlates fairly well in stellarators if the gradients in density and temperature are allowed to vary, and it can furthermore make fairly accurate predictions of stability in stellarators although both effects can largely be explained by the ITG: the Spearman correlation coefficient between the ITG and nonlinear energy flux is $0.80$ (not included in table 1), and the high correlation of the available energy is thus largely explained by its dependence on it, increasing to $0.84$ due to its additional dependencies. Similarly, performing XGBoost regression and classification analyses on the varying-gradient subset where the only features retained are the gradients in density and temperature (also not included in table 1), gives $R^2 \approx 0.71$ and $\mathrm{accuracy} \approx 0.92$ , respectively. The available energy increases these values to $0.76$ and $0.94$ , respectively, showing that predictive capabilities increase, though only slightly. The most surprising finding is that available energy can predict stability in tokamaks fairly well too. This is likely a consequence of the anti-correlation found (with higher available energy corresponding to short connection lengths and, possibly, stability), and is thus likely strongly dependent on the probability space in which the tokamaks are sampled (e.g. sampling at constant connection length may lead to different results). We note in passing that similar implications on sampling likely exist in the stellarator database too (e.g. sampling geometries with high mirror force will likely give different outcomes).