2.1 Allowable domains

In the calculation of the ground state, we restrict our attention to a subregion of the plasma which has the shape of a slender flux tube aligned with the magnetic field, allowing us to approximate the distribution function by its first-order Taylor expansion in the directions perpendicular to the field. Only then does it seem possible to solve the integro-differential ground-state equation (1.7) in general. Somewhat surprisingly, it turns out that the calculation is only possible if the cross-section of the flux tube is elliptical, as we shall now see.

We first introduce the following re-scaled coordinates:

(2.1a)\begin{gather} x = (\psi - \psi_0)/{\mathrm \Delta} \psi, \end{gather}

(2.1b)\begin{gather}y = (\alpha - \alpha_0)/{\mathrm \Delta} \alpha. \end{gather}

Here, $\psi _0$ and $\alpha _0$ are the flux-surface label and field-line label of the magnetic field line defining the centre of the flux tube, and ${\mathrm \Delta} \psi$ and ${\mathrm \Delta} \alpha$ define its scale lengths in the $\psi$- and $\alpha$-directions, respectively. In contrast to the situation in most gyrokinetic simulations, the cross-section of the flux tube is not assumed to be rectangular. In these coordinates, we denote the first-order expansion of the distribution function $f$ for fixed values of $\mu$, $\mathcal {J}$, and $t$ by

(2.2)\begin{equation} f(x,y,\mu,\mathcal{J},t) \equiv f_0 + f_{x} x + f_{y} y. \end{equation}

For a given smooth distribution function, this representation should always be sufficiently accurate in the limit of small ${\mathrm \Delta} \psi$ and ${\mathrm \Delta} \alpha$ (i.e. for small enough $\psi - \psi _0$ and $\alpha - \alpha _0$ one can approximate $f \approx f_0 + (\psi - \psi _0) \partial _\psi f + (\alpha - \alpha _0) \partial _\alpha f$), but there is no guarantee that the ground state corresponding to such a function can be similarly approximated by its linear expansion. Indeed, this turns out to be true, in general, only if the domain $\varOmega$ in the $(x,y)$-plane under consideration is elliptical. As already mentioned, phase-space volume conservation dictates that the following integral is independent of time for all values of $\phi$:

(2.3)\begin{equation} \frac{\mathrm{d}}{\mathrm{d}t}\int_{\varOmega} \varTheta \left[\phi - f(x,y,\mu,\mathcal{J},t) \right] \mathrm{d}x\,\mathrm{d}y = 0. \end{equation}

In particular, this integral must be the same for the initial (given) distribution function and its corresponding ground state. We assume that both can be described by their linear approximations (2.2), whose gradients will, however, in general point in different directions (i.e. $\boldsymbol {\nabla }_\perp f \propto \hat {\boldsymbol {x}} f_x +\hat {\boldsymbol {y}} f_y$ can take on any value, with the hatted quantities denoting unit vectors). We thus require that the gradient of the ground state may point in an arbitrary direction and derive a constraint on the shape of $\varOmega$ from this condition.

It is useful to define a rotated coordinate system $(x,y)\mapsto (\tilde {x},\tilde {y})$

(2.4)\begin{equation} \left.\begin{gathered} \tilde{x} = x \cos \vartheta + y \sin \vartheta + x_{\rm min}(\vartheta), \\ \tilde{y} ={-} x \sin \vartheta + y \cos \vartheta + y_{\rm min}(\vartheta), \end{gathered}\right\} \end{equation}

such that $f$ assumes its smallest value, $f_\textrm {min}$, at $\tilde x = 0$ and the gradient of $f$ points in the direction of $\tilde x$

(2.5)\begin{equation} f = f_\mathrm{min} + f_{\tilde{x}} \tilde{x}. \end{equation}

In such a coordinate system the domain rotates as the function evolves, whilst fixing the minimal values of $x$ and $y$ to $\tilde {x}=\tilde {y}=0$, and this coordinate system will aid us in making statements about this domain shape. The initial state and the ground state thus correspond to two different values of $\vartheta$. The width of the domain $\varOmega$ in the $\tilde x$-direction is denoted by $2D(\vartheta )$, so that the maximum value of $f$ becomes $f_\mathrm {max} = f_\mathrm {min} + 2 D f_{\tilde {x}}$, and the distribution function may be written as

(2.6)\begin{equation} f = f_\mathrm{min} + \frac{\tilde{x}}{2D} \left( f_\mathrm{max} - f_\mathrm{min} \right), \end{equation}

where we note that $f_\mathrm {min}$ and $f_\mathrm {max}$ depend only on $\mu$ and $\mathcal {J}$ and are thus independent of $\vartheta$.

We now return to the problem at hand, codified by (2.3). The argument of the Heaviside function vanishes in the new coordinate system when

(2.7)\begin{equation} \tilde{x}(\phi) = 2D \frac{\phi - f_\mathrm{min}}{f_\mathrm{max}-f_\mathrm{min}}. \end{equation}

Let us now consider a value of $\phi$ very close to $f_\mathrm {min}$, i.e. $\phi = f_\mathrm {min} + |\epsilon |$ with $|\epsilon | \rightarrow 0^+$. Equation (2.7) now reduces to $\tilde {x}_\epsilon = 2D|\epsilon |/(f_\mathrm {max}-f_\mathrm {min})$. The integral in (2.3) is equal to the area of the subdomain of $\varOmega$ over which $f \le \phi$, i.e. the coloured region in figure 1. In the limit $\tilde {x}_\epsilon \ll D$, the radius of curvature, $R$, of the domain boundary in the region $\tilde x < \tilde x_\epsilon$ is approximately constant (if the boundary is sufficiently smooth). Hence, we may approximate the integral by instead considering the area of a circular segment with the same radius of curvature $R$

(2.8)\begin{align} \mathcal{A} & = 2\int_0^{\tilde{x}_\epsilon} \sqrt{2R\tilde{x} - \tilde{x}^2}\,\mathrm{d}\tilde{x} \nonumber\\ & \approx 2 \int_0^{\tilde{x}_\epsilon} \sqrt{2 R\tilde{x}}\,\mathrm{d}\tilde{x} = \frac{4}{3} \sqrt{2 R \tilde{x}_\epsilon^3}, \end{align}

where we have used the equation for a circle $\tilde {y}^2 = R^2 - (\tilde {x} - R)^2$. The next step is to realize that, as the distribution function evolves, its gradient may point in any direction in the $(x,y)$-plane, so that the angle $\vartheta$ may assume any value. A sketch of the domain at a different orientation is also given in figure 1. Both $D$ and $R$ depend on $\vartheta$, but they are not independent. Indeed, from (2.3) we have

(2.9)\begin{equation} R(\vartheta)D(\vartheta)^3 = C_\mathcal{A} = \mathrm{constant}, \end{equation}

where we have used the equation for the area found in (2.8). In particular, the radius of curvature is the same at the two extremal locations $\tilde {x}=0$ and $\tilde {x}=2D$, since a rotation of ${\rm \pi}$ leaves $D$ unchanged. Moreover, on purely geometrical grounds there is a relation between the second derivative of the domain width and the curvature, as shown in Appendix A,

(2.10)\begin{equation} \frac{\mathrm{d}^2 D}{\mathrm{d} \vartheta^2} = R - D. \end{equation}

Using the constancy of $RD^3$, we find that

(2.11)\begin{equation} \frac{\mathrm{d}^2 D}{\mathrm{d} \vartheta^2} + D = \frac{C_\mathcal{A}}{D^3} \implies \frac{\mathrm{d}}{\mathrm{d} \vartheta} \left[\left(\frac{\mathrm{d} D}{\mathrm{d} \vartheta}\right)^2 + D^2 + \frac{C_\mathcal{A}}{D^2} \right]=0. \end{equation}

As shown in Appendix A, this differential equation has the general solution

(2.12)\begin{equation} D(\vartheta) = (C^2-C_\mathcal{A})^{1/4}\sqrt{C + \cos 2 \vartheta}, \end{equation}

implying that the radius of curvature becomes

(2.13)\begin{equation} R(\vartheta) = \frac{C_\mathcal{A}}{(C^2-C_\mathcal{A})^{3/4}(C + \cos 2 \vartheta)^{3/2}}. \end{equation}

The radius of curvature for an ellipse is of the exact same form as (2.13) (as shown in Appendix A), implying that $\varOmega$ must be elliptical. We have thus found that the only domain shape in which Liouville's theorem can generally be satisfied is a slanted ellipse, if the distribution function is to be approximated by its first-order Taylor expansion in the coordinates perpendicular to the magnetic field.

Figure 1. A sketch of the area, coloured red, as a subset of the entire domain $\varOmega$. Done at two different times, so that the orientation of the domain $\varOmega$ is different.

2.2 The ground state and the available energy

As shown in the previous section the domain has to be a slanted ellipse, and here we shall derive what the ground state and Æ is in such a system. To this end, let us first realize that it is sufficient to calculate the Æ in an unslanted ellipse, which in the coordinates $x$ and $y$ of the previous section has the equation

(2.14)\begin{equation} \frac{x^2}{a^2} + \frac{y^2}{b^2} = 1, \end{equation}

with $a \leq b$. We may then apply a rotation operator to find the result of any slanted ellipse. Let us start, then, by considering (1.7) and one must evaluate integrals of the form

(2.15)\begin{equation} I[h] = \int_\varOmega \delta[h]\,\mathrm{d}\psi\,\mathrm{d}\alpha, \end{equation}

where $h(\psi,\alpha )$ is an arbitrary smooth function that vanishes in the centre of the domain, $h(0,0)=0$. This condition has to do with the fact that a constant Maxwellian is the ground state corresponding to an initial condition without gradients. As such, to lowest order the difference between the Maxwellian and the ground state should vanish, implying that $h(0,0)=0$ for the functions appearing in (1.7). Since we take the flux tube to be slender, we approximate the function $h$ as

(2.16)\begin{equation} h(\psi,\alpha) \approx x {\mathrm \Delta} \psi \partial_\psi h + y {\mathrm \Delta} \alpha \partial_\alpha h. \end{equation}

We next define the vector $\boldsymbol {p} = ( x, y )$, so that one can write the linear expansion of $h$ as

(2.17)\begin{equation} h \approx \boldsymbol{p} \boldsymbol{\cdot} \boldsymbol{n} |\partial_{\boldsymbol{p}} h|, \end{equation}

where we have defined $\boldsymbol {n} \equiv \partial _{\boldsymbol {p}} h/ |\partial _{\boldsymbol {p}} h|$, and derivatives of $h$ are evaluated at the centre of the flux tube. Equation (2.15) then becomes

(2.18)\begin{equation} I[h] = \frac{{\mathrm \Delta} \psi {\mathrm \Delta} \alpha }{|\partial_{\boldsymbol{p}} h|} \int_\varOmega \delta[\boldsymbol{p} \boldsymbol{\cdot} \boldsymbol{n}]\,\mathrm{d}x\,\mathrm{d}y. \end{equation}

The integral measures the length of the line where $h=0$ in the domain, which may be found as follows. We realize that this line satisfies the equation

(2.19)\begin{equation} y ={-}\frac{h_x}{h_y} x, \end{equation}

where we have used the notation $h_x \equiv \partial _x h$. We may now find its intersection with the ellipse given by (2.14), and one may readily verify that the points of intersection $(x,y)$ satisfy the equation

(2.20)\begin{equation} x^2 + y^2 = (ab)^2 \frac{(h_x)^2+( h_y)^2}{(a h_x)^2+(b h_y)^2}, \end{equation}

and the integral thus becomes

(2.21)\begin{align} I[h] & = \frac{{\mathrm \Delta} \psi {\mathrm \Delta} \alpha }{|\partial_{\boldsymbol{p}} h|} 2ab\sqrt{\frac{h_x^2 + h_y^2}{(a h_x)^2+(b h_y)^2}} \nonumber\\ & = {\mathrm \Delta} \psi {\mathrm \Delta} \alpha \frac{ 2 ab }{\sqrt{(ah_x)^2+(bh_y)^2}}. \end{align}

Our final step is to rotate the domain, which is equivalent to rotating the vector $\partial _{\boldsymbol {p}} h$, for which one can use the rotation mapping

(2.22a)\begin{gather} h_x \mapsto h_x \cos \vartheta - h_y \sin \vartheta, \end{gather}

(2.22b)\begin{gather}h_y \mapsto h_x \sin \vartheta + h_y \cos \vartheta, \end{gather}

where $\vartheta$ measures the angle of rotation. Combining results of (2.21) and (2.22), we find

(2.23)\begin{equation} \left.\begin{gathered} I[h] = \frac{{\mathrm \Delta} \psi {\mathrm \Delta} \alpha}{ b \sqrt{ h_x^2 + h_y^2 - \mathcal{E}_\vartheta^2 } }, \\ \mathcal{E}_\vartheta(h_x,h_y) \equiv e \left( h_x \cos \vartheta - h_y \sin \vartheta \right), \end{gathered}\right\} \end{equation}

where the eccentricity $e \in [0,1)$ is defined as $e^2 \equiv 1 - a^2/b^2$. We are now in a position to evaluate the ground state; first choose the initial distribution function to be a Maxwellian (as usually done in gyrokinetic simulations)

(2.24)\begin{equation} f_M = n(\psi) \left(\frac{m}{2 {\rm \pi}T(\psi)} \right)^{3/2} \exp({-\epsilon/T(\psi)}). \end{equation}

Here, $n(\psi )$ is the number density, $m$ is the electron mass and $T(\psi )$ is the electron temperature. The spatial derivatives of the distribution and energy function can be related to the bounce-averaged precession frequencies of the electron orbits (as shown in Appendix B)

(2.25a)\begin{gather} \left(\frac{\partial \epsilon}{\partial \psi}\right)_{\mu,\mathcal{J},\alpha} ={+} q \omega_\alpha, \end{gather}

(2.25b)\begin{gather}\left(\frac{\partial \epsilon}{\partial \alpha}\right)_{\mu,\mathcal{J},\psi} ={-} q \omega_\psi, \end{gather}

(2.25c)\begin{gather}\left(\frac{\partial f_{M}}{\partial \psi}\right)_{\mu,\mathcal{J},\alpha} ={+} q \frac{f_{M,0}}{T_0} \left( \omega_*^T - \omega_\alpha \right), \end{gather}

(2.25d)\begin{gather}\left(\frac{\partial f_{M}}{\partial \alpha}\right)_{\mu,\mathcal{J},\psi} ={+} q \frac{f_{M,0}}{T_0} \omega_\psi. \end{gather}

Here, quantities with the subscript zero are quantities that are evaluated at the centre of the flux tube, $q$ is the electron charge, $\omega _\psi$ is the drift (precession) frequency in the $\psi$ direction, $\omega _\alpha$ is the corresponding frequency in the $\alpha$ direction and $\omega _*^T$ is the electron diamagnetic frequency. With these expressions, the numerator in (1.7) for the ground state becomes

(2.26)\begin{equation} \iint \delta \left[w - \epsilon_0(\mu,\mathcal{J}) - q \omega_\alpha \psi + q \omega_\psi \alpha \right] \mathrm{d}\psi \,\mathrm{d}\alpha= \frac{{\mathrm \Delta} \psi {\mathrm \Delta} \alpha }{|qb| \sqrt{(\omega_\alpha {\mathrm \Delta} \psi )^2 + (\omega_\psi {\mathrm \Delta} \alpha )^2 - (\mathcal{E}_\vartheta)^2 }}. \end{equation}

The denominator of (1.7) can be calculated in a similar manner. We conclude that the ground-state distribution function has the following derivative:

(2.27)\begin{equation} \left.\begin{gathered} \left(\frac{\partial f_g(w,\mu,\mathcal{J})}{\partial w}\right)_{\mu,\mathcal{J}} ={-} \frac{f_{M,0}}{T_0} F, \\ F \equiv \frac{\sqrt{(\omega_*^T - \omega_\alpha)^2 ({\mathrm \Delta} \psi)^2 + \omega_\psi^2 ({\mathrm \Delta} \alpha)^2 - \mathcal{E}_\vartheta ({\mathrm \Delta} \psi [\omega_*^T-\omega_\alpha], {\mathrm \Delta} \alpha \omega_\psi )^2 }}{\sqrt{\omega_\alpha^2 ({\mathrm \Delta} \psi)^2 + \omega_\psi^2 ({\mathrm \Delta} \alpha)^2- \mathcal{E}_\vartheta ({\mathrm \Delta} \psi \omega_\alpha, -{\mathrm \Delta} \alpha \omega_\psi )^2}}. \end{gathered}\right\} \end{equation}

Since the right-hand side is negative, the solution decreases with increasing energy, which is a requirement for any ground state. In the limit of omnigeneity (thus $\omega _\psi = 0$) one may verify that the terms involving $\mathcal {E}_\vartheta$ cancel, and one retrieves the expression previously found by Helander (Reference Helander2020). To evaluate the Æ, we must find the difference in energy between the initial distribution function and the ground state

(2.28)\begin{equation} A = \int \epsilon \left( f_i - f_g \right) \sqrt{g}\,\mathrm{d}\psi\,\mathrm{d}\alpha\,\mathrm{d}\mu\,\mathrm{d}\mathcal{J}. \end{equation}

Here, $\sqrt {g}$ denotes the phase-space Jacobian, which reduces to $\sqrt {g}= 4 {\rm \pi}/ m^2$ (Helander Reference Helander2017). Moreover, $f_i$ and $f_g$ can be approximated by their Taylor expansions

(2.29)\begin{equation} f_g\left[\epsilon(\psi,\alpha,\mu,\mathcal{J}),\mu,\mathcal{J}\right] = f_{i,0} + \left(\frac{\partial \epsilon}{\partial \psi} \psi + \frac{\partial \epsilon}{\partial \alpha} \alpha\right) \left.\frac{\partial f_g(w,\mu,\mathcal{J})}{\partial w}\right|_{w = \epsilon_0}+ \cdots. \end{equation}

It is useful to realize that the total number of particles within the domain for each pair of $\mu$ and $\mathcal {J}$ is conserved, giving

(2.30)\begin{equation} \iint \left( f_i - f_g \right) \mathrm{d}\psi \,\mathrm{d}\alpha = 0, \end{equation}

which in turn implies that the Æ can be calculated, to leading order, as

(2.31)\begin{equation} A = \frac{4 {\rm \pi}{\mathrm \Delta} \psi {\mathrm \Delta} \alpha}{m^2} \int \left( {\mathrm \Delta} \psi^2 \frac{\partial \epsilon}{\partial \psi} \frac{\partial (f_i - f_g)}{\partial \psi} x^2 + {\mathrm \Delta} \alpha^2 \frac{\partial \epsilon}{\partial \alpha} \frac{\partial (f_i - f_g)}{\partial \alpha} y^2 \right) \mathrm{d}x\,\mathrm{d}y\,\mathrm{d}\mu\,\mathrm{d}\mathcal{J}. \end{equation}

The integration across $x$ and $y$ is readily carried out by rotating the function $x^2$ or $y^2$ by some angle, and integrating this function over the domain whose boundary is given in (2.14). This may be readily calculated as

(2.32a)\begin{gather} \int (x \cos \vartheta - y \sin \vartheta )^2 \,\mathrm{d}x \,\mathrm{d}y = \frac{{\rm \pi} a^2b^2}{4} \left( \cos^2 \vartheta \sqrt{1-e^2} + \frac{\sin^2 \vartheta}{\sqrt{1 - e^2}} \right) \equiv \frac{{\rm \pi} a^2 b^2}{4} \hat{\mathcal{E}}_{\vartheta}^2, \end{gather}

(2.32b)\begin{gather}\int (y \cos \vartheta + x \sin \vartheta )^2 \,\mathrm{d}x \,\mathrm{d}y = \frac{{\rm \pi} a^2b^2}{4} \left( \sin^2 \vartheta \sqrt{1-e^2} + \frac{\cos^2 \vartheta}{\sqrt{1 - e^2}} \right) \equiv \frac{{\rm \pi} a^2 b^2}{4} \check{\mathcal{E}}_{\vartheta}^2. \end{gather}

Our final step is to impose that this area of the ellipse is the same area as the unit circle ($a b = 1$), which may be achieved by an appropriate choice of ${\mathrm \Delta} \psi$ and/or ${\mathrm \Delta} \alpha$. We thus have

(2.33)\begin{align} A & = {\rm \pi}^2 \left( \frac{ q {\mathrm \Delta} \psi {\mathrm \Delta} \alpha }{m}\right)^2 \iint \frac{f_{M,0}}{T_0} \nonumber\\ & \quad \times \left[ \hat{\mathcal{E}}_{\vartheta}^2 \omega_\alpha^2 \left( \frac{\omega_*^T}{\omega_\alpha} - 1 + F \right) \frac{{\mathrm \Delta} \psi}{{\mathrm \Delta} \alpha} + \check{\mathcal{E}}_{\vartheta}^2\omega_\psi^2 \left({-}1 + F \right) \frac{{\mathrm \Delta} \alpha}{{\mathrm \Delta} \psi} \right] \mathrm{d}\mu \,\mathrm{d}\mathcal{J} . \end{align}

Again, in the limit $\omega _\psi \rightarrow 0$, the Æ reduces to the previously found expression (Helander Reference Helander2020).Footnote ¹ The above equation is the central result of this paper, and may be interpreted as the most general result of the local available energy of trapped particles. To simplify further steps of the calculation, we shall specialize to the case in which the cross-section is circular in $(x,y)$, which implies $e = 0$. Note that this does not imply that the cross-section is also circular in $(\psi,\alpha )$-space: in these coordinates it is still an ellipse whose semi-major and semi-minor axes lie on lines of constant $\psi$ and $\alpha$. This furthermore has as a consequence that the various functions dependent on the angle of rotation and eccentricity simplify to $\mathcal {E}_\vartheta ^2 = 0$, $\hat {\mathcal {E}}_{\vartheta }^2=\check {\mathcal {E}}_{\vartheta }^2 = 1$, simplifying the calculations significantly. If the eccentricity is non-zero but small, the leading-order correction is of order $\mathcal {O}(e^2)$.

2.3 Making the integral dimensionless

We proceed by making the expression for Æ dimensionless. We first normalize the magnetic field to some reference strength $B_0$

(2.34)\begin{equation} B(\ell) = B_0\hat{B}. \end{equation}

Here, $B(\ell )$ is the magnetic-field strength as a function of the arc length $\ell$ along the flux tube. As in neoclassical transport theory (see Helander & Sigmar (Reference Helander and Sigmar2005), chapter 7), it is useful to perform a change of variables, $(\mu,\mathcal {J}) \mapsto (\lambda, z)$, with

(2.35a,b)\begin{equation} \lambda = \frac{\mu B_0}{\epsilon_0},\quad z = \frac{\epsilon_0}{T_0}. \end{equation}

In these variables, the second adiabatic invariant becomes

(2.36)\begin{equation} \mathcal{J} = \sqrt{2 m T_0} \sqrt{z} \int_{\{\ell_{b}\}} \sqrt{1 - \lambda \hat{B}}\,\mathrm{d} \ell, \end{equation}

if the electric field parallel to $\boldsymbol {B}$ vanishes. Here, we have introduced a shorthand for the integration domain, which is given by the interval in $\ell$ between two consecutive bounce points

(2.37)\begin{equation} \{\ell_{b}\}; \quad 1- \lambda \hat{B} (\ell_b) = 0. \end{equation}

The derivatives of $\mathcal {J}$ can be used to find bounce-averaged precession frequencies, as shown in Appendix B. We now define dimensionless precession frequencies, which are the dimensionless equivalents of the precession and drift frequencies in (2.25),

(2.38a)\begin{gather} \hat{\omega}_\alpha(\lambda,{\mathrm \Delta} \psi) \equiv \frac{q {\mathrm \Delta} \psi}{\epsilon_0}\omega_\alpha(\lambda), \end{gather}

(2.38b)\begin{gather}\hat{\omega}_\psi(\lambda,{\mathrm \Delta} \alpha) \equiv \frac{q{\mathrm \Delta} \alpha}{\epsilon_0} \omega_\psi(\lambda), \end{gather}

(2.38c)\begin{gather}\hat{\omega}_*^T(z, {\mathrm \Delta} \psi) \equiv \frac{q {\mathrm \Delta} \psi}{\epsilon_0} \omega_*^T(z). \end{gather}

Finally, one needs to account for the change in the volume element from $(\mu,\mathcal {J}) \mapsto (z,\lambda )$. The Jacobian of this transformation, which we denote by $G^{1/2}$, is

(2.39)\begin{align} G^{1/2}(z,\lambda) & = \frac{L}{\bar{B}} \frac{\sqrt{m} T_0^{3/2}}{\sqrt{2} } \sqrt{z} \int_{\{\ell_{b}\}} \frac{1}{\sqrt{1-\lambda\hat{B}}} \frac{ \mathrm{d}\ell}{L} \nonumber\\ & \equiv \frac{L}{\bar{B}} \frac{\sqrt{m }T_0^{3/2}}{\sqrt{2} } \sqrt{z} \hat{G}^{1/2}(\lambda), \end{align}

where $L$ denotes the length of the magnetic-field line. We are now in a position to express the Æ in terms of these new variables

(2.40)\begin{align} A & = \frac{1}{4 \sqrt{\rm \pi}} \frac{{\rm \pi} {\mathrm \Delta} \psi {\mathrm \Delta} \alpha L}{\bar{B}} n_0 T_0 \iint \sum_{\text{wells}(\lambda)} e^{{-}z} z^{5/2} \nonumber\\ & \quad \times\left[ \hat{\omega}_\alpha^2 \left( \frac{\hat{\omega}_*^T}{\hat{\omega}_\alpha} - 1 + \hat{F} \right) + \hat{\omega}_\psi^2 \left({-}1 + \hat{F} \right) \right] \hat{G}^{1/2} \,\mathrm{d}z \,\mathrm{d}\lambda . \end{align}

Here, we have introduced a summation over all bounce wells corresponding to a given value of $\lambda$, and we have introduced the nonlinear function

(2.41)\begin{equation} \hat{F} = \frac{\sqrt{(\hat{\omega}_*^T - \hat{\omega}_\alpha)^2 + \hat{\omega}_\psi^2 }}{\sqrt{\hat{\omega}_\alpha^2 + \hat{\omega}_\psi^2 } }. \end{equation}

The prefactor in front of the integral in (2.40) is readily interpreted. The factor ${\rm \pi} n_0 {\mathrm \Delta} \psi {\mathrm \Delta} \alpha L / \bar {B}$ is roughly the amount of particles residing in our flux tube, which we will call $N$, and hence the Æ is proportional to the total energy of these particles, $NT_0$. It is furthermore relatively straightforward to show that the integrand of (2.40) is always positive definite, as follows from the expression

(2.42)\begin{align} I(\hat{\omega}_\alpha, \hat{\omega}_\psi, \hat{\omega}_*^T) & = \hat{\omega}_\alpha^2 \left( \frac{\hat{\omega}_*^T}{\hat{\omega}_\alpha} - 1 + \frac{\sqrt{(\hat{\omega}_*^T - \hat{\omega}_\alpha)^2 + (\hat{\omega}_\psi)^2}}{\sqrt{(\hat{\omega}_\alpha)^2 + (\hat{\omega}_\psi)^2}} \right) \nonumber\\ & \quad +\hat{\omega}_\psi^2 \left({-}1 + \frac{\sqrt{(\hat{\omega}_*^T - \hat{\omega}_\alpha)^2 + (\hat{\omega}_\psi)^2}}{\sqrt{(\hat{\omega}_\alpha)^2 + (\hat{\omega}_\psi)^2}} \right). \end{align}

Since

(2.43)\begin{align} \frac{I(\hat{\omega}_\alpha,\hat{\omega}_\psi,\hat{\omega}_*^T)}{\hat{\omega}_\alpha^2 + \hat{\omega}_\psi^2} & = \frac{\hat{\omega}_*^T \hat{\omega}_\alpha}{\hat{\omega}_\alpha^2 + \hat{\omega}_\psi^2} + \sqrt{\frac{(\hat{\omega}_\alpha-\hat{\omega}_*^T)^2 + \hat{\omega}_\psi^2}{\hat{\omega}_\alpha^2 + \hat{\omega}_\psi^2}} - 1 \nonumber\\ & = \frac{\hat{\omega}_*^T \hat{\omega}_\alpha}{\hat{\omega}_\alpha^2 + \hat{\omega}_\psi^2} -1 + \sqrt{\left( \frac{\hat{\omega}_*^T \hat{\omega}_\alpha}{\hat{\omega}_\alpha^2 + \hat{\omega}_\psi^2} - 1 \right)^2 + \left( \frac{\hat{\omega}_*^T \hat{\omega}_\psi}{ \hat{\omega}_\alpha^2 + \hat{\omega}_\psi^2 } \right)^2 } \geq 0, \end{align}

the integrand in (2.40) must be positive definite, and so is therefore the Æ. It is also clear from this argument that even a small degree of non-omnigeneity in the magnetic-field geometry endows available energy to a plasma which otherwise has none, for the right-hand side of (2.43) is always positive if $\hat {\omega }_\psi \ne 0$. In other words, a Maxwellian can only be a ground state if the magnetic field is omnigenous (Helander Reference Helander2020).

Our final step is to find the fraction of energy that is available. The thermal energy of a plasma in a flux tube can readily be calculated to leading order in the perpendicular coordinates as

(2.44)\begin{equation} E_t = \frac{3}{2} \int \frac{n T}{B} \mathrm{d}\psi \,\mathrm{d}\alpha \,\mathrm{d}\ell \approx \frac{3}{2} n_0 T_0 \frac{{\rm \pi} {\mathrm \Delta} \psi {\mathrm \Delta} \alpha L}{B_0} \int \frac{1}{\hat{B}} \frac{\mathrm{d} \ell}{L}. \end{equation}

Hence, the fraction of energy that is available is equal to

(2.45)\begin{align} \frac{A}{E_t} & = \frac{2}{3} \left( \int \frac{4 \sqrt {\rm \pi}}{\hat{B}} \frac{ \mathrm{d} \ell}{L} \right)^{{-}1} \iint \sum_{\text{wells}(\lambda)} e^{{-}z} z^{5/2} \nonumber\\ & \quad \times\left[ \hat{\omega}_\alpha^2 \left( \frac{\hat{\omega}_*^T}{\hat{\omega}_\alpha} - 1 + \hat{F} \right) + \hat{\omega}_\psi^2 \left({-}1 + \hat{F} \right) \right] \hat{G}^{1/2} \,\mathrm{d}z \,\mathrm{d}\lambda. \end{align}

2.4 Relating available energy to turbulence

Our next step is to relate Æ to typical turbulence quantities, but there is basic difficulty having to do with the question of how these scale with the size of the system. There exist types of turbulence which depend on the size and shape of the domain in which it takes place, but we are mainly interested in turbulence for which this is not the case if the domain is large enough. We will refer to such turbulence as ‘local’. In a sufficiently resolved simulation, the correlation length of local turbulence will be smaller than the computational domain. Furthermore, if we expect Æ to encapsulate information about local turbulence, it should act as an extensive thermodynamic variable and thus scale linearly with the simulation volume. In expression (2.40), however, we see that under the replacement $({\mathrm \Delta} \psi, {\mathrm \Delta} \alpha ) \mapsto C({\mathrm \Delta} \psi, {\mathrm \Delta} \alpha )$ the Æ transforms as $A \mapsto C^4 A$ (as $\hat {\omega }_\psi$ and $\hat {\omega }_\alpha$ scale linearly with $C$). The underlying physical reason is that the Æ measures the maximum amount of energy that can be released by redistributing plasma over the entire domain, which scales as the volume of the domain (${\sim }C^2$) multiplied by the square of the variation of the distribution function over the domain (${\sim }C^2$), see Helander (Reference Helander2017). Local turbulence, however, is only able to redistribute plasma over some finite length scale comparable to the correlation length. The latter can be different in the two directions perpendicular to the magnetic field, and the appropriate choice for the domain size over which the system can redistribute particles are these length scales in the $\psi$- and $\alpha$-directions.Footnote ² Let us denote these length scales by ${\mathrm \Delta} \psi _{A}$ and ${\mathrm \Delta} \alpha _{A}$, respectively. For fixed ${\mathrm \Delta} \psi _{A}$ and ${\mathrm \Delta} \alpha _{A}$ we see that (2.45) is indeed independent of the domain size, and hence the total Æ acts as a thermodynamic variable.

To choose these length scales appropriately, we define radial and binormal coordinates in units of length by

(2.46a,b)\begin{equation} r = a \sqrt{\frac{\psi}{\psi_\mathrm{tot}}},\quad s = r_0 \alpha, \end{equation}

where $\psi _\mathrm {tot}$ is the total toroidal flux passing through the last closed flux surface, and $r_0 = a \sqrt {\psi _0/\psi _\mathrm {tot}}$. In terms of these coordinates, we require that the length scales are of the form

(2.47a,b)\begin{equation} {\mathrm \Delta} r_A = C_r \rho,\quad {\mathrm \Delta} s_A = C_s \rho, \end{equation}

where $\rho$ is the Larmor radius, and $C_r$ and $C_y$ are functions of order $\mathcal {O}(\rho ^0)$, which can depend on the equilibrium parameters such as the rotational transform $\iota$ or the magnetic shear $s$. We note that, by choosing proportionality to the gyroradius, the Æ becomes proportional to the expansion parameter $(\rho _*)^2 \equiv (\rho /L_\mathrm {ref})^2$, with $L_\mathrm {ref}$ being some equilibrium length scale. This is because $\omega _\alpha$ and $\omega _\psi$ are set by the equilibrium and exhibit an $L_\mathrm {ref}^{-1}$ dependence. The ratio of $C_r/C_s$ is a measure of anisotropy; if $C_r/C_s=1$ the correlation length is similar in the radial and binormal direction. Conversely, if there are large radial streamers present in a system one could reasonably expect $C_r/C_s \gg 1$. As such, it is unclear a priori what a proper length scaleis for both $C_r$ and $C_s$, as it will in general depend on the specific structure of the turbulence. For the rest of the investigation, we shall make use of the simplest possible ansatz for $C_r$ and $C_s$: we shall take them to be equal and independent of equilibrium parameters, i.e.

(2.48)\begin{equation} C_r = C_s = 1. \end{equation}

Finally, to facilitate a comparison with nonlinear gyrokinetic turbulence simulations, we define a dimensionless Æ

(2.49)\begin{equation} \hat{A} = \frac{A}{E_t \rho_*^2}. \end{equation}

Note that this choice for $\hat {A}$ differs from the one used in Mackenbach et al. (Reference Mackenbach, Proll and Helander2022), (up to a constant prefactor) by the factor $(\int \hat {B}^{-1} \,\mathrm {d} \ell / L)^{-1}$. In the results which are to be presented the inclusion or exclusion of this factor does not alter the results significantly. In devices with large variation in the magnetic-field strength, however, this factor can change $\hat {A}$ to a greater extent. Since the above expression for $\hat {A}$ is the fraction of energy that is available in the flux tube, we find that this definition is more natural. Finally, we denote some important dependencies of $\hat {A}$. Firstly, upon rescaling $(C_r,C_s) \mapsto C (C_r,C_s)$ by a factor $C$, the Æ scales as $\hat {A} \mapsto C^2 \hat {A}$. A rescaling of $C_r$ and $C_s$ may thus have a significant impact on the Æ. Furthermore, if one increases the magnitude of $C_s$ for fixed $C_r$, (2.43) implies that the Æ will increase as well. In general, one should thus expect results to be dependent on the choice of $C_r$ and $C_s$, and the choice of $C_r=C_s=1$ should be interpreted as the simplest possible model for these functions.