2. Loss cones and analytic models

Before discussing the loss-cone instabilities themselves, it is important to discuss the structures of the loss cones themselves. Consider a system with magnetic field strength $B_\text {mid}$ at the midplane and field strength $B_\text {end}$ at the mirror ends. Let $R$ denote the mirror ratio, which is defined by

(2.1)\begin{equation} R \doteq \frac{B_\text{end}}{B_\text{mid}} . \end{equation}

Let $\boldsymbol {v}$ denote the velocity of a particle and $m$ denote its mass. Define $B \doteq |\boldsymbol {B}|$, $\hat b \doteq \boldsymbol {B} / B$, $v_{\|} \doteq \hat b \boldsymbol {\cdot } \boldsymbol {v}$, $\boldsymbol {v}_\perp \doteq \boldsymbol {v} - v_{\|} \hat b$ and $v_\perp \doteq |\boldsymbol {v}_\perp |$. Furthermore, let $\mu$ denote the first adiabatic invariant,

(2.2)\begin{equation} \mu \doteq \frac{m v_\perp^2}{2 B} . \end{equation}

In the limit in which collisions can be neglected, and in which the rotation frequency is small compared with the gyrofrequency (Volosov et al. Reference Volosov, Pal'chikov and Tsel'nik1969; Thyagaraja & McClements Reference Thyagaraja and McClements2009), this $\mu$ is conserved. For a particle subject to some potential $\varPhi$ (including the electrostatic potential and, in a rotating system, the centrifugal potential) the conservation of energy can be written as

(2.3)\begin{equation} \tfrac{1}{2} m v_{\|}^2 + \mu B + \varPhi = \text{constant}. \end{equation}

Then if $\Delta \varPhi$ is the difference in $\varPhi$ between the point of maximum $B$ and the midplane (and assuming that $\varPhi$ does not have interior extrema away from the midplane) the condition for a particle to be trapped is

(2.4)\begin{equation} v_{\|}^2 \leq (R-1) v_\perp^2 + \frac{2 \Delta \varPhi}{m} . \end{equation}

This condition is written in terms of the particle's midplane velocity. If (2.4) is rewritten in terms of spherical $(v,\theta,\varphi )$ velocity coordinates, this becomes

(2.5)\begin{equation} \sin^2 \theta \geq \frac{1}{R} \left( 1 - \frac{2 \Delta \varPhi}{m v^2} \right), \end{equation}

which is the well-known generalization of the mirror loss cone to the case with a potential $\varPhi$ (Pastukhov Reference Pastukhov1987). The trapping regions defined by (2.5) are shown in figure 1. When $\Delta \varPhi = 0$, the trapped region is a cone in phase space. When $\Delta \varPhi \neq 0$ and $R > 1$, the trapped region becomes a hyperboloid of revolution. Finally, when $\Delta \varPhi > 0$ and $R < 1$, the trapped region is ellipsoidal. This last case is not a focus of the present paper but it has advantages in certain scenarios (Volosov Reference Volosov2006; Mlodik et al. Reference Mlodik, Munirov, Rubin and Fisch2023; Munirov & Fisch Reference Munirov and Fisch2023). Other cases with $R < 1$ are not shown because these cases do not have a trapped region.

Over the course of this paper, it will often be important to write down a model for a distribution function consistent with this trapping condition for a given $R$ and $\Delta \varPhi$ (neglecting, for simplicity, any spatial variation within the trapped region). Perhaps the simplest choice is a truncated Maxwellian,

(2.6)\begin{equation} f_T(\boldsymbol{v}) = A {\rm e}^{- m v^2 / 2T} \varTheta \left[ \frac{2 \Delta \varPhi}{m} + (R-1) v_\perp^2 - v_{\|}^2 \right]. \end{equation}

Here $\varTheta$ is the Heaviside step function and the normalization constant $A$ is chosen such that

(2.7)\begin{equation} \int \mathrm{d}^3 \boldsymbol{v} \, f(\boldsymbol{v}) = 1. \end{equation}

This model is simple and easy to understand, but does have the disadvantage of having (arguably unphysically) sharp boundaries around the trapping region.

Another choice, essentially equivalent to one used by Volosov et al. (Reference Volosov, Pal'chikov and Tsel'nik1969) and Turikov (Reference Turikov1973), is

(2.8)\begin{equation} f_S(\boldsymbol{v}) = A \sqrt{ \frac{2 \Delta \varPhi}{m} + (R-1) v_\perp^2 - v_{\|}^2} \, {\rm e}^{{-}m v_\perp^2 / 2T} \varTheta \left[ \frac{2 \Delta \varPhi}{m} + (R-1) v_\perp^2 - v_{\|}^2 \right] . \end{equation}

This model has a smoother transition to the loss cone. Its original use by Volosov et al. (Reference Volosov, Pal'chikov and Tsel'nik1969) was motivated by the fact that, when $\Delta \varPhi = 0$ and $R = 2$, it reduces to an analytic model of a non-rotating mirror used by Rosenbluth & Post (Reference Rosenbluth and Post1965). However, it comes with a major disadvantage: it does not reduce to a Maxwellian when $\Delta \varPhi$ is large. This suggests that $f_S$ might be a reasonable model in the limit of slow rotation, but that it is unphysical for larger $\Delta \varPhi$.

For practical purposes, as well as to facilitate easier translation between our results and those elsewhere in the literature, it may be helpful to explicitly show the mapping between $\Delta \varPhi$ and the rotational Mach number. Suppose $\Delta \varPhi$ is entirely due to a centrifugal potential (that is, suppose electrostatic effects can be neglected); then if the angular frequency of rotation $\varOmega$ is constant along a flux surface, and if $r_\text {min}$ and $r_\text {max}$ are the minimum and maximum radii along the flux surface, a particle of mass $m_i$ sees

(2.9)\begin{equation} \Delta \varPhi = \frac{1}{2} m_i (r_\text{max}^2-r_\text{min}^2) \varOmega^2 . \end{equation}

If $R = r_\text {max}^2 / r_\text {min}^2$, then this can be rewritten in terms of the mirror ratio $R$ as

(2.10)\begin{equation} \Delta \varPhi = \frac{R-1}{2 R} m_i r_\text{max}^2 \varOmega^2 . \end{equation}

Then if the rotational Mach number is defined as

(2.11)\begin{equation} {Ma} \doteq \sqrt{ \frac{m_i r_\text{max}^2 \varOmega^2}{2T} } \end{equation}

for a temperature $T$, we have

(2.12)\begin{equation} \frac{\Delta \varPhi}{T} = \frac{R-1}{R} \, {Ma}^2. \end{equation}

Of course, (2.12) can (and, for most real devices, will) be modified by electrostatic potentials, which need to be included in $\Delta \varPhi$ alongside the centrifugal potential. Quasineutrality generally requires the combined electrostatic and centrifugal potentials for electrons and ions to be roughly equal. In the case of singly charged ions and equal ion and electron temperatures, this leads to an electrostatic potential with a magnitude half that of the centrifugal potential, equalizing the species’ confining potentials by trapping the electrons and reducing the trapping of the ions. It will also be modified in cases where $R \neq r_\text {max}^2 / r_\text {min}^2$. This can happen if the plasma occupies a volume with an annular rather than circular cross-section. It can also result from plasma diamagnetism. As such, (2.12) should be understood as a reasonable scaling but not as a precise mapping for all devices. In any case, an important takeaway from (2.12) is that when considering stability thresholds described in the subsequent sections of this paper, one should be aware that $\Delta \varPhi$ can be understood largely as a proxy for the rotation speed but that the mapping between $\Delta \varPhi$ and ${Ma}$ does also depend on the mirror ratio $R$.

3. The HFCLC instability

The first loss-cone mode to be considered is the HFCLC instability. In the regime in which it is unstable, the HFCLC is characterized by travelling waves that grow as they propagate. As we will see, the typical wavelengths are on the scale of the ion Debye length.

The HFCLC is derived using a number of assumptions. First, the ion plasma frequency $\omega _{pi}$ and ion cyclotron frequency $\omega _{ci}$ satisfy $\omega _{pi}\gg \omega _{ci}$. This is a statement about the density regime of interest. Second, the parallel and perpendicular components of the wavenumber $\boldsymbol {k}$ satisfy $k_{\|} \ll k_\perp$. This constraint is often imposed in mirror-type systems due to the large electron mobility. Third, the mode frequency $\omega$ satisfies $\omega _{ci} \ll \omega \ll \omega _{ce}$. This means that the plasma response can be calculated under the assumption that the Larmor gyration of the ions can be ignored but that the electron motion is almost entirely along field lines (Post Reference Post1987). Then, not taking into account the stabilizing effects of finite electron temperature, the electrostatic dispersion relation can be written as (Rosenbluth & Post Reference Rosenbluth and Post1965; Post & Rosenbluth Reference Post and Rosenbluth1966; Post Reference Post1987)

(3.1)\begin{equation} 1 + \frac{\omega_{pe}^2}{\omega_{ce}^2} = \frac{\omega_{pe}^2}{\omega^2} \frac{k_{\|}^2}{k^2} + \frac{\omega_{pi}^2}{k^2 \bar{v}_i^2} F \left( \frac{\omega}{k_\perp \bar{v}_i} \right), \end{equation}

where $\bar {v}_i$ is the ion thermal velocity and

(3.2)\begin{gather} F(\kern0.7pt y) \doteq{-} 2 \int_0^\infty \left( 1 - \frac{x}{y^2} \right)^{{-}1/2} \frac{\partial \psi}{\partial x} \, \mathrm{d} x , \end{gather}

(3.3)\begin{gather}\psi\left( \frac{v_\perp^2}{\bar v_i^2} \right) \doteq \bar\psi \int_{-\infty}^{+\infty} f(v_\perp^2,v_{\|}^2) \, \mathrm{d} v_{\|}, \end{gather}

in which $f$ is the ion distribution function, $x = v_\perp ^2 / \bar v_i^2$, and the normalization constant $\bar \psi$ is chosen such that

(3.4)\begin{equation} \int_0^\infty \psi(x) \, \mathrm{d} x = 1. \end{equation}

Post & Rosenbluth (Reference Post and Rosenbluth1966) showed that the most unstable $\boldsymbol {k}$ satisfies

(3.5)\begin{equation} | \text{Im}(k_{\|}) |_\text{max} \approx \frac{1}{2} \sqrt{\frac{m_e}{m_i}} \frac{\omega_{pi}}{\bar{v}_i} \left( 1 + \frac{\omega_{pe}^2}{\omega_{ce}^2} \right)^{{-}1/2} \left[ - y \text{Im} F(\kern0.7pt y) \right]_\text{max}, \end{equation}

where $m_e / m_i$ is the electron–ion mass ratio; $\omega _{ps}$ is the plasma frequency of species $s$; $\omega _{cs}$ is the cyclotron frequency of species $s$; $-y \text {Im} F(\kern0.7pt y)$ is maximized over $y$.

Note that

(3.6)\begin{equation} -y \text{Im} F(\kern0.7pt y) = 2 y^2 \int_{y^2}^\infty \left( x - y^2 \right)^{{-}1/2} \frac{\partial \psi}{\partial x} \, \mathrm{d} x . \end{equation}

It is clear from (3.6) that the mode cannot be unstable if $\psi$ is a monotonically decreasing function of $x$.Footnote ¹ Similarly, the mode cannot be unstable if there is no choice of $y$ for which $y \text {Im}F(\kern0.7pt y)$ is negative. The second condition turns out to be more easily satisfied than the first. It is also possible to define a third, even more easily satisfied condition by requiring that the mode must have $\text {Im}(\omega ) > \omega _{ci}$, since this is an additional requirement of the HFCLC (Mikhailovskii Reference Mikhailovskii1974). This leads to a requirement that $-y \text {Im} F(\kern0.7pt y)$ be larger than a finite threshold rather than simply that it be positive. This was the approach of Turikov. This third condition is a density-dependent correction to the second, and depends on the value of $\omega _{pi} / k \bar {v}_i$. In the limit where $\omega _{pi} \gg k \bar {v}$, the second and third conditions are equivalent. For the sake of simplicity (and a smaller parameter space), we will focus on the first and second sufficient conditions for stability, which can be summarized as follows.

Perpendicular monotonicity condition:

(3.7)\begin{equation} \frac{\partial \psi}{\partial x} \leq0 \quad \forall x \geq 0; \end{equation}

HFCLC integral condition:

(3.8)\begin{equation} \int_{y^2}^\infty \left( x - y^2 \right)^{{-}1/2} \frac{\partial \psi}{\partial x} \, \mathrm{d} x \leq 0 \quad \forall y > 0. \end{equation}

Either condition serves independently as a sufficient condition for HFCLC stability; the distributions that satisfy the former condition are a subset of the distributions that satisfy the latter.

These conditions can be evaluated for each of the distributions discussed in § 2. Going forward, it will be convenient to work with the dimensionless confinement potential

(3.9)\begin{equation} \phi \doteq \frac{\Delta \varPhi}{T} . \end{equation}

When $\phi > 0$, the potential is confining. First, integrating $f_T$ from (2.6), we find that the corresponding $\psi$ is

(3.10)\begin{equation} \psi_T(x) = \psi_0 {\rm e}^{{-}x} \text{erf} \left[ \sqrt{ (R-1) x + \phi } \right] \varTheta \left[ (R-1) x + \phi \right], \end{equation}

where $\varTheta$ is again the Heaviside step function. The normalization constant $\psi _0$ is

(3.11)\begin{equation} \psi_0^{{-}1} = \text{erf} \left[ \sqrt{\phi} \right] + {\rm e}^{\phi/(R-1)} \text{erfc} \left[ \sqrt{ \frac{R \phi}{R-1} } \right], \quad \text{for } R>1, \phi \geq 0 \end{equation}

and

(3.12)\begin{equation} \psi_0^{{-}1} = {\rm e}^{\phi / (R-1)} \sqrt{\frac{R-1}{R}}, \quad \text{for } R>1, \phi < 0. \end{equation}

Integrating $f_S$ from (2.8) with respect to $v_{\|}$,

(3.13)\begin{equation} \psi_S(x) = \psi_0 \left[ (R-1) x + \phi \right] {\rm e}^{{-}x} \varTheta \left[ (R-1) x + \phi \right] \end{equation}

where

(3.14)\begin{equation} \psi_0^{{-}1} = \phi + R - 1, \quad \text{for } R>1, \phi \geq 0 \end{equation}

and

(3.15)\begin{equation} \psi_0^{{-}1} = (R-1) {\rm e}^{\phi/(R-1)}, \quad \text{for } R>1, \phi < 0. \end{equation}

In figure 2, $\psi _T$ and $\psi _S$ are plotted for several choices of $R$ and $\phi$.

Figure 2. The projected perpendicular-energy distribution functions $\psi _T$ and $\psi _S$ for several choices of $R$ and $\phi$.

The first sufficient condition for stabilization – that is, the monotonicity of $\psi$ – is relatively easy to check. Assuming $R > 1$, neither $\psi _T$ nor $\psi _S$ will meet this condition when $\phi < 0$. Now $\psi _T$ satisfies the perpendicular monotonicity condition when

(3.16)\begin{equation} \sqrt{{\rm \pi} \phi} e^\phi \text{erf}\sqrt{\phi} \geq R-1 , \end{equation}

and $\psi _S$ satisfies the same monotonicity condition when

(3.17)\begin{equation} \phi \geq R - 1. \end{equation}

Note that in both cases, higher $R$ pushes the stability threshold to higher $\phi$.

The second sufficient condition is somewhat more involved to evaluate. It is convenient to change integration variables to $u \doteq \sqrt {x-y^2}$, so that

(3.18)\begin{equation} \int_{y^2}^\infty \left( x - y^2 \right)^{{-}1/2} \frac{\partial \psi}{\partial x} \, \mathrm{d} x = 2 \int_0^\infty \frac{\partial \psi}{\partial x} \, \mathrm{d} u . \end{equation}

When $\phi \geq 0$, the second sufficient condition for $\psi _T$ can be written as

(3.19)\begin{align} & \int_0^\infty {\rm e}^{{-}u^2-y^2} \text{erf} \left[ \sqrt{(R-1) (u^2+y^2) + \phi} \right] \mathrm{d} u \nonumber\\ & \qquad \geq \int_0^\infty \frac{(R-1) {\rm e}^{{-}R(u^2+y^2)-\phi}}{\sqrt{\rm \pi} \sqrt{(R-1)(u^2+y^2)+\phi}} \, \mathrm{d} u \quad \forall y > 0. \end{align}

The right-hand side integral can be evaluated, so that this is equivalent to

(3.20)\begin{align} & \int_0^\infty {\rm e}^{{-}u^2-y^2} \text{erf} \left[ \sqrt{(R-1)(u^2+y^2)+\phi} \, \right] \mathrm{d} u \nonumber\\ & \qquad \geq \frac{\sqrt{R-1}}{2\sqrt{\rm \pi}} \exp \left[ \frac{1}{2} \left({-}R y^2 - 2 \phi + \frac{R \phi}{R-1} \right) \right] {\rm K}_0 \left[ \frac{R}{2} \left( y^2 + \frac{\phi}{R-1} \right) \right] \quad \forall y > 0. \end{align}

Here ${\rm K}_0$ is the modified Bessel function of the second kind. This stability condition is calculated numerically in figure 3.

Figure 3. The first and second HFCLC stability conditions for two choices of model distribution. Either condition is sufficient for stability; the integral condition is more easily satisfied than the monotonicity condition.

The second condition for $\psi _S$ is

(3.21)\begin{align} & \int_0^\infty {\rm e}^{{-}x} \left\{ (R-1) \left[ (R-1) x + \phi \right] \delta \left[ (R-1) x + \phi \right] \right.\nonumber\\ & \qquad \left. - \left[ (R-1)(x-1) + \phi \right] \varTheta \left[ (R-1)x + \phi \right] \right\} \mathrm{d} u \leq 0 \quad \forall y > 0. \end{align}

When $\phi \geq 0$, the condition becomes

(3.22)\begin{equation} \int_0^\infty {\rm e}^{{-}u^2-y^2} \left[ (R-1) (u^2+y^2-1) + \phi \right] \mathrm{d} u \geq 0 \quad \forall y > 0. \end{equation}

This integral can be evaluated explicitly,

(3.23)\begin{equation} \int_0^\infty {\rm e}^{{-}u^2-y^2} \left[ (R-1) (u^2+y^2-1) + \phi \right] \mathrm{d} u = \frac{\sqrt{\rm \pi}}{4} {\rm e}^{{-}y^2} \left[ (R-1) (2 y^2 - 1) + 2 \phi \right] , \end{equation}

so the condition is satisfied when

(3.24)\begin{equation} \phi \geq \frac{R-1}{2} . \end{equation}

To check the condition when $\phi < 0$, it is sufficient to evaluate the integral when $y$ is less than $\sqrt {-\phi /(R-1)}$,

(3.25)\begin{align} & \int_{\sqrt{-\phi/(R-1)}}^\infty {\rm e}^{{-}u^2-y^2} \left[ (R-1) (u^2+y^2-1)+\phi \right] \mathrm{d} u \nonumber\\ & \quad = {\rm e}^{{-}y^2} \left\{ \frac{1}{2 \sqrt{2}} {\rm e}^{\phi/(R-1)} \sqrt{-2 \phi(R-1)} \right.\nonumber\\ & \left.\qquad - \frac{\sqrt{\rm \pi}}{4} \left[ (R-1)(1-2y^2)-2\phi\right] \text{erfc} \left[ \sqrt{\frac{-\phi}{R-1}} \, \right] \right\} . \end{align}

Note that for $R > 1$ and $\phi < 0$,

(3.26)\begin{gather} \frac{1}{2 \sqrt{2}} < \frac{\sqrt{\rm \pi}}{4} , \end{gather}

(3.27)\begin{gather}\sqrt{-(R-1) 2 \phi} < R - 1 - 2 \phi, \end{gather}

(3.28)\begin{gather}{\rm e}^{\phi/(R-1)} < \text{erfc} \sqrt{\frac{-\phi}{R-1}} , \end{gather}

so the second sufficient condition is never met for $\psi _S$ when $\phi < 0$.

These conditions are shown for both model distributions in figure 3. A higher mirror ratio $R$ consistently makes the HFCLC mode more difficult to stabilize. This makes sense; the population inversion appears through the projected perpendicular distribution $\psi (x)$, and higher $R$ makes the loss-cone structure steeper in that distribution.

Especially at higher $R$, the truncated Maxwellian model $f_T$ is stabilized much more easily than the smoothed polynomial prefactor model $f_S$. This makes sense; as $\phi$ becomes larger, $f_S$ remains anisotropic even far from the loss-cone boundary itself, whereas $f_T$ reverts to a Maxwellian everywhere away from the boundary.

There are a variety of additional corrections to the stability condition not considered here. Turikov (Reference Turikov1973) discusses a density-dependent correction to the centrifugal stabilization effect that becomes more important at lower densities. Non-zero electron temperature makes the HFCLC mode easier to stabilize (Post Reference Post1987). So does any other modification of the ion distribution that helps to fill in the population inversion, such as the ‘warm plasma stabilization’ concept (Post Reference Post1967, Reference Post1987). Finite-geometry effects may also help to stabilize the mode (Rosenbluth & Post Reference Rosenbluth and Post1965; Aamodt & Book Reference Aamodt and Book1966; Gerver Reference Gerver1979).

3.1. Stability of distributions from Fokker–Planck simulations

To compare the validity of the two models, we can also generate equilibria using Fokker–Planck simulations. We consider the thermally normalized, sourced, steady-state diffusion equation in pitch angle $\theta$ and normalized velocity $z \doteq |v|/\sqrt {2 T / m}$ (Pastukhov Reference Pastukhov1974; Najmabadi, Conn & Cohen Reference Najmabadi, Conn and Cohen1984; Ochs, Munirov & Fisch Reference Ochs, Munirov and Fisch2023),

(3.29)\begin{gather} 0 = \frac{\partial }{\partial z} \left[\sqrt{g}\left( \frac{z}{z^3+z_0^3} f + \frac{1}{2 (z^3+z_0^3)} \frac{\partial f}{\partial z}\right)\right] + \frac{\partial }{\partial \theta} \left[\sqrt{g} \left(\frac{Z_\perp }{z^3+z_0^3} \frac{\partial f}{\partial \theta} \right)\right] + \sqrt{g} s(z,\theta) , \end{gather}

(3.30)\begin{gather}\sqrt{g} \doteq 4{\rm \pi} z^2 \sin \theta, \qquad s \doteq {\rm \pi}^{{-}3/2} {\rm e}^{{-}z^2/T_s}, \end{gather}

with reflecting boundary conditions everywhere except the loss cone, determined by

(3.31)\begin{equation} R \sin^2 \theta = 1- \frac{\phi}{z^2}. \end{equation}

For $Z = 1$ (singly charged) ions, $Z_\perp = 1/2$, and we additionally take $z_0 = 0.01$.Footnote ² We also allow for different source temperatures $T_s$ relative to the normalized temperature towards which collisions drive the distribution function.

This collision model is a simplified form of the asymptotic high-velocity limit of the Rosenbluth potentials. Although not the most accurate model for thermal ion–ion collisions, it has the nice feature that it very rapidly drives the low-energy part of the distribution to a Maxwellian, due to the high diffusion coefficient, while preserving the diffusion-influenced structure of the solution near the loss cone. For implementation, we use the Fokker–Planck code based on DolfinX (Scroggs et al. Reference Scroggs, Dokken, Richardson and Wells2022) and gmsh (Geuzaine & Remacle Reference Geuzaine and Remacle2009) that was developed in Ochs et al. (Reference Ochs, Munirov and Fisch2023). We use a mesh size of $\Delta z = 0.025$.

The stability boundary results of these simulations for the HFCLC stability criterion (3.8) are shown in figure 4, alongside the truncated Maxwellian and Volosov models, for several values of the normalized source temperature $T_s$. The simulations show the same general trend of decreasing stability with greater $R$ as the truncated Maxwellian model, but the $R$ dependence in the Fokker–Planck simulations is noticeably weaker. To some extent, this is to be expected. In the truncated Maxwellian model, the effect of the loss cone scales with the volume of the loss cone in phase space (since that region is simply empty whereas the rest of the distribution is unaffected). The numerical results include the effects that the loss cone has on the population outside of the loss cone itself. If there is any loss region at a given $v$, pitch-angle scattering into that region will have some tendency to suppress the whole distribution at $v$. So, although the stability boundary predicted from the truncated Maxwellian is qualitatively similar to the simulated results, they are not a perfect match. Note that this dependence is on $R$ for fixed $\phi$, not fixed Mach number. Apart from that, the simulations make it clear that a higher-temperature source somewhat destabilizes the distribution (as expected). Finally, it is clear that the Volosov model does not remotely match the results of the simulations at high $R$.

Figure 4. The HFCLC stability boundary for the Fokker–Planck simulations, with sources at different temperatures, plotted alongside the truncated Maxwellian model threshold $\phi _T$ (corresponding to distribution $f_T$) and Volosov's model $\phi _S$ (corresponding to distribution $f_S$). Estimated error bars based on pseudoerror scaling are plotted, but fall within the square markers. The pseudoerror is calculated by varying the mesh size $\Delta z$, calculating the rate of convergence of the results, and using this to infer the numerical error at the chosen resolution. The qualitative trend of decreasing stability at greater $R$ remains true in the model, though with a less steep slope than for $\phi _T$. Meanwhile, $\phi _S$ very quickly diverges from the simulation lines. See (2.12), and the discussion of the ambipolar potential that follows it, for more on the mapping between $\phi$ and the rotational Mach number.

As before, we can also examine the projected perpendicular energy distribution $\psi (x)$, where $x = v_\perp ^2/\bar v_i^2$, for the simulation results for set values of the mirror ratio and confining potential. This analysis is shown in figure 5. Comparing with figure 2, we see that once again, the curves qualitatively match the behaviour of the truncated Maxwellian model $f_T$, but not of the Volosov model $f_S$.

Figure 5. Projected perpendicular distribution $\psi (x)$ (unnormalized, arbitrary units) from Fokker–Planck simulations for several values of the mirror ratio and confining potential. Comparing with figure 2, we see that the behaviour much more closely matches the truncated Maxwellian model $f_T$ than the Volosov model $f_S$.

4. The DCLC mode

The DCLC mode is another loss-cone instability. In its simplest form, it is electrostatic and has vanishing $k_{\|}$ (Post & Rosenbluth Reference Post and Rosenbluth1966). As was the case for the HFCLC, small $k_{\|}$ is the result of the electron mobility. The DCLC, like the HFCLC, draws free energy from the loss-cone population inversion. Unlike the HFCLC, it also draws free energy from the presence of a (spatial) density gradient. Both factors must be present for the DCLC to occur as an unstable mode. Instability requires that the gradient scale length $|\boldsymbol {\nabla } n / n|^{-1}$ not be very much larger than the ion Larmor radius (Kotelnikov et al. Reference Kotelnikov, Chernoshtanov and Prikhodko2017). The DCLC is related to the drift cyclotron instability, which draws only on the density gradient (Mikhailovsky Reference Mikhailovsky1965). Intuitively, the DCLC can be understood as a kind of drift wave that is destabilized by the presence of the loss cone.

The short-wavelength electrostatic $k_{\|}=0$ DCLC for a slab is probably best known in terms of the following dispersion relation (Post & Rosenbluth Reference Post and Rosenbluth1966; Post Reference Post1987; Kotelnikov et al. Reference Kotelnikov, Chernoshtanov and Prikhodko2017):

(4.1)\begin{equation} w^2 \cot w + \beta w = \beta^{2/3} (\varepsilon a_i) {\rm \pi}^{4/3} \left( \frac{Z m_e}{m_i} + \frac{\omega_{ci}^2}{\omega_{pi}^2} \right)^{{-}2/3} \end{equation}

where $w \doteq {\rm \pi}\omega / \omega _{ci}$, $Z$ is the ion charge state and

(4.2)\begin{equation} \beta \doteq {\rm \pi}(k a_i)^3 \left( \frac{Z m_e}{m_i} + \frac{\omega_{ci}^2}{\omega_{pi}^2} \right) . \end{equation}

In (4.1), $\varepsilon$ is the gradient scale length of the density $n$; in Cartesian $(x,y,z)$ coordinates, if $\boldsymbol {\nabla } n$ is directed in the $y$ direction, it can be written as

(4.3)\begin{equation} \varepsilon \doteq \frac{1}{n} \frac{\partial n}{\partial y} . \end{equation}

Finally, $a_i$ is typically given by

(4.4)\begin{equation} a_i = \frac{1}{\varOmega_i} \left[ \int \frac{f_i}{v_\perp^3} \, \mathrm{d}^3 \boldsymbol{v} \right]^{{-}1/3} . \end{equation}

As written, nothing about this dispersion relation suggests that the monotonicity of $f_i$ should play any special role in the stability of the mode. The shape of the distribution comes into the problem through $a_i$, and it turns out that the value of $a_i$ determines the minimal unstable gradient scale length $\varepsilon$ (Post & Rosenbluth Reference Post and Rosenbluth1966),

(4.5)\begin{equation} (\varepsilon \, a_i)_\text{crit} = 0.38 \left( \frac{Z m_e}{m_i} + \frac{\omega_{ci}^2}{\omega_{pi}^2} \right)^{2/3} . \end{equation}

Suppression of the loss cone does not appear to be stabilizing except in this limited sense.

However, it turns out that (4.1) is derived assuming that $f_i$ must vanish as $v_\perp \rightarrow 0$. The derivation of (4.1) includes two steps in which boundary terms involving $f_i |_{v_\perp = 0}$ are discarded. Equation (4.1) is obtained by integrating over unperturbed particle trajectories; see, for example, the discussions in Post & Rosenbluth (Reference Post and Rosenbluth1966), Swanson (Reference Swanson2003) or Kotelnikov et al. (Reference Kotelnikov, Chernoshtanov and Prikhodko2017). The dispersion relation can be written in terms of the electron and ion susceptibilities $\chi _e$ and $\chi _i$ as

(4.6)\begin{equation} 0 = 1 + \chi_e + \chi_i. \end{equation}

The electron susceptibility is unchanged from what can be found in previous sources,

(4.7)\begin{equation} \chi_e ={-} \frac{\omega_{pe}^2}{\omega_{ce}^2} - \frac{\omega_{pe}^2 \varepsilon}{k_\perp \omega \omega_{ce}}, \end{equation}

where $k = k_\perp$ is the (purely perpendicular) wavenumber. The ion contribution comes from the following:

(4.8)\begin{equation} \chi_i = \frac{2 {\rm \pi}\omega_{pi}^2}{k^2} \int_{-\infty}^{+\infty} \mathrm{d} v_{\|} \int_0^\infty v_\perp \, \mathrm{d} v_\perp \sum_{n={-}\infty}^{\infty} \frac{{\rm J}_n^2(k v_\perp{/} \omega_{ci})}{\omega + n \omega_{ci}} \frac{n \omega_{ci}}{v_\perp} \frac{\partial f_i}{\partial v_\perp} \end{equation}

where ${\rm J}_n$ is a Bessel function of the first kind. Using the identity that $\sum _{n=-\infty }^{\infty } {\rm J}_n^2(x) = 1$, this can be rewritten as

(4.9)\begin{equation} \chi_i = \frac{2 {\rm \pi}\omega_{pi}^2}{k^2} \int_{-\infty}^{+\infty} \mathrm{d} v_{\|} \int_0^\infty \mathrm{d} v_\perp \, \frac{\partial f_i}{\partial v_\perp} \left[ 1 - \sum_{n={-}\infty}^{\infty} \frac{\omega {\rm J}_n^2(k v_\perp{/} \omega_{ci})}{\omega + n \omega_{ci}} \right] . \end{equation}

Then, following Post & Rosenbluth (Reference Post and Rosenbluth1966), the sum can be evaluated in the short-wavelength limit by considering the asymptotic behaviour of the Bessel functions,

(4.10)\begin{align} \sum_{n={-}\infty}^{\infty} \frac{\omega {\rm J}_n^2(k v_\perp{/} \omega_{ci})}{\omega + n \omega_{ci}} & \approx \frac{\omega}{{\rm \pi} k v_\perp} \sum_{n={-}\infty}^{\infty} \frac{1}{n + \omega / \omega_{ci}} \end{align}

(4.11)\begin{align} & = \frac{\omega}{k v_\perp} \cot \left( \frac{{\rm \pi} \omega}{\omega_{ci}} \right). \end{align}

(Here we have taken $k$ to be positive.) This leads to

(4.12)\begin{align} \chi_i & = \frac{2 {\rm \pi}\omega_{pi}^2}{k^2} \int_{-\infty}^{+\infty} \mathrm{d} v_{\|} \int_0^\infty \mathrm{d} v_\perp \frac{\partial f_i}{\partial v_\perp} \left[ 1 - \frac{\omega}{k v_\perp} \cot \left( \frac{{\rm \pi} \omega}{\omega_{ci}} \right) \right] \end{align}

(4.13)\begin{align} & ={-} \frac{2 {\rm \pi}\omega_{pi}^2}{k^2} \int_{-\infty}^{+\infty} \mathrm{d} v_{\|} \left[ f_i \big|_{v_\perp{=} 0} + \frac{\omega}{k} \cot \left( \frac{{\rm \pi} \omega}{\omega_{ci}} \right) \int_0^\infty \mathrm{d} v_\perp \frac{1}{v_\perp} \frac{\partial f_i}{\partial v_\perp} \right]. \end{align}

For cases in which $f_i$ vanishes when $v_\perp$ is small, the first term in square brackets vanishes and the second can be integrated by parts to get the form of $a_i$ given in (4.4). For more general distributions, (4.1) becomes

(4.14)\begin{equation} w^2 \cot w + \beta w = \beta^{2/3} (\varepsilon a_i) {\rm \pi}^{4/3} \left( \frac{Z m_e}{m_i} + \frac{\omega_{ci}^2}{\omega_{pi}^2} + \frac{2 {\rm \pi}\omega_{ci}^2}{k^2} \int_{-\infty}^{+\infty} \mathrm{d} v_{\|} \, f_i \big|_{v_\perp{=} 0} \right)^{{-}2/3} \end{equation}

with

(4.15)\begin{equation} \beta = {\rm \pi}(k a_i)^3 \left( \frac{Z m_e}{m_i} + \frac{\omega_{ci}^2}{\omega_{pi}^2} + \frac{2 {\rm \pi}\omega_{ci}^2}{k^2} \int_{-\infty}^{+\infty} \mathrm{d} v_{\|} \, f_i \big|_{v_\perp{=} 0} \right) \end{equation}

and

(4.16)\begin{equation} a_i = \frac{1}{\omega_{ci}} \left[ \int \frac{1}{v_\perp^2} \frac{\partial f_i}{\partial v_\perp} \, \mathrm{d}^3 \boldsymbol{v} \right]^{{-}1/3}. \end{equation}

This can be understood as a generalization of the original DCLC dispersion relation to include more general distributions. Equation (4.14) can also be written as

(4.17)\begin{equation} w^2 \cot w + \beta w = {\rm \pi}^2 k^2 \varepsilon a_i^3. \end{equation}

Equations (4.15) and (4.16) are the same coefficients $\beta$ and $a_i$ described by (4.2) and (4.4), evaluated in the more general case where $f_i$ need not vanish as $v_\perp \rightarrow 0$. The form of $a_i$ described by (4.4) is always positive, but the more general $a_i$ in (4.16) can be negative, and is guaranteed to be negative any time $\int f_i \, \mathrm {d} v_{\|}$ is a monotonically decreasing function of $v_\perp$. When $a_i > 0$, there is always a threshold value of $\varepsilon$ beyond which this dispersion relation yields an unstable DCLC mode. This is no longer guaranteed when $a_i$ is negative.

To see this, it is helpful to go back to the original argument for the DCLC critical gradient given by Post & Rosenbluth (Reference Post and Rosenbluth1966). This goes as follows. Suppose $\beta > 0$ (as is always the case when $a_i > 0$). Then $w^2 \cot w + \beta w$ has a local maximum between $w = 0$ and $w = {\rm \pi}$. This implies that there is a value of $\varepsilon$ beyond which two real $w$ solutions vanish. The vanishing of these two real solutions corresponds to the appearance of two complex solutions (one of which will have an imaginary part of each sign), and the onset of instability.

Such a critical gradient can always be found so long as $\beta > 0$, and $\beta > 0$ so long as $a_i > 0$. However, for more general $\beta$, the local maximum persists only if $\beta > -1$. For $\beta < -1$, there is no critical gradient at which the number of unstable modes changes, and we can consider the DCLC stabilized. This constitutes a sufficient condition for stability, which can be written as follows.

DCLC integral condition:

(4.18)\begin{equation} \frac{{\rm \pi} k^3}{\omega_{ci}^3} \left[ \int \frac{1}{v_\perp^2} \frac{\partial f_i}{\partial v_\perp} \, \mathrm{d}^3 \boldsymbol{v} \right]^{{-}1} \left( \frac{Z m_e}{m_i} + \frac{\omega_{ci}^2}{\omega_{pi}^2} + \frac{2 {\rm \pi}\omega_{ci}^2}{k^2} \int_{-\infty}^{+\infty} \mathrm{d} v_{\|} \, f_i \big|_{v_\perp{=} 0} \right) \leq{-}1 . \end{equation}

Recall that $k > 0$ by assumption. The condition can be rewritten as

(4.19)\begin{equation} 0 \leq{-} \int \frac{1}{v_\perp^2} \frac{\partial f_i}{\partial v_\perp} \, \mathrm{d}^3 \boldsymbol{v} \leq \frac{{\rm \pi} k^3}{\omega_{ci}^3} \left( \frac{Z m_e}{m_i} + \frac{\omega_{ci}^2}{\omega_{pi}^2} + \frac{2 {\rm \pi}\omega_{ci}^2}{k^2} \int_{-\infty}^{+\infty} \mathrm{d} v_{\|} \, f_i \big|_{v_\perp{=} 0} \right) . \end{equation}

As we saw for the HFCLC, it is the population inversion in the projection $\int f \, \mathrm {d} v_{\|}$ that matters. In terms of the perpendicular energy distribution $\psi$ defined by (3.3), this can be rewritten as

(4.20)\begin{equation} 0 \geq \int_0^\infty x^{{-}1/2} \frac{\partial \psi}{\partial x} \, \mathrm{d} x \geq{-} \frac{{\rm \pi} k^3 \bar v_i^3}{2 \omega_{ci}^3} \left( \frac{Z m_e}{m_i} + \frac{\omega_{ci}^2}{\omega_{pi}^2} + \frac{2 \omega_{ci}^2}{k^2 \bar v_i^2} \, \psi(0) \right) . \end{equation}

The first inequality in (4.20) is the same as the HFCLC integral condition when $y \rightarrow 0$. For the HFCLC integral condition, the choice of $y$ adjusts which parts of $\partial \psi / \partial x$ are weighted most heavily. For the model distributions $f_S$ and $f_T$, positive $\partial \psi / \partial x$ (when it happens) is concentrated at the smallest values of $x$. Therefore, the HFCLC condition and the first inequality in the DCLC condition yield the same results for $f_S$ and $f_T$. However, one could construct a distribution for which this would not be the case (for example, a distribution for which $\partial \psi / \partial x$ went from negative to positive to negative again as $x$ increased).

The second inequality in the DCLC integral condition is less intuitive. However, it may be a less serious limitation than it initially appears. Note that the DCLC integral condition is satisfied any time

(4.21)\begin{equation} 0 \leq{-} \int_{-\infty}^{+\infty} \mathrm{d} v_{\|} \int_0^\infty \frac{\omega_{ci}}{k v_\perp} \frac{\partial f_i}{\partial v_\perp} \, \mathrm{d} v_\perp{\leq} - {\rm \pi}\int_{-\infty}^{+\infty} \mathrm{d} v_{\|} \int_0^\infty \frac{\partial f_i}{\partial v_\perp} \, \mathrm{d} v_\perp . \end{equation}

Recall that the dispersion relation used in this section was derived using a short-wavelength assumption. Specifically, in (4.9), the integrand was simplified by taking the asymptotic limit $k v_\perp / \omega _{ci} \gg 1$. The two integrands in (4.21) differ by this same factor.

This tells us two things. First, within the limits of applicability of the dispersion relation, the same perpendicular monotonicity condition discussed for the HFCLC guarantees that the DCLC integral condition will also be satisfied. Second, even in cases where $\psi$ may not be monotonic, the second inequality in the DCLC integral condition is not a very strong constraint. A distribution can satisfy the first inequality in (4.20) but fail to satisfy the second only if $k v_\perp / \omega _{ci}$ is not too large for $v_\perp$ on the scale on which $f_i$ varies; it is not clear that the underlying dispersion relation will still be valid in such a case.

Numerically, it turns out that the stability boundaries from the HFCLC integral condition and the DCLC integral condition coincide for the cases shown in figures 5 and 5, for the same reason that they coincide for the analytic models $f_T$ and $f_S$: the regions in which $\psi$ has the most positive slope are concentrated at low $v_\perp$. Therefore, figure 4 can also be read as the stability boundary for the DCLC integral condition.

As was the case for the HFCLC analysis, there are a variety of additional effects not considered here, some of which can be stabilizing for practical devices. These include finite-$\beta$ effects, finite system size and the same ‘warm plasma stabilization’ mentioned in § 3 (Tang, Pearlstein & Berk Reference Tang, Pearlstein and Berk1972; Berk & Gerver Reference Berk and Gerver1976; Gerver Reference Gerver1976; Kotelnikov et al. Reference Kotelnikov, Chernoshtanov and Prikhodko2017). We have also not considered the effects of multiple-ion-species distributions, which can be non-trivial (Kotelnikov & Chernoshtanov Reference Kotelnikov and Chernoshtanov2018).

5. The DGH mode

The last mode we will discuss in detail is the DGH mode (Dory et al. Reference Dory, Guest and Harris1965). In particular, we will focus on the electrostatic $k_{\|}=0$ ‘zero-frequency’ DGH mode, in which $\text {Re}(\omega )$ vanishes. This is not the most general form of the mode; see, for example, Callen & Guest (Reference Callen and Guest1971). However, it has the advantage of simplicity, and it was the focus of much of the original work on the subject. Physically, this mode can be understood in terms of single-particle drifts in an inhomogeneous electric field (Post Reference Post1987). An appropriate wave field can cause corrections to the homogeneous $\boldsymbol {E} \times \boldsymbol {B}$ drift that cause particles with different charges to move in opposite directions. Depending on the kinetic distribution of particles, this can amplify these waves and lead to instability. As was the case for the HFCLC and DCLC, loss-cone distributions can have the necessary characteristics to drive the instability.

If the mode is driven by an ion loss cone, the ion contribution follows from (4.9); then the relevant dispersion relation is

(5.1)\begin{equation} 1 + \frac{\omega_{pe}^2}{\varOmega_e^2} = \frac{2 {\rm \pi}\omega_{pi}^2}{k^2} \int_{-\infty}^{+\infty} \mathrm{d} v_{\|} \int_0^\infty \mathrm{d} v_\perp \, \frac{\partial f_i}{\partial v_\perp} \left[ 1 - \sum_{n={-}\infty}^{\infty} \frac{\omega {\rm J}_n^2(k v_\perp{/} \omega_{ci})}{\omega + n \omega_{ci}} \right] . \end{equation}

A marginally stable distribution with $\text {Re}(\omega ) = 0$ satisfies

(5.2)\begin{equation} 1 + \frac{\omega_{pe}^2}{\varOmega_e^2} = \frac{2 {\rm \pi}\omega_{pi}^2}{k^2} \int_{-\infty}^{+\infty} \mathrm{d} v_{\|} \int_0^\infty \mathrm{d} v_\perp \, \frac{\partial f_i}{\partial v_\perp} \left[ 1 - {\rm J}_0^2(k v_\perp{/} \omega_{ci}) \right] . \end{equation}

An essentially equivalent expression was found by Post & Rosenbluth (Reference Post and Rosenbluth1966). If the right-hand side of (5.2) is positive, then the zero-frequency DGH mode is stable (Post & Rosenbluth Reference Post and Rosenbluth1966). Note that ${\rm J}_0^2(x) \leq 1$ for all $x$. It follows immediately that the same monotonicity condition we have discussed for the HFCLC and DCLC modes – that is, (3.7) – is a sufficient condition for the stability of this mode as well. A more easily met sufficient condition for stability is to have

(5.3)\begin{equation} \int_{-\infty}^{+\infty} \mathrm{d} v_{\|} \int_0^\infty \mathrm{d} v_\perp \frac{\partial f_i}{\partial v_\perp} \left[ 1 - {\rm J}_0^2(k v_\perp{/} \omega_{ci}) \right] \leq 0. \end{equation}

In terms of the perpendicular energy distribution $\psi$, this stability condition can be written as follows.

DGH integral condition:

(5.4)\begin{equation} \int_0^\infty \frac{\partial \psi}{\partial x} \left[ 1 - {\rm J}_0^2 (\xi x^{1/2}) \right] \,\mathrm{d} x \leq0 \end{equation}

where

(5.5)\begin{equation} \xi \doteq \frac{k \bar v_i}{\omega_{ci}} . \end{equation}

The same stability condition can be written equivalently as

(5.6)\begin{equation} \int_0^\infty \xi x^{{-}1/2} {\rm J}_0 (\xi x^{1/2}) {\rm J}_1(\xi x^{1/2}) \psi(x) \, \mathrm{d} x \geq 0. \end{equation}

This condition was also discussed by Post & Rosenbluth (Reference Post and Rosenbluth1966).

The DGH mode is much easier to stabilize than the HFCLC or DCLC modes. Consider the implications of the Bessel function weighting factor in (5.4). One minus the square of the Bessel function ${\rm J}_0^2$ vanishes when its argument is small, then exhibits damping oscillations with a mean that tends towards unity as the argument increases. Mirror-like equilibria, including both $f_T$ and $f_S$, typically have their most destabilizing parts (that is, the regions of velocity space in which the population inversion is steepest) at small $v_\perp$; this can be seen in figure 2. The exception, for $f_S$ and $f_T$, is when $\phi$ is negative (when the potential is deconfining rather than confining). When $\phi < 0$, the regions with the most positive $\partial \psi / \partial x$ are shifted to larger values of $x$, below which the distribution entirely vanishes, and it is possible to find a choice of $R$ and $\xi$ for which the DGH integral condition is not satisfied. However, $\phi < 0$ is not the case that is relevant for centrifugal confinement.

The numerical solutions described in § 3 are also stable to the DGH in most cases. However, it is possible to find DGH-unstable equilibria in cases with relatively high-temperature sources (for example, $T_s = 2$) and little rotation ($\phi < 1/2$). Because the numerical instability threshold is strongly determined by the source term, and because the DGH is consistently stabilized at a $\phi$ threshold below that shown in figure 4, we will not focus on this effect here. The two main things that are worth noting, for present purposes, are (i) that the DGH is stable for $\phi \geq 0$ for the analytic models but not always for the numerical simulations and (ii) that even in the numerical simulations the DGH still appears to be more easily stabilized than either the HFCLC or the DCLC.

6. Free energy and stability of flute-like modes

A common thread across all of these modes is that stability is guaranteed whenever $\psi$ is a monotonically decreasing function: that is, whenever

(6.1)\begin{equation} \frac{\partial}{\partial v_\perp^2} \int_{-\infty}^{+\infty} f(v_\perp^2, v_{\|}^2) \, \mathrm{d} v_{\|} \leq 0. \end{equation}

As we have seen, this condition is sufficient but not necessary for each of the modes in question. Nonetheless, it is interesting because of its simplicity and because of the fact that it appears in the stability analyses for a range of different modes.

Another commonality (which applies to all loss-cone modes, not just the ones considered here) is that the physical intuition for how they work is typically posed in terms of free energy and phase-space rearrangements. Systems with loss cones typically have population inversions. This means that energy can be released by rearranging the contents of phase space such that high-energy particles move to lower energy. This energy then becomes available to drive instabilities.

This intuition is closely related to the idea underlying the theory of plasma free energy, wherein the capacity of a system for instability is set by the amount of energy that can be released by certain classes of phase-space rearrangements. Given a set of rules for which phase-space rearrangements are allowed, it is possible to calculate the free energy (sometimes called the ‘available energy’) by determining the maximal energy that can be extracted from the system using these rearrangements. This energy could then drive instabilities, or one could imagine extracting it intentionally, as in the case of alpha-channelling (Fisch & Rax Reference Fisch and Rax1992).

Several different classes of rearrangement have been proposed, each corresponding to a different way of defining the free energy. Gardner (Reference Gardner1963) suggested a formulation that is now sometimes called ‘Gardner restacking’, in which any operation that preserves phase space densities is permitted (Dodin & Fisch Reference Dodin and Fisch2005; Helander Reference Helander2017, Reference Helander2020). An alternative proposed by Fisch & Rax (Reference Fisch and Rax1993), sometimes called the ‘diffusively accessible free energy’, instead allows operations that mix or average the populations between phase space elements (Hay, Schiff & Fisch Reference Hay, Schiff and Fisch2015, Reference Hay, Schiff and Fisch2017; Kolmes, Helander & Fisch Reference Kolmes, Helander and Fisch2020a; Kolmes & Fisch Reference Kolmes and Fisch2020, Reference Kolmes and Fisch2022). There are variants of both formulations that allow for additional restrictions to be imposed on the allowed rearrangements, enforcing a conservation law – for example, on an adiabatic invariant $\mu$ – by requiring that exchange or mixing operations must act only between elements of phase space with the same value of $\mu$ (Helander Reference Helander2017, Reference Helander2020; Kolmes et al. Reference Kolmes, Helander and Fisch2020a). Recent evidence suggests that adding suitable adiabatic invariants to the Gardner free energy leads to a metric that can predict the saturation levels of some types of turbulence (Mackenbach, Proll & Helander Reference Mackenbach, Proll and Helander2022; Mackenbach et al. Reference Mackenbach, Proll, Snoep and Helander2023a, Reference Mackenbach, Proll, Wakelkamp and Helanderb).

When we try to understand the behaviour of flute-like loss-cone modes in terms of free energy theory, we arrive at an apparent paradox. On the one hand, free energy/rearrangement theories appear to match our intuitions for the physics behind loss-cone instabilities. Loss-cone instabilities happen because it is possible to release energy by rearranging the distribution in velocity space; this is exactly the kind of situation that free energy theories should be able to describe. On the other hand, the mirror ratio dependence found in §§ 3 and 4, and shown in figure 3, is not at all consistent with what would have been predicted from either the Gardner or the diffusive exchange theories. A larger mirror ratio means a smaller loss cone, with a correspondingly smaller energy release if that loss cone is filled in, but the stabilization thresholds for the HFCLC and DCLC suggest that configurations with larger mirror ratios are actually more difficult to stabilize.

To resolve this contradiction, it is helpful to begin by considering the condition for a distribution to be a ‘ground state’ in different versions of free energy theory. In the absence of any additional conservation laws, a system is in a ground state with vanishing free energy only when the distribution $f$ is a monotonically decreasing function of energy. That is, if $f = f(v^2)$, and if the energy is $m v^2 / 2$, being in a ground state requires that

(6.2)\begin{equation} \frac{\mathrm{d} f(v^2)}{\mathrm{d} v^2} \leq 0 \quad \forall v. \end{equation}

This holds for both the Gardner and the diffusive theories. If $\mu$ invariance is enforced, a distribution of the form $f = f(v^2, \mu )$ is a ground state only if (Helander Reference Helander2017)Footnote ³

(6.3)\begin{equation} \frac{\partial f(v^2, \mu)}{\partial v^2} \leq 0 \quad \forall v. \end{equation}

Considering the resemblance of these conditions to (6.1), and the similarity of the intuitions behind the free-energy theories and the loss-cone modes, it is natural to wonder whether the behaviour of these loss-cone instabilities can be understood in terms of the free energy.

To motivate a free energy with ground states corresponding to (6.1), consider the standard quasilinear diffusion equation describing the velocity-space diffusion induced by wave–particle interactions,

(6.4)\begin{equation} \frac{\partial f}{\partial t} = \frac{\partial}{\partial \boldsymbol{v}} \boldsymbol{\cdot} {\mathsf{D}} \boldsymbol{\cdot} \frac{\partial}{\partial \boldsymbol{v}} f, \end{equation}

where $ {\mathsf {D}}$ is a rank-2 tensor given by

(6.5)\begin{equation} {\mathsf{D}} \doteq D_0 \int \frac{\omega_i \mathcal{E}_{\boldsymbol{k}}}{(\boldsymbol{k} \boldsymbol{\cdot} \boldsymbol{v} - \omega_r)^2 + \omega_i^2} \frac{\boldsymbol{k} \boldsymbol{k}}{k^2} \, \mathrm{d} \boldsymbol{k}, \end{equation}

for some species-dependent constant $D_0$ and spectral wave energy density $\mathcal {E}_{\boldsymbol {k}}$, and with $\omega _r$ and $\omega _i$ denoting the real and imaginary parts of the wave frequency $\omega$ (Krall & Trivelpiece Reference Krall and Trivelpiece1973). In the limit where $k_{\|}$ is vanishingly small, this diffusion operator acts only in perpendicular directions, and does not distinguish between different values of $v_{\|}$.

If we expect phase space to be rearranged by interactions with small-$k_{\|}$ modes, then it is sensible to consider the free energy associated with mixing operations that can act only on the projected distribution function $\int _{-\infty }^{+\infty } f \, \mathrm {d} v_{\|}$. Then the ground-state condition immediately recovers (6.1).

This rule bears some similarity to the restacking conservation laws introduced by Helander (Helander Reference Helander2017, Reference Helander2020). Indeed, it leads to the conservation of $v_{\|}$, but it is more restrictive than Helander's rule. The difference comes from the assumption that these flute-like modes not only cannot move particles in the parallel velocity-space direction, but they also cannot perform rearrangements that distinguish between populations with the same $v_\perp$ but different $v_{\|}$. This distinction is illustrated in figure 6.

Figure 6. Two different kinds of restriction on the allowed rearrangement operations, each enforcing a conservation law. (a) $\mu$ conservation is enforced by allowing any rearrangement that acts on two elements of phase space with the same value of $\mu$. (b) $v_{\|}$ is conserved by instead enforcing that rearrangement operations must act to exchange populations between the entire region of phase space with one value of $v_\perp$ and the entire region with another $v_\perp$.

7. Rearrangements in loss-cone systems

Thus far, the discussion of the allowed rearrangements has not taken into account the existence of the loss cone. The preceding sections have considered loss-cone instabilities as those instabilities which result from the particle distributions that are characteristic of systems with loss cones, without considering the effects of the phase-space loss region itself. In one sense, there is nothing wrong with this; there is no reason why loss-cone-like distributions, and their associated instabilities, cannot exist without the presence of loss regions in phase space. (In other words, it would be self-consistent to imagine a mirror-like distribution of particles that happens to exist in an infinite homogeneous space, and to calculate its stability properties). However, if we want to apply free energy theory to cases in which there are regions of phase space from which particles are promptly lost, then we need to modify our accessible states accordingly.

Once a particle enters the loss cone, it exits the system. This means that it should not be possible to ‘fill up’ a loss cone. It is possible to model the energetic fate of these vanishing particles in more than one way, but the simplest approach is to assert that they leave the system without any further change in energy. Then a particle starting at energy $\varepsilon _i$ and leaving the system with energy $\varepsilon _f$ leaves behind energy $\varepsilon _i - \varepsilon _f$ within the system. This energy difference is thus ‘available’ in the same sense that is usually measured by the Gardner free energy. The situation would be very different if there were some other pathway by which the liberated energy could exit the system.

When describing rearrangement processes in the presence of a loss cone, it is notationally convenient (and physically equivalent) to instead say that any region of phase space within the loss cone interacts with other parts of phase space as if it had zero population. As an illustrative example, consider the four-state discrete system in which the states have dimensionless energies $\varepsilon _0 = 0$, $\varepsilon _1 = 1$, $\varepsilon _2 = 2$ and $\varepsilon _3 = 3$. Suppose the second-lowest-energy state is part of a region of phase space from which particles are promptly lost; we will denote this by shading it light blue. Then the state

has energy $\sum _i \varepsilon _i \, f_i$. If not for the presence of the loss cone, Gardner restacking would lead to a ground state simply by sorting the $f_i$ from largest to smallest. For example, it would take

\begin{align*} \begin{array}{|c|c|c|c|} \hline 4 & 0 & 2 & 1 \\ \hline \end{array} & \rightarrow \begin{array}{|c|c|c|c|} \hline 4 & 2 & 0 & 1 \\ \hline \end{array} \rightarrow \begin{array}{|c|c|c|c|} \hline 4 & 2 & 1 & 0 \\ \hline \end{array} \end{align*}

for a total energy release of $(2 \varepsilon _2 + \varepsilon _3) - (2 \varepsilon _1 + \varepsilon _2) = 3$. On the other hand, if the second cell is part of a loss cone, the mapping would instead be

for a release of $(2 \varepsilon _2 + \varepsilon _3) - 3 \varepsilon _1 = 4$. The loss-cone region never runs out of space; we still imagine that phase-space densities are being conserved, because this notation is shorthand for the idea that there are many copies of that box sitting at the same energy (corresponding to a spatial coordinate as particles exit the confinement device). In other words, the notation describes a projection of phase space in which rearrangements into loss cones can appear to produce condensate-like accumulations of particles, even though phase space volumes are still conserved.

Now consider the corresponding continuous problem for a simple loss cone with no added trapping or detrapping potential. This is the case shown in figure 1(a). In this scenario, one can see that the Gardner free energy is equal to the entire kinetic energy of the system, since any populated region of phase space can exit the system through the part of the loss cone with zero energy. In the more general case where the lowest-energy region of the loss cone has energy $\varepsilon _\text {min} > 0$ (for example, figure 1b), any particle that starts with energy $\varepsilon > \varepsilon _\text {min}$ can be removed from the system at $\varepsilon _\text {min}$, though the details of the ground state will also depend on the initial conditions for the parts of phase space with $\varepsilon < \varepsilon _\text {min}$. The situation is similar if we consider diffusive rearrangements rather than Gardner restacking; an element of phase space with initial energy above the lowest-energy region of the loss cone may be repeatedly averaged against that empty region until it has entirely been removed from the system at the minimal loss-cone energy.

This is straightforward enough so far. Things become more complicated if we consider the flute-like rearrangements discussed in § 6. Consider, for example, a discrete system in which we distinguish between the parallel and perpendicular energy in different boxes as follows:

Here the state with $(\varepsilon_\perp, \varepsilon _{\|}) = (0,1)$ is taken to be a loss region. This is roughly the kind of loss-cone structure we would expect for the loss cone of a system with $R > 1$ and some confining potential.

Suppose we insist, following the discussion in § 6, that rearrangements must act between whole columns of states with given values of $\varepsilon _\perp$. One example of the Gardner restacking procedure might be

This produces a ground state subject to the given constraints. But then consider the following initial condition:

This system has only one allowed exchange operation, and it leads to a higher-energy state,

This transforms the system from a state with energy 2 to a state with energy 3. However, consider the following sequence:

This sends the system from energy 2 to energy 3, and then to energy 1. In other words, this initial condition cannot be transformed into any state with lower energy with any single allowed operation, but because of the interaction between the flute-like constraint and the loss-cone, it is possible to map it to a lower-energy state by using a sequence of two operations. The intermediate step, in which the energy of the system is raised, is sometimes called an ‘annealing operation’ (Hay et al. Reference Hay, Schiff and Fisch2015; Kolmes & Fisch Reference Kolmes and Fisch2022). A very similar scenario can be found if diffusive exchange operations are used in place of Gardner restacking.

Annealing operations are never necessary to reach the ground state in the original versions of the restacking or diffusive exchange theories (that is, without this interaction between the flute-like-mode constraint and the loss cone) (Hay et al. Reference Hay, Schiff and Fisch2015). Their appearance in this version of the theory suggests that we should make a distinction between two different kinds of ground states. The first, which we might call a weak ground state, is any state which cannot be mapped to a lower-energy state using any single allowed operation. The second, which we might call a strong ground state, is any state which cannot be mapped to a lower-energy state using any sequence of allowed operations. Of course, every strong ground state is also a weak ground state. The perpendicular monotonicity condition found in the context of the HFCLC, DCLC and DGH modes corresponds to the weak ground states of the constrained system.