3.1. The set-up

Let us consider the simplest possible example in which to compute the available energy $A$ with no invariant coordinates ( $\boldsymbol{y} = \varnothing$ ): an electrostatic plasma consisting of a single species, similar to Chen (Reference Chen1966a ). The potential energy of such a system is

(3.1) \begin{equation} E_{\phi } = \int \mathrm{d}{\boldsymbol{r}}\, \frac {|\boldsymbol{\nabla }\phi |^{2}}{8\pi }, \end{equation}

where the electrostatic potential $\phi$ is related to density fluctuations $\delta n_e$ by Poisson’s equation, $\boldsymbol{\nabla} ^2 \phi \propto \delta n_e$ , or quasineutrality $\phi \propto \delta n_e$ . To make progress in solving (2.7), we assume that the initial distribution function $f(t=0)$ consists of a single waterbag: a distribution that can take on only two values; $\eta$ or zero. Under such conditions, the generally infinite number of Casimir invariants reduce to a single constraint for the simple reason that there is only a single level set (i.e. phase-volume conservation becomes equivalent to particle-number conservation; Kadomtsev & Pogutse Reference Kadomtsev and Pogutse1970; Ewart et al. Reference Ewart, Brown, Adkins and Schekochihin2022).

Consider for the sake of illustration the case depicted in figure 1. It should be evident that the Gardner ground state will collapse to the lowest velocity coordinate, whilst maintaining its total volume to satisfy incompressibility. In reality such a system will release its energy into electrostatic fluctuations, coupled to the density fluctuations, which forbid the system from reaching this state. Imposing that the ground state has some spatial density variation, $n(\boldsymbol{r})$ , will result in this ground state having a varying thermal width to accommodate the density variation (which will be symmetric about $v_{x}=0$ to minimise kinetic energy). Such considerations carry over to arbitrary dimensions: the ground state will collapse onto a sphere around $\boldsymbol{v} = 0$ (since the isocontours of particle energy are spheres), whose radius may vary spatially in order to support density variations. Expressing (2.5) in a convenient form (which in the case of the waterbag is either $\eta$ or zero, as captured by a Heaviside function), we thus have

(3.2) \begin{equation} F\left [\frac {1}{2}m|\boldsymbol{v}|^{2} + \kappa (\boldsymbol{r}) \right ] = \eta \varTheta [ a(\boldsymbol{r})^2 - |\boldsymbol{v}|^{2}], \end{equation}

where $a(\boldsymbol{r})$ is chosen such that $\int \mathrm{d}{\boldsymbol{v}} F = n(\boldsymbol{r})$ , i.e.

(3.3) \begin{equation} a(\boldsymbol{r}) = \frac {1}{\sqrt {\pi }}\left [\varGamma \left (1 + \frac {d}{2}\right ) \frac {n(\boldsymbol{r})}{\eta }\right ]^{1/d}, \end{equation}

where the gamma function $\varGamma (x) = (x-1)!$ has emerged to give an expression for the volumes of a ball in $d$ dimensions. The constrained ground state now readily gives the system’s kinetic energy

(3.4) \begin{equation} E_{\mathrm{K}} \geqslant \iint \mathrm{d}{\boldsymbol{r}}\mathrm{d}{\boldsymbol{v}}\frac {1}{2}m|\boldsymbol{v}|^{2}F = E_{\mathrm{G}} \left \langle \left [\frac {n(\boldsymbol{r})}{n_{0}} \right ]^{({d+2})/{d}} \right \rangle , \end{equation}

where $\langle \ldots \rangle \equiv \int \mathrm{d}{\boldsymbol{r}} \ldots / V$ is a volume average ( $V$ being the system’s volume), and the ground-state energy with an unconstrained density profile ( $\kappa =0$ ), $E_{\mathrm{G}}$ , is given as

(3.5) \begin{equation} E_{\mathrm{G}} = \frac {mV}{ \pi \eta ^{2/d} (2+d)\varGamma \left ({d}/{2}\right )} \left [ \varGamma \left ( 1 + \frac {d}{2} \right ) n_0 \right ]^{({d+2})/{d}}. \end{equation}

While (3.4) is specific to the single-waterbag model it serves to analytically illustrate the cost of a density perturbation: in order to conserve phase-space volume while enhancing density, more particles must be elevated to larger velocities, which costs more energy.

Figure 1.

A schematic showing the distinction between the Gardner state and an accessible ground state. In the leftmost panel we show a possible initial condition with energy available (the two-stream instability); its Gardner state (without a constraint on the density), $F_{\mathrm{G}}$ , is given in the middle panel, and it may be seen that it is spatially homogenous. However, in order to support an electric field containing the energy released by the instability, the ground state must have spatial variation, and an example of such a spatially varying ground state, $F$ , is given in the rightmost panel.

Figure 1. Long description

The image consists of three panels side by side, each depicting a different state of a plasma system. The leftmost panel shows an initial condition with energy available, specifically illustrating the two-stream instability. This panel features two horizontal shaded regions centered around positive and negative velocities. The middle panel represents the Gardner state without a constraint on the density. This state is spatially homogeneous, as indicated by the uniform shading across the entire panel. The rightmost panel shows an example of a spatially varying ground state, which is necessary to support an electric field containing the energy released by the instability. This panel features a wavy, non-uniform shading, indicating spatial variation. The x-axis in all panels ranges from negative pi to positive pi, and the y-axis ranges from negative v0 to positive v0, representing velocity. The panels collectively illustrate the mappin from ant initial state with available energy to a spatially homogeneous Gardner state and to a spatially varying ground state to support the released energy.

The Gardner bound is given as $A_{\mathrm{G}} = E_{\mathrm{init}} - E_{\mathrm{G}}$ , so that the available energy subject to a density perturbation is

(3.6) \begin{equation} A = A_{\mathrm{G}} - E_{\mathrm{G}}\langle [1 + \nu (\boldsymbol{r}) ]^{({d+2})/{d}} - 1 \rangle , \end{equation}

where we have set $n(\boldsymbol{r}) = n_0[1+\nu (\boldsymbol{r})]$ . Of course, all this energy must be stored in the energy of the fields $E_{\phi }$ which, for an electrostatic plasma, gives the condition

(3.7) \begin{equation} E_{\mathrm{\phi }}[\nu (\boldsymbol{r})] = \int \mathrm{d}{\boldsymbol{r}}\,\frac {|\boldsymbol{\nabla }\phi |^{2}}{8\pi } = A_{\mathrm{G}} - E_{\mathrm{G}}\langle [1 + \nu (\boldsymbol{r}) ]^{({d+2})/{d}} -1 \rangle , \end{equation}

where $\phi [\nu (\boldsymbol{r})]$ . Note that there are two additional constraints: one to conserve the total number of particles, $\int \mathrm{d}{\boldsymbol{r}} \nu = 0$ , and one to have positive number density, namely $\nu (\boldsymbol{r}) \geqslant -1, \ \text{for all }\boldsymbol{r}$ . In order to close these equations one needs to specify how $\phi$ and $\nu$ are related. For a fixed background of ions, the relationship is simply

(3.8) \begin{equation} -\lambda _{D,{\mathrm{ref}}}^2 \boldsymbol{\nabla} ^2 \hat {\phi } = \nu , \end{equation}

where we have introduced the Debye length $\lambda _{D,{\mathrm{s}}} = \sqrt {T_{\mathrm{s}}/4\pi \mathrm{e}^2 n_0}$ as a reference length scale, and the electric potential is normalised as $\hat {\phi } = \mathrm{e} \phi /T_{\mathrm{ref}}$ , where $T_{\mathrm{ref}}$ is some reference temperature (here an appropriate choice would be, e.g. $T_{\mathrm{ref}} \propto E_{\mathrm{G}}$ ). Equations (3.7) and (3.8) now constitute an implicit relation for the density constraint $\nu (\boldsymbol{x})$ . Before solving them, it is instructive to quickly arrive at an intuitive result. Letting the density (and thus potential) be non-zero only at some predefined wavenumber, $k_0$ , gives a field energy, via (3.7) and (3.8), that scales as $E_\phi \propto k_0^2 \phi _k^2 \propto \nu _{k}/k_0^2$ . Assuming the fluctuations are small, the available energy (3.6) becomes $A \approx A_{\mathrm{G}} - C \nu _k^2$ , where $C$ is some positive constant. By equating $A$ to $E_\phi$ we see that $E_\phi =A$ is possible in the limit $k_0 \rightarrow 0$ : purely electrostatic fluctuations maximise their released energy by producing the large-scale modes (i.e. long wavelengths can store large amounts of potential energy with relatively little density perturbation). For neutralising backgrounds, then, the Gardner bound is always attainable unless the lowest accessible wavenumber is non-zero (although whether this mode is dynamically activated is not answered by this framework).

However, the appearance of structures much larger than the Debye length is at odds with the assumption of a static background, as for sufficiently long wavelengths (i.e. longer than the Debye length) the electric field is shielded. With two species, Poisson’s equation becomes

(3.9) \begin{equation} -\lambda _{D,{\mathrm{bg}}}^2 \boldsymbol{\nabla} ^2 \hat {\phi } = \nu - \nu _{\mathrm{bg}}, \end{equation}

where the subscript ‘bg’ refers to a background species. Assume this species has a Boltzmann response $\nu _{\mathrm{bg}} = B ( \exp [\hat {\phi }] -1 )$ , having set $T_{\mathrm{ref}} = T_{\mathrm{bg}}$ , which is the temperature associated with the background. The potential is described by the Poisson–Boltzmann equation

(3.10) \begin{equation} - \lambda _{D,{\mathrm{bg}}}^2 \boldsymbol{\nabla} ^2 \hat {\phi } + B ({\rm e}^{\hat {\phi }} -1 ) = \nu , \end{equation}

where $B \in [0,1]$ is the parameter that controls if the background is fully neutralising ( $B=0$ ) or adheres to the Boltzmann distribution ( $B=1$ ), introduced so that both cases are captured in the calculation that follows.Footnote ³ The Boltzmann response introduces an additional effect that should be accounted for: if $B=1$ , the background species will inherit the perturbations of $\nu$ , which brings an associated energetic cost. Assuming this energetic cost scales like the waterbag solution, it may be set to $A_{\mathrm{bg}} = - E_{\mathrm{bg}} \langle [1 +\nu _{\mathrm{bg}}]^{({d+2})/{d}} - 1 \rangle ,$ where we have assumed that the background has no free energy in its initial state (i.e. the background does not contribute to $A_{\mathrm{G}}$ ).Footnote ⁴ The total available energy is now the sum of the two species, resulting in

(3.11) \begin{equation} A = A_{\mathrm{G}} - E_{\mathrm{G}}\big\langle [1 + \nu ]^{({d+2})/{d}} - 1 \big\rangle - E_{\mathrm{bg}} \big\langle [1 + \nu _{\mathrm{bg}} ]^{({d+2})/{d}} - 1 \big\rangle . \end{equation}

Our task is to maximise this available energy, subject to $A = E_\phi$ and a positive number density. In general, these equations do not permit analytical solutions, but they may be solved perturbatively when the initial condition is close to the ground state, and scalings may be derived when it is very far from it.

3.2. Limit of small ground-state density fluctuations

Let $|\nu |^2 \sim |\nu _{\mathrm{bg}}|^2 \sim A_{\mathrm{G}}/E_{\mathrm{G}}\ll 1$ and, in order to preserve the number of particles, we set $\int \mathrm{d}{\boldsymbol{r}} \nu = 0$ . Taylor expanding the Poisson–Boltzmann equation and taking its Fourier transform results in

(3.12) \begin{equation} \big(\lambda _{D,{\mathrm{bg}}}^2 k^2 + B \big) \hat {\phi }_{\boldsymbol{k}} = \nu _{\boldsymbol{k}}, \end{equation}

where the subscript $\boldsymbol{k}$ denotes the $\boldsymbol{k}$ th spectral component (e.g. $\nu (\boldsymbol{r}) = \int \mathrm{d}{\boldsymbol{k}} \nu _{\boldsymbol{k}} {\rm e}^{\mathrm{i} \boldsymbol{k} \boldsymbol{\cdot }\boldsymbol{r}}$ ). Expanding (3.11), invoking Plancherel’s theorem, and employing Poisson’s equation gives

(3.13) \begin{equation} A = A_{\mathrm{G}} - \frac {d+2}{d^2} \frac {(2 \pi )^d}{V} \int \mathrm{d}{\boldsymbol{k}}\, |\nu _{\boldsymbol{k}}|^2 \left (E_{\mathrm{G}} + E_{\mathrm{bg}} \frac {B^2}{\big(B+\lambda _{D,{\mathrm{bg}}}^2 k^2\big)^2}\right )\!. \end{equation}

Similarly manipulating the field energy yields

(3.14) \begin{equation} E_{\phi } = E_{\phi ,0} \frac {(2 \pi )^d}{V} \int \mathrm{d}{\boldsymbol{k}}\, \frac { \lambda _{D,{\mathrm{bg}}}^2 k^2 |\nu _{\boldsymbol{k}} | ^2}{\big(B + \lambda _{D,{\mathrm{bg}}}^2 k^2 \big)^2}, \end{equation}

where $E_{\phi ,0} = V n_0 T_{\mathrm{bg}}/2$ . Equating these expressions, $A=E_\phi$ , is equivalent to

(3.15) \begin{equation} \int \mathrm{d}{\boldsymbol{k}}\, \frac {(2\pi )^d|\nu _{\boldsymbol{k}}|^2}{V A_{\mathrm{G}}} \left [ E_{\phi ,0} \frac { \lambda _{D,{\mathrm{bg}}}^2 k^2 }{\left (B + \lambda _{D,{\mathrm{bg}}}^2 k^2 \right )^2} + E_{\mathrm{G}} \frac {d+2}{d^2} + E_{\mathrm{bg}} \frac {d+2}{d^2} \frac {B^2}{\big(B+ \lambda _{D,{\mathrm{bg}}}^2 k^2\big)^2} \right ] = 1, \end{equation}

and defining the integrand of (3.15) to be the distribution $\upsilon (\boldsymbol{k})$ with unit density automatically satisfies the equation. The field energy may now be expressed as

(3.16) \begin{equation} E_{\phi } = A_{\mathrm{G}}\! \int\! \mathrm{d}{\boldsymbol{k}}\, \upsilon (\boldsymbol{k}) \frac {E_{\phi ,0} \lambda _{D,{\mathrm{bg}}}^2 k^2 }{E_{\phi ,0} \lambda _{D,{\mathrm{bg}}}^2 k^2 + E_{\mathrm{G}} (({d+2})/{d^2}) \big(B+\lambda _{D,{\mathrm{bg}}}^2 k^2\big)^2 + B^2 E_{\mathrm{bg}} (({d+2})/{d^2})}, \end{equation}

and this quantity is maximised by letting $\upsilon (\boldsymbol{k})$ be a Dirac-delta distribution, non-zero at a wavenumber that maximises the fraction. This optimal wavenumber is given as

(3.17) \begin{equation} \lambda _{D,{\mathrm{bg}}}^2 k_0^2 = B \sqrt {1 + \frac {E_{\mathrm{bg}}}{E_{\mathrm{G}}}}, \end{equation}

and the available energy becomes

(3.18) \begin{align} A &= A_{\mathrm{G}} \left [1+ \frac {d+2}{d^2} \frac {E_{\mathrm{G}}\big(B+ \lambda _{D,{\mathrm{bg}}}^2 k_0^2\big)^2 + E_{\mathrm{bg}} B^2}{E_{\phi ,0} \lambda _{D,{\mathrm{bg}}}^2 k_0^2} \right ]^{-1} \nonumber\\[4pt] &=A_{\mathrm{G}} \left [ 1 + 2 B \frac {2+d}{d^2} \frac {E_{\mathrm{G}}}{E_{\phi ,0}} \left ( 1 + \sqrt {1+\frac {E_{\mathrm{bg}}}{E_{\mathrm{G}}} }\right ) \right ]^{-1}. \end{align}

An interesting result is found by considering the limiting case where $E_{\mathrm{bg}} \sim$ $E_{ \phi ,0} \gg E_{\mathrm{G}}$ , in which case the Gardner bound is approached asymptotically, $A \approx A_{\mathrm{G}}$ . This behaviour can be motivated physically: in order to circumvent the shielding response provided by the background species, it is beneficial to be at scales much smaller than the Debye length $k_0 \gg \lambda _{D,{\mathrm{bg}}}^{-1}$ . The background then behaves neutralising, and, as previously shown, the Gardner bound is approached if the potential has structure sizes that are much larger than the Debye length associated with the Gardner state (i.e. $k_0 \ll \lambda _{D,\mathrm{G}}^{-1} \propto \sqrt {E_{\mathrm{G}}}$ ). As $E_{\mathrm{bg}} \sim E_{ \phi ,0} \rightarrow \infty$ , the Debye length of the background and kinetic species asymptotically separate. By choosing a wavenumber that lies in between these two Debye lengths, the shielding response is evaded whilst simultaneously being at sufficiently long wavelength to approach the Gardner bound. Indeed, (3.17) satisfies these conditions and the Gardner bound is approached. Conversely, if $E_{\mathrm{bg}} \sim E_{ \phi ,0} \ll E_{\mathrm{G}}$ then maximising the energy release requires going to large length scales, but these are energetically disfavoured by the Debye shielding, so that to leading order $A \approx 0$ in this limit. Finally, we note that the scaling of the available energy with respect to the wavenumber, $k_0$ , is in line with Chen (Reference Chen1966a ) for $B=0$ (see Appendix A for comparison).

3.3. The constraint of non-negative number density: a second-order phase transition

As suggested by (3.15), for increasing initial energies, the variation in the density profile increases to a point where the constraint of a positive number density ( $\nu \geqslant -1$ ) becomes relevant (due to the factor $|\nu _{\boldsymbol{k}}|^2/A_{\mathrm{G}}$ ). In order to circumvent this constraint for as long as possible, it is thermodynamically beneficial to have a one-dimensional profile of density meaning the isotropic ground state must spontaneously break spatial symmetry (but remains isotropic in velocity space). To see why this is the case, consider an example in dimensionality $d=2$ of two profiles with the same energetic cost to the ground-state energy: $n_1(x,y)/n_0 =1 + \cos (k_x x)$ and $n_2(x,y)/n_0 = 1 + \cos (k_x x)/\sqrt {2} + \cos (k_y y)/\sqrt {2}$ (where $k_x = k_y = k_0$ ). It is clear that $n_1/n_0 \geqslant 0$ , whereas $\min [n_2/n_0] = 1- \sqrt {2} \lt 0$ , thus breaking positive definiteness of the density, so that this solution is prohibited. More generally, when a discrete set of modes $\{\nu _i \}$ are present, the density is bounded from below by $1 - \sum |\nu _i|$ , which, for fixed $\sum \nu _i^2$ , is maximised by having only one non-zero $\nu _i$ , thus favouring one-dimensional profiles. Let us then assume that only a single mode is present, and we wish to estimate the Gardner energy at which the constraint of positive definiteness starts limiting the ground-state energy, which we denote as $A_{\mathrm{ph}}$ . Setting $\nu = \nu _0 \cos (k_0 x)$ , and finding the critical initial energy for which $\nu _0^2 = 1$ (still employing the small- $\nu$ equation) gives

(3.19) \begin{align} A_{\mathrm{ph}} & \approx \frac{d+2}{2d^2} E_{\mathrm{G}} + \frac{ (({d+2})/{d^2}) B^2 E_{\mathrm{bg}}+ \lambda_{D,{\mathrm{bg}}}^2 k_0^2 E_{\phi,0}}{2\big(B + \lambda_{D,{\mathrm{bg}}}^2k_0^2 \big)^2} \nonumber \\ & \approx \frac{d+2}{2d^2}E_{\mathrm{G}} + \frac{(({d+2})/{d^2}) B E_{\mathrm{bg}} + E_{\phi,0} \sqrt{1+({d^2}/({d+2}))({E_{\mathrm{bg}}}/{E_{\mathrm{G}}}) }}{2B( 1 + \sqrt{1+({d^2}/({d+2}))({E_{\mathrm{bg}}}/{E_{\mathrm{G}}}) } )^2}. \end{align}

It is important to stress that this is an estimate – except for $d=2$ where the result is exact – since, for $\nu _0 \sim 1$ , the nonlinear terms in (3.11) will affect the mode distribution.

When the critical energy is reached, new modes are activated in order to maintain positive definiteness. It may be postulated that near this energy the ground state exhibits a phase transition,Footnote ⁵ where the mode amplitude scales as

(3.20) \begin{equation} \nu _{\boldsymbol{k}} \approx \mathfrak{F}_0(\boldsymbol{k})A_{\mathrm{G}} - \mathfrak{F}(\boldsymbol{k}) \varTheta (A_{\mathrm{G}}-A_{\mathrm{ph}}) (A_{\mathrm{G}}-A_{\mathrm{ph}})^{\alpha }, \end{equation}

where $\mathfrak{F}_0(\boldsymbol{k})$ contains the mode distribution at $A_{\mathrm{G}} = A_{\mathrm{ph}}$ and $\mathfrak{F}(\boldsymbol{k})$ is some function describing the relative amplitude of the newly activated modes. Substituting this expression into (3.14) and imposing that it is increasing in $A_{\mathrm{G}}$ yields $\alpha \geqslant 1$ . In Helander & Mackenbach (Reference Helander and Mackenbach2024) it is shown that the Helmholtz free energy $H$ can be made equivalent to the available energy, at least in a certain asymptotic regime. Positing that $A \sim H$ here too, and that $A_{\mathrm{G}} \sim T$ , we conclude that $\partial _T^{\lceil \alpha \rceil } H$ is discontinuous and is thus a $\lceil \alpha \rceil$ -order phase transition (where $\lceil \alpha \rceil$ rounds $\alpha$ up to the nearest integer). Though we have postulated and argued for the general case here as a consequence of avoiding a negative number density and maximising the minimal number density (processes that are likely unaffected by additional details), one may calculate this phase-transition behaviour exactly in certain limiting cases, as we shall see in the proceeding paragraph.

The calculation of the ground state and mode spectrum can be carried out exactly for $d=2$ , in the case of a system of finite size with a neutralising background, where $A_{\mathrm{ph}} = E_{\mathrm{G}}/2 + E_{\phi ,0}/2\lambda _{D,{\mathrm{bg}}}^2k_0^2$ . For ease of notation, we set $\nu (\boldsymbol{x}) = \sum _{\boldsymbol{k}} \nu _{\boldsymbol{k}} e^{\mathrm{i} \boldsymbol{k} \boldsymbol{\cdot }\boldsymbol{x}}$ , so that $A=E_\phi$ yields

(3.21) \begin{equation} \sum _{\boldsymbol{k}} \frac { E_{\mathrm{G}} + \big({E_{\phi ,0}}/{\lambda _{D,{\mathrm{bg}}}^2 k^2}\big) }{A_{\mathrm{G}}} |\nu _{\boldsymbol{k}}|^2 = 1. \end{equation}

Analogously to (3.15) and the discussion below it, this is achieved in (3.21) by defining the summand to be $\mathcal{U}_{\boldsymbol{k}}$ whose sum evaluates to unity. The field energy is then

(3.22) \begin{equation} E_\phi = A_{\mathrm{G}} \sum _{\boldsymbol{k}} \frac {\mathcal{U}_{\boldsymbol{k}}}{E_{\mathrm{G}} + \big({E_{\phi ,0}}/{\lambda _{D,{\mathrm{bg}}}^2 k^2}\big)} \frac {E_{\phi ,0}}{\lambda _{D,{\mathrm{bg}}}^2 k_0^2}, \end{equation}

maximised by setting $\mathcal{U}_{\boldsymbol{k}} = 1/2$ for $\boldsymbol{k} = (k_x,k_y) = (\pm k_0,0)$ (or rotations thereof), below the critical energy. At the critical energy we have $\nu _{\boldsymbol{k}} = 1/2$ for $\boldsymbol{k} = (\pm k_0,0)$ , so that just above it new modes are activated as

(3.23) \begin{equation} \nu _{\boldsymbol{k}} \approx \begin{cases} \dfrac{1}{2} + \tilde {C}_1, & \boldsymbol{k} = (\pm k_0,0), \\[10pt] \tilde {C}_2, & \boldsymbol{k} = (\pm 2 k_0,0), \\[3pt] 0, & \text{otherwise}, \end{cases} \end{equation}

where $\tilde {C}_1 = C_1 (A_{\mathrm{G}} - A_{\mathrm{ph}})^{\alpha _1}$ and $\tilde {C}_2 = C_2 (A_{\mathrm{G}} - A_{\mathrm{ph}})^{\alpha _2}$ . Positive definiteness of the density then requires $\tilde {C}_1 \leqslant \tilde {C}_2$ , and the requirement that summing $\mathcal{U}_{\boldsymbol{k}}$ gives unity necessitates

(3.24) \begin{equation} A_{\mathrm{G}} - A_{\mathrm{ph}} = 4 A_{\mathrm{ph}} \tilde {C}_1 + 4 A_{\mathrm{ph}} \tilde {C}_1^2 + A_{\mathrm{ph}} \tilde {C}_2^2 + \frac {3 E_{\mathrm{G}}}{2} \tilde {C}_2^2, \end{equation}

and the field energy can similarly be calculated:

(3.25) \begin{align} E_{\phi } &= \frac {E_{\phi ,0}}{2\lambda _{D,{\mathrm{bg}}}^2 k_0^2} \big( 1 + 4 \tilde {C}_1 + 4\tilde {C}_1^2 + \tilde {C}_2^2 \big)\nonumber \\ &= \frac {E_{\phi ,0}}{2\lambda _{D,{\mathrm{bg}}}^2 k_0^2} \left ( \frac {A_{\mathrm{G}} }{A_{\mathrm{ph}}} - \frac {3E_{\mathrm{G}}}{2 A_{\mathrm{ph}}} \tilde {C}_2^2 \right )\!. \end{align}

Balancing of (3.24) can be achieved in the following four distinct ways.

1. (i) If $2\alpha _2 \gt 1$ , implying that $\alpha _1 = 1$ and $C_1 =1/(4 A_{\mathrm{ph}})$ . Positive definiteness requires $\alpha _1 \geqslant \alpha _2$ , so that $\alpha _2 \in (1/2,1]$ . The term decreasing the field energy goes as ${\sim} C_2^2 (A_{\mathrm{G}} - A_{\mathrm{ph}})^{2\alpha _2}$ , where setting $\alpha _2 = 1$ yields the slowest decrease. Since $C_2 \geqslant 1/4A_{\mathrm{ph}}$ , the available energy is maximal when $C_2$ touches its lower bound yielding

   (3.26) \begin{equation} E_{\phi } = \frac {E_{\phi ,0}}{2\lambda _{D,{\mathrm{bg}}}^2 k_0^2} \left [ \frac {A_{\mathrm{G}} }{A_{\mathrm{ph}}} - \frac {3E_{\mathrm{G}}( A_{\mathrm{G}} - A_{\mathrm{ph}} )^2}{32 A_{\mathrm{ph}}^3} \right ]\!. \end{equation}

2. (ii) If $2 \alpha _2 \lt 1$ , necessitating $\alpha _1=2\alpha _2$ . We conclude immediately that the field energy decreases more rapidly with increasing $A_{\mathrm{G}}$ in this case than in the previous, and we thus disregard it.

3. (iii) If $2 \alpha _2 = 1$ and $\alpha _1 \gt 1$ , we conclude again that penalty is subquadratic and this solution may be set aside.

4. (iv) If $2\alpha _2=\alpha _1=1$ , the solution can similarly be dismissed.

The available energy is hence given as

(3.27) \begin{equation} A \approx \frac {E_{\phi ,0}}{2\lambda _{D,{\mathrm{bg}}}^2 k_0^2} \left [ \frac {A_{\mathrm{G}}}{A_{\mathrm{ph}}} - \frac {3E_{\mathrm{G}}( A_{\mathrm{G}} - A_{\mathrm{ph}} )^2}{32 A_{\mathrm{ph}}^3} \varTheta (A_{\mathrm{G}} - A_{\mathrm{ph}}) \right ]\!, \end{equation}

and we conclude that the ground state exhibits a second-order phase transition. A plot of the density profile for different values of $A_{\mathrm{G}}$ is given in figure 2, where it can be seen that the adjustment of the modes keep the density profile positive even beyond the critical energy. However, a second transition point arises at $A_{\mathrm{G}} = 5A_{\mathrm{ph}}/3$ , since $n^{\prime\prime}(x)=0$ at the minimum, and one has to perform a similar exercise by including a new mode tailored to keep the density profile positive for higher values of $A_{\mathrm{G}}$ . Ultimately, this is indicating a form of diminishing returns: a large $A_{\mathrm{G}}$ (free energy) gives a greater release of energy, but an increasing amount of that additional energy is spent on density perturbations than on the electric field. Physically, this phase transition is difficult to realise in any numerical simulation as the available energy can generally only be achieved instantaneously (as we have not imposed stationarity of the ground state), and this particular transition corresponds to ever more energy being channelled into the density fluctuations. Qualitatively, however, it agrees with the intuition that large amplitude turbulence tends to ‘cellularise’.

Figure 2.

Density profile for the Fourier modes given in (3.23) for various values of the Gardner bound, where $A_{\mathrm{ph}} = \lambda _{D,{\mathrm{bg}}} = k_0 = 1$ . It may be seen that the density profile narrows beyond the critical energy in order to maintain positive definiteness. If the Gardner free energy becomes too large, negative densities arise since $n^{\prime\prime}(\pi )\lt 0$ , and these regions are shaded red.

Figure 2. Long description

The line graph displays the density profile for the Fourier modes for various values of the Gardner bound. The x-axis ranges from negative pi to positive pi, and the y-axis represents the normalized density n(x)/n0, ranging from 0 to 4. Multiple lines, each representing different values of the Gardner bound, show how the density profile changes. The lines are color-coded according to the Gardner bound values, with a color bar on the right indicating the range from 0.5 to 3.0. The graph shows that the density profile narrows beyond the critical energy to maintain positive definiteness. Regions where the Gardner free energy becomes too large, leading to negative densities, are shaded red.

As one keeps increasing $A_{\mathrm{G}}$ , one enters another asymptotic regime where the density distribution can be approximated by a Dirac-delta peak. This case is treated in detail in Appendix B, where different background responses are considered for one-dimensional profiles, and all scalings found (including the linear regime) are numerically verified.

3.4. Comparison to simulations

We assess the validity of the analytical results in the case where $d=1$ , for which a set of $32$ nonlinear PIC simulations have been constructed, performed with the osiris-code (Fonseca et al. Reference Fonseca2002). The initial condition of the simulation is a simple two-stream instability consisting of ions, with beams modelled by waterbags, as depicted in figure 1. Furthermore, 16 of these simulations were performed with a neutralising background (i.e. $B=0$ in the notation of (3.10)) and 16 have been performed with a kinetic background consisting of heavy electrons ( $m_e/m_i=1/100$ ), but still light enough that we may take them to be approximately adiabatic. In each set of 16 simulations the width of the beam (in the notation of figure 1, $\varDelta$ ) is varied so that the background temperature remains fixed, whereas the equivalent ‘temperature’ of the ground state, ${\sim} \varDelta ^2$ , changes, altering both $E_{\mathrm{G}}/E_{\phi ,0}$ and $E_{\mathrm{G}}/E_{\mathrm{bg}}$ .Footnote ⁶ A depiction of one such simulation with a kinetic background species is given in figure 3, where one can see the initial condition consisting of two streams, where at a later time Bernstein–Greene–Kruskal (BGK) modes are present (Bernstein, Greene & Kruskal Reference Bernstein, Greene and Kruskal1957).

Figure 3.

A PIC simulation of the evolution of an ion-two-stream instability with a Maxwellian background of electrons for two times $t=0\, \omega _{\mathrm{pe}}^{-1}$ (left-hand column) and $t=400\, \omega _{\mathrm{pe}}^{-1}$ where the electric-field energy has attained ${\sim}80 \,\%$ of its maximum value in this simulation (right-hand column). While the density fluctuations in the ions naturally occur on $\lambda _{D}$ scales, and the electron response is approximately Boltzmann-like, the ion response is significantly different from the ground state (cf. figure 1).

Figure 3. Long description

A PIC simulation of ion-two-stream instability with Maxwellian electrons at two times. The left column shows the initial state at time t equals 0 inverse omega pe, with ions forming horizontal bands and electrons displaying a uniform distribution. The right column shows the state at time t equals 400 inverse omega pe, where the turbulence is fully developed. Ions exhibit complex, wave-like patterns (BGK-holes) with labeled wavelength lambda D, while electrons show a more diffuse, wavy structure. The color scale on the right indicates the normalized distribution function f over f max, ranging from 0 to 1.

The proportionality constant $E_{\mathrm{bg}}$ is set by the assumed energetic cost of density perturbations to the background $A_{\mathrm{bg}} = - E_{\mathrm{bg}} \langle [1 + \nu _{\mathrm{bg}}]^{3} -1 \rangle$ and, as will be shown in § 4.1, this relationship is exactly true for a Maxwellian background in the small fluctuation limit, where the prefactor may be calculated analytically. Assuming the prefactor from the linear regime holds true generally, we set

(3.28) \begin{equation} E_{\mathrm{bg}} = \frac {n_0 T_{\mathrm{bg}} V}{6}. \end{equation}

Similarly for the beams, using the relationship $n_0 = 2 \Delta \eta$ gives the initial energy $E_i = m v_0^2V/2 + m \varDelta ^2 V/24$ and ground-state energy $E_{\mathrm{G}} = m n_0 \varDelta ^2 V/6$ . The Gardner bound becomes

(3.29) \begin{equation} A_{\mathrm{G}} = 3 E_{\mathrm{G}} \left ( \frac {v_0^2}{\varDelta ^2} - \frac {1}{4} \right ), \end{equation}

where it is required that $v_0/\varDelta \geqslant 1/2$ . Let us reiterate that the prefactor for the field energy is defined to be $E_{\phi ,0} = n_0 T_{\mathrm{bg}} V/2$ . Given $m$ , $\varDelta$ , $v_0$ , $T_{\mathrm{bg}}$ and the above relationships, one can calculate $A$ . This value is compared with the maximal field energy in the simulation, denoted as $E_{\phi ,{\mathrm{max}}}$ , seeing what fraction of the total available energy is actually released. Let us note that, for the neutralising case with finite box size, using the above relations, the available energy in the linear regime is given as

(3.30) \begin{equation} A_{\mathrm{linear}} = \frac {A_{\mathrm{G}}}{1 + \big({d_e^2 k_{\mathrm{min}}^2 \varDelta ^2}/{c^2}\big)}, \end{equation}

where $c$ is the speed of light and $d_e$ is the electron skin depth (which is larger than the Debye length by a factor $c/v_T$ , with $v_T$ being the thermal velocity). Since $\varDelta /c$ is chosen small in all simulations, whereas the domain size is chosen to be one skin depth, the correction factor in the denominator is in all cases small and the available energy is close to the Gardner bound.

The main result is displayed in figure 4, where one sees how much of the available energy is released (including comparisons to the Gardner bound, the linearised result and the fully nonlinear result; see Appendix D for details on how the latter is calculated numerically). As more and more physics is included (going from red to green), one can see the estimates get increasingly better, as expected. However, it can also be seen that the available energy overestimates the total energy release by a factor of roughly $10$ , and this discrepancy is caused by two complications. Firstly, the dominant wavenumber, ${\sim} \lambda _{D,{\mathrm{bg}}} k_0$ , that the simulation settles in is in general not the same as given in (3.17) (e.g. in the neutralising case there are modes present with wavelengths smaller than the box size). Secondly, from (2.5) and (2.6), a necessary condition for a density-constrained ground state is that

(3.31) \begin{equation} \langle F^{\prime} [\epsilon +\kappa (\boldsymbol{r}),\boldsymbol{y}] \rangle \leqslant 0, \end{equation}

where a prime denotes a derivative with respect to the first argument. As such, we can investigate the velocity-space structure of the volume-averaged distribution function in the simulation, where we additionally average over time in the ‘saturated’ (i.e. turbulent) state to extract an average structure (further reducing the noise typical in PIC simulations). The result is displayed in figure 5, where the time-and-volume-averaged distribution function is plotted, normalised by its maximal value, denoted as $\langle \overline {f}(x)\rangle$ . For $\varDelta /v_0 = 1.4$ , the distribution function clearly breaks the condition given in (3.31), indicating that more energy may be released by restacking, which is likely inhibited due to the presence of BGK modes (seen in figure 3). In a more nonlinear regime, $\varDelta /v_0=0.1$ , the distribution function gets closer to a density-constrained ground state, explaining the improved prediction. We note that plotting at $\langle f(t=t_{\mathrm{max}},x) \rangle$ , where $E_{\phi }(t_{\mathrm{max}}) = E_{\phi , \mathrm{max}}$ , shows similar trends, though somewhat obfuscated by the higher levels of noise.

Figure 4.

The fraction of the maximal field energy measured in the simulation, over the available energy. Three available energies are added for the comparison: the Gardner value $A_{\mathrm{G}}$ , the linear estimate $A_{\mathrm{linear}}$ (see (3.18)), where for the neutralising case, $B=0$ and $\lambda _{D,{\mathrm{bg}}} k_0 = \lambda _{D,{\mathrm{bg}}} k_{\mathrm{min}}$ ), and a full numerical calculation with no additional assumptions denoted by $A$ .

Figure 4. Long description

The scatter plot compares the fraction of the maximal field energy measured in the simulation over the available energy using three different methods. The x axis represents the variable delta over v naught, ranging from 0 to 1.5. The y axis represents the percentage, ranging from 0 to 30 percent. The graph is divided into two sections: Neutralising and Kinetic. Three data series are plotted: E phi max over A, E phi max over A linear, and E phi max over A G. Each data series is represented by different markers: green circles for E phi max over A, orange crosses for E phi max over A linear, and red plus signs for E phi max over A G. In the Neutralising section, all three data series show a decreasing trend as delta over v naught increases. In the Kinetic section, the green circles and orange crosses show a peak around delta over v naught equals 0.5 and then decrease, while the red plus signs show a decreasing trend throughout.

Figure 5.

Velocity-space structure of the waterbag distribution function averaged over the simulation volume and the time in the ‘saturated’ state, normalised by its maximal value. Blue lines denote simulations with a neutralising background, whereas the red lines have a kinetic background species. The initial distribution at $t=0$ , normalised by its maximal value, is included in grey. The top row has $\varDelta /v_0=1.4$ , whereas the bottom row has $\varDelta /v_0=0.1$ . The left column has a neutralising background whereas the right has a kinetic background.

Figure 5. Long description

The image contains four line graphs comparing the velocity-space structure of the waterbag distribution function. The graphs are arranged in a 2x2 grid. The top row represents a higher value of Δ/v0, while the bottom row represents a lower value. The left column shows data for a neutralising background, and the right column shows data for a kinetic background. Each graph plots the normalised distribution function against velocity (v/c). Blue lines denote simulations with a neutralising background, and red lines denote simulations with a kinetic background. Grey lines represent the initial distribution function. The graphs illustrate how the distribution function evolves over time in the saturated state, highlighting differences between neutralising and kinetic backgrounds, and, crucially, showing that the distribution function do not always fully relax in an averaged sense.

Finally, scalings of how much energy goes into perturbing the density of the ground state are given in (3.11), where we remind ourselves that the scaling for the background was posited to be of the same form as the waterbag species, employing the result from § 4.1 to set $E_{\mathrm{bg}} = n_0 T_{\mathrm{bg}} V/6$ . It is thus worthwhile to plot the thermal energy as a function of the density-fluctuation level, seeing if the predictions hold. For the waterbag species, we have that the thermal energy increases as $E_{\mathrm{T,wb}} - E_{\mathrm{G}} = E_{\mathrm{G}}\langle (1 + \nu )^3 - 1 \rangle$ , whereas for the background, we have $E_{\mathrm{T,bg}} - E_{\mathrm{T,bg,0}} = E_{\mathrm{bg}} \langle (1 + \nu )^3 - 1 \rangle$ , having defined $E_{\mathrm{T,bg}}(t=0)=E_{\mathrm{T,bg,0}}$ and $E_{\mathrm{T}}$ as the thermal energy. To investigate how closely these scalings are followed, we have scattered the thermal energies against the density-fluctuation level for all the simulations, displayed in figure 6. It can be seen that the background species adheres to the proposed scaling fairly closely, giving confidence in the posited scaling, with some deviations arising due to non-adiabatic effects. The same, however, can not be said for the waterbag distributions. At the start of the simulation one can observe that the thermal energy decreases with increasing fluctuation level. This is to be expected from energy conservation considerations: as energy is being released by the instability, it must be transferred to the field via the density fluctuations, explaining the observed trend. However, as one approaches the ground state in which no more energy can be released, density fluctuations instead cost the electric-field energy. Such an uptick in thermal energy with increasing fluctuation level is not seen for the waterbags in figure 6, giving further evidence that this distribution is unable to relax to a regime where it behaves similarly to a ground state. This suggests that a linear description of the plasma dynamics may be sufficient: we perform such an analysis in Appendix C where it is indeed shown that the linear dynamics explain the structures seen in figure 6 fairly well.

Figure 6.

Scatter of the deviation in thermal energy against the density-fluctuation level for the waterbag distribution in simulations with a neutralising background (left), with a kinetic background that is approximately Boltzmannian (centre) and for the background species itself (right). The expected scaling for a distribution that is in its ground state is included as a black dashed line in all plots, and the points are furthermore coloured according to their simulation-time value (normalised by the total simulated time in each simulation).

Figure 6. Long description

Three scatter plots compare the deviation in thermal energy against the density-fluctuation level for different simulation scenarios. The left plot shows data for simulations with a neutralising background, the center plot for simulations with a kinetic background that is approximately Boltzmannian, and the right plot for the background species itself. Each plot includes a black dashed line representing the expected scaling for a distribution in its ground state. The points in the plots are colored according to their simulation-time value, normalized by the total simulated time in each simulation. The x-axis for all plots represents the density-fluctuation level, while the y-axis represents the deviation in thermal energy. The plots illustrate how the energy deviation varies with the density-fluctuation level under different background conditions.

All in all, it is clear that by including more physics the available energy bound gets sharper and closer to what is found in simulations. However, fairly large mismatches remain, which we partly attribute to BGK modes making it difficult for the waterbags to relax to a state akin to the ground state. Since the behaviour, and stability, of such modes has special dependence on the dimensionality of the system (see, e.g. Ng & Bhattacharjee Reference Ng and Bhattacharjee2005), it may be interesting to do a similar exercise in higher dimensionality, though we make no such attempt here.

4.1. The ground state and available energy

We turn our attention to the ground-state equation, which can be found by equating the Casimirs for each $\boldsymbol{y}$ . The volume of phase space in which $w$ exceeds $\epsilon +\kappa$ at fixed $\boldsymbol{y}$ is given by

(4.10) \begin{equation} \varOmega (w,\boldsymbol{y}) = \int \mathrm{d}{\boldsymbol{z}}\, \varTheta [ w - \epsilon - \kappa ] \sqrt {g}. \end{equation}

The ground-state equation (2.7) then becomes

(4.11) \begin{equation} \int \mathrm{d}{\boldsymbol{z}}\, \varTheta [ f(\boldsymbol{x}) - F(w,\boldsymbol{y}) ] \sqrt {g} = \varOmega (w,\boldsymbol{y}), \end{equation}

and we proceed to solve this equation in orders of $\varepsilon$ . Since $f_0 = F_0$ , to leading order we have

(4.12) \begin{equation} \varTheta [ f_0(\boldsymbol{x}) - F_0(w,\boldsymbol{y}) ] = \varTheta [w -\epsilon ], \end{equation}

and taking the first and second derivative of this equation with respect to $w$ , one finds that

(4.13) \begin{equation} \delta [f_0(\boldsymbol{x}) - F_0(w,\boldsymbol{y})] = - \frac {\delta [w -\epsilon ]}{ F_0^{\prime} (w,\boldsymbol{y})} \end{equation}

and

(4.14) \begin{equation} \delta ^{\prime} [f_0(\boldsymbol{x}) - F_0(w,\boldsymbol{y})] = \frac {1}{F_0^{\prime} (w,\boldsymbol{y})} \frac {\partial }{\partial w} \left ( \frac {\delta [w -\epsilon ]}{F_0^{\prime} (w,\boldsymbol{y})} \right )\!, \end{equation}

respectively.

The next order, $\mathcal{O}(\varepsilon )$ , gives an equation involving derivatives of the Heaviside function, which is the Dirac-delta distribution $\delta [x]$ , giving rise to

(4.15) \begin{equation} \int \mathrm{d}{\boldsymbol{z}}\, [ f_1(\boldsymbol{x}) - F_1(w,\boldsymbol{y}) ]\, \delta [f_0(\boldsymbol{x}) - F_0(w,\boldsymbol{y})] \sqrt {g} = - \int \mathrm{d}{\boldsymbol{z}}\, \kappa _1(\boldsymbol{r}) \delta [w -\epsilon ] \sqrt {g} , \end{equation}

and using (4.13) this may be written as

(4.16) \begin{equation} \int \mathrm{d}{\boldsymbol{z}} ( f_1 - F_1) \delta [w -\epsilon ] \sqrt {g} = F_0^{\prime} (w,\boldsymbol{y}) \int \mathrm{d}{\boldsymbol{z}}\, \kappa _1 \delta [w -\epsilon ] \sqrt {g}. \end{equation}

Filling this into the leading-order result of the available energy (4.7), writing $\epsilon = w \delta [w-\epsilon ]$ , gives

(4.17) \begin{align} A &\approx \varepsilon \int \mathrm{d}{\boldsymbol{y}}\, \mathrm{d}{\boldsymbol{z}}\, \mathrm{d}{w} [ f_1(\boldsymbol{x}) - F_1(w,\boldsymbol{y}) - \kappa _1(\boldsymbol{r}) F_0^{\prime} (w,\boldsymbol{y}) ] w \delta [w-\epsilon ] \sqrt {g}\nonumber \\ &\approx 0. \end{align}

The available energy to order $\mathcal{O}(\varepsilon )$ thus vanishes and one is required to go to second order. It shall prove useful to rewrite (4.16) in terms of $F_1$ alone, using (4.9). One finds the following integral equation for $F_1$ ,

(4.18) \begin{equation} \int \mathrm{d}{\boldsymbol{z}} \left (\frac {F_1}{F_0^{\prime} } - \frac {\langle F_1 \rangle _{\boldsymbol{v}}}{\langle F_0^{\prime} \rangle _{\boldsymbol{v}}} \right ) \delta [w-\epsilon ] \sqrt {g} = \int \mathrm{d}{\boldsymbol{z}} \left (\frac {f_1}{F_0^{\prime} } - \frac {n_1}{\langle F_0^{\prime} \rangle _{\boldsymbol{v}}} \right ) \delta [w-\epsilon ] \sqrt {g}, \end{equation}

having defined $\langle \ldots \rangle _{\boldsymbol{v}} = \int \mathrm{d}{\boldsymbol{v}} \ldots$ .

To second order $\mathcal{O}(\varepsilon ^2)$ , the ground-state equation becomes

(4.19) \begin{align} &\!\int\! \mathrm{d}{\boldsymbol{z}}\, \{ [f_1(\boldsymbol{x}) - F_1(w,\boldsymbol{y})]^2 \delta ^{\prime} [f_0(\boldsymbol{x}) - F_0(w,\boldsymbol{y})] -F_2(w,\boldsymbol{y}) \delta [f_0(\boldsymbol{x}) - F_0(w,\boldsymbol{y})] \} \sqrt {g} \nonumber\\[4pt] &\quad = \int \mathrm{d}{\boldsymbol{z}}\, ( \kappa _1(\boldsymbol{r})^2 \delta ^{\prime} [ w - \epsilon ] - \kappa _2(\boldsymbol{r}) \delta [w - \epsilon ] ) \sqrt {g} . \end{align}

Using (4.13) and (4.14) we find that

(4.20) \begin{align} &\int \mathrm{d}{\boldsymbol{z}} \left \{ (f_1 - F_1)^2\frac {\partial }{\partial w} \left ( \frac {\delta [w -\epsilon ]}{F_0^{\prime} } \right ) - \kappa _1^2 \frac {\partial F_0^{\prime} \delta [ w - \epsilon ]}{\partial w}\right \} \sqrt {g} \nonumber\\[4pt] &\quad = -\int \mathrm{d}{\boldsymbol{z}} [ F_2 + \kappa _1^2 F_0^{\prime\prime} + \kappa _2 F_0^{\prime} ]\, \delta [w - \epsilon ] \sqrt {g} . \end{align}

Substituting this into the available energy, it may be written as

(4.21) \begin{align} A &\approx -\frac {\varepsilon ^2}{2} \int \mathrm{d}{\boldsymbol{y}}\, \mathrm{d}{\boldsymbol{z}}\, \mathrm{d}{w} [ F_2(w,\boldsymbol{y}) + \kappa _1^2 F_0^{\prime\prime}(w,\boldsymbol{y}) + \kappa _2 F_0^{\prime} (w,\boldsymbol{y})\nonumber\\&\quad\qquad\quad\qquad \quad\qquad + 2 \kappa _1 F_1^{\prime} (w,\boldsymbol{y}) ] w \delta [w - \epsilon ] \sqrt {g} \nonumber\\ &\approx \frac {\varepsilon ^2}{2} \int \mathrm{d}{\boldsymbol{y}}\, \mathrm{d}{\boldsymbol{z}}\, \mathrm{d}{w} \bigg ( [f_1 - F_1]^2\frac {\partial }{\partial w} \left ( \frac {\delta [w -\epsilon ]}{F_0^{\prime} } \right ) - \kappa _1^2 \frac {\partial (F_0^{\prime} \delta [ w - \epsilon ])}{\partial w}\nonumber\\&\quad\qquad \quad\qquad\quad\qquad - 2 \kappa _1 F_1^{\prime} \delta [w - \epsilon ] \bigg ) w \sqrt {g} . \end{align}

We integrate the terms that have a derivative with respect to $w$ by parts, giving

(4.22) \begin{align} \int\! \mathrm{d}{w} w (f_1 - F_1)^2\frac {\partial }{\partial w} \left ( \frac {\delta [w -\epsilon ]}{F_0^{\prime} } \right ) = -\!\int\! \mathrm{d}{w} [(f_1 - F_1)^2 - 2 w (f_1 - F_1) F_1^{\prime} ]\frac {\delta [w -\epsilon ]}{F_0^{\prime} }\nonumber\\[-10pt]\nonumber\\ \end{align}

and

(4.23) \begin{equation} \int \mathrm{d}{w} w \frac {\partial (F_0^{\prime} \delta [ w - \epsilon ])}{\partial w} = -\int \mathrm{d}{w} F_0^{\prime} \delta [ w - \epsilon ] . \end{equation}

Combining all results, we find that

(4.24) \begin{align} A \approx \frac {\varepsilon ^2}{2} \int \mathrm{d}{\boldsymbol{y}}\, \mathrm{d}{\boldsymbol{z}}\, \mathrm{d}{w} \left [ -\frac {(f_1-F_1)^2}{F_0^{\prime} } + \kappa _1^2 F_0^{\prime} + 2w F_1^{\prime} \left ( \frac {f_1-F_1}{F_0^{\prime} } - \kappa _1\right ) \right ]\delta [w-\epsilon ] \sqrt {g} , \end{align}

where the last term vanishes due to (4.16).

All in all, we have thus shown that

(4.25) \begin{align} A \approx \frac {\varepsilon ^2}{2} \int \mathrm{d}{\boldsymbol{x}}\, \frac {(f_1(\boldsymbol{x})-F_1[\epsilon ,\boldsymbol{y}])^2}{|F_0^{\prime} [\epsilon ,\boldsymbol{y}]|} - \frac {\varepsilon ^2}{2} \int \mathrm{d}{\boldsymbol{x}}\, \kappa _1^2(\boldsymbol{r}) |F_0^{\prime} [\epsilon ,\boldsymbol{y}]| +\mathcal{O}(\varepsilon ^3). \end{align}

The first term is identical to the result in Helander & Mackenbach (Reference Helander and Mackenbach2024) and is closely related to the notion of Helmholtz free energy. This term is furthermore always positive definite, so that if $\kappa =0$ , the available energy is always greater than or equal to zero. The second term encoding the constraint on the density can only serve to decrease the available energy, which is manifestly true as a more constrained calculation of it can only lower the value. It thus encodes what the energetic penalty is of some perturbation to the ground-state density, $n_1$ , which couples to $\kappa _1$ via (4.8) and (4.16).

Another useful form of the available energy may be obtained by substituting in (4.8), so that the result becomes

(4.26) \begin{align} A \approx \frac {\varepsilon ^2}{2} \int \mathrm{d}{\boldsymbol{x}} \left [ \frac {(f_1(\boldsymbol{x})-F_1[\epsilon ,\boldsymbol{y}])^2}{|F_0^{\prime} [\epsilon ,\boldsymbol{y}]|} - \left (\frac {n_1(\boldsymbol{r}) - \langle F_1 \rangle _{\boldsymbol{v}} }{|\langle F_0^{\prime} \rangle _{\boldsymbol{v}} |} \right )^2 F_0^{\prime} [\epsilon ,\boldsymbol{y}] \right ]\!. \end{align}

This can be simplified further by noting that $\int F_0^{\prime} [\epsilon ,\boldsymbol{y}] \mathrm{d}{\boldsymbol{x}} = \int \langle F_0^{\prime} \rangle _{\boldsymbol{v}} \mathrm{d}{\boldsymbol{r}}$ , so that we haveFootnote ⁷

(4.27) \begin{align} A \approx \frac {\varepsilon ^2}{2} \int \mathrm{d}{\boldsymbol{x}}\, \frac {(f_1(\boldsymbol{x})-F_1[\epsilon ,\boldsymbol{y}])^2}{|F_0^{\prime} [\epsilon ,\boldsymbol{y}]|} - \frac {\varepsilon ^2}{2} \int \mathrm{d}{\boldsymbol{r}}\, \frac {(n_1(\boldsymbol{r}) - \langle F_1 \rangle _{\boldsymbol{v}})^2}{|\langle F_0^{\prime} \rangle _{\boldsymbol{v}}|} . \end{align}

In the above form, we require only to solve for $F_1$ using (4.18) and we may calculate the available energy. We now solve (4.18) and (4.27), or equivalently (4.16), (4.9) and (4.25) in a number of explicit cases.

4.2. Maxwellian initial condition with no invariants

Assume that one is given a distribution function that is initially Maxwellian with slightly varying density and temperature across the domain, and we have no additional invariants (i.e. $\boldsymbol{y} = \varnothing$ ). As such, set $n(\boldsymbol{r}) = n_0(1+\nu )$ and $T(\boldsymbol{r})=T_0(1+\tau )$ , where we assume that $\nu (\boldsymbol{r}) \sim \tau (\boldsymbol{r}) \sim \varepsilon \ll 1$ and $\int \mathrm{d}{\boldsymbol{r}} \nu = \int \mathrm{d}{\boldsymbol{r}} \tau = 0$ . Expanding the Maxwellian distribution function, $f_M = n(m/2\pi T)^{3/2}{\rm e}^{-\epsilon /T}$ , gives

(4.28) \begin{equation} f_1(\boldsymbol{x}) = f_{M,0} \left ( \nu + \left [ \frac {\epsilon }{T_0} - \frac {3}{2} \right ] \tau \right )\!. \end{equation}

Let the density perturbation to the ground state be denoted as $\tilde {\nu } = n_1/n_0$ , so that (4.18) becomes

(4.29) \begin{equation} \frac {F_1(w)}{f_{M,0}} - \frac {\langle F_1 \rangle _{\boldsymbol{v}} }{n_0 } = \frac {1}{ \varOmega ^{\prime} (w) } \int \mathrm{d}{\boldsymbol{x}}\left ( \nu - \tilde {\nu } + \left [ \frac {\epsilon }{T_0} - \frac {3}{2} \right ] \tau \right ) \delta [w - \epsilon ] . \end{equation}

Choosing the density perturbation so that $\int \mathrm{d}{\boldsymbol{r}} \tilde {\nu } = 0$ (conserving the number of particles) causes the right-hand side terms to vanish, giving

(4.30) \begin{equation} \frac {F_1[\epsilon ]}{f_{M,0}[\epsilon ]} = \frac {\int \mathrm{d}{\boldsymbol{v}} F_1[\epsilon ] }{n_0}, \end{equation}

and we thus conclude that $F_1 = C f_{M,0}$ , where $C \sim \mathcal{O}(\varepsilon )$ is some constant.Footnote ⁸ Filling this result into the available energy results in

(4.31) \begin{align} A &= \frac {T_0}{2} \int \mathrm{d}{\boldsymbol{x}} \left ( \nu + \left [ \frac {\epsilon }{T_0} - \frac {3}{2}\right ] \tau - C \right )^2 f_{M,0} - \frac {n_0T_0}{2} \int \mathrm{d}{\boldsymbol{r}} \left ( \tilde {\nu } - C \right )^2 \nonumber\\ &= \frac { n_0 T_0 V}{2} \left \langle \nu ^2 - \tilde {\nu }^2 + \frac {3}{2} \tau ^2 \right \rangle , \end{align}

where angular brackets denote a volume average. Setting $\tilde {\nu } = \nu$ we recover the result from Helander (Reference Helander2020), who found that for a fixed density profile (i.e. the initial and final density profile are equal), the available energy is independent of $\nu$ and is equal to zero if there are no temperature perturbations. It turns out that the Maxwellian is actually a special case: any other choice of the initial distribution function with only density fluctuations that are furthermore held fixed would have resulted in non-zero available energy, as is proven in Appendix E. Two such examples shall be discussed in §§ 4.6 and 4.8.

Expression (4.31) can furthermore be used to bound the root-mean-square amplitude of the density fluctuations, which is maximal if the entire unconstrained available energy goes into density fluctuations of the ground state,

(4.32) \begin{equation} \langle \tilde {\nu }^2 \rangle ^{1/2} \leqslant \left \langle \nu ^2 + \frac {3}{2} \tau ^2 \right \rangle ^{1/2}. \end{equation}

In order to couple the ground state to the fields one can take both the ion and electron response into account, which we do in the next section.

4.3. Multiple Maxwellian species in free space

Suppose that we are given a ion–electron plasma with differing temperatures, sitting in free space so that the only invariants are the Casimirs. The available energy is given as the sum of the individual species in (4.31), thus,

(4.33) \begin{equation} A = A_{\mathrm{G}} - \frac { n_0 T_i V}{2} \langle \tilde {\nu }_i^2 \rangle - \frac { n_0 T_e V}{2} \langle \tilde {\nu }_e^2 \rangle , \end{equation}

where the subscript $i$ and $e$ denotes the particle species (ion or electron) and we remind ourselves that $A_{\mathrm{G}}$ is the ‘Gardner’ available energy where $\tilde {\nu }_i = \tilde {\nu }_e = 0$ . The electric potential generated by the ground states is given as

(4.34) \begin{equation} -\lambda _{D,e}^2 \boldsymbol{\nabla} ^2 \hat {\phi } = \tilde {\nu }_{i} - \hat {\phi }, \end{equation}

where we have assumed an adiabatic electron response $\tilde {\nu }_e \approx \mathrm{e} \phi / T_e = \hat {\phi }$ . Substituting (4.34) for $\tilde {\nu }_i$ and employing the adiabatic electrons response gives the following equation to ensure energetic consistency, i.e. $A = E_{\phi }$ :

(4.35) \begin{equation} A_{\mathrm{G}} - \frac { n_0 T_i V}{2} \langle ( \hat {\phi } - \lambda _{D,e}^2 \boldsymbol{\nabla} ^2 \hat {\phi } )^2 \rangle - \frac { n_0 T_e V}{2} \langle \hat {\phi }^2 \rangle = \frac {n_0 T_e V}{2} \langle ( \lambda _{D,e} \boldsymbol{\nabla }\hat {\phi } )^2 \rangle . \end{equation}

Furthermore, we wish to maximise the left-hand side subject to this equality constraint. This is most readily done in Fourier space, defined via $\hat {\phi }=\int \mathrm{d}{\boldsymbol{k}}\, \hat {\phi }_{\boldsymbol{k}} {\rm e}^{\mathrm{i} \boldsymbol{k} \boldsymbol{\cdot }\boldsymbol{r}}$ , so that invoking Plancherel’s theorem gives the constraint

(4.36) \begin{align} A_{\mathrm{G}} &= \frac {n_0 T_i (2 \pi )^d}{2} \int \mathrm{d}{\boldsymbol{k}}\, | \hat {\phi }_{\boldsymbol{k}} |^2 \left [\big(1 + \lambda _{D,e}^2 k^2\big)^2 + \frac {T_e}{T_i} \big( 1 + \lambda _{D,e}^2k^2 \big) \right ]\!, \end{align}

where $d$ is the dimensionality. As with (3.15) and (3.16), (4.36) is automatically satisfied by any distribution

(4.37) \begin{equation} \varphi (\boldsymbol{k}) = \frac {n_0 T_i (2 \pi )^d}{2 A_{\mathrm{G}}} | \hat {\phi }_{\boldsymbol{k}} |^2 \left [\big(1 + \lambda _{D,e}^2 k^2\big)^2 + \frac {T_e}{T_i} \big( 1 + \lambda _{D,e}^2k^2 \big) \right ] \end{equation}

that has unit density, i.e. $\int \mathrm{d}{\boldsymbol{k}} \varphi (\boldsymbol{k}) = 1$ . The available energy is equal to

(4.38) \begin{equation} \frac {A}{A_{\mathrm{G}}} = \frac {T_e}{T_i} \int \mathrm{d}{\boldsymbol{k}}\, \varphi (\boldsymbol{k}) \frac { \lambda _{D,e}^2k^2}{\big(1 + \lambda _{D,e}^2 k^2\big)^2 + ({T_e}/{T_i}) \big( 1 + \lambda _{D,e}^2k^2 \big) } \leqslant 1. \end{equation}

The available energy may be maximised by condensing all modes onto a single wavenumber, $k=k_0$ , that maximises the fraction. The wavenumber is given as

(4.39) \begin{equation} \lambda _{D,e}^2k_0^2 = \sqrt {1 + \frac {T_e}{T_i}}, \end{equation}

so that the available energy becomes

(4.40) \begin{equation} A = \left ( 1 - \frac {4}{2 + 2 \sqrt {1 + T_e/T_i}} \right ) A_{\mathrm{G}} , \end{equation}

as seen in figure 7. For equal temperatures, we find that the Gardner bound is multiplied by a factor $3 - 2\sqrt {2} \approx 0.17$ , a significant decrease highlighting the importance of taking the field equations into consideration. In the limit $T_e/T_i \rightarrow \infty$ , the available energy tends to the Gardner bound, whereas in the limit $T_e/T_i \rightarrow 0$ the available energy tends to zero, for the same reasons as discussed in § 3. However, electrons can hardly be considered adiabatic in the limit as the ratio of thermal velocities, $v_e/v_i = \sqrt {m_iT_e/m_e T_i}$ , tends to zero.

Figure 7.

The factor by which the available energy is reduced in a two-species plasma with adiabatic electrons, as a function of the ratio of temperatures.

Figure 7. Long description

A line graph showing the factor by which the available energy is reduced in a two-species plasma with adiabatic electrons, as a function of the ratio of temperatures. The x-axis represents the ratio of electron temperature to ion temperature, ranging from 0 to 10. The y-axis represents the ratio of available energy to the initial available energy, ranging from 0 to 1. The graph shows a curve that starts at the origin and increases gradually, indicating that as the ratio of temperatures increases, the available energy is reduced by a smaller factor.

One can similarly bound the fluctuation amplitude of a consistent ground state by combining (4.34) and (4.37). The ion-density fluctuation can be written as

(4.41) \begin{equation} |\tilde {\nu }_{i,\boldsymbol{k}} |^2 = \varphi (\boldsymbol{k}) \frac {2A_{\mathrm{G}}}{n_0 T_i (2 \pi )^d} \frac {\big(1+\lambda _{D,e}^2k^2\big)}{\big(1 + \lambda _{D,e}^2 k^2\big)^2 + ({T_e}/{T_i}) \big( 1 + \lambda _{D,e}^2k^2 \big)}, \end{equation}

which is maximal for a Dirac-delta $\varphi$ , centred on $k^2 \rightarrow 0$ . This yields a maximal fluctuation amplitude of

(4.42) \begin{equation} \langle \tilde {\nu }_{i}^2 \rangle \leqslant \frac {2 A_{\mathrm{G}}}{n_0 T_i V} \frac {1}{1+ T_e/T_i}, \end{equation}

closely related to the result of (4.32), where the electric field was not accounted for. Including constraints imposed by the field reduces the bound by a factor dependent on the ratio of temperatures.

If instead both electrons and ions are kinetic, one gains an additional degree of freedom in order to maximise the available energy, namely $\tilde {\nu }_e$ . This may be employed to interpret $T_e/T_i$ as an additional parameter to optimise over where the maximal value of (4.33) is $A = A_{\mathrm{G}}$ , implying that for kinetic electrons the Gardner bound is attainable (where dynamical accessibility is, of course, not guaranteed).

4.4. Charged particles with a given ground-state density profile in a magnetic field

Here the available energy of ions in the presence of a magnetic field is calculated by assuming that $\boldsymbol{y} = \{\mu \}$ is conserved, the density profile is furthermore kept fixed. Such a model may, for example, be useful in the presence of adiabatic electrons (i.e. $n \propto \phi$ ) as there can be no advection of electrons down the gradient, $\boldsymbol{E} \times \boldsymbol{B} \boldsymbol{\cdot }\boldsymbol{\nabla }n = 0$ . When $\mu$ is a conserved quantity, the phase-space volume element is given as

(4.43) \begin{equation} \mathrm{d} \boldsymbol{x} = \frac {2 \pi B}{m} \; \mathrm{d} \mu \mathrm{d} v_\| \mathrm{d} \boldsymbol{r}, \end{equation}

so that (4.16) becomes

(4.44) \begin{equation} \int _{w\gt \mu B} \mathrm{d}{\boldsymbol{r}} \frac {B(f_1 - F_1 - \kappa _1 F_0^{\prime} )}{\sqrt {w-\mu B}} = 0. \end{equation}

Equation (4.44) does not permit analytical solutions in general, but may be solved if the magnetic field is assumed to be constant, giving

(4.45) \begin{equation} F_1 = \langle f_1 - \kappa _1 F_0^{\prime} \rangle . \end{equation}

The Lagrange multiplier $\kappa _1$ is, in turn, determined via (4.9), yielding

(4.46) \begin{equation} \kappa _1 = \frac {n_1- C n_0 }{\int \mathrm{d}{\boldsymbol{v}} F_0^{\prime} }, \end{equation}

where we have introduced some unknown constant of order $\varepsilon$ , $\int \mathrm{d}{\boldsymbol{v}} F_1 = C n_0$ .

Taking the initial distribution function to be bi-Maxwellian,

(4.47) \begin{equation} f_{B} = n \left ( \frac {m}{2 \pi T_\perp ^{2/3}T_\|^{1/3}} \right )^{3/2} \exp \left [ -\frac {\epsilon }{T_{\|}} + \mu B \left ( \frac {1}{T_{\|} } - \frac {1}{T_{\perp }} \right )\right ]\!, \end{equation}

where $T_{\|}$ and $T_{\perp }$ measure the parallel and perpendicular temperatures. Setting $n=n_0(1+\nu )$ , $T_\|=T_{\|,0}(1+\tau _\|)$ , $T_\perp =T_{\perp ,0}(1+\tau _\perp )$ and expanding around the smallness of $\nu \sim \tau _\perp \sim \tau _\| \sim \varepsilon$ gives

(4.48) \begin{equation} f_1 = f_{B,0} \left [ \nu + \left ( \frac {\mu B}{T_{\perp ,0}}-1 \right )\tau _{\perp } + \left ( \frac {\epsilon - \mu B}{T_{\|,0}}- \frac {1}{2}\right )\tau _{\|} \right ]\!. \end{equation}

Noting that $\partial _\epsilon f_{B,0} = - f_{B,0}/T_{\|,0}$ , it follows that the ground state is described as

(4.49a,b) \begin{align} \kappa _1 = T_{\|,0} \left ( C - \tilde {\nu } \right ), \quad F_1 = f_{B,0} \left [ C + \left ( \frac {\mu B}{T_{\perp ,0}}-1 \right )\langle \tau _{\perp } \rangle + \left ( \frac {\epsilon - \mu B}{T_{\|,0}}- \frac {1}{2}\right ) \langle \tau _{\|} \rangle \right ]\!,\nonumber\\[-10pt]\nonumber\\ \end{align}

where we have used $\langle \nu \rangle = \langle \tilde {\nu } \rangle = 0$ . Filling these results into (4.25) gives

(4.50) \begin{equation} A = \frac {V n_0 T_{\|,0}}{2} \left \langle \nu ^2 - \tilde {\nu }^2 +\frac {1}{2} \left ( \tau _\|^2 - \langle \tau _\|^2 \rangle \right ) + \tau _\perp ^2 - \langle \tau _\perp ^2 \rangle \right \rangle , \end{equation}

which generalises the result of Helander & Mackenbach (Reference Helander and Mackenbach2024) to $\tilde {\nu } \neq 0$ . Somewhat surprisingly, we see that for a plasma with isotropic temperatures, the result reduces to the case without invariants. Thus, an isotropic plasma that is in its ground state with a given density profile frees up the same amount of available energy to leading order in closeness to a ground state, regardless if a constant magnetic field is present or not. A spatially varying magnetic field and/or the introduction of anisotropy reduces the available energy below its unmagnetised counterpart. However, Poisson’s equation is adjusted when a sufficiently strong magnetic field is present so that a consistent coupling of ground states leads to differing behaviours of it, as we see in the next section.

4.5. Ions and electrons in a constant magnetic field

We consider again an ion–electron plasma now in the presence of a constant magnetic field, a so-called slab. To make further progress, we will take finite-Larmor-radius effects to be small, i.e. $ \rho _i k_\perp \ll 1$ , where $k_\perp$ is the wave vector perpendicular to the magnetic field and $\rho _i = \sqrt {m_i T_i}/(\mathrm{e} B )$ is the ion Larmor radius ( $m_i$ being the ion mass). The invariant energy in the electrostatic, long wavelength limit of gyrokinetics is given as (Dubin et al. Reference Dubin, Krommes, Oberman and Lee1983, see also Helander Reference Helander2017, Appendix A)

(4.51) \begin{align} E_{\mathrm{tot}} \approx & \int \mathrm{d}{\boldsymbol{x}} \epsilon f_i + \int \mathrm{d}{\boldsymbol{x}} \epsilon f_e + \int \mathrm{d}{\boldsymbol{r}} \frac {|\boldsymbol{\nabla }\phi |^2}{8 \pi } + \frac {m_i n_0}{2 B^2} \int \mathrm{d}{\boldsymbol{r}} | \boldsymbol{b} \times \boldsymbol{\nabla }\phi |^2, \end{align}

where the unit magnetic-field vector is given as $\boldsymbol{b} = \boldsymbol{B}/B$ and $\epsilon = \mu B + mv_\|^2/2$ . Besides the field energy, a kinetic ‘sloshing’ energy (the last term in (4.51)) is now present also, so that the available energy bounds the sum. Furthermore, Poisson’s equation with adiabatic electrons and small $\rho _i k_\perp$ becomes

(4.52) \begin{align} \tilde {\nu }_{i,\boldsymbol{k}} &\approx \left ( \frac {T_i}{T_e} + \rho _i^2 k_\perp ^2 + \lambda _{D,i}^2 k^2 \right ) \hat {\phi }_{\boldsymbol{k}} \nonumber\\ &\approx \varPhi _i(\boldsymbol{k}) \hat {\phi }_{\boldsymbol{k}}, \end{align}

where $\tilde {\nu }_{i,\boldsymbol{k}}$ denotes the $\boldsymbol{k}$ th Fourier component of density fluctuations associated with the ground state of $f_i$ , the electrostatic potential is normalised as $\hat {\phi }_{\boldsymbol{k}} = \mathrm{e} \phi _{\boldsymbol{k}} / T_i$ and the shorthand $\varPhi (\boldsymbol{k})$ is used for the factor in brackets.Footnote ⁹ Assume that the initial distribution $f_i$ is given as a Maxwellian (as is $f_e$ ) with slight variations in temperature and density, so that the sum of the two available energies gives the total

(4.53) \begin{equation} A = A_{\mathrm{G}} - \frac {n_0 T_i (2\pi )^3}{2} \int \mathrm{d}{\boldsymbol{k}}\, | \varPhi _i(\boldsymbol{k})\hat {\phi }_{\boldsymbol{k}} |^2 - \frac {n_0 T_i (2\pi )^3}{2} \int \mathrm{d}{\boldsymbol{k}}\, \frac {T_i}{T_e} |\hat {\phi }_{\boldsymbol{k}}|^2. \end{equation}

In order to satisfy the conservation of (4.51), one must have

(4.54) \begin{equation} A_{\mathrm{G}} = \frac {n_0 T_i (2 \pi )^3}{2} \int \mathrm{d}{\boldsymbol{k}}\, | \hat {\phi }_{\boldsymbol{k}} |^2 \left ( |\varPhi _i(\boldsymbol{k})|^2 + \frac {T_i}{T_e} +\rho _i^2 k_\perp ^2 + \lambda _{D,i}^2 k^2 \right )\!, \end{equation}

which – as in (3.15) and (3.16) – may be achieved by defining the distribution

(4.55) \begin{equation} \varphi (\boldsymbol{k}) = \frac {n_0 T_i (2\pi )^3}{2A_{\mathrm{G}}} | \hat {\phi }_{\boldsymbol{k}} |^2 \left ( |\varPhi _i(\boldsymbol{k})|^2 + \frac {T_i}{T_e} + \rho _i^2 k_\perp ^2 + \lambda _{D,i}^2 k^2 \right )\!, \end{equation}

whose integral over wavenumber space satisfies $\int \mathrm{d}{\boldsymbol{k}} \varphi ({\boldsymbol{k}}) = 1$ . The available energy then becomes

(4.56) \begin{equation} \frac {A}{A_{\mathrm{G}}} = \int \mathrm{d}{\boldsymbol{k}}\, \varphi (\boldsymbol{k}) \frac { \rho _i^2 k_\perp ^2 + \lambda _{D,i}^2 k^2 }{ |\varPhi _i(\boldsymbol{k})|^2 + ({T_i}/{T_e}) + \rho _i^2 k_\perp ^2 + \lambda _{D,i}^2 k^2} \leqslant 1 , \end{equation}

which again may be maximised by finding the wave vector $\boldsymbol{k}=\boldsymbol{k}_0$ that does so. Setting the magnitude of the wavenumber $k^2 = k_\perp ^2 + k_z^2$ where the perpendicular component is $k_\perp ^2 = k_x^2 + k_y^2$ (i.e. the metric tensor is diagonal and constant) yields a maximising perpendicular wavenumber

(4.57) \begin{equation} \big(\rho _i^2 + \lambda _{D,i}^2\big)k_\perp ^2 =-\lambda _{D,i}^2 k_z^2 + \sqrt {\frac {T_i^2}{T_e^2} + \frac {T_i}{T_e}}, \end{equation}

where solubility requires that $\lambda _{D,i}^2k_z^2 \leqslant \sqrt {({T_i^2}/{T_e^2}) + ({T_i}/{T_e})}$ . Note the close correspondence to the result given in (4.39), as setting $\rho _i = 0$ gives the same maximising wavenumber magnitude. Indeed, the available energy is unchanged and becomes

(4.58) \begin{equation} A = \left ( 1 - \frac {4}{2 + 2 \sqrt {1 + T_e/T_i}} \right ) A_{\mathrm{G}}, \end{equation}

the exact same result as in (4.38). This is perhaps unsurprising, as one can still deposit the energy in electrostatic modes in this regime, though it is no longer solely in the form of the field energy. The energy is distributed between the field and sloshing terms as

(4.59) \begin{equation} \frac {\int \mathrm{d}{\boldsymbol{r}}\, |\boldsymbol{\nabla }\phi |^2/8\pi }{\int \mathrm{d}{\boldsymbol{r}}\, m_i n_0 |\boldsymbol{b} \times \boldsymbol{\nabla }\phi |/2B^2} = \frac {\lambda _{D,i}^2}{\rho _i^2} \left (\frac { \sqrt {({T_i^2}/{T_e^2}) + ({T_i}/{T_e})} + \rho _i^2 k_z^2}{\sqrt { ({T_i^2}/{T_e^2}) + ({T_i}/{T_e}}) - \lambda _{D,i}^2 k_z^2} \right )\!. \end{equation}

Since strongly magnetised plasmas tend to be very anisotropic, $\rho _i k_z \sim \rho _i/L \ll 1$ and $\lambda _{D,i} k_z \sim \lambda _{D,i}/L \ll 1$ , where $L$ is some global length scale (in a magnetic confinement device $L$ is of the order of the major radius), the ratio of the field-to-sloshing energy is approximately $\lambda _{D,i}^2/\rho _i^2$ . In most practical applications $\lambda _{D,i} \ll \rho _i$ , meaning that one can safely neglect $\lambda _{D,i}$ effects in the ground state.

Increasing $T_e/T_i$ in its large-value limit we see two behaviours of interest: firstly, for $\lambda _{D,i}= 0$ , the maximising wavenumber (4.57) goes as $\rho _i^2 k_\perp ^2 \approx \sqrt {T_i/T_e}$ , confirming that the long-wavelength assumption is satisfied here, and one could furthermore expect that the transport is greatly enhanced by being concentrated in long wavelengths. A similar behaviour is present in the simple Vlasov–Poisson plasma, as for fixed ion-temperature and large $T_e/T_i$ , (4.39) goes as $\lambda _{D,i}^2k^2 \approx \sqrt {T_i/T_e}$ . Secondly, the transport could be exacerbated further when increasing $T_e/T_i$ as the available energy (4.58) goes as $A \approx A_{\mathrm{G}}(1 - 2\sqrt {T_i/T_e})$ , approaching the Gardner bound as the electron temperature is raised at fixed ion temperature. The amplification of transport with increasing $T_e$ at fixed $T_i$ is a known feature of gyrokinetic turbulence having been found both in dispersion relations and simulations (Biglari, Diamond & Rosenbluth Reference Biglari, Diamond and Rosenbluth1989; Romanelli Reference Romanelli1989; Plunk et al. Reference Plunk, Helander, Xanthopoulos and Connor2014; Zocco et al. Reference Zocco, Xanthopoulos, Doerk, Connor and Helander2018; Helander & Plunk Reference Helander and Plunk2022), and may be a possible culprit degrading energy confinement in both tokamaks and stellarators that are electron heated (Beurskens et al. Reference Beurskens2021). In the simple model considered, we may understand it simply as the upper bound on the combined field and kinetic energy growing, which can lead to increased transport if the system approaches the bound (in value or in scaling), where an increase in the maximising wavelength furthermore contributes.

4.5.1. Zonal flows

The fact that the ground state derived in the preceding section allows for radial streamers (that is, $k_x=0$ modes) is reminiscent of results found by Hammett et al. (Reference Hammett, Beer, Dorland, Cowley and Smith1993), who simulated ion-temperature-gradient-driven turbulence using a nonlinear gyrofluid model assuming that $\delta n_e \propto \phi$ , and found that the turbulent state does not saturate due to these streamers. Saturation is achieved by correcting the adiabatic response,

(4.60) \begin{equation} \frac {\delta n_e}{n_0} = \frac {\mathrm{e} (\phi - \langle \phi \rangle _{\mathrm{fs}})}{T_e}, \end{equation}

where the flux-surface average of the electrostatic potential, denoted as $\langle \phi \rangle _{\mathrm{fs}}$ , is subtracted due to the rapid circulation of electrons around the flux surface. It is worthwhile to compare this behaviour to the results found in § 3: for purely zonal potentials ( $\phi = \langle \phi \rangle _{\mathrm{fs}}$ ), the background density is unperturbed and acts like the neutralising case of the waterbags. In the language of zonal-flow literature, this is because such flows are modes of minimal inertia (Diamond et al. Reference Diamond, Itoh, Itoh and Hahm2005): they are not screened by the Boltzmann response of the electrons. Focusing on the ground state, by concentrating the energy into arbitrarily long wavelength zonal flows we thus expect that the Gardner bound can be reached, which is indeed what we find.

The flux-surface average in Fourier space is equivalent to setting $\hat {\phi }_{\mathrm{Z}} = \hat {\phi }_{\boldsymbol{k}}(k_x,$ $k_y=0,k_z=0)$ (where the subscript $\mathrm{Z}$ stands for the zonal component), thus retaining radial dependence alone. Similarly, using this response in our calculation modifies Poisson’s equation to

(4.61) \begin{equation} \tilde {\nu }_{i,\boldsymbol{k}} \approx \begin{cases} \left (\varPhi _i(\boldsymbol{k}) - \dfrac {T_i}{T_e} \right )\hat {\phi }_{\mathrm{Z}} & \text{if } k_y=k_z=0, \\[3pt] \varPhi _i(\boldsymbol{k}) \hat {\phi }_{\boldsymbol{k}} & \text{otherwise}. \end{cases} \end{equation}

Supposing that the potential is purely zonal ( $k_y=k_z=0$ , consequently $\hat {\phi }_{\boldsymbol{k}}=\hat {\phi }_{\mathrm{Z}}$ ), we may proceed with the calculation as before, with the total available energy (4.53) being adjusted to

(4.62) \begin{equation} A = A_{\mathrm{G}} - \frac {n_0 T_i (2\pi )}{2} \int \mathrm{d}{k_x} \left ( \varPhi _i(k_x) -\frac {T_i}{T_e} \right )^2 |\hat {\phi }_{\mathrm{Z}} |^2, \end{equation}

whereas conservation of energy (4.54) now gives

(4.63) \begin{equation} A_{\mathrm{G}} = \frac {n_0 T_i (2 \pi )}{2} \int \mathrm{d}{k_x} \left [ \left (\varPhi _i(k_x) - \frac {T_i}{T_e} \right )^2 +\rho _i^2 k_x^2 + \lambda _{D,i}^2 k_x^2 \right ] | \hat {\phi }_{\mathrm{Z}} |^2. \end{equation}

As before, we define the distribution $\varphi (k_x)$ ,

(4.64) \begin{equation} \varphi (k_x) = \frac {n_0 T_i (2\pi )}{2A_{\mathrm{G}}} \left [ \left (\varPhi _i(k_x) - \frac {T_i}{T_e} \right )^2 +\rho _i^2 k_x^2 + \lambda _{D,i}^2 k_x^2 \right ]| \hat {\phi }_{\mathrm{Z}} |^2, \end{equation}

so that the part of the Gardner bound that may be released becomes

(4.65) \begin{equation} \frac {A}{A_{\mathrm{G}}} = \int \mathrm{d}{k_x} \frac { \varphi (k_x)}{ 1 + \big(\rho _i^2 + \lambda _{D,i}^2\big) k_x^2 } \leqslant 1. \end{equation}

For a Dirac-delta $\varphi$ at $k_x = k_{x,0}$ , we thus find that

(4.66) \begin{equation} \frac {A}{A_{\mathrm{G}}} = \frac {1}{1+\big(\rho _i^2+\lambda _{D,i}^2\big)k_{x,0}^2} \end{equation}

is completely analogously to the waterbag calculation with a neutralising background, where the Gardner bound is reached for $k_{x,0}\rightarrow 0$ . This gives an interesting example of a state that has no transport due to the stabilising zonal structure, despite having reached the Gardner bound. Furthermore, the ground state is a zonally dominated state, as is the case in the Dimits regime (Dimits et al. Reference Dimits2000).

4.6. Impurities subject to quasineutrality

We consider a crude model to capture the available energy of three different Maxwellian species of particles in a toroidal geometry, namely electrons (trapped and passing), hydrogen and an additional impurity species with charge number $Z$ . We model the trapped electrons and both ion species as having no invariants aside from the Casimirs, whereas the passing electron species due to its fast motion along the field line (averaging the radial drift to zero) conserves its flux-surface label $\psi$ . This is, of course, a gross oversimplification. In reality, all species obey conservation of $\mu$ , and the trapped electron species furthermore conserves the parallel invariant as treated in Mackenbach et al. (Reference Mackenbach, Proll, Wakelkamp and Helander2023). However, we are not too interested in capturing the Gardner value of the available energy accurately here, and we instead wish to focus on how the constraint of quasineutrality reduces this Gardner bound, continuing as such.

In a toroidal geometry the density and temperature are flux functions, $\nu _{e,p} = \nu _{e,p}(\psi )$ and $\tau _{e,p} = \tau _{e,p}(\psi )$ , rendering it impossible for the passing electron species to liberate available energy (as the plasma distribution function is a ground state at fixed $\psi$ , consequently keeping the density and temperature profiles of these particles fixed). In order to have quasineutrality, one thus requires that

(4.67) \begin{equation} f_p\nu _{e,p}(\psi ) + (1-f_p) \tilde {\nu }_{e,t}(\psi ) = \frac {n_H}{n_e}\tilde {\nu }_{H}(\psi ) + \frac {Zn_Z}{n_e} \tilde {\nu }_{Z}(\psi ), \end{equation}

where the subscripts $e$ , $H$ and $Z$ denote whether we are dealing with electrons, hydrogen or the remaining ionic species, respectively. Furthermore, subscripts $t$ and $p$ indicates if the electrons are trapped or passing and $f_p \in [0,1]$ denotes the passing fraction (related to the trapping fraction $f_t$ as $f_t = 1-f_p$ ). The available energy, assuming initially Maxwellian distribution functions with slightly varying density and temperature, becomes

(4.68) \begin{equation} A = A_{\mathrm{G} } - \frac {n_e T_{e} V (1-f_p)^2}{2} \langle \tilde {\nu }_{e,t}^2 \rangle - \frac {n_H T_{H} V}{2} \langle \tilde {\nu }_{H}^2 \rangle - \frac {n_Z T_{Z} V}{2} \langle \tilde {\nu }_{Z}^2 \rangle . \end{equation}

It is evident that in order to maximise (4.68) the trapped electron density variation can be set to zero, leaving only the ion densities to optimise over. Solving (4.67) for $\tilde {\nu }_Z$ and substituting this expression into (4.68) gives

(4.69) \begin{equation} A = A_{\mathrm{G} } - \frac {n_e T_{H} V}{2} \langle \tilde {\nu }_{H}^2 \rangle - \frac {n_Z T_{Z} V}{2Z^2} \left \langle \left ( \frac { n_e f_p \nu _{e,p}}{n_Z} - \frac {n_H \tilde {\nu }_{H}}{n_Z} \right )^2 \right \rangle , \end{equation}

which is maximal for

(4.70) \begin{equation} \tilde {\nu }_{H}(\psi ) = \frac {f_p n_0 T_{Z}}{n_HT_{Z} + Z^2 n_Z T_{H}} \nu _{e,p}(\psi ). \end{equation}

Using this equation for the available energy and using the zeroth-order solution of the quasineutrality condition $n_e = n_H + Zn_Z$ , we obtain the result

(4.71) \begin{equation} A = A_{\mathrm{G}} - \frac {f_p^2 n_e^2 T_H T_Z V}{2( n_H T_Z + Z^2n_Z T_H ) } \langle \nu _{e,p}^2 \rangle . \end{equation}

It is clear that having a high charge-number species present in the plasma (at fixed number densities) is disadvantageous not only in terms of having higher radiation losses or rarefying the fusion fuel, it also increases the available energy and may thus enhance turbulence in this very crude model. The reason for this is simple: if the charge number is very high, one can satisfy (4.67) with small values of $\tilde {\nu }_H$ and $\tilde {\nu }_Z$ , requiring little density variation in the ground state and thus approaching the Gardner value. However, this argument ignores the fact that the Gardner bound itself may depend on $Z$ and other parameters, and we wish to analyse such dependencies in the following paragraph.

To investigate dependence on impurity content $n_Z$ , let us assume that the initial states all have variation in density alone, further assuming that $\nu _{e,t}=\nu _{e,p}=\nu _e$ , so that the Gardner bound becomes

(4.72) \begin{equation} A_{\mathrm{G}} = \frac {n_H T_H V}{2} \langle \nu _{H}^2 \rangle + \frac {n_e T_e f_p V}{2} \langle \nu _{e}^2 \rangle + \frac {n_Z T_Z V}{2} \langle \nu _{Z}^2 \rangle . \end{equation}

The initial states also have to be quasineutral ( $n_e=n_H+Zn_Z$ ), meaning that

(4.73) \begin{equation} \nu _{e}(\psi ) = \frac {n_H}{n_0}\nu _{H}(\psi ) + \frac {Zn_Z}{n_0} \nu _{Z}(\psi ). \end{equation}

The precise functional form of $A(n_Z)$ depends, of course, on what is held constant as $n_Z$ is varied, and we make the following choices: we use the quasineutrality equations to eliminate dependence on $n_e$ and $\nu _e$ , set $T_H=T_Z=T_e=T$ and express the impurity density in terms of the effective charge $Z_{\mathrm{eff}} = n_H/n_e + Z^2 n_Z/n_e$ , and define a relative available energy $\overline {A}=A/A(Z_{\mathrm{eff}}=0)$ . Further specialising to the case where $\langle \nu _H^2 \rangle = 1$ , $\langle \nu _H \nu _Z \rangle = 1/4$ , $\langle \nu _Z^2 \rangle = 1/16$ , $f_p=1/3$ and $Z=6$ gives a relative available energy dependent on only $Z_{\mathrm{eff}}$ (which itself depends on $n_H$ ), namely

(4.74) \begin{equation} \overline {A} = \frac {-486 Z_{\mathrm{eff}}-({1350}/{Z_{\mathrm{eff}}-6})-({529}/{Z_{\mathrm{eff}}})+5065}{4320}. \end{equation}

A plot of this function is given in figure 8, where we see a pronounced minimum in $\overline {A}$ for a certain $Z_{\mathrm{eff}}$ . Similar dependence was found in García-Regaña et al. (Reference García-Regaña, Calvo, Parra and Thienpondt2024) for the turbulent transport, though we stress that the parameters chosen here are not the same – we are more interested in seeing if qualitative behaviours are captured, since the model we have used is crude. The fact that qualitatively such trends arise may indicate that the intricate dependencies of turbulent transport on impurities can be, in part, captured by considerations of energy budgets in conjunction with quasineutrality.

Figure 8.

Plot of the relative available energy’s dependence on $Z_{\mathrm{eff}}$ .

4.7. A kappa distribution in free space

Let us now consider a non-Maxwellian distribution function that is nonetheless a ground state, namely the kappa distribution given as

(4.75) \begin{equation} f_\chi (\boldsymbol{x}) = n_0 \left ( \frac {m }{2 \pi T} \right )^{d/2} \frac {\varGamma ( \chi + 1 + {d}/{2} )}{ \chi ^{d/2}\varGamma ( \chi + 1)} \left ( 1 + \frac {1}{\chi } \frac {\epsilon }{T} \right )^{-\chi - 1 - d/2}, \end{equation}

where $d$ denotes the dimensionality (having replaced the typical $\kappa$ parameter with $\chi$ , since the former is used as a Lagrange multiplier). Perturbing the density and temperature slightly, $n = n_0(1+\nu )$ and $T = T_0 (1 + \tau )$ , gives the leading-order derivative

(4.76) \begin{equation} f_{\chi ,0}^{\prime} [\epsilon ] = - f_{\chi ,0}[\epsilon ] \frac {1 + \chi + {d}/{2}}{\epsilon +T_0 \chi }, \end{equation}

which is indeed negative definite for any $\chi \gt 0$ , implying it is a ground state. Going to the next order gives the perturbed distribution function

(4.77) \begin{equation} f_1(\boldsymbol{x}) = f_{\chi ,0} \left [ \nu + \tau \frac {{\epsilon }/{T_0}(1+\chi ) - {\chi d}/{2}}{{\epsilon }/{T_0} + \chi } \right ]\!. \end{equation}

As in § 4.1, using $\langle \tilde {\nu } \rangle =\langle \nu \rangle =\langle \tau \rangle = 0$ simplifies the ground-state equation (4.18) to

(4.78) \begin{equation} F_1[\epsilon ] = \frac { C \chi (1 + \chi + {d}/{2} )}{(1+\chi ) ( ({mv^2}/{2T_0}) + \chi )} f_{\chi ,0}, \end{equation}

having defined $\int \mathrm{d}{\boldsymbol{v}} F_1 =C n_0$ and having used $\int \mathrm{d}{\boldsymbol{v}} F_0^{\prime} = -n_0(1+\chi )/T_0\chi$ . The available energy simplifies to

(4.79) \begin{equation} A = \frac {n_0 T_0 V}{2} \left \langle \frac {d(\nu + \tau )^2 + \chi (2 \nu ^2 + d \tau ^2)}{2+d+2 \chi } \right \rangle - \frac {n_0 T_0 V}{2} \left \langle \frac {\chi \tilde {\nu }^2}{1 + \chi } \right \rangle\! . \end{equation}

The first term is the Gardner bound as derived in Helander & Mackenbach (Reference Helander and Mackenbach2024), whereas the second term encodes the penalty of the ground-state density fluctuation. As shown in Appendix E, even if the density profile is fixed and there are no temperature fluctuations ( $\tilde {\nu } = \nu$ and $\tau =0$ ) there is a non-zero available energy, i.e.

(4.80) \begin{equation} A = \frac {n_0 T_0 V}{2} \left \langle \frac {d\nu ^2 }{2(1+\chi )(2+d+2\chi )} \right \rangle , \end{equation}

meaning that energy can be extracted from density fluctuations even if it cannot flatten, indicating that there are phase-space rearrangements possible that release energy. (One such rearrangement for a simplified scenario is shown in figure 14 of Appendix E).

This may have implications for the stability properties of plasmas with non-Maxwellian features, adiabatic electrons and a pure density gradient. In Helander (Reference Helander2020) it is proven that an isothermal Maxwellian with a fixed density profile is nonlinearly stable due to its vanishing available energy. However, as we see here, for any other distribution function, there are phase-space rearrangements possible that liberate energy, meaning some instability can arise.

4.8. A three-dimensional slowing-down distribution in a constant magnetic field

Consider now a distribution function used often to model the fast alpha particles born from fusion reactions, the so-called slowing-down distribution. This distribution function is derived by imposing a Dirac-delta source of particles at a single velocity and solving the steady-state Fokker–Planck equation where the fast particles collide with a Maxwellian background species of electrons (Gaffey Jr Reference Gaffey1976; Moseev & Salewski Reference Moseev and Salewski2019), yielding

(4.81) \begin{equation} f_s(\boldsymbol{x}) = \frac {3n}{4 \pi \ln \big( 1 + {v_\alpha ^3}/{v_c^3} \big)} \frac {\varTheta ( v_\alpha - v )}{v^3 + v_c^3}. \end{equation}

Here $n(\boldsymbol{x})$ is the fast-particle density, $v_c(\boldsymbol{x})$ is the ‘cross-over’ velocity where the drag on fast ions and electrons is equal and $v_\alpha$ is the velocity of the born alpha particles (so that $m_\alpha v_\alpha ^2/2 \approx 3.5 \; \text{MeV}$ ). Note that $\partial _\epsilon f_s \leqslant 0$ , so that a spatially homogeneous slowing-down distribution is a ground state. We wish to find the corresponding ground state of the spatially varying slowing-down distribution function (though the turbulence itself may alter the distribution function Wilkie (Reference Wilkie2018) – an effect that we will ignore). A comment is warranted here, as the available energy tells us how much energy can be released collisionlessly, whereas the slowing-down distribution is created by a collisional process. As such, one can best interpret the results as a measure of how unstable such a collisional distribution is to collisionless instabilities. A proper analysis should take into account the ordering of the different time scales (slowing-down time, growth rate of collisionless instabilities, etc.), we shall not attempt to do so here.

Set $n = n_0( 1 + \nu )$ and $v_c = v_{c,0}(1 + \tau /2)$ since $v_c \propto \sqrt {T_e}$ , though $\tau$ may include spatial variation in the main ion and electron densities too (see Wilkie Reference Wilkie2015, Equation (1.22)). The perturbed distribution function is

(4.82) \begin{equation} f_{s,1} = f_{s,0} \left \{ \nu + \frac {3 \tau }{2} \left [ \frac {v_\alpha ^3}{\big(v_\alpha ^3 + v_{c,0}^3\big) \ln \big( 1 + {v_\alpha ^3}/{v_{c,0}^3} \big)} - \frac {v_{c,0}^3}{v_{c,0}^3 + v^3} \right ] \right \}\!. \end{equation}

The ground state in a constant magnetic field is given by (4.45), (4.46), and for these equations, we require $\partial _\epsilon f_{s,0}$ and its integral over velocity space. These may be evaluated analytically as

(4.83) \begin{equation} \partial _\epsilon f_{s,0} = -\frac {9 n_0}{4 \pi m \ln \big( 1 + {v_\alpha ^3}/{v_{c,0}^3} \big)} \left \{\frac {v \varTheta [v_\alpha - v]}{\big( v_{c,0}^3 + v^3 \big)^2} + \frac {\delta [v_\alpha - v]}{3v\big( v_{c,0}^3 + v^3 \big)} \right \} \end{equation}

and

(4.84) \begin{align} \int \mathrm{d}{\boldsymbol{v}} \partial _\epsilon f_{s,0} = -&\frac {n_0}{2 m v_{c,0}^2 \ln \big( 1 + {v_\alpha ^3}/{v_{c,0}^3} \big)} \times \Bigg ( \frac {\pi }{\sqrt {3}} - 2\sqrt {3} \cot ^{-1} \left [ \frac {v_{c,0} \sqrt {3}}{v_{c,0} - 2 v_{\alpha }} \right ] \nonumber\\ & +2 \ln [ v_{c,0} + v_{\alpha } ] - \ln \big[ v_{c,0}^2 + v_\alpha ^2 - v_{c,0} v_\alpha \big] \Bigg ) \nonumber\\ \equiv & -\frac {n_0}{2m v_{c,0}^2 \ln \big( 1 +{v_\alpha ^3}/{v_{c,0}^3} \big)} \times \mathcal{E}(v_{c,0},v_\alpha ), \end{align}

where $\cot ^{-1}(x)$ is the inverse of $\cot (x)=1/\tan (x)$ , and having defined the factor in brackets as $\mathcal{E}(v_{c,0},v_\alpha )$ . As such, the Lagrange multiplier satisfies

(4.85) \begin{equation} \kappa _1 = 2 m v_{c,0}^2 \frac { C - \tilde {\nu } }{ \mathcal{E}(v_{c,0},v_\alpha ) } \ln \left ( 1 + \frac {v_\alpha ^3}{v_{c,0}^3} \right )\!, \end{equation}

so that the ground state becomes

(4.86) \begin{equation} F_1 = \frac {9 n_0 v_{c,0}^2 C}{2 \pi \mathcal{E}(v_{c,0},v_\alpha ) } \left \{\frac {v \varTheta [v_\alpha - v]}{\big( v_{c,0}^3 + v^3 \big)^2} + \frac {\delta [v_\alpha - v]}{3v\big( v_{c,0}^3 + v^3 \big)} \right \}\! ,\end{equation}

having used $\langle f_1 \rangle = 0$ . One can now calculate the available energy, realising special care is required for $v=v_\alpha$ , as the integrand is dominated by the Dirac-delta peak here:

(4.87) \begin{align} \lim _{v\rightarrow v_\alpha }\frac {(f_1-F_1)^2}{2 \partial _{\epsilon }f_{s,0}} = \frac {6 m n_0 v_{c,0}^4 v_\alpha (C-\tilde {\nu })^2 \delta (0) \ln \big(1 + {v_\alpha ^3}/{v_{c,0}^3} \big)}{\big(v_{c,0}^3 + v_\alpha ^3\big) \mathcal{E}(v_{c,0},v_\alpha )^2}. \end{align}

Combining all results yields

(4.88) \begin{align} A & = \frac {n_0 m v_\alpha ^2 V}{4} \Bigg \langle \frac {2}{\mathcal{L}_1} \nu ^2 - \frac {8 p^2 \mathcal{L}_1}{\mathcal{E}} \tilde {\nu }^2 + \frac {{3}/{1 + p^3} -p^2 \mathcal{L}_1 (\mathcal{E}-2\mathcal{L}_2)}{\mathcal{L}_1^2} \nu \tau \nonumber\\ & \quad +\frac {(p^3+1) p^2 \mathcal{L}_1 \left [\mathcal{E} (p^3 \mathcal{L}_1+\mathcal{L}_1-6)-2 \mathcal{L}_2 (p^3 \mathcal{L}_1+\mathcal{L}_1-6 )+6 p \mathcal{L}_1\right ]+9}{2 (p^3+1)^2 \mathcal{L}_1^3} \tau ^2 \Bigg \rangle ,\nonumber\\[-10pt]\nonumber\\ \end{align}

having defined $p = v_{c,0}/v_\alpha$ , $\mathcal{L}_1 = \ln (1 + v_\alpha ^3/v_{c,0}^3)$ and $\mathcal{L}_2 = 2 \ln (v_{c,0} + v_\alpha ) - \ln ( v_{c,0}^2 + v_\alpha ^2 - v_{c,0} v_\alpha )$ . As before, it is interesting to keep the density profile fixed with no temperature variation, $\tilde {\nu } = \nu$ and $\tau = 0$ , yielding

(4.89) \begin{equation} A(\tilde {\nu }=\nu ,\tau =0) = \frac {n_0 m v_\alpha ^2 V}{4} \frac {\mathcal{E} - 4 p^2 \mathcal{L}_1^2}{\mathcal{E}\mathcal{L}_1} \langle \nu ^2 \rangle , \end{equation}

thus, energy is still available. A plot showing how $A$ varies with respect to the relative cross-over velocity $p=v_{c,0}/v_\alpha$ , as in (4.89), is shown in figure 9. A local minimum is found around $p\approx 1/3$ (correct to $\lt 1\,\%$ ) and, hence, this may be an optimal operating point for the given conditions.

Figure 9.

Plot of the dependence of the available energy with fixed density profile and no temperature gradients on the relative cross-over velocity $p = v_{c,0}/v_\alpha$ .

Figure 9. Long description

A line graph depicts the relationship between the available energy and the relative cross-over velocity. The x-axis represents the relative cross-over velocity on a logarithmic scale ranging from 10^-2 to 10^0. The y-axis represents the available energy, also on a logarithmic scale ranging from 10^-1 to 10^1. The graph shows a curve that starts low, rises slightly, dips, and then sharply increases as the relative cross-over velocity approaches 10^0.

4.9. A fast-ion species subject to quasineutrality

In this section multiple species are coupled via quasineutrality as in § 4.5, but with the impurities replaced by a slowing-down distribution of mean density $n_\alpha$ and with $Z=Z_\alpha$ , which may be taken to be a model for a fast-ion population. We note that formally the slowing-down distribution is derived under the assumption of small $n_\alpha$ compared with the main ion and electron densities, but we continue nonetheless to assess how quasineutrality alters the available energy. The condition for quasineutrality is simply (as in (4.67), already setting the variation in trapped electron density to zero as this maximises available energy)

(4.90) \begin{equation} f_p \nu _{e,p} = \frac {n_H}{n_0} \tilde {\nu }_H + Z_\alpha \frac {n_\alpha }{n_0} \tilde {\nu }_\alpha , \end{equation}

whereas the available energy becomes

(4.91) \begin{equation} A = A_{\mathrm{G}} - \frac {n_H T_H V}{2} \langle \tilde {\nu }_H^2 \rangle - \frac {n_\alpha T_\alpha V}{2} \frac {8 p^2 \mathcal{L}_1}{\mathcal{E}} \langle \tilde {\nu }^2_\alpha \rangle, \end{equation}

having defined $m_\alpha v_\alpha ^2/2 = T_\alpha$ . Expressing $\tilde {\nu }_\alpha$ in terms of $\tilde {\nu }_H$ using (4.80) gives the optimal profile

(4.92) \begin{equation} \tilde {\nu }_H = \frac {f_p n_e T_\alpha }{n_H T_\alpha + n_\alpha T_H ( {Z_\alpha ^2\mathcal{E}}/{8 p^2 \mathcal{L}_1})} \nu _{e,p}, \end{equation}

so that the available energy becomes

(4.93) \begin{equation} A = A_{\mathrm{G}} - \frac {f_p^2 n_e^2 T_H T_\alpha V}{2(n_HT_\alpha + n_\alpha T_H ({Z_\alpha ^2\mathcal{E}}/{8 p^2 \mathcal{L}_1}) )} \langle \nu _{e,p}^2 \rangle . \end{equation}

This is equivalent to (4.71) in § 4.5, where the neutral density and temperature is replaced with that of the $\alpha$ particles, and $Z^2 \rightarrow Z_\alpha ^2 \mathcal{E}/8p^2\mathcal{L}_1$ . A plot of the dependence of this factor on $p$ is given in figure 10, where we see that small values of $p$ lead to large $\mathcal{E}/8p^2\mathcal{L}_1$ and, thus, a smaller penalty due to quasineutrality. However, changing $p$ will, in general, also affect $A_{\mathrm{G}}$ (analogous to the discussion in § 4.5) and an optimal operating $p$ is not clear a priori.

Figure 10.

Plot of $\mathcal{E}/8p^2\mathcal{L}_1$ as a function of the relative cross-over velocity $p=v_{c,0}/v_\alpha$ .

4.10. Parallelly local ions in a flux tube

Let us now instead assume that there is an ion species that is parallelly local as in Mackenbach et al. (Reference Mackenbach, Helander, Landreman, Brunner and Proll2025), i.e. we assume that $k_\| v_T \ll \boldsymbol{k}_\perp \boldsymbol{\cdot }\boldsymbol{v}_D$ , with $v_T$ and $\boldsymbol{v}_D$ being the thermal velocity and the particle drifts and $k_\|$ and $\boldsymbol{k}_\perp$ being the parallel and perpendicular wavenumbers, respectively. We do so in the geometry of a flux tube: a slender domain around a given magnetic-field line that is everywhere parallel to the magnetic field, and we consequently take $\mu$ , $v_\|$ and $\ell$ to be conserved. This is worked out with the asymptotic framework in Appendix F, where we find a different result than in the aforementioned publication. The reason for this is the same as discussed in Helander & Mackenbach (Reference Helander and Mackenbach2024): the ground state breaks the presented asymptotic ordering when only one coordinate is left for the Gardner restacking, as this leads to sharp derivatives in the distribution function that we do not account for. As such, when calculating the ground state and available energy with the asymptotic framework considered in this paper, one finds the erroneous result of Appendix F.