Non-thermal particle acceleration and power-law tails via relaxation to universal Lynden-Bell equilibria

Abstract

Collisionless and weakly collisional plasmas often exhibit non-thermal quasi-equilibria. Among these quasi-equilibria, distributions with power-law tails are ubiquitous. It is shown that the statistical-mechanical approach originally suggested by Lynden-Bell (Mon. Not. R. Astron. Soc., vol. 136, 1967, p. 101) can easily recover such power-law tails. Moreover, we show that, despite the apparent diversity of Lynden-Bell equilibria, a generic form of the equilibrium distribution at high energies is a ‘hard’ power-law tail $\propto \varepsilon ^{-2}$, where $\varepsilon$ is the particle energy. The shape of the ‘core’ of the distribution, located at low energies, retains some dependence on the initial condition but it is the tail (or ‘halo’) that contains most of the energy. Thus, a degree of universality exists in collisionless plasmas.

1. Introduction

It is well known that the ultimate fate of a homogeneous collisional plasma is to become a Maxwellian. This result was first inferred for neutral particles by Maxwell (1860) on statistical grounds and given solid dynamical foundation by Boltzmann (1896) with his collision integral. Plasma physics was to wait for Landau (1936) and later Balescu (1960) and Lenard (1960) to be equipped with its own collision integral, and the resulting universality. Nevertheless, distributions with power-law tails, a far cry from Maxwellian equilibria, are observed in a myriad of plasma systems including cosmic rays (Becker Tjus & Merten 2020; Amato & Casanova 2021), the solar corona and solar flares (Dudík et al. 2017; Oka et al. 2018), the solar wind (Gloeckler et al. 2008; Fisk & Gloeckler 2014; Livadiotis, Desai & Wilson 2018; Moncuquet et al. 2020), the Earth's magnetosheath (Birn et al. 2012; Ergun et al. 2020), and laser plasmas (Cruz et al. 2018; Hartouni et al. 2022).

That such non-Maxwellian distributions emerge should perhaps come as no surprise. In a plasma, the time scale of relaxation to Maxwellian equilibrium is associated with two-body Coulomb collisions, but, due to the long-range nature of the forces involved, the plasma may evolve, exchanging energy between fields and particles, on much shorter time scales. Indeed, in the absolute absence of collisions, the Vlasov equation has an infinite set of nonlinearly stable equilibria: all distributions that are monotonically decreasing functions of particle energy are certainly stable (Gardner 1963). However, the set of all monotonically decreasing functions of energy is very large, and is certainly not exhausted by the Maxwellian equilibrium, which depends on only two parameters. It is, therefore, an outstanding challenge to determine whether any of these stable collisionless equilibria are naturally favoured by the dynamics of the system in a way that is not sensitive to initial conditions, i.e. whether a degree of universality exists in collisionless plasmas. It is certainly the case that nature appears to prefer distributions with power-law tails, and direct numerical simulations have indicated that power-law tails are the natural result of a number of dynamical processes including relativistic and non-relativistic shocks (Sironi & Spitkovsky 2010; Caprioli & Spitkovsky 2014; Crumley et al. 2019), magnetic reconnection (Sironi & Spitkovsky 2014; Werner & Uzdensky 2021; Uzdensky 2022), and various types of plasma turbulence (Kunz, Stone & Quataert 2016; Zhdankin et al. 2017, 2019; Comisso & Sironi 2018, 2022; Zhdankin 2021, 2022b).

In addition to suggesting dynamical paths towards distributions with power-law tails, there have been multiple attempts to justify the ubiquity of such distributions from a thermodynamic point of view. This is, however, entangled with the question of whether the standard Gibbs–Shannon entropy is applicable to systems with long-range interactions and, if not, then what entropy should be used. Naturally, many entropies have emerged to fill this niche. A popular contender is the Tsallis (1988) entropy (or $\alpha$-structural entropy: see Havrda & Charvát 1967). The Tsallis entropy was designed to be a non-additive version of the Gibbs–Shannon entropy as a way to model systems with correlations that therefore should be non-extensive (see, e.g. Livadiotis & McComas (2009), Pierrard & Lazar (2010), and references therein). While this model produces good fits to observed distributions, it has a free parameter that is needed to quantify the degree of the non-extensivity and cannot be determined without fitting data, or additional input of physics currently lacking (note some recent progress suggesting that this additional physics might be deducible from free-energy considerations: Zhdankin 2022a,b).

An early attempt to tackle the question of entropy in collisionless systems was made by Lynden-Bell (1967). Let us consider a system of $N$ particles with canonical positions $\boldsymbol {r}_{i}$ and momenta $\boldsymbol {p}_{i}$ that evolve subject to a Hamiltonian $\mathcal {H}(\boldsymbol {r}_{i},\boldsymbol {p}_{i})$. Such a system can be said to be ‘collisionless’ if the evolution equation for the single-particle distribution function $f(\boldsymbol {r},\boldsymbol {p})$ is well approximated by an effective Hamiltonian $\mathcal {H}^{\mathrm {eff}}(\boldsymbol {r},\boldsymbol {p})$ acting on a single particle (i.e. if the mean-field dynamics are a sufficiently good approximation to the true dynamics), viz.,

(1.1)\begin{equation} \frac{\partial{f}}{\partial{t}} + \frac{\partial{\mathcal{H}^{\mathrm{eff}}}}{\partial{\boldsymbol{p}}}\cdot\frac{\partial{f}}{\partial{\boldsymbol{r}}} -\frac{\partial{\mathcal{H}^{\mathrm{eff}}}}{\partial{\boldsymbol{r}}}\cdot\frac{\partial{f}}{\partial{\boldsymbol{p}}} = 0. \end{equation}

In his original treatment, Lynden-Bell focused on relaxation of stellar systems, but the spirit of his statistical mechanics is the same for all collisionless systems, including plasmas. While keeping the calculations as general as possible, one can think of (1.1) as the collisionless Vlasov equation for a plasma, which could be electrostatic or electromagnetic in a non-relativistic or relativistic regime. The collisionless dynamics described by (1.1) conserve an infinite number of invariants, equivalent to conserving the volume of level sets of the distribution function $f(\boldsymbol {r},\boldsymbol {p})$ in phase space. Thus, the dynamics can be viewed as an extremely complicated rearranging of the elements of phase space, which, however much they are distorted and stirred, will keep the same level sets (often referred to as ‘waterbags’, in analogy with parcels of incompressible fluid). Lynden-Bell posited that, after a short time, the exact phase-space density $f(\boldsymbol {r},\boldsymbol {p})$ would become so chaotic that it could be treated as a random field and that any measurement of it – in practice, of a coarse-grained version of it – was in fact a measurement of the mean phase-space density. This allowed the construction of a statistical mechanics, with an entropy closely related to the Gibbs–Shannon entropy, that encoded an infinite number of invariants and thus predicted the steady states from a given initial condition. These steady states are the Lynden-Bell equilibria.

Since its genesis, Lynden-Bell's theory (often referred to as the theory of ‘violent relaxation’) has received continued attention both thermodynamically (Chavanis, Sommeria & Robert 1996; Arad & Johansson 2005; Chavanis 2006a,b; Levin, Pakter & Teles 2008; Levin et al. 2014) and dynamically, viz., effective ‘collisionless collision integrals’ have been proposed that recovered Lynden-Bell equilibria as their fixed points (Kadomtsev & Pogutse 1970; Severne & Luwel 1980; Chavanis 2004, 2022; Ewart et al. 2022). However, the main strength of the theory is also its weakness. Unlike in the non-extensive entropy formulations, there is no ad hoc parameter in the Lynden-Bell theory: equilibria are uniquely determined by the ‘waterbag content’ of the initial conditions. However, this necessarily means that the equilibria depend (seemingly, in a complicated way) on an infinite family of invariants (sometimes referred to as ‘Casimirs’). This has limited any actual calculations with Lynden-Bell equilibria to simplified situations with only a small number of level sets (in practice, between one and three, e.g. Assllani et al. 2012). At any rate, the intricate dependence on an infinite family of invariants might not appear to be a step towards general power-law tails or any other meaningful form of universality.

To see just how non-universal Lynden-Bell equilibria can be, one only needs to consider the relation between the Lynden-Bell equilibria and the aforementioned Gardner distributions. Should the initial distribution be a monotonically decreasing function of particle energy, then it is a Gardner distribution and there are no possible rearrangements of the phase volume that do not increase energy. Hence, the only state available via collisionless dynamics is this Gardner distribution, which must therefore be its own Lynden-Bell equilibrium.¹ But since any monotonically decreasing function of energy is a Gardner distribution, these minimum-energy states are clearly highly non-universal. However, this is only a good intuition for systems where the number of level sets is small or where phase-volume conservation conspires with energy conservation to render much of the phase space inaccessible to the system (as is the case for Gardner distributions). In this paper, by solving for the full Lynden-Bell equilibria numerically (as well as analytically, in a tractable limit), we will show that most Lynden-Bell equilibria are much more generic. Namely, we will show that, in the limit of a continuum of level sets, and for energies sufficiently greater than the ground-state energy (the energy of the corresponding Gardner distribution), the Lynden-Bell equilibria exhibit power-law tails at high energies, typically with a scaling of $\varepsilon ^{-2}$, where $\varepsilon$ is the particle energy.

The physical argument for these power-law tails is as follows. Phase-volume conservation effectively makes the particles occupying each waterbag (level set of the distribution) behave as if they were members of a separate species, which can communicate with the other waterbags only via the equilibration of some effective ‘temperature’ subject to competition for the same volumes of phase space. In essence, this turns the system into a ensemble of many different fermionic species, all of which exclude each other. When the system is in its minimum-energy (ground) state, the Gardner distribution, the competition for phase volume is the overpowering factor, giving the distribution a highly non-universal shape. However, when the energy of the system is increased, more of the phase space becomes accessible, so, as in Fermi–Dirac statistics, the competition for the same phase volume becomes weaker. For sufficiently large energies, the competition for any volume of phase space is minuscule. In this limit, each waterbag will form its own Maxwellian distribution, in thermal equilibrium with all other waterbags. However, despite these Maxwellian equilibria having the same effective temperature, they will have different thermal spreads because waterbags of larger phase-space density ‘cost’ more energy to be placed at a given momentum $p$ in phase space. The true distribution function is recovered by summing up (in the limit of many waterbags, integrating) the contributions from each of these Maxwellians to the mean phase-space density. This procedure naturally gives rise to a power-law tail that depends on the relative weighting for each Maxwellian (this is qualitatively similar to the formalism of ‘superstatistics’; cf. Beck & Cohen 2003; Chavanis 2006a; Davis et al. 2023). This weighting turns out to depend only weakly on the level sets of the initial condition. For a wide class of initial conditions, the resulting Lynden-Bell equilibria turn out to have the same universal power-law tail, $\varepsilon ^{-2}$.²

The rest of this paper is organised as follows. In § 2, we will review briefly the Lynden-Bell formalism, to state the problem and establish notation. We will then proceed to perform a systematic exploration of the nature of the Lynden-Bell equilibria. In § 3.1, we will argue that to each initial condition, one can uniquely assign a Gardner distribution function with the same Casimir invariants (waterbag content). All Lynden-Bell equilibria can thus be viewed as the result of adding some amount of energy to a Gardner distribution with the same waterbag content and letting it reach a maximum-entropy state (one can think of this approach as describing how adding energy to a population of collisionless particles causes them to form a ‘non-thermal’ distribution with a tail). In § 3.2, we will show that the function describing the waterbag content of a large class of Gardner distributions has a relatively generic form, which will contribute to the universality of the resulting equilibria. In § 3.3, we will solve for the Lynden-Bell equilibria in the limit where the energy of the system far exceeds the energy of the corresponding Gardner distribution. This will ensure that competition for volumes of phase space can be neglected. This makes the problem analytically tractable, and the resulting analytical solution will exhibit the universal power-law tail $\propto \varepsilon ^{-2}$ at high energies. In § 4, by solving for the Lynden-Bell equilibria numerically, we will show that the qualitative features of this analytical solution are retained even for energies that are of the same order as the energy of the Gardner distribution. Therefore, a large class of Lynden-Bell equilibria display a universal power-law tail. This tail contains much of the distribution's energy, whereas the low-energy ‘core’ retains some dependence on the initial conditions. In § 5, we summarise our findings and discuss their implications for real (observed) plasmas.

2. Lynden-Bell's statistical mechanics

In this section, we present a brief re-derivation of Lynden-Bell's equilibria as applied to a homogeneous system. We begin with the collisionless Vlasov equation (1.1) evolving a single species of particles. As well as particle number, momentum, and total energy (i.e. of fields and particles), (1.1) conserves an infinite number of ‘Casimir’ invariants, e.g. the volume of phase space where the exact phase-space density is greater than a given value:

(2.1)\begin{equation} \Gamma(\eta) = \iint\mathrm{d}\boldsymbol{r}\, \mathrm{d}\boldsymbol{p}H(f(\boldsymbol{r},\boldsymbol{p}) - \eta) = \mathrm{const.}, \end{equation}

where $H(x)$ is the Heaviside function (unity for $x>0$ and zero otherwise). As discussed above, despite the existence of these invariants, the system's evolution can still be highly chaotic, which prompted Lynden-Bell to consider the exact phase-space density $f(\boldsymbol {r},\boldsymbol {p})$ as a random field. Therefore, one may introduce the probability density $P(\boldsymbol {r},\boldsymbol {p},\eta )$ for the exact phase-space density $f(\boldsymbol {r},\boldsymbol {p})$ to take the value $\eta$ at position $(\boldsymbol {r},\boldsymbol {p})$ (Robert & Sommeria 1991; Chavanis et al. 1996). The distinct values of $\eta$ will be referred to as ‘waterbags’ since this term conjures up the correct mental image: parcels of phase space of a certain density that can be distorted and moved, but not rarefied, compressed, or superimposed. The mean phase-space density is then

(2.2)\begin{equation} \langle f \rangle(\boldsymbol{p}) = \int \mathrm{d}\eta \eta P(\boldsymbol{p},\eta). \end{equation}

Here, we have applied the intuition that the steady-state distribution function $P(\boldsymbol {p},\eta )$ will be homogeneous in space (this contrasts with Lynden-Bell's original treatment, which focused on gravitationally bound, and therefore inhomogeneous, systems). Lynden-Bell's statistical mechanics amounts to positing that, before the onset of ‘true’ collisions, $P(\boldsymbol {p},\eta )$ will maximise the Gibbs–Shannon entropy

(2.3)\begin{equation} S ={-}\iint\mathrm{d}\boldsymbol{p}\, \mathrm{d}\eta P(\boldsymbol{p},\eta)\ln P(\boldsymbol{p},\eta). \end{equation}

Note that the integral in $\eta$ must run over all the possible values, including $\eta = 0$ (the empty waterbag).

Equipped with an entropy ripe for maximisation, we must decide upon a set of reasonable constraints under which to maximise it. Naturally, since $P(\boldsymbol {p},\eta )$ is a probability-density function in $\eta$ at a given $\boldsymbol {p}$, its integral in $\eta$ at any $\boldsymbol {p}$ must equal unity:

(2.4)\begin{equation} \int\mathrm{d}\eta P(\boldsymbol{p},\eta) = 1. \end{equation}

As well as this, we fix the energy density of the system:

(2.5)\begin{equation} \iint \mathrm{d}\boldsymbol{p}\, \mathrm{d}\eta \varepsilon(\boldsymbol{p})\eta P(\boldsymbol{p},\eta) = E = \mathrm{const.}, \end{equation}

where $\varepsilon (\boldsymbol {p})$ is the energy of a particle as a function of its momentum $\boldsymbol {p}$. Note that, within this formalism, one could include the interaction energy of particles with fields (electromagnetic, gravitational, etc.), so that in its most general form $\varepsilon$ would be a function of both position and momentum, which could need to be solved self-consistently with $P$. Here, we will neglect this rich complexity, assuming instead that, in the relaxed state, the energy of the fields has decayed to a negligible fraction of the total energy and, in the process of decaying, has mediated the relaxation of the distribution function.

Next, we enforce the conservation of the Casimir invariants (2.1) by requiring that the volume-integrated probability of each waterbag stays constant, viz.,

(2.6)\begin{equation} \int \mathrm{d}\boldsymbol{p} P(\boldsymbol{p},\eta) = \rho(\eta) = \mathrm{const.} \end{equation}

The function $\rho (\eta )$ will be referred to as the ‘waterbag content’ and is determined by initial conditions. The waterbag content of the initial condition can be read off by integrating over all portions of phase space where the initial exact phase-space density is equal to a particular value, viz.,

(2.7)\begin{equation} \rho(\eta) = \frac{1}{V}\iint\mathrm{d}\boldsymbol{r}\, \mathrm{d}\boldsymbol{p}\delta(\eta - f(\boldsymbol{r},\boldsymbol{p},t = 0)) ={-}\frac{1}{V}\frac{\mathrm{d}\Gamma}{\mathrm{d}\eta}, \end{equation}

where $V$ is the system's spatial volume. A priori, in Lynden-Bell's statistical mechanics, the degree of universality of the equilibrium distribution is determined by $\rho (\eta )$: all initial conditions with the same waterbag content and energy lead to the same equilibrium.

Maximising the entropy (2.3) subject to the constraints (2.4)–(2.6)³ is equivalent to maximising, unconditionally, the functional

(2.8)\begin{align} & S[P(\boldsymbol{p},\eta)] - \int\mathrm{d}\boldsymbol{p}\lambda(\boldsymbol{p})\left[\int \mathrm{d}\eta P(\boldsymbol{p},\eta) - 1\right] - \beta \left[\iint \mathrm{d}\boldsymbol{p}\, \mathrm{d}\eta \varepsilon(\boldsymbol{p}) \eta P(\boldsymbol{p},\eta) - E \right]\nonumber\\ & \quad + \beta\int \mathrm{d}\eta \eta \mu(\eta)\left[ \int \mathrm{d}\boldsymbol{p} P(\boldsymbol{p},\eta) - \rho(\eta)\right], \end{align}

where $\lambda (\boldsymbol {p})$, $\beta$ and $-\beta \eta \mu (\eta )$ are Lagrange multipliers. By analogy with textbook statistical mechanics, we will sometimes refer to $\mu (\eta )$ as the ‘chemical potential’ (which it is, being the Lagrange multiplier that fixes the number of particles in waterbag $\eta$). Doing so, we find the Lynden-Bell equilibria

(2.9)\begin{equation} P(\boldsymbol{p},\eta) = \frac{{\rm e}^{-\beta \eta[\varepsilon(\boldsymbol{p}) - \mu(\eta) ]}}{\displaystyle \int \mathrm{d}\eta'\, {\rm e}^{-\beta \eta'[\varepsilon(\boldsymbol{p})- \mu(\eta') ]}}, \end{equation}

where $\lambda (\boldsymbol {p})$ has been computed explicitly to arrange for the correct normalisation (2.4), whereas $\beta$ and $\mu (\eta )$ must be chosen in such a way as to satisfy the constraints (2.5) and (2.6). We note that, despite (2.6), the mean phase-space density $\langle f \rangle$, given by (2.2), will, in general, not have the same level sets (2.1) as the exact one $f$ (since $\langle f \rangle$ is an averaged quantity). The equilibria (2.9) are both homogeneous and isotropic: an inevitable consequence of the system having no preferred position or direction.

The similarity between the Lynden-Bell equilibria (2.9) and the Fermi–Dirac distribution is immediately apparent.⁴ This should come as no surprise, because phase-volume conservation functions analogously to Pauli's exclusion principle: pieces of the same waterbag, or different waterbags, cannot cohabit in phase space. The equilibria, therefore, have degeneracy effects incorporated within them.

The prescription for computing Lynden-Bell equilibria is now clear: given an initial condition, with the initial energy density $E$ and waterbag content $\rho (\eta )$ (determined by (2.7)), solve two coupled integral equations (2.5) and (2.6) with $P(\boldsymbol {p},\eta )$ given by (2.9), determine $\beta$ and $\mu (\eta )$, and substitute back into (2.9). Before considering the numerical solutions of this problem in § 4, we will first seek to understand the system analytically.

3. Theory: degenerate and non-degenerate equilibria

3.1. Lynden-Bell equilibria as excited Gardner distributions

Just as it is only meaningful to consider a Maxwellian with a positive energy, it is only meaningful to solve the Lynden-Bell equilibria (2.9) subject to reasonable choices of the constraints (2.5) and (2.6). It is therefore instructive to understand the properties of the waterbag-content function (2.6). To get a feel for waterbag contents of typical initial conditions, one could compute the integral (2.7) for a range of examples (this is relatively simple due to the presence of the delta function). One quickly discovers that many different initial conditions have similar waterbag contents, just as many different initial conditions can have the same energy. To see this, we note that, from (2.7), any volume-preserving transformation of the coordinates $(\boldsymbol {r},\boldsymbol {p})$, including those transformations that ‘splice’ the phase space discontinuously, will leave the waterbag content unchanged. This is unsurprising because the true, incompressible, flow of probability in phase space is precisely one such volume-preserving transformation. It is this freedom that implies that vastly different initial conditions can possess identical, or similar, waterbag contents. A cartoon illustrating this is given in figure 1, showing how seemingly complex distributions have the same waterbag content as very simple distributions. There will be families of initial conditions that have the same waterbag content, but different energies.

Figure 1. A cartoon contour plot in phase space of three possible distribution functions, all of which possess identical waterbag contents. Panel (a) shows the Gardner distribution function corresponding to this waterbag content. Panels (b) and (c) show distributions, at different (higher) energies, which can be reduced to the Gardner distribution by deforming and splicing the phase space incompressibly. A small patch of phase space is highlighted in red between plots to show the effect of the deformation.

For every family of initial conditions possessing the same waterbag content, there will be a unique distribution function that has that waterbag content but is a monotonically decreasing function of energy and, therefore, has the minimum possible energy associated with that waterbag content. Such a distribution function, for which the exact phase-space density satisfies $f(\boldsymbol {r},\boldsymbol {p}) = f_{G}(\varepsilon (\boldsymbol {p}))$, is known as the Gardner distribution (Gardner 1963; Helander 2017). The sequence of deformations of the distribution function to map an initial condition to its Gardner distribution function is often referred to as a ‘restacking’, as it amounts to a reordering of phase-space elements into their minimum-energy configuration (Dodin & Fisch 2005; Kolmes & Fisch 2020; Kolmes, Helander & Fisch 2020). Gardner distributions can be viewed as ‘ground states’ associated with a given waterbag content (e.g. Helander 2017), since no more energy can be extracted from such a distribution without violating phase-volume conservation. This fact intuitively guarantees that any initial condition that is a Gardner distribution is its own Lynden-Bell equilibrium since no other states are available to the system. Indeed, one can show that any Gardner distribution can be reconstructed from (2.9) for a particular choice of $\beta$ and $\mu (\eta )$, although the proof is technical and left to Appendix A.

A generic initial condition can then be viewed as equivalent to taking some Gardner distribution and driving it out of equilibrium by the injection of some energy without changing the waterbag content. The Lynden-Bell equilibria are then simply the collisionless, phase-volume preserving, entropy-maximising equilibria of these higher-energy states, making them the natural excited states of Gardner distributions. Therefore, to capture the set of all possible waterbag contents, we need only study the set of all these ‘ground states’, to which we would then add energy – the first step in the direction of universal outcomes. Physically, this approach is equivalent to asking to what distribution a population of collisionless particles will relax once a certain amount of energy is injected into it – in a manner of speaking, a ‘thermodynamical’ approach to ‘non-thermal’ particle acceleration.

3.2. Waterbag content of Gardner distributions

Having stated the problem in this way, we now consider the waterbag content associated with Gardner distributions. In what follows, we will consider Gardner distribution functions that are truncated at some minimum phase-space density $\eta _{\min }$. This mathematical convenience will turn out to be a physical necessity. Thankfully, while it is mathematically and physically important that the cutoff $\eta _{\min }$ be finite, it will only appear logarithmically in the outcomes of our calculations, making them highly insensitive to its actual value – yet another theory where the need for a cutoff is unavoidable but non-lethal.

As a prototypical example, we compute $\rho (\eta )$ for a particular Gardner distribution: a truncated Maxwellian, viz.,

(3.1)\begin{equation} f_{G}(\boldsymbol{p}) = \begin{cases} \eta_{\max} \, {\rm e}^{-\varepsilon(\boldsymbol{p})/ \varepsilon_{0}} \quad & \text{for} \ \varepsilon(\boldsymbol{p}) < \varepsilon_{0}\ln\displaystyle\frac{\eta_{\max}}{\eta_{\min}},\\ 0 & \text{for} \ \varepsilon(\boldsymbol{p}) >\varepsilon_{0}\ln \displaystyle\frac{\eta_{\max}}{\eta_{\min}}. \end{cases} \end{equation}

Besides $\eta _{\min }$, the parameters of this distribution are the energy scale $\varepsilon _{0}$ and the maximum phase-space density $\eta _{\max }$. The latter is straightforwardly related to the particle's spatial density, e.g. ${n_{0} = \eta _{\max }(2{\rm \pi} m \varepsilon _{0})^{3/2}}$ in the limit ${\eta _{\min }/ \eta _{\max } \to 0}$ for a three-dimensional, non-relativistic plasma, where $\varepsilon (\boldsymbol {p}) = p^{2}/2m$. The waterbag content of the Gardner distribution for such a plasma is then, from (2.7),

(3.2)\begin{equation} \rho(\eta) = {\Gamma}_{\mathrm{free}}\delta(\eta) + \begin{cases} \displaystyle\frac{2n_{0}}{\sqrt{\rm \pi}\eta_{\max}\eta} \left(\ln\displaystyle\frac{\eta_{\max}}{\eta} \right)^{1/2} & \text{for } \eta_{\min} < \eta <\eta_{\max},\\ 0 & \text{otherwise}, \end{cases} \end{equation}

where $\Gamma _{\mathrm {free}}$ is the total volume of the momentum space that is unoccupied (i.e. where the exact phase-space density is zero). Of course, in reality, momentum space is unbounded and so $\Gamma _{\mathrm {free}}$ is infinite. Formally, we are solving for the waterbag content and Lynden-Bell equilibrium in a momentum space of large, but finite, volume and will take this volume to infinity at the end of the calculation – of course, nothing physical will depend on $\Gamma _{\mathrm {free}}$ as it becomes large.

The presence of a $\delta (\eta )$ term in the expression for $\rho (\eta )$ is a generic feature, not restricted to the specific example (3.2). When it comes to solving (2.6) with $P(\boldsymbol {p},\eta )$ given by (2.9), in order to find $\mu (\eta )$, this delta function can be accommodated by writing the chemical potential as

(3.3)\begin{equation} {\rm e}^{\beta \eta \mu(\eta)} = \eta_{\mathrm{ref}}\delta(\eta) + \begin{cases} \eta_{\mathrm{ref}}F(\eta) & \text{for } \eta_{\min} < \eta < \eta_{\max}, \\ 0 & \text{otherwise}, \end{cases} \end{equation}

where $\eta _{\mathrm {ref}}$ is some reference constant that must have dimensions of phase-space density. Its value is unimportant because $e^{\beta \eta \mu (\eta )}$ can always be rescaled by a constant without changing the Lynden-Bell equilibrium (2.9) – in essence, $\eta _{\mathrm {ref}}$ is a gauge choice for the function $\mu (\eta )$. By analogy to textbook statistical mechanics, the function $F(\eta )$ will be referred to as the ‘fugacity’ of the distribution. The form (3.3) results in the following expression for the Lynden-Bell equilibrium (2.9):

(3.4)\begin{equation} P(\boldsymbol{p},\eta) = \frac{\delta(\eta) + {\rm e}^{-\beta \eta \varepsilon(\boldsymbol{p})}F(\eta)}{\displaystyle 1 + \int_{\eta_{\min}}^{\eta_{\max}} \mathrm{d}\eta'\, {\rm e}^{-\beta \eta'\varepsilon(\boldsymbol{p})}F(\eta')}. \end{equation}

In (3.4), the first, $\delta (\eta )$, term in the numerator accounts for the probability density of finding phase space to be empty at a given location (this part of the phase space is referred to, aptly, as the ‘vacuum’ by Chavanis 2006a), whereas the second term accounts for non-empty waterbags. Already $\Gamma _{\mathrm {free}}$ has dropped out of the calculation, as it must, and it is safe to let $\Gamma _{\mathrm {free}} \to \infty$. Likewise, it is immediately obvious that the reference phase-space density $\eta _{\mathrm {ref}}$ has cancelled, as it also must. Thus, the Lynden-Bell equilibrium distribution (3.4) depends only on the Lagrange multiplier $\beta$ and the fugacity $F(\eta )$ (which themselves depend on $\rho (\eta )$ for $\eta > 0$ and the energy density $E$).

The key feature of the example (3.2) is the $\eta ^{-1}$ power-law behaviour. While the specific form of the logarithmic factor in (3.2) was set by our (non-universal) choice of a Maxwellian Gardner distribution, the $\eta ^{-1}$ behaviour at small $\eta$ is relatively universal. For any exponentially decaying Gardner distribution, i.e. for any distribution that at large energies can be approximated by $\propto \exp [-(\varepsilon /\varepsilon _{0} )^{\sigma } ]$, for some $\sigma > 0$ (or indeed can be bounded between two such functions), one will find a waterbag content with an $\eta ^{-1}$ asymptotic at $\eta \ll \eta _{\max }$.

To see this, we note that, since $f_{G}(\varepsilon )$ is monotonically decreasing with energy, it has a well-defined inverse $f_{G}^{-1}(\eta )$. In terms of this inverse, one can explicitly express the waterbag content (2.7) as

(3.5)\begin{equation} \rho(\eta) = \int\mathrm{d}\varepsilon g(\varepsilon)\delta(\eta - f_{G}(\varepsilon) ) ={-}g(f_{G}^{{-}1}(\eta))\frac{\mathrm{d}f_{G}^{{-}1}}{\mathrm{d}\eta} \implies \frac{\mathrm{d}f_{G}}{\mathrm{d}\varepsilon} ={-}\frac{g(\varepsilon)}{\rho(f_{G}(\varepsilon))}, \end{equation}

where $g(\varepsilon )$ is the density of states in energy, defined by the equation

(3.6)\begin{equation} \int \mathrm{d}\boldsymbol{p} (\cdots) = \int \mathrm{d}\varepsilon g(\varepsilon) (\cdots). \end{equation}

The first equality in (3.5) is a straightforward way to calculate the waterbag content of a given Gardner distribution, while the second is an equation from which the Gardner distribution can be constructed given knowledge of the system's waterbag content (cf. Dodin & Fisch 2005; Helander 2017). It is immediately clear why the $\eta ^{-1}$ scaling should arise in $\rho (\eta )$. For any exponentially decaying $f_{G}(\varepsilon )$, the inverse function will be some logarithmic function of $\eta$, which, after differentiation in (3.5), will give an $\eta ^{-1}$ asymptotic multiplied by some logarithmic function of $\eta$. To illustrate this, in figure 2, we give three examples of starkly different Gardner distributions that, despite their differences, all possess waterbag contents which scale as $\eta ^{-1}$ at low $\eta$.

Figure 2. (a) Three example Gardner distribution functions (the phase-space density here is plotted as a function of energy) and (b) their corresponding waterbag contents. All three distributions were chosen to have the same particle density $n_{0}$ and energy density $E_{G}$. The maximum phase-space density $\eta _{\max }$ of the distribution sets the upper cutoff of the waterbag content in $\eta$, shown by the dashed vertical lines in (b). The lower cutoff $\eta _{\min }$ is justified in § 3.5. We see that large differences at low $\varepsilon$ only change the behaviour of $\rho (\eta )$ significantly at $\eta \sim \eta _{\max }$. For $\eta \ll \eta _{\max }$, all three waterbag contents asymptote to a universal $\eta ^{-1}$ scaling.

To convince a doubtful reader, we consider an alternative argument in support of the $\eta ^{-1}$ scaling of $\rho (\eta )$. First, let us imagine a system in which we are allowed to vary $\eta _{\min }$ freely while leaving the exact phase-space density $f(\boldsymbol {v})$ otherwise unchanged, as if there were some true distribution that had $\eta _{\min } = 0$ and we were examining successive approximations to it (e.g. the difference between a truncated Maxwellian (3.1) and a true Maxwellian). Let us consider the following integrals of such distribution functions (similar to the Casimir invariants considered by Zhdankin 2022a):

(3.7)\begin{equation} \frac{1}{V}\iint_{f(\boldsymbol{r},\boldsymbol{p}) > \eta_{\min}}\, \mathrm{d}\boldsymbol{r}\, \mathrm{d}\boldsymbol{p} [f(\boldsymbol{r},\boldsymbol{p})]^{\gamma} = \int_{\eta_{\min}}^{\eta_{\max}} \mathrm{d}\eta \eta^{\gamma}\rho(\eta) , \end{equation}

where we have used (2.7), and the phase-space integral is taken over the volumes where $\eta > \eta _{\min }$. We have also used the property that the process of varying $\eta _{\min }$ while otherwise leaving $f(\boldsymbol {r},\boldsymbol {p})$ unchanged only changes the integration limits of $\rho (\eta )$, without changing $\rho (\eta )$ itself. The idea now is to vary the values of $\eta _{\min }$ and $\gamma$ and use what we know about $f(\boldsymbol {r},\boldsymbol {p})$ to deduce the form of $\rho (\eta )$. Clearly, for $\gamma = 1$, (3.7) is the particle density of the truncated $f(\boldsymbol {r},\boldsymbol {p})$, which must be finite. This tells us that $\rho (\eta )$ must integrate to a finite value when multiplied by $\eta$. Furthermore, should $f$ be any exponentially decaying function, then there would be a characteristic momentum scale above which the distribution is suppressed. This means that, as $\eta _{\min }$ is taken to zero, both sides of (3.7) must converge for $\gamma = 1$. This is effectively a statement that the amount of probability contained beyond a few standard deviations is small. However, for an exponentially decaying phase-space density $f$ and any positive power $\gamma > 0$, $f^{\gamma }$ is also exponentially decaying, so the same argument applies. Therefore, for exponentially decaying phase-space densities, $\rho (\eta )$ must be such that, when multiplied by any positive power of $\eta$, it integrates to some finite value and is largely independent of the choice of lower cutoff $\eta _{\min }$ of the integral. However, for $\gamma = 0$, (3.7) becomes the (spatially averaged) momentum-space volume occupied by the truncated distribution: its support. For exponentially decaying distributions, which do not have compact support without truncation, this quantity will continue to grow without bound as $\eta _{\min }$ is decreased. Therefore, $\rho (\eta )$ integrated with no powers of $\eta$ must diverge as $\eta _{\min } \to 0$. It is obvious that $\eta ^{-1}$ is a function that has all these properties, but, more generally, $\rho (\eta )$ could be any function of the form

(3.8)\begin{equation} \rho(\eta) = \frac{n_{0}}{\eta_{\max}\eta}G\left(\frac{\eta}{\eta_{\max}} \right) \quad \text{for } \eta_{\min} < \eta < \eta_{\max}, \end{equation}

where $G(x)$ is a dimensionless function whose dependence on $x$ is weaker than any power law, viz.,

(3.9)\begin{equation} \lim_{x \to 0} x^{\gamma - 1} G(x) = \begin{cases} 0 & \text{if } \gamma > 1, \\ \infty & \text{if } \gamma < 1. \end{cases} \end{equation}

Note that the exact limit at $\gamma = 1$ cannot be determined by this argument; it is, in fact, dependent on the exact details of the exponential decay, but we will not require it for any calculations. This concludes the argument that Gardner distributions with exponential tails have waterbag contents with a universal low-$\eta$ asymptotic.

It is a straightforward extension of this argument to work out what the waterbag content will be for Gardner distributions with non-exponential tails at high energies (low $\eta$). Suppose that, instead of being exponentially decaying, the Gardner distribution behaves as a power law at large momenta. Then there is a choice of $\gamma > 0$ for which $f^{\gamma }$, multiplied by the density of states in $|\boldsymbol {p}|$, decays slower than $|\boldsymbol {p}|^{-1}$, implying that the integral (3.7) will diverge as $\eta _{\min }$ is decreased. This means that $\rho (\eta )$ multiplied by $\eta ^{\gamma }$ for some $\gamma > 0$ must still have a divergent integral, so it must take the form

(3.10)\begin{equation} \rho(\eta) = \frac{n_{0}}{\eta_{\max}^{2-\delta}\eta^{\delta}}G \left(\frac{\eta}{\eta_{\max}} \right) \quad \text{for } \eta_{\min} < \eta < \eta_{\max}, \end{equation}

with some $\delta > 1$. It turns out that (3.10) also holds, but with $\delta < 1$, for phase-space densities that go to zero algebraically even when $\eta _{\min } = 0$.⁵ This is because the integral (3.7) will be over a finite momentum-space volume even as $\eta _{\min }$ is taken to zero, so, while $\gamma < 0$ will make $f^{\gamma }$ diverge in this finite domain, that divergence will still be integrable if $\gamma$ is chosen sufficiently small and negative.

The explicit value of $\delta$ can be found from (3.5). To enable an explicit calculation, we will assume henceforth that the density of states $g(\varepsilon )$ is related to energy by a simple power law, viz.,

(3.11)\begin{equation} g(\varepsilon) = A\varepsilon^{a}, \end{equation}

where $A$ is an appropriate constant with dimensions $[n_0/\eta _{\max }\varepsilon ^{1+a}]$, and $a$ is a real number. The assumption (3.11) is not too restrictive as it can capture both non-relativistic and ultra-relativistic systems of any dimensionality. With the assumption (3.11), we may now use (3.5) to link the value of $\delta$ to the high-energy asymptotic of the Gardner distribution. If the Gardner distribution has the power-law tail $f_{G}(\boldsymbol {p}) \propto \varepsilon (\boldsymbol {p})^{-(\chi + a)}$ at high energies, then, via (3.5), one finds $\delta = 1 + (a+1)/(a+\chi )$. If, instead, the Gardner distribution goes to zero at some finite energy $\varepsilon _{\max }$ so that $f_{G}(\boldsymbol {p}) \propto [\varepsilon _{\max } - \varepsilon (\boldsymbol {p}) ]^{\chi }$ near $\varepsilon = \varepsilon _{\max }$, then $\delta = 1 - 1/\chi$.

Having catalogued the possible Gardner distributions and their corresponding $\rho (\eta )$, an important open question is now what Gardner distributions are the most common. Within the domain of numerical experiments, this is clearly decided by the whims of the numerical experimenter. However, it would seem reasonable to conjecture that in nature, the most common Gardner distributions should be ones with exponential tails. The reason for this is that the only processes that can change the Gardner distribution are, by definition, collisional, and collisional processes naturally relax the system to a Maxwell–Boltzmann distribution. More concretely, one often thinks of the sources of energy for violent relaxation as being large-scale inhomogeneities (e.g. counter-propagating flows, collapsing distributions of matter, etc.) of a system that is locally collisionally relaxed, and, therefore, has an exponential Gardner distribution.

We are now safe in the knowledge that the exponentially decaying Gardner distributions have $\rho (\eta ) \propto \eta ^{-1}$, or, in more exotic cases, $\rho (\eta ) \propto \eta ^{-\delta }$. We may now return to the task of solving (2.5) and (2.6) in light of these facts.

3.3. Non-degenerate Lynden-Bell equilibria

Taking inspiration from the similarity between Fermi–Dirac and Lynden-Bell statistics, we may expect that there are two analytically tractable limits: degenerate (‘cold’) and non-degenerate (‘hot’). In this section, we will explore the latter limit, which will turn out to be far more useful than the former (which is, nevertheless, also treated, for completeness, in Appendix A).

We define the non-degenerate limit as one in which the probability of finding the exact phase-space density to be non-zero is small, viz.,

(3.12)\begin{equation} D(\varepsilon) = \int_{\eta_{\min}}^{\eta_{\max}} \mathrm{d}\eta\, {\rm e}^{-\beta \eta\varepsilon}F(\eta) \ll 1. \end{equation}

We shall call $D(\varepsilon )$ the ‘degeneracy parameter’ since, from (3.4), the probability that a position in phase space with a given energy is non-empty is given by the quotient $D(\varepsilon )/[1+D(\varepsilon )]$. In the limit (3.12), the distribution function (3.4) can be approximated by

(3.13)\begin{equation} P(\boldsymbol{p},\eta) \simeq \delta(\eta)[1 - D(\varepsilon(\boldsymbol{p})) ] + {\rm e}^{-\beta \eta\varepsilon(\boldsymbol{p})}F(\eta). \end{equation}

The effect of this simplification is that the competition for any particular sub-volume of phase space is so weak that the waterbags are free to arrange themselves as Maxwellians $\eta$ by $\eta$. The waterbags with lower $\eta$ cost less energy to be placed at larger momenta – therefore, they have a larger thermal spread. In the non-relativistic limit, this is equivalent to the intuition that particles belonging to the waterbags with higher phase-space densities behave as though they have larger masses.

The approximate form (3.13) of the Lynden-Bell distribution makes computing the momentum-space integral in (2.6) and determining the fugacity $F(\eta )$ a simple matter. Substituting (3.13) into (2.6) and using the explicit form (3.11) of the density of states, we find the fugacity in the non-degenerate limit to be

(3.14)\begin{equation} F(\eta) = \frac{(\beta \eta)^{1+a}}{A \varGamma(a+1)}\rho(\eta), \quad \eta_{\min} \leq \eta \leq \eta_{\max}, \end{equation}

where $\varGamma (a+1)$ is the gamma function. Substituting (3.14) back into (3.13) and using (2.2) finally gives us an expression for the mean phase-space density (although still in terms of the as yet unspecified parameter $\beta$):

(3.15)\begin{equation} \langle f \rangle(\boldsymbol{p}) = \frac{\beta^{1+a}}{A \varGamma(a+1)}\int_{\eta_{\min}}^{\eta_{\max}}\, \mathrm{d}\eta \eta^{2+a}\rho(\eta)\, {\rm e}^{-\beta \eta \varepsilon(\boldsymbol{p})}. \end{equation}

From (3.15), it might seem as though, by diverse choices of $\rho (\eta )$, a wide variety of distribution functions $\langle f \rangle$ can be obtained. However, as we showed in § 3.2, a diversity of choices of waterbag contents is exactly what we do not have. Instead, fairly generic Gardner distributions with any form of exponential tails possess waterbag contents that are highly universal at low $\eta$. Using (3.8), we can make a convenient change of variables in (3.15), $x = \beta \eta \varepsilon (\boldsymbol {p})$, to find

(3.16)\begin{equation} \langle f \rangle(\boldsymbol{p}) = \frac{n_{0}}{A\beta\varGamma(a+1) \eta_{\max}\varepsilon(\boldsymbol{p})^{2+a}} \int_{\beta \eta_{\min}\varepsilon(\boldsymbol{p})}^{\beta \eta_{\max}\varepsilon(\boldsymbol{p})}\, \mathrm{d}xx^{1+a}G\left(\frac{x}{\beta \eta_{\max}\varepsilon(\boldsymbol{p})} \right) {\rm e}^{{-}x}. \end{equation}

This form exposes the fact that there is a natural power-law behaviour at energies such that

(3.17)\begin{equation} \frac{1}{\beta \eta_{\max}} \ll \varepsilon(\boldsymbol{p}) \ll \frac{1}{\beta \eta_{\min}}. \end{equation}

This corresponds to the range of energies that are well within the thermal spread of the least dense waterbags, but far outside the thermal spread of the densest ones. At ${\varepsilon (\boldsymbol {p}) \gg 1/\beta \eta _{\min }}$, the lower limit of the integral in (3.16) imposes an exponential cutoff on $\langle f \rangle$. The distribution function $N(\varepsilon )$ of particle energies corresponding to (3.16) is obtained by multiplying the mean phase-space density $\langle f \rangle (\boldsymbol {p})$ by the density of states:

(3.18)\begin{equation} N(\varepsilon) = g(\varepsilon)\langle f \rangle(\boldsymbol{p}) = \frac{n_{0}}{ \beta \varGamma(a+1) \eta_{\max} \varepsilon^{2}}\int_{\beta \eta_{\min}\varepsilon}^{{\beta \eta_{\max}\varepsilon}}\, \mathrm{d}x x^{1+a}G\left(\frac{x}{\beta \eta_{\max}\varepsilon} \right) {\rm e}^{{-}x}. \end{equation}

This shows that the non-degenerate Lynden-Bell equilibria express a natural power-law tail and, furthermore, that this power law is independent of the type of plasma system under consideration. Note that the origin of the $\varepsilon ^{-2}$ scaling found here is entirely different than the $\varepsilon ^{-2}$ arising from particle acceleration in shocks (Bell 1978).

Consider now what happens if the Gardner distribution does not have an exponential tail. Then $\rho (\eta )$ can be written as (3.10). Following all the same steps as before from (3.15) onwards, but using (3.10) in place of (3.8), one arrives at

(3.19)\begin{equation} N(\varepsilon) = \frac{n_{0}}{ \beta^{2-\delta} \varGamma(a+1) \eta_{\max}^{2-\delta} \varepsilon^{3-\delta}}\int_{\beta \eta_{\min}\varepsilon}^{{\beta \eta_{\max}\varepsilon}}\, \mathrm{d}x x^{2+ a -\delta}G\left(\frac{x}{\beta \eta_{\max}\varepsilon} \right) {\rm e}^{{-}x}. \end{equation}

Thus, the resulting Lynden-Bell equilibrium again displays a power-law tail. The power law's exponent is set by the particular value of $\delta$ in (3.10), which is related to the Gardner distribution of that Lynden-Bell equilibrium via (3.5), as explained at the end of § 3.2. For Gardner distributions that already have power-law tails, the resulting Lynden-Bell equilibria have shallower (ultra-‘hard’) power-law tails that strongly diverge in energy, giving the cutoff $\eta _{\min }$ pivotal importance. For Gardner distributions that decay faster than any exponential, the Lynden-Bell equilibria have ‘soft’ power-law tails, with total energy depending only very weakly on the cutoff.⁶

Presently, however, we shall return to the Lynden-Bell equilibria (3.18) arising from exponential Gardner distributions, for which we complete the calculation.

3.4. Calculation of $\beta$ and inevitability of partial degeneracy

Both (3.17) and (3.18) still depend on the as yet unknown parameter $\beta$, which, as well as fixing the energy, will determine the accuracy of the non-degeneracy approximation (3.12). To find $\beta$, we must compute the energy of our Lynden-Bell equilibrium (3.13) according to (2.5). Equivalently, from (3.15),

(3.20)\begin{equation} E = \int\mathrm{d}\boldsymbol{p} \varepsilon(\boldsymbol{p})\langle f \rangle(\boldsymbol{p}) = \frac{a+1}{\beta}\int_{\eta_{\min}}^{\eta_{\max}} \mathrm{d}\eta \rho(\eta). \end{equation}

Therefore, in the non-degenerate limit,

(3.21)\begin{equation} \beta = \frac{a+1}{E}\int_{\eta_{\min}}^{\eta_{\max}} \mathrm{d}\eta\rho(\eta). \end{equation}

We see that $\beta$ decreases with increasing total energy of the distribution function; this is natural if $\beta$ is viewed as an inverse thermodynamic temperature.

To see how this affects the underlying assumption (3.12) of non-degeneracy, we evaluate $D(\varepsilon )$ at $\varepsilon \to 0$, where the degeneracy effect is strongest, since $D(\varepsilon ) \leq D(0)$. Requiring $D(0) \ll 1$ gives the following condition on the total energy of the distribution, via (3.12), (3.14) and (3.21):

(3.22)\begin{equation} E \gg (a+1)\left[\int_{\eta_{\min}}^{\eta_{\max}} \mathrm{d}\eta \frac{\eta^{a+1}\rho(\eta)}{A\varGamma(a+1)} \right]^{1/(a+1)}\int_{\eta_{\min}}^{\eta_{\max}} \mathrm{d}\eta\rho(\eta). \end{equation}

Since $\rho (\eta )$ is a function only of the Gardner distribution, the right-hand side of (3.22) must scale with the energy density $E_{G}$ of the Gardner distribution that has the same waterbag content $\rho (\eta )$, but will always have $E_{G} \leq E$. For example, for the Gardner distribution (3.1), the right-hand side of (3.22) should be proportional to $n_{0}\varepsilon _{0}$. Thus, just like in Fermi–Dirac statistics, the non-degeneracy approximation becomes more accurate as the distribution's energy density $E$ begins to dwarf the energy density $E_{G}$ of the ground state. Note, however, that (3.22) contains an integral of the waterbag content $\rho (\eta )$ with no weighting by $\eta$, which, by (2.4) and (2.6), is the total phase volume occupied by non-empty waterbags. Since $\rho (\eta ) \propto \eta ^{-1}$ at small $\eta$, this will be large, depending, albeit logarithmically, on the minimum waterbag density $\eta _{\min }$. Indeed, for our example (3.1), (3.22) becomes

(3.23)\begin{equation} E \gg \frac{4}{3\sqrt{\rm \pi}}n_{0}\varepsilon_{0} \left(\ln\frac{\eta_{\max}}{\eta_{\min}} \right)^{3/2}\left[ 1 + \mathcal{O}\left(\frac{\eta_{\min}}{\eta_{\max}} \right) \right]. \end{equation}

This means that the condition (3.22) requires the energy of the distribution to be much greater not just than the energy of the corresponding Gardner distribution, but than the Gardner energy multiplied by a polylogarithmic function of $\eta _{\min }$. This is a manifestation of the fact that, as $\eta _{\min }$ is taken to zero, more of the phase space is pervaded by low-density waterbags, making true non-degeneracy harder to achieve. It is thus impossible to achieve a non-degenerate limit unless $\eta _{\min }$ is kept finite.

This is not the only place where the finiteness of the cutoff $\eta _{\min }$ has raged against the dying of the light. The same effect is manifest in the power law of $\varepsilon ^{-2}$ appearing in (3.18), which would have led to a logarithmically divergent mean particle energy were it not for the exponential cutoff at $\varepsilon \sim 1/\beta \eta _{\min }$. This is obvious in (3.20), where the integral of $\rho (\eta )$ has the same logarithmic divergence with $\eta _{\min } \to 0$ as it did in the right-hand side of (3.22). This then makes its way into the expression (3.21) for $\beta$, so, formally in the limit $\eta _{\min } \to 0$, $\beta \to \infty$ always!

Since we must keep $\eta _{\min }$ finite, it is of substantial importance to understand the physical significance of it and the extent to which one should be prepared to accept one's equilibrium's dependence on its value.

3.5. Physical meaning of minimum waterbag density

A Lynden-Bell equilibrium is essentially the thermal equilibrium of a collection of correlated blobs in phase space. This is to say that, inherent to the idea of computing the mean phase-space density, we have assumed that an exact phase-space density, i.e. a finite value of $\eta$, is a meaningful concept. But, of course, an exact phase-space density is a fiction, since a plasma is composed of many discrete particles, and a phase-space density is only an average occupation number of particles’ positions in phase space. The only sense in which an exact, continuous, phase-space density can be meaningful in a collisionless plasma then is if, within a small enough phase-space volume $\Delta \Gamma$, many particles can be considered to move as a collective entity: a waterbag. Then, on the scale of $\Delta \Gamma$, the system is composed of many waterbags with some ‘exact’ phase-space density, whereas on scales much larger than $\Delta \Gamma$, the system can attain a mean phase-space density. This ‘correlation volume’ provides a natural way to introduce the minimum non-zero phase-space density: clearly that should correspond to a single particle sitting in $\Delta \Gamma$, giving ${\eta _{\min } = \Delta \Gamma ^{-1}}$.

Determining the value of $\Delta \Gamma$, is, however, a non-trivial challenge (see, e.g. discussions in Kadomtsev & Pogutse 1970; Chavanis 2022; Ewart et al. 2022). Ewart et al. (2022) argued, on the grounds that any meaningful collisionless-relaxation rate must be smaller than the plasma frequency but larger than the rate at which collisions break phase-volume conservation, that a reasonable constraint to place on the correlation volume is

(3.24)\begin{equation} \Delta\Gamma \sim \frac{1}{\eta_{\mathrm{eff}}}(n_{0} \lambda_{D}^{3} )^{\alpha}, \quad \frac{2}{3} < \alpha < 1, \end{equation}

where $\eta _{\mathrm {eff}}$ is some typical phase-space density, which we can here estimate by $\eta _{\max }$, and $\lambda _{D}$ is the Debye length. This gives us an estimate for the minimum waterbag density

(3.25)\begin{equation} \eta_{\min} = \Delta\Gamma^{{-}1} \sim \eta_{\max} (n_{0} \lambda_{D}^{3} )^{-\alpha}. \end{equation}

Since the typical number of particles in a Debye sphere can be as large as $10^{6}$ to $10^{8}$ in collisionless plasma environments (such as the solar wind or interstellar gas: see, e.g. Ferrière 2019; Verscharen, Klein & Maruca 2019), the estimate (3.25) might seem damningly small. It is, in fact, ideal. The existence of a broad power-law tail in (3.18) required a scale separation between $\eta _{\max }$ and $\eta _{\min }$. The estimate (3.25) certainly provides this separation, tied to the plasma parameter. As for the breakdown of the non-degenerate approximation and the marginal divergence of the mean particle energy of a distribution with an $\varepsilon ^{-2}$ power law, we are saved by the fact that only the logarithm of the ratio $\eta _{\max }/\eta _{\min }$ will appear, which, while large, can only ever be in the range of 10–30. This is somewhat reminiscent of the situations in the conventional theory of Coulomb collisions in plasmas, where forcible introduction of a phase-space cutoff into the collision integral results in only a weak dependence on the value of this cutoff, via the so-called Coulomb logarithm (see, e.g. Helander & Sigmar 2005).

Nevertheless, the presence of the logarithmic divergence with $\eta _{\min }$ in the expression for $\beta$, signposted at the end of § 3.4, will make it difficult to satisfy the non-degeneracy approximation.⁷ There is good reason to suppose, however, that its breakdown may only be partial. We note that evaluating (3.12) at $\varepsilon = 0$ is tantamount to requesting non-degeneracy everywhere in the distribution. Since the degeneracy parameter $D(\varepsilon )$ decreases with increasing energy, it is reasonable to expect (and indeed we will see) some degeneracy at low energies, which gives way to non-degeneracy at higher energies, where our power-law tails will be recovered. As ever, the true solution lies on the cusp of asymptotic theory, and to go any further we must resort to numerical computation.

4. Numerical results: partially degenerate equilibria

In this section, we shall recover the analytically predicted power-law tail in (3.18) by solving the constraint equations (2.5) and (2.6) for the Lynden-Bell equilibria (3.4). The numerical scheme for this is documented in detail in Appendix B, amounting to an iterative method coupled to a one-dimensional root finder. For these numerical results, we have restricted ourselves to a three-dimensional, non-relativistic plasma with $\varepsilon (\boldsymbol {p}) = p^{2}/2m$, although we anticipate from § 3 that these results extend, qualitatively, to general regimes.

To capture a broad range of initial conditions, we consider a family of waterbag contents defined by

(4.1)\begin{equation} \rho(\eta) = \frac{4{\rm \pi} p_{\mathrm{th,}\sigma}^{3}}{\sigma \eta}\left(\ln \frac{\eta_{\mathrm{max,}\sigma}}{\eta} \right)^{(3-\sigma)/\sigma}, \quad \eta_{\min} < \eta < \eta_{\mathrm{max,}\sigma}, \end{equation}

with $\sigma > 0$. This defines the family of Gardner distributions

(4.2)\begin{equation} f_{\mathrm{G,}\sigma}(\boldsymbol{p}) = \begin{cases} \eta_{\max,\sigma}\, {\rm e}^{-(p/p_{\mathrm{th,}\sigma})^{\sigma}} & \text{for } p < p_{\mathrm{th,}\sigma}\left(\ln \displaystyle\frac{\eta_{\mathrm{max,}\sigma}}{\eta_{\min}} \right)^{1/\sigma}, \\ 0 & \text{for } p > p_{\mathrm{th,}\sigma}\left(\ln \displaystyle\frac{\eta_{\mathrm{max,}\sigma}}{\eta_{\min}} \right)^{1/\sigma}, \end{cases} \end{equation}

which, in the limit of $\eta _{\min }/\eta _{\mathrm {max,}\sigma } \to 0$, have particle densities $n_{0}$ and energy densities $E_{G}$ that satisfy

(4.3a,b)\begin{equation} n_{0} = \frac{4{\rm \pi}}{\sigma}\varGamma\left(\frac{3}{\sigma}\right)p_{\mathrm{th,}\sigma}^{3} \eta_{\mathrm{max,}\sigma}, \quad E_{G} = \frac{\varGamma(5/\sigma)}{\varGamma(3/\sigma)}n_{0} \frac{p_{\mathrm{th,}\sigma}^{2}}{2m}. \end{equation}

Since these Gardner distributions represent minimum-energy states, we will be able to scan in the energy density of the system for all $E \geq E_{G}$, imagining some initial distribution of particles, with waterbag content $\rho (\eta )$ and energy density $E_{G}$, being accelerated to the energy density $E$, and then seeking a maximum-entropy state (see § 3.1). We can also vary $\sigma$ in order to see the effects of the shape of the underlying Gardner distribution on the resulting Lynden-Bell equilibria.

4.1. Degrees of degeneracy

Let us first scan in $\eta _{\mathrm {max,}\sigma }/\eta _{\min }$ in order to show that we can indeed recover the fully non-degenerate limit solved in § 3. Figure 3 shows the result of such a scan for $\sigma = 2$ and $E = 10E_{G}$. By comparing the exact (numerically calculated) fugacity with the theoretical prediction (3.14) obtained in the absence of phase-space degeneracy, we see that the agreement is nearly perfect when $\eta _{\max }$ is close to $\eta _{\min }$, e.g. when $\eta _{\max }/\eta _{\min } = 10$. This is as expected, since the non-degenerate limit is valid when (3.22) holds, which it does, as can be confirmed from the solid red colour in panel (b), showing the probability $D(\varepsilon )/[1 + D(\varepsilon )]$ (with $D(\varepsilon )$ defined in (3.12)) that a region of phase space is occupied. However, at $\eta _{\max }/\eta _{\min } = 10$, there is an insufficient range of waterbag levels to achieve the scale separation (3.17) necessary to resolve a power-law tail in energies. To see a power-law tail, one must increase $\eta _{\max }/\eta _{\min }$ to higher values, but this comes at the price of increasing the degeneracy of the system and, hence, undermining the asymptotic regime in which the tail was derived in the first place. For the values of $\eta _{\max }/\eta _{\min }$ that we argued in § 3.5 to be realistic, the system becomes strongly degenerate at low energies, which can again be seen from the solid colours in panel (b). In spite of this, at high energies, the degeneracy falls away, and, correspondingly, $F(\eta )$ at low $\eta$ still agrees well with the non-degenerate approximation (3.14). All this conspires to ensure that, even formally outside the non-degenerate limit, the power-law tail $N(\varepsilon ) \propto \varepsilon ^{-2}$ is still manifestly present.

Figure 3. Numerically computed Lynden-Bell equilibria for a range of $\eta _{\mathrm {max,}\sigma }/\eta _{\min }$ and with $\rho (\eta )$ given by (4.1) with $\sigma = 2$. The energy density is equal to $10 E_{G}$ in all cases. (a) The numerically computed fugacity $F(\eta )$ (solid lines) compared with the analytical solution (3.14) obtained in the non-degenerate limit (dashed lines). (b) The resulting distributions $N(\varepsilon )$ of particle energies, with the universal power law $\propto \varepsilon ^{-2}$ shown for reference, cf. (3.18). Overplotted in solid colour (with the value range shown on the right) is the level of degeneracy $D(\varepsilon )/[1+ D(\varepsilon )]$ (the probability that a given energy is occupied by a non-empty waterbag) as a function of energy; $D(\varepsilon )$ is defined in (3.12).

4.2. Energisation of particles and power-law tails

Let us now scan in the energies densities $E$ of the distribution and again look for power-law tails and assess the effect of degeneracy. Figure 4 shows the results of such a scan, again with ${\rho (\eta )}$ specified by (4.1) with $\sigma = 2$. Plotted underneath the mean phase-space densities are the contributions from a number of finite ranges of exact phase-space density defined by

(4.4)\begin{equation} \langle f \rangle_{\eta_{1} < \eta < \eta_{2}}(\boldsymbol{p}) = \int_{\eta_{1}}^{\eta_{2}} \, \mathrm{d}\eta \eta P(\boldsymbol{p},\eta). \end{equation}

As one would anticipate, when $E$ is only slightly larger than $E_{G}$, the effects of phase-space degeneracy are most prominent: the densest portions of the phase space clog up the lowest energies, forcing less dense portions to higher energies. As the total energy density is increased, we see that the contribution from each range of waterbags spreads out. This is because the increased energy allows dense portions of phase space to be promoted to larger energies, making room at lower energies for less dense portions of phase space to fill.

Figure 4. Numerically computed Lynden-Bell equilibria for a range of energy densities $E$ (in multiples of the energy density $E_{G}$ of the underlying Gardner distribution) with $\rho (\eta )$ given by (4.1) with $\sigma = 2$ and $\eta _{\max }/\eta _{\min } = 10^{6}$. In each plot, the dashed line is the mean phase-space density, while the underplotted solid lines are the contributions from four distinct ranges of exact phase-space density as functions of energy. Note that, while the exact phase-space densities have been grouped into four, this is not the same as solving for a four-waterbag Lynden-Bell equilibrium, as each grouping is still composed of a continuum of waterbags.

While the solutions plotted in figure 4 might appear qualitatively similar to the Lynden-Bell equilibria obtained in numerical experiments with a small discrete number of level sets (see, e.g. Assllani et al. 2012; Ewart et al. 2022), this hides key universal features of systems with a continuum of level sets. To showcase this universality, in figure 5, we plot the Lynden-Bell equilibria for two different waterbag contents, $\sigma = 2$ and $\sigma = 8$ in (4.1), and a range of energy densities $E$. The $\varepsilon ^{-2}$ power-law tails of these equilibria are immediately apparent, as predicted in (3.18).

Figure 5. Numerically computed Lynden-Bell equilibria with waterbag content given by the $\sigma = 2$ (top) and $\sigma = 8$ (bottom) cases of (4.1), $\eta _{\mathrm {max,}\sigma } / \eta _{\min } = 10^{6}$. (a) The phase-space densities shown in linear scale, (b) the corresponding distributions of particle energies in logarithmic scale, for a range of ratios of $E/E_{G}$. The small deviations from the $\varepsilon ^{-2}$ tail can be attributed to the logarithmic corrections arising from the $x$ integral in (3.18).

Figure 5 shows how these power-law tails become more prominent as one adds more energy to the Gardner distribution. At $E = E_{G}$, one has a highly non-universal Gardner equilibrium. As a small amount of energy $E - E_{G} \ll E_{G}$ is added to $E_{G}$, the mean phase-space density at low energies is largely unaffected, while the power-law tail grows from the lowest-density waterbags. Thus, for energies close to the energy of the underlying Gardner distribution, the Lynden-Bell equilibria have a ‘core-halo’ structure: the ‘core’, which has energy density $\sim E_{G}$, is comprised of dense waterbags, which the system does not have sufficient energy to excite, whereas the ‘halo’ is the tail comprised of those less dense waterbags that are capable of sampling a larger portion of phase space and thus arrange themselves into a universal $\varepsilon ^{-2}$ power law, containing the excess energy density $\sim E - E_{G}$. As the energy of the distribution is further increased, more waterbags have sufficient energy to sample a larger range of phase space, the halo continues to eat into the core, but both thermally broaden. At $E \gg E_{G}$, the asymptotically non-degenerate solution (3.18) with a power law in the energy range (3.17) suggests that the system strives for a state in which the halo has much more energy than the core. In this limit, one expects the transition between the core and halo to occur at $\varepsilon \sim 1/\beta \eta _{\max }$. Since the halo should be exponentially suppressed at $\varepsilon \gtrsim 1/\beta \eta _{\min }$, one can compute the ratio of the core energy to the halo energy. Owing to the $\varepsilon ^{-2}$ tail, this ratio will be proportional to $\ln (\eta _{\max }/\eta _{\min })$ raised to some power, which depends on the specific functional form of the Gardner distribution. Whatever this power is, the vast majority of the total energy will be contained in the universal power-law tail for any system whose energy is much larger than that of its Gardner state.

5. Conclusion

5.1. Summary

The Lynden-Bell (1967) equilibria are the natural maximum-entropy states for systems, such as a plasma described by the collisionless Vlasov equation (1.1), which conserve not only density, momentum, and energy, but also an infinite family of further invariants (2.1). These additional invariants are due to the conservation of phase volume, encoded by the ‘waterbag content’ function $\rho (\eta )$ given by (2.7) (equivalent to the Casimir invariants (3.7)), which measures the amount of phase volume where the exact phase-space density takes the value $\eta$ (a ‘waterbag’), per unit $\eta$. Maximising entropy subject to all these conservation laws then gives the mean phase-space density (2.2) in the form of the Lynden-Bell equilibrium (2.9) coupled with the constraints (2.5) and (2.6). In this paper, we have solved these constraint equations numerically in the general case, as well as analytically in a tractable limit (which turned out to be the practically relevant one). We have been able to show that, despite their apparent dependence on non-universal initial conditions, Lynden-Bell equilibria generically exhibit power-law tails, and, in particular, that a broad class of initial conditions will give rise to the distribution of particles’ energies scaling as $\varepsilon ^{-2}$ at high $\varepsilon$.

To study the Lynden-Bell equilibria systematically, we first considered what values the invariants of the system, the energy density $E$ and waterbag content $\rho (\eta )$, could take. This led us to the concept of a Gardner distribution function, which is any monotonically decreasing function of the particle's energy. In § 3.1, we argued that to each possible initial condition one could assign a unique Gardner distribution, with the same waterbag content (and, therefore, the same Casimir invariants) as the initial condition, but a different energy, the Gardner distribution having, by definition, the lowest possible energy of all distributions with a given waterbag content $\rho (\eta )$. The Lynden-Bell equilibria at higher energies and the same $\rho (\eta )$ can, therefore, be thought of as excited states of this Gardner distribution. In § 3.2, we argued that the typical $\rho (\eta )$ would have a fairly generic power-law form at low $\eta$, and, in particular, that it would scale as $\eta ^{-1}$ for a wide class of initial conditions (see (3.8)). In § 3.3, we were able to find the Lynden-Bell equilibria analytically provided the energy of the system was sufficiently large for the competition of waterbags for phase space to be ignorable. In this ‘non-degenerate’ limit, particles belonging to each waterbag arrange themselves into a separate Maxwellian equilibrium with an effective ‘temperature’ inversely proportional to the phase-space density of that waterbag: denser portions of phase space are energetically costlier to move to higher energies. The resulting mean phase-space density (3.15) was found by integrating the contributions of all waterbags, each with their own Maxwellian distribution weighted by the amount of phase space which that waterbag occupied. Since the amount of each waterbag had a universal form (3.8), this gave rise to a universal power-law tail (3.18) scaling as $\varepsilon ^{-2}$ at high particle energies.

However, the non-degenerate limit required the system's energy to be asymptotically larger than the energy of the corresponding Gardner distribution. Formally, this turned out to be a very stringent limitation, and indeed one that could hardly ever be strictly fulfilled. Our analytical results were rescued by the argument, confirmed by the numerical solutions presented in § 4, that the effects of phase-space degeneracy were confined to the low-energy part of the distribution. The universal $\varepsilon ^{-2}$ tail was numerically confirmed to be a robust feature of the generic Lynden-Bell equilibria. As well as ascertaining that a range of different initial conditions (4.1) gave rise to the same power-law tail, the numerical solution also showed how this power-law tail formed. We found (figure 5) that at energies comparable to the Gardner energy, the Lynden-Bell equilibria had a ‘core-halo’ structure. The halo, consisting of the $\varepsilon ^{-2}$ tail, was formed from low-density waterbags, which had sufficient energy to explore large portions of phase space, while the non-universal core was made up of denser waterbags, which did not have sufficient energies to be excited. As the energy of the Lynden-Bell equilibrium was increased, the halo ate its way into the core, making the distribution less and less degenerate and more universal in its shape. This behaviour is perhaps reminiscent of the measurements of the ‘non-thermal fraction’ of particles in the solar wind (see, e.g. Pierrard & Lazar 2010; Oka et al. 2015).

5.2. Limitations and applications

The Lynden-Bell statistical mechanics presents an attractive scenario for the universal generation of power-law tails, which could perhaps offer insight into the distribution of particles in such astrophysical systems as cosmic rays and the solar wind. Indeed, a power law of $\varepsilon ^{-2}$ in energy is roughly consistent with the observed power laws typical in the quiet-time solar wind (e.g. Gloeckler et al. 2008; Fisk & Gloeckler 2014; Yang et al. 2020) and close to, but distinct from, the inferred value for cosmic-ray sources (e.g. Ormes & Freier 1978; Reichherzer et al. 2021). However, the universality of this predicted power law fails to capture the wide range of power laws seen both in numerical simulations (e.g. Sironi & Spitkovsky 2014; Werner & Uzdensky 2017; Zhdankin et al. 2017) and observationally nearer to the Sun (see Oka et al. (2018) and references therein). Such systems are usually turbulent, possibly inhomogeneous, and invariably magnetised. In contrast, the theory that we have proposed for a universal power-law tail assumes a homogeneous system, in which all the fluctuating field's energy has decayed to a negligible fraction of the total energy. It is therefore an intriguing question whether our theory can be adjusted to apply to these cases. It is this question that we will address in this section, speculatively.

One common limitation in the applicability of Lynden-Bell's theory is that the fluctuating fields may decay away faster than the equilibrium state is reached – this is a well-understood feature in galactic dynamics, where it is referred to as ‘incomplete relaxation’ (see Chavanis (2006b) and references therein). Here we have ignored such a possibility, assuming effectively that there will always be a sufficient amount of fluctuations to see the plasma through to its maximum-entropy state.

Perhaps an even more pressing concern is the possible reliance of the theory on precise phase-volume conservation. The validity of the Lynden-Bell statistics rests on the assumption that such an equilibrium can be reached long before the conservation of phase volume is broken. Conventional (linear) estimates for the time scale on which it would be broken, scaling as an inverse fractional power of the true Coulomb-collision frequency (Su & Oberman 1968), indicate that such a relaxation should be possible. However, recent progress (Beraldo e Silva et al. 2017; Zhdankin 2022a; Nastac et al. 2023) has shown that in turbulent systems, the conservation of phase volume may be broken on fast, collision-frequency-independent time scales. While it is possible that this is damning evidence against the existence (or persistence) of Lynden-Bell equilibria, it is also possible that it points to interesting interplay between the Lynden-Bell relaxation and collisional effects. For instance, one can imagine the possibility that the effect of collisions is to evolve $\rho (\eta )$ without immediately pinning the system to a Maxwellian equilibrium. This can happen if collisions are already acting to erase small-scale phase-space structure of the exact phase-space density $f$ but not yet to change its mean $\langle f \rangle$ directly. In such a situation, one's aim would be to compute the evolution of $\rho (\eta )$ or, equivalently, of the underlying Gardner ground state $f_{G}(\varepsilon )$. In particular, one can imagine a regime in which the underlying Gardner distribution evolves slower than the Lynden-Bell equilibrium is reached, causing the distribution to go through a sequence of Lynden-Bell equilibria. If the effect of collisions is to smooth the fluctuations $f - \langle f \rangle$ diffusively, it is clear that this can only cause the energy $E_{G}$ of the Gardner distribution to increase or remain constant (cf. Kolmes & Fisch 2020). One would therefore expect that, as the system evolves through a sequence of Lynden-Bell equilibria, it will steadily become more degenerate as $E_{G}$ approaches the system's energy $E$, with the core of the Lynden-Bell distribution eating into its tail (halo; see § 4.2). This partially collisionless evolution would finally freeze once all fluctuations are diffused away, leaving a degenerate Lynden-Bell equilibrium doomed to further gradual erosion by weak Coulomb collisions on $\langle f \rangle$. Since we have argued that an $\varepsilon ^{-2}$ tail is generic for any Lynden-Bell equilibrium whose underlying Gardner distribution has an exponential tail, this scenario suggests that the universal power-law tail could persist as long as $E_{G} < E$, despite the breaking of phase-volume conservation due to weak collisionality.

Finally, there is the question of whether any of this formalism can be ported smoothly to magnetised equilibria. Here, the most obvious straw to grasp at is that there is no guarantee that the invariants (2.5) and (2.6) are the only invariants respected on the relaxation time scale. As previously mentioned, in drift-kinetic plasmas, each particle conserves its magnetic moment $\mu _{b}$. This implies a new conserved function

(5.1)\begin{equation} \frac{2{\rm \pi}}{V}\int\mathrm{d}\boldsymbol{r}\int\mathrm{d}v_{{\parallel}}B(\boldsymbol{r})\delta(\eta - f(\boldsymbol{r},v_{{\parallel}},\mu_{b})) = \rho(\eta, \mu_{b}), \end{equation}

which would supersede the conservation of the now mundane $\rho (\eta )$. While some Gardner distributions have been studied for such systems (see Helander 2020; Mackenbach, Proll & Helander 2022), the Lynden-Bell equilibria in them are unexplored and may contain a wealth of interesting physics. This said, in turbulent systems, the conservation law (5.1) may be just as fragile as the conservation of phase volume. Indeed it has been suggested that the breaking of adiabatic invariance may be essential to understanding the transport properties of non-thermal particles (see Ruszkowski & Pfrommer (2023) and references therein).

From the previous discussion it becomes clear that perhaps the most relevant limitation of Lynden-Bell's statistical mechanics – or indeed of any equilibrium statistical mechanics – in application to observed plasma phenomena is that much of real plasma dynamics are out of equilibrium in a physically essential way: any local relaxation processes, collisional or collisionless, tend to have to be taken into account alongside various ‘sources’ and ‘sinks’ of particles and/or energy, e.g. the energisation and escape of cosmic rays (Schlickeiser 1989; Chandran 2000; Becker Tjus & Merten 2020; Hopkins et al. 2022; Kempski & Quataert 2022), the turbulent heating and radiative cooling of the intracluster medium (e.g. Zhuravleva et al. (2014) and references therein) or accretion-disc plasmas (e.g. Lesur (2021), Kawazura et al. (2022) and references therein), a veritable zoo of such processes in the solar wind (e.g. Verscharen et al. 2019; Chen et al. 2020) and the Earth's magnetosphere (e.g. Lucek et al. 2005), the birth of energetic $\alpha$-particles in fusion reactions and their subsequent slowing down and escape from confined plasmas (e.g. Helander & Sigmar 2005; Mailloux et al. 2022), etc. In plasmas where Coulomb collisions can be assumed to relax the particle distribution quickly to a local Maxwellian, we have a robust analytical framework for handling all these non-equilibrium processes in terms of the evolution of the density, momentum and temperature of that Maxwellian and a separation of the dynamics into that evolution plus the turbulence of small fluctuations around the local equilibrium (e.g. Schekochihin et al. 2009; Abel et al. 2013). In collisionless plasmas, such a framework is lacking as both the turbulence and the nature of the underlying equilibrium are mysterious and indeed it is not even guaranteed that they can be understood without detailed reference to each system's particular initial circumstances. If Lynden-Bell's statistical mechanics proves to be a viable collisionless substitute for Maxwell's, a path could be charted towards a theory of the dynamics and thermodynamics of collisionless plasmas possessing a modicum of universality.