2.1. Linear theory for $f_{k\alpha }$

To solve for the evolution of $f_{\boldsymbol {k}\alpha }$ from (2.5) and (2.6), we follow Landau (Reference Landau1936) and introduce the Laplace transform

(2.7)\begin{equation} \hat{\varphi}(p) = \int_{0}^{\infty}{\rm e}^{{-}pt}\varphi(t)\,\mathrm{d}t \,, \end{equation}

and similarly for $f_{\boldsymbol {k}\alpha }(t)$. We now take the Laplace transform of (2.5) and (2.6) to get

(2.8)\begin{gather} \hat{f}_{\boldsymbol{k}\alpha}(p) = {\rm i}\frac{q_{\alpha}}{m_{\alpha}}\frac{\hat{\varphi}_{\boldsymbol{k}}(p)}{p+{\rm i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}}\boldsymbol{k}\boldsymbol{\cdot}\frac{\partial f_{0\alpha}}{\partial \boldsymbol{v}} + \hat{h}_{\boldsymbol{k}\alpha}(p), \end{gather}

(2.9)\begin{gather} \hat{\varphi}_{\boldsymbol{k}}(p) = \sum_{\alpha'}\frac{4{\rm \pi} q_{\alpha'}}{k^{2}\epsilon_{\boldsymbol{k}}(p)}\int\,\mathrm{d}\boldsymbol{v}' \, \hat{h}_{\boldsymbol{k}\alpha'}(p), \end{gather}

where the dielectric function has emerged, defined by

(2.10)\begin{equation} \epsilon_{\boldsymbol{k}}(p) = 1 - {\rm i}\sum_{\alpha'}\frac{4{\rm \pi} q_{\alpha'}^{2}}{m_{\alpha'}k^{2}}\int\,\mathrm{d}\boldsymbol{v}' \,\frac{1}{p+{\rm i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}'}\boldsymbol{k}\boldsymbol{\cdot}\frac{\partial f_{0\alpha'}}{\partial \boldsymbol{v}'}, \end{equation}

while information about the initial condition enters via

(2.11)\begin{equation} \hat{h}_{\boldsymbol{k}\alpha}(p) = \frac{g_{\boldsymbol{k}\alpha}(\boldsymbol{v})}{p+{\rm i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}}. \end{equation}

where $g_{\boldsymbol {k}\alpha }(\boldsymbol {v}) = f_{\boldsymbol {k}\alpha }(t=0,\boldsymbol {v})$.

The time-dependent solution is given by the inverse Laplace transforms:

(2.12)\begin{gather} \varphi_{\boldsymbol{k}}(t) = \frac{1}{2{\rm \pi} {\rm i} }\int_{-{\rm i}\infty + \sigma}^{{\rm i}\infty + \sigma}\,\mathrm{d}p \, {\rm e}^{pt}\hat{\varphi}_{\boldsymbol{k}}(p), \end{gather}

(2.13)\begin{gather} f_{\boldsymbol{k}\alpha}(t) = \int_{-{\rm i}\infty + \sigma}^{{\rm i}\infty + \sigma}\frac{\mathrm{d}p \,}{2{\rm \pi} {\rm i} } {\rm e}^{pt}\left[{\rm i}\frac{q_{\alpha}}{m_{\alpha}}\frac{\hat{\varphi}_{\boldsymbol{k}}(p)}{p+{\rm i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}}\boldsymbol{k}\boldsymbol{\cdot}\frac{\partial f_{0\alpha}}{\partial \boldsymbol{v}} + \hat{h}_{\boldsymbol{k}\alpha}(p) \right], \end{gather}

where $\sigma$ must be chosen so that for all $p$ with $\mathrm {Re}(p)> \sigma$ the integrands are analytic functions.

2.2. Quasilinear evolution of $f_{0\alpha }$

Having computed $\varphi (t)$ and $f_{\boldsymbol {k}\alpha }(t)$ we are now in a position to determine the evolution of $f_{0\alpha }$. Substituting (2.12) and (2.13) into (2.4) we obtain the earliest form of our collision integral for the evolution of $f_{0\alpha}$:

(2.14)\begin{align} \frac{\partial f_{0\alpha}}{\partial t} ={-}\frac{\partial }{\partial \boldsymbol{v}}\boldsymbol{\cdot} \frac{q_{\alpha}}{m_{\alpha}}\sum_{\boldsymbol{k}}\boldsymbol{k}\mathrm{Im}\iint\frac{\mathrm{d}p \,\,\mathrm{d}p' \,}{(2{\rm \pi})^{2}} {\rm e}^{(p+p')t}\left[{\rm i}\frac{q_{\alpha}}{m_{\alpha}}\frac{\langle \hat{\varphi}_{\boldsymbol{k}}(p)\hat{\varphi}^{*}_{\boldsymbol{k}}(p'^{*}) \rangle}{p + {\rm i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}}\boldsymbol{k}\boldsymbol{\cdot}\frac{\partial f_{0\alpha}}{\partial \boldsymbol{v}}\right.\nonumber\\ \left.\quad\quad +\, \langle \hat{h}_{\boldsymbol{k}\alpha}(p,\boldsymbol{v})\hat{\varphi}^{*}_{\boldsymbol{k}}(p'^{*}) \rangle \vphantom{\frac{q_{\alpha}}{m_{\alpha}}}\right]. \end{align}

Note that, to avoid confusion in (2.14), both inverse Laplace transforms are taken along the contour running from $-{\rm i}\infty +\sigma$ to ${\rm i}\infty + \sigma$. Since the $p'$ contour in (2.14) comes from a complex conjugation, this results in an overall minus sign appearing.

Presently, (2.14) is still in need of information, as it still depends on the averages of $\hat {\varphi }_{\boldsymbol {k}}$ and $\hat {h}_{\boldsymbol {k}}$. The eventual aim is to ‘close’ this collision integral: to replace these correlators with expressions involving only the mean distribution $f_{0\alpha }$. This will give a differential equation for each $f_{0\alpha }$ in terms of other $f_{0\alpha '}$ only, which, in principle, can then be solved. Such a closure will come from some model of the averages that we have yet to compute, and will be enabled by assumptions introduced in §§ 2.5 and 3.

First we will manipulate (2.14) into a more agreeable form by rewriting the averages $\langle \hat {\varphi }_{\boldsymbol {k}}(p)\hat {\varphi }^{*}_{\boldsymbol {k}}(p'^{*}) \rangle$ and $\langle \hat {h}_{\boldsymbol {k}\alpha }(p,\boldsymbol {v})\hat {\varphi }^{*}_{\boldsymbol {k}}(p'^{*}) \rangle$ as correlators solely of $\hat {h}_{\boldsymbol {k}\alpha }$. For the first term in (2.14), we get, using (2.9),

(2.15)\begin{align} {\rm i}\frac{q_{\alpha}^{2}}{m_{\alpha}^{2}}& \frac{\langle \varphi_{\boldsymbol{k}}(p)\varphi^{*}_{\boldsymbol{k}}(p'^{*}) \rangle}{p+ {\rm i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}}\boldsymbol{k}\boldsymbol{\cdot}\frac{\partial f_{0\alpha}}{\partial \boldsymbol{v}} \nonumber= \\ & {\rm i}\frac{q_{\alpha}^{2}}{m_{\alpha}^{2}}\sum_{\alpha'\alpha''}\frac{16{\rm \pi}^{2}q_{\alpha'}q_{\alpha''}}{k^{4}\epsilon_{\boldsymbol{k}}(p)\epsilon^{*}_{\boldsymbol{k}}(p'^{*})} \iint \,\mathrm{d}\boldsymbol{v}' \,\,\mathrm{d}\boldsymbol{v}'' \,\frac{\langle \hat{h}_{\boldsymbol{k}\alpha''}(p,\boldsymbol{v}'')\hat{h}^{*}_{\boldsymbol{k}\alpha'}(p'^{*},\boldsymbol{v}') \rangle}{p+{\rm i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}}\boldsymbol{k}\boldsymbol{\cdot}\frac{\partial f_{0\alpha}}{\partial \boldsymbol{v}}. \end{align}

The second term in (2.14), again via (2.9), becomes

(2.16)\begin{align} \frac{q_{\alpha}}{m_{\alpha}}\langle \hat{h}_{\boldsymbol{k}\alpha}(p)\varphi_{\boldsymbol{k}}^{*}(p'^{*}) \rangle & = \sum_{\alpha'}\frac{4{\rm \pi} q_{\alpha}q_{\alpha'}}{m_{\alpha}k^{2}\epsilon_{\boldsymbol{k}}^{*}(p'^{*})}\int\,\mathrm{d}\boldsymbol{v} \,\langle \hat{h}_{\boldsymbol{k}\alpha}(p,\boldsymbol{v})\hat{h}^{*}_{\boldsymbol{k}\alpha'}(p'^{*},\boldsymbol{v}') \rangle \nonumber\\ & = \sum_{\alpha'}\frac{4{\rm \pi} q_{\alpha}q_{\alpha'}}{m_{\alpha}k^{2}\epsilon_{\boldsymbol{k}}(p)\epsilon^{*}_{\boldsymbol{k}}(p'^{*})}\int\,\mathrm{d}\boldsymbol{v}' \,\left[\vphantom{\sum_{\alpha'}\frac{4{\rm \pi} q_{\alpha}q_{\alpha'}}{m_{\alpha}k^{2}\epsilon_{\boldsymbol{k}}(p)\epsilon^{*}_{\boldsymbol{k}}(p'^{*})}}\langle \hat{h}_{\boldsymbol{k}\alpha}(p,\boldsymbol{v})\hat{h}^{*}_{\boldsymbol{k}\alpha'}(p'^{*},\boldsymbol{v}') \rangle \right.\nonumber\\ & \left.\quad-\,{\rm i} \sum_{\alpha''}\frac{4{\rm \pi} q_{\alpha''}^{2}}{m_{\alpha''}k^{2}}\int\,\mathrm{d}\boldsymbol{v}'' \,\frac{\langle \hat{h}_{\boldsymbol{k}\alpha}(p,\boldsymbol{v})\hat{h}^{*}_{\boldsymbol{k}\alpha'}(p'^{*},\boldsymbol{v}') \rangle}{p + {\rm i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}''}\boldsymbol{k}\boldsymbol{\cdot}\frac{\partial f_{0\alpha''}}{\partial \boldsymbol{v}''}\right]. \end{align}

At the last step, seemingly gratuitously, we multiplied and divided by the dielectric function (2.10). This will prove to be a useful tactic, as it separates this correlator into two parts: the second, which bears strong resemblance to (2.15), and the first, which will later vanish under certain assumptions. Since this first term is destined to vanish, we will continue as though it has already done so and confirm its disappearance at the end of § 2.5. The full evolution equation for the mean distribution function of species $\alpha$ is then

(2.17)\begin{gather} \frac{\partial f_{0\alpha}}{\partial t} ={-}\frac{q_{\alpha}}{m_{\alpha}}\frac{\partial }{\partial \boldsymbol{v}}\boldsymbol{\cdot}\sum_{\boldsymbol{k}}\frac{\boldsymbol{k}\boldsymbol{k}}{k^{4}}\boldsymbol{\cdot}\mathrm{Re}\iint\frac{\mathrm{d}p \,\,\mathrm{d}p' \,}{(2{\rm \pi})^{2}} {\rm e}^{(p+p')t}\sum_{\alpha'\alpha''}\frac{16{\rm \pi}^{2}q_{\alpha'}q_{\alpha''}}{\epsilon_{\boldsymbol{k}}(p)\epsilon^{*}_{\boldsymbol{k}}(p'^{*})} \iint \,\mathrm{d}\boldsymbol{v}' \,\,\mathrm{d}\boldsymbol{v}'' \,\nonumber\\ \left[\frac{\langle \hat{h}_{\boldsymbol{k}\alpha''}(p,\boldsymbol{v}'')\hat{h}^{*}_{\boldsymbol{k}\alpha'}(p'^{*},\boldsymbol{v}') \rangle}{p+{\rm i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}} \frac{q_{\alpha}}{m_{\alpha}}\frac{\partial f_{0\alpha}}{\partial \boldsymbol{v}} - \frac{\langle \hat{h}_{\boldsymbol{k}\alpha}(p,\boldsymbol{v})\hat{h}^{*}_{\boldsymbol{k}\alpha'}(p'^{*},\boldsymbol{v}') \rangle}{p+{\rm i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}''} \frac{q_{\alpha''}}{m_{\alpha''}}\frac{\partial f_{0\alpha''}}{\partial \boldsymbol{v}''} \right]. \end{gather}

This can be written more compactly as

(2.18)\begin{equation} \frac{\partial f_{0\alpha}}{\partial t} = \frac{\partial }{\partial \boldsymbol{v}}\boldsymbol{\cdot}\sum_{\alpha''}\int\,\mathrm{d}\boldsymbol{v}'' \,\left[\boldsymbol{\mathsf{D}}^{\alpha}_{\alpha\alpha''}(\boldsymbol{v},\boldsymbol{v}'')\boldsymbol{\cdot}\frac{\partial f_{0\alpha}}{\partial \boldsymbol{v}} - \boldsymbol{\mathsf{D}}^{\alpha}_{\alpha''\alpha}(\boldsymbol{v}'',\boldsymbol{v})\boldsymbol{\cdot}\frac{\partial f_{0\alpha''}}{\partial \boldsymbol{v}''}\right], \end{equation}

where the ‘diffusion kernel’ is

(2.19)\begin{align} \boldsymbol{\mathsf{D}}^{\alpha}_{\mu \nu}(\boldsymbol{w},\boldsymbol{v}) & ={-}\sum_{\nu'}\frac{16{\rm \pi}^{2}q_{\mu}^{2}q_{\nu}q_{\nu'}}{m_{\alpha}m_{\mu}}\mathrm{Re}\sum_{\boldsymbol{k}}\frac{\boldsymbol{k}\boldsymbol{k}}{k^{4}}\iint\frac{\mathrm{d}p \,\,\mathrm{d}p' \,}{(2{\rm \pi})^{2}} \frac{{\rm e}^{(p + p')t}}{\epsilon_{\boldsymbol{k}}(p)\epsilon^{*}_{\boldsymbol{k}}(p'^{*})}\nonumber\\& \quad \times\int\,\mathrm{d}\boldsymbol{v}' \,\frac{\langle \hat{h}_{\boldsymbol{k}\nu}(p,\boldsymbol{v})\hat{h}^{*}_{\boldsymbol{k}\nu'}(p'^{*},\boldsymbol{v}') \rangle}{p + {\rm i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{w}}. \end{align}

2.3. Simplification of the diffusion kernel

Even before we make any assumptions about the nature of the initial condition to decompose the averages, it is possible to simplify the diffusion kernel by appealing to the separation of timescales between the mean and the fluctuations. To do so, we will carry out the $p$ and $p'$ integrals. First we rewrite (2.19) as follows:

(2.20)\begin{equation} \boldsymbol{\mathsf{D}}^{\alpha}_{\mu\nu}(\boldsymbol{w},\boldsymbol{v}) = \sum_{\nu'}\frac{16{\rm \pi}^{2}q_{\mu}^{2}q_{\nu}q_{\nu'}}{m_{\alpha}m_{\mu}}\mathrm{Re}\sum_{\boldsymbol{k}}\frac{\boldsymbol{k}\boldsymbol{k}}{k^{4}}\int \,\mathrm{d}\boldsymbol{v} \,'\langle g_{\boldsymbol{k}\nu}(\boldsymbol{v})g^{*}_{\boldsymbol{k}\nu'}(\boldsymbol{v}') \rangle I_{\boldsymbol{k}}(\boldsymbol{v},\boldsymbol{v}',\boldsymbol{w}), \end{equation}

where

(2.21)\begin{equation} I_{\boldsymbol{k}}(\boldsymbol{v},\boldsymbol{v}',\boldsymbol{w}) ={-}\iint\frac{\mathrm{d}p \,\,\mathrm{d}p' \,}{(2{\rm \pi})^{2}}\frac{{\rm e}^{(p+p')t}}{\epsilon_{\boldsymbol{k}}(p)\epsilon^{*}_{\boldsymbol{k}}(p'^{*})}\frac{1}{(p+ {\rm i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v})(p+{\rm i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{w})(p'-{\rm i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}')}. \end{equation}

Since this is a contour integral of a holomorphic function [note that $\epsilon _{\boldsymbol {k}}^{*}(p'^{*})$ is a holomorphic function of $p'$], these integrals can be computed by deforming the $p$ and $p'$ contours far into the left-hand plane, where the real part of ${\rm e}^{pt}$ will suppress the integral. All that will remain from this is the contribution from the poles of the integral. Generally, as well as the ‘ballistic poles’ on the imaginary line, the dielectric function can have poles in the left and right halves of the complex plane. Poles of the dielectric function in the left half of the complex plane, however, correspond to decaying modes, which we will neglect (see figure 1). Poles of the dielectric function in the right half of the complex plane correspond to linear instabilities of the distribution function and in principle cannot be neglected. However, we may restrict ourselves to linearly stable distribution functions with the proviso that we take our initial condition to be the distribution function after all instabilities have vanished and that the initial evolution with instabilities growing and saturating must be treated by a different theory. Within these assumptions the result of the contour integration is

(2.22)\begin{gather} I_{\boldsymbol{k}}(\boldsymbol{v},\boldsymbol{v}',\boldsymbol{w}) = \frac{{\rm i}}{\epsilon_{\boldsymbol{k}}(-{\rm i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v})\epsilon^{*}_{\boldsymbol{k}}(-{\rm i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}')} \frac{{\rm e}^{-{\rm i}\boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{v}-\boldsymbol{v}')t}}{\boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{v}-\boldsymbol{w})}\left[1 - \frac{\epsilon_{\boldsymbol{k}}(-{\rm i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v})}{\epsilon_{\boldsymbol{k}}(-{\rm i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{w})} {\rm e}^{{\rm i}\boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{v}-\boldsymbol{w})t} \right]. \end{gather}

We have previously assumed a separation of timescales between the mean distribution function and the perturbed distribution function (known generally as the Bogoliubov ansatz: see Swanson Reference Swanson2008). Therefore, the perturbation can be allowed to evolve to large $t$ before we consider its effect on the mean distribution function. Accordingly, we take the limit $t\to \infty$ in (2.21):

(2.23)\begin{gather} I_{\boldsymbol{k}}(\boldsymbol{v},\boldsymbol{v}',\boldsymbol{w}) = {\rm \pi}\frac{{\rm e}^{-{\rm i}\boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{v}-\boldsymbol{v}')t}}{\epsilon_{\boldsymbol{k},\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}}\epsilon^{*}_{\boldsymbol{k},\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}'}}\delta\left(\boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{v}-\boldsymbol{w}) \right), \end{gather}

where $\epsilon _{\boldsymbol {k},\boldsymbol {k}\boldsymbol {\cdot }\boldsymbol {v}} \equiv \epsilon _{\boldsymbol {k}}(-{\rm i}\boldsymbol {k}\boldsymbol {\cdot }\boldsymbol {v})$. Substituting (2.23) back into (2.20) gives a simplified expression for the diffusion kernel:

(2.24)\begin{align} \boldsymbol{\mathsf{D}}^{\alpha}_{\mu\nu}(\boldsymbol{w},\boldsymbol{v}) & = \sum_{\nu'}\frac{16{\rm \pi}^{3}q_{\mu}^{2}q_{\nu}q_{\nu'}}{m_{\alpha}m_{\mu}}\mathrm{Re}\sum_{\boldsymbol{k}}\frac{\boldsymbol{k}\boldsymbol{k}}{k^{4}}\delta(\boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{w}-\boldsymbol{v}))\nonumber\\& \quad \times\int\,\mathrm{d}\boldsymbol{v}' \, \frac{{\rm e}^{-{\rm i}\boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{v}- \boldsymbol{v}')t}\langle g_{\boldsymbol{k}\nu}(\boldsymbol{v})g^{*}_{\boldsymbol{k}\nu'}(\boldsymbol{v}') \rangle}{\epsilon_{\boldsymbol{k},\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}}\epsilon_{\boldsymbol{k},\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}'}^{*}}. \end{align}

Note that it is in obtaining this expression that the quasilinear approximation has first been truly used. In the more general expression (2.19), $\hat {h}_{\boldsymbol {k}\nu }(p,\boldsymbol {v})$, instead of being given by (2.11), can always be assumed to contain the nonlinear contributions neglected in (2.5).

Figure 1. The contours used for the $p$ integration in (2.21). The right-hand contour is the original contour, which is calculated by deforming it around the poles of (2.22) onto the left-hand contour. Poles of the dielectric function would correspond to decaying modes but are ignored (shown in the dotted line).

2.4. Conservation laws of the general quasilinear collision integral

Due to the symmetries that are possessed by (2.24), it is now possible to show that the collision integral (2.18) with the diffusion kernel (2.24) conserves the particle number, total mean momentum, and total mean energy between the species.

The number of particles of species $\alpha$

(2.25)\begin{equation} N_{\alpha} = V\int\,\mathrm{d}\boldsymbol{v} \,f_{0\alpha}, \end{equation}

is trivially conserved because the right-hand side of (2.18) is a full derivative with respect to $\boldsymbol {v}$.

The total momentum of the system is given by

(2.26)\begin{equation} P = V\sum_{\alpha}\int\,\mathrm{d}\boldsymbol{v} \,m_{\alpha}\boldsymbol{v}f_{0\alpha}. \end{equation}

Taking the time derivative of (2.26) using (2.18), we get

(2.27)\begin{align} \frac{\mathrm{d}P}{\mathrm{d}t} & = V\sum_{\alpha\alpha''}\iint\,\mathrm{d}\boldsymbol{v} \,\,\mathrm{d}\boldsymbol{v}'' \,m_{\alpha}\boldsymbol{v} \frac{\partial }{\partial \boldsymbol{v}}\boldsymbol{\cdot}\left[\boldsymbol{\mathsf{D}}^{\alpha}_{\alpha\alpha''}(\boldsymbol{v},\boldsymbol{v}'')\boldsymbol{\cdot}\frac{\partial f_{0\alpha}}{\partial \boldsymbol{v}}- \boldsymbol{\mathsf{D}}^{\alpha}_{\alpha''\alpha}(\boldsymbol{v}'',\boldsymbol{v})\boldsymbol{\cdot}\frac{\partial f_{0\alpha''}}{\partial \boldsymbol{v}''} \right]\nonumber\\ & ={-}V\sum_{\alpha\alpha''}\iint\,\mathrm{d}\boldsymbol{v} \,\,\mathrm{d}\boldsymbol{v}'' \,\left[m_{\alpha}\boldsymbol{\mathsf{D}}^{\alpha}_{\alpha\alpha''}(\boldsymbol{v},\boldsymbol{v}'')- m_{\alpha''}\boldsymbol{\mathsf{D}}^{\alpha''}_{\alpha\alpha''}(\boldsymbol{v},\boldsymbol{v}'') \right]\boldsymbol{\cdot}\frac{\partial f_{0\alpha}}{\partial \boldsymbol{v}} \nonumber\\ & = 0. \end{align}

In the second equality, we integrated by parts and swapped the $\alpha$ and $\alpha ''$ indices, as well as the arguments $\boldsymbol {v}$ and $\boldsymbol {v}''$ in the second term of the integral. This expresses the condition for momentum conservation as the symmetry

(2.28)\begin{equation} m_{\alpha}\boldsymbol{\mathsf{D}}^{\alpha}_{\alpha\alpha''}(\boldsymbol{v},\boldsymbol{v}'')= m_{\alpha''}\boldsymbol{\mathsf{D}}^{\alpha''}_{\alpha\alpha''}(\boldsymbol{v},\boldsymbol{v}''). \end{equation}

From (2.24), and indeed already from (2.20), we see that this is manifestly satisfied, hence the total momentum is conserved but particles of different species may exchange momentum.

We apply much the same procedure to the total energy

(2.29)\begin{equation} E = V\sum_{\alpha}\int \,\mathrm{d}\boldsymbol{v} \,\frac{1}{2}m_{\alpha}|\boldsymbol{v}|^{2}f_{0\alpha}(\boldsymbol{v}). \end{equation}

Taking the time derivative of (2.29), we get

(2.30)\begin{align} \frac{\mathrm{d}E}{\mathrm{d}t} & = V\sum_{\alpha\alpha''}\iint\,\mathrm{d}\boldsymbol{v} \,\,\mathrm{d}\boldsymbol{v}'' \,\frac{1}{2}m_{\alpha}|\boldsymbol{v}|^{2}\frac{\partial }{\partial \boldsymbol{v}}\boldsymbol{\cdot}\left[\boldsymbol{\mathsf{D}}^{\alpha}_{\alpha\alpha''}(\boldsymbol{v},\boldsymbol{v}'')\boldsymbol{\cdot}\frac{\partial f_{0\alpha}}{\partial \boldsymbol{v}}- \boldsymbol{\mathsf{D}}^{\alpha}_{\alpha''\alpha}(\boldsymbol{v}'',\boldsymbol{v})\boldsymbol{\cdot}\frac{\partial f_{0\alpha''}}{\partial \boldsymbol{v}''} \right]\nonumber\\ & ={-}V\sum_{\alpha\alpha''}\iint\,\mathrm{d}\boldsymbol{v} \,\,\mathrm{d}\boldsymbol{v}'' \,\left[m_{\alpha}\boldsymbol{v}\boldsymbol{\cdot} \boldsymbol{\mathsf{D}}^{\alpha}_{\alpha\alpha''}(\boldsymbol{v},\boldsymbol{v}'')- m_{\alpha''}\boldsymbol{v}''\boldsymbol{\cdot} \boldsymbol{\mathsf{D}}^{\alpha''}_{\alpha\alpha''}(\boldsymbol{v},\boldsymbol{v}'') \right]\boldsymbol{\cdot}\frac{\partial f_{0\alpha}}{\partial \boldsymbol{v}} \nonumber\\ & = 0. \end{align}

Again, the swapping of the indices and velocities allows the condition of total-energy conservation to be cast as the symmetry

(2.31)\begin{equation} m_{\alpha}\boldsymbol{v}\boldsymbol{\cdot} \boldsymbol{\mathsf{D}}^{\alpha}_{\alpha\alpha''}(\boldsymbol{v},\boldsymbol{v}'') = m_{\alpha''}\boldsymbol{v}''\boldsymbol{\cdot} \boldsymbol{\mathsf{D}}^{\alpha''}_{\alpha\alpha''}(\boldsymbol{v},\boldsymbol{v}''), \end{equation}

which is satisfied by (2.24) due to the delta function that emerged from the approximation relating to the separation of timescales. Essentially, what this means is that, when the mean distribution function is considered to be linearly stable and to change slowly compared to the fluctuations, the energy $E$ of this distribution function cannot change.

While it is encouraging that we have a collision integral in a general form into which conservation laws are hard-wired, without a general form of the correlator $\langle g_{\alpha \boldsymbol {k}}(\boldsymbol {v})g^{*}_{\nu \boldsymbol {k}}(\boldsymbol {v}') \rangle$ it is still not possible to determine the evolution of $f_{0\alpha }$. We will now introduce the approximations necessary to arrive at collision integrals in a closed form.

2.5. Microgranulation ansatz

The collision integral (2.18), with $\boldsymbol {\mathsf {D}}$ given by (2.24), expresses the evolution of the mean distribution function for each species in terms of correlators of the form $\langle g_{\boldsymbol {k}\alpha }(\boldsymbol {v})g^{*}_{\boldsymbol {k}\alpha '}(\boldsymbol {v}') \rangle$. Despite $g_{\boldsymbol {k}\alpha }$ first entering the calculation as an initial condition, we must now interpret it in a slightly different way.

Suppose we begin at time $t_{0}$ with a fluctuation in the distribution function in the form of an initial condition $g_{\boldsymbol {k}\alpha } = f_{\boldsymbol {k}\alpha }(t = t_{0})$. We then evolve this perturbation, $f_{\boldsymbol {k}\alpha }(t)$, according to the linearised Vlasov equation (2.5). This evolution is fast, and we assume that the mean distribution, $f_{0\alpha }(t_0)$, remains constant during it. But of course, the mean distribution does evolve, albeit on a longer time scale. So if we pick some later time $t_{1}$, such that $t_{1}-t_{0}$ is comparable to that longer time scale, we will find a slightly altered mean distribution function, $f_{0\alpha }(t_{1})$. It would then make sense to evolve the perturbation $f_{\boldsymbol {k}\alpha }$ using this new mean distribution. We therefore ‘reset’ the initial condition to $g_{\boldsymbol {k}\alpha } = f_{\boldsymbol {k}\alpha }(t=t_{1})$ and restart the evolution of $f_{\boldsymbol {k}\alpha }(t)$ from this ‘new’ initial state.

This is very similar to how one might do this numerically: we evolve linearly under a given $f_{0\alpha }$ and initial condition; then, after a small amount of time we, ratchet-like, restart the system, declaring the new initial conditions of this next iteration to be the final state of the previous iteration. The ideal model for the correlation function of $g_{\boldsymbol {k}\alpha }$ would then be one that continually updated to be the correct one for a given mean distribution function. The determination of the steady-state phase-space correlation function associated with a given $f_{0\alpha }$ is a difficult problem, not in general solved (cf. Adkins & Schekochihin Reference Adkins and Schekochihin2018). Instead, following Kadomtsev & Pogutse (Reference Kadomtsev and Pogutse1970), Balescu (Reference Balescu1960) and Lenard (Reference Lenard1960), we will make some simplifying assumptions about the correlations of $g_{\boldsymbol {k}\alpha }$. The first step towards doing this is to assume that these correlations are very short ranged in phase space, i.e., that the correlator $\langle g_{\alpha }(\boldsymbol {r},\boldsymbol {v})g_{\alpha '}(\boldsymbol {r}',\boldsymbol {v}') \rangle$ is zero unless the points $(\boldsymbol {r},\boldsymbol {v})$ and $(\boldsymbol {r}',\boldsymbol {v}')$ lie very close to each other. We also assume that the perturbed distribution functions of different species are uncorrelated. Mathematically, we express these assumptions in the form of the microgranulation ansatz

(2.32)\begin{equation} \langle g_{\nu}(\boldsymbol{r},\boldsymbol{v})g_{\nu'}(\boldsymbol{r}',\boldsymbol{v}') \rangle = \Delta\Gamma_{\nu}\delta_{\nu\nu'}\delta(\boldsymbol{r}-\boldsymbol{r}')\delta(\boldsymbol{v}-\boldsymbol{v}')\langle g_{\nu}^{2} \rangle(\boldsymbol{v}). \end{equation}

Here, the remaining correlator $\langle g_{\nu }^{2} \rangle (\boldsymbol {v})$ is assumed to be spatially independent due to the statistical homogeneity of the system. The new parameter $\Delta \Gamma _{\nu }$ is the ‘volume’ of phase space over which the distribution function of a given species has correlated fluctuations. The ansatz (2.32) does not yet constitute a closure as we have not specified the correlator $\langle g_{\nu }^{2} \rangle (\boldsymbol {v})$ in terms of $f_{0\nu }(\boldsymbol {v})$. But before this is done, let us implement the ansatz (2.32) to simplify the diffusion kernel (2.24) again.

The Fourier-space correlation function that appears in (2.24) is easily computed from (2.32):

(2.33)\begin{align} \langle g_{\boldsymbol{k}\nu}(\boldsymbol{v})g_{\boldsymbol{k}\nu'}^{*}(\boldsymbol{v}') \rangle & = \iint\frac{\mathrm{d}\boldsymbol{r} \,\,\mathrm{d}\boldsymbol{r}' \,}{V^{2}} {\rm e}^{-{\rm i}\boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{r}-\boldsymbol{r}')}\langle g_{\nu}(\boldsymbol{r},\boldsymbol{v})g_{\nu'}(\boldsymbol{r}',\boldsymbol{v}') \rangle\nonumber\\ & = \frac{\Delta\Gamma_{\nu}}{V} \delta_{\nu \nu'}\langle g_{\nu}^{2} \rangle(\boldsymbol{v})\delta(\boldsymbol{v}-\boldsymbol{v}'). \end{align}

Thus, the microgranulation ansatz simplifies (2.24) to

(2.34)\begin{equation} \boldsymbol{\mathsf{D}}^{\alpha}_{\mu \nu}(\boldsymbol{w},\boldsymbol{v}) = \frac{16{\rm \pi}^{3}q_{\mu}^{2}q_{\nu}^{2}\Delta\Gamma_{\nu}}{m_{\alpha}m_{\mu}}\sum_{\boldsymbol{k}}\frac{\boldsymbol{k}\boldsymbol{k}}{k^{4}}\frac{\delta(\boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{w}-\boldsymbol{v}))}{|\epsilon_{\boldsymbol{k},\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}}|^{2}}\langle g_{\nu}^{2} \rangle(\boldsymbol{v}). \end{equation}

This turns the collision integral (2.18) into the following form, written in terms of an as yet undetermined correlator $\langle g_{\alpha }^{2} \rangle (\boldsymbol {v})$:

(2.35)\begin{gather} \frac{\partial f_{0\alpha}}{\partial t} = \sum_{\alpha''}\frac{16{\rm \pi}^{3}q_{\alpha}^{2}q_{\alpha''}^{2}}{m_{\alpha}V}\frac{\partial }{\partial \boldsymbol{v}}\boldsymbol{\cdot}\sum_{\boldsymbol{k}}\frac{\boldsymbol{k}\boldsymbol{k}}{k^{4}}\boldsymbol{\cdot}\int \,\mathrm{d}\boldsymbol{v}'' \,\frac{\delta(\boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{v}-\boldsymbol{v}''))}{\left|\epsilon_{\boldsymbol{k},\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}} \right|^{2}} \nonumber\\ \left[\frac{\Delta\Gamma_{\alpha''}}{m_{\alpha}}\langle g_{\alpha''}^{2} \rangle(\boldsymbol{v}'')\frac{\partial f_{0\alpha}}{\partial \boldsymbol{v}} - \frac{\Delta\Gamma_{\alpha}}{m_{\alpha''}}\langle g_{\alpha}^{2} \rangle(\boldsymbol{v})\frac{\partial f_{0\alpha''}}{\partial \boldsymbol{v}''}\right]. \end{gather}

In the next section we will show that one can link $\langle g_{\alpha }^{2} \rangle (\boldsymbol {v})$ to assumptions about the nature of the exact distribution function, leading finally to a closure in terms of $f_{0\alpha }$.

However, first, let us take care of a piece of unfinished business: we are now in a position to confirm that the first term in (2.16) does indeed vanish. The contribution of that term to (2.14) is

(2.36)\begin{equation} -\frac{\partial }{\partial \boldsymbol{v}}\boldsymbol{\cdot}\sum_{\boldsymbol{k}}\boldsymbol{k}\mathrm{Im}\sum_{\alpha'}\frac{4{\rm \pi} q_{\alpha}q_{\alpha'}}{m_{\alpha}k^{2}} \iint\frac{\mathrm{d}p \,\,\mathrm{d}p' \,}{(2{\rm \pi})^{2}} \frac{{\rm e}^{(p+p')t}}{\epsilon_{\boldsymbol{k}}(p)\epsilon_{\boldsymbol{k}}^{*}(p'^{*})} \int\,\mathrm{d}\boldsymbol{v}' \,\langle \hat{h}_{\boldsymbol{k}\alpha}(p,\boldsymbol{v})\hat{h}_{\boldsymbol{k}\alpha'}^{*}(p'^{*},\boldsymbol{v}') \rangle. \end{equation}

Therefore, this term will be zero if the quantity

(2.37)\begin{equation} \sum_{\alpha'}\frac{4{\rm \pi} q_{\alpha}q_{\alpha'}}{m_{\alpha}k^{2}} \iint\frac{\mathrm{d}p \,\,\mathrm{d}p' \,}{(2{\rm \pi})^{2}} \frac{{\rm e}^{(p+p')t}}{\epsilon_{\boldsymbol{k}}(p)\epsilon_{\boldsymbol{k}}^{*}(p'^{*})} \int\,\mathrm{d}\boldsymbol{v}' \,\frac{\langle g_{\boldsymbol{k}\alpha}(\boldsymbol{v})g^{*}_{\boldsymbol{k}\alpha'}(\boldsymbol{v}') \rangle}{\left(p+{\rm i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}\right)\left(p'-{\rm i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}'\right)}, \end{equation}

is real. Taking the complex conjugate of (2.37) and permuting $p\leftrightarrow p'$, we find

(2.38)\begin{equation} \sum_{\alpha'}\frac{4{\rm \pi} q_{\alpha}q_{\alpha'}}{m_{\alpha}k^{2}}\iint\frac{\mathrm{d}p \,\,\mathrm{d}p' \,}{(2{\rm \pi})^{2}} \frac{{\rm e}^{(p+p')t}}{\epsilon_{\boldsymbol{k}}(p)\epsilon^{*}_{\boldsymbol{k}}(p'^{*})}\int\,\mathrm{d}\boldsymbol{v}' \,\frac{\langle g^{*}_{\boldsymbol{k}\alpha}(\boldsymbol{v})g_{\boldsymbol{k}\alpha'}(\boldsymbol{v}') \rangle}{\left(p' - {\rm i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}\right)\left(p + {\rm i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}'\right)}, \end{equation}

which is the same as (2.37) if the correlation function satisfies

(2.39)\begin{equation} \int\,\mathrm{d}\boldsymbol{v}' \,\frac{\langle g_{\boldsymbol{k}\alpha}(\boldsymbol{v})g^{*}_{\boldsymbol{k}\alpha'}(\boldsymbol{v}') \rangle}{\left(p+{\rm i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}\right)\left(p'-{\rm i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}'\right)} =\int\,\mathrm{d}\boldsymbol{v}' \,\frac{\langle g^{*}_{\boldsymbol{k}\alpha}(\boldsymbol{v})g_{\boldsymbol{k}\alpha'}(\boldsymbol{v}') \rangle}{\left(p' - {\rm i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}\right)\left(p + {\rm i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}'\right)}. \end{equation}

The microgranulation ansatz (2.32) manifestly does satisfy this, so our earlier neglect of the first term in (2.16) is vindicated. Of course, the microgranulation ansatz is quite a simplification of the correlation function of the phase-space density. It is possible that the true correlation function would not have this symmetry and, consequently, give rise to a qualitatively different evolution of the mean distribution function. Such a possibility will be explored in a separate publication.

3.1. Single-waterbag closure

When only one waterbag density $\eta _{\alpha }$ for each species $\alpha$ is assumed, this has immediate implications for the mean phase-space density. If the exact phase-space density $f_{\alpha }(\boldsymbol {r},\boldsymbol {v})$ is only ever $\eta _{\alpha }$ or zero, then its average value is directly related to the probability of a portion of phase space being occupied. We will therefore define $p_{\alpha }(\boldsymbol {v})$ as the probability that, at a given point $\boldsymbol {v}$, the phase-space density is $\eta _{\alpha }$ for particles of species $\alpha$. This can also be said as ‘$p_{\alpha }(\boldsymbol {v})$ is the probability that there is a waterbag of species $\alpha$ at $\boldsymbol {v}$’ with the understanding that a waterbag is a patch of phase space of a certain density.

Written in terms of this $p_{\alpha }(\boldsymbol {v})$, the mean phase-space density $f_{0\alpha }(\boldsymbol {v})$ is then

(3.2)\begin{equation} f_{0\alpha}(\boldsymbol{v}) = \langle f_{\alpha}(\boldsymbol{r},\boldsymbol{v}) \rangle =\eta_{\alpha} p_{\alpha}(\boldsymbol{v}). \end{equation}

Note that there is no $\boldsymbol {r}$ dependence because our system is assumed to be statistically homogeneous in the position space. Likewise, the squared phase-space density will be $\eta _{\alpha }^{2}$ with probability $p_{\alpha }$ or zero otherwise. Therefore,

(3.3)\begin{equation} \langle f_{\alpha}^{2} \rangle(\boldsymbol{v}) = \eta_{\alpha}^{2}p_{\alpha}(\boldsymbol{v}) = \eta_{\alpha} f_{0\alpha}(\boldsymbol{v}). \end{equation}

The correlator $\langle g_{\alpha }^{2} \rangle (\boldsymbol {v})$ can then be determined immediately:

(3.4)\begin{equation} \langle g_{\alpha}^{2} \rangle(\boldsymbol{v}) = \left\langle \left[f_{\alpha}(\boldsymbol{r},\boldsymbol{v}) - f_{0\alpha}(\boldsymbol{v})\right]^{2} \right\rangle = \langle f_{\alpha}^{2} \rangle(\boldsymbol{v})-f_{0\alpha}^{2}(\boldsymbol{v}) = \left[\eta_\alpha - f_{0\alpha}(\boldsymbol{v}) \right]f_{0\alpha}(\boldsymbol{v}). \end{equation}

With this closure, (2.35) becomes the multi-species generalisation of Kadomtsev and Pogutse's collision integral:

(3.5)\begin{align} \frac{\partial f_{0\alpha}}{\partial t} & = \sum_{\alpha''}\frac{16{\rm \pi}^{3}q_{\alpha}^{2}q_{\alpha''}^{2}}{m_{\alpha}V}\frac{\partial }{\partial \boldsymbol{v}}\boldsymbol{\cdot}\sum_{\boldsymbol{k}}\frac{\boldsymbol{k}\boldsymbol{k}}{k^{4}}\boldsymbol{\cdot}\int\,\mathrm{d}\boldsymbol{v}'' \,\frac{\delta(\boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{v}-\boldsymbol{v}''))}{\left|\epsilon_{\boldsymbol{k},\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}} \right|^{2}}\nonumber\\ & \qquad \left\{\frac{\Delta\Gamma_{\alpha''}}{m_{\alpha}}\left[\eta_{\alpha''}-f_{0\alpha''}(\boldsymbol{v}'')\right]f_{0\alpha''}(\boldsymbol{v}'')\frac{\partial f_{0\alpha}}{\partial \boldsymbol{v}} - \frac{\Delta\Gamma_{\alpha}}{m_{\alpha''}}\left[\eta_{\alpha} - f_{0\alpha}(\boldsymbol{v})\right]f_{0\alpha}(\boldsymbol{v})\frac{\partial f_{0\alpha''}}{\partial \boldsymbol{v}''}\right\}. \end{align}

The fixed points of this collision integral are Fermi–Dirac distributions (see figure 2):

(3.6)\begin{equation} f_{0\alpha}(\boldsymbol{v}) = \frac{\eta_{\alpha}}{1 + \exp\left[\beta\Delta\Gamma_{\alpha}\eta_{\alpha}\left(\frac{1}{2}m_{\alpha}|\boldsymbol{v}|^{2} - \mu_{\alpha}\right) \right]}. \end{equation}

We can see that these are indeed fixed points by noting that, for (3.6),

(3.7)\begin{equation} \frac{\Delta\Gamma_{\alpha''}}{m_{\alpha}}\frac{\partial f_{0\alpha}}{\partial \boldsymbol{v}} ={-}\beta\Delta\Gamma_{\alpha}\Delta\Gamma_{\alpha''}\boldsymbol{v}\left[\eta_{\alpha} - f_{0\alpha}(\boldsymbol{v}) \right]f_{0\alpha}(\boldsymbol{v}). \end{equation}

Thus, the two bracketed terms in (3.5) become identical except for a factor of $\boldsymbol {v}$ or $\boldsymbol {v}''$. This difference vanishes when the bracket is dotted with $\delta (\boldsymbol {k}\boldsymbol {\cdot }(\boldsymbol {v}-\boldsymbol {v}''))\boldsymbol {k}$, setting the collision integral to zero. The demonstration that the fixed points (3.6) are stable will come in § 3.3, where we prove the general H-theorem for the multi-waterbag model, which reduces to the single-waterbag model in the appropriate limit. This H-theorem will guarantee that the mean phase-space density tends towards the maximum of a certain entropy (exactly the entropy proposed by Lynden-Bell Reference Lynden-Bell1967) subject to the conservation laws proved in § 2.4. These conservation laws will therefore enforce a particular choice of the parameters $\beta$ and $\mu _{\alpha }$ in 3.6, which we call the ‘thermodynamic beta’ and ‘chemical potentials’, respectively. This is the correct number of free parameters, one for each conservation law. Note that we have, without loss of generality, moved into the frame where the plasma has zero net momentum.

Figure 2. A one-dimensional cartoon of the exact and mean distribution functions (phase-space densities) for the single-waterbag model. As discussed in § 2.5, the mean distribution should be considered in the ensemble-averaged sense, so that the many realisations of systems like the exact one given in panel (a) evolve on average to look like the Fermi–Dirac distribution in panel (b).

Physically, the Fermi–Dirac distribution has emerged because it is the maximiser of an entropy, subject to the conservation of energy, particle number, and phase volume. Phase-volume conservation manifests as a Pauli-like exclusion effect: the waterbags cannot overlap, so they are forced to different points of phase space, removing the possibility of a Maxwellian equilibrium. Instead, the qualitative shape of the distribution function is set by the relation between the constant volume occupied by the waterbags and the constant energy of those waterbags. Clearly, for an initial condition with waterbags of species $\alpha$ occupying a phase volume $\varGamma _{\alpha }$ in phase space, there will be some finite minimum energy that the initial condition must possess. This minimum energy is non-zero because the possibility of all waterbags sinking to $\boldsymbol {v} = 0$ is forbidden by the exclusion principle. Instead, some waterbags will be forced to higher energies, giving a constant mean phase-space density for all energies below a certain value, the Fermi energy $\epsilon _{\mathrm {F}\alpha }$. The Fermi energy is then set by the initial volume of the waterbags with

(3.8)\begin{equation} \frac{\varGamma_{\alpha}}{V} = \frac{4{\rm \pi}}{3}\left(\frac{2\epsilon_{\mathrm{F}\alpha}}{m_{\alpha}} \right)^{3/2}, \end{equation}

where $V$ is the volume of the position space. The minimum possible energy $E_{\mathrm {min}}$ for such an initial condition is

(3.9)\begin{equation} \frac{E_{\mathrm{min}}}{V} = \sum_{\alpha}\frac{2{\rm \pi} m_{\alpha} \eta_{\alpha}}{5}\left(\frac{2\epsilon_{\mathrm{F}\alpha}}{m_{\alpha}} \right)^{5/2} = \sum_{\alpha}\frac{2{\rm \pi} m_{\alpha} \eta_{\alpha}}{5} \left(\frac{3\varGamma_{\alpha}}{4{\rm \pi} V} \right)^{5/3}. \end{equation}

When the energy of the system far exceeds this minimum, ${E \gg E_{\mathrm {min}}}$, the volume of phase space available to a given waterbag will be much greater than the volume occupied by other waterbags. In this limit, the exclusion effect will be unimportant and the distribution functions will be approximately Maxwellian with the parameters $\mu _{\alpha }$ and $\beta$ in (3.6) chosen to give the correct particle number for each species and the correct total energy. In the opposite limit, where the energy of the system begins very close to the minimum possible energy, $E-E_{\mathrm {min}} \ll E_{\mathrm {min}}$, the exclusion effect will be paramount, with waterbags tightly packed at low energies. This gives rise to a mean phase-space density that is a smoothed step function, taking a value nearly $\eta _{\alpha }$ at energies below $\epsilon _{\mathrm {F}\alpha }$ and nearly zero above it. The values of $\beta$ and $\mu _{\alpha }$ can then be determined via a Sommerfeld expansion in the same way as it is done in standard statistical mechanics.

3.2. Multi-waterbag closure

By assuming, as we did in § 3.1, that the exact phase-space density was a single-waterbag distribution, we made the closure for $\langle g_{\alpha }^{2} \rangle (\boldsymbol {v})$ very simple. The reason for this boils down to the fact that, for a single-waterbag distribution, knowledge of the mean phase-space density was sufficient to determine all the requisite statistical information about the exact one, namely the probabilities $p_{\alpha }(\boldsymbol {v})$. In a similar vein, one can ask what information there is to be determined in an $N$-waterbag system, where the exact phase-space density can take the values $\{ \eta _{J\alpha } \}_{J =1,2,\ldots,N}$. As a direct generalisation of the single-waterbag model, we then define $p_{J \alpha }(\boldsymbol {v})$ to be the probability that the phase-space density of species $\alpha$ will be equal to $\eta _{J\alpha }$ at velocity $\boldsymbol {v}$. In other words, ‘the $\eta _{J\alpha }$ waterbag of species $\alpha$ has probability $p_{J\alpha }(\boldsymbol {v})$ of being present at $\boldsymbol {v}$’. Then we can write the mean phase-space density as

(3.10)\begin{align} f_{0\alpha}(\boldsymbol{v}) = \langle f_{\alpha}(\boldsymbol{r},\boldsymbol{v}) \rangle = \sum_{J}\eta_{J\alpha}p_{J\alpha}(\boldsymbol{v}), \end{align}

and the mean square one as

(3.11)\begin{equation} \langle f_{\alpha}^{2} \rangle(\boldsymbol{v}) = \sum_{J}\eta_{J\alpha}^{2}p_{J\alpha}(\boldsymbol{v}). \end{equation}

Now the complication introduced by multiple waterbags becomes obvious. When ${N>1}$, knowledge of just $f_{0\alpha }$ cannot uniquely determine the probabilities $p_{J\alpha }(\boldsymbol {v})$ and so does not exactly determine $\langle g_{\alpha }^{2} \rangle$ in (2.35). Therefore, it is not possible to close (2.35) without further information. In principle, this information could be extracted. One could construct the evolution equations for the higher-order moments of the exact phase-space density $\langle f_{\alpha }^{2} \rangle$, $\langle f_{\alpha }^{3} \rangle$, …, $\langle f_{\alpha }^{N} \rangle$. The evolution of the $i$-th such moment would generally depend on the $(i+1)$-st moment. In that way, one would have $N$ analogues of (3.10) and (3.11) for the $N$ unknowns $p_{J\alpha }(\boldsymbol {v})$. This would then allow one to write $\langle f_{\alpha }^{N+1} \rangle$ as a function of all the lower-order moments and finally close the system, featuring $N$ coupled collision integrals.

In § 4, we will show that such a scheme can be made tractable by calculating all such moments in one fell swoop by increasing the dimension of the phase space. First, however, we consider a simple closure that will prove illuminating in understanding the relaxation of collisionless systems. Instead of hoping to gain any further knowledge of the system, we ask what is the most likely assignment of probabilities $p_{J\alpha }(\boldsymbol {v})$ given that the mean phase-space density is $f_{0\alpha }(\boldsymbol {v})$ and nothing else is known. This can be answered in the spirit of statistical inference by maximising the standard Shannon (Reference Shannon1948) entropy:

(3.12)\begin{equation} S_{\alpha} ={-}\int\,\mathrm{d}\boldsymbol{v} \,\sum_{J = 0}^{N}p_{J\alpha}(\boldsymbol{v})\ln p_{J\alpha}(\boldsymbol{v}), \end{equation}

where the sum now includes $J=0$ as ‘the empty waterbag’, representing the probability that a given point in phase space is has zero density of particles. The Shannon entropy (3.12) must be maximised subject to the condition that the mean phase-space density is $f_{0\alpha }(\boldsymbol {v})$, i.e., that the $p_{J\alpha }(\boldsymbol {v})$'s obey (3.10), and that the probabilities $p_{J\alpha }(\boldsymbol {v})$ sum to unity at each $\boldsymbol {v}$. Another constraint on these probabilities is that the total number of particles per unit volume contained within each waterbag,

(3.13)\begin{equation} n_{J\alpha} = \eta_{J\alpha}\int\,\mathrm{d}\boldsymbol{v} \,p_{J\alpha}(\boldsymbol{v}), \end{equation}

is fixed for all $J \neq 0$, since the phase volume corresponding to each $\eta _{J\alpha }$ must be conserved by the evolution under the Vlasov equation. Such invariants are often called (or are equivalent to) Casimir invariants (cf. Chavanis Reference Chavanis2004, Reference Chavanis2005; Zhdankin Reference Zhdankin2021).

Thus, we maximise the following functional for each species:

(3.14)\begin{align} -\int& \,\mathrm{d}\boldsymbol{v} \,\sum_{J}p_{J\alpha}(\boldsymbol{v})\ln p_{J\alpha}(\boldsymbol{v}) - \int\,\mathrm{d}\boldsymbol{v} \,\psi_{\alpha}(\boldsymbol{v})\left[\sum_{J}\eta_{J\alpha}p_{J\alpha}(\boldsymbol{v}) - f_{0\alpha}(\boldsymbol{v}) \right] \nonumber\\ & - \sum_{J\neq 0}\gamma_{J\alpha}\left[\int \,\mathrm{d}\boldsymbol{v} \,p_{J\alpha}(\boldsymbol{v}) - \frac{n_{J\alpha}}{\eta_{J\alpha}}\right] - \int\,\mathrm{d}\boldsymbol{v} \,\lambda_{\alpha}(\boldsymbol{v}) \left[\sum_{J}p_{J\alpha}(\boldsymbol{v}) - 1 \right] \to \max, \end{align}

where $\psi _{\alpha }(\boldsymbol {v})$, $\gamma _{J\alpha }$, and $\lambda _{\alpha }(\boldsymbol {v})$ are Lagrange multipliers. The result is

(3.15)\begin{gather} p_{J\alpha}(\boldsymbol{v}) = \frac{1}{Z_{\alpha}\left(\psi_{\alpha}(\boldsymbol{v})\right)} {\rm e}^{-\psi_{\alpha}(\boldsymbol{v})\eta_{J\alpha}- \gamma_{J\alpha}}, \end{gather}

where the ‘partition function’ of species $\alpha$ is

(3.16)\begin{gather} Z_{\alpha}(\psi_{\alpha}(\boldsymbol{v})) = 1 + \sum_{J\neq 0} {\rm e}^{-\psi_{\alpha}(\boldsymbol{v})\eta_{J\alpha} - \gamma_{J\alpha}}, \end{gather}

and the Lagrange multipliers $\psi _{\alpha }(\boldsymbol {v})$ and $\gamma _{J\alpha }$ must be chosen to enforce the constraints (3.10) and (3.13). Analogously to the standard Gibbs (Reference Gibbs1902) statistical mechanics, (3.10) becomes

(3.17)\begin{gather} f_{0\alpha}(\boldsymbol{v}) = \frac{1}{Z_{\alpha}\left(\psi_{\alpha}(\boldsymbol{v})\right)}\sum_{J}\eta_{J\alpha} {\rm e}^{-\psi_{\alpha}(\boldsymbol{v})\eta_{J\alpha} - \gamma_{J\alpha}} ={-}\frac{\partial \ln Z_{\alpha}}{\partial \psi_{\alpha}}. \end{gather}

Thus, the mean phase-space density plays the role that energy does in the regular prescription of statistical mechanics, and $\psi _{\alpha }(\boldsymbol {v})$ that of inverse temperature, which is local in $\boldsymbol {v}$. In this formalism, therefore,

(3.18)\begin{align} \langle g_{\alpha}^{2} \rangle(\boldsymbol{v}) & = \frac{1}{Z_{\alpha}(\psi_{\alpha}(\boldsymbol{v}))}\sum_{J}\eta_{J\alpha}^{2} {\rm e}^{-\psi_{\alpha}(\boldsymbol{v})\eta_{J\alpha} - \gamma_{J\alpha}} - f_{0\alpha}^{2}(\boldsymbol{v}) \nonumber\\ & = \frac{1}{Z_{\alpha}}\frac{\partial^{2} Z_{\alpha}}{\partial \psi_{\alpha}^{2}} - \frac{1}{Z_{\alpha}^{2}}\left(\frac{\partial Z_{\alpha}}{\partial \psi_{\alpha}}\right)^{2} = \frac{\partial^{2} \ln Z_{\alpha}}{\partial \psi_{\alpha}^{2}} ={-}\frac{\partial f_{0\alpha}}{\partial \psi_{\alpha}}, \end{align}

reminiscent of the heat capacity of a system in the Gibbs ensemble. Such a closure appears to have been first proposed by Chavanis (Reference Chavanis2005), in the context of geophysical turbulence.Footnote ¹

With $\langle g_{\alpha }^{2} \rangle (\boldsymbol {v})$ thus specified, we may substitute (3.18) into (2.35) to get

(3.19)\begin{gather} \frac{\partial f_{0\alpha}}{\partial t} ={-}\sum_{\alpha''}\frac{16{\rm \pi}^{3}q_{\alpha}^{2}q_{\alpha''}^{2}}{m_{\alpha}V}\frac{\partial }{\partial \boldsymbol{v}}\boldsymbol{\cdot}\sum_{\boldsymbol{k}}\frac{\boldsymbol{k}\boldsymbol{k}}{k^{4}}\boldsymbol{\cdot}\int \,\mathrm{d}\boldsymbol{v}'' \,\frac{\delta(\boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{v}-\boldsymbol{v}''))}{\left|\epsilon_{\boldsymbol{k},\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}} \right|^{2}}\nonumber\\ \left(\frac{\Delta\Gamma_{\alpha''}}{m_{\alpha}}\frac{\partial \psi_{\alpha}}{\partial \boldsymbol{v}} - \frac{\Delta\Gamma_{\alpha}}{m_{\alpha''}}\frac{\partial \psi_{\alpha''}}{\partial \boldsymbol{v}''} \right)\frac{\partial f_{0\alpha}}{\partial \psi_{\alpha}}(\boldsymbol{v})\frac{\partial f_{0\alpha''}}{\partial \psi_{\alpha''}}(\boldsymbol{v}''). \end{gather}

This is the collision integral for a multi-waterbag Lynden-Bell plasma. The instantaneous relationship between $f_{0\alpha }(\boldsymbol {v})$ and $\psi _{\alpha }(\boldsymbol {v})$ is given by (3.17) with $Z_{\alpha }$ defined by (3.16) and $\gamma _{J\alpha }$'s set by (3.13) and (3.15). The set of constants $n_{J\alpha }$ in (3.13) is fixed by the initial condition and cannot change during the evolution of $f_{0\alpha }$. This is the way in which phase-volume conservation in a collisionless plasma imprints a signature of the initial distribution (its ‘waterbag content’) on its otherwise universal evolution towards the Lynden-Bell equilibria.

Note that the closure proposed above amounts to assuming that the system always quickly attains a ‘local’ equilibrium in phase space given by (3.15) with a $\boldsymbol {v}$-dependent ‘inverse phase temperature’ $\psi _{\alpha }(\boldsymbol {v})$. The integral (3.19) then describes the evolution toward a ‘global’ equilibrium. We shall discuss the plausibility of this assumption in § 4.3.

3.3. Properties of the multi-waterbag collision integral

Having derived the collision integral (3.19), we now proceed to study its properties. First, we will determine its fixed points, then confirm that they are stable attractors by proving an H-theorem for our collision integral.

3.3.1. Multi-waterbag equilibria

Since $\partial f_{0\alpha } / \partial \psi _{\alpha } \leq 0$ (and only zero in pathological cases), the integral on the right-hand side of (3.19) can vanish only if its integrand vanishes. Manifestly, it does so if

(3.20)\begin{equation} \psi_{\alpha}(\boldsymbol{v}) = \beta\Delta\Gamma_{\alpha}\frac{1}{2}m_{\alpha}|\boldsymbol{v}|^{2} \equiv \beta \Delta\Gamma_{\alpha}\epsilon_{\alpha}(\boldsymbol{v}), \end{equation}

where $\epsilon _{\alpha }(\boldsymbol {v})$ is the energy of a particle of species $\alpha$ with velocity $\boldsymbol {v}$. Therefore, from (3.15),

(3.21)\begin{equation} p_{J\alpha}(\boldsymbol{v}) = \frac{\exp\left\{-\beta \Delta\Gamma_{\alpha}\eta_{J\alpha}\left[\epsilon_{\alpha}(\boldsymbol{v}) - \mu_{J\alpha}\right]\right\}}{1+ \sum_{J'\neq 0}\exp\left\{-\beta \Delta\Gamma_{\alpha}\eta_{J'\alpha}\left[\epsilon_{\alpha}(\boldsymbol{v}) - \mu_{J'\alpha}\right]\right\}}. \end{equation}

We have brought this into a form pleasingly similar to the Fermi–Dirac distribution by rewriting $\gamma _{J\alpha } = -\beta \Delta \Gamma _{\alpha }\eta _{J\alpha }\mu _{J\alpha }$ to define the chemical potential $\mu _{J\alpha }$ of the waterbag $J$ of species $\alpha$. These are the equilibria derived by Lynden-Bell (Reference Lynden-Bell1967) by a statistical-mechanical entropy-maximisation method applied to a multi-waterbag system. In our derivation, they have emerged as fixed points of a plasma's dynamical evolution. In § 3.3.2, we will prove that these fixed points of our collision integral (3.19) are stable and that the system relaxes towards them. First, however, let us discuss the qualitative nature of these solutions.

The similarity between the multi-waterbag equilibria (3.21) and the Fermi–Dirac distribution is no accident. Just like the Fermi–Dirac distribution, the Lynden-Bell equilibria maximise an entropy subject to the conservation of energy, momentum, and phase volume. Conservation of phase volume in the single-waterbag case meant waterbags excluded each other. In the multi-waterbag case, the same is true and it must be stressed that waterbags of different phase-space densities exclude each other indiscriminately. The result is a very broad class of possible distributions, whose general features, however, can be understood from the relation between the total energy of the system and the phase volume occupied by each waterbag.

Just like in the case of the Fermi–Dirac distribution, for an initial condition with waterbags of phase-space densities $\eta _{J\alpha }$ of species $\alpha$ each of which occupies a phase volume $\varGamma _{J\alpha }=Vn_{J\alpha }/\eta _{J\alpha }$, there is a minimum possible energy of the system, $E_{\mathrm {min}}$. A state with this minimum energy is the state that has the waterbag of the highest phase-space density at the lowest energy with waterbags of progressively lower phase-space density forced to higher energies by the exclusion of phase volume. This is analogous to how a suspension of liquids of different densities will arrange themselves in a glass, with the densest taking the lowest energy (cf. Lorenz Reference Lorenz1955; Gardner Reference Gardner1963). Let us order $\eta _{1\alpha }<\eta _{2\alpha }<\cdots <\eta _{N\alpha }$ by increasing waterbag density. The densest waterbag $\eta _{N\alpha }$ will occupy the sphere in velocity space with maximum energy $\epsilon _{\mathrm {F}J\alpha }$ given by

(3.22)\begin{equation} \frac{\varGamma_{N\alpha}}{V} = \frac{4{\rm \pi}}{3}\left(\frac{2\epsilon_{\mathrm{F}N\alpha}}{m_{\alpha}} \right)^{3/2}, \end{equation}

while the subsequent waterbags will be forced to higher energies, filling shells in velocity space. The waterbag of density $\eta _{J\alpha }$ will then occupy the shell between energies $\epsilon _{\mathrm {F}(J+1)\alpha }$ and $\epsilon _{\mathrm {F}J\alpha }$ given by

(3.23)\begin{equation} \frac{\varGamma_{J\alpha}}{V} = \frac{4{\rm \pi}}{3}\left[\left(\frac{2\epsilon_{\mathrm{F}J\alpha}}{m_{\alpha}} \right)^{3/2} - \left(\frac{2\epsilon_{\mathrm{F}(J+1)\alpha}}{m_{\alpha}} \right)^{3/2} \right], \end{equation}

or, equivalently,

(3.24)\begin{equation} \frac{4{\rm \pi}}{3}\left(\frac{2\epsilon_{\mathrm{F}J\alpha}}{m_{\alpha}} \right)^{3/2} = \sum_{n= J}^{N}\frac{\varGamma_{n\alpha}}{V}. \end{equation}

The minimum energy is then

(3.25)\begin{align} E_{\mathrm{min}} & = V\sum_{\alpha}\sum_{J=1}^{N}\frac{2{\rm \pi} m_{\alpha} \eta_{J\alpha}}{5}\left[\left(\frac{2\epsilon_{\mathrm{F}J\alpha}}{m_{\alpha}} \right)^{5/2} - \left(\frac{2\epsilon_{\mathrm{F}(J+1)\alpha}}{m_{\alpha}}\right)^{5/2} \right] \nonumber\\ & = V\sum_{\alpha}\sum_{J=1}^{N}\frac{2{\rm \pi} m_{\alpha} \eta_{J\alpha}}{5}\left(\frac{3}{4{\rm \pi} V}\right)^{5/3}\left[\left(\sum_{n=J}^{N}\varGamma_{n\alpha}\right)^{5/3} - \left(\sum_{n=J+1}^{N}\varGamma_{n\alpha}\right)^{5/3} \right]. \end{align}

At precisely this energy, the distribution as a function of velocity will look like a ziggurat, with sharp steps corresponding to the different waterbag densities. This state is known generically as a Gardner (Reference Gardner1963), or ‘Gardner-restacked’ distribution (normally thought of in the continuous limit of infinitely many infinitesimally separated waterbags – see § 3.4). It is a generic ‘ground state’ for a Vlasov plasma (cf. Helander Reference Helander2017).

However, this state is not at all representative of the full range of possible Lynden-Bell equilibria. Since the energy of the system remains constant under the collision integrals that we are considering (see § 2.4), for the Gardner distribution to reachable, the system must have started with the minimum energy possible, but that necessarily implies that it already was in the Gardner state. As we did in § 3.1 for the Fermi–Dirac distribution, let us ask for what range of initial energies the Lynden-Bell distribution (3.21) will be qualitatively similar to the Gardner distribution. To answer this question, we must consider the energy required to allow two waterbags to intermingle freely in velocity space as opposed to being cleanly separated. Should the $(J+1)$-st and $J$-th waterbags freely mix, they would have an average phase-space density

(3.26)\begin{equation} \bar{\eta} = \frac{\eta_{(J+1)\alpha} \varGamma_{(J+1)\alpha}+ \eta_{J\alpha} \varGamma_{J\alpha}}{\varGamma_{(J+1)\alpha} + \varGamma_{J\alpha}} \end{equation}

between the energies $\epsilon _{(J+2)\alpha }$ and $\epsilon _{J \alpha }$. The difference between such a state's energy and the original energy is

(3.27)\begin{align} {\Delta E_{J\alpha}} & = \frac{2{\rm \pi} m_{\alpha}V}{5}\left(\frac{2}{m_{\alpha}} \right)^{5/2}\left[\bar{\eta} \left(\epsilon_{\mathrm{F}J\alpha}^{5/2} - \epsilon_{\mathrm{F}(J+2)\alpha}^{5/2} \right)\right.\nonumber\\ & \left. \quad -\, \eta_{J\alpha} \left(\epsilon_{\mathrm{F}J\alpha}^{5/2} - \epsilon_{\mathrm{F}(J+1)\alpha}^{5/2} \right) - \eta_{(J+1)\alpha} \left(\epsilon_{\mathrm{F}(J+1)\alpha}^{5/2} - \epsilon_{\mathrm{F}(J+2)\alpha}^{5/2} \right)\right] \nonumber\\ & = \frac{2^{7/2}{\rm \pi} V}{5m_{\alpha}^{3/2}}(\eta_{(J+1)\alpha} - \eta_{J\alpha})\nonumber\\ & \quad\times \left(\frac{\varGamma_{(J+1)\alpha}}{\varGamma_{(J+1)\alpha}+ \varGamma_{J\alpha}}\epsilon_{\mathrm{F}J\alpha}^{5/2} + \frac{\varGamma_{J\alpha}}{\varGamma_{(J+1)\alpha}+ \varGamma_{J\alpha}}\epsilon_{\mathrm{F}(J+2)\alpha}^{5/2} -\epsilon_{\mathrm{F}(J+1)\alpha}^{5/2} \right). \end{align}

The somewhat complex form of (3.27) is due to the fact that the Lynden-Bell equilibria allow for a huge number of possible initial conditions. Nevertheless, some important, and simple, features can be gleaned from (3.27). First, $\Delta E_{J\alpha }$ can be very small if either $\eta _{(J+1)\alpha } - \eta _{J\alpha }$ is small or, more subtly, if one of $\varGamma _{J\alpha }$ or $\varGamma _{(J+1)\alpha }$ is small compared to the other. The first of these possibilities is a manifestation of the fact that, if the phase-space densities of two waterbags are similar, the distribution function does not change significantly by allowing them to intermingle. The second possibility makes the energy required to intermingle the waterbags small because only a small portion of phase space needs to be ‘excited’.

The corollary of (3.27) is, therefore, that the initial energy of the system must only be greater than its minimum possible energy by an amount $E-E_{\mathrm {min}} \sim \mathrm {min}_{J}\Delta E_{J\alpha }$ for the Lynden-Bell equilibrium distribution of species $\alpha$ to be considerably distinct from the Gardner distribution. If, for instance, $\Delta E_{J\alpha }$ were minimal for a particular $J = \bar {J}$, then at energies such that $E-E_{\mathrm {min}} \sim \Delta E_{\bar {J}\alpha }$, the equilibrium would allow the waterbags $\bar {J}$ and $\bar {J}+1$ to be intermingled without the phase space occupied by the other waterbags being ‘unlocked.’ An example of this is shown in figure 3 for a one-dimensional three-waterbag system: panels (a,c) show the exact phase-space densities $f(v)$, panels (b,d) the corresponding mean ones $f_{0}(v)$. For the latter, we also show the contributions $\eta _{J}p_{J}(v)$ from each of the three waterbags. In panels (a,b), where the energy of the system was chosen to be very close to $E_{\mathrm {min}}$, the equilibrium Lynden-Bell looks like a ziggurat. More interestingly, when the system has more energy, as it does in panels (c,d), it first ‘unlocks’ mixing between the two densest waterbags. Thus, the second densest waterbag makes a significant contribution to the phase-space density at zero velocity. The least dense waterbag is still very localised in phase space, indicating that more energy would need to be present in the system to make all waterbags fully non-degenerate. What all this means for Lynden-Bell equilibria is that, despite the relative simplicity of solving a set of transcendental equations enforcing waterbag-content and energy conservation, the shape of these equilibria can be significantly varied for different initial conditions.

Figure 3. A one-dimensional cartoon depicting exact and mean phase-space densities for a three-waterbag system. Panels (b,d) show $f_{0}(v)$ corresponding to the exact $f(v)$ shown in panels (a,c), respectively. The dotted, dashed and dot-dashed lines in panels (b,d) show the contributions $\eta _{J}p_{J}(v)$ of each of the three waterbag densities to the mean phase-space density. The system in the top two panels has a total energy equal to $E = 1.003 E_{\mathrm {min}}$, where $E_{\mathrm {min}}$ is the minimal energy that the system could have subject to phase-volume conservation. The system in the bottom two panels has $E = 1.05 E_{\mathrm {min}}$. This allows the second most dense waterbag (dashed line) to intermingle with the densest waterbag (dotted line) giving a substantial contribution to the mean phase-space density at $v=0$. In contrast, in the top panels, the little energy that the system has above $E_{\mathrm {min}}$ permits very little mixing between waterbags, making the mean phase-space density appear very close to a step function (a ‘ziggurat’).

3.3.2. The H-theorem

To prove the stability of the Lynden-Bell equilibria (3.21), we will now prove an H-theorem for the collision integral (3.19). Namely, we will prove that there is a functional of $f_{0\alpha }$ that can only be increased by evolution under (3.19). This functional is then the entropy of our system, and if it has a maximum, this is a stable attractor of the evolution, since the system could not depart from this state without lowering the entropy.

Consider the following obvious candidate for entropy:

(3.28)\begin{equation} S ={-}\sum_{\alpha}\frac{V}{\Delta\Gamma_{\alpha}}\int\,\mathrm{d}\boldsymbol{v} \,\sum_{J}p_{J\alpha}(\boldsymbol{v})\ln p_{J\alpha}(\boldsymbol{v}). \end{equation}

As before, the sum in this definition includes the empty waterbag $J=0$, whose probability is $p_{0\alpha } = 1 - \sum _{J\neq 0}p_{J\alpha }$. Therefore, (3.28) reduces to the well-known Fermi–Dirac entropy in the single-waterbag case. According to our closure scheme, $p_{J\alpha }(\boldsymbol {v})$ can be written as (3.15). The entropy (3.28) then becomes

(3.29)\begin{align} S & = \sum_{\alpha}\frac{V}{\Delta\Gamma_{\alpha}}\int\,\mathrm{d}\boldsymbol{v} \,\sum_{J}p_{J\alpha}(\boldsymbol{v})\left[\psi_{\alpha}(\boldsymbol{v})\eta_{J\alpha} + \gamma_{J\alpha} + \ln Z_{\alpha}(\psi_{\alpha}(\boldsymbol{v}))\right] \nonumber\\ & = \sum_{\alpha}\frac{V}{\Delta\Gamma_{\alpha}}\left\{\int\,\mathrm{d}\boldsymbol{v} \, \left[f_{0\alpha}(\boldsymbol{v})\psi_{\alpha}(\boldsymbol{v}) + \ln Z_{\alpha}(\psi_{\alpha}(\boldsymbol{v}))\right] + \sum_{J}\gamma_{J\alpha}\frac{n_{J\alpha}}{\eta_{J\alpha}}\right\}. \end{align}

We now take the time derivative of (3.29):

(3.30)\begin{equation} \frac{\mathrm{d}S}{\mathrm{d}t} = \sum_{\alpha}\frac{V}{\Delta\Gamma_{\alpha}}\left[\int\,\mathrm{d}\boldsymbol{v} \, \left(\frac{\partial f_{0\alpha}}{\partial t}\psi_{\alpha} + f_{0\alpha}\frac{\partial \psi_{\alpha}}{\partial t} + \frac{\partial \ln Z_{\alpha}}{\partial t}\right) + \sum_{J}\frac{\partial \gamma_{J\alpha}}{\partial t}\frac{n_{J\alpha}}{\eta_{J\alpha}}\right]. \end{equation}

This can be simplified by taking the time derivative of (3.17) and noting that

(3.31)\begin{align} \frac{\partial \ln Z_{\alpha}}{\partial t} & ={-}\frac{1}{Z_{\alpha}}\sum_{J} {\rm e}^{-\psi_{\alpha}\eta_{J\alpha} - \gamma_{J\alpha}}\left(\frac{\partial \psi_{\alpha}}{\partial t}\eta_{J\alpha} + \frac{\partial \gamma_{J\alpha}}{\partial t} \right)\nonumber\\ & ={-} f_{0\alpha}\frac{\partial \psi_{\alpha}}{\partial t} - \sum_{J}p_{J\alpha}\frac{\partial \gamma_{J\alpha}}{\partial t}. \end{align}

Substituting this into (3.30) and using (3.13), we get

(3.32)\begin{equation} \frac{\mathrm{d}S}{\mathrm{d}t} = \sum_{\alpha}\frac{V}{\Delta\Gamma_{\alpha}}\int\,\mathrm{d}\boldsymbol{v} \,\frac{\partial f_{0\alpha}}{\partial t}\psi_{\alpha}. \end{equation}

In this form, we may finally use the collision integral (3.19) and integrate by parts:

(3.33)\begin{align} {\frac{\mathrm{d}S}{\mathrm{d}t}} & = \sum_{\alpha\alpha''}\frac{16{\rm \pi}^{3}q_{\alpha}^{2}q_{\alpha''}^{2}}{\Delta\Gamma_{\alpha}\Delta\Gamma_{\alpha''}}\sum_{\boldsymbol{k}}\iint \,\mathrm{d}\boldsymbol{v} \,\,\mathrm{d}\boldsymbol{v}'' \,\frac{\delta(\boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{v}-\boldsymbol{v}''))}{k^{4}|\epsilon_{\boldsymbol{k},\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}}|^{2}}\frac{\partial f_{0\alpha}}{\partial \psi_{\alpha}}(\boldsymbol{v})\frac{\partial f_{0\alpha''}}{\partial \psi_{\alpha''}}(\boldsymbol{v}'') \nonumber\\ & \qquad \left[\frac{\Delta\Gamma_{\alpha''}^{2}}{m_{\alpha}^{2}}\left(\boldsymbol{k}\boldsymbol{\cdot}\frac{\partial \psi_{\alpha}}{\partial \boldsymbol{v}} \right)^{2} - \frac{\Delta\Gamma_{\alpha}\Delta\Gamma_{\alpha''}}{m_{\alpha}m_{\alpha''}}\boldsymbol{k}\boldsymbol{\cdot}\frac{\partial \psi_{\alpha''}}{\partial \boldsymbol{v}''}\boldsymbol{k}\boldsymbol{\cdot}\frac{\partial \psi_{\alpha}}{\partial \boldsymbol{v}}\right] \nonumber\\ & = \sum_{\alpha\alpha''}\frac{8{\rm \pi}^{3}q_{\alpha}^{2}q_{\alpha''}^{2}}{\Delta\Gamma_{\alpha}\Delta\Gamma_{\alpha''}}\sum_{\boldsymbol{k}}\iint \,\mathrm{d}\boldsymbol{v} \,\,\mathrm{d}\boldsymbol{v}'' \, \frac{\delta(\boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{v}-\boldsymbol{v}''))}{k^{4}|\epsilon_{\boldsymbol{k},\boldsymbol{k}\boldsymbol{\cdot}{\boldsymbol{v}}}|^{2}}\frac{\partial f_{0\alpha}}{\partial \psi_{\alpha}}(\boldsymbol{v})\frac{\partial f_{0\alpha''}}{\partial \psi_{\alpha''}}(\boldsymbol{v}'')\nonumber\\ & \qquad \left(\frac{\Delta\Gamma_{\alpha''}}{m_{\alpha}}\boldsymbol{k}\boldsymbol{\cdot}\frac{\partial \psi_{\alpha}}{\partial \boldsymbol{v}}- \frac{\Delta\Gamma_{\alpha}}{m_{\alpha''}}\boldsymbol{k}\boldsymbol{\cdot}\frac{\partial \psi_{\alpha''}}{\partial \boldsymbol{v}''} \right)^{2} \geq 0. \end{align}

Here, to get to the second equality, we have symmetrised the expression by swapping ${\alpha \leftrightarrow \alpha ''}$ and $\boldsymbol {v} \leftrightarrow \boldsymbol {v}''$. From (3.18), we note that $\partial f_{0\alpha }/\partial \psi _{\alpha }$ is a negative semi-definite quantity so the integrand in (3.33) is certainly non-negative. This proves that the entropy (3.29) is increased by the collision integral (3.19) unless the integrand is zero, in which case the entropy is conserved. For the latter possibility to be realised, the squared expression in (3.33) must vanish, which it does only for $\psi _{\alpha }(\boldsymbol {v})$ given by (3.20). Thus, the collision integral (3.19) relaxes the mean phase-space density of the system towards the Lynden-Bell multi-waterbag equilibria (3.21).

3.4. Continuous limit of the multi-waterbag formalism

While conceptually enlightening, using a discrete number of waterbag densities is somewhat problematic. At face value, it restricts one to considering initial conditions that are piecewise-constant functions. In fact, a continuous phase-space density $f_{\alpha }(\boldsymbol {v})$ can easily be accommodated by making the grid of $\eta$ arbitrarily fine. We do this by choosing some large number of waterbags densities to fill the interval $[0,\eta _{\alpha \max }]$ with spacing $\Delta \eta \to 0$. We can then rewrite what were sums over $J$ as integrals with respect to the waterbag density $\eta$:

(3.34)\begin{equation} \sum_{J} \to \frac{1}{\Delta\eta}\int_{0}^{\eta_{\alpha \max}}\,\mathrm{d}\eta \,. \end{equation}

Note that the sum above includes the empty waterbag, and so care must be taken to ensure that, when the sums are transformed in this way, the empty waterbag has been included. Now all discrete quantities indexed by $J$ are to be upgraded to continuous functions of $\eta$. Namely, we let

(3.35)\begin{equation} \gamma_{\alpha}(\eta) \equiv \gamma_{J\alpha}, \quad P_{0\alpha}(\boldsymbol{v},\eta) \equiv \frac{p_{J\alpha}(\boldsymbol{v})}{\Delta\eta}, \end{equation}

the latter function being a probability density with respect to $\eta$ (hence the normalisation to $\Delta \eta$). Then (3.15) becomes

(3.36)\begin{gather} P_{0\alpha}(\boldsymbol{v},\eta) = \frac{1}{Z_{\alpha}(\psi_{\alpha}(\boldsymbol{v}))} {\rm e}^{-\psi_{\alpha}(\boldsymbol{v})\eta - \gamma_{\alpha}(\eta)}, \end{gather}

with the partition function (3.16) redefined as

(3.37)\begin{gather} Z_{\alpha}(\psi_{\alpha}(\boldsymbol{v})) = \int_{0}^{\eta_{\alpha\mathrm{max}}}\,\mathrm{d}\eta \, {\rm e}^{-\psi_{\alpha}(\boldsymbol{v})\eta - \gamma_{\alpha}(\eta)}. \end{gather}

The Lagrange multiplier $\psi _{\alpha }(\boldsymbol {v})$ is now determined by the continuous version of (3.10) and (3.17):

(3.38)\begin{equation} f_{0\alpha}(\boldsymbol{v}) = \int\,\mathrm{d}\eta \, \eta P_{0\alpha}(\boldsymbol{v},\eta) ={-}\frac{\partial \ln Z_{\alpha}}{\partial \psi_{\alpha}}. \end{equation}

To determine $\gamma _{\alpha }(\eta )$, one must solve the continuous version of the constraint (3.13) fixing the ‘waterbag content’ of the distribution:

(3.39)\begin{equation} \int \,\mathrm{d}\boldsymbol{v} \, P_{0\alpha}(\boldsymbol{v},\eta) = \frac{n_{\alpha}(\eta)}{\eta \Delta\eta} \equiv \rho_{\alpha}(\eta), \end{equation}

where $n_{\alpha }(\eta )$ is the continuous generalisation of $n_{J\alpha }$, and the newly defined function $\rho _{\alpha }(\eta )$ encodes all the information about the species $\alpha$ that must be retained from its initial distribution (an infinite set of Casimir invariants).

Thus, the price of the generalisation to a continuous-waterbag model is having to solve the two coupled integral equations (3.38) and (3.39) for the functions $\psi _{\alpha }(\boldsymbol {v})$ and $\gamma _{\alpha }(\eta )$. The collision integral is then again (3.19), its fixed points are (3.38) with $\psi _{\alpha }(\boldsymbol {v})$ given by (3.20), viz.,

(3.40)\begin{gather} f_{0\alpha}(\boldsymbol{v}) = \frac{\int\,\mathrm{d}\eta \,\eta {\rm e}^{-\psi_{\alpha}(\boldsymbol{v})\eta - \gamma_{\alpha}(\eta)}}{\int \,\mathrm{d}\eta \, {\rm e}^{-\psi_{\alpha}(\boldsymbol{v})\eta - \gamma_{\alpha}(\eta)}}, \end{gather}

and the H-theorem (3.33) continues to hold, with the continuous limit of the entropy (3.28):

(3.41)\begin{equation} S ={-}\sum_{\alpha}\frac{V}{\Delta\Gamma_{\alpha}}\iint\,\mathrm{d}\boldsymbol{v} \,\,\mathrm{d}\eta \,P_{0\alpha}(\boldsymbol{v},\eta)\ln P_{0\alpha}(\boldsymbol{v},\eta). \end{equation}

4.1. Hyperkinetic collision integral

To construct the kinetics of the individual waterbags, we first define the ‘waterbag distribution function’

(4.1)\begin{equation} P_{\alpha}(\boldsymbol{r},\boldsymbol{v},\eta) = \delta(f_{\alpha}(\boldsymbol{r},\boldsymbol{v}) - \eta), \end{equation}

which is the probability density of finding the exact phase-space density $f_{\alpha }(\boldsymbol {r},\boldsymbol {v})$ to have the value $\eta$ at the phase-space position $(\boldsymbol {r},\boldsymbol {v})$. The evolution equation for $P_{\alpha }$ takes the same form as the Vlasov equation:

(4.2)\begin{align} \frac{\partial P_{\alpha}}{\partial t} & = \frac{\partial f_{\alpha}}{\partial t}\delta'(f_{\alpha}(\boldsymbol{r},\boldsymbol{v}) - \eta) \nonumber\\ & = \left(-\boldsymbol{v}\boldsymbol{\cdot}\nabla f_{\alpha} + \frac{q_{\alpha}}{m_{\alpha}}\nabla \varphi \boldsymbol{\cdot} \frac{\partial f_{\alpha}}{\partial \boldsymbol{v}}\right)\delta'(f_{\alpha}(\boldsymbol{r},\boldsymbol{v}) - \eta) \nonumber\\ & ={-}\boldsymbol{v}\boldsymbol{\cdot}\nabla P_{\alpha} + \frac{q_{\alpha}}{m_{\alpha}}\nabla \varphi \boldsymbol{\cdot} \frac{\partial P_{\alpha}}{\partial \boldsymbol{v}}. \end{align}

Since we are still concerned with the relaxation to equilibrium, we now wish to find the collision integral for the evolution of $P_{\alpha }$ in much the same way as we did for $f_{\alpha }$ in § 2.

As before, we Fourier decompose

(4.3)\begin{gather} P_{\alpha}(\boldsymbol{r},\boldsymbol{v},\eta) = P_{0\alpha}(\boldsymbol{v},\eta) + \sum_{\boldsymbol{k}} {\rm e}^{{\rm i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{r}}P_{\boldsymbol{k}\alpha}(\boldsymbol{v},\eta), \end{gather}

and expect the ensemble average of the homogeneous part of the hyperkinetic distribution, $P_{0\alpha }$, to be much greater in size than the fluctuating part, sanctioning the quasilinear approach. Therefore, we again linearise our hyperkinetic equation (4.2) to get the fluctuating part of the distribution:

(4.4)\begin{equation} \frac{\partial P_{\boldsymbol{k}\alpha}}{\partial t} + {\rm i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}P_{0\alpha} = {\rm i}\frac{q_{\alpha}}{m_{\alpha}}\varphi_{\boldsymbol{k}}\boldsymbol{k}\boldsymbol{\cdot}\frac{\partial P_{0\alpha}}{\partial \boldsymbol{v}}, \end{equation}

which in turn determines the evolution of the homogeneous part:

(4.5)\begin{equation} \frac{\partial P_{0\alpha}}{\partial t} = \frac{\partial }{\partial \boldsymbol{v}}\boldsymbol{\cdot}\left(\frac{q_{\alpha}}{m_{\alpha}}\sum_{\boldsymbol{k}}\boldsymbol{k}\mathrm{Im}\langle \varphi_{\boldsymbol{k}}^{*}P_{\boldsymbol{k}\alpha} \rangle \right). \end{equation}

The only difference with the calculation in § 2 comes from the fact that electric-field perturbations are still made by charge-density perturbations, which are a cumulative effect of all the waterbags. In other words, Poisson's equation (2.6) is

(4.6)\begin{equation} \varphi_{\boldsymbol{k}} = \sum_{\alpha}\frac{4{\rm \pi} q_{\alpha}}{k^{2}}\int \,\mathrm{d}\boldsymbol{v} \,f_{\boldsymbol{k}\alpha}(\boldsymbol{v}) =\sum_{\alpha}\frac{4{\rm \pi} q_{\alpha}}{k^{2}}\iint \,\mathrm{d}\boldsymbol{v} \,\,\mathrm{d}\eta \, \eta P_{\boldsymbol{k}\alpha}(\boldsymbol{v},\eta). \end{equation}

The derivation of the collision integral for $P_{0\alpha }$ can now be ported over from § 2 near verbatim. The only difference comes from Poisson's equation (4.6): instead of velocity integrals, we now have integrals over both velocities and waterbag densities

(4.7)\begin{equation} \int\,\mathrm{d}\boldsymbol{v} \,\left(\cdots\right) \to \iint\,\mathrm{d}\boldsymbol{v} \,\,\mathrm{d}\eta \, \eta\left(\cdots\right). \end{equation}

To avoid repetition, we will just note the key equations that are reached throughout the derivation, with reference to the corresponding equations in § 2.

The dielectric function, previously (2.10), is now

(4.8)\begin{equation} \epsilon_{\boldsymbol{k}}(p) = 1- {\rm i}\sum_{\alpha'}\frac{4{\rm \pi} q_{\alpha'}^{2}}{m_{\alpha'}k^{2}}\iint\,\mathrm{d}\boldsymbol{v}' \,\,\mathrm{d}\eta' \,\frac{\eta'}{p+{\rm i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}'}\boldsymbol{k}\boldsymbol{\cdot}\frac{\partial P_{0\alpha'}}{\partial \boldsymbol{v}'}. \end{equation}

The resultant general form of the collision integral, previously (2.18), becomes

(4.9)\begin{gather} \frac{\partial P_{0\alpha}}{\partial t} = \frac{\partial }{\partial \boldsymbol{v}}\boldsymbol{\cdot}\sum_{\alpha''}\!\iint\,\mathrm{d}\boldsymbol{v}'' \,\,\mathrm{d}\eta'' \,\eta'' \left[\boldsymbol{\mathsf{D}}^{\alpha}_{\alpha\alpha''}(\boldsymbol{v},\boldsymbol{v}'',\eta'')\,\boldsymbol{\cdot}\,\left.\frac{\partial P_{0\alpha}}{\partial \boldsymbol{v}}\right|_{\eta} \!-\! \boldsymbol{\mathsf{D}}^{\alpha}_{\alpha''\alpha}(\boldsymbol{v}'',\boldsymbol{v},\eta)\boldsymbol{\cdot} \left. \frac{\partial P_{0\alpha''}}{\partial \boldsymbol{v}''}\right|_{\eta''} \right]\!. \end{gather}

The diffusion kernel, via an expression analogous to (2.19), can again be simplified by assuming that the fluctuations evolve much more rapidly than the mean: a calculation identical to that done in § 2.3 returns the following analogue of (2.24):

(4.10)\begin{align} & \boldsymbol{\mathsf{D}}^{\alpha}_{\mu\nu}(\boldsymbol{w},\boldsymbol{v},\eta) = \sum_{\nu'}\frac{16{\rm \pi}^{3}q_{\mu}^{2}q_{\nu}q_{\nu'}}{m_{\alpha}m_{\mu}}\nonumber\\ & \qquad \mathrm{Re}\sum_{\boldsymbol{k}}\frac{\boldsymbol{k}\boldsymbol{k}}{k^{4}}\delta(\boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{w} - \boldsymbol{v}))\iint\,\mathrm{d}\boldsymbol{v}' \,\,\mathrm{d}\eta' \,\eta' {\rm e}^{-{\rm i}\boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{v}-\boldsymbol{v}')t}\frac{\langle g_{\boldsymbol{k}\nu}(\boldsymbol{v},\eta)g^{*}_{\boldsymbol{k}\nu'}(\boldsymbol{v}',\eta') \rangle}{\epsilon_{\boldsymbol{k},\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}}\epsilon^{*}_{\boldsymbol{k},\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}'}}, \end{align}

where $g_{\boldsymbol {k}\nu }(\boldsymbol {v},\eta ) = P_{\boldsymbol {k}\nu }(t=0,\boldsymbol {v},\eta )$ is the initial distribution and, as in § 2, $\epsilon _{\boldsymbol {k},\boldsymbol {k}\boldsymbol {\cdot }\boldsymbol {v}} \equiv \epsilon _{\boldsymbol {k}}(-{\rm i}\boldsymbol {k}\boldsymbol {\cdot }\boldsymbol {v})$.

Now we need a closure for the correlation function $\langle g_{\boldsymbol {k}\nu }(\boldsymbol {v},\eta )g^{*}_{\boldsymbol {k}\nu '}(\boldsymbol {v}',\eta ') \rangle$ in terms of $P_{0\nu }(\boldsymbol {v},\eta )$. The first step is again the microgranulation ansatz introduced in § 2.5, viz., the assumption that only very near points in phase space are correlated with each other, over a phase-space volume $\Delta \Gamma _{\nu }$:

(4.11)\begin{equation} \langle g_{\nu}(\boldsymbol{r},\boldsymbol{v},\eta)g_{\nu'}(\boldsymbol{r}',\boldsymbol{v}',\eta') \rangle = \Delta\Gamma_{\nu}\delta_{\nu\nu'}\langle g_{\nu}(\eta)g_{\nu}(\eta') \rangle(\boldsymbol{v})\delta(\boldsymbol{v} - \boldsymbol{v}')\delta(\boldsymbol{r} - \boldsymbol{r}'), \end{equation}

or, in Fourier space,

(4.12)\begin{equation} \langle g_{\boldsymbol{k}\nu}(\boldsymbol{v},\eta)g^{*}_{\boldsymbol{k}\nu'}(\boldsymbol{v}',\eta') \rangle = \frac{\Delta\Gamma_{\nu}}{V}\delta_{\nu\nu'}\langle g_{\nu}(\eta)g_{\nu}(\eta') \rangle(\boldsymbol{v})\delta(\boldsymbol{v} - \boldsymbol{v}'). \end{equation}

This is the generalisation of (2.32). Note that, while localising the correlator in the $(\boldsymbol {r},\boldsymbol {v})$ space, we allow for correlations between different values of $\eta$. To determine this remaining correlator in terms of $P_{0\nu }$, we use (4.1) to find

(4.13)\begin{align} \langle g_{\nu}(\eta)g_{\nu}(\eta') \rangle(\boldsymbol{v}) & = \langle [\delta(f_{\nu}(\boldsymbol{r},\boldsymbol{v}) - \eta) - P_{0\nu}(\boldsymbol{v},\eta)][\delta(f_{\nu}(\boldsymbol{r},\boldsymbol{v}) - \eta') - P_{0\nu}(\boldsymbol{v},\eta')] \rangle \nonumber\\ & = \langle \delta(f_{\nu}(\boldsymbol{r},\boldsymbol{v}) - \eta)\delta(f_{\nu}(\boldsymbol{r},\boldsymbol{v}) - \eta') \rangle - P_{0\nu}(\boldsymbol{v},\eta)P_{0\nu}(\boldsymbol{v},\eta')\nonumber\\ & = \delta(\eta - \eta')P_{0\nu}(\boldsymbol{v},\eta) - P_{0\nu}(\boldsymbol{v},\eta)P_{0\nu}(\boldsymbol{v},\eta'). \end{align}

This formula is a direct generalisation of the single-waterbag closure (3.4), but of course there is no longer any need to assume a single-waterbag distribution. Neither is there any need to assume a local maximisation of entropy, as we did for our multi-waterbag distribution in § 3.2: working in the extend phase space $(\boldsymbol {r},\boldsymbol {v},\eta )$ has resolved the closure problem for $\langle f_{\alpha }^{2} \rangle$ vs $f_{0\alpha }$ automatically.

Using (4.12) and (4.13), we can now write the diffusion kernel (4.10) in a closed form:

(4.14)\begin{equation} \boldsymbol{\mathsf{D}}^{\,\alpha}_{\mu\nu}(\boldsymbol{w},\boldsymbol{v},\eta) = \frac{16{\rm \pi}^{3}q_{\mu}^{2}q_{\nu}^{2}\Delta\Gamma_{\nu}}{m_{\alpha}m_{\mu}V} \sum_{\boldsymbol{k}}\frac{\boldsymbol{k}\boldsymbol{k}}{k^{4}}\frac{\delta(\boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{w}- \boldsymbol{v}))}{|\epsilon_{\boldsymbol{k},\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}}|^{2}} \left[\eta - f_{0\nu}(\boldsymbol{v}) \right]P_{0\nu}(\boldsymbol{v},\eta), \end{equation}

where $f_{0\nu }(\boldsymbol {v}) = \int \,\mathrm {d}\eta ' \,\eta 'P_{0\nu }(\boldsymbol {v},\eta ')$. Finally, substituting (4.14) into (4.9) gives us the generalised hyperkinetic collision integral:

(4.15)\begin{gather} \frac{\partial P_{0\alpha}}{\partial t} = \sum_{\alpha''}\frac{16{\rm \pi}^{3}q_{\alpha}^{2}q_{\alpha''}^{2}}{m_{\alpha}V}\frac{\partial }{\partial \boldsymbol{v}}\boldsymbol{\cdot}\sum_{\boldsymbol{k}}\frac{\boldsymbol{k}\boldsymbol{k}}{k^{4}}\boldsymbol{\cdot}\int\,\mathrm{d}\boldsymbol{v}'' \,\frac{\delta(\boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{v}-\boldsymbol{v}''))}{|\epsilon_{\boldsymbol{k},\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}}|^{2}}\int\,\mathrm{d}\eta'' \,\eta''\nonumber\\ \left\{ \frac{\Delta\Gamma_{\alpha''}}{m_{\alpha}}\left[ \eta'' - f_{0\alpha''}(\boldsymbol{v}'') \right]P_{0\alpha''}(\boldsymbol{v}'',\eta'')\left.\frac{\partial P_{0\alpha}}{\partial \boldsymbol{v}}\right|_{\eta} - \frac{\Delta\Gamma_{\alpha}}{m_{\alpha''}}\left[ \eta - f_{0\alpha}(\boldsymbol{v}) \right]P_{0\alpha}(\boldsymbol{v},\eta)\left.\frac{\partial P_{0\alpha''}}{\partial \boldsymbol{v}''}\right|_{\eta''}\right\}. \end{gather}

To reiterate, the only closure required in the derivation of this collision integral was the microgranulation ansatz (4.11). No need to break higher-order correlators has arisen, because, within the hyperkinetic formalism, all moments of the phase-space density can be derived from the distribution function $P_{0\alpha }$:$\langle f_{\alpha }^{n} \rangle = \int \,\mathrm {d}\eta \,\eta ^{n}P_{0\alpha }(\boldsymbol {v},\eta )$. Indeed, if one takes the first moment of (4.15) with respect to $\eta$, it is easy to show that one recovers (2.35), which was the general collision integral before a choice of waterbag closure was made in § 3.

As well as conserving total energy, momentum, and particle number, (4.15) has two further invariants, which represent the conservation of phase volume,

(4.16)\begin{equation} \int \,\mathrm{d}\boldsymbol{v} \, P_{0\alpha}(\boldsymbol{v},\eta) = \rho_{\alpha}(\eta), \end{equation}

and the conservation of probability

(4.17)\begin{equation} \int\,\mathrm{d}\eta \,P_{0\alpha}(\boldsymbol{v},\eta) = 1. \end{equation}

As in (3.39), (4.16) distils to a single function $\rho _{\alpha }(\eta )$ (the ‘waterbag content’ of the distribution, or the infinite set of its Casimir invariants) all the information from the initial condition that must be preserved by collisionless evolution. Given these invariants, it is unsurprising that the fixed points of the collision integral (4.15) will again be the Lynden-Bell equilibria, viz.,

(4.18)\begin{gather} P_{0\alpha}(\boldsymbol{v},\eta) = \frac{{\rm e}^{-\beta\Delta\Gamma_{\alpha}\eta\left[\epsilon_{\alpha} - \mu_{\alpha}(\eta)\right]}}{\int\,\mathrm{d}\eta \, {\rm e}^{-\beta\Delta\Gamma_{\alpha}\eta\left[\epsilon_{\alpha} - \mu_{\alpha}(\eta)\right]}}, \end{gather}

where $\epsilon _{\alpha } = m_{\alpha }|\boldsymbol {v}|^{2}/2$, and the normalisation (4.17) has been enforced. The parameters $\beta$ and $\mu (\eta )$ are determined from (4.16) and the energy conservation,

(4.19)\begin{equation} V\iint\,\mathrm{d}\boldsymbol{v} \,\,\mathrm{d}\eta \,\eta \epsilon_{\alpha}(\boldsymbol{v})P_{0\alpha}(\boldsymbol{v},\eta) = E. \end{equation}

That (4.18) are indeed fixed points of (4.15) can be confirmed by direct substitution. Just as in § 3.3.2, we must provide an H-theorem to prove their stability.

4.2. Hyperkinetic H-theorem

We define the Shannon entropy as before but now upgraded to its continuous variant (3.41). Taking the time derivative of $S$ using (4.15) and integrating by parts in $\boldsymbol {v}$, we get

(4.20)\begin{align} {\frac{\mathrm{d}S}{\mathrm{d}t}} & ={-}\sum_{\alpha}\frac{V}{\Delta\Gamma_{\alpha}}\iint\,\mathrm{d}\boldsymbol{v} \,\,\mathrm{d}\eta \,\left(1 + \ln P_{0\alpha} \right)\frac{\partial P_{0\alpha}}{\partial t} \nonumber\\ & = \sum_{\alpha\alpha''}\frac{16{\rm \pi}^{3}q_{\alpha}^{2}q_{\alpha''}^{2}}{\Delta\Gamma_{\alpha}\Delta\Gamma_{\alpha''}}\iint\,\mathrm{d}\boldsymbol{v} \,\,\mathrm{d}\boldsymbol{v}'' \,\sum_{\boldsymbol{k}}\frac{\delta(\boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{v}-\boldsymbol{v}''))}{k^{4}|\epsilon_{\boldsymbol{k},\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}}|^{2}}\nonumber\\ & \quad \iint\,\mathrm{d}\eta \,\,\mathrm{d}\eta'' \,\eta'' \left\{ \frac{\Delta\Gamma_{\alpha''}^{2}}{m_{\alpha}^{2}}\left[\eta''- f_{0\alpha''}(\boldsymbol{v}'')\right]\frac{P_{0\alpha''}(\boldsymbol{v}'',\eta'')}{P_{0\alpha}(\boldsymbol{v},\eta)}\left(\boldsymbol{k}\boldsymbol{\cdot}\left.\frac{\partial P_{0\alpha}}{\partial \boldsymbol{v}}\right|_{\eta} \right)^{2}\right. \nonumber\\ & \left.\qquad - \,\frac{\Delta\Gamma_{\alpha}\Delta\Gamma_{\alpha''}}{m_{\alpha}m_{\alpha''}}\left[\eta - f_{0\alpha}(\boldsymbol{v})\right]\boldsymbol{k}\boldsymbol{\cdot}\left. \frac{\partial P_{0\alpha}}{\partial \boldsymbol{v}} \right|_{\eta}\boldsymbol{k}\boldsymbol{\cdot} \left. \frac{\partial P_{0\alpha}}{\partial \boldsymbol{v}''}\right|_{\eta''} \vphantom{\left(\boldsymbol{k}\boldsymbol{\cdot}\left.\frac{\partial P_{0\alpha}}{\partial \boldsymbol{v}}\right|_{\eta} \right)^{2}}\right\}. \end{align}

First, we notice that the $\eta ''$ integral of the first bracketed term is

(4.21)\begin{equation} \int \,\mathrm{d}\eta'' \,\eta'' \left[\eta'' - f_{0\alpha''}(\boldsymbol{v}'')\right]P_{0\alpha''}(\boldsymbol{v}'',\eta'') = \int\,\mathrm{d}\eta'' \,\left[\eta'' - f_{0\alpha''}(\boldsymbol{v}'') \right]^{2}P_{0\alpha''}(\boldsymbol{v}'',\eta''). \end{equation}

Secondly, in the second bracketed term, anything that is not multiplied by both $\eta$ and $\eta ''$ integrates to zero by (4.17):

(4.22)\begin{equation} \int\,\mathrm{d}\eta \,\left.\frac{\partial P_{0\alpha}}{\partial \boldsymbol{v}}\right|_{\eta} = \frac{\partial }{\partial \boldsymbol{v}}\int \,\mathrm{d}\eta \,P_{0\alpha}(\boldsymbol{v},\eta) = 0. \end{equation}

With these insights, the rate of change of entropy becomes

(4.23)\begin{align} {\frac{\mathrm{d}S}{\mathrm{d}t}} & = \sum_{\alpha\alpha''}\frac{16{\rm \pi}^{3}q_{\alpha}^{2}q_{\alpha''}^{2}}{\Delta\Gamma_{\alpha}\Delta\Gamma_{\alpha''}}\iint\,\mathrm{d}\boldsymbol{v} \,\,\mathrm{d}\boldsymbol{v}'' \,\sum_{\boldsymbol{k}}\frac{\delta(\boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{v}-\boldsymbol{v}''))}{k^{4}|\epsilon_{\boldsymbol{k},\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}}|^{2}}\nonumber\\ & \qquad \iint\,\mathrm{d}\eta \,\,\mathrm{d}\eta'' \, \left\{ \frac{\Delta\Gamma_{\alpha''}^{2}}{m_{\alpha}^{2}}\left[\eta''- f_{0\alpha''}(\boldsymbol{v}'')\right]^{2}\frac{P_{0\alpha''}(\boldsymbol{v}'',\eta'')}{P_{0\alpha}(\boldsymbol{v},\eta)}\left(\boldsymbol{k}\boldsymbol{\cdot}\left.\frac{\partial P_{0\alpha}}{\partial \boldsymbol{v}}\right|_{\eta} \right)^{2} \right.\nonumber\\ & \left.\qquad - \frac{\Delta\Gamma_{\alpha}\Delta\Gamma_{\alpha''}}{m_{\alpha}m_{\alpha''}}\eta \eta''\boldsymbol{k}\boldsymbol{\cdot}\left. \frac{\partial P_{0\alpha}}{\partial \boldsymbol{v}} \right|_{\eta}\boldsymbol{k}\boldsymbol{\cdot} \left. \frac{\partial P_{0\alpha''}}{\partial \boldsymbol{v}''}\right|_{\eta''} \right\}. \end{align}

Finally, we symmetrise the entire expression by swapping $\alpha \leftrightarrow \alpha ''$, $\eta \leftrightarrow \eta ''$, $\boldsymbol {v} \leftrightarrow \boldsymbol {v}''$, which allows us to write (4.23) in an explicitly positive semi-definite form:

(4.24)\begin{align} {\frac{\mathrm{d}S}{\mathrm{d}t}} & = \sum_{\alpha\alpha''}\frac{8{\rm \pi}^{3}q_{\alpha}^{2}q_{\alpha''}^{2}}{\Delta\Gamma_{\alpha}\Delta\Gamma_{\alpha''}}\iint\,\mathrm{d}\boldsymbol{v} \,\,\mathrm{d}\boldsymbol{v}'' \,\sum_{\boldsymbol{k}}\frac{\delta(\boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{v}-\boldsymbol{v}''))}{k^{4}|\epsilon_{\boldsymbol{k},\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}}|^{2}}\nonumber\\ & \qquad \iint\,\mathrm{d}\eta \,\,\mathrm{d}\eta'' \,\left\{\frac{\Delta\Gamma_{\alpha''}}{m_{\alpha}}\left[\eta'' - f_{0\alpha''}(\boldsymbol{v}'')\right]\sqrt{\frac{P_{0\alpha''}(\boldsymbol{v}'',\eta'')}{P_{0\alpha}(\boldsymbol{v},\eta)}}\boldsymbol{k}\boldsymbol{\cdot} \left.\frac{\partial P_{0\alpha}}{\partial \boldsymbol{v}}\right|_{\eta} \right.\nonumber\\ & \left.\qquad -\,\frac{\Delta\Gamma_{\alpha}}{m_{\alpha''}}\left[\eta - f_{0\alpha}(\boldsymbol{v})\right]\sqrt{\frac{P_{0\alpha}(\boldsymbol{v},\eta)}{P_{0\alpha''}(\boldsymbol{v}'',\eta'')}}\boldsymbol{k}\boldsymbol{\cdot} \left.\frac{\partial P_{0\alpha''}}{\partial \boldsymbol{v}''}\right|_{\eta''} \right\}^{2} \geq 0. \end{align}

Expanding the squared expression in (4.24) does indeed recover (4.23) as all excess terms cancel. The collision integral (4.15) therefore never decreases the Shannon entropy (3.41).

To prove that all initial conditions will eventually reach their corresponding Lynden-Bell equilibrium (4.18), we must show that the rate of entropy growth (4.24) will equal zero if and only if the Lynden-Bell equilibrium has been reached. Owing to the squared expression in (4.24), the entropy growth will equal zero when

(4.25)\begin{align} & \frac{1}{\Delta\Gamma_{\alpha}m_{\alpha}\left[\eta - f_{0\alpha}(\boldsymbol{v}) \right]P_{0\alpha}(\boldsymbol{v},\eta)}\boldsymbol{k}\boldsymbol{\cdot}\left.\frac{\partial P_{0\alpha}}{\partial \boldsymbol{v}}\right|_{\eta} \nonumber=\\& \qquad\qquad\qquad\qquad\qquad \frac{1}{\Delta\Gamma_{\alpha''}m_{\alpha''}\left[\eta'' - f_{0\alpha''}(\boldsymbol{v}'') \right]P_{0\alpha''}(\boldsymbol{v}'',\eta'')}\boldsymbol{k}\boldsymbol{\cdot}\left.\frac{\partial P_{0\alpha''}}{\partial \boldsymbol{v}''}\right|_{\eta''}, \end{align}

for all $\alpha$, $\alpha ''$, $\eta$ and $\eta ''$ when $\boldsymbol {k}\boldsymbol {\cdot }(\boldsymbol {v}- \boldsymbol {v}'' ) = 0$. Save for $\boldsymbol {v}$ and $\boldsymbol {v}''$, the left- and right-hand sides of (4.25) are functions of distinct, independent variables. The only solution to (4.25) must, therefore, satisfy

(4.26)\begin{equation} \boldsymbol{k}\boldsymbol{\cdot}\left.\frac{\partial P_{0\alpha}}{\partial \boldsymbol{v}}\right|_{\eta} ={-}\beta\Delta\Gamma_{\alpha}m_{\alpha}\left[\eta - f_{0\alpha}(\boldsymbol{v}) \right]P_{0\alpha}(\boldsymbol{v},\eta)\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}, \end{equation}

where $\beta$ is an as yet undetermined constant, although it is clear it will come to mean the thermodynamic beta shortly. Note that by writing $\boldsymbol {k}\boldsymbol {\cdot }\boldsymbol {v}$ we have assumed, without loss of generality, that we are in the zero-net-momentum frame. It is clear that the Lynden-Bell equilibria (4.18) satisfy the condition (4.26) [as they must, being maximisers of the entropy (3.41)].

To prove that these are the only solutions, we rewrite (4.26) as

(4.27)\begin{equation} \left.\frac{\partial \ln P_{0\alpha}}{\partial \boldsymbol{v}}\right|_{\eta} ={-}\beta\Delta\Gamma_{\alpha}\eta m_{\alpha}\boldsymbol{v} + \beta m_{\alpha}\Delta\Gamma_{\alpha}\int\,\mathrm{d}\eta \,\eta P_{0\alpha}(\boldsymbol{v},\eta)\boldsymbol{v}. \end{equation}

Then, without loss of generality,

(4.28)\begin{equation} P_{0\alpha}(\boldsymbol{v},\eta) = \frac{C_{\alpha}(\boldsymbol{v},\eta)}{\int \,\mathrm{d}\eta \,C_{\alpha}(\boldsymbol{v},\eta)}, \end{equation}

for which (4.27) becomes

(4.29)\begin{align} \left.\frac{\partial \ln C_{\alpha}}{\partial \boldsymbol{v}}\right|_{\boldsymbol{v}} + \beta \Delta\Gamma_{\alpha}\eta m_{\alpha}\boldsymbol{v} & = \beta\Delta\Gamma_{\alpha}m_{\alpha}\frac{\int\,\mathrm{d}\eta \,\eta C_{\alpha}(\boldsymbol{v},\eta)}{\int\,\mathrm{d}\eta \,C_{\alpha}(\boldsymbol{v},\eta)}\boldsymbol{v} + \frac{\partial }{\partial \boldsymbol{v}}\ln\int\,\mathrm{d}\eta \,C_{\alpha}(\boldsymbol{v},\eta)\boldsymbol{v}\nonumber\\ & = \frac{\partial }{\partial \boldsymbol{v}}\varPhi_{\alpha}(\boldsymbol{v}). \end{align}

In (4.29), we have collected all $\eta$-dependent terms on the right-hand side and further declared that, since the curl of the left-hand side is zero, these terms can be written as the gradient of some function $\varPhi _{\alpha }(\boldsymbol {v})$, which will be determined self-consistently after determining $C_{\alpha }(\boldsymbol {v})$. With this sleight of hand, (4.29) is solved by

(4.30)\begin{gather} C_{\alpha}(\boldsymbol{v},\eta) = {\rm e}^{\varPhi_{\alpha}(\boldsymbol{v})} {\rm e}^{-\beta\Delta\Gamma_{\alpha}\eta\left[\epsilon_{\alpha}(\boldsymbol{v}) - \mu_{\alpha}(\eta)\right]}, \end{gather}

where $\epsilon _{\alpha }(\boldsymbol {v}) = m_{\alpha }|\boldsymbol {v}|^{2}/2$ and the chemical potential $\mu _{\alpha }(\eta )$ has emerged as an integration constant. Now, substituting this into (4.28), we find that $P_{0\alpha }(\boldsymbol {v},\eta )$ is the Lynden-Bell equilibrium (4.18), proving that it is the only solution for which the entropy does not grow. Note that the determination of $\varPhi _{\alpha }(\boldsymbol {v})$ is unimportant for the calculation of $P_{0\alpha }(\boldsymbol {v},\eta )$ because $C_{\alpha }(\boldsymbol {v})$ can be freely multiplied by any function of $\boldsymbol {v}$ without changing $P_{0\alpha }(\boldsymbol {v},\eta )$; also, given (4.30), the solvability condition

(4.31)\begin{equation} \frac{\partial }{\partial \boldsymbol{v}}\varPhi_{\alpha}(\boldsymbol{v}) = \beta\Delta\Gamma_{\alpha}m_{\alpha}\frac{\int\,\mathrm{d}\eta \,\eta C_{\alpha}(\boldsymbol{v},\eta)}{\int\,\mathrm{d}\eta \,C_{\alpha}(\boldsymbol{v},\eta)}\boldsymbol{v} + \frac{\partial }{\partial \boldsymbol{v}}\ln\int\,\mathrm{d}\eta \,C_{\alpha}(\boldsymbol{v},\eta)\boldsymbol{v}, \end{equation}

is satisfied for all functions $\varPhi _{\alpha }(\boldsymbol {v})$. This completes the proof that all initial conditions will reach their Lynden-Bell equilibria.

4.3. Relation between multi-waterbag and hyperkinetic collision integrals

In this section, we aim to draw a relation between the two collision integrals derived above: the multi-waterbag collision integral (3.19) and the hyperkinetic collision integral (4.15). The apparent distinction between these collision integrals is the need for an artificial closure in the multi-waterbag collision integral that is not present in the hyperkinetic collision integral. We will begin to shed light on this by stating the following exact relation between the two collision integrals:

(i) The hyperkinetic collision integral recovers the multi-waterbag collision integral if the waterbag distribution function $P_{0\alpha }(\boldsymbol {v},\eta )$ satisfies the closure relation (3.36). To prove this, we observe, as we did already in section (4.1), that the first moment of the hyperkinetic collision integral (4.15) with respect to $\eta$ is (2.35) with $\langle g_{\alpha }^{2} \rangle (\boldsymbol {v})$ given by

(4.32)\begin{equation} \langle g_{\alpha}^{2} \rangle(\boldsymbol{v}) = \int\,\mathrm{d}\eta \,\left[\eta - f_{0\alpha}(\boldsymbol{v})\right]^{2}P_{0\alpha}(\boldsymbol{v},\eta) = \int \,\mathrm{d}\eta \, \eta^{2}P_{0\alpha}(\boldsymbol{v},\eta) - \left[\int \,\mathrm{d}\eta \, \eta P_{0\alpha}(\boldsymbol{v},\eta) \right]^{2}. \end{equation}

Substituting $P_{0\alpha }(\boldsymbol {v},\eta )$ from (3.36) into (4.32) gives

(4.33)\begin{gather} \langle g_{\alpha}^{2} \rangle(\boldsymbol{v}) = \frac{1}{Z_{\alpha}(\psi_{\alpha})}\int\,\mathrm{d}\eta \,\eta^{2} {\rm e}^{-\psi_{\alpha}(\boldsymbol{v})\eta - \gamma_{\alpha}(\eta)} - f_{0\alpha}^{2}(\boldsymbol{v}), \end{gather}

which is the continuous version of (3.18) used to arrive at the multi-waterbag collision integral. Therefore, if the waterbag distribution function satisfies (3.36), the hyperkinetic collision integral reduces to the multi-waterbag collision integral. However, this is not an automatic guarantee that the waterbag distribution function $P_{0\alpha }$ will ‘stay on the closure’. In this case, however, it is possible to prove that it will.

(ii) If the waterbag distribution function $P_{0\alpha }(\boldsymbol {v},\eta )$ evolving under the hyperkinetic collision integral (4.15) ever satisfies the closure (3.36), then it will continue to satisfy it for all future times. For a given waterbag distribution function $P_{0\alpha }(\boldsymbol {v},\eta )$, let us define $\bar {P}_{0\alpha }(\boldsymbol {v},\eta )$ to be the waterbag distribution function that maximises the entropy (3.41) subject to having the same mean phase-space density $f_{0\alpha }(\boldsymbol {v})$ and waterbag content $\rho _{\alpha }(\eta )$ as $P_{0\alpha }(\boldsymbol {v},\eta )$, i.e., $\bar {P}_{0\alpha }(\boldsymbol {v},\eta )$ is the waterbag distribution function (3.36). Now, we examine the difference between the entropies of $\bar {P}_{0\alpha }(\boldsymbol {v},\eta )$ and $P_{0\alpha }(\boldsymbol {v},\eta )$

(4.34)\begin{equation} \bar{S} - S = \sum_{\alpha}\frac{V}{\Delta\Gamma_{\alpha}}\iint\,\mathrm{d}\boldsymbol{v} \,\,\mathrm{d}\eta \, \left(P_{0\alpha}\ln P_{0\alpha} - \bar{P}_{0\alpha}\ln\bar{P}_{0\alpha}\right). \end{equation}

We first note that this relative entropy is, by definition of $\bar {P}_{0\alpha }$, non-negative and only zero if $P_{0\alpha } = \bar {P}_{0\alpha }$ for all $\alpha$. Furthermore, since $\bar {P}_{0\alpha }$ is determined by $P_{0\alpha }$ and $\bar {S}$ by $\bar {P}_{0\alpha }$, we may calculate this relative entropy at each time, and take its time derivative. Taking the time derivative of $S$ is straightforward given the collision integral (4.15). The time derivative of $\bar {S}$ can be determined in terms of only $\psi _{\alpha }$ and $\partial f_{0\alpha }/ \partial t$, as the entropy has the form (3.30), the only difference being that the evolution of $f_{0\alpha }$ is now governed by the hyperkinetic collision integral, not the multi-waterbag one. Using (3.32), we may, therefore, write the evolution of the relative entropy as

(4.35)\begin{align} \frac{\mathrm{d}}{\mathrm{d}t}\left(\bar{S} - S \right) & = \sum_{\alpha}\frac{V}{\Delta\Gamma_{\alpha}}\int\,\mathrm{d}\boldsymbol{v} \,\left[\psi_{\alpha}\frac{\partial f_{0\alpha}}{\partial t} + \int\,\mathrm{d}\eta \,\left(1 + \ln P_{0\alpha} \right)\frac{\partial P_{0\alpha}}{\partial t}\right]\nonumber\\ & = \sum_{\alpha}\frac{V}{\Delta\Gamma_{\alpha}}\iint\,\mathrm{d}\boldsymbol{v} \,\,\mathrm{d}\eta \, \left(\eta\psi_{\alpha} + \ln P_{0\alpha}\right)\frac{\partial P_{0\alpha}}{\partial t}. \end{align}

Consequently, when the waterbag distribution function satisfies the closure (3.36), the relative entropy will be zero by definition of $\bar {P}_{0\alpha }$ and, by (4.35), its time derivative will also be zero, so the relative entropy will stay zero for all future times. Since the relative entropy is zero only for $P_{0\alpha } = \bar {P}_{0\alpha }$ this proves our second statement. Therefore, it is possible to construct initial conditions for which the multi-waterbag collision integral (3.19) correctly describes the evolution due to the hyperkinetic collision integral for all future times.

Given that the distribution function $P_{0\alpha }(\boldsymbol {v},\eta )$ remains on the closure if it begins exactly on the closure, it is natural to ask if the evolution of the distribution function is, in fact, forced towards the closure by the hyperkinetic collision integral. If this were true, it would make the multi-waterbag collision integral (3.19) a valid approximation for the hyperkinetic collision integral after an initial transient period. However, for this to be true, there would have to be two timescales hidden within the hyperkinetic collision operator: a shorter timescale on which the hyperkinetic distribution function approached the closure and a longer timescale over which the closure evolved according to the multi-waterbag collision integral. It is trivially clear that this separation of timescales exists in a one-dimensional system, since the one-dimensional hyperkinetic collision integral forbids the mean phase-space density to change. In this case, the only evolution is towards satisfying the closure (3.36). In a higher number of dimensions, it may therefore be useful to consider two types of collisions: those which do and do not alter the mean phase-space density $f_{0\alpha }(\boldsymbol {v})$. This artificial divide would then describe two processes by which the system raises its entropy: one by making $f_{0\alpha }(\boldsymbol {v})$ approach a Lynden-Bell equilibrium and the other by reordering waterbags without altering $f_{0\alpha }(\boldsymbol {v})$. The latter process is one in which entropy is maximised locally in phase space (i.e., for each $\boldsymbol {v}$). A useful analogy might be the conventional collisional dynamics of gases, which relaxes quickly, at the collision rate, to a local Maxwellian, and slowly, at the diffusion rate, to the global one.

5.1. Balescu–Lenard collision integral

In reality, a plasma is not a phase fluid but a collection of $N$ particles. The true, exact phase-space density of these particles is the Klimontovich distribution (see, e.g., Klimontovich Reference Klimontovich1967):

(5.1)\begin{equation} f(\boldsymbol{r},\boldsymbol{v}) = \sum_{i=1}^{N}\delta(\boldsymbol{r}-\boldsymbol{r}_{i})\delta(\boldsymbol{v}-\boldsymbol{v}_{i}), \end{equation}

where $\boldsymbol {r}_{i}$ and $\boldsymbol {v}_{i}$ are the particles’ instantaneous positions and velocities, respectively. Since each particle thus occupies precisely zero phase volume, it might seem as though phase-volume conservation were a meaningless idea. The reason one can talk about phase-volume conservation at all is that one assumes that particles that are neighbours in phase space move in a similar way, implying that replacing the Klimontovich distribution with a smoothed phase-space density is a reasonable approximation. It is then the phase volumes associated with this smoothed function that are conserved. This is closely related to the microgranulation ansatz (2.32), which posits that, within a phase volume $\Delta \Gamma _{\alpha }$, the fluctuations of the phase-space density are correlated, implying that particles are moving collectively. A ‘true collision’ occurs when a single particle is not correlated at all with its neighbouring particles. This can be accommodated within the microgranulation ansatz by assuming that each particle is its own waterbag. Mathematically this is just the single-waterbag model of § 3.1 in the limit where the correlation volume $\Delta \Gamma _{\alpha }$ is made smaller that the inter-particle separation in phase space. Then, by assumption, there is only one particle in the correlation volume, so

(5.2)\begin{equation} \Delta\Gamma_{\alpha}\eta_{\alpha} = 1. \end{equation}

The smoothed distribution function is given by the average occurrence of single particles, hence $f_{0\alpha }(\boldsymbol {v}) \ll \eta _{\alpha }$. Under these assumptions, the Kadomtsev–Pogutse collision integral (3.5) becomes the Balescu–Lenard collision integral:

(5.3)\begin{gather} \frac{\partial f_{0\alpha}}{\partial t} = \sum_{\alpha''}\frac{16{\rm \pi}^{3}q_{\alpha}^{2}q_{\alpha''}^{2}}{m_{\alpha}V}\frac{\partial }{\partial \boldsymbol{v}}\boldsymbol{\cdot}\sum_{\boldsymbol{k}}\frac{\boldsymbol{k}\boldsymbol{k}}{k^{4}}\boldsymbol{\cdot}\int\,\mathrm{d}\boldsymbol{v}'' \,\frac{\delta(\boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{v}-\boldsymbol{v}''))}{\left|\epsilon_{\boldsymbol{k},\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}} \right|^{2}}\nonumber\\ \left[\frac{1}{m_{\alpha}}f_{0\alpha''}(\boldsymbol{v}'')\frac{\partial f_{0\alpha}}{\partial \boldsymbol{v}} - \frac{1}{m_{\alpha''}}f_{0\alpha}(\boldsymbol{v})\frac{\partial f_{0\alpha''}}{\partial \boldsymbol{v}''}\right]. \end{gather}

This reduction allows one to make a useful comparison between the ‘true’ collision rate, which is the typical collision rate associated with the Balescu–Lenard integral (5.3), and the effective collision rate associated with the hyperkinetic ‘collisionless collision integral’ (4.15). Before proceeding to do that, we observe in passing that the formalism developed in §§ 2 and 3 also allows one to derive very efficiently true collision operators for quantum plasmas consisting of fermionic and bosonic species. This is done in Appendix A.

5.2. Effective collision rates

To keep this discussion as transparent as possible, let us consider only the like-particle effective collision rates. Consider the generic quasilinear collision integral (2.35), with $\alpha '' = \alpha$ and

(5.4)\begin{gather} \langle g_{\alpha}^{2} \rangle(\boldsymbol{v}) = \int \,\mathrm{d}\eta \,\left[\eta - f_{0\alpha}(\boldsymbol{v}) \right]^{2} P_{0\alpha}(\boldsymbol{v},\eta), \quad f_{0\alpha}(\boldsymbol{v}) = \int\,\mathrm{d}\eta \,\eta P_{0\alpha}(\boldsymbol{v},\eta), \end{gather}

where $P_{0\alpha }(\boldsymbol {v},\eta )$ is evolved by (4.15). This integral should be compared to its Balescu–Lenard counterpart (5.3), also with $\alpha '' = \alpha$. It is immediately apparent that the difference in the rates of the true and effective collisions is

(5.5)\begin{equation} \frac{\nu_{\alpha\alpha}^{\mathrm{eff}}}{\nu_{\alpha\alpha}^{\mathrm{true}}} \sim \frac{\Delta\Gamma_{\alpha}\langle g_{\alpha}^{2} \rangle}{f_{0\alpha}} \equiv \Delta\Gamma_{\alpha}\eta_{\alpha}^{\mathrm{eff}}, \end{equation}

where $\eta _{\alpha }^{\mathrm {eff}}$ is the typical deviation of the phase-space density from its mean. The quantity on the right-hand side of (5.5) is the typical variation in the number of particles contained in the correlation volume $\Delta \Gamma _{\alpha }$. Since this collisionless theory is built upon the assumption that the correlation volume is sufficiently large for its mean phase-space density $\eta$ to be meaningfully specified, we have inherently assumed that the number of particles contained in a correlation volume is large. However, this does not tell us immediately about the typical deviation of this number from its mean. We therefore aim to compare $\eta _{\alpha }^{\mathrm {eff}}$ to the typical value of the phase-space density, $n_{\alpha }/v_{\mathrm {th}\alpha }^{3}$. To make this comparison quantitative, we estimate (5.5) close to a Lynden-Bell equilibrium. In this case, using (3.18) with ${\psi _{\alpha } = \beta \Delta \Gamma _{\alpha }\epsilon _{\alpha }}$, we find

(5.6)\begin{equation} \Delta\Gamma_{\alpha}\eta_{\alpha}^{\mathrm{eff}} ={-}\frac{1}{\beta}\frac{\partial \ln f_{0\alpha}}{\partial \epsilon_{\alpha}} \sim \frac{1}{\beta m_{\alpha} v_{\mathrm{th}\alpha}^{2}}, \end{equation}

where $v_{\mathrm {th}\alpha }$ is the typical velocity scale of $f_{0\mathrm {\alpha }}$. Since $\beta$ must be determined from the constraints of energy conservation (2.29) and phase-volume conservation (3.39), (5.6) provides an implicit expression for the $\eta _{\alpha }^{\mathrm {eff}}$ in terms of $E$ and $\rho _{\alpha }(\eta )$. For the purposes of order-of-magnitude estimates, the exact calculation of $\beta$ is unnecessary and it will suffice to consider two relevant limits. In the degenerate case, when the initial condition is very close to the minimum-energy state (the Gardner state; see § 3.3.1), $\beta$ will be very large and can be computed by a Sommerfeld-like expansion. In this case, naively ordering ${\eta _{\mathrm {\alpha }}^{\mathrm {eff}} \sim n_{\alpha }/v_{\mathrm {th}\alpha }^{3}}$ would be an overestimation. This is because deviations from the mean phase-space density require the system to have sufficient energy to allow two species of waterbag to intermingle (see discussion around § 3.3.1). Thankfully, as discussed in § 3.3.1, the system does not need to be energetically very far from its Gardner distribution to become, at least partially, non-degenerate. In the non-degenerate limit, which is by far the most common, waterbags of different phase-space density can freely intermingle, and so we can estimate $\eta _{\mathrm {\alpha }}^{\mathrm {eff}} \sim n_{\alpha }/v_{\mathrm {th}\alpha }^{3}$, implying $\Delta \Gamma _{\alpha }\eta _{\alpha }^{\mathrm {eff}}$ is on the order of the number of particles in a correlation volume, which we have assumed to be large.

Thus, the effective collision rate is generally much larger than the true collision rate. Note, however, that this fast collisionless relaxation must still be slower that the evolution of the fluctuations. The most straightforward estimate of the typical rate of the latter is the plasma frequency, $\omega _{\mathrm {pe}} = (4{\rm \pi} {\rm e}^{2}n_{\mathrm {e}}/m_{\mathrm {e}})^{1/2}$ (specialising to electrons for the purposes of this estimate). Thus, we require (and expect) the ordering

(5.7)\begin{equation} \nu_{\mathrm{ee}}^{\mathrm{true}} \sim \frac{\omega_{\mathrm{pe}}}{n_{\mathrm{e}}\lambda_{\mathrm{De}}^{3}} \ll \nu_{\mathrm{ee}}^{\mathrm{eff}} \ll \omega_{\mathrm{pe}}, \end{equation}

where $\lambda _{\mathrm {De}}\sim \omega _{\mathrm {pe}}v_{\mathrm {the}}$ is the Debye length. This places a constraint on the correlation volume:

(5.8)\begin{equation} 1 \ll \Delta\Gamma_{\mathrm{e}}\eta_{\mathrm{e}}^{\mathrm{eff}} \ll n_{\mathrm{e}}\lambda_{\mathrm{De}}^{3}. \end{equation}

If the plasma parameter $n_{\mathrm {e}}\lambda _{\mathrm {De}}^{3}$ is large (i.e., if the ideal-gas approximation applies), this constraint is not very stringent and still allows the effective collision frequency to be much larger than $\nu _{\mathrm {ee}}^{\mathrm {true}}$.

A more stringent constraint emerges if one works out the time $\tau _{\mathrm {c}}$ that it takes for exact phase-volume conservation to be broken. This timescale is substantially less well understood as it depends on the exact rate at which the fluctuating part $\delta f_{\alpha }$ of the exact phase-space density is altered irrevocably by collisions. A typical estimate (see, e.g., Su & Oberman Reference Su and Oberman1968) is

(5.9)\begin{equation} \tau_{\mathrm{c}}^{{-}1} \sim (\nu_{ee}^\mathrm{true})^{1/3}\omega_{\mathrm{pe}}^{2/3} \sim \frac{\omega_{\mathrm{pe}}}{(n_{\mathrm{e}}\lambda_{\mathrm{De}}^{3})^{1/3}}. \end{equation}

We may then argue that, for the collisionless relaxation to the Lynden-Bell equilibria to be of any importance, it must happen long before phase-volume conservation is broken: (5.8) is then revised to

(5.10)\begin{equation} (n_{\mathrm{e}}\lambda_{\mathrm{De}}^{3})^{2/3}\ll \Delta\Gamma_{\mathrm{e}}\eta_{\mathrm{e}}^{\mathrm{eff}} \ll n_{\mathrm{e}}\lambda_{\mathrm{De}}^{3}. \end{equation}

Depending on just how large the plasma parameter is, this could be a more difficult ordering to satisfy. However, if satisfied, it would make the collisionless relaxation rate far larger than the rate of relaxation due to ‘true’ collisions. To know just how much larger, we must be able to calculate $\Delta \Gamma _{\alpha }$ independently. Obviously this, along with a quantitative assessment of the validity of the microgranulation ansatz, requires a full theory of the two-point correlation function of the phase-space density. Without such a theory, certain revealing estimates can, nevertheless, be made.

5.3. Energy of fluctuations and the correlation volume

In all the above, the major shortcoming of the theory is the lack of clarity about the size of $\Delta \Gamma _{\alpha }$. In order to relate $\Delta \Gamma _{\alpha }$ to something measurable within a plasma, let us calculate the energy $E_{\varphi }$ stored in the electric-field fluctuations under the microgranulation ansatz (4.12). This energy is given by

(5.11)\begin{equation} E_{\varphi} = \int\,\mathrm{d}\boldsymbol{r} \,\frac{|\nabla\varphi|^{2}}{8{\rm \pi}} = \frac{1}{2}V\sum_{\alpha}q_{\alpha}\sum_{\boldsymbol{k}}\iint\,\mathrm{d}\boldsymbol{v} \,\,\mathrm{d}\eta \,\eta\mathrm{Re}\langle \varphi_{\boldsymbol{k}}^{*}P_{\boldsymbol{k}\alpha}(\boldsymbol{v},\eta) \rangle, \end{equation}

where we have integrated by parts and used Poisson's equation (4.6) in the second equality. Only the real part has survived because the summation in $\boldsymbol {k}$ is even. While one could again apply Poisson's equation to the remaining $\varphi _{\boldsymbol {k}}$ in (5.11) and then naively utilise the microgranulation ansatz (4.12), it is worth noting that the correlator in (5.11) is the real part of the same correlator the imaginary part of which appeared in (4.5). The real part of this correlator is equivalent to the hyperkinetic generalisation of the discarded term of (2.16). Recovering the hyperkinetic generalisation of this term, from (2.16) via (2.37), we get

(5.12)\begin{gather} \mathrm{Re}\langle \varphi_{\boldsymbol{k}}^{*}P_{\boldsymbol{k}\alpha}(\boldsymbol{v},\eta) \rangle = \mathrm{Re}\sum_{\alpha}\frac{4{\rm \pi} q_{\alpha'}}{k^{2}}\iint \frac{\mathrm{d}p \,\,\mathrm{d}p' \,}{(2{\rm \pi})^{2}}\frac{{\rm e}^{(p+p')t}}{\epsilon_{\boldsymbol{k}}(p)\epsilon_{\boldsymbol{k}}^{*}(p'^{*})} \nonumber\\ \times\iint\,\mathrm{d}\boldsymbol{v}' \,\,\mathrm{d}\eta' \,\eta' \frac{\langle g_{\boldsymbol{k}\alpha}(\boldsymbol{v},\eta)g_{\boldsymbol{k}\alpha'}^{*}(\boldsymbol{v}',\eta') \rangle}{(p+{\rm i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v})(p'-{\rm i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}')}, \end{gather}

with the understanding that, after the microgranulation ansatz, the real part of the other contributions will vanish. Carrying out the $p$ and $p'$ contour integration (see § 2.3 for details) and applying the microgranulation ansatz (4.12), we find

(5.13)\begin{equation} \mathrm{Re}\langle \varphi_{\boldsymbol{k}}^{*}P_{\boldsymbol{k}\alpha}(\boldsymbol{v},\eta) \rangle = \frac{4{\rm \pi} q_{\alpha}\Delta\Gamma_{\alpha}}{k^{2}V|\epsilon_{\boldsymbol{k},\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}}|^{2}}\int \,\mathrm{d}\eta' \,\eta' \langle g_{\boldsymbol{k}\alpha}(\boldsymbol{v},\eta)g_{\boldsymbol{k}\alpha}^{*}(\boldsymbol{v},\eta') \rangle. \end{equation}

Substituting (5.13) into (5.11) and using the formula (4.13) for the correlation function of $g_{\boldsymbol {k}\alpha }(\boldsymbol {v},\eta )$, we get

(5.14)\begin{align} {E_{\varphi}} & = \sum_{\alpha}\sum_{\boldsymbol{k}}\frac{2{\rm \pi} q_{\alpha}^{2}\Delta\Gamma_{\alpha}}{k^{2}}\int\,\mathrm{d}\boldsymbol{v} \,\frac{1}{{|\epsilon_{\boldsymbol{k},\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}}|^{2}}}\iint\,\mathrm{d}\eta \,\,\mathrm{d}\eta' \,\eta\eta' \left[\delta(\eta-\eta')P_{0\alpha}(\boldsymbol{v},\eta)\right.\nonumber\\ & \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad\left. \quad -\,P_{0\alpha}(\boldsymbol{v},\eta)P_{0\alpha}(\boldsymbol{v},\eta') \right] \nonumber\\ & = \sum_{\alpha}\sum_{\boldsymbol{k}}\frac{2{\rm \pi} q_{\alpha}^{2}\Delta\Gamma_{\alpha}}{k^{2}}\int\,\mathrm{d}\boldsymbol{v} \,\frac{\langle g_{\alpha}^{2} \rangle(\boldsymbol{v})}{|\epsilon_{\boldsymbol{k},\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}}|^{2}}, \end{align}

where $\langle g_{\alpha }^{2} \rangle (\boldsymbol {v})$ is given by (5.4). This formula describes the energy in the fluctuating electric field for a plasma obeying the microgranulation ansatz.

In a Lynden-Bell equilibrium, (5.14) can be further simplified by calculating $\Delta \Gamma _{\alpha }\langle g_{\alpha }^{2} \rangle$ via (5.5) and (5.6). Approximating also ${|\epsilon _{\boldsymbol {k},\boldsymbol {k}\boldsymbol {\cdot }\boldsymbol {v}}|^{2}\approx 1}$, we have

(5.15)\begin{equation} E_{\varphi} \approx \sum_{\alpha}\sum_{\boldsymbol{k}}\frac{2{\rm \pi} q_{\alpha}^{2}}{k^{2}}\int \,\mathrm{d}\boldsymbol{v} \,\left(-\frac{1}{\beta}\frac{\partial f_{0\alpha}}{\partial \epsilon_{\alpha}} \right) = \sum_{\alpha}\frac{4Vq_{\alpha}^{2}k_{\mathrm{max}}}{\beta m_{\alpha}}\int_{0}^{\infty} \,\mathrm{d}v \,f_{0\alpha}(v), \end{equation}

where $k_{\mathrm {max}}$ is the ultraviolet cutoff for the wavenumber integral (in three dimensions). Using (5.6) and making a rough estimate of everything, we have

(5.16)\begin{equation} E_{\varphi} \sim \sum_{\alpha}\Delta\Gamma_{\alpha}\eta_{\alpha}^{\mathrm{eff}}q_{\alpha}^{2}n_{\alpha}k_{\mathrm{max}}V. \end{equation}

Let us compare these estimates with the energy stored in electric fluctuations associated with true, Coulomb collisions (i.e., with discrete particle noise). The latter can be calculated in a manner analogous to the above but using the Balescu–Lenard collision operator (5.3). To do this, we set

(5.17)\begin{equation} P_{0\alpha}(\boldsymbol{v},\eta) = \left[1-\frac{f_{0\alpha}(\boldsymbol{v})}{\eta_{\alpha}}\right]\delta(\eta) + \frac{f_{0\alpha}(\boldsymbol{v})}{\eta_{\alpha}}\delta(\eta - \eta_{\alpha}), \end{equation}

where $\eta _{\alpha } = \Delta \Gamma _{\alpha }^{-1}$ according to (5.2). Inserting (5.17) into (5.14) and collecting only the lowest-order terms in the limit of $f_{0\alpha }/\eta _{\alpha } \to 0$ gives

(5.18)\begin{equation} E_{\varphi,\mathrm{BL}} = \sum_{\alpha}\sum_{\boldsymbol{k}}\frac{2{\rm \pi} q_{\alpha}^{2}}{k^{2}}\int\,\mathrm{d}\boldsymbol{v} \,\frac{f_{0\alpha}(\boldsymbol{v})}{|\epsilon_{\boldsymbol{k},\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}}|^{2}}\approx \sum_{\alpha}4V q_{\alpha}^{2}k_{\mathrm{max}}n_{\alpha}, \end{equation}

where the last estimate has been obtained in the same manner as (5.15). From a direct comparison of (5.14) and (5.16) with (5.18), it is clear that, both $\boldsymbol {k}$ by $\boldsymbol {k}$ and overall, a plasma obeying the microgranulation ansatz has more energy stored in the electric field than a collisional plasma, by a factor of $\Delta \Gamma _{\alpha }\eta _{\alpha }^{\mathrm {eff}} \gg 1$, the typical number of particles in a correlation volume.Footnote ²

Note that, despite the fluctuation energy being large compared to particle noise, it will still be small compared to the kinetic energy of the particles

(5.19)\begin{equation} K = V\sum_{\alpha}\int\,\mathrm{d}\boldsymbol{v} \,\frac{m_{\alpha}|\boldsymbol{v}|^{2}}{2}f_{0\alpha} \sim V\sum_{\alpha}m_{\alpha}v_{\mathrm{th}\alpha}^{2}n_{\alpha}. \end{equation}

Therefore, using (5.16), we find

(5.20)\begin{equation} \frac{E_{\varphi}}{K} \sim \frac{\sum_{\alpha}\Delta\Gamma_{\alpha}^{\mathrm{eff}}\eta_{\alpha}^{\mathrm{eff}}q_{\alpha}^{2}n_{\alpha}k_{\mathrm{max}}}{\sum_{\alpha}m_{\alpha}v_{\mathrm{th}\alpha}^{2}n_{\alpha}} \sim \frac{\Delta\Gamma_{\mathrm{e}}\eta_{\mathrm{e}}^{\mathrm{eff}}}{n_{\mathrm{e}}\lambda_{\mathrm{De}}^{3}}k_{\mathrm{max}}\lambda_{\mathrm{De}} \ll 1, \end{equation}

ignoring at the last step any potential disparities between contributions from different species. Taking a cue from our discussion in § 5.2, we conclude that this ratio will be small in all conceivable cases of interest (at least as long as this theory is valid).

Since the distributions evolved by the quasilinear collision integrals are assumed to be linearly stable (see § 2.3), their kinetic energy does not change [see (2.30)] and, therefore, neither can the fluctuation energy $E_{\varphi }$ change. It is then useful to think of $E_{\varphi }$ as a feature of the system that tells us about the size of the correlation volume $\Delta \Gamma _{\alpha }$ and of the applicability of the collisionless relaxation, via the estimate (5.16).

This calculation does not determine $\Delta \Gamma _{\alpha }$, however, together with the discussion of the collision timescales, it sets feasible limits on what the allowed values of $\Delta \Gamma _{\alpha }$ can be. Crucially, we see that a larger value of $\Delta \Gamma _{\alpha }$ corresponds to a larger energy in the fluctuating electric field. This further suggests that the validity of the Balescu–Lenard collision integral should be called into question when the fluctuations in the electric field are anomalously large compared to (5.18). The Balescu–Lenard collision integral is designed to describe slow relaxation mediated by discrete particle noise. The hyperkinetic collision integral, on the other hand, describes a faster relaxation mediated by correlated volumes of phase space.

5.4. Caveats on the existence of a Lynden-Bell plasma

Despite promising universality, previous attempts at the numerical demonstration of systems reaching their Lynden-Bell equilibria have not been met with universal success. Instead, only certain initial conditions will cause a system to reach its Lynden-Bell equilibrium (see, e.g., Arad & Johansson Reference Arad and Johansson2005; Levin et al. Reference Levin, Pakter, Rizzato, Teles and Benetti2014 and references therein), which naturally forces one to consider the validity of this theory, in particular, its chief assumption: the microgranulation ansatz. From §§ 5.1–5.3, we are now equipped with a better understanding of how the microgranulation ansatz has affected not only the collision integral, but also the fluctuations of the electric field. It is therefore now prudent, in spite of the comfort provided by the estimates (5.10) and (5.16), to question what has been lost and, therefore, what the validity of this closure is.

The defining, but also the most fragile, feature of the collisionless collision integrals derived above is phase-volume conservation by the evolution of the exact phase-space density $f_{\alpha }$. We argued in § 5.2 that, to be in any way important, our relaxation had to be faster than phase-volume conservation could be broken. While estimates (5.8) and (5.9) serve as guesses for how fast phase-volume conservation is broken, there is numerical evidence (Zhdankin Reference Zhdankin2021) that in turbulent systems Casimir invariants can be broken on timescales that are independent of the true collision rate – analogously perhaps to the way in which ‘dissipative anomalies’ arise generically in turbulent environments (cf. Schekochihin et al. Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009; Eyink Reference Eyink2018). If this is the case, then one would not be justified in writing a collisionless Vlasov equation (2.1) for $f_{\alpha }$ or (4.2) for $P_{\alpha }$, rendering the subsequent calculations formally invalid (cf. Beraldo e Silva et al. Reference Beraldo e Silva, de Siqueira Pedra, Sodré, Perico and Lima2017). In such systems, it may nevertheless be a meaningful question to ask whether the waterbag density $\rho (\eta )$ (equivalently the Casimir invariants) will evolve towards anything universal. If so, then despite its non-conservation, waterbag content could still be a relevant constraint.

Should phase-volume conservation survive, we must further ask if it is reasonable to expect that the correlation volume $\Delta \Gamma _{\alpha }$ should be a constant independent of time. Naturally one might assume that, as fluctuations travel through phase space, they gradually break up and form smaller and smaller structures. Another factor that indicates that $\Delta \Gamma _{\alpha }$ should be dependent on time is (5.14) itself. Since the kinetic energy of the mean distribution function is also constant, then to conserve the total energy, the electric-field energy must also be constant. By (5.14), this implies that $\Delta \Gamma _{\alpha }$ must vary in concert with the integral of $\langle g_{\alpha }^{2} \rangle$ over velocities. Alternatively one could take the view that the mean energy of the distribution function is only approximately constant, justified by the smallness of the electric-field energy compared to the kinetic energy as computed in (5.20).

Only one of our results requires that $\Delta \Gamma _{\alpha }$ not be a function of time: the H-theorem. If $\Delta \Gamma _{\alpha }$ evolves with time, then the H-theorem is broken because our entropy (3.41) has a prefactor proportional to $\Delta \Gamma _{\alpha }$. This could be remedied if $\Delta \Gamma _{\alpha }$ were independent of the species, because it would be an overall multiplicative prefactor that the entropy need not include. This would result in a working, but weaker, H-theorem under which any state could be a steady state if $\Delta \Gamma _{\alpha }$ decayed sufficiently fast to halt its evolution. Such relaxation has been called ‘incomplete violent relaxation’, indicating that the system tried to reach its Lynden-Bell equilibrium, but stalled before the relaxation could be completed (Chavanis Reference Chavanis2006).

Thus, for the phase-volume conservation and the correlation volume to be meaningful features of any complete theory, they must earn their place in it. Since the microgranulation ansatz grants them a privileged position without question, it cannot be trusted without question. Nevertheless, it provides an insight into what interesting effects could be contained in a theory where these features are valid and meaningful. We give one such interesting example in the next section.

6.1. Preview of strange relaxation

To study the interaction between species we will write our collision operator (4.15) as

(6.3)\begin{equation} \frac{\partial P_{\alpha}}{\partial t} = \sum_{\alpha'}C_{\alpha\alpha'}[P_{\alpha},P_{\alpha'}], \end{equation}

in terms of the interspecies collision operator $C_{\alpha \alpha '}$, which can be easily read from elements of the species sum in (4.15). For compactness of notation, we shall henceforth drop the zeros from the subscripts of the mean distribution function. For electrons interacting with ions, the dominant effect is the diffusion-like first term of (4.15). Comparing this to the true collision operator (5.3), we see that (aside from acting on $P_{\mathrm {e}}$ instead of $f_{\mathrm {e}}$) the new, ‘collisionless’ feature in this diffusion term is that $f_{\mathrm {i}}$ is replaced by $\Delta \Gamma _{\mathrm {i}}\langle g_{\mathrm {i}}^{2} \rangle$ with $\langle g_{\mathrm {i}}^{2} \rangle$ defined by (4.32). Therefore, any effects on the electrons relating to $f_{\mathrm {i}}$ in the standard collisional theory will now be replaced by the analogous effects relating to $\Delta \Gamma _{\mathrm {i}}\langle g_{\mathrm {i}}^{2} \rangle$ (see figure 4 for an illustration of this in a single-waterbag distribution). In particular, instead of being isotropised by the ions of density $n_{\mathrm {i}}$ and dragged towards the ion velocity $\boldsymbol {u}_{\mathrm {i}}$, the electrons will see an ‘anomalous density’ $n_{\mathrm {i}}^{\mathrm {a}}$ and find themselves dragged towards an ‘anomalous velocity’ $\boldsymbol {u}_{\mathrm {i}}^{\mathrm {a}}$, given by

(6.4)\begin{gather} n_{\mathrm{i}}^{\mathrm{a}} = \Delta\Gamma_{\mathrm{i}}\int\,\mathrm{d}\boldsymbol{v} \,'\langle g_{\mathrm{i}}^{2} \rangle(\boldsymbol{v}'), \quad \boldsymbol{u}_{\mathrm{i}}^{\mathrm{a}} = \frac{\int\,\mathrm{d}\boldsymbol{v}' \,\boldsymbol{v}'\langle g_{\mathrm{i}}^{2} \rangle(\boldsymbol{v}')}{\int\,\mathrm{d}\boldsymbol{v}' \, \langle g_{\mathrm{i}}^{2} \rangle(\boldsymbol{v}')}. \end{gather}

Since ion relaxation is slow compared to the electron one, the ion distribution need not be isotropic, so, in general, $\boldsymbol {u}_{\mathrm {i}}^{\mathrm {a}} \neq 0$. Furthermore, there is no need for $\boldsymbol {u}_{\mathrm {i}}^{\mathrm {a}}$ to point in the same direction as the ions’ mean velocity $\boldsymbol {u}_{\mathrm {i}}$, and indeed it is even possible that $\boldsymbol {u}_{\mathrm {i}} = 0$ while $\boldsymbol {u}_{\mathrm {i}}^{\mathrm {a}} \neq 0$. This means that collisionless relaxation can lead to spontaneous generation of current – and, therefore, of magnetic field – from ion anisotropies (e.g., from an ion heat flux).

Figure 4. A cartoon of a possible anisotropic ion distribution function that possesses an anomalous velocity $\boldsymbol {u}_{\mathrm {i}}^{\mathrm {a}}$. The solid line represents the mean ion distribution function. If the exact ion phase-space density is viewed as a single waterbag [as in the Kadomtsev–Pogutse collision integral (3.5)], then the ‘effective distribution’ of ions with which the electrons will interact is shown by the dashed black line. The mean ion velocity is zero, but the anomalous velocity is not, as shown by the dashed red line. Note that the height of the effective distribution would be rescaled if $\eta$ and $\Delta \Gamma _{\mathrm {i}}$ took different values, but this would only affect the rate of relaxation to the anomalous velocity, not the anomalous velocity itself.

The physical significance of this becomes especially clear for a single-waterbag example shown in figure 4. The phase-space exclusion effect makes the ions behave as though they were fermions. Therefore, the probability for an electron to have a ‘collision’ with an ion whose velocity is $\boldsymbol {v}$ is not just the probability $f_{\mathrm {i}}(\boldsymbol {v})/\eta$ that an ion can be found at that velocity but the probability that an ion can be found at that velocity $\boldsymbol {v}$ and that the velocity $\boldsymbol {v} + \Delta \boldsymbol {v}$ into which that ion will be scattered by the interaction is not already occupied. Since an ion is deflected by a very small amount in a collision with an electron ($\Delta \boldsymbol {v} \ll \boldsymbol {v}$), this probability is approximately proportional to $f_{\mathrm {i}}(\boldsymbol {v})[\eta - f_{\mathrm {i}}(\boldsymbol {v})]$, which is precisely $\langle g_{\mathrm {i}}^{2} \rangle (\boldsymbol {v})$ from which the anomalous ion velocity and density are defined in (6.4). This means that densely packed portions of the ion phase space become effectively invisible to the electrons. Thus, the mean velocity towards which the electrons are dragged is not the true ion mean velocity but the mean velocity of those ions that move in the less densely occupied regions of the phase space.

After an initial period of this ‘strange relaxation’ the electrons will converge to their own Lynden-Bell equilibrium moving at the anomalous ion velocity $\boldsymbol {u}_{\mathrm {i}}^{\mathrm {a}}$. The ions, however, will not thus far have relaxed significantly except to ensure the conservation of total momentum (and thus alter the mean velocity by a mass-ratio-small amount ${\sim -\boldsymbol {u}_{\mathrm {i}}^{\mathrm {a}}m_{\mathrm {e}}n_{\mathrm {e}}/m_{\mathrm {i}}n_{\mathrm {i}}}$). Once the ion–ion relaxation timescale is reached, the ions will begin to erase any anisotropy in their distribution function as they proceed towards their own Lynden-Bell equilibrium. As their anisotropy vanishes, the anomalous mean ion velocity will tend to the true mean ion velocity, $\boldsymbol {u}_{\mathrm {i}}^{\mathrm {a}}\to \boldsymbol {u}_{\mathrm {i}}$, and the electrons’ mean velocity will doggedly follow it, relaxing towards two Lynden-Bell equilibria of equal velocities but distinct temperatures. Finally, on the longest (ion–electron) interaction timescale, the distributions will equalise their temperatures, reaching the overall maximum-entropy state and completing the relaxation.

This is the physics of the strange relaxation process. In the remainder of this section, we demonstrate formally that this is indeed what happens by carrying out the mass-ratio expansion of the electron–ion (§ 6.3.2) and the ion–electron (§ 6.4) collision operators. This calculation follows the standard path well trodden for true collisions (specifically, as presented, e.g., in Parra Reference Parra2019), but with a few important adjustments.

6.2. Landau form of the hyperkinetic collision operator

We will first make a further simplification by reducing our hyperkinetic collision operator (4.15) to the so-called ‘Landau form’. Doing so amounts to finding an approximate expression for the $\boldsymbol {k}$ sum

(6.5)\begin{equation} \frac{1}{V}\sum_{\boldsymbol{k}}\frac{\boldsymbol{k}\boldsymbol{k}}{k^{4}}\frac{\delta\left(\boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{v}-\boldsymbol{v}') \right)}{|\epsilon_{\boldsymbol{k},\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}}|^{2}}, \end{equation}

which is complicated by the presence of the dielectric function

(6.6)\begin{equation} \epsilon_{\boldsymbol{k},\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}} = 1 + \sum_{\alpha'}\frac{4{\rm \pi} q_{\alpha'}^{2}}{m_{\alpha'}k^{2}}\iint \,\mathrm{d}\boldsymbol{v}' \,\,\mathrm{d}\eta' \,\frac{\eta'}{\boldsymbol{k}\boldsymbol{\cdot}\left(\boldsymbol{v}-\boldsymbol{v}' \right)}\boldsymbol{k}\boldsymbol{\cdot}\frac{\partial P_{\alpha'}}{\partial \boldsymbol{v}}. \end{equation}

Here, the $\boldsymbol {v}'$ integral is taken along the Landau contour. To simplify (6.5), the important feature to note is that, for length scales shorter than the Debye length and in the absence of instabilities, the second term in (6.6) will be small and $\epsilon _{\boldsymbol {k},\boldsymbol {k}\boldsymbol {\cdot }\boldsymbol {v}}$ can be approximated by unity. For scales significantly longer than the Debye length, $\epsilon _{\boldsymbol {k},\boldsymbol {k}\boldsymbol {\cdot }\boldsymbol {v}} \sim (\lambda _{\mathrm {De}}k)^{-2}$, which will make their contribution to the sum in (6.5) small in $(\lambda _{\mathrm {De}}k)^{4}$. Therefore, following Landau (Reference Landau1936), we truncate the $\boldsymbol {k}$ sum at $(k\lambda _{\mathrm {De}})\gtrsim 1$, and approximate the dielectric function by unity. Denoting $\boldsymbol {w} = \boldsymbol {v}-\boldsymbol {v}'$ and integrating in cylindrical coordinates such that $\boldsymbol {k} = k_{\parallel }\hat {\boldsymbol {w}} + \boldsymbol {k}_{\perp }$, we get

(6.7)\begin{equation} \frac{1}{V}\sum_{\boldsymbol{k}}\frac{\boldsymbol{k}\boldsymbol{k}}{k^{4}}\frac{\delta(\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{w})}{|\epsilon_{\boldsymbol{k},\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}}|^{2}} \approx \iiint\frac{\mathrm{d}\boldsymbol{k} \,}{\left(2{\rm \pi} \right)^{3}}\frac{\boldsymbol{k}_{{\perp}}\boldsymbol{k}_{{\perp}}}{k_{{\perp}}^{4}}\delta(k_{{\parallel}}w) = \frac{1}{8{\rm \pi}^{2}}\frac{1}{w}\left(\mathbb{I} -\frac{\boldsymbol{w}\boldsymbol{w}}{w^{2}} \right)\ln\left(\frac{k_{\mathrm{max}}}{k_{\mathrm{min}}}\right). \end{equation}

Here, we have truncated the $k_{\perp }$ integral, not only at $k_{\mathrm {min}}^{-1} \sim \lambda _{\mathrm {De}}$ as discussed above, but also at scales shorter than $k_{\mathrm {max}}^{-1}$. In the standard collisional theory, this corresponds to cutting off the integral at the distance of closest approach of two particles and gives rise to the Coulomb logarithm. In our treatment, such a high-$k$ cutoff should instead be associated with the physical length scale of the correlations, which may be thought of as the distance of closest approach for two waterbags. Henceforth, we will denote this generalised Coulomb-logarithm-like quantity by $\varLambda _{\alpha \alpha '}$.

With the approximation (6.7), the interspecies collision operator in (4.15) becomes

(6.8)\begin{align} & C_{\alpha\alpha'}[P_{\alpha},P_{\alpha'}] = \frac{\gamma_{\alpha\alpha'}}{m_{\alpha}}\frac{\partial }{\partial \boldsymbol{v}}\boldsymbol{\cdot}\int\frac{\mathrm{d}\boldsymbol{v} \,'}{w}\left( \mathbb{I} -\frac{\boldsymbol{w}\boldsymbol{w}}{w^{2}}\right)\boldsymbol{\cdot}\int\,\mathrm{d}\eta' \,\eta'\nonumber\\ & \qquad \left\{ \frac{\Delta\Gamma_{\alpha'}}{m_{\alpha}}\left[\eta' - f_{\alpha'}(\boldsymbol{v}') \right]P_{\alpha'}(\boldsymbol{v}',\eta')\left.\frac{\partial P_{\alpha}}{\partial \boldsymbol{v}}\right|_{\eta} - \frac{\Delta\Gamma_{\alpha}}{m_{\alpha'}}\left[\eta - f_{\alpha}(\boldsymbol{v}) \right]P_{\alpha}(\boldsymbol{v},\eta)\left.\frac{\partial P_{\alpha'}}{\partial \boldsymbol{v}'}\right|_{\eta'}\right\}, \end{align}

where $\gamma _{\alpha \alpha '} = 2{\rm \pi} q_{\alpha }^{2}q_{\alpha '}^{2}\varLambda _{\alpha \alpha '}$.

6.3. Electron–ion relaxation

6.3.1. Zeroth order: isotropisation of electrons

To lowest order in $\sqrt {m_{\mathrm {e}}/m_{\mathrm {i}}}$, the second term in brackets in (6.8) vanishes entirely and the velocity difference between the electrons and ions, $\boldsymbol {w}$, becomes approximately equal to the electron velocity $\boldsymbol {v}$. This leaves the collision operator in the immensely simple form

(6.9)\begin{equation} C_{\mathrm{ei}} = \frac{\gamma_{\mathrm{ei}}n_{\mathrm{i}}^{\mathrm{a}}}{m_{\mathrm{e}}^{2}}\frac{\partial }{\partial \boldsymbol{v}}\boldsymbol{\cdot}\left[\frac{1}{v}\left(\mathbb{I} - \frac{\boldsymbol{v}\boldsymbol{v}}{v^{2}}\right)\boldsymbol{\cdot} \left. \frac{\partial P_{\mathrm{e}}}{\partial \boldsymbol{v}}\right|_{\eta} \right], \end{equation}

where $n_{\mathrm {i}}^{\mathrm {a}}$ is the anomalous ion density given by the first expression in (6.4). This operator is a pitch-angle-scattering (Lorentz) operator (see, e.g., Helander & Sigmar Reference Helander and Sigmar2005), which causes relaxation on timescales comparable to that of electron–electron interactions. Its effect is to isotropise $P_{\mathrm {e}}$. Physically this is a consequence of the ions’ high mass. Like ping-pong balls bouncing off bowling balls, the electrons bounce off the ions without exchanging any energy and the only effect is to isotropise their distribution. However, to this lowest order, we have only retained terms ordered with the electron thermal velocity, losing any speeds ordered with the ion velocity. To retain velocities of that size, viz., the electrons’ mean velocity, we must go to next order.

6.3.2. First order: anomalous drag

To next order, determining the collision integral becomes moderately less trivial. First, we will expand the electron hyperkinetic distribution function in mass ratio and assume that its lowest-order part, denoted by $P_{\mathrm {e}}^{\mathrm {I}}(v,\eta )$, is isotropic. We further assume that any deviations from isotropy are due to velocities ordered with the ion thermal velocity, so the anisotropic correction, denoted $P_{\mathrm {e}}^{\mathrm {A}}(\boldsymbol {v},\eta )$, will be smaller than $P_{\mathrm {e}}^{\mathrm {I}}(v,\eta )$ by a factor of $\sqrt {m_{\mathrm {e}}/m_{\mathrm {i}}}$. As well as this, since the $\boldsymbol {v}'$ integral in (6.8) is over ion velocities, $\boldsymbol {v}'$ will be small in $\sqrt {m_{\mathrm {e}}/m_{\mathrm {i}}}$ compared to $\boldsymbol {v}$. Therefore, we can now expand the $\boldsymbol {w}$-dependent tensor in (6.8) to first order in $\sqrt {m_{\mathrm {e}}/m_{\mathrm {i}}}$, giving

(6.10)\begin{align} \frac{1}{w}\left(\mathbb{I} - \frac{\boldsymbol{w}\boldsymbol{w}}{w^{2}} \right) & \approx \frac{1}{v}\left( \mathbb{I} - \frac{\boldsymbol{v}\boldsymbol{v}}{v^{2}} \right) - \boldsymbol{v}'\boldsymbol{\cdot}\frac{\partial }{\partial \boldsymbol{v}}\left[\frac{1}{v}\left(\mathbb{I} - \frac{\boldsymbol{v}\boldsymbol{v}}{v^{2}} \right) \right] \nonumber\\ & = \frac{1}{v}\left(1 + \frac{3\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{v}'}{v^{2}} \right)\left(\mathbb{I} - \frac{\boldsymbol{v}\boldsymbol{v}}{v^{2}} \right) + \frac{\boldsymbol{v}'\boldsymbol{v} + \boldsymbol{v}\boldsymbol{v}' - 2(\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{v}')\mathbb{I}}{v^{3}}. \end{align}

Since the anisotropic part of the hyperkinetic distribution function is already small, we only need the lowest-order contribution from (6.10) to act on $P_{\mathrm {e}}^{\mathrm {A}}(\boldsymbol {v},\eta )$. For the isotropic part $P_{\mathrm {e}}^{\mathrm {I}}(v,\eta)$, the lowest-order contribution vanishes (as it must, since any isotropic distribution is a solution of the lowest-order problem) and we must keep the next-order terms. These terms are

(6.11)\begin{align} \frac{1}{w}\left(\mathbb{I} - \frac{\boldsymbol{w}\boldsymbol{w}}{w^{2}} \right)\boldsymbol{\cdot} \frac{\partial P_{\mathrm{e}}^{\mathrm{I}}}{\partial \boldsymbol{v}} & \approx \frac{\boldsymbol{v}'\boldsymbol{v} + \boldsymbol{v}\boldsymbol{v}' - 2(\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{v}')\mathbb{I}}{v^{3}}\boldsymbol{\cdot} \frac{\boldsymbol{v}}{v}\frac{\partial P_{\mathrm{e}}^{\mathrm{I}}}{\partial v} \nonumber\\ & = \frac{1}{v}\left(\mathbb{I} -\frac{\boldsymbol{v}\boldsymbol{v}}{v^{2}} \right)\boldsymbol{\cdot}\frac{\boldsymbol{v}'}{v}\frac{\partial P_{\mathrm{e}}^{\mathrm{I}}}{\partial v} = \frac{1}{v}\left(\mathbb{I} - \frac{\boldsymbol{v}\boldsymbol{v}}{v^{2}} \right)\boldsymbol{\cdot} \frac{\partial }{\partial \boldsymbol{v}}\left(\frac{\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{v}'}{v}\frac{\partial P_{\mathrm{e}}^{\mathrm{I}}}{\partial v} \right). \end{align}

In the final equality, we made use of the isotropy of $P_{\mathrm {e}}^{\mathrm {I}}$ in order to insert an additional derivative into the product. The effect of this is to cast this term in the form of a pitch-angle-scattering operator. The resulting collision operator is

(6.12)\begin{equation} C_{\mathrm{ei}}[P_{\mathrm{e}},P_{\mathrm{i}}] = \frac{\gamma_{\mathrm{ei}}}{m_{\mathrm{e}}^{2}}n_{\mathrm{i}}^{\mathrm{a}}\frac{\partial }{\partial \boldsymbol{v}}\boldsymbol{\cdot} \left\{\frac{1}{v}\left(\mathbb{I} -\frac{\boldsymbol{v}\boldsymbol{v}}{v^{2}} \right)\boldsymbol{\cdot} \frac{\partial }{\partial \boldsymbol{v}}\left[P_{\mathrm{e}}^{\mathrm{A}}(\boldsymbol{v},\eta) + \frac{\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{u}_{\mathrm{i}}^{\mathrm{a}}}{v}\frac{\partial P_{\mathrm{e}}^{\mathrm{I}}}{\partial v} \right]\right\}, \end{equation}

where $n_{\mathrm {i}}^{\mathrm {a}}$ and $\boldsymbol {u}_{\mathrm {i}}^{\mathrm {a}}$ are given by (6.4).

Despite not knowing the exact evolution due to (6.12), we may use the fact that it is a pitch-angle-scattering operator to read off the steady state. This will occur when the expression in the square brackets is isotropic. Without loss of generality, we may define $P_{\mathrm {e}}^{\mathrm {A}}$ to have zero spherical average (isotropic part), in which case the only way for the term in square brackets in (6.12) to be isotropic is for it to be identically zero. Therefore the fixed point of the collision operator (6.12) is

(6.13)\begin{equation} P_{\mathrm{e}}(\boldsymbol{v},\eta) = P_{\mathrm{e}}^{\mathrm{I}}(v,\eta) + P_{\mathrm{e}}^{\mathrm{A}}(\boldsymbol{v},\eta) = P_{\mathrm{e}}^{\mathrm{I}}(v,\eta) - \frac{\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{u}_{\mathrm{i}}^{\mathrm{a}}}{v}\frac{\partial P_{\mathrm{e}}^{\mathrm{I}}}{\partial v} \approx P_{\mathrm{e}}^{\mathrm{I}}(|\boldsymbol{v}-\boldsymbol{u}_{\mathrm{i}}^{\mathrm{a}}|,\eta), \end{equation}

the last equality holding with $O(m_{\mathrm {e}}/m_{\mathrm {i}})$ precision.

Thus, as prophesied in § 6.1, the electrons become isotropic around the anomalous ion velocity $\boldsymbol {u}_{\mathrm {i}}^{\mathrm {a}}$ rather than the true mean ion velocity. Dynamically, this manifests itself as a drag on the electrons with the rate of change of the electron momentum $m_{\mathrm {e}}n_{\mathrm {e}}\boldsymbol {u}_{\mathrm {e}}$ being

(6.14)\begin{align} m_{\mathrm{e}}n_{\mathrm{e}}\frac{\mathrm{d}\boldsymbol{u}_{\mathrm{e}}}{\mathrm{d}t}& = \frac{\mathrm{d}}{\mathrm{d}t}\int \,\mathrm{d}\boldsymbol{v} \,m_{\mathrm{e}}\boldsymbol{v}\int \,\mathrm{d}\eta \,\eta P_{\mathrm{e}}(\boldsymbol{v},\eta) = \int\,\mathrm{d}\boldsymbol{v} \,m_{\mathrm{e}}\boldsymbol{v}\int \,\mathrm{d}\eta \,\eta C_{\mathrm{ei}}[P_{\mathrm{e}},P_{\mathrm{i}}] \nonumber\\ & ={-}\frac{\gamma_{\mathrm{ei}}n_{\mathrm{i}}^{\mathrm{a}}}{m_{\mathrm{e}}}\int\frac{\mathrm{d}\boldsymbol{v} \,}{v}\left(\mathbb{I} -\frac{\boldsymbol{v}\boldsymbol{v}}{v^{2}} \right)\boldsymbol{\cdot}\frac{\partial }{\partial \boldsymbol{v}}\left(f_{\mathrm{e}}^{\mathrm{A}} + \frac{\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{u}_{\mathrm{i}}^{\mathrm{a}}}{v}\frac{\partial f_{\mathrm{e}}^{\mathrm{I}}}{\partial v} \right) \nonumber\\ & ={-}\frac{\gamma_{\mathrm{ei}}n_{\mathrm{i}}^{\mathrm{a}}}{m_{\mathrm{e}}}\left[\int \frac{\mathrm{d}\boldsymbol{v} \,}{v} \left(\mathbb{I} -\frac{\boldsymbol{v}\boldsymbol{v}}{v^{2}} \right)\frac{\partial f_{\mathrm{e}}^{\mathrm{I}}}{\partial v} \right]\boldsymbol{u}_{\mathrm{i}}^{\mathrm{a}} - \frac{2\gamma_{\mathrm{ei}}n_{\mathrm{i}}^{\mathrm{a}}}{m_{\mathrm{e}}}\int \,\mathrm{d}\boldsymbol{v} \,\frac{\boldsymbol{v}}{v^{3}}f_{\mathrm{e}}^{\mathrm{\mathrm{A}}} \nonumber\\ & = \frac{8{\rm \pi} \gamma_{\mathrm{ei}}n_{\mathrm{i}}^{\mathrm{a}}f_{\mathrm{e}}^{\mathrm{I}}(0)}{3m_{\mathrm{e}}}\boldsymbol{u}_{\mathrm{i}}^{\mathrm{a}} - \frac{2\gamma_{\mathrm{ei}}n_{\mathrm{i}}^{\mathrm{a}}}{m_{\mathrm{e}}}\int\,\mathrm{d}\boldsymbol{v} \,\frac{\boldsymbol{v}}{v}f_{\mathrm{e}}^{\mathrm{A}}, \end{align}

where $f_{\mathrm {e}}^{\mathrm {I}}(v)$ and $f_{\mathrm {e}}^{\mathrm {A}}(\boldsymbol {v})$ are, respectively, the isotropic and anisotropic parts of the electron phase-space density. Going from the third to the fourth line we used the isotropy of $f_{\mathrm {e}}^{\mathrm {I}}$ to compute the angle integral explicitly. The same cannot be done for the integral over the anisotropic part of the phase-space density. However, if we assume that, to lowest order, the electron phase-space density is isotropic around the mean electron velocity $\boldsymbol {u}_{\mathrm {e}}$, then

(6.15)\begin{equation} f_{\mathrm{e}}(\boldsymbol{v}) = f_{\mathrm{e}}^{\mathrm{I}}(|\boldsymbol{v}-\boldsymbol{u}_{\mathrm{e}}|) \approx f_{\mathrm{e}}^{\mathrm{I}} - \frac{\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{u}_{\mathrm{e}}}{v}\frac{\partial f_{\mathrm{e}}^{\mathrm{I}}}{\partial v} \implies f_{\mathrm{e}}^{\mathrm{A}}(\boldsymbol{v}) ={-}\frac{\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{u}_{\mathrm{e}}}{v}\frac{\partial f_{\mathrm{e}}^{\mathrm{I}}}{\partial v}, \end{equation}

which allows the final integral to be computed, giving the anomalous drag force

(6.16)\begin{equation} \boldsymbol{F}_{\mathrm{ei}} = m_{\mathrm{e}}n_{\mathrm{e}}\frac{\mathrm{d}\boldsymbol{u}_{\mathrm{e}}}{\mathrm{d}t} ={-}\frac{8{\rm \pi} \gamma_{\mathrm{ei}}n_{\mathrm{i}}^{\mathrm{a}}f_{\mathrm{e}}^{\mathrm{I}}(0)}{3m_{\mathrm{e}}}\left(\boldsymbol{u}_{\mathrm{e}} - \boldsymbol{u}_{\mathrm{i}}^{\mathrm{a}} \right). \end{equation}

This is the same as the standard expression for the collisional drag (Spitzer Reference Spitzer1967), but replacing the ion density and mean velocity with the corresponding anomalous variants (6.4).

We note that neither the lowest- nor the first-order approximations of the electron–ion collision operator yet describe the transfer of energy between the two species. To achieve this, and obtain the expected relaxation of temperature, one must go to higher order, which is easiest to do via the ion–electron collision operator.

6.4. Ion–electron relaxation: temperature equilibration

The ion–electron collision operator (6.8) is

(6.17)\begin{align} C_{\mathrm{ie}} & = \dfrac{\gamma_{\mathrm{ei}}}{m_{\mathrm{i}}}\frac{\partial }{\partial \boldsymbol{v}}\int \dfrac{\mathrm{d}\boldsymbol{v}' \,}{w}\left(\mathbb{I} -\dfrac{\boldsymbol{w}\boldsymbol{w}}{w^{2}} \right)\boldsymbol{\cdot} \int\,\mathrm{d}\eta' \,\eta'\nonumber\\ & \qquad \Bigg\{ \underbrace{\dfrac{\Delta\Gamma_{\mathrm{e}}}{m_{\mathrm{i}}}\left[\eta' - f_{\mathrm{e}}(\boldsymbol{v}') \right]P_{\mathrm{e}}(\boldsymbol{v}',\eta')\frac{\partial P_{\mathrm{i}}}{\partial \boldsymbol{v}}}_{{\sim} \dfrac{\Delta\Gamma_{\mathrm{e}}\eta_{\mathrm{e}}^{\mathrm{eff}}}{m_{\mathrm{i}}v_{\mathrm{thi}}}P_{\mathrm{e}}P_{\mathrm{i}}} - \underbrace{\dfrac{\Delta\Gamma_{\mathrm{i}}}{m_{\mathrm{e}}}\left[\eta - f_{\mathrm{i}}(\boldsymbol{v}) \right]P_{\mathrm{i}}(\boldsymbol{v},\eta)\frac{\partial P_{\mathrm{e}}}{\partial \boldsymbol{v}'}}_{{\sim} \dfrac{\Delta\Gamma_{\mathrm{i}}\eta_{\mathrm{i}}^{\mathrm{eff}}}{m_{\mathrm{e}}v_{\mathrm{the}}}P_{\mathrm{e}}P_{\mathrm{i}}} \Bigg\}. \end{align}

Naively, this gives an ion–electron relaxation rate that is comparable to the ion–ion relaxation rate, which is smaller than the electron–ion one by a factor of $\sqrt {m_{\mathrm {e}}/m_{\mathrm {i}}}$. However, this neglects the fact that the electron distribution will isotropise itself on the electron–ion relaxation timescale. Thus, the lowest-order term in (6.17) will vanish before it ever gets a chance to participate. This renders the ion–electron relaxation rate smaller than the ion–ion one by another factor of $\sqrt {m_{\mathrm {e}}/m_{\mathrm {i}}}$. As in the case of true collisions, this does not preclude the conservation of total momentum because an order-unity velocity change of the electrons only requires an order-mass-ratio adjustment to the ion velocity in order to conserve momentum, and the ion–electron collision operator is still fully capable of achieving a change this small on the electron–ion collision timescale.

Using the fact that the electron distribution will long since have become isotropic to lowest order, we can again carry out a mass-ratio expansion of the ion–electron collision operator. Since the first bracketed term in (6.17) is nominally smaller, the lowest-order contribution only requires the isotropic part of the electron distribution function. This becomes (noting that now the integration is over electron velocities so, $\boldsymbol {v}'' \gg \boldsymbol {v}$)

(6.18)\begin{equation} \int\dfrac{\mathrm{d}\boldsymbol{v}' \,}{v'}\left(\mathbb{I} - \frac{\boldsymbol{v}'\boldsymbol{v}'}{v'^{2}} \right)\int \,\mathrm{d}\eta' \,\eta'\left[\eta' - f_{\mathrm{e}}(\boldsymbol{v}') \right]P_{\mathrm{e}}(\boldsymbol{v}',\eta) = \frac{8{\rm \pi}}{3}\mathbb{I} \int_{0}^{\infty}\langle (g_{\mathrm{e}}^{\mathrm{I}})^{2} \rangle(\boldsymbol{v}')v'\,\mathrm{d}v' \,. \end{equation}

The second bracketed term in (6.17) requires both the isotropic and anisotropic contributions, leading us to evaluate the integral

(6.19)\begin{align} \int\frac{\mathrm{d}\boldsymbol{v}' \,}{w}\left(\mathbb{I} - \frac{\boldsymbol{w}\boldsymbol{w}}{w^{2}} \right)\boldsymbol{\cdot} \frac{\partial f_{\mathrm{e}}}{\partial \boldsymbol{v}'} & = \int \frac{\mathrm{d}\boldsymbol{v}' \,}{w}\left(\mathbb{I} - \frac{\boldsymbol{w}\boldsymbol{w}}{w^{2}} \right)\boldsymbol{\cdot} \left(\frac{\partial f_{\mathrm{e}}^{\mathrm{A}}}{\partial \boldsymbol{v}'} + \frac{\boldsymbol{v}'}{v'}\frac{\partial f_{\mathrm{i}}^{\mathrm{I}}}{\partial v'} \right) \nonumber\\ & = \int \frac{\mathrm{d}\boldsymbol{v}' \,}{v'}\left(\mathbb{I} - \frac{\boldsymbol{v}'\boldsymbol{v}'}{v'^{2}} \right)\boldsymbol{\cdot} \frac{\partial f_{\mathrm{e}}^{\mathrm{A}}}{\partial \boldsymbol{v}'} + \left[\int\frac{\boldsymbol{v}'}{v'^{2}}\left(\mathbb{I} -\frac{\boldsymbol{v}'\boldsymbol{v}'}{v'^{2}} \right)\frac{\partial f_{\mathrm{e}}^{\mathrm{I}}}{\partial v'} \right]\boldsymbol{\cdot}\boldsymbol{v} \nonumber\\ & = 2\int \,\mathrm{d}\boldsymbol{v}' \,\frac{\boldsymbol{v}'}{v'^{2}}f_{\mathrm{e}}^{\mathrm{A}}(\boldsymbol{v}') - \frac{8{\rm \pi} f_{\mathrm{e}}^{\mathrm{I}}(0)}{3}\boldsymbol{v}. \end{align}

Collecting the contributions (6.18) and (6.19) together, we find that the electron–ion collision operator is, to lowest order,

(6.20)\begin{gather} C_{\mathrm{ie}} ={-}\frac{\gamma_{\mathrm{ei}}\Delta\Gamma_{\mathrm{i}}}{m_{\mathrm{e}}m_{\mathrm{i}}}\frac{\partial }{\partial \boldsymbol{v}}\boldsymbol{\cdot}\left\{ \left[\eta - f_{\mathrm{i}}(\boldsymbol{v}) \right]P_{\mathrm{i}}(\boldsymbol{v},\eta)\left[2\int \,\mathrm{d}\boldsymbol{v}' \,\frac{\boldsymbol{v}'}{v'^{2}}f_{\mathrm{e}}^{\mathrm{A}}(\boldsymbol{v}') - \frac{8{\rm \pi} f_{\mathrm{e}}^{\mathrm{I}}(0)}{3}\boldsymbol{v} \right]\right\} \nonumber\\ + \,\frac{8{\rm \pi} \gamma_{\mathrm{ei}}\Delta\Gamma_{\mathrm{e}}}{3m_{\mathrm{i}}^{2}}\int\langle (g_{\mathrm{e}}^{\mathrm{I}})^{2} \rangle(v')v'\,\mathrm{d}v' \,\frac{\partial^{2} P_{\mathrm{i}}}{\partial \boldsymbol{v}^{2}}. \end{gather}

While not amazingly insightful, it is at least obvious how to verify that, together with the electron–ion collision operator (6.12), this conserves the total momentum of the system, giving the opposite momentum change to (6.14), as it must.

Of course, we need not have stopped at the simplest assumption that the electron distribution is, to lowest order, isotropic. The electron–electron interactions, which are just as frequent as the electron–ion interactions, will push the electron distribution function towards a Lynden-Bell equilibrium, which will have some mean velocity $\boldsymbol {u}_{\mathrm {e}}$ that lies close to the anomalous ion velocity $\boldsymbol {u}_{\mathrm {i}}^{\mathrm {a}}$ given by (6.4). Such a Lynden-Bell equilibrium for the electrons has the form

(6.21)\begin{gather} f_{\mathrm{e}}(\boldsymbol{v}) = \frac{\int \,\mathrm{d}\eta \,\eta {\rm e}^{-\beta_{\mathrm{e}}\Delta\Gamma_{\mathrm{e}}\eta\left[\frac{1}{2}m_{\mathrm{e}}|\boldsymbol{v}- \boldsymbol{u}_{\mathrm{e}}|^{2} - \mu_{\mathrm{e}}(\eta) \right]}}{\int \,\mathrm{d}\eta \, {\rm e}^{-\beta_{\mathrm{e}}\Delta\Gamma_{\mathrm{e}}\eta\left[\frac{1}{2}m_{\mathrm{e}}|\boldsymbol{v}- \boldsymbol{u}_{\mathrm{e}}|^{2} - \mu_{\mathrm{e}}(\eta) \right]}}. \end{gather}

This immediately allows $f_{\mathrm {e}}^{\mathrm {A}}$ in (6.20) to be computed as in (6.15). Using the property of Lynden-Bell equilibria expressed by (5.5) and (5.6), we also find

(6.22)\begin{equation} \langle (g_{\mathrm{e}}^{\mathrm{I}})^{2} \rangle ={-}\frac{1}{\beta_{\mathrm{e}}\Delta\Gamma_{\mathrm{e}}m_{\mathrm{e}}v}\frac{\partial f_{\mathrm{e}}^{\mathrm{I}}}{\partial v}. \end{equation}

The ion–electron collision operator (6.20) then reduces to

(6.23)\begin{gather} C_{\mathrm{ie}}[P_{\mathrm{i}},P_{\mathrm{e}}] = \frac{8{\rm \pi} \gamma_{\mathrm{ei}}\Delta\Gamma_{\mathrm{i}}f_{\mathrm{e}}^{\mathrm{I}}(0)}{3m_{\mathrm{i}}m_{\mathrm{e}}}\frac{\partial }{\partial \boldsymbol{v}}\boldsymbol{\cdot}\left\{ \frac{1}{\beta_{\mathrm{e}}\Delta\Gamma_{\mathrm{i}}m_{\mathrm{i}}}\frac{\partial P_{\mathrm{i}}}{\partial \boldsymbol{v}} + \left(\boldsymbol{v}-\boldsymbol{u}_{\mathrm{e}} \right)\left[\eta - f_{\mathrm{i}}(\boldsymbol{v}) \right]P_{\mathrm{i}}(\boldsymbol{v},\eta)\right\}. \end{gather}

Clearly, this will be stationary for a Lynden-Bell equilibrium if the mean electron velocity $\boldsymbol {u}_{\mathrm {e}}$ is equal to the mean ion velocity (i.e., if neither distribution has a mean velocity in the frame moving with $\boldsymbol {u}_{\mathrm {e}}$), and the thermodynamic beta $\beta _{\mathrm {i}}$ of the ions’ Lynden-Bell equilibrium is equal to that of the electrons, $\beta _{\mathrm {e}}$.

Because the ion–ion interactions increase the entropy of the ion hyperkinetic distribution function and are faster than the ion–electron interactions, it is certainly true that the ion distribution will be a Lynden-Bell equilibrium. Furthermore, because the anomalous ion velocity $\boldsymbol {u}_{\mathrm {i}}^{\mathrm {a}}$ is equal to the mean ion velocity for isotropic distributions like the Lynden-Bell equilibria, the electron–ion collision operator will have seen to it that the mean ion and electron velocities match. However, since the electron–electron, electron–ion and ion–ion collision operators to lowest order do not alter the energies of the two distributions, it is not guaranteed that the electrons and ions will have the same thermodynamic beta. Thus, as anticipated, the final piece of the relaxation puzzle, the relaxation of temperatures, occurs on the ion–electron interaction timescale. By assuming a Lynden-Bell equilibrium for ions with mean velocity equal to the electron mean velocity (i.e., both species stationary in the zero momentum frame), we can calculate the rate of change of the kinetic energy of the ion distribution function:

(6.24)\begin{align} \frac{\mathrm{d}E_{\mathrm{i}}}{\mathrm{d}t} & = \frac{\mathrm{d}}{\mathrm{d}t}\int\,\mathrm{d}\boldsymbol{v} \,\frac{1}{2}m_{\mathrm{i}}v^{2}\int \,\mathrm{d}\eta \,\eta P_{\mathrm{i}}(\boldsymbol{v},\eta) \nonumber\\ & = \frac{8{\rm \pi} \gamma_{\mathrm{ei}}f_{\mathrm{e}}^{\mathrm{I}}(0)}{3m_{\mathrm{i}}m_{\mathrm{e}}}\left[ \frac{3}{\beta_{\mathrm{e}}}\int \,\mathrm{d}\boldsymbol{v} \,f_{\mathrm{i}}(\boldsymbol{v}) - m_{\mathrm{i}}\Delta\Gamma_{\mathrm{i}}\int \,\mathrm{d}\boldsymbol{v} \,|\boldsymbol{v}|^{2}\langle g_{\mathrm{i}}^{2} \rangle(\boldsymbol{v}) \right] \nonumber\\ & = \frac{8{\rm \pi} \gamma_{\mathrm{ei}}f_{\mathrm{e}}^{\mathrm{I}}(0)}{3m_{\mathrm{i}}m_{\mathrm{e}}}\left(\frac{3n_{\mathrm{i}}}{\beta_{\mathrm{e}}} + \frac{1}{\beta_{\mathrm{i}}}\int \,\mathrm{d}\boldsymbol{v} \,\boldsymbol{v}\boldsymbol{\cdot}\frac{\partial f_{\mathrm{i}}}{\partial \boldsymbol{v}} \right) \nonumber\\ & = \frac{8{\rm \pi} \gamma_{\mathrm{ei}}f_{\mathrm{e}}^{\mathrm{I}}(0)n_{\mathrm{i}}}{3m_{\mathrm{i}}m_{\mathrm{e}}} \left(\frac{1}{\beta_{\mathrm{e}}} - \frac{1}{\beta_{\mathrm{i}}} \right), \end{align}

where, in going to the second line we integrated by parts, in going to the third line we exploited the relationship (6.22) between $\langle g_{\mathrm {i}}^{2} \rangle$ and $f_{\mathrm {i}}$, and finally integrated by parts again in going to the last line. This is simply another manifestation of the second law of thermodynamics (guaranteed by the H-theorem): energy will flow from the thermodynamically hotter species to the thermodynamically colder one until temperatures equalise. Note, however, that this energy flow need not produce the equilibration of kinetic energies, which are no longer directly tied to the thermodynamic temperatures, the equilibria being non-Maxwellian.