q_{\alpha}^{2}\Delta\Gamma_{\alpha}}{k^{2}}\int\mathrm{d}\boldsymbol{v}\,\frac{1}{{|\epsilon_{\boldsymbol{k},\boldsymbol{k}\cdot\boldsymbol{v}}|^{2}}}\iint\mathrm{d}\eta\,\mathrm{d}\eta^{\prime}\,\eta\eta^{\prime}&\Big{[}\delta(\eta-\eta^{\prime})P_{0\alpha}(\boldsymbol{v},\eta)\\\
&\quad-
P_{0\alpha}(\boldsymbol{v},\eta)P_{0\alpha}(\boldsymbol{v},\eta^{\prime})\Big{]}\end{aligned}$
$\displaystyle\begin{aligned}
\mathllap{}=\sum_{\alpha}\sum_{\boldsymbol{k}}\frac{2\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}\cdot\boldsymbol{v}}|^{2}},\end{aligned}$
(129)

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

In a Lynden-Bell equilibrium, (5.3) can be further simplified by calculating
$\Delta\Gamma_{\alpha}\langle g_{\alpha}^{2}\rangle$ via (120) and (121).
Approximating also
${|\epsilon_{\boldsymbol{k},\boldsymbol{k}\cdot\boldsymbol{v}}|^{2}\approx
1}$, we have

$E_{\varphi}\approx\sum_{\alpha}\sum_{\boldsymbol{k}}\frac{2\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),$ (130)

where $k_{\mathrm{max}}$ is the UV cutoff for the wavenumber integal (in 3D).
Using (121) and making a rough estimate of everything, we have

$E_{\varphi}\sim\sum_{\alpha}\Delta\Gamma_{\alpha}\eta_{\alpha}^{\mathrm{eff}}q_{\alpha}^{2}n_{\alpha}k_{\mathrm{max}}V.$
(131)

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 (118). To do this, we set

$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}),$
(132)

where $\eta_{\alpha}=\Delta\Gamma_{\alpha}^{-1}$ according to (117). Inserting
(132) into (5.3) and collecting only the lowest-order terms in the limit of
$f_{0\alpha}/\eta_{\alpha}\to 0$ gives

$E_{\varphi,\mathrm{BL}}=\sum_{\alpha}\sum_{\boldsymbol{k}}\frac{2\pi
q_{\alpha}^{2}}{k^{2}}\int\mathrm{d}\boldsymbol{v}\,\frac{f_{0\alpha}(\boldsymbol{v})}{|\epsilon_{\boldsymbol{k},\boldsymbol{k}\cdot\boldsymbol{v}}|^{2}}\approx\sum_{\alpha}4Vq_{\alpha}^{2}k_{\mathrm{max}}n_{\alpha},$
(133)

where the last estimate has been obtained in the same manner as (130). From a
direct comparison of (5.3) and (131) with (133), 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.222As with the estimates of the collision
timescale, we do not expect this factor to be large for distributions
sufficiently close to their ground states (Gardner distributions) when
${\eta_{\alpha}^{\mathrm{eff}}\ll n_{\alpha}/v_{\mathrm{th}\alpha}^{3}}$, but
such systems will barely evolve anyway.

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

$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}.$ (134)

Therefore, using (131), we find

$\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,$ (135)

ignoring at the last step any potential disparities between contributions from
different species. Taking cue from our discussion in section 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 section 2.3), their kinetic energy does not
change [see (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 (131).

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 (133).
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

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 (125) and (131), 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 section
5.2 that, to be in any way important, our relaxation had to be faster than
phase-volume conservation could be broken. While estimates (123) and (124)
serve as guesses for how fast phase-volume conservation is broken, there is
numerical evidence (Zhdankin, 2021) 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. Eyink 2018,
Schekochihin et al. 2009). If this is the case, then one would not be
justified in writing a collisionless Vlasov equation (1) for $f_{\alpha}$ or
(82) for $P_{\alpha}$, rendering the subsequent calculations formally invalid.
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.3) 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.3), 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 (135).

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 (80) 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,
2006).

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 Strange relaxation in multispecies Lynden-Bell plasma

Having derived the hyperkinetic collision integral, its conservation laws,
steady states and H-theorem, it would seem that relaxation to equilibrium has
successfully been turned into an app. Like a figure of Greek myth, the
collision integral is then destined to be made redundant by its own offspring:
the H-theorem. This is because the H-theorem prescribes for each initial
condition a steady-state distribution function, which makes evolving the
hyperkinetic collision integral in time seem moot. In this section, however,
we will show that there are certain fairly general initial conditions for
which the process of relaxation has interesting and non-trivial properties.

For this to be the case, it is clear that there must be multiple timescales
within the problem. Based on our estimate (120) of effective collision
frequency relative to the true collision frequency, it is clear that one way
to make this possible is to consider two species with a large mass ratio:
electrons and ions. Borrowing intuition from the conventional collisional
theory, this will allow for a separation of timescales whereby particles of
the lighter, faster, species will collide most frequently while the heavier,
slower, species will collide, and therefore evolve, at a lower rate. This
intuitive picture is, however, complicated by the fact that, from our
knowledge of the steady states (98), we expect the thermal velocities of the
particles to be such that

$\beta\Delta\Gamma_{\alpha}\eta_{\alpha}^{\mathrm{eff}}\frac{1}{2}m_{\alpha}v_{\mathrm{th}\alpha}^{2}\sim
1\implies
v_{\mathrm{th}\alpha}\sim\sqrt{\frac{2}{\beta\Delta\Gamma_{\alpha}^{\mathrm{eff}}\eta_{\alpha}^{\mathrm{eff}}m_{\alpha}}}.$
(136)

This is to say that, relative to other species, particles of a given species
behave as though they had an effective mass
$\Delta\Gamma_{\alpha}\eta_{\alpha}^{\mathrm{eff}}$ times greater than their
true mass. To avoid this complexity, we will restrict ourselves to the case
where the typical number of particles within a correlation volume is
comparable for electrons and ions, viz.,

$\Delta\Gamma_{\mathrm{e}}\eta_{\mathrm{e}}^{\mathrm{eff}}\sim\Delta\Gamma_{\mathrm{i}}\eta_{\mathrm{i}}^{\mathrm{eff}}.$
(137)

This amounts to assuming that collisionless effects do not override the scale
separation imposed by the mass ratio, which is what creates the interesting
effects in the conventional collisional theory. From (95) and the orderings
(136) and (137), it is apparent that the same-species collision operators only
contain a single timescale. This is simply a statement that the only purpose
of the same-species collision operator is to increase the entropy of that
species independent of the others, and that it can only do that on a single
timescale. In contrast, owing to the mass ratio, the inter-species collisions
can have multiple timescales. We know from standard collisional theory that
these give rise to the rates at which momentum and temperature are equalised
between species. To extract similar, but distinct, features from our new
collision integral, we would therefore have to expand it in small mass ratio.
Before we do this, however, it is possible to anticipate from simple
observations what the interesting new physics will be.

###  6.1 Preview of strange relaxation

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

$\frac{\partial P_{\alpha}}{\partial
t}=\sum_{\alpha^{\prime}}C_{\alpha\alpha^{\prime}}[P_{\alpha},P_{\alpha^{\prime}}],$
(138)

in terms of the interspecies collision operator $C_{\alpha\alpha^{\prime}}$,
which can be easily read from elements of the species sum in (95). For
compactness of notation, we shall henceforth drop 0’s from the subscripts of
the mean distribution function. For electrons interacting with ions, the
dominant effect is the diffusion-like first term of (95). Comparing this to
the true-collision operator (118), 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 (112). 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

$n_{\mathrm{i}}^{\mathrm{a}}=\Delta\Gamma_{\mathrm{i}}\int\mathrm{d}\boldsymbol{v}\,^{\prime}\langle
g_{\mathrm{i}}^{2}\rangle(\boldsymbol{v}^{\prime}),\quad\quad\boldsymbol{u}_{\mathrm{i}}^{\mathrm{a}}=\frac{\int\mathrm{d}\boldsymbol{v}^{\prime}\,\boldsymbol{v}^{\prime}\langle
g_{\mathrm{i}}^{2}\rangle(\boldsymbol{v}^{\prime})}{\int\mathrm{d}\boldsymbol{v}^{\prime}\,\langle
g_{\mathrm{i}}^{2}\rangle(\boldsymbol{v}^{\prime})}.$ (139)

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).

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})\left[\eta-
f_{\mathrm{i}}(\boldsymbol{v})\right]$, which is precisely $\langle
g_{\mathrm{i}}^{2}\rangle(\boldsymbol{v})$ from which the anomalous ion
velocity and density are defined in (139). 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 (section 6.3.2) and
the ion-electron (section 6.4) collision operators. This calculation follows
the standard path well trodden for true collisions, but with a few important
adjustments.

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 (44)], 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.

###  6.2 Landau form of the hyperkinetic collision operator

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

$\frac{1}{V}\sum_{\boldsymbol{k}}\frac{\boldsymbol{k}\boldsymbol{k}}{k^{4}}\frac{\delta\left(\boldsymbol{k}\cdot(\boldsymbol{v}-\boldsymbol{v}^{\prime})\right)}{|\epsilon_{\boldsymbol{k},\boldsymbol{k}\cdot\boldsymbol{v}}|^{2}},$
(140)

which is complicated by the presence of the dielectric function

$\epsilon_{\boldsymbol{k},\boldsymbol{k}\cdot\boldsymbol{v}}=1+\sum_{\alpha^{\prime}}\frac{4\pi
q_{\alpha^{\prime}}^{2}}{m_{\alpha^{\prime}}k^{2}}\iint\mathrm{d}\boldsymbol{v}^{\prime}\,\mathrm{d}\eta^{\prime}\,\frac{\eta^{\prime}}{\boldsymbol{k}\cdot\left(\boldsymbol{v}-\boldsymbol{v}^{\prime}\right)}\boldsymbol{k}\cdot\frac{\partial
P_{\alpha^{\prime}}}{\partial\boldsymbol{v}}.$ (141)

Here the $\boldsymbol{v}^{\prime}$ integral is taken along the Landau contour.
To simplify (140), 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 (141) will be small and
$\epsilon_{\boldsymbol{k},\boldsymbol{k}\cdot\boldsymbol{v}}$ can be
approximated by unity. For scales significantly longer than the Debye length,
$\epsilon_{\boldsymbol{k},\boldsymbol{k}\cdot\boldsymbol{v}}\sim(\lambda_{\mathrm{De}}k)^{-2}$,
which will make their contribution to the sum in (140) small in
$(\lambda_{\mathrm{De}}k)^{4}$. Therefore, following Landau (1936), 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}^{\prime}$ and integrating in
cylindrical coordinates such that
$\boldsymbol{k}=k_{\parallel}\hat{\boldsymbol{w}}+\boldsymbol{k}_{\perp}$, we
get

$\frac{1}{V}\sum_{\boldsymbol{k}}\frac{\boldsymbol{k}\boldsymbol{k}}{k^{4}}\frac{\delta(\boldsymbol{k}\cdot\boldsymbol{w})}{|\epsilon_{\boldsymbol{k},\boldsymbol{k}\cdot\boldsymbol{v}}|^{2}}\approx\iiint\frac{\mathrm{d}\boldsymbol{k}\,}{\left(2\pi\right)^{3}}\frac{\boldsymbol{k}_{\perp}\boldsymbol{k}_{\perp}}{k_{\perp}^{4}}\delta(k_{\parallel}w)=\frac{1}{8\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).$
(142)

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
$\Lambda_{\alpha\alpha^{\prime}}$.

With the approximation (142), the interspecies collision operator in (95)
becomes

$C_{\alpha\alpha^{\prime}}[P_{\alpha},P_{\alpha^{\prime}}]=\frac{\gamma_{\alpha\alpha^{\prime}}}{m_{\alpha}}\frac{\partial}{\partial\boldsymbol{v}}\cdot\int\frac{\mathrm{d}\boldsymbol{v}\,^{\prime}}{w}\left(\mathbb{I}-\frac{\boldsymbol{w}\boldsymbol{w}}{w^{2}}\right)\cdot\int\mathrm{d}\eta^{\prime}\,\eta^{\prime}\\\
\Bigg{\\{}\frac{\Delta\Gamma_{\alpha^{\prime}}}{m_{\alpha}}\left[\eta^{\prime}-f_{\alpha^{\prime}}(\boldsymbol{v}^{\prime})\right]P_{\alpha^{\prime}}(\boldsymbol{v}^{\prime},\eta^{\prime})\left.\frac{\partial
P_{\alpha}}{\partial\boldsymbol{v}}\right|_{\eta}-\frac{\Delta\Gamma_{\alpha}}{m_{\alpha^{\prime}}}\left[\eta-
f_{\alpha}(\boldsymbol{v})\right]P_{\alpha}(\boldsymbol{v},\eta)\left.\frac{\partial
P_{\alpha^{\prime}}}{\partial\boldsymbol{v}^{\prime}}\right|_{\eta^{\prime}}\Bigg{\\}},$
(143)

where $\gamma_{\alpha\alpha^{\prime}}=2\pi
q_{\alpha}^{2}q_{\alpha^{\prime}}^{2}\Lambda_{\alpha\alpha^{\prime}}$.

###  6.3 Electron-ion relaxation

####  6.3.1 Zeroth order: istropisation of electrons

To lowest order in $\sqrt{m_{\mathrm{e}}/m_{\mathrm{i}}}$, the second term in
brackets in (143) 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

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

where $n_{\mathrm{i}}^{\mathrm{a}}$ is the anomalous ion density given by the
first expression in (139). This operator is a pitch-angle-scattering (Lorentz)
operator (see, e.g., Helander & Sigmar 2005), 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}^{\prime}$ integral in (143) is over ion velocities,
$\boldsymbol{v}^{\prime}$ 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 (143) to
first order in $\sqrt{m_{\mathrm{e}}/m_{\mathrm{i}}}$, giving

$\begin{split}\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}^{\prime}\cdot\frac{\partial}{\partial\boldsymbol{v}}\left[\frac{1}{v}\left(\mathbb{I}-\frac{\boldsymbol{v}\boldsymbol{v}}{v^{2}}\right)\right]\\\
&=\frac{1}{v}\left(1+\frac{3\boldsymbol{v}\cdot\boldsymbol{v}^{\prime}}{v^{2}}\right)\left(\mathbb{I}-\frac{\boldsymbol{v}\boldsymbol{v}}{v^{2}}\right)+\frac{\boldsymbol{v}^{\prime}\boldsymbol{v}+\boldsymbol{v}\boldsymbol{v}^{\prime}-2(\boldsymbol{v}\cdot\boldsymbol{v}^{\prime})\mathbb{I}}{v^{3}}.\end{split}$
(145)

Since the anisotropic part of the hyperkinetic distribution function is
already small, we only need the lowest-order contribution from (145) to act on
$P_{\mathrm{e}}^{\mathrm{A}}(\boldsymbol{v},\eta)$. For the isotropic part
$P_{\mathrm{e}}^{\mathrm{I}}(\boldsymbol{v})$, 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

$\begin{split}\frac{1}{w}\left(\mathbb{I}-\frac{\boldsymbol{w}\boldsymbol{w}}{w^{2}}\right)\cdot\frac{\partial
P_{\mathrm{e}}^{\mathrm{I}}}{\partial\boldsymbol{v}}&\approx\frac{\boldsymbol{v}^{\prime}\boldsymbol{v}+\boldsymbol{v}\boldsymbol{v}^{\prime}-2(\boldsymbol{v}\cdot\boldsymbol{v}^{\prime})\mathbb{I}}{v^{3}}\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)\cdot\frac{\boldsymbol{v}^{\prime}}{v}\frac{\partial
P_{\mathrm{e}}^{\mathrm{I}}}{\partial
v}=\frac{1}{v}\left(\mathbb{I}-\frac{\boldsymbol{v}\boldsymbol{v}}{v^{2}}\right)\cdot\frac{\partial}{\partial\boldsymbol{v}}\left(\frac{\boldsymbol{v}\cdot\boldsymbol{v}^{\prime}}{v}\frac{\partial
P_{\mathrm{e}}^{\mathrm{I}}}{\partial v}\right).\end{split}$ (146)

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

$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}}\cdot\left\\{\frac{1}{v}\left(\mathbb{I}-\frac{\boldsymbol{v}\boldsymbol{v}}{v^{2}}\right)\cdot\frac{\partial}{\partial\boldsymbol{v}}\left[P_{\mathrm{e}}^{\mathrm{A}}(\boldsymbol{v},\eta)+\frac{\boldsymbol{v}\cdot\boldsymbol{u}_{\mathrm{i}}^{\mathrm{a}}}{v}\frac{\partial
P_{\mathrm{e}}^{\mathrm{I}}}{\partial v}\right]\right\\},$ (147)

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

Despite not knowing the exact evolution due to (147), 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 (147) to be isotropic is for it to be identically zero.
Therefore the fixed point of the collision operator (147) is

$\begin{split}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}\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{split}$
(148)

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

Thus, as prophesied in section 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

$\begin{split}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}}]\\\
&=-\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)\cdot\frac{\partial}{\partial\boldsymbol{v}}\left(f_{\mathrm{e}}^{\mathrm{A}}+\frac{\boldsymbol{v}\cdot\boldsymbol{u}_{\mathrm{i}}^{\mathrm{a}}}{v}\frac{\partial
f_{\mathrm{e}}^{\mathrm{I}}}{\partial v}\right)\\\
&=-\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}}}\\\
&=\frac{8\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{split}$
(149)

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

$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}\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}\cdot\boldsymbol{u}_{\mathrm{e}}}{v}\frac{\partial
f_{\mathrm{e}}^{\mathrm{I}}}{\partial v},$ (150)

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

$\boldsymbol{F}_{\mathrm{ei}}=m_{\mathrm{e}}n_{\mathrm{e}}\frac{\mathrm{d}\boldsymbol{u}_{\mathrm{e}}}{\mathrm{d}t}=-\frac{8\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).$
(151)

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

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 (143) is

$C_{\mathrm{ie}}=\frac{\gamma_{\mathrm{ei}}}{m_{\mathrm{i}}}\frac{\partial}{\partial\boldsymbol{v}}\int\frac{\mathrm{d}\boldsymbol{v}^{\prime}\,}{w}\left(\mathbb{I}-\frac{\boldsymbol{w}\boldsymbol{w}}{w^{2}}\right)\cdot\int\mathrm{d}\eta^{\prime}\,\eta^{\prime}\\\
\Bigg{\\{}\underbrace{\frac{\Delta\Gamma_{\mathrm{e}}}{m_{\mathrm{i}}}\left[\eta^{\prime}-f_{\mathrm{e}}(\boldsymbol{v}^{\prime})\right]P_{\mathrm{e}}(\boldsymbol{v}^{\prime},\eta^{\prime})\frac{\partial
P_{\mathrm{i}}}{\partial\boldsymbol{v}}}_{\sim\frac{\Delta\Gamma_{\mathrm{e}}\eta_{\mathrm{e}}^{\mathrm{eff}}}{m_{\mathrm{i}}v_{\mathrm{thi}}}P_{\mathrm{e}}P_{\mathrm{i}}}-\underbrace{\frac{\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}^{\prime}}}_{\sim\frac{\Delta\Gamma_{\mathrm{i}}\eta_{\mathrm{i}}^{\mathrm{eff}}}{m_{\mathrm{e}}v_{\mathrm{the}}}P_{\mathrm{e}}P_{\mathrm{i}}}\Bigg{\\}}.$
(152)

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 (152) 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 (152)
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}^{\prime\prime}\gg\boldsymbol{v}$)

$\int\frac{\mathrm{d}\boldsymbol{v}^{\prime}\,}{v^{\prime}}\left(\mathbb{I}-\frac{\boldsymbol{v}^{\prime}\boldsymbol{v}^{\prime}}{v^{\prime
2}}\right)\int\mathrm{d}\eta^{\prime}\,\eta^{\prime}\left[\eta^{\prime}-f_{\mathrm{e}}(\boldsymbol{v}^{\prime})\right]P_{\mathrm{e}}(\boldsymbol{v}^{\prime},\eta)=\frac{8\pi}{3}\mathbb{I}\int_{0}^{\infty}\langle(g_{\mathrm{e}}^{\mathrm{I}})^{2}\rangle(\boldsymbol{v}^{\prime})v^{\prime}\mathrm{d}v^{\prime}\,.$
(153)

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

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

Collecting the contributions (153) and (154) together, we find that the
electron-ion collision operator is, to lowest order,

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

While not amazingly insightful, it is at least obvious how to verify that,
together with the electron-ion collision operator (147), this conserves the
total momentum of the system, giving the opposite momentum change to (149), 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 (139). Such a Lynden-Bell
equilibrium for the electrons has the form

$f_{\mathrm{e}}(\boldsymbol{v})=\frac{\int\mathrm{d}\eta\,\eta
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\,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]}}.$
(156)

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

$\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}.$ (157)

The ion-electron collision operator (155) then reduces to

$C_{\mathrm{ie}}[P_{\mathrm{i}},P_{\mathrm{e}}]=\frac{8\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}}\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\\}.$
(158)

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:

$\begin{split}\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)\\\
&=\frac{8\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]\\\
&=\frac{8\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}\cdot\frac{\partial
f_{\mathrm{i}}}{\partial\boldsymbol{v}}\right)\\\
&=\frac{8\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{split}$
(159)

where, in going to the second line we integrated by parts, in going to the
third line we exploited the relationship (157) 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.

##  7 Summary and discussion

The existence of collision integrals that have the Lynden-Bell (1967)
equilibria as their steady-state solutions implies that these equilibria can
be reached dynamically. Two collision integrals have arisen: the multi-
waterbag collision integral (58), and the hyperkinetic collision integral
(95), both generalising the single-waterbag collision integral first derived
by Kadomtsev & Pogutse (1970). It is a key ingredient that these collision
integrals are equipped with H-theorems, confirming that the dynamical
evolution towards the Lynden-Bell equilibria is made inevitable by the
requirement to increase the Lynden-Bell entropy. The derivation of both
collision integrals had to contend with a closure problem: the evolution of
the mean phase-space density $f_{0\alpha}$ requires knowledge of the
correlation function of the exact phase-space density between separate parts
of phase space: ${\langle
f_{\alpha}(\boldsymbol{r},\boldsymbol{v})f_{\alpha^{\prime}}(\boldsymbol{r}^{\prime},\boldsymbol{v}^{\prime})\rangle}$,
which is not, in general, known. For both collision integrals, this problem is
partially resolved by the microgranulation ansatz (see section 2.5), a version
of which was first used by Kadomtsev & Pogutse (1970) to derive their
collision integral. The microgranulation ansatz posits that the exact phase-
space density of particles of species $\alpha$ is correlated over a small, but
non-zero, volume in phase space $\Delta\Gamma_{\alpha}$, but is perfectly
mixed over larger volumes. In contrast, one can derive the Balescu-Lenard
integral (118) describing ‘true’ particle collisions by assuming that all
particles are statistically independent, i.e., that a particle is only
correlated with itself. The effect of the microgranulation ansatz is to reduce
the problem of calculating the worrisome two point correlator $\langle
f_{\alpha}(\boldsymbol{r},\boldsymbol{v})f_{\alpha^{\prime}}(\boldsymbol{r}^{\prime},\boldsymbol{v}^{\prime})\rangle$
to one of calculating a more manageable one-point correlator $\langle
f_{\alpha}^{2}\rangle(\boldsymbol{v})$.

A closure is still required because the variance of a random quantity cannot,
in general, be determined by its mean. Kadomtsev & Pogutse (1970) restricted
the exact phase-space density $f_{0\alpha}(\boldsymbol{r},\boldsymbol{v})$ to
only two possible values, $\eta_{\alpha}$ or $0$, which immediately implied
${\langle f_{\alpha}^{2}\rangle=\eta_{\alpha}f_{0\alpha}}$, removing the
closure problem (see section 3.1). However, this single-waterbag model is
obviously extremely non-general. To move past it, in section 3.2, a scheme
motivated by statistical mechanics is employed, which can be traced back to
the treatment of geophysical turbulence by Chavanis (2005). By maximising the
Shannon entropy (51) subject to a continuum of constraints (78) (the ‘waterbag
content’ of the distribution function, or its Casimir invariants), one finds
that the correlator $\langle f_{\alpha}^{2}\rangle$ can be written implicitly
in terms of $f_{0\alpha}$. This leads to the multi-waterbag collision integral
(58), which grows the Lynden-Bell entropy (67) and has the Lynden-Bell
equilibria (60) as its only fixed points. Notably, if laterally to our main
purpose here, this formalism allows one to recover very efficiently the ‘true’
collision integrals: Balescu-Lenard (section 5.1) and the collision integrals
for fermionic and bosonic quantum plasmas (Appendix A).

An alternative route to solving the closure problem for $\langle
f_{\alpha}^{2}\rangle$ is to dodge it entirely. In section 4, instead of
trying to determine the evolution of the mean phase-space density
$f_{0\alpha}(\boldsymbol{v})$, a ‘hyperkinetic’ approach is introduced,
treating the exact phase-space density
$f_{\alpha}(\boldsymbol{r},\boldsymbol{v})$ as a random field and asking for
the probability of finding it to have a particular value $\eta$ at a
particular position in phase space. This method, pioneered by Severne & Luwel
(1980) (for discrete waterbags) and having its origins in vortex kinetics and
galactic dynamics (Chavanis et al., 1996), resolves the closure problem by
calculating $f_{0\alpha}$ and $\langle f_{\alpha}^{2}\rangle$ as the first and
second moments, respectively, of the ‘hyperkinetic distribution function’
$P_{0\alpha}(\boldsymbol{v},\eta)=\langle\delta(f_{\alpha}(\boldsymbol{r},\boldsymbol{v})-\eta)\rangle$.
This function is evolved by the hyperkinetic collision integral (95), which
also grows the Lynden-Bell entropy (95) and has the Lynden-Bell equilibria
(98) as its only fixed points.

It is immediately apparent that the statistical mechanical closure (75) used
to derive the multi-waterbag collision integral (58) is simply a special case
of the hyperkinetic distribution function. In section 4.3, we show that in
fact their connection is deeper: we prove that, should the hyperkinetic
distribution be placed precisely on the closure (75), the two collision
integrals would then be equivalent for all future times. This could imply that
entropy maximisation local in phase space is an inherent feature of the
hyperkinetic collision integral (95). For this to be so, there would have to
be a shorter effective collision timescale on which the hyperkinetic
distribution function relaxed towards the closure (75) before continuing its
relaxation towards a Lynden-Bell equilibrium on a longer timescale.

While conceptually fascinating, to be put on solid ground, the ‘collisionless
collision integals’ must be verified. While there is substantial evidence of
the validity (within certain regimes) of the Balescu-Lenard integral (118)
describing ‘true’ Coulomb collisions, it remains to be seen whether the
hyperkinetic collision integral, and, more generally and especially, the
microgranulation ansatz, are valid. We therefore outline a number potential
observable quantities that would be indicative of the collisionless
relaxation. Most obviously, the verification of Lynden-Bell equilibria
themselves would vindicate the assumption of phase-volume conservation despite
the fairly stringent formal limitation on it imposed by (125) (the requirement
that collisionless relaxation occur before the exact phase-space density is
filamented down to collisional scales), and potentially an even more stringent
limitation in systems where a dissipative anomaly is present and the phase-
volume conservation is broken on a time scale independent of ‘true’
collisionality (see discussion and references in section 5.4). As discussed in
section 3.3.1, the Lynden-Bell equilibria can be extremely varied, but there
is a direct correspondence between the waterbag content of the initial
condition and the ultimate equilibrium, which can clearly be tested. The
relaxation process itself would bear telltale signs, should it obey the
microgranulation ansatz. In section 6, we showed that a plasma obeying the
microgranulation ansatz has an effective collision rate much higher than the
true-collision rate associated with the Balescu-Lenard collision integral, and
that the typical energy stored in the electric fluctuations in such a plasma
is much greater than the fluctuation level arising from Coulomb collisions.
This creates a picture of collisionless relaxation mediated not by single-
particle collisions but by the effective collisions of larger correlated
volumes of phase space, which seed larger perturbations to the electric field.

In section 6, this collisionless relaxation is shown to have some curious
consequences due to inter-species interactions. Most interestingly, electrons
are dragged towards a mean velocity that is not necessarily the mean velocity
of the ions, as the case of ‘true’ collisions, but a certain anomalous
velocity associated with an anisotropic ion distribution, which can be non-
zero even if the mean ion velocity is zero (e.g., for an ion distribution that
carries a heat flux but no net momentum). A current and, therefore, magnetic
field, could be spontaneously generated by this mechanism.

Despite accommodating larger fluctuations in the electric field than the
collisional theory, the collisionless relaxation described above is still
assumed to take place in the quasilinear regime. It seems likely that there
exist regimes where the fluctuations are larger still, and therefore a full
nonlinear solution to the two-point correlation function of the phase-space
density is required. What these regimes are, and whether there is a limit in
which they recover the microgranulation ansatz, will be the subject of future
work.

Note. Shortly before this paper was submitted, a preprint by Chavanis (2021)
appeared, where he independently develops a theory of collisionless relaxation
along the lines that are, in certain respects, parallel to our reasoning, and
reaches some of the same (or similar) conclusions. We refer the reader to his
paper also for a detailed historical review of his and others’ work on
collisionless relaxation in a range of physical systems—primarily gravitating
and fluid-dynamical.

## Acknowledgements

We would like to thank Michael Barnes, Ben Chandran, Ilya Dodin, Bill Dorland,
Jean-Baptiste Fouvry, Chris Hamilton, Per Helander, Sid Parameswaran, Luis
Silva, and Dmitri Uzdensky for illuminating discussions. RJE would like to
further thank Felix Parra, whose lectures on collisional plasma physics
inspired the treatment of strange relaxation in section 6. We also thank
Pierre-Henri Chavanis for bringing to our attention his early work on
statistical-mechanical closures for collision integrals (Chavanis, 2004,
2005). RJE’s and TA’s work was supported by UK EPSRC studentships. TA’s work
was also supported by the Euratom research and training programme within the
framework of the EUROfusion consortium (grant agreement No. 633053) and by the
UKRI Energy Programme (grant EP/T012250/1). The views and opinions expressed
herein do not necessarily reflect those of the European Commission. AB was
supported by a Marshall Scholarship. The work of AAS was supported in part by
a UK EPSRC grant EP/R034737/1.

##  Appendix A Collision integrals in quantum plasmas

In this appendix, we show how the scheme used in sections 2 and 3 can be
applied to work out swiftly the ‘true’ collision integrals for fermions and
bosons. It is perhaps unsruprising that this is possible given the close
analogy between phase-volume conservation and the exclusion principle.
However, the procedure is non-rigorous because sections 2 and 3 assume point
particles occupying definite positions in the
($\boldsymbol{r},\boldsymbol{v}$) phase space which is not the case for
quantum gases. This issue can be handled rigorously by the derivation and
expansion of the Uehling & Uhlenbeck (1933) collision operator leading to the
quantum Balescu-Lenard collision operator and its Landau siblings (see, e.g.,
Danielewicz 1980). We instead opt for the cheap and fast route to arrive at
the (not-yet-closed) collision operator (35). This is written in terms of the
unknown correlation volume $\Delta\Gamma_{\alpha}$ and the correlator $\langle
g_{\alpha}^{2}\rangle(\boldsymbol{v})$, which we can resolve by falling back
on the quantum nature of our particles. To consider true collisions, we assume
that particles are uncorrelated wherever this is compatible with quantum
constraints (the exclusion principle). Crucially, however, due to the
quantisation of their positions and momenta, particles occupy – and,
therefore, must be correlated over – a finite phase space volume,
$\Delta\Gamma_{\alpha}$. We can infer this volume by considering particles in
a box of volume $V$, for which the possible momenta are

$\boldsymbol{p}=\frac{2\pi\hbar}{V^{1/3}}\left(i_{x},i_{y},i_{z}\right),$
(160)

where $i_{x}$, $i_{y}$, $i_{z}$ are integers. The phase-volume per particle
can then be read off from the spacing of their momenta:

$\Delta\Gamma_{\alpha}=\left(\frac{2\pi\hbar}{m_{\alpha}}\right)^{3}\equiv\eta_{0\alpha}^{-1},$
(161)

where we have defined the ‘quantum density’ $\eta_{0\alpha}$ (equivalent in
spirit to the ‘waterbag density’ considered in section 3) to be the inverse of
this correlation volume. The possible values of the ‘exact phase-space
density’ $f_{\alpha}$ are then simply $\eta_{0\alpha}$ multiplied by the
possible occupation numbers.

To calculate the correlator $\langle g_{\alpha}^{2}\rangle(\boldsymbol{v})$
(the variance of the initial condition, understood in the sense discussed in
section 2.5), we will appeal to the closure scheme in section 3.2, which will
account for the particle statistics. We consider all possible occupation
numbers at velocity $\boldsymbol{v}$, and maximise entropy subject to knowing
the mean phase-space density $f_{0\alpha}$ (equal to $\eta_{0\alpha}$ times
the mean occupation number). The correlator $\langle
g_{\alpha}^{2}\rangle(\boldsymbol{v})$ is then $\eta_{0\alpha}^{2}$ times the
variance of the occupation number given this maximum-entropy assignment. We
will consider the case where the bosons or fermions have $\sigma_{\alpha}$
possible spins (or indeed any quantum number giving rise to degeneracy). Then
the partition function becomes

$Z_{\alpha}=\sum_{n_{1},n_{2}...n_{\sigma_{\alpha}}}e^{-\psi_{\alpha}(\boldsymbol{v})\eta_{0\alpha}\sum_{j}n_{j}}=\left[\sum_{n}e^{-\psi_{\alpha}(\boldsymbol{v})\eta_{0\alpha}n}\right]^{\sigma_{\alpha}},$
(162)

where $\psi_{\alpha}(\boldsymbol{v})$ is a Lagrange-multiplier function chosen
as in (53) to guarantee the correct mean phase-space density
$f_{0\alpha}(\boldsymbol{v})$. The first sum in (162) is over all possible
occupation numbers $n_{j}$ of each of the possible spins. Since the effect on
the phase-space density due to occupation numbers from different spins is
additive, this is then split into the product of $\sigma_{\alpha}$ identical
sums. The remaining sum is taken over the allowed values of the occupation
number at each spin: $n\in\\{0,1\\}$ for fermions and $n\in\\{0,1,2,...\\}$
for bosons. The result is

$Z_{\alpha}=\left[1\pm
e^{-\psi_{\alpha}(\boldsymbol{v})\eta_{0\alpha}}\right]^{\pm\sigma_{\alpha}},$
(163)

with the ‘$+$’ sign for fermions and the ‘$-$’ sign for bosons. As in (57),
$\langle g_{\alpha}^{2}\rangle$ can now be computed thus:

$\langle g_{\alpha}^{2}\rangle=\frac{\partial^{2}\ln
Z_{\alpha}}{\partial\psi_{\alpha}^{2}}=\left(\eta_{0\alpha}\mp\frac{f_{0\alpha}}{\sigma_{\alpha}}\right)f_{0\alpha},$
(164)

with the ‘$-$’ sign for fermions and the ‘$+$’ sign for bosons.

Using (164) in (35) gives the desired collision operator for charged particles
obeying Bose-Einstein or Fermi-Dirac statistics:

$\frac{\partial f_{0\alpha}}{\partial
t}=\sum_{\alpha^{\prime}}\frac{16\pi^{3}q_{\alpha}^{2}q_{\alpha^{\prime}}^{2}}{m_{\alpha}V}\frac{\partial}{\partial\boldsymbol{v}}\cdot\sum_{\boldsymbol{k}}\frac{\boldsymbol{k}\boldsymbol{k}}{k^{4}}\cdot\int\mathrm{d}\boldsymbol{v}^{\prime}\,\frac{\delta(\boldsymbol{k}\cdot(\boldsymbol{v}-\boldsymbol{v}^{\prime}))}{\left|\epsilon_{\boldsymbol{k},\boldsymbol{k}\cdot\boldsymbol{v}}\right|^{2}}\\\
\left\\{\frac{1}{m_{\alpha}}\left[1\mp\frac{f_{0\alpha^{\prime}}(\boldsymbol{v}^{\prime})}{\sigma_{\alpha^{\prime}}\eta_{0\alpha^{\prime}}}\right]f_{0\alpha^{\prime}}(\boldsymbol{v}^{\prime})\frac{\partial
f_{0\alpha}}{\partial\boldsymbol{v}}-\frac{1}{m_{\alpha^{\prime}}}\left[1\mp\frac{f_{0\alpha}(\boldsymbol{v})}{\sigma_{\alpha}\eta_{0\alpha}}\right]f_{0\alpha}(\boldsymbol{v})\frac{\partial
f_{0\alpha^{\prime}}}{\partial\boldsymbol{v}^{\prime}}\right\\}.$ (165)

This clearly recovers the Balescu-Lenard collision operator (118) when the
system is nowhere near degeneracy (i.e., when
$f_{\alpha}\ll\sigma_{\alpha}\eta_{0\alpha}$). We further note that the
results regarding strange relaxation in section 7 are equally applicable (and
easier to interpret physically) for the collision operators (165). This
implies an anomalous resistivity of quantum plasmas resulting from degeneracy
prohibiting certain collisions. The modification to the electron-ion friction
force in quantum plasmas has been studied before (e.g., see Daligault 2016;
Rightley & Baalrud 2021) in systems close to the relevant equilibrium (i.e.,
after isotropisation has removed the possibility of
$\boldsymbol{u}_{\mathrm{i}}^{\mathrm{a}}\neq\boldsymbol{u}_{\mathrm{i}}$). We
note, however, that, owing to the scaling of (161) with species mass, ions are
conventionally less degenerate than electrons, so achieving
$\boldsymbol{u}_{\mathrm{i}}^{\mathrm{a}}\neq\boldsymbol{u}_{\mathrm{i}}$
requires extreme densities.

## References

  * Adkins & Schekochihin (2018) Adkins, T. & Schekochihin, A. A. 2018 A solvable model of Vlasov-kinetic plasma turbulence in Fourier–Hermite phase space. J. Plasma Phys. 84, 905840107.
  * Balescu (1960) Balescu, R. 1960 Irreversible processes in ionized gases. Phys. Fluids 3, 52.
  * Boltzmann (1896) Boltzmann, L. 1896 Vorlesugnen über Gastheorie. Leipzig: J. A. Barth.
  * Chavanis (2004) Chavanis, P.-H. 2004 Generalized thermodynamics and kinetic equations: Boltzmann, Landau, Kramers and Smoluchowski. Physica A 332, 89.
  * Chavanis (2005) Chavanis, P.-H. 2005 Statistical mechanics of geophysical turbulence: application to jovian flows and Jupiter’s great red spot. Physica D 200, 257.
  * Chavanis (2006) Chavanis, P.-H. 2006 Quasi-stationary states and incomplete violent relaxation in systems with long-range interactions. Physica A 365, 102.
  * Chavanis (2021) Chavanis, P.-H. 2021 Kinetic theory of collisionless relaxation for systems with long-range interactions, arXiv: 2112.13664.
  * Chavanis et al. (1996) Chavanis, P. H., Sommeria, J. & Robert, R. 1996 Statistical mechanics of two-dimensional vortices and collisionless stellar systems. Astrophys. J. 471, 385.
  * Coonin (2014) Coonin, A. 2014 From Marble to Flesh: The Biography of Michelangelo’s David. B’Gruppo.
  * Daligault (2016) Daligault, J. 2016 On the quantum Landau collision operator and electron collisions in dense plasmas. Phys. Plasmas 23, 032706.
  * Danielewicz (1980) Danielewicz, P. 1980 Nonrelativistic and relativistic Landau/Fokker-Planck equation for arbitrary statistics. Physica A 100, 167.
  * Eyink (2018) Eyink, G. L. 2018 Cascades and dissipative anomalies in nearly collisionless plasma turbulence. Phys. Rev. X 8, 041020\.
  * Gardner (1963) Gardner, C. S. 1963 Bound on the energy available from a plasma. Phys. Fluids 6, 839.
  * Gibbs (1902) Gibbs, J. W. 1902 Elementary Principles In Statistical Mechanics Developed With Especial Reference To The Rational Foundation Of Thermodynamics. New York: Charles Scribner.
  * Helander (2017) Helander, P. 2017 Available energy and ground states of collisionless plasmas. J. Plasma Phys. 83, 715830401.
  * Helander & Sigmar (2005) Helander, P. & Sigmar, D. J. 2005 Collisional Transport in Magnetised Plasmas. Cambridge: Cambridge University Press.
  * Kadomtsev (1965) Kadomtsev, B. B. 1965 Plasma Turbulence. New York: Academic Press.
  * Kadomtsev & Pogutse (1970) Kadomtsev, B. B. & Pogutse, O. P. 1970 Collisionless relaxation in systems with Coulomb interactions. Phys. Rev. Lett. 25, 1155.
  * Klimontovich (1967) Klimontovich, Y. L. 1967 The Statistical Theory of Non-Equilibrium Processes in a Plasma. MIT Press.
  * Landau (1936) Landau, L. D. 1936 Transport equation in the case of Coulomb interaction. Zh. Eksp. Teor. Fiz. 7, 203.
  * Lenard (1960) Lenard, A. 1960 On Bogoliubov’s kinetic equation for a spatially homogeneous plasma. Ann. Phys 10, 390.
  * Lorenz (1955) Lorenz, E. N. 1955 Available potential energy and the maintenance of the general circulation. Tellus 7, 157.
  * Lynden-Bell (1967) Lynden-Bell, D. 1967 Statistical mechanics of violent relaxation in stellar systems. Mon. Not. R. Astron. Soc. 136, 101\.
  * Maxwell (1860) Maxwell, J. C. 1860 V. Illustrations of the dynamical theory of gases.—Part I. On the motions and collisions of perfectly elastic spheres. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 19, 19.
  * Rightley & Baalrud (2021) Rightley, S. & Baalrud, S. D. 2021 Kinetic model for electron-ion transport in warm dense matter. Phys. Rev. E 103, 063206\.
  * Schekochihin et al. (2009) Schekochihin, A. A., Cowley, S. C., Dorland, W., Hammett, G. W., Howes, G. G., Quataert, E. & Tatsuno, T. 2009 Astrophysical gyrokinetics: kinetic and fluid turbulent cascades in magnetized weakly collisional plasmas. Astrophys. J. Suppl. Ser. 182, 310.
  * Severne & Luwel (1980) Severne, G. & Luwel, M. 1980 Dynamical theory of collisionless relaxation. Astrophys. Space Sci. 72, 293.
  * Shannon (1948) Shannon, C. E. 1948 A mathematical theory of communication. The Bell System Technical Journal 27, 379\.
  * Spitzer (1967) Spitzer, L. 1967 Physics of Fully Ionized Gases. New York: Interscience Publishers.
  * Su & Oberman (1968) Su, C. H. & Oberman, C. 1968 Collisional damping of a plasma echo. Phys. Rev. Lett. 20, 427.
  * Swanson (2008) Swanson, D. G. 2008 Plasma Kinetic Theory. New York: Chapman and Hall.
  * Uehling & Uhlenbeck (1933) Uehling, E. A. & Uhlenbeck, G. E. 1933 Transport phenomena in Einstein-Bose and Fermi-Dirac gases. I. Phys. Rev. 43, 552.
  * Zhdankin (2021) Zhdankin, V. 2021 Generalized entropy production in collisionless plasma flows and turbulence, arXiv: 2110.07025.


Copyright for this paper by its authors. Use permitted under Creative Commons
License Attribution 4.0 International (CC BY 4.0).

Ital-IA 2024: 4th National Conference on Artificial Intelligence, organized by
CINI, May 29-30, 2024, Naples, Italy

<EMAIL_ADDRESS>] [1] [2]
<EMAIL_ADDRESS>]<EMAIL_ADDRESS>]
<EMAIL_ADDRESS>]

[1]Corresponding author. [1]Authors are listed in alphabetical order. [2]Phd
student enrolled in the National PhD in Artificial Intelligence for Health and
Life Sciences, XXXVIII cycle, Università Campus Bio-Medico, Roma.

# Graph Neural Networks for Gut Microbiome Metaomic data: A preliminary work

Christopher Irwin  Flavio Mignone  Stefania Montani  Luigi Portinale  DISIT,
Computer Science Institute, University of Piemonte Orientale, Alessandria,
Italy  SmartSeq s.r.l., Alessandria, Italy  Interdepartmental Research Center
on Artificial Intelligence (AI@UPO), University of Piemonte Orientale, 15121
Alessandria, Italy

(2024)

###### Abstract

The gut microbiome, crucial for human health, presents challenges in analyzing
its complex metaomic data due to high dimensionality and sparsity. Traditional
methods struggle to capture its intricate relationships. We investigate graph
neural networks (GNNs) for this task, aiming to derive meaningful
representations of individual gut microbiomes. Unlike methods relying solely
on taxa abundance, we directly leverage phylogenetic relationships, in order
to obtain a generalized encoder for taxa networks. The representation learnt
from the encoder are then used to train a model for phenotype prediction such
as Inflammatory Bowel Disease (IBD).

###### keywords:

Multiomics Microbiome Graph Representation Learning

##  1 Introduction

The gut microbiome, a community of microorganisms residing within the human
gut, plays a significant role in health and disease. Advancements in
sequencing have yielded vast amounts of metaomic data, detailing the gut
microbiome’s composition and function. However, effectively analyzing this
high-dimensional data remains a challenge. Traditional methods often fail to
capture the intricate relationships between the diverse microbial species
within the microbiome.

This work explores the potential of graph neural networks (GNNs) for analyzing
gut microbiome metaomic data. GNNs, a powerful class of deep learning models
adept at handling graph-structured data, hold promise for unraveling the
complex relationships within the microbiome. Our primary objective is to
develop a method that utilizes metaomic data to obtain a meaningful
representation of a patient’s gut microbiome. In particular, we begin by
constructing a network that captures the relationships between genes, species,
and genus. This network is built using the gene expression levels obtained
from all patient samples. By leveraging GNN techniques, we then learn latent
representations (embeddings) for each entity within the network. These
embeddings capture the underlying functional relationships between genes,
species, and genus. Finally, we aggregate these embeddings based on a
patient’s specific gene expression profile to obtain a global representation
of their unique gut microbiome. This patient-specific microbiome
representation can then be used by a classifier to predict a particular
phenotype, such as the presence or absence of Inflammatory Bowel Disease (IBD)
as done in the present work.

Recently, a wide rage of methods have been proposed to analyse microbiome data
and use it to obtain predictive models for biomarkers discovery and
classifying patients based on some phenotype. Some of these methods rely on
the use of the microbiome features like the relative abundance of taxa.
However, such approaches can be limited by the inherent nature of microbiome
datasets: sparsity and high dimensionality (e.g., the presence of many
features). This may lead to overfitting, with models failing to generalize to
other datasets [1].

In order to address these kind of limitations, many works try to leverage the
intrinsic phylogenetic information between taxa [2]. For example PopPhy-CNN
[3] uses a CNN-based model that works on the 2D matrix derived from the
phylogenetic tree representation. MIOSTONE [1] uses a gated phylogenetic
encoded neural network. TaxoNN [4] leverages an ensemble of CNNs each
specializing to particular taxonomic phylum.

Our method directly builds upon the concept of leveraging phylogenetic
relationships to create a meaningful representation of the microbiome. We
achieve this by employing a GNN-based method to learn a generalized encoder
specifically for the network of taxa. This encoder is trained independently
from the downstream classification task, which is handled by a separate model.

##  2 Dataset

The dataset in use is the IBDMDB which is described in [5] related to a IBD
cohort of patients. The database contains observations coming from 2 different
omic levels: microbiome metagenomics (MGX) and metatranscriptomics (MTX) gene
expression levels (expressed in count per milion - CPMs) along with some other
information related to the patient state and the presence of absence of IBD.

Different omic levels have different number of samples and different number of
features, and in particular:

  * •

Metagenomics: 1635 samples with 108433 features.

  * •

Metatranscriptomics: 736 samples with 70711 features.

Regarding the patient labels, these are related to the presence of absence of
IDB (binary). Notably, there are 1594 patients with labels. The dataset has a
quite substantial skew, since there are 83% negative and 17% positive samples.

The dataset has a large dimensionality and a different rate of missing
samples. To this end, our approach is to use methods which can deal with
missing samples for a particular omic-level and a different number of
features. Moreover, as already expressed in the previous section, when dealing
with metaomics, it can potentially be useful to being able to model an
interaction between the different microorganisms.

Figure 1: Overview of the proposed method: input data characterized by
microbiome gene expression across different omic levels (left). Extracted
graph connecting genes based on species and genus (center). Application of the
GNN module to contextualize node representations: the set of node embeddings
representing a patient are aggregated by taking in to account the relative
abundance expressed in CPMs and fed to the final classifier predicting IBD
presence (bottom)

##  3 Proposed Method

The main task we want to carry out is to obtain a meaningful representation of
a patient using the meta-omic data coming from the analysis of gut microbiome
samples. To do this, the idea is to establish relationships among the various
microorganisms present in the microbiome, whose presence is determined by the
abundance of enzymes/genes that characterize them, using the relational
information of species and genus encoded within a graph. Moreover, since it is
common to have incomplete datasets (either patients may have different gene
expression or missing data may come from a particular omic level), we model
our task and the way of treating the data in a way that allows transparency
from this point of view. Figure 1 gives an overview of the proposed pipeline.

The following subsections will describe how to get the representation of a
patient starting from the metagenomic level. The same technique can be applied
in the same way to other omic levels and seamlessly integrated in this
pipeline.

###  3.1 Graph Construction

Given a set of patients with gene expression levels of gut microbiome
microorganisms, and associated data on enzymes, species, and genus
(phylogenetic information) for each microorganism, our aim is to construct a
graph $\mathcal{G}=(U,V)$ where:

  * •

$U$ denotes the set of nodes, encompassing enzymes, species, and genus.

  * •

$V$ represents the set of edges, categorized into two types: (enzyme, species)
and (species, genus).

Once the phylogenetic graph is constructed, we can represent each patient
$P_{i}$ as a subset of nodes $U_{i}\subseteq U$, where the patient exhibits a
gene expression level greater than 0.

###  3.2 Graph Representation Learning Module

Given a graph $\mathcal{G}$ that encodes the relationships between genes, we
aim to learn a representation (embedding) of the nodes that reflects their
connectivity within the network. To achieve this, we leverage established
methods from the field of graph representation learning, specifically: Graph
Laplacian Eigenvector Positional Encoding, Random Walk Positional Encoding,
and Node2Vec.

Here, we provide a concise description of each method.

Graph Laplacian Eigenvector Positional Encoding (LPE): this technique relies
on the factorization of the graph’s Laplacian matrix. The embedding for a node
is defined by the k-smallest, non-trivial eigenvectors associated with that
node [6].

Random Walk Positional Encoding (RWPE): This approach, as described in [6],
utilizes the random walk matrix to generate node embeddings.

Node2Vec (N2V): This method, presented in [7], is a semi-supervised algorithm
that leverages random walks on the graph to learn meaningful features for the
nodes.

###  3.3 Aggregation Function

Having obtained embeddings for all the nodes in the graph, we now require a
methodology to aggregate a subset of these embeddings into a single vector
that represents a patient. Two aggregation levels are employed for this
purpose.

The first level of aggregation aims to generate a flattened representation for
each gene. Here, a gene’s embedding is calculated as the mean of the
embeddings associated with its corresponding phylogenetic subgraph.

The second level of aggregation focuses on converting the set of genes
expressed by a patient into a single embedding representation. This is
achieved by utilizing a hyperparameter that determines the number of most
highly expressed genes considered for the aggregation process. The optimal
value of this hyperparameter will be investigated in Section 4.3.

###  3.4 Different Omic Levels Integration

Previous sections have outlined the methodology for transforming gene
abundance data (obtained from a metagenomic analysis) of a patient’s
microbiome into a single patient representation. However, incorporating
metatranscriptomic features is straightforward. We simply include genes
identified at the transcript level during the construction of the phylogenetic
graph, these additional nodes are then considered when generating the patient
representation.

##  4 Experimental results

To evaluate the performance of our proposed method, we conducted preliminary
experiments using the IBD dataset described in Section 2. These experiments
aimed to investigate the method’s effectiveness in predicting the presence or
absence of IBD under various parameter settings. We addressed the following
key questions:

  1. Q1.

Node Embedding Comparison: which node embedding technique (LPE, RWPE, N2V)
yields the best performance?

  2. Q2.

Impact of Multi-Omics Integration: how does the model’s performance vary when
using different omic levels as input: metagenomics only versus metagenomics
combined with metatranscriptomics?

  3. Q3.

Gene Selection Strategy: how does the number of genes considered during
patient representation generation impact the final model performance?

For the final classification task utilizing patient representations as input
features, we employed a Support Vector Machine (SVM) classifier using a Radial
Basis Function (RBF). Our choice of a relatively simple model reflects our
primary focus on evaluating the efficacy of the graph representation
techniques.

###  4.1 Evaluation of Node Embedding Techniques (Q1)

To assess the performance of different node embedding techniques (LPE, RWPE,
N2V), we employed an 80/10/10 train-validation-test split. For each technique,
we optimized model hyperparameters using the Optuna framework [8]. Following
hyperparameter tuning, each optimal model configuration was evaluated over
independent runs, utilizing randomized train-validation-test splits and random
initialization. Table 1 reports the results related to Q1. It is clear from
such results that the best choices are LPE and N2V, which perform very
similarly on both datasets. While RWPE significantly underperforms with
repsect to other methods.

###  4.2 Impact of Omic Levels (Q2)

To assess the influence of data derived from different omic levels on model
performance (Q2), we trained two distinct models. The first model used only
data from the MGX level, while the second incorporated data from both
metagenomics (MGX) and metatrascriptomics (MTX) levels. Importantly, the
evaluation for Q2 employed the same experimental setup established for Q1
(detailed in the previous section). The corresponding results for both models
are presented in Table 1. From the results it is possible to notice that the
model trained using only MGX data sligthly underperforms with respect to the
model that includes MTX data.

Table 1: Results obtained by the models using different methods for encoding
the phylogenetic graph

Dataset | F1 score | ROC AUC
---|---|---
LPE | RWPE | N2V | LPE | RWPE | N2V
MGX only | 0.839 | 0.664 | 0.837 | 0.908 | 0.709 | 0.914
MGX+MTX | 0.871 | 0.681 | 0.873 | 0.929 | 0.721 | 0.923

Figure 2: The plot compares the ROC AUC achieved by the model trained on two
different datasets: MGX and MGX+MTX. In this case the graph encoder used the
LPE tecnhique for the node embedding representation.

###  4.3 Gene Selection Strategy (Q3)

Relatively to Q3, we focused on the model trained using the LPE encoding for
the phylogenetic graph. To obtain the representation of a patient, we trained
the model using as features a different number $k$ of genes. In particular, we
considered the top-$k$ most expressed genes. Moreover, we also compared the
results between the two datasets (MGX-only, MGX+MTX). As shown in Figure 2,
the model’s ROC-AUC initially increases with the number of features, peaking
at a certain point before a slight decrease and stabilization. This likely
reflects the importance of highly expressed genes for class separation. Genes
with lower expression levels are likely to carry less information.
Interestingly, despite a small difference, the optimal number of genes differs
between the datasets. The MGX-only model achieves its peak ROC-AUC with 3568
genes, while the MGX+MTX model peaks at 3366. This might suggest that the MTX
data provides additional useful information, resulting in a smaller number of
necessary genes.

##  5 Conclusion and future work

This paper proposes a method leveraging Graph Neural Networks (GNNs) to
extract phylogenetic information from microbiome metagenomics data. The method
aims to create a generalizable encoder for the microbiome, enabling its
application to phenotype classification tasks. The approach achieves good
performance on real-world datasets even with simple classification models,
although it falls short of models specifically designed for classification
tasks [9]. Future work will involve:

  * •

Expanding comparisons with other methods.

  * •

Evaluating performance on additional datasets.

  * •

Exploring the use of clustering for phenotype classification using the learned
encoder.

###### Acknowledgements.

We acknowledge the use of the Chameleon Cloud platform [10] for providing the
essential computational resources that facilitated the execution of our
experiments and model training.

## References

  * Jiang et al. [2024] Y. Jiang, M. Aton, Q. Zhu, Y. Y. Lu,  Miostone: Modeling microbiome-trait associations with taxonomy-adaptive neural networks,  bioRxiv (2024). URL: https://www.biorxiv.org/content/early/2024/02/04/2023.11.04.565596. doi:10.1101/2023.11.04.565596.
  * Hernández Medina et al. [2022] R. Hernández Medina, S. Kutuzova, K. N. Nielsen, J. Johansen, L. H. Hansen, M. Nielsen, S. Rasmussen,  Machine learning and deep learning applications in microbiome research,  ISME Communications 2 (2022) 98. URL: https://doi.org/10.1038/s43705-022-00182-9. doi:10.1038/s43705-022-00182-9.
  * Reiman et al. [2020] D. Reiman, A. A. Metwally, J. Sun, Y. Dai,  Popphy-cnn: A phylogenetic tree embedded architecture for convolutional neural networks to predict host phenotype from metagenomic data,  IEEE Journal of Biomedical and Health Informatics 24 (2020) 2993–3001. doi:10.1109/JBHI.2020.2993761.
  * Sharma et al. [2020] D. Sharma, A. D. Paterson, W. Xu,  TaxoNN: ensemble of neural networks on stratified microbiome data for disease prediction,  Bioinformatics 36 (2020) 4544–4550. URL: https://doi.org/10.1093/bioinformatics/btaa542. doi:10.1093/bioinformatics/btaa542.
  * Lloyd-Price et al. [2019] J. Lloyd-Price, C. Arze, A. N. Ananthakrishnan, M. Schirmer, J. Avila-Pacheco, T. W. Poon, E. Andrews, N. J. Ajami, K. S. Bonham, C. J. Brislawn, et al.,  Multi-omics of the gut microbial ecosystem in inflammatory bowel diseases,  Nature 569 (2019) 655–662.
  * Dwivedi et al. [2021] V. P. Dwivedi, A. T. Luu, T. Laurent, Y. Bengio, X. Bresson,  Graph neural networks with learnable structural and positional representations,  arXiv preprint arXiv:2110.07875 (2021).
  * Grover and Leskovec [2016] A. Grover, J. Leskovec,  node2vec: Scalable feature learning for networks,  in: Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD ’16, Association for Computing Machinery, New York, NY, USA, 2016, p. 855–864. URL: https://doi.org/10.1145/2939672.2939754. doi:10.1145/2939672.2939754.
  * Akiba et al. [2019] T. Akiba, S. Sano, T. Yanase, T. Ohta, M. Koyama,  Optuna: A next-generation hyperparameter optimization framework,  in: Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, 2019.
  * Hu et al. [2023] M. Hu, J. Zhu, G. Peng, W. Lu, H. Wang, Z. Xie,  IMOVNN: incomplete multi-omics data integration variational neural networks for gut microbiome disease prediction and biomarker identification,  Briefings in Bioinformatics 24 (2023) bbad394. URL: https://doi.org/10.1093/bib/bbad394. doi:10.1093/bib/bbad394.
  * Keahey et al. [2020] K. Keahey, J. Anderson, Z. Zhen, P. Riteau, P. Ruth, D. Stanzione, M. Cevik, J. Colleran, H. S. Gunawi, C. Hammock, J. Mambretti, A. Barnes, F. Halbach, A. Rocha, J. Stubbs,  Lessons learned from the chameleon testbed,  in: Proceedings of the 2020 USENIX Annual Technical Conference (USENIX ATC ’20), USENIX Association, 2020.