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 Reference Gardner1963; Helander Reference Helander2017). 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 Reference Dodin and Fisch2005; Kolmes & Fisch Reference Kolmes and Fisch2020; Kolmes, Helander & Fisch Reference Kolmes, Helander and Fisch2020). Gardner distributions can be viewed as ‘ground states’ associated with a given waterbag content (e.g. Helander Reference Helander2017), 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 Reference Chavanis2006a), 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 Reference Dodin and Fisch2005; Helander Reference Helander2017). 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 Reference Zhdankin2022a):

(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$.Footnote ⁵ 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 Reference Bell1978).

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.Footnote ⁶

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 Reference Kadomtsev and Pogutse1970; Chavanis Reference Chavanis2022; Ewart et al. Reference Ewart, Brown, Adkins and Schekochihin2022). Ewart et al. (Reference Ewart, Brown, Adkins and Schekochihin2022) 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 Reference Ferrière2019; Verscharen, Klein & Maruca Reference Verscharen, Klein and Maruca2019), 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 Reference Helander and Sigmar2005).

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.Footnote ⁷ 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.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. Reference Assllani, Fanelli, Turchi, Carletti and Leoncini2012; Ewart et al. Reference Ewart, Brown, Adkins and Schekochihin2022), 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.