1 Introduction
Magnetised plasmas are subject to a plethora of instabilities. Some of these are catastrophic in the sense that they lead to a sudden loss of confinement or plasma termination, as in tokamak disruptions (Schuller Reference Schuller1995) or solar flares (Fletcher et al. Reference Fletcher, Dennis, Hudson, Krucker, Phillips, Veronig, Battaglia, Bone, Caspi and Chen2011), but most instabilities are of a more benign nature. In particular those with short wavelengths (microinstabilities) usually saturate at a low level and merely cause turbulent fluctuations and degradation of the confinement due to transport (Garbet et al. Reference Garbet, Mantica, Angioni, Asp, Baranov, Bourdelle, Budny, Crisanti, Cordey and Garzotti2004).
The traditional way of mathematically investigating plasma stability is to identify an equilibrium state and then calculate the various linear eigenmodes and their growth rates. In simple cases, this can be done analytically, but accurate mode structures and growth rates can often only be found numerically. Great efforts have gone into the construction of computer codes for this purpose, see e.g. Kotschenreuther, Rewoldt & Tang (Reference Kotschenreuther, Rewoldt and Tang1995).
In reality, however, linear microinstabilities are rarely observed in laboratory and space plasmas. These are instead maintained in a nonlinearly fluctuating state, and the relevance of linear stability theory often seems questionable. The amplitude of the observed fluctuations do indeed seem to be related to the ‘drives’ identified by linear stability theory, such as density and temperature gradients or unfavourable magnetic curvature (Garbet et al. Reference Garbet, Idomura, Villard and Watanabe2010), but the assumptions made in linear stability theory are so restrictive that they rarely apply in practice. A question that arises, then, is whether there is some other way of characterising how ‘unstable’ a plasma is, without resorting to any linear approximation.
It is the purpose of the present paper to suggest a quantity that naturally lends itself to this purpose, which we shall call the ‘available energy’ of the plasma. It is related to the concept of ‘available potential energy’ in meteorology, which was introduced by Lorenz (Reference Lorenz1955) and denotes the excess of potential energy of the atmosphere above the minimum that could result from any adiabatic redistribution of mass. Analogously, we shall consider the excess of thermal energy of a magnetically confined plasma over the minimum accessible to it under the various constraints that limit the possible forms of plasma motion. This excess measures how much thermal energy is ‘available’ for conversion into linear and nonlinear plasma instabilities, and therefore constitutes an intuitive measure of plasma stability. Of particular interest are plasmas with zero available energy. Such a plasma is in its lowest possible energy state, the ‘ground state’, and is therefore stable to linear and nonlinear perturbations.
The available energy is in general much smaller than the total thermal energy, and depends on what type of instability is considered – or rather on what constraints are imposed on the motion of the plasma. For instance, if we consider collisionless, electrostatic instabilities with frequencies below the ion cyclotron frequency in a plasma embedded in a constant magnetic field of strength
$B$
, the magnetic moment
$\unicode[STIX]{x1D707}=mv_{\bot }^{2}/2B$
is conserved for all plasma particles, and the average perpendicular pressure
$$\begin{eqnarray}\langle p_{\bot }\rangle =\left\langle \int \unicode[STIX]{x1D707}Bfd^{3}v\right\rangle\end{eqnarray}$$
therefore cannot change. (Here
$f$
is the distribution function,
$m$
the mass,
$v_{\bot }$
the velocity perpendicular to the magnetic field, and angular brackets denote a volume average.) The part of the thermal energy that is associated with particle motion perpendicular to the magnetic field is therefore not available to drive such instabilities. This condition was recently used to derive upper bounds on magnetic field generation by dynamo action in collisionless plasmas (Helander, Strumik & Schekochihin Reference Helander, Strumik and Schekochihin2016).
Another important constraint was identified by Gardner (Reference Gardner1963) and applies if the distribution function is governed by the Vlasov equation,
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}(\sqrt{g}f)}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\sqrt{g}\dot{\boldsymbol{x}}f)=0,\end{eqnarray}$$
where
$\boldsymbol{x}$
denotes arbitrary phase-space coordinates,
$\unicode[STIX]{x1D735}=\unicode[STIX]{x2202}/\unicode[STIX]{x2202}\boldsymbol{x}$
and
$\sqrt{g}$
is the Jacobian, which satisfies Liouville’s theorem,
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}(\sqrt{g})}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\sqrt{g}\dot{\boldsymbol{x}})=0.\end{eqnarray}$$
Multiplying the Vlasov equation by
$G^{\prime }(f)$
and integrating over
$\boldsymbol{x}$
gives
$$\begin{eqnarray}\frac{\text{d}}{\text{d}t}\left\langle \int G(f)d^{3}v\right\rangle =0\end{eqnarray}$$
for any function
$G$
. There is thus an uncountable infinity of constraints on the evolution of
$f$
, and Gardner went on to determine the distribution function
$f_{0}(v)$
that minimises the energy under these constraints. The purpose of the present work is to extend his calculation by accounting for additional constraints and apply it to magnetically confined plasmas.
In the following section, we define the available energy for any kinetic plasma model that satisfies a Liouville theorem. In § 3 we calculate the available energy of a Vlasov plasma and show that it is larger than the corresponding quantity in ideal and resistive magnetohydrodynamics (MHD) in § 4. In the next two sections, we consider the effect of conservation laws such as magnetic-moment conservation, and find that the minimum-energy state in general is inhomogeneous unless the magnetic field is constant over the plasma volume. Finally, in § 7, we explore the effect of the second adiabatic invariant and show that its conservation implies that plasmas confined in certain magnetic configurations have particularly small available energy.
2 The ground state
Let us consider the evolution of the distribution function of some plasma species, assuming that it satisfies the Vlasov equation or any kinetic equation of the form (1.2) with the property that Liouville’s theorem (1.3) holds. The question we are asking is: what is the minimum energy
$$\begin{eqnarray}E=\int \unicode[STIX]{x1D716}f\sqrt{g}\,\text{d}\boldsymbol{x}\end{eqnarray}$$
accessible to the distribution function? Here
$\unicode[STIX]{x1D716}(\boldsymbol{x})$
denotes the energy (usually
$mv^{2}/2$
) of an individual particle at position
$\boldsymbol{x}$
in phase space. Clearly, the lowest possible value for
$E$
would be attained if the support of
$f$
(i.e. the region in which it assumes non-zero values) were concentrated to the point(s) where
$\unicode[STIX]{x1D716}(\boldsymbol{x})$
vanishes. The energy
$E$
would then also vanish, but this state is usually not accessible from the initial condition.
If we denote the energy-minimising distribution function (the ‘ground state’) by
$f_{0}(\boldsymbol{x})$
, then it should not be possible for the system to evolve away from this state in such a way that
$E$
decreases. Now, if
$f(\boldsymbol{x},t_{0})=f_{0}(\boldsymbol{x})$
, then for small
$\unicode[STIX]{x1D6FF}t$
we have
$$\begin{eqnarray}f(\boldsymbol{x},t_{0}+\unicode[STIX]{x1D6FF}t)-f(\boldsymbol{x},t_{0})\simeq \unicode[STIX]{x1D6FF}t\frac{\unicode[STIX]{x2202}f}{\unicode[STIX]{x2202}t}=\unicode[STIX]{x1D6FF}\boldsymbol{y}\boldsymbol{\cdot }\unicode[STIX]{x1D735}f_{0},\end{eqnarray}$$
where
$\unicode[STIX]{x1D6FF}\boldsymbol{y}=-\dot{\boldsymbol{x}}(\boldsymbol{x},t_{0})\unicode[STIX]{x1D6FF}t$
. The variation in energy is thus
$$\begin{eqnarray}\unicode[STIX]{x1D6FF}E=E(t_{0}+\unicode[STIX]{x1D6FF}t)-E(t_{0})=\int \unicode[STIX]{x1D716}(\unicode[STIX]{x1D6FF}\boldsymbol{y}\boldsymbol{\cdot }\unicode[STIX]{x1D735}f_{0})\sqrt{g}\,\text{d}\boldsymbol{x}+O(\unicode[STIX]{x1D6FF}t^{2}),\end{eqnarray}$$
and should vanish to first order in
$\unicode[STIX]{x1D6FF}t$
. This will clearly be the case if the functional
$\unicode[STIX]{x1D6FF}E[\unicode[STIX]{x1D6FF}\boldsymbol{y}]$
vanishes for all trial functions
$\unicode[STIX]{x1D6FF}\boldsymbol{y}$
such that
$\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FF}\boldsymbol{y}\sqrt{g})=0$
. We thus introduce a Lagrange multiplier
$\unicode[STIX]{x1D706}(\boldsymbol{x})$
and demand from
$f_{0}$
that it should make the functional
$$\begin{eqnarray}\unicode[STIX]{x1D6FF}W=\int [\unicode[STIX]{x1D716}(\unicode[STIX]{x1D6FF}\boldsymbol{y}\boldsymbol{\cdot }\unicode[STIX]{x1D735}f_{0})\sqrt{g}+\unicode[STIX]{x1D706}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FF}\boldsymbol{y}\sqrt{g})]\,\text{d}\boldsymbol{x},\end{eqnarray}$$
vanish for all functions
$\unicode[STIX]{x1D6FF}\boldsymbol{y}$
that vanish at infinity (or other boundaries of phase space). This is the case if, and only if,
$$\begin{eqnarray}\unicode[STIX]{x1D716}\unicode[STIX]{x1D735}f_{0}=\unicode[STIX]{x1D735}\unicode[STIX]{x1D706},\end{eqnarray}$$
which implies that
$f_{0}$
must be a function of
$\unicode[STIX]{x1D716}$
alone,
$f_{0}(\boldsymbol{x})=F_{0}[\unicode[STIX]{x1D716}(\boldsymbol{x})]$
. This can, for instance, be seen by differentiating the last equation, giving
$$\begin{eqnarray}(\unicode[STIX]{x2202}_{i}\unicode[STIX]{x1D716})(\unicode[STIX]{x2202}_{j}\hspace{0.6pt}f_{0})-(\unicode[STIX]{x2202}_{j}\unicode[STIX]{x1D716})(\unicode[STIX]{x2202}_{i}\hspace{0.6pt}f_{0})=0\end{eqnarray}$$
for all index pairs
$(i,j)$
, where we have written
$\unicode[STIX]{x2202}_{i}=\unicode[STIX]{x2202}/\unicode[STIX]{x2202}x_{i}$
.
To determine the function
$F_{0}(\unicode[STIX]{x1D716})$
, we appeal to (1.4) with
$G(f)=\unicode[STIX]{x1D6E9}(f-\unicode[STIX]{x1D719})$
, where
$\unicode[STIX]{x1D6E9}$
denotes the Heaviside step function and
$\unicode[STIX]{x1D719}$
an arbitrary constant, and we thus conclude that the function
$$\begin{eqnarray}H(\unicode[STIX]{x1D719})=\int \unicode[STIX]{x1D6E9}[f(\boldsymbol{x})-\unicode[STIX]{x1D719}]\sqrt{g}\,\text{d}\boldsymbol{x}\end{eqnarray}$$
must equal
$$\begin{eqnarray}\int \unicode[STIX]{x1D6E9}[F_{0}(\unicode[STIX]{x1D716}(\boldsymbol{x}))-\unicode[STIX]{x1D719}]\sqrt{g}\,\text{d}\boldsymbol{x}.\end{eqnarray}$$
This statement is just a reformulation of Liouville’s theorem: since the motion in phase space is incompressible, the volume in which the distribution exceeds a certain value,
$\unicode[STIX]{x1D719}$
, cannot change with time, whatever the choice of
$\unicode[STIX]{x1D719}$
. Furthermore, since
$F_{0}(\unicode[STIX]{x1D716})$
is a decreasing function of energy (see below), equation (2.8) must be equal to
$$\begin{eqnarray}\int \unicode[STIX]{x1D6E9}[\unicode[STIX]{x1D716}_{\unicode[STIX]{x1D719}}-\unicode[STIX]{x1D716}(\boldsymbol{x})]\sqrt{g}\,\text{d}\boldsymbol{x}=\unicode[STIX]{x1D6FA}(\unicode[STIX]{x1D716}_{\unicode[STIX]{x1D719}}),\end{eqnarray}$$
where
$\unicode[STIX]{x1D716}_{\unicode[STIX]{x1D719}}$
denotes the energy for which
$F_{0}(\unicode[STIX]{x1D716}_{\unicode[STIX]{x1D719}})=\unicode[STIX]{x1D719}$
. In other words,
$\unicode[STIX]{x1D716}_{\unicode[STIX]{x1D719}}$
is the inverse of the function
$F_{0}(\unicode[STIX]{x1D716})$
and
$$\begin{eqnarray}\unicode[STIX]{x1D6FA}^{\prime }(y)=\int \unicode[STIX]{x1D6FF}[\unicode[STIX]{x1D716}(\boldsymbol{x})-y]\sqrt{g}\,\text{d}\boldsymbol{x}\end{eqnarray}$$
is the ‘density of states’ of energy
$y$
. We thus conclude that the ground state is determined by the integral equation
$$\begin{eqnarray}H[F_{0}(\unicode[STIX]{x1D716})]=\unicode[STIX]{x1D6FA}(\unicode[STIX]{x1D716}),\end{eqnarray}$$
which is fundamental to what follows. Any monotonically decreasing function of energy,
$f(\boldsymbol{x})=F[\unicode[STIX]{x1D716}(\boldsymbol{x})]$
with
$\text{d}F/\text{d}\unicode[STIX]{x1D716}<0$
, is a possible ground state, because such a function trivially satisfies (2.11). If this equation is differentiated, we obtain an integro–differential equation
$$\begin{eqnarray}\frac{\text{d}F_{0}}{\text{d}\unicode[STIX]{x1D716}}=\frac{\unicode[STIX]{x1D6FA}^{\prime }(\unicode[STIX]{x1D716})}{H^{\prime }[F_{0}(\unicode[STIX]{x1D716})]},\end{eqnarray}$$
with the boundary condition
$F_{0}(\unicode[STIX]{x1D716})\rightarrow 0$
as
$\unicode[STIX]{x1D716}\rightarrow \infty$
. This equation has earlier been derived by Dodin & Fisch (Reference Dodin and Fisch2005) for the slightly more special case that the trajectories are determined by Hamiltonian equations of motion. We prefer to work with (2.11) since it represents the integral of the Dodin–Fisch equation.Footnote
1
Once either of these equations has been solved, the energy of the ground state can be calculated from
$$\begin{eqnarray}E_{0}=\int \unicode[STIX]{x1D716}F_{0}\sqrt{g}\,\text{d}\boldsymbol{x},\end{eqnarray}$$
and the available energy,
$A$
, defined to be the difference between the energy (2.1) of the initial state
$f$
and that of the ground state,
$$\begin{eqnarray}A=E-E_{0}.\end{eqnarray}$$
This is the maximum energy available for conversion into nonlinear fluctuations in the plasma. Strictly speaking, it is however only an upper bound on this energy, because we have not accounted for constraints other than that implied by Liouville’s theorem. For instance, we have not considered whether the final state, with kinetic energy in fluctuations, is actually accessible from the initial condition. Fisch & Rax (Reference Fisch and Rax1993) and Hay, Schiff & Fisch (Reference Hay, Schiff and Fisch2015) explored the maximum extractable energy under diffusion by waves, and found that it is generally less than that given by (2.14).
3 Available energy of Vlasov plasmas
Having defined the concepts of ground states and available energy, we now consider a couple of particularly simple examples where the evolution of the plasma is governed by the Vlasov equation with no additional constraints. The ground state will then always have constant density and temperature, since
$f_{0}$
is function of energy alone. Conversely, if the initial plasma density or temperature, defined by
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle n(\boldsymbol{r})=\int f(\boldsymbol{r},\boldsymbol{v})\,\text{d}\boldsymbol{v},\\ \displaystyle T(\boldsymbol{r})=\frac{2}{3n(\boldsymbol{r})}\int \unicode[STIX]{x1D716}(\boldsymbol{v})f(\boldsymbol{r},\boldsymbol{v})\,\text{d}\boldsymbol{v},\end{array}\right\}\end{eqnarray}$$
in Cartesian phase-space coordinates
$\boldsymbol{x}=(\boldsymbol{r},\boldsymbol{v})$
, varies over the plasma volume, then the plasma is not in a ground state and the available energy will in general be non-zero.
Since
$\unicode[STIX]{x1D716}=mv^{2}/2$
, the integrated density of states (2.9) is
$$\begin{eqnarray}\unicode[STIX]{x1D6FA}(y)=\int \text{d}\boldsymbol{r}\int _{0}^{\infty }\unicode[STIX]{x1D6E9}\left(y-\frac{mv^{2}}{2}\right)4\unicode[STIX]{x03C0}v^{2}\,\text{d}v=\frac{4\unicode[STIX]{x03C0}V}{3}\left(\frac{2y}{m}\right)^{3/2},\end{eqnarray}$$
where
$V$
is the volume of the spatial domain. This function plays a fundamental role in statistical mechanics and has recently stood in the centre of a debate on the correct definition of entropy (Dunkel & Hilbert Reference Dunkel and Hilbert2014).
3.1 Bi-Maxwellian initial condition
Our first example is a bi-Maxwellian distribution function, for simplicity taken to be constant in space,
$$\begin{eqnarray}f(\boldsymbol{v})=M\text{e}^{-mv_{\bot }^{2}/2T_{\bot }-mv_{\Vert }^{2}/2T_{\Vert }},\end{eqnarray}$$
where
$M=n(m/2\unicode[STIX]{x03C0}\bar{T})^{3/2}$
with
$\bar{T}=T_{\bot }^{2/3}T_{\Vert }^{1/3}$
. The energy density associated with this distribution function is
$$\begin{eqnarray}\frac{E}{V}=\frac{3nT}{2},\end{eqnarray}$$
where
$T=(2T_{\bot }+T_{\Vert })/3$
denotes the temperature defined in (3.1).
To calculate the ground state, we note that the surface
$f(\boldsymbol{v})=\unicode[STIX]{x1D719}$
in velocity space is given by
$$\begin{eqnarray}\frac{mv_{\bot }^{2}}{2T_{\bot }}+\frac{mv_{\Vert }^{2}}{2T_{\Vert }}=\ln \left(\frac{M}{\unicode[STIX]{x1D719}}\right),\end{eqnarray}$$
and represents an ellipsoid that encloses the volume
$$\begin{eqnarray}H(\unicode[STIX]{x1D719})=\frac{4\unicode[STIX]{x03C0}V}{3}\left(\frac{2\bar{T}}{m}\ln \frac{M}{\unicode[STIX]{x1D719}}\right)^{3/2}.\end{eqnarray}$$
The ground-state equation (2.11) thus becomes
$$\begin{eqnarray}\frac{4\unicode[STIX]{x03C0}V}{3}\left(\frac{2\bar{T}}{m}\ln \frac{M}{F_{0}}\right)^{3/2}=\frac{4\unicode[STIX]{x03C0}V}{3}\left(\frac{2\unicode[STIX]{x1D716}}{m}\right)^{3/2},\end{eqnarray}$$
where we have used (3.2), and the solution is simply a Maxwellian,
$$\begin{eqnarray}F_{0}(\unicode[STIX]{x1D716})=M\text{e}^{-\unicode[STIX]{x1D716}/\bar{T}}.\end{eqnarray}$$
The available energy is thus
$$\begin{eqnarray}A=\frac{3nV}{2}(T-\bar{T}),\end{eqnarray}$$
and is always positive, since
$T$
denotes the arithmetic mean between the temperatures
$(T_{\bot },T_{\bot },T_{\Vert })$
in the three directions of velocity space, which always exceeds the geometric mean
$\bar{T}$
.
This example clearly illustrates the difference between equilibrium states and ground states of the Vlasov equation. A spatially constant bi-Maxwellian is an equilibrium solution of this equation, but the ground state is Maxwellian, even in the absence of collisions.
Fowler (Reference Fowler, Simon and Thompson1968) used thermodynamic arguments, specifically the Helmholtz free energy, to derive bounds on the fluctuation amplitude in turbulent plasmas. Schekochihin (Reference Schekochihin2017) has shown that when his technique is applied to a bi-Maxwellian, the result is the same as (3.9).
3.2 Maxwellian initial condition
Our second example is an initially Maxwellian distribution function,
$$\begin{eqnarray}f_{M}(\boldsymbol{r},\boldsymbol{v})=n(\boldsymbol{r})\left(\frac{m}{2\unicode[STIX]{x03C0}T(\boldsymbol{r})}\right)^{3/2}\text{e}^{-mv^{2}/2T(\boldsymbol{r})},\end{eqnarray}$$
whose density and temperature are now allowed to vary with
$\boldsymbol{r}$
. The total thermal energy is
$$\begin{eqnarray}E_{M}=\frac{3V}{2}\langle nT\rangle ,\end{eqnarray}$$
where angular brackets again denote a spatial average. If we denote the maximum (over
$\boldsymbol{v}$
) of the Maxwellian (3.10) by
$$\begin{eqnarray}M(\boldsymbol{r})=f_{M}(\boldsymbol{r},0)=n(\boldsymbol{r})\left(\frac{m}{2\unicode[STIX]{x03C0}T(\boldsymbol{r})}\right)^{3/2},\end{eqnarray}$$
the function (2.7) becomes
$$\begin{eqnarray}H(\unicode[STIX]{x1D719})=\frac{4\unicode[STIX]{x03C0}}{3}\int \left(\frac{2T}{m}\ln \frac{M}{\unicode[STIX]{x1D719}}\right)^{3/2}\unicode[STIX]{x1D6E9}(M-\unicode[STIX]{x1D719})\,\text{d}\boldsymbol{r},\end{eqnarray}$$
and the ground-state equation (2.11)
$$\begin{eqnarray}\unicode[STIX]{x1D716}^{3/2}=\left\langle \left[T(\boldsymbol{r})\ln \frac{M(\boldsymbol{r})}{F_{0}(\unicode[STIX]{x1D716})}\right]^{3/2}\unicode[STIX]{x1D6E9}[M(\boldsymbol{r})-F_{0}(\unicode[STIX]{x1D716})]\right\rangle .\end{eqnarray}$$
Note that the ground state
$F_{0}(\unicode[STIX]{x1D716})$
determined by this nonlinear integral equation is generally not Maxwellian.
We are now in a position to calculate the ground-state energy
$$\begin{eqnarray}\frac{E_{0}}{V}=\int _{0}^{\infty }\unicode[STIX]{x1D716}F_{0}4\unicode[STIX]{x03C0}v^{2}\,\text{d}v=\frac{4\unicode[STIX]{x03C0}}{5}\left(\frac{2}{m}\right)^{3/2}\int _{0}^{maxF_{0}}\unicode[STIX]{x1D716}^{5/2}\,\text{d}F_{0},\end{eqnarray}$$
where we substitute the solution (3.14) and obtain
$$\begin{eqnarray}\frac{E_{0}}{V}=\frac{4\unicode[STIX]{x03C0}}{5}\left(\frac{2}{m}\right)^{3/2}\int _{0}^{maxF_{0}}\left\langle \left(T\ln \frac{M}{F_{0}}\right)^{3/2}\unicode[STIX]{x1D6E9}(M-F_{0})\right\rangle ^{5/3}\,\text{d}F_{0}.\end{eqnarray}$$
According to Jensen’s inequality, we have
$$\begin{eqnarray}\left\langle \left(T\ln \frac{M}{F_{0}}\right)^{3/2}\unicode[STIX]{x1D6E9}(M-F_{0})\right\rangle ^{5/3}\leqslant \left\langle \left(T\ln \frac{M}{F_{0}}\right)^{5/2}\unicode[STIX]{x1D6E9}(M-F_{0})\right\rangle ,\end{eqnarray}$$
and it follows that
$$\begin{eqnarray}E_{0}\leqslant E_{M},\end{eqnarray}$$
with equality if, and only if,
$n(\boldsymbol{r})$
and
$T(\boldsymbol{r})$
are constant over the domain in question. As expected, a Maxwellian (3.10) with spatially varying density or temperature thus has higher energy (3.11) than the ground state.
A general formula for the difference – the available energy – can be derived if the density and temperature only vary slightly, so that we can write
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}n(\boldsymbol{r})=\langle n\rangle [1+\unicode[STIX]{x1D708}(\boldsymbol{r})],\\ T(\boldsymbol{r})=\langle T\rangle [1+\unicode[STIX]{x1D70F}(\boldsymbol{r})],\end{array}\right\}\end{eqnarray}$$
with
$\unicode[STIX]{x1D708}\sim \unicode[STIX]{x1D70F}\ll 1$
. Then the energy of the Maxwellian (3.11) becomes
$$\begin{eqnarray}E_{M}=\frac{3V}{2}\langle n\rangle \langle T\rangle \langle 1+\unicode[STIX]{x1D708}\unicode[STIX]{x1D70F}\rangle ,\end{eqnarray}$$
and the function (3.12)
$$\begin{eqnarray}M=\bar{M}(1+\unicode[STIX]{x1D708})(1+\unicode[STIX]{x1D70F})^{-3/2},\end{eqnarray}$$
where
$$\begin{eqnarray}\bar{M}=\langle n\rangle \left(\frac{m}{2\unicode[STIX]{x03C0}\langle T\rangle }\right)^{3/2}.\end{eqnarray}$$
The ground-state energy (3.16) can now be expressed as
$$\begin{eqnarray}\frac{E_{0}}{V}=\frac{4}{5\sqrt{\unicode[STIX]{x03C0}}}\langle n\rangle \langle T\rangle \int _{x_{0}}^{\infty }\left\langle (1+\unicode[STIX]{x1D70F})^{3/2}\left[x+\ln \frac{1+\unicode[STIX]{x1D708}}{(1+\unicode[STIX]{x1D70F})^{3/2}}\right]^{3/2}\unicode[STIX]{x1D6E9}\right\rangle ^{5/3}\text{e}^{-x}\,\text{d}x,\end{eqnarray}$$
where
$x=\ln (\bar{M}/F_{0})$
and
$x_{0}=\ln (\bar{M}/\max F_{0})$
. Most of the contribution to the integral on the right-hand side comes from values of
$x$
much larger than
$\unicode[STIX]{x1D708}\sim \unicode[STIX]{x1D70F}\ll 1$
. The integral can therefore be calculated by expanding the integrand to second order in
$\unicode[STIX]{x1D708}$
and
$\unicode[STIX]{x1D70F}$
, and then integrating over
$x$
from 0 to
$\infty$
. Since
$$\begin{eqnarray}\ln \frac{1+\unicode[STIX]{x1D708}}{(1+\unicode[STIX]{x1D70F})^{3/2}}\simeq \unicode[STIX]{x1D708}-\frac{\unicode[STIX]{x1D708}^{2}}{2}-\frac{3\unicode[STIX]{x1D70F}}{2}+\frac{3\unicode[STIX]{x1D70F}^{2}}{4},\end{eqnarray}$$
and
$\langle \unicode[STIX]{x1D708}\rangle =\langle \unicode[STIX]{x1D70F}\rangle =0$
, we have
$$\begin{eqnarray}\displaystyle & & \displaystyle \left\langle (1+\unicode[STIX]{x1D70F})^{3/2}\left[x+\ln \frac{1+\unicode[STIX]{x1D708}}{(1+\unicode[STIX]{x1D70F})^{3/2}}\right]^{3/2}\right\rangle ^{5/3}\nonumber\\ \displaystyle & & \displaystyle \quad \simeq x^{5/2}+\frac{5x^{3/2}}{4}\langle 3\unicode[STIX]{x1D708}\unicode[STIX]{x1D70F}-\unicode[STIX]{x1D708}^{2}-3\unicode[STIX]{x1D70F}^{2}\rangle +\frac{5x^{1/2}}{8}\left\langle \unicode[STIX]{x1D708}^{2}-3\unicode[STIX]{x1D708}\unicode[STIX]{x1D70F}+\frac{9\unicode[STIX]{x1D70F}^{2}}{4}\right\rangle ,\end{eqnarray}$$
which can be inserted into the integral (3.23), giving
$$\begin{eqnarray}\frac{E_{0}}{V}=\frac{3V}{2}\langle n\rangle \langle T\rangle \left\langle 1-\frac{\unicode[STIX]{x1D708}^{2}}{3}-\frac{\unicode[STIX]{x1D70F}^{2}}{2}+\unicode[STIX]{x1D708}\unicode[STIX]{x1D70F}\right\rangle .\end{eqnarray}$$
Comparing with (3.20), we can thus express the available energy asFootnote 2
$$\begin{eqnarray}\frac{A}{E_{M}}=\left\langle \frac{\unicode[STIX]{x1D708}^{2}}{3}+\frac{\unicode[STIX]{x1D70F}^{2}}{2}\right\rangle .\end{eqnarray}$$
3.3 Several particle species
A simple consequence of (3.27) is that it implies a limit on the amount of energy that can be transferred between different kinds of particles in a multi-species plasma. A common situation in astrophysics is that one particle species, usually the ions, has a significantly higher temperature than the other(s). A legitimate question, then, is whether there is any mechanism, linear or nonlinear, by which energy can be transferred from the hotter to the colder species.
Equation (3.27) implies an essentially negative answer to this question. In the absence of collisions, no energy can be extracted from a Maxwellian (or any other decreasing function of energy) unless there are spatial density and temperature variations. Only the available energy associated with the latter is accessible for extraction and transfer to other species. The temperature difference between different species does not itself contribute to the available energy. On the other hand, density and temperature gradients may drive instabilities that indirectly lead to energy exchange between different species – even if their temperatures happen to be equal (Barnes, Abiuso & Dorland Reference Barnes, Abiuso and Dorland2017).
4 Comparison with MHD
It is interesting to compare the result (3.27) with the corresponding one obtained in magnetohydrodynamics by Helander et al. (Reference Helander, Strumik and Schekochihin2016). From the continuity equation and the entropy conservation law of ideal MHD,
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \frac{\text{d}\unicode[STIX]{x1D70C}}{\text{d}t}+\unicode[STIX]{x1D70C}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{V}=0,\\ \displaystyle \frac{\text{d}}{\text{d}t}\left(\frac{p}{\unicode[STIX]{x1D70C}^{\unicode[STIX]{x1D6FE}}}\right)=0,\end{array}\right\}\end{eqnarray}$$
where
$\unicode[STIX]{x1D70C}$
denotes the density,
$p=nT$
the pressure,
$\unicode[STIX]{x1D6FE}=5/3$
the adiabatic index and
$\boldsymbol{V}$
the fluid velocity, it is straightforward to derive the equation
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}p^{1/\unicode[STIX]{x1D6FE}}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(p^{1/\unicode[STIX]{x1D6FE}}\boldsymbol{V})=0.\end{eqnarray}$$
If there is no flow across the boundary, it follows that
$$\begin{eqnarray}S=\langle p^{1/\unicode[STIX]{x1D6FE}}\rangle\end{eqnarray}$$
is a conserved quantity, which is a manifestation of the entropy law in the present context. If other constraints (e.g. topological ones) are ignored, the lowest-energy state accessible to an MHD plasma is obtained by minimising the thermal energy
$$\begin{eqnarray}E=\frac{\langle p\rangle }{\unicode[STIX]{x1D6FE}-1}V,\end{eqnarray}$$
at fixed
$S$
. The resulting state is one of constant pressure
$p_{0}=\langle p^{1/\unicode[STIX]{x1D6FE}}\rangle ^{\unicode[STIX]{x1D6FE}}$
, and the available energy is thus given by (Helander et al.
Reference Helander, Strumik and Schekochihin2016)
$$\begin{eqnarray}A_{MHD}=\frac{\langle p\rangle -\langle p^{1/\unicode[STIX]{x1D6FE}}\rangle ^{\unicode[STIX]{x1D6FE}}}{\unicode[STIX]{x1D6FE}-1}V.\end{eqnarray}$$
If the initial pressure perturbations are small,
$p(\boldsymbol{r})=\langle p\rangle +\unicode[STIX]{x1D6FF}p(\boldsymbol{r})$
, the available energy becomes
$$\begin{eqnarray}A_{MHD}=\frac{\langle p\rangle }{2\unicode[STIX]{x1D6FE}}\left\langle \left(\frac{\unicode[STIX]{x1D6FF}p}{p}\right)^{2}\right\rangle V\end{eqnarray}$$
to second order in the fluctuations. By making the identification
$p=nT=\langle n\rangle \langle T\rangle (1+\unicode[STIX]{x1D708}+\unicode[STIX]{x1D70F}+\unicode[STIX]{x1D708}\unicode[STIX]{x1D70F})$
, we can compare this result with (3.27). The difference in available energies becomes
$$\begin{eqnarray}A-A_{MHD}=\frac{\langle p\rangle }{5}\left\langle \left(\unicode[STIX]{x1D708}-\frac{3\unicode[STIX]{x1D70F}}{2}\right)^{2}\right\rangle V\geqslant 0,\end{eqnarray}$$
showing that more energy is available in a Maxwellian Vlasov plasma than in MHD. This conclusion does not change if resistivity and viscosity are introduced into the latter. The entropy then increases, causing
$S$
to grow and
$A_{MHD}$
to drop. Note that
$A=A_{MHD}$
exactly when the relative fluctuations in density are everywhere 50 % larger than those in temperature,
$\unicode[STIX]{x1D708}=3\unicode[STIX]{x1D70F}/2$
, implying that the linearised entropy perturbation vanishes,
$$\begin{eqnarray}\unicode[STIX]{x1D6FF}\ln \left(\frac{p}{\unicode[STIX]{x1D70C}^{\unicode[STIX]{x1D6FE}}}\right)=\frac{\unicode[STIX]{x1D6FF}p}{p}-\frac{\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x1D70C}}=\unicode[STIX]{x1D70F}-\frac{2\unicode[STIX]{x1D708}}{3}=0.\end{eqnarray}$$
5 Conserved quantities
The derivation given in § 2 can easily be extended to cases where the motion is constrained to satisfy certain conservation laws. The ground states then become less trivial and the concept of available energy correspondingly more interesting. If some quantity such as the magnetic moment
$\unicode[STIX]{x1D707}(\boldsymbol{x})$
is conserved, then the trial function
$\unicode[STIX]{x1D6FF}\boldsymbol{y}$
in (2.4) must be required to satisfy
$\unicode[STIX]{x1D6FF}\boldsymbol{y}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D707}=0$
, which can be accounted for by using
$\unicode[STIX]{x1D707}$
as one of the phase-space coordinates. We thus write
$\boldsymbol{x}=(\boldsymbol{z},\unicode[STIX]{x1D707})$
and demand that the integral
$$\begin{eqnarray}\unicode[STIX]{x1D6FF}W[\unicode[STIX]{x1D6FF}\boldsymbol{y};\unicode[STIX]{x1D707}]=\int [\unicode[STIX]{x1D716}(\unicode[STIX]{x1D6FF}\boldsymbol{y}\boldsymbol{\cdot }\unicode[STIX]{x1D735}f_{0})\sqrt{g}+\unicode[STIX]{x1D706}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FF}\boldsymbol{y}\sqrt{g})]\,\text{d}\boldsymbol{z},\end{eqnarray}$$
should vanish for every value of
$\unicode[STIX]{x1D707}$
, which is to be held fixed in the integration. This condition implies that
$f_{0}(\boldsymbol{x})$
depends on the phase-space coordinates only through
$\unicode[STIX]{x1D716}$
and
$\unicode[STIX]{x1D707}$
, i.e. it must be a function of the form
$f_{0}(\boldsymbol{z},\unicode[STIX]{x1D707})=F_{0}[\unicode[STIX]{x1D716}(\boldsymbol{z},\unicode[STIX]{x1D707}),\unicode[STIX]{x1D707}]$
. In analogy with the constraint (1.4) we now have
$$\begin{eqnarray}\frac{\text{d}}{\text{d}t}\int G(f)\sqrt{g}\,\text{d}\boldsymbol{z}=0\end{eqnarray}$$
for all
$\unicode[STIX]{x1D707}$
, which is again to kept constant in the integral. As in (2.7)–(2.9), we define
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle H(\unicode[STIX]{x1D719},\unicode[STIX]{x1D707})=\int \unicode[STIX]{x1D6E9}[f(\boldsymbol{x})-\unicode[STIX]{x1D719}]\sqrt{g}\,\text{d}\boldsymbol{z},\\ \displaystyle \unicode[STIX]{x1D6FA}(w,\unicode[STIX]{x1D707})=\int \unicode[STIX]{x1D6E9}[w-\unicode[STIX]{x1D716}(\boldsymbol{x})]\sqrt{g}\,\text{d}\boldsymbol{z},\end{array}\right\}\end{eqnarray}$$
and conclude that the ground state
$F_{0}(\unicode[STIX]{x1D716},\unicode[STIX]{x1D707})$
is determined by the equation
$$\begin{eqnarray}H[F_{0}(\unicode[STIX]{x1D716},\unicode[STIX]{x1D707}),\unicode[STIX]{x1D707}]=\unicode[STIX]{x1D6FA}(\unicode[STIX]{x1D716},\unicode[STIX]{x1D707}).\end{eqnarray}$$
Of course, it is trivial to generalise these equations to the case of multiple conservation laws. One merely needs to replace the scalar
$\unicode[STIX]{x1D707}$
by a vector representing all the conserved quantities. Note that any distribution function of the form
$F_{0}(\unicode[STIX]{x1D716},\unicode[STIX]{x1D707})$
is a ground state if
$\unicode[STIX]{x2202}F_{0}/\unicode[STIX]{x2202}\unicode[STIX]{x1D716}\leqslant 0$
everywhere.
6 Conservation of the magnetic moment
We now consider the case where
$\unicode[STIX]{x1D707}$
specifically denotes the magnetic moment,
$\unicode[STIX]{x1D707}=mv_{\bot }^{2}/(2B)$
, and demonstrate that the density and temperature of the ground state will in general not be constant but depend on the strength of the magnetic field
$B(\boldsymbol{r})$
.
A particularly simple example is furnished by a two-dimensional plasma with a magnetic field
$\boldsymbol{B}=B(x,y)\unicode[STIX]{x1D735}z$
pointing in the ignorable direction and varying in the two other directions. Any function
$F_{0}(\unicode[STIX]{x1D707})$
of
$\unicode[STIX]{x1D707}$
alone is then a ground state, with density
$$\begin{eqnarray}n=\int _{0}^{\infty }F_{0}2\unicode[STIX]{x03C0}v_{\bot }\,\text{d}v_{\bot }=\frac{2\unicode[STIX]{x03C0}B}{m}\int _{0}^{\infty }F_{0}\,\text{d}\unicode[STIX]{x1D707},\end{eqnarray}$$
proportional to
$B(x,y)$
. In other words, there is no available energy associated with density variations proportional to
$B$
if the magnetic moment is conserved. As a consequence, if the plasma is initially in a state of constant density but varying magnetic field strength, turbulence may be expected to spontaneously cause the density to become non-uniform (Yankov & Nycander Reference Yankov and Nycander1997). The temperature is
$$\begin{eqnarray}T=\frac{1}{n}\int _{0}^{\infty }\frac{mv_{\bot }^{2}}{2}F_{0}2\unicode[STIX]{x03C0}v_{\bot }\,\text{d}v_{\bot }=\frac{2\unicode[STIX]{x03C0}B^{2}}{mn}\int _{0}^{\infty }\unicode[STIX]{x1D707}F_{0}\,\text{d}\unicode[STIX]{x1D707},\end{eqnarray}$$
and will then also tend to a state proportional to
$B(x,y)$
.
In three dimensions, any function
$F_{0}(\unicode[STIX]{x1D716},\unicode[STIX]{x1D707})$
of energy
$\unicode[STIX]{x1D716}=mv^{2}/2$
and magnetic moment is a ground state if
$\unicode[STIX]{x2202}F_{0}/\unicode[STIX]{x2202}\unicode[STIX]{x1D716}\leqslant 0$
(Taylor Reference Taylor1963). A simple example is the spatially homogeneous bi-Maxwellian (3.3), which is a ground state if the magnetic field strength is constant but, as we have seen, is not a ground state if the magnetic moment is allowed to vary. The conservation of
$\unicode[STIX]{x1D707}$
thus prevents access to the lower-energy state (3.8) otherwise available to the plasma.
More generally, the density and temperature of a ground state
$F_{0}(\unicode[STIX]{x1D716},\unicode[STIX]{x1D707})$
,
$$\begin{eqnarray}\binom{n}{T}=\unicode[STIX]{x03C0}B\left(\frac{2}{m}\right)^{3/2}\int _{0}^{\infty }\left(\begin{array}{@{}c@{}}1\\ {\displaystyle \frac{2\unicode[STIX]{x1D716}}{3n}}\end{array}\right)\text{d}\unicode[STIX]{x1D716}\int _{0}^{\unicode[STIX]{x1D716}/B}\frac{F_{0}(\unicode[STIX]{x1D716},\unicode[STIX]{x1D707})\,\text{d}\unicode[STIX]{x1D707}}{\sqrt{\unicode[STIX]{x1D716}-\unicode[STIX]{x1D707}B}},\end{eqnarray}$$
will be non-uniform (unless
$\unicode[STIX]{x2202}F_{0}/\unicode[STIX]{x2202}\unicode[STIX]{x1D707}=0$
) if the magnetic field is inhomogeneous.
7 The adiabatic invariant of parallel motion
We now turn our attention to plasmas where the frequencies of collisions and any fluctuations are lower than that of particle motion along the magnetic field, so that not only the magnetic moment, but also the parallel adiabatic invariant,
$$\begin{eqnarray}J=\int mv_{\Vert }\,\text{d}l,\end{eqnarray}$$
is conserved. Here
$l$
denotes the arc length along
$\boldsymbol{B}$
, and the integration is carried out between two consecutive turning points of the motion. (
$J$
is usually conserved for electrons in hot fusion plasmas, since the frequency of the turbulence is comparable to the diamagnetic frequency, which is smaller than that of the electron transit frequency.) Any distribution function of the form
$F_{0}(\unicode[STIX]{x1D716},\unicode[STIX]{x1D707},J)$
is now a ground state if
$$\begin{eqnarray}\left(\frac{\unicode[STIX]{x2202}F_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D716}}\right)_{\unicode[STIX]{x1D707},J}\leqslant 0\end{eqnarray}$$
for every choice of
$\unicode[STIX]{x1D707}$
and
$J$
. This condition coincides with the linear stability criterion derived by Taylor (Reference Taylor1964) for electrostatic flute modes in mirror machines in the zero-gyroradius limit, but is of much wider significance (Schmidt Reference Schmidt1965; Taylor & Newton Reference Taylor and Newton2015). If this criterion is satisfied for all particle species, the plasma is not only linearly stable, but is in a global lowest-energy state and is therefore also nonlinearly stable. A detailed argument using gyrokinetics is given in the appendix A.
In a tokamak or stellarator, it is useful to write the magnetic field as
$\boldsymbol{B}=\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}$
, where
$\unicode[STIX]{x1D713}$
labels the toroidally nested flux surfaces and
$\unicode[STIX]{x1D6FC}$
the different field lines on each such surface.Footnote
3
The distribution function of each species is usually of the form of a Maxwellian whose density and temperature are constant on flux surfaces,
$$\begin{eqnarray}f_{M}(\unicode[STIX]{x1D713},\unicode[STIX]{x1D716})=n(\unicode[STIX]{x1D713})\left[\frac{m}{2\unicode[STIX]{x03C0}T(\unicode[STIX]{x1D713})}\right]^{3/2}\text{e}^{-\unicode[STIX]{x1D716}/T(\unicode[STIX]{x1D713})}.\end{eqnarray}$$
Expressing
$\unicode[STIX]{x1D713}$
as a function
$\unicode[STIX]{x1D713}=\unicode[STIX]{x1D713}(\unicode[STIX]{x1D707},J,\unicode[STIX]{x1D716})$
, we find
$$\begin{eqnarray}\left(\frac{\unicode[STIX]{x2202}f_{M}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D716}}\right)_{\unicode[STIX]{x1D707},J}=\left(\frac{\unicode[STIX]{x2202}f_{M}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}\right)_{\unicode[STIX]{x1D716}}\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D716}}\right)_{\unicode[STIX]{x1D707},J}-\frac{f_{M}}{T},\end{eqnarray}$$
where (Helander Reference Helander2014b )
$$\begin{eqnarray}\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D716}}\right)_{\unicode[STIX]{x1D707},J}=\left.-\left(\frac{\unicode[STIX]{x2202}J}{\unicode[STIX]{x2202}\unicode[STIX]{x1D716}}\right)_{\unicode[STIX]{x1D707},\unicode[STIX]{x1D713}}\right/\left(\frac{\unicode[STIX]{x2202}J}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}\right)_{\unicode[STIX]{x1D707},\unicode[STIX]{x1D716}}=\frac{1}{q\unicode[STIX]{x1D714}_{\unicode[STIX]{x1D6FC}}},\end{eqnarray}$$
$q$
denotes the electric charge and
$\unicode[STIX]{x1D714}_{\unicode[STIX]{x1D6FC}}=\overline{\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}}$
the frequency of trapped-particle precession in the
$\unicode[STIX]{x1D6FC}$
direction. Here
$\boldsymbol{v}_{d}$
is the drift velocity and an overbar indicates an orbit average. Thus, in standard gyrokinetic notation,
$$\begin{eqnarray}\left(\frac{\unicode[STIX]{x2202}f_{M}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D716}}\right)_{\unicode[STIX]{x1D707},J}=\frac{f_{M}}{T}\left(\frac{\unicode[STIX]{x1D714}_{\ast }^{T}}{\overline{\unicode[STIX]{x1D714}}_{d}}-1\right),\end{eqnarray}$$
where
$\unicode[STIX]{x1D714}_{d}=\boldsymbol{k}_{\bot }\boldsymbol{\cdot }\boldsymbol{v}_{d}$
with
$\boldsymbol{k}_{\bot }=k_{\unicode[STIX]{x1D6FC}}\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}$
,
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D714}_{\ast }=\frac{k_{\unicode[STIX]{x1D6FC}}T}{q}\frac{\text{d}\ln n}{\text{d}\unicode[STIX]{x1D713}},\\ \displaystyle \unicode[STIX]{x1D714}_{\ast }^{T}=\unicode[STIX]{x1D714}_{\ast }\left[1+\unicode[STIX]{x1D702}\left(\frac{\unicode[STIX]{x1D716}}{T}-\frac{3}{2}\right)\right],\end{array}\right\}\end{eqnarray}$$
and
$\unicode[STIX]{x1D702}=(\text{d}\ln T/\text{d}\unicode[STIX]{x1D713})/(\text{d}\ln n/\text{d}\unicode[STIX]{x1D713})$
. (The wavenumber
$k_{\unicode[STIX]{x1D6FC}}$
drops out and is only introduced for reasons of convention.) According to (7.6), the Maxwellian (7.3) is a ground state if
$$\begin{eqnarray}\frac{\unicode[STIX]{x1D714}_{\ast }^{T}}{\overline{\unicode[STIX]{x1D714}}_{d}}<1\end{eqnarray}$$
for all values of
$\unicode[STIX]{x1D716}$
,
$\unicode[STIX]{x1D707}$
and
$J$
. As we shall see in the next two subsections, there are two different ways to satisfy this condition.
7.1 Maximum-
$J$
devices
The criterion (7.8) holds if the following two conditions are both satisfied:
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D714}_{\ast }\overline{\unicode[STIX]{x1D714}}_{d}\leqslant 0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & 0\leqslant \unicode[STIX]{x1D702}<{\textstyle \frac{2}{3}}. & \displaystyle\end{eqnarray}$$
The first of these conditions is the maximum-
$J$
criterion,
$\unicode[STIX]{x2202}J/\unicode[STIX]{x2202}\unicode[STIX]{x1D713}<0$
if
$\text{d}n/\text{d}\unicode[STIX]{x1D713}<0$
, which was identified by Rosenbluth (Reference Rosenbluth1968) as stabilising to low-frequency interchange modes. It has recently been applied to gyrokinetic instabilities in quasi-isodynamic stellarators by Proll et al. (Reference Proll, Helander, Connor and Plunk2012) and Helander, Proll & Plunk (Reference Helander, Proll and Plunk2013), who showed that the two criteria (7.9) and (7.10) imply linear stability of collisionless trapped-particle modes. The physical significance of (7.9) is that it guarantees that all trapped particles experience favourable magnetic curvature on a time average along their orbits.
This is never the case in tokamaks. The magnetic curvature is unfavourable on the outboard side of the torus, where deeply trapped particles reside, whereas more shallowly trapped particles spend most of their time on the inboard side of the torus, where they experience good curvature. The drift frequency
$\overline{\unicode[STIX]{x1D714}}_{d}$
is thus positive for some orbits and negative for others, so the distribution function (7.3) cannot possibly correspond to a ground state.
The situation is different in certain stellarators, where the regions of particle trapping and bad magnetic curvature are separated from each other. This tends to be the case, at least to some approximation, in quasi-isodynamic stellarators at high beta (normalised plasma pressure). Thanks to the diamagnetic property of the plasma, its pressure digs a magnetic well surrounding the magnetic axis, causing the poloidal precession caused by the grad-
$B$
drift to reverse relative to the diamagnetic frequency, so that
$\unicode[STIX]{x1D714}_{\ast }\overline{\unicode[STIX]{x1D714}}_{d}$
becomes negative for all trapped orbits. If
$0\leqslant \unicode[STIX]{x1D702}<2/3$
, the plasma is then in a ground state relative to fluctuations that conserve
$\unicode[STIX]{x1D707}$
and
$J$
. This tends to be the case for the plasma electrons, since the typical turbulence frequency is of order
$\unicode[STIX]{x1D714}_{\ast }$
, which is much smaller than the electron bounce frequency if with
$k_{\bot }$
is comparable to the inverse ion gyroradius (Isichenko et al.
Reference Isichenko, Gruzinov, Diamond and Yushmanov1996). For the ions, however, the bounce frequency does not exceed the turbulence frequency, and
$J$
is not expected to be conserved.
An interesting question that arises, then, is the nature of the lowest-energy state in a plasma where
$J$
is conserved for the electrons but not for the ions. Different constraints then apply to the different particle species, which are however coupled to each other by the requirement of quasi-neutrality. The ground state of such a plasma should be sought by minimising the energy of both species simultaneously, each subject to its own particular constraints and additionally to that of quasi-neutrality. This problem is beyond the scope of the present publication.
7.2 Unfavourable magnetic curvature
Remarkably, there is another way of satisfying the stability criterion (7.8) if the magnetic curvature is unfavourable for all orbits. The bounce-averaged drift frequency is proportional to the kinetic energy,
$$\begin{eqnarray}\overline{\unicode[STIX]{x1D714}}_{d}=\tilde{\unicode[STIX]{x1D714}}_{d}(\unicode[STIX]{x1D706},\unicode[STIX]{x1D713},\unicode[STIX]{x1D6FC})\frac{\unicode[STIX]{x1D716}}{T(\unicode[STIX]{x1D713})},\end{eqnarray}$$
where
$\unicode[STIX]{x1D706}=\unicode[STIX]{x1D707}/\unicode[STIX]{x1D716}$
, and if
$\tilde{\unicode[STIX]{x1D714}}_{d}\unicode[STIX]{x1D714}_{\ast }$
is positive for all trapped orbits, so that all such trajectories on average experience bad magnetic curvature, then (7.8) is satisfied if simultaneously
$$\begin{eqnarray}\unicode[STIX]{x1D702}>{\textstyle \frac{2}{3}}\end{eqnarray}$$
and
$$\begin{eqnarray}\frac{\unicode[STIX]{x1D702}\unicode[STIX]{x1D714}_{\ast }}{\tilde{\unicode[STIX]{x1D714}}_{d}}<1\end{eqnarray}$$
for all particles. It is impossible to satisfy the second of these criteria in a tokamak, because
$\tilde{\unicode[STIX]{x1D714}}_{d}$
has different signs for deeply and shallowly trapped orbits. It therefore goes through zero, making
$\unicode[STIX]{x1D702}\unicode[STIX]{x1D714}_{\ast }/\tilde{\unicode[STIX]{x1D714}}_{d}$
arbitrarily large.
In a dipole magnetic field, however, the magnetic curvature is bad everywhere, and the precession frequency is an increasing function of
$\unicode[STIX]{x1D706}$
(Kesner & Hastie Reference Kesner and Hastie2002), bounded from below by
$$\begin{eqnarray}\tilde{\unicode[STIX]{x1D714}}_{d}(\unicode[STIX]{x1D706})\geqslant \tilde{\unicode[STIX]{x1D714}}_{d}(\unicode[STIX]{x1D706}=0).\end{eqnarray}$$
To calculate this bound, we note that for
$\unicode[STIX]{x1D706}=0$
,
$$\begin{eqnarray}J=mvL(\unicode[STIX]{x1D713}),\end{eqnarray}$$
where
$L(\unicode[STIX]{x1D713})$
denotes the length of the field lines on the flux surface labelled by
$\unicode[STIX]{x1D713}$
. For a point dipole,
$\boldsymbol{B}=\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}$
, where
$\unicode[STIX]{x1D6FC}$
denotes the toroidal angle, we have
$$\begin{eqnarray}\unicode[STIX]{x1D713}=\frac{M}{r}\sin ^{2}\unicode[STIX]{x1D703}\end{eqnarray}$$
in spherical coordinates. It follows that
$L(\unicode[STIX]{x1D713})$
is inversely proportional to
$\unicode[STIX]{x1D713}$
and thus
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \frac{\unicode[STIX]{x2202}J}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}=-\frac{J}{\unicode[STIX]{x1D713}},\\ \displaystyle \frac{\unicode[STIX]{x2202}J}{\unicode[STIX]{x2202}\unicode[STIX]{x1D716}}=\frac{J}{2\unicode[STIX]{x1D716}},\end{array}\right\}\end{eqnarray}$$
so that the drift frequency from (7.5) becomes
$$\begin{eqnarray}\unicode[STIX]{x1D714}_{d}(\unicode[STIX]{x1D706}=0)=\frac{2k_{\unicode[STIX]{x1D6FC}}\unicode[STIX]{x1D716}}{q\unicode[STIX]{x1D713}}.\end{eqnarray}$$
According to (7.12) and (7.13), the plasma is thus in a ground state if
$$\begin{eqnarray}\frac{2}{3}\frac{\text{d}\ln n}{\text{d}\ln \unicode[STIX]{x1D713}}<\frac{\text{d}\ln T}{\text{d}\ln \unicode[STIX]{x1D713}}<2,\end{eqnarray}$$
and nonlinear stability is then guaranteed. This region in parameter space occupies a subset of the linear stability domain identified in earlier works (Kesner & Hastie Reference Kesner and Hastie2002; Helander Reference Helander2014a ; Helander & Connor Reference Helander and Connor2016), see figure 1.
Stability diagram for low-frequency electrostatic perturbations of a collisionless plasma in the magnetic field of a point dipole. The linear stability region was calculated under the simplifying assumption that the precession frequency is independent of the trapping parameter
$\unicode[STIX]{x1D706}$
, which is not quite true but was shown by Kesner & Hastie (Reference Kesner and Hastie2002) to be a good approximation. The nonlinear stability domain is calculated without this approximation.
7.3 Partly unfavourable curvature
It is, of course, impossible to realise a pure dipole magnetic field, and we therefore briefly turn our attention to the field produced by a circular coil of finite radius. In this case, the two ways of satisfying the ground-state criterion (7.8) we have just discussed can be combined in such a way that (7.9)–(7.10) hold in some parts of the plasma and (7.12)–(7.13) elsewhere. If
$I$
denotes the current and
$a$
the radius of the coil, the magnetic flux function becomes (Jackson Reference Jackson1975)
$$\begin{eqnarray}\unicode[STIX]{x1D713}(r,\unicode[STIX]{x1D703})=\frac{\unicode[STIX]{x1D707}_{0}Iar\sin \unicode[STIX]{x1D703}}{\unicode[STIX]{x03C0}\sqrt{a^{2}+r^{2}+2ar\sin \unicode[STIX]{x1D703}}}\frac{(2-k^{2})K(k)-2E(k)}{k^{2}},\end{eqnarray}$$
where
$K$
and
$E$
denote complete elliptic integrals of the argument
$$\begin{eqnarray}k^{2}=\frac{4ar\sin \unicode[STIX]{x1D703}}{a^{2}+r^{2}+2ar\sin \unicode[STIX]{x1D703}}.\end{eqnarray}$$
The magnetic drift is everywhere in the same toroidal direction, so
$\unicode[STIX]{x1D714}_{d}$
never changes sign. The plasma density and temperature, however, should go to zero both close to the coil, where
$\unicode[STIX]{x1D713}$
becomes large, and far away from it, where
$\unicode[STIX]{x1D713}\rightarrow 0$
. The density and temperature will therefore peak on some intermediate surface, at
$\unicode[STIX]{x1D713}=\unicode[STIX]{x1D713}_{0}$
say, where the diamagnetic frequency thus changes sign, see figure 2. The product
$\tilde{\unicode[STIX]{x1D714}}_{d}\unicode[STIX]{x1D714}_{\ast }^{T}$
is negative inside this surface (
$\unicode[STIX]{x1D713}>\unicode[STIX]{x1D713}_{0}$
), reflecting favourable magnetic curvature, whereas
$\tilde{\unicode[STIX]{x1D714}}_{d}\unicode[STIX]{x1D714}_{\ast }^{T}>0$
and the curvature is unfavourable in the outer region (
$\unicode[STIX]{x1D713}<\unicode[STIX]{x1D713}_{0}$
). The ground-state condition (7.8) can therefore be satisfied throughout the entire plasma by choosing the density and temperature profiles so that
$0<\unicode[STIX]{x1D702}<2/3$
in the inner region and
$\unicode[STIX]{x1D702}>2/3$
in the outer region. In the latter region, the temperature profile should also satisfy the condition
$\unicode[STIX]{x1D702}\unicode[STIX]{x1D714}_{\ast }<\min \tilde{\unicode[STIX]{x1D714}}_{d}$
, which far from the coil (
$r\gg a$
) reduces to (7.19). It thus appears possible to use a single (levitated) coil to confine a collisionless plasma in such a way that it is nonlinearly stable to all low-frequency electrostatic perturbations.
Work is currently underway to construct a such a device for confining electron–positron plasmas (T. S. Pedersen, private communication (2017)). The density will be so low that the Debye length exceeds the gyroradius by a large factor, and for this reason modes with frequencies above the bounce frequency are predicted to be linearly stable (Helander Reference Helander2014a ). Given this linear stability at high frequencies and the nonlinear stability at lower ones, there is every reason to hope for good confinement.
Excellent plasma performance in dipole magnetic fields has indeed been observed experimentally (Boxer et al. Reference Boxer, Bergmann, Ellsworth, Garnier, Kesner, Mauel and Woskov2010; Yoshida et al. Reference Yoshida, Saitoh, Morikawa, Yano, Watanabe and Ogawa2010) and predicted theoretically on the basis of linear stability calculations (Simakov, Hastie & Catto Reference Simakov, Hastie and Catto2002) or relaxation arguments (Hasegawa, Chen & Mauel Reference Hasegawa, Chen and Mauel1990). The latter, however, imply non-Maxwellian distribution functions that are unlikely to apply if the confinement time exceeds the collision time.
Field lines produced by a circular coil at
$R=1$
,
$z=0$
. If the density and temperature peak on the field line labelled by
$\unicode[STIX]{x1D713}=\unicode[STIX]{x1D713}_{0}$
, the magnetic curvature is favourable inside this line and unfavourable outside it. The plasma is in a ground state if
$0<\unicode[STIX]{x1D702}<2/3$
in the inner region and the gradients in the outer region lie in the nonlinearly stable region corresponding to figure 1.
8 Conclusions
We have defined a ground state of a plasma species as a state with the lowest possible energy subject to whatever constraints limit the possible motion. The available energy is the difference between the actual energy of the plasma and the ground-state energy. It provides an upper bound on how much thermal energy can be converted into instabilities or turbulence.
If no constraint other than the Liouville theorem is imposed, the ground state of a Vlasov plasma is homogeneous, and the available energy of an initially Maxwellian plasma is proportional to the mean square of the density and temperature fluctuations, if these are assumed to be small. This energy exceeds the corresponding quantity calculated from ideal or visco-resistive MHD. The ground state of a spatially constant bi-Maxwellian distribution function is Maxwellian, even in the absence of collisions.
If the first or second adiabatic invariant (
$\unicode[STIX]{x1D707}$
or
$J$
) is conserved, the plasma parameters of most ground states vary over the plasma volume. If instabilities or turbulence bring the plasma toward such a state, spontaneous density and temperature peaking may thus result.
In maximum-
$J$
devices such as quasi-isodynamic stellarators, the average magnetic curvature is favourable for all particle trajectories and the plasma is in a ground state with respect to fluctuations that conserve
$\unicode[STIX]{x1D707}$
and
$J$
if the temperature profile is less than two thirds as steep as the density profile. The plasma is thus stable, even nonlinearly, to electrostatic instabilities with frequencies smaller than the bounce frequency of all species. Collisionless trapped-ion modes fall into this category, but the more common trapped-electron modes have frequencies that are comparable to the ion bounce frequency, so that
$J$
-conservation is broken for the ions. Since it still holds for the electrons, these act stabilising and any instability must draw its energy from the ions. This result has earlier been found through linear stability analysis by Proll et al. (Reference Proll, Helander, Connor and Plunk2012), but is here found to be true even for large perturbations.
Perhaps more surprisingly, the plasma can be in a ground state even in magnetic fields with unfavourable curvature, if the temperature profile is more than two thirds as steep as the density profile. For instance, in the field of a magnetic dipole, nonlinear stability to low-frequency electrostatic modes should prevail if (7.19) is satisfied. As a result, it appears that a low-density plasma could be confined extremely well by the field from a single levitated coil.
We have limited our discussion to collisionless and electrostatic instabilities, but it would be interesting also to explore the concepts of available energy and ground states in more general settings.
Acknowledgements
I would like to thank G. Plunk, N. Fisch, A. Schekochihin and I. Abel for many insightful comments and interesting discussions.
Appendix A. Nonlinear stability
On several occasions in this paper, it has been stated that ground states are stable, not only linearly but also to finite-sized perturbations. The precise meaning of this statement is that if the plasma is prepared in such a way that all distribution functions initially are in ground states, then the energy of the electric field fluctuations, suitably defined, cannot grow with time, independently of their initial amplitude. In order to prove this assertion, it is necessary to consider the full system of equations describing the electric field as well as all plasma particle species. We shall do so for the Vlasov–Poisson and gyrokinetic systems of equations, respectively. In both cases, it is sufficient to show that a quantity of the form
$$\begin{eqnarray}W=E+\unicode[STIX]{x1D6F7}\end{eqnarray}$$
is conserved, where
$\unicode[STIX]{x1D6F7}$
denotes the energy of fluctuations, which must be required to be positive definite, and
$E$
denotes the sum of all particle energies, taking the form
$$\begin{eqnarray}E=\mathop{\sum }_{a}\int \unicode[STIX]{x1D716}f_{a}\sqrt{g}\,\text{d}\boldsymbol{x}\end{eqnarray}$$
for some suitable function
$\unicode[STIX]{x1D716}(\boldsymbol{x})$
. It then follows that
$$\begin{eqnarray}\unicode[STIX]{x1D6F7}=W-E\leqslant W-E_{0}\end{eqnarray}$$
if
$E\geqslant E_{0}$
. Thus, if the system is initiated in a state with total energy
$W$
and the particles in a ground state, the fluctuation energy
$\unicode[STIX]{x1D6F7}$
cannot grow.
A.1 The Vlasov–Poisson system
The Vlasov–Poisson system of equations for an arbitrary number of species (distinguished by a subscript
$a$
) in an electric field
$\boldsymbol{E}=-\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}$
, is
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \frac{\unicode[STIX]{x2202}f_{a}}{\unicode[STIX]{x2202}t}+\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}f_{a}-\frac{q_{a}}{m_{a}}\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}f_{a}}{\unicode[STIX]{x2202}\boldsymbol{v}}=0,\\ \displaystyle \unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D719}=-4\unicode[STIX]{x03C0}\mathop{\sum }_{a}q_{a}\int f_{a}\,\text{d}\boldsymbol{v}.\end{array}\right\}\end{eqnarray}$$
If there is no flux to infinity or across other boundaries, this system of equations conserves the energy
$$\begin{eqnarray}W=E+\unicode[STIX]{x1D6F7},\end{eqnarray}$$
where
$$\begin{eqnarray}E=\mathop{\sum }_{a}\left\langle \int \frac{m_{a}v^{2}}{2}f_{a}\,\text{d}\boldsymbol{v}\right\rangle\end{eqnarray}$$
denotes the energy (2.1) of the particles, and
$$\begin{eqnarray}\unicode[STIX]{x1D6F7}=\left\langle \frac{|\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}|^{2}}{8\unicode[STIX]{x03C0}}\right\rangle\end{eqnarray}$$
that of the electric field. Since both energies are positive definite, the stability criterion formulated above is satisfied.
A.2 Gyrokinetic system of equations
We now demonstrate nonlinear stability of a plasma with respect to low-frequency fluctuations if the gyro-averaged distribution functions of all species satisfy
$$\begin{eqnarray}\left(\frac{\unicode[STIX]{x2202}F_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D716}}\right)_{\unicode[STIX]{x1D707},J}\leqslant 0\end{eqnarray}$$
everywhere in phase space. To this end, we note that the nonlinear, electrostatic system of gyrokinetic equations conserves the energy (Dubin et al. Reference Dubin, Krommes, Oberman and Lee1983)
$$\begin{eqnarray}W=E+\unicode[STIX]{x1D6F7},\end{eqnarray}$$
where
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle E=\mathop{\sum }_{a}\int f_{a}(\boldsymbol{R},\unicode[STIX]{x1D707},v_{\Vert })\unicode[STIX]{x1D705}\sqrt{g}\,\text{d}\boldsymbol{x},\\ \displaystyle \unicode[STIX]{x1D705}=\unicode[STIX]{x1D707}B+\frac{m_{a}v_{\Vert }^{2}}{2},\end{array}\right\}\end{eqnarray}$$
and
$$\begin{eqnarray}\unicode[STIX]{x1D6F7}=\int \left(\frac{|\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}|^{2}}{8\unicode[STIX]{x03C0}}+\mathop{\sum }_{a}\frac{m_{a}n_{a}}{2}\left|\frac{\boldsymbol{B}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D719}}{B^{2}}\right|^{2}\right)\text{d}\boldsymbol{r},\end{eqnarray}$$
denotes the fluctuation energy in the long wavelength limit, which is appropriate for low-frequency instabilities. Here
$\boldsymbol{R}$
denotes the gyro-centre position and
$\sqrt{g}\,\text{d}\boldsymbol{x}=(2\unicode[STIX]{x03C0}B_{\Vert }^{\ast }/m_{a})\,\text{d}\boldsymbol{R}\,\text{d}v_{\Vert }\,\text{d}\unicode[STIX]{x1D707}$
the phase-space volume element, which can also be written as
$$\begin{eqnarray}\sqrt{g}\,\text{d}\boldsymbol{x}=\frac{4\unicode[STIX]{x03C0}}{m_{a}^{2}|v_{\Vert }|}\,\text{d}H\,\text{d}\unicode[STIX]{x1D707}\,\text{d}\unicode[STIX]{x1D713}\,\text{d}\unicode[STIX]{x1D6FC}\,\text{d}l,\end{eqnarray}$$
where
$$\begin{eqnarray}H=\unicode[STIX]{x1D705}+e_{a}\unicode[STIX]{x1D719}(\boldsymbol{R})\end{eqnarray}$$
is the total (kinetic
$+$
potential) guiding-centre energy. The parallel adiabatic invariant is
$$\begin{eqnarray}J(\unicode[STIX]{x1D707},H,\unicode[STIX]{x1D713},\unicode[STIX]{x1D6FC})=\int \sqrt{2m_{a}(H-\unicode[STIX]{x1D707}B-e_{a}\unicode[STIX]{x1D719})}\,\text{d}l,\end{eqnarray}$$
and its derivative with respect to
$H$
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}J}{\unicode[STIX]{x2202}H}=\int \frac{\text{d}l}{|v_{\Vert }|}=\unicode[STIX]{x1D70F}_{b},\end{eqnarray}$$
the bounce time. Hence we can write
$$\begin{eqnarray}E=\mathop{\sum }_{a}\frac{4\unicode[STIX]{x03C0}}{m_{a}^{2}}\int \frac{\unicode[STIX]{x1D705}f_{a}\,\text{d}l}{|v_{\Vert }|\unicode[STIX]{x1D70F}_{b}}\,\text{d}\unicode[STIX]{x1D707}\,\text{d}J\,\text{d}\unicode[STIX]{x1D713}\,\text{d}\unicode[STIX]{x1D6FC}.\end{eqnarray}$$
For slow fluctuations (compared with the bounce frequency), the distribution function is nearly independent of
$l$
, and thus
$$\begin{eqnarray}W=E+\unicode[STIX]{x1D6F7},\end{eqnarray}$$
where
$$\begin{eqnarray}E=\mathop{\sum }_{a}\frac{4\unicode[STIX]{x03C0}}{m_{a}^{2}}\int \unicode[STIX]{x1D716}f_{a}\,\text{d}\unicode[STIX]{x1D707}\,\text{d}J\,\text{d}\unicode[STIX]{x1D713}\,\text{d}\unicode[STIX]{x1D6FC},\end{eqnarray}$$
and
$\unicode[STIX]{x1D716}(\unicode[STIX]{x1D707},J,\unicode[STIX]{x1D713},\unicode[STIX]{x1D6FC})$
is defined as the bounce average of the kinetic gyro-centre energy,
$$\begin{eqnarray}\unicode[STIX]{x1D716}=\frac{1}{\unicode[STIX]{x1D70F}_{b}}\int \unicode[STIX]{x1D705}\frac{\text{d}l}{|v_{\Vert }|}.\end{eqnarray}$$
With this definition of the particle energy
$\unicode[STIX]{x1D716}(\boldsymbol{x})$
in (A 2), which is equal to the local kinetic energy if the electrostatic potential is a flux function (i.e. only depends on
$\unicode[STIX]{x1D713}$
), the conditions for nonlinear stability are thus satisfied.
