# Parallel expansion of a fuel pellet plasmoid

Alistair M. Arnold1<EMAIL_ADDRESS>Pavel Aleynikov1  Boris N.
Breizman2  1Max-Planck-Institut für Plasmaphysik, Greifswald, Germany
2Institute for Fusion Studies, University of Texas at Austin, Austin, USA

## Abstract

The problem of the expansion and assimilation of a cryogenic fuel pellet
injected into a hot plasma is considered. Due to the transparency of the
plasmoid to ambient particles, it is found that electrons reach a ‘quasi-
equilibrium’ (QE) which is characterised by a steady-state on the fastest
collisional timescale. The simplified electron kinetic equation of the quasi-
equilibrium state is solved. Taking a velocity moment of the electron kinetic
equation permits a fluid closure, yielding an evolution equation for the
parameters describing the QE distribution function. In contrast to the
Braginskii equations, the closure does not require that electrons have a short
mean free path compared to the size of density perturbations and permits an
anisotropic and highly non-Maxwellian distribution function. Since the QE
electron distribution function accounts for both trapped and passing
electrons, the self-consistent electric potential that causes the expansion
can be properly described, in contrast to earlier models of pellet plasmoid
expansion with an unbounded potential. The plasmoid expansion is simulated
using both a Vlasov model and a cold fluid model for the ions. During the
expansion plasmoid ions and electrons obtain a nearly equal amount of energy;
as hot ambient electrons provide this energy in the form of collisional
heating of plasmoid electrons, the expansion of a pellet plasmoid is expected
to be a potent mechanism for the transfer of energy from electrons to ions on
a timescale shorter than that of ion-electron thermalisation.

##  1 Introduction

During a recent experimental campaign of the W7-X stellarator, fuel pellet
injection was found to be associated with a large transfer of energy from
electrons to ions (Baldzuhn et al.,, 2019, 2020; Bozhenkov et al.,, 2020).
Such an energy transfer is generally desirable, since it acts to bring the ion
temperature up to the electron temperature (the ion temperature being lower
than the electron temperature during normal operation), and a higher ion
temperature results in a larger fusion cross-section. Subsequent investigation
of the dynamics of the injection and assimilation of fuel pellets have
suggested that a possible mechanism for the energy transfer is the rapid
ambipolar expansion of the ionised pellet material – the pellet plasmoid –
along magnetic field lines (Aleynikov et al.,, 2019; Runov et al.,, 2021;
Arnold et al.,, 2021). The aim of this paper is to resolve the inconsistencies
present in these models and provide a rigorous model of the parallel plasmoid
expansion. What follows is a brief recapitulation of the processes by which
the pellet plasmoid is formed, the reason for its parallel expansion and
concomitant electron-ion energy transfer, and a summary of the approaches and
pitfalls of earlier models. An outline of a new approach which does not suffer
from these pitfalls is provided before being realised mathematically.

When a fuel pellet is injected into an MCF (Magnetic Confinement Fusion)
device, the incoming energy flux from the ambient plasma ablates the surface
of the pellet and produces a gas cloud (Parks et al.,, 1977). The pellet and
gas are composed of electrically neutral molecules, but plasma is continuously
generated within the gas cloud by the collisions of the high-energy ions and
electrons composing the multi $\mathrm{keV}$ ambient plasma with gas
molecules. Subsequently, the pellet and gas cloud continue to cross magnetic
field lines at the speed at which they were injected, but some of the newly-
ionised plasma is left behind; that which was not collisionally ‘dragged’
along with the moving gas cloud. This is because the plasma constituents are
charged particles and follow Larmor orbits that ‘pin’ the particles to the
field line. The result is that a plasmoid, a localised excess density of
plasma, is deposited on field lines that intersected the gas cloud as it
traversed the device.

Since the plasmoid is a localised density perturbation and electrons have a
much higher thermal velocity than ions, the electric potential required to
maintain quasineutrality acts to trap electrons inside the plasmoid and
accelerate ions away from the plasmoid. Figure 1 shows a schematic of the
plasmoid and electric potential. Since the potential acts to trap electrons,
we will use the names ‘well’ and ‘potential’ interchangeably. With regards to
pellet plasmoids the density is such that the the electric potential drives
parallel dynamics much more quickly than transverse dynamics occur, the latter
being due to drifts. Therefore, as with previous investigations, we consider
only the parallel expansion of the pellet plasmoid on one given field line.

We stress that not all of the plasma produced from the pellet ablatant fits
the description of the previous paragraph. Naturally, plasma that is ‘dragged
along’ with the gas cloud exhibits quite different dynamics. However, for
‘fast’ pellet injection devices proportionally more of the plasmoid dynamics
occur on field lines where the gas cloud and pellet have departed (Arnold et
al.,, 2021). The injection devices fitting the criteria of being ‘fast’ are
becoming the norm in MCF experiments, so the conclusions drawn from the
parallel plasmoid expansion in the absence of gas can be expected to
reasonably well apply to pellet injection in future MCF devices.

The dynamics of any plasmoid immersed in an ambient plasma depend greatly upon
the plasma and plasmoid parameters, such as their relative temperatures,
densities, the plasmoid size, and so on. Naturally, it is difficult to
describe plasmoid dynamics with a too wide-ranging choice of parameters, so
our attention must be restricted to plasmoids broadly corresponding to a those
produced by pellet injection in a state-of-the-art MCF device. We take as a
reference point the W7-X stellarator, since the success of its pellet
injection campaign provides motivation for studying pellet plasmoids. Further,
the temperatures and densities in the core of W7-X are generally comparable to
other high-performance MCF devices.

For the purpose of untangling the different phenomena involved in plasmoid
expansion it is helpful to provide concrete plasma parameters We consider an
ambient plasma of electron density $n_{a}=5\times 10^{19}\,\mathrm{m}^{-3}$ at
a temperature $T_{a}=5\,\mathrm{keV}$. A typical line-integrated density along
the field line of a fuel pellet plasmoid in W7-X is
$N_{p}=10^{22}\,\mathrm{m}^{-2}$ (Arnold et al.,, 2021).

The fuel pellets in W7-X contain approximately $10^{20}$ electrons and
penetrate roughly $0.1\,\mathrm{m}$ into the plasma, resulting in the average
line-density in the radial direction of $10^{21}\,\mathrm{m}^{-1}$ (Baldzuhn
et al.,, 2019). In W7-X the flux surface located at minor radius
$r=0.3\,\mathrm{m}$, given the major radius $R=5\,\mathrm{m}$, has a flux-
surface integrated density of $3\times 10^{21}\,\mathrm{m}^{-1}$ if the
density of this flux surface is $5\times 10^{19}\,\mathrm{m}^{-3}$. Hence for
such a flux surface the temperature is not strongly quenched after
assimilation. For high-performance scenarios in W7-X the quenching is even
weaker due to the higher plasma densities. Therefore, unlike killer pellets,
which quench the temperature completely, on large flux surfaces fuelling
pellets only slightly affect the temperature. We therefore neglect any change
in $T_{a}$ during the plasmoid expansion.

We consider irrational flux surfaces where individual field lines have a
connection length $L_{F}$ that is, in principle, infinite. However, given that
the plasmoid has a transverse size $r_{I}$, the connection length of the flux
tube containing the plasmoid is in practice $(2\pi R)(2\pi r)/r_{I}$, since
this is the length after which the flux tube of diameter $r_{I}$ self-
intersects. For W7-X pellets $r_{I}\approx 0.1\,\mathrm{m}$ shortly after
injection (Baldzuhn et al.,, 2019; Arnold et al.,, 2021), giving a maximum
connection length of $600\,\mathrm{m}$. Since the plasmoid is expected to
reach this size after its density has dropped to practically the value of the
ambient plasma for $N_{p}=10^{22}\,\mathrm{m}^{-2}$ (Aleynikov et al.,, 2019),
we formally take the connection length to be infinite. Arnold et al., (2023)
in contrast treated the electron kinetics in a high-Z plasmoid on a field line
of finite connection length, accounting for the quenching of the ambient
plasma.

Since plasmoid electrons are ‘born’ at energies comparable to the ionisation
energy, of order tens of $\mathrm{eV}$, but are immersed in an ambient plasma
with a temperature on the order of several $\mathrm{keV}$, the electron
distribution function as a whole will consist of a cold, dense core of
plasmoid electrons and a hot, sparse tail of ambient electrons. The
distribution function is only close to a Maxwellian after the plasmoid
electrons have been sufficiently heated by the ambient electrons, which
happens after the plasmoid has significantly expanded with the plasma
parameters we use here. The primary concern with previous models of the
expansion is that they treated only the cold plasmoid electrons, assuming that
they have a near-Maxwellian distribution function, but did not treat ambient
electrons. These electrons were simply assumed to be of a constant density,
hence providing collisional heating to the plasmoid electrons. Further,
ambient ions were not considered at all in the fluid model for the ions. The
consequence of this approach is that the electric potential decreases without
bound as the plasmoid density vanishes. Clearly, one cannot use this approach
as a basis for treating both trapped and passing electrons. The fact that the
electric potential was unphysical also called into question the result of the
electron-ion energy transfer and other aspects of the expansion.

A more sophisticated approach to electrons is required to resolve these
issues. We will consider electron kinetics in the variables of parallel energy
$\mathcal{E}_{\parallel}=m_{e}v_{\parallel}^{2}/2-e\phi(z)$, where
$v_{\parallel}$ is the velocity parallel to the field line, perpendicular
energy $\mathcal{E}_{\perp}=m_{e}v_{\perp}^{2}/2$, where $v_{\perp}$ is the
speed perpendicular to the field lines, $z$, the position along the field
line, and $t$, time. In anticipation of the form of the distribution function
for different energies we split the phase space into regions I, II, and III,
respectively corresponding to the deeply trapped electrons, hot trapped
electrons, and hot passing electrons (Fig. 2).

In each region we employ the separation of timescales appropriate for the
plasmoid and plasma parameters mentioned earlier in this section in order to
obtain a simplified kinetic equations for the electrons. We find that trapped
electrons collide with the cold, dense plasmoid electrons (and the plasmoid
ions) much more quickly than with the passing electrons. At the same time,
owing to the high temperature of the ambient plasma, the mean free path of
passing and hot trapped electrons exceeds the length of the plasmoid; the
plasmoid appears essentially transparent to the ambient electrons, and hot
trapped electrons bounce inside the well many times before colliding. The
latter effect means that trapped electrons behave ‘adiabatically’ as the
potential well expands.

We will show that, except at the very earliest stage of expansion, the
ordering of timescales leads to the electrons reaching a ‘quasi-equilibrium’
(QE) state which is characterised by rapid electron collisions within the
plasmoid causing the electron distribution to exhibit a steady-state on the
timescale on which trapped electrons collide with the plasmoid. The steady-
state is established with no net flux of electrons into the trapped region of
phase-space to prevent the ‘charging up’ of the plasmoid (and hence the
violation of quasineutrality) on this timescale. The QE electron distribution
function is analogous to an equilibrium distribution function, which is indeed
attained if electron-electron collisions are the fastest effect. In our case,
however, the bounce period of hot trapped electrons is considerably shorter
than the collision timescale. We note that the QE state and the equilibrium
state (which has a Maxwellian energy distribution) differ conceptually; a
Maxwellian distribution exhibits no collisional flux, but the QE state is
characterised by a vanishing divergence of collisional flux; in this sense QE
is a ‘dynamical’ steady-state.

It will be shown that the QE distribution is specified in terms of the ‘deeply
trapped’ distribution function occupying region I, a Maxwellian with
homogeneous temperature $T$, which is uniquely defined by two parameters.
These two parameters must be such that there is no net flux of electrons into
the trapped region on the timescales on which QE is established. This allows
us to express one parameter of the lowest-order distribution in terms of the
other. Once the electron distribution is known in terms of this remaining
parameter, which we choose to be the temperature, its zeroth moment may be
taken to obtain an expression for electron density, which, combined with the
quasineutrality condition provides an implicit expression for the self-
consistent electric potential $\phi$ in terms of the temperature. The velocity
moment corresponding to line-integrated energy density is then taken over the
electron kinetic equation to obtain an energy conservation law, which is
practically used as the evolution equation for $T$.

A description of the expansion requires a model for ion motion. Two models
were considered: a cold-fluid model with a single flow velocity and a
collisionless kinetic (Vlasov) model. The first model is pragmatically
justified by a possible application of this work being to provide a simplified
model for pellet plasmoid expansion in an established fluid code. The second
is justified by the long mean-free-path of hot ambient ions. These models
represent opposite collisionality regimes for ions, and we therefore expect
the shared qualitative properties to remain in a more sophisticated and
accurate model for the ions. The qualitative property of greatest concern is
the electron to ion energy transfer during the expansion.

With each ion model the system was evolved until the plasmoid and ambient
densities were similar and the plasmoid electron temperature $T$ had reached
$T_{a}$; the plasmoid assimilated with the ambient plasma. After this point
the electric field does not provide much energy to the ions, so the energy
transfer from electrons to ions may be considered complete. With the given
choice of plasma parameters the densities and temperatures equilibrate at
approximately the same time. With a larger line-integrated density $N_{p}$ the
temperature equilibration would occur well before the densities are
comparable. With a much smaller line-integrated density the densities become
similar well before the temperatures have equilibrated. More discussion of how
the QE formalism fits into the larger topic of plasmoid expansion is given in
a later section.

Late in the plasmoid expansion with the cold-fluid model for ions, a
steepening of the density profile results from the fast plasmoid moving into
the ambient plasma causing a shock near the extremities of the plasmoid. This
shock which will cause the generation of sound waves and solitons that will
propagate into the ambient plasma. The wave propagation will sap the kinetic
energy of the plasmoid and possibly the electron-ion energy balance. However,
the shock and wave dynamics may only be properly accounted for with Poisson’s
equation for the electric potential since a deviation from quasineutrality
will occur near the shock. As we neglect any deviation from quasineutrality,
sound wave and soliton generation is suppressed.

In the Vlasov ion model, which accounts for the ambient ion temperature, no
shock is observed; the density profile smoothly decreases to the ambient
density. This is because many ambient ions can now traverse the entire
plasmoid or are rapidly reflected from the potential, hence do not pile up at
the expanding edge of the plasmoid. It is expected that as the ion temperature
decreases the system is more prone to forming a shock.

###  1.1 Self-similar solution to plasmoid expansion

Aleynikov et al., (2019) provided the self-similar solution to plasmoid
expansion given that the plasmoid is transparent to the ambient plasma.
Although we seek to rectify the issues in the model therein, the density and
temperature profiles obtained with the plasma parameters given earlier will be
used to justify the ordering which will be used to simplify the electron
kinetic problem. We require only order-of-magnitude estimates to find this
ordering, so we deem the self-similar solution good enough for this purpose.
However, as we wish to model the long-term expansion of the plasmoid, well
past the applicability of the self-similar model, we must modify the profiles
in Aleynikov et al., (2019) so that they are valid (i.e. giving a correct
order of magnitude estimate) for the long-term expansion.

Firstly, we note that in the self-similar expansion the plasmoid electron
temperature is given by $T=\nu_{h}T_{a}t$ for $t\ll\nu_{h}^{-1}$, where

$\nu_{h}=\frac{n_{a}e^{4}\ln\Lambda}{6\sqrt{2}\pi^{\frac{3}{2}}\varepsilon_{0}^{2}m_{e}^{\frac{1}{2}}T_{a}^{\frac{3}{2}}}$
(1)

is the inverse heating time of the cold plasmoid electrons by a hot population
of density $n_{a}$ and temperature $T_{a}$. The expression

$T\sim T_{a}\left(1-\mathrm{e}^{-\nu_{h}t}\right)$ (2)

agrees with this linearly increasing temperature at early times, but
exponentially approaches $T_{a}$ as time advances, which is the characteristic
behaviour of a cold Maxwellian being heated by a hotter one. Therefore we
expect the above expression to be adequate in describing the plasmoid electron
temperature in both the short and long-term evolution of the plasmoid.

In the self-similar solution, the density becomes infinite as $t\rightarrow 0$
and vanishes for $t\rightarrow\infty$, neglecting the fact that the electron
density approaches $n_{a}$ as time advances. Therefore the expression

$n_{m}\sim N_{p}\nu_{h}\sqrt{\frac{3m_{i}}{8\pi(\nu_{h}t)^{3}T_{a}}}+n_{a},$
(3)

for the peak plasma density, which simply adds the ambient electron density
$n_{a}$ to the self-similar solution, is a plausible expression for both long
and short term evolution.

Although the electric potential in Aleynikov et al., (2019) diverges as
$|z|\rightarrow\infty$, the Boltzmann relation

$e\phi_{m}\sim T\ln\left(\frac{n_{m}}{n_{a}}\right)$ (4)

provides a good estimate for the height of the potential. This estimate is
also supported by solution to the self-consistent electron kinetic problem in
Arnold et al., (2023), which showed that the potential height is the same
order of magnitude as that suggested by the Boltzmann relation, its exact
value being somewhat larger when $T\ll T_{a}$.

Equation (3) expresses is the peak plasma density, which we will subsequently
use in the ordering. Of course, electrons move throughout the plasmoid, more
for passing electrons and less for trapped electrons, so the most rigorous
approach would be to consider ‘average’ quantities throughout the orbit. This
would be very complicated, and relies on a detailed knowledge of the shape of
the potential which we have not yet obtained. Therefore, we apply the ordering
using the quantities at the peak of the plasmoid with the reasoning that the
plasmoid density is very large at its peak and decreases rapidly as one moves
away from it. Hence, when considering trapped electrons the orbit-average of
any quantity will be heavily weighted by the value at the peak. When consider
passing electrons, we will actually obtain the same expressions as those
obtained by rigorous consideration.

Figure 1: Schematic of the electric potential induced by the presence of the
plasmoid (Arnold et al.,, 2023). Example trapped (with turning points $\pm
z_{c}$) and passing electron trajectories are included. The profiles are
assumed to be even and monotonically decreasing in $z$ with the electron
density and potential reaching their maxima $n_{m}$ and $\phi_{m}$ at $z=0$.
Figure 2: Schematic of the phase-space domain of the electron kinetic problem
at $z=0$. This is also the domain for the bounce-averaged kinetic problems.
The dotted line indicates
$\mathcal{E}=\mathcal{E}_{\parallel}+\mathcal{E}_{\perp}=0$. The diagonal
dashed line indicates $\mathcal{E}=\mathcal{E}_{\mathrm{I/II}}$, which
separates regions I and II. The vertical dashed line indicates
$\mathcal{E}_{\parallel}=0$, the trapped-passing separatrix.

##  2 Electron kinetics

The kinetic equation for the electron distribution function $f$ is given by

$\frac{\partial f}{\partial t}+v_{\parallel}\frac{\partial f}{\partial
z}+\frac{e}{m_{e}}\frac{\partial\phi}{\partial z}\frac{\partial f}{\partial
v_{\parallel}}=C(f,f)+\sum_{k}C_{e,ik}(f)$ (5)

in the variables $(v_{\parallel},v_{\perp},z,t)$, where $v_{\parallel}$ is the
velocity parallel to the field line, $v_{\perp}$ is the speed perpendicular to
the field line, $z$ is the coordinate along the field line, $t$ is time,
$C(f,f)$ is the electron self-collision operator and $C_{e,ik}(f)$ is the
collision operator for electrons colliding against ion species $k$.

We change to the independent variables
$(\mathcal{E}_{\parallel}=m_{e}v_{\parallel}^{2}/2-e\phi,\mathcal{E}_{\perp}=m_{e}v_{\perp}^{2}/2,z,t)$:

$\frac{\partial f}{\partial t}+v_{\parallel}\frac{\partial f}{\partial
z}-e\frac{\partial\phi}{\partial t}\frac{\partial
f}{\partial\mathcal{E}_{\parallel}}=C(f,f)+\sum_{k}C_{e,ik}(f),$ (6)

showing that collisionless change in electron energy is associated with time-
variation of the electric potential. We note that with a stationary potential
both $\mathcal{E}_{\parallel}$ and $\mathcal{E}_{\perp}$ are constants of
motion in the absence of collisions. Passing electrons have
$\mathcal{E}_{\parallel}\geq 0$ and trapped $\mathcal{E}_{\parallel}<0$. The
minimum parallel energy of an electron is
$\mathcal{E}_{\parallel}=-e\phi_{m}$, which corresponds to an electron with
$v_{\parallel}=0$ which remains at $z=0$.

We now solve the kinetic equations in regions I, II, and III. We split up the
distribution function into its representation in each region,
$f_{\mathrm{I}}$, $f_{\mathrm{II}}$, and $f_{\mathrm{III}}$. That is, when we
are in region I, for example, we solve the kinetic equation for
$f_{\mathrm{I}}$. No matter where we are in phase-space, collisions are
experienced with every other region of phase-space via
$C(\cdot,f_{\mathrm{I}}+f_{\mathrm{II}}+f_{\mathrm{III}})$. In this sense it
is understood that $f_{\mathrm{I}}$ (or $f_{\mathrm{II}}$ or
$f_{\mathrm{III}}$) is the entire distribution function in region I (or II or
III), but is zero outside of its own region.

###  2.1 Solving the kinetic equation for passing electrons (region III)

The kinetic equation for electrons in region III is given by

$\frac{\partial f_{\mathrm{III}}}{\partial t}+v_{\parallel}\frac{\partial
f_{\mathrm{III}}}{\partial z}-e\frac{\partial\phi}{\partial t}\frac{\partial
f_{\mathrm{III}}}{\partial\mathcal{E}_{\parallel}}=C(f_{\mathrm{III}},f)+\sum_{k}C_{e,ik}(f_{\mathrm{III}}).$
(7)

The collision frequency of an electron inside the plasmoid with the plasmoid
electrons and ions is approximately given by

$\nu_{p}(v)=\frac{n_{p}(1+Z_{\mathrm{eff}})e^{4}\ln\Lambda}{8\pi\varepsilon_{0}^{2}m_{e}^{2}v^{3}},$
(8)

where $v$ is the electron speed, $n_{p}$ is the plasmoid electron density, and

$Z_{\mathrm{eff}}=\frac{\sum_{k}n_{ik}Z_{k}^{2}}{\sum_{k}n_{ik}Z_{k}}$ (9)

is the effective charge of the ions. Equation (8) is derived from the
frequency with which an electron experiences pitch-angle scattering, assuming
that quasineutrality

$n_{e}\approx n_{p}=\sum_{k}n_{ik}Z_{k},$ (10)

where $n_{e}$ is the electron density, holds inside the plasmoid. The typical
velocity of a passing electron is given by the ambient electron thermal
velocity $v_{T_{a}}=\sqrt{2T_{a}/m_{e}}$. The inverse of the time taken for a
passing electron to transit the plasmoid is $\nu_{T}=v_{T_{a}}/L_{p}$ for
plasmoid length $L_{p}=N_{p}/n_{p}(z=0)$. Hence, for a hydrogenic plasma, in
region III the ratio of the collision frequency with the plasmoid and the
inverse transit time is given by

$\frac{\nu_{p}(v_{T_{a}})}{\nu_{T}}:=\mu=\frac{N_{p}e^{4}\ln\Lambda}{16\pi\varepsilon_{0}^{2}T_{a}^{2}},.$
(11)

With the plasma and plasmoid parameters given in the Introduction this ratio
is much less than unity:

$\mu\ll 1.$ (12)

where for the purpose of calculation the Coulomb logarithm $\ln\Lambda$ is
assumed to be equal to $15$. This implies that the mean free path of passing
electrons is much longer than the plasmoid, so the plasmoid appears
essentially transparent to ambient electrons. We therefore refer to $\mu$ as
the opacity of the plasmoid. $\mu$ is independent of any parameters that
change during the expansion, hence it is small throughout the entire
expansion.

The terms in the kinetic equation (7) containing time derivatives correspond
to collisionless changes in the energy of passing electrons, which certainly
occur on a longer timescale than the transit time. Therefore the shortest
timescale in the region III is the transit time, which is associated with the
convective term in the kinetic equation; the lowest-order kinetic equation for
passing electrons is

$v_{\parallel}\frac{\partial f_{\mathrm{III}}}{\partial z}=0,$ (13)

which implies that $f_{\mathrm{III}}$ is independent of $z$. The kinetic
equation for passing electrons (7) then reduces to

$\frac{\partial f_{\mathrm{III}}}{\partial t}-e\frac{\partial\phi}{\partial
t}\frac{\partial
f_{\mathrm{III}}}{\partial\mathcal{E}_{\parallel}}=C(f_{\mathrm{III}},f)+\sum_{k}C_{e,ik}(f_{\mathrm{III}}),$
(14)

the solution to which can be obtained immediately by bounce-averaging. We
define the bounce integral of a function
$g(\mathcal{E}_{\parallel},\mathcal{E}_{\perp},z,t)$ to be

$\oint g\,\mathrm{d}z:=2\int^{z_{c}}_{-z_{c}}g\,\mathrm{d}z,$ (15)

where $z_{c}(\mathcal{E}_{\parallel},t)>0$ is the turning point such that

$\mathcal{E}_{\parallel}+e\phi(z_{c},t)=0$ (16)

(cf. Fig. 1). The bounce average of $g$ is given by

$\left<g\right>=\frac{1}{\tau}\oint g\,\mathrm{d}z$ (17)

for bounce period

$\tau=\oint\frac{\mathrm{d}z}{v_{\parallel}}.$ (18)

We note that on an infinitely long magnetic field line the bounce-average of
any function $g$ for $\mathcal{E}_{\parallel}>0$ is given by

$\left<g\right>\big{|}_{\mathcal{E}_{\parallel}>0}=\lim_{|z|\rightarrow\infty}g.$
(19)

Since the potential vanishes at infinity and we expect the distribution
function to be constant at infinity, the bounce average of Eq (14) is solved
by the Maxwellian defining the ambient plasma:

$f_{\mathrm{III}}=f_{a}=n_{a}\left(\frac{m_{e}}{2\pi
T_{a}}\right)^{\frac{3}{2}}\mathrm{e}^{-\frac{\mathcal{E}}{T_{a}}}.$ (20)

We note that, as mentioned earlier, we have used the expression for peak
plasmoid density to deduce the ordering (12) and obtain Eq. (20), despite the
fact that plasmoid electrons spend little time inside the plasmoid. The same
result can be obtained by a slightly different approach: Eq. (11) is identical
to the expression for the number of passing electrons slowed down by
collisions with the plasmoid as they traverse it (Arnold et al.,, 2021). If
the number of passing electrons slowed by the plasmoid is small, then the
passing electron distribution is modified little by collisions and is simply
replenished from flux emerging from $|z|\rightarrow\infty$, resulting in the
distribution identical to that at infnity; hence Eq. (20).

###  2.2 Solving the kinetic equation for deeply trapped electrons (region I)

The kinetic equation for electrons in region I is given by

$\frac{\partial f_{\mathrm{I}}}{\partial t}+v_{\parallel}\frac{\partial
f_{\mathrm{I}}}{\partial z}-e\frac{\partial\phi}{\partial t}\frac{\partial
f_{\mathrm{I}}}{\partial\mathcal{E}_{\parallel}}=C(f_{\mathrm{I}},f_{\mathrm{I}})+C(f_{\mathrm{I}},f_{\mathrm{II}}+f_{\mathrm{III}})+\sum_{k}C_{e,ik}(f_{\mathrm{I}}).$
(21)

We expect the potential well to be parabolic at its peak; deeply trapped
electrons bounce inside this parabola and collide with the plasmoid when its
density is near its peak. For a parabolic potential well of height $\phi_{m}$
and width $L_{p}$, the bounce frequency is given by

$\nu_{B}\sim\frac{L_{p}}{\sqrt{2e\phi_{m}/m_{e}}},$ (22)

and is associated with the convective term, so we write

$v_{\parallel}\frac{\partial f_{\mathrm{I}}}{\partial
z}\sim\nu_{B}f_{\mathrm{I}}.$ (23)

Substituting temperature (2) and density (3) from the modified self-similar
solution into the Boltzmann relation (4) yields

$\frac{e\phi_{m}}{T_{a}}\sim(1-\mathrm{e}^{-\nu_{h}t})\ln\left(1+\frac{2}{\pi\sqrt{3}}\sqrt{\frac{m_{i}}{m_{e}}}\mu(\nu_{h}t)^{-\frac{3}{2}}\right).$
(24)

Given that $(2/(\pi\sqrt{3}))\sqrt{m_{i}/m_{e}}\mu$ is at least order unity,
the above is order unity for several heating times $\nu_{h}^{-1}$, with the
exception of the very early stage of the expansion. This is the case for the
plasma parameters given in the Introduction. We also note that while $T<T_{a}$
the Boltzmann relation provides something of an underestimate, strengthening
the argument that $e\phi_{m}\sim T_{a}$ for smaller values of $\mu$.

Owing to the height of the potential, the bounce frequency in region I is of
the same order as the transit frequency of region III:

$\nu_{B}\sim\frac{L_{p}}{\sqrt{2T_{a}/m_{e}}}.$ (25)

Since we expect $f_{\mathrm{I}}$ to correspond to a dense population of cold
electrons, we associate the collision terms against $f_{\mathrm{I}}$ with the
frequency of collisions with the plasmoid; we write

$C(f_{\mathrm{I}},f_{\mathrm{I}})\sim\sum_{k}C_{e,ik}(f_{\mathrm{I}})\sim\nu_{p}\left(v_{T}\right)f_{\mathrm{I}},$
(26)

where $v_{T}=\sqrt{2T/m_{e}}$ the typical electron speed within region I.

Since $f_{\mathrm{II}}$ and $f_{\mathrm{III}}$ represent the hot tail, we
associate the collision terms against these with the heating rate. Since $T$
approaches exponentially $T_{a}$ as time advances, this heating rate decreases
exponentially as time advances: we define the heating rate to be

$\frac{1}{T_{a}}\frac{\mathrm{d}T}{\mathrm{d}t}=\nu_{h}\mathrm{e}^{-\nu_{h}t}$
(27)

and write

$C(f_{\mathrm{I}},f_{\mathrm{II}}+f_{\mathrm{III}})\sim\nu_{h}\mathrm{e}^{-\nu_{h}t}f_{\mathrm{I}}.$
(28)

Comparing the frequency of collisions with the plasmoid with the heating using
the modified self-similar temperature (2) and density (3) gives

$\frac{\nu_{h}\mathrm{e}^{-\nu_{h}t}}{\nu_{p}(v_{T})}\sim\frac{4}{3\sqrt{\pi}}\mathrm{e}^{-\nu_{h}t}\left(1-\mathrm{e}^{-\nu_{h}t}\right)^{\frac{3}{2}}\left(1+\frac{2}{\pi\sqrt{3}}\sqrt{\frac{m_{i}}{m_{e}}}\mu(\nu_{h}t)^{-\frac{3}{2}}\right)^{-1}.$
(29)

With the plasma parameters used in the Introduction the above is much smaller
than unity for all times: collisions with the plasmoid occur much more
frequently than collisions that cause heating. This effect is also enhanced by
the fact that in a parabolic potential well we expect the heating rate to be
slightly reduced, due to the reduction in the density of passing electrons and
their reduction in collisionality (Arnold et al.,, 2023). Hence we can write

$\frac{\nu_{h}\mathrm{e}^{-\nu_{h}t}}{\nu_{p}(v_{T})}\ll 1$ (30)

Since region I represents the cold electrons, we expect the time-dependent
terms in Eq. (21) to act on a much longer timescale than the collision time.
In region III we noted that the transit frequency greatly exceeds the
collision frequency with the plasmoid. However, as we move into region I,
which contains lower-velocity electrons, the bounce frequency (which is
comparable to the transit frequency, cf. Eq. (25)), remains the same as the
collision frequency increases. Therefore, collisions with the plasmoid and
bounce motion are associated with the two shortest timescales in region I.
Accordingly, the lowest-order kinetic equation in region I is

$v_{\parallel}\frac{\partial f_{\mathrm{I}}}{\partial
z}=C(f_{\mathrm{I}},f_{\mathrm{I}})+\sum_{k}C_{e,ik}(f_{\mathrm{I}}).$ (31)

When $T\ll T_{a}$, $T\ll e\phi_{m}$, which allows us to define
$\mathcal{E}_{\mathrm{I/II}}$ such that region I extends for several times $T$
in both the parallel and perpendicular direction. Therefore, the solution to
the above, assuming that collisions with ions are well-approximated by pitch-
angle scattering, is a Maxwellian in energy

$f_{\mathrm{I}}=f_{0}=\eta\left(\frac{m_{e}}{2\pi
T}\right)\mathrm{e}^{-\frac{\mathcal{E}}{T}}$ (32)

for parameters $\eta(t)$, $T(t)$ (Aleynikov et al.,, 2019). We see now that in
the earlier work Aleynikov et al., (2019); Runov et al., (2021); Arnold et
al., (2021) only electrons in region I, where a purely Maxwellian electron
distribution function is exhibited, are treated, whereas in this investigation
we continue the analysis in regions II and III.

The remaining kinetic equation in region I, corresponding to the heating
timescale, is given by

$\frac{\partial f_{\mathrm{I}}}{\partial t}-e\frac{\partial\phi}{\partial
t}\frac{\partial
f_{\mathrm{I}}}{\partial\mathcal{E}_{\parallel}}=C(f_{\mathrm{I}},f_{\mathrm{II}}+f_{\mathrm{III}}),$
(33)

which captures the collisionless change in electron energy due to the
expanding well and the heating of the cold Maxwellian by the hot electrons.

###  2.3 Choosing $\mathcal{E}_{\mathrm{I/II}}$

Owing to the ordering, when $T\ll T_{a}$ we understand that the potential well
is deep enough for a Maxwellian to reside in region I, provided
$\mathcal{E}_{\mathrm{I/II}}$ is chosen close enough to zero. We now decide
upon an explicit definition for $\mathcal{E}_{\mathrm{I/II}}$ which is
consistent with the distribution in region I: we require that collisions with
the cut-off Maxwellian in region I are well-approximated by collisions with
the full Maxwellian that extends to arbitrarily large energies. This will
allow the linearisation of the kinetic problem in region II in terms of
collisions with the Maxwellian. However, we must not artificially extend
region I past the point where the distribution function would be Maxwellian;
this would result in an incorrect distribution function.

From the Boltzmann relation (4) and the density of the full Maxwellian
distribution (32) being $n_{p}=\eta\exp(e\phi/T)$ we find that

$n_{p}\sim n_{a}\mathrm{e}^{\frac{e\phi}{T}}.$ (34)

Considering an electron with with parallel energy
$\mathcal{E}_{\mathrm{I/II}}$, it has a turning point $z_{c}$ such that
$e\phi(z_{c})+\mathcal{E}_{\mathrm{I/II}}=0$. Therefore, at this turning
point,

$n_{p}(z_{c})\sim n_{a}\mathrm{e}^{-\frac{\mathcal{E}_{\mathrm{I/II}}}{T}}.$
(35)

During its orbit, this electron collides with a plasmoid density that strictly
larger than the above. Therefore, if we choose $\mathcal{E}_{\mathrm{I/II}}$
such that $n_{p}(z_{c})>an_{a}$ for $a\gg 1$, then the collisions the electron
experiences are completely dominated by collisions with the cold Maxwellian
throughout its entire orbit. Above this energy, the electron collides
considerably with the ambient electrons in the extremities of the plasmoid as
well as with the plasmoid electrons in the core, so the distribution function
at these higher energies is not necessarily Maxwellian. Hence, the upper bound
for $\mathcal{E}_{\mathrm{I/II}}$ is expressed as

$\mathcal{E}_{\mathrm{I/II}}<-T\ln a.$ (36)

The lower bound is fixed by collisions with the cut-off Maxwellian in region I
being well-approximated by collisions with the full Maxwellian. The simplest
way to guarantee this is to have

$\frac{f_{0}(\mathcal{E}_{\mathrm{I/II}})}{f_{0}(-e\phi_{m})}<\frac{1}{a}$
(37)

for $a\gg 1$. Then, the lower bound for $\mathcal{E}_{\mathrm{I/II}}$ is given
by

$\mathcal{E}_{\mathrm{I/II}}>-e\phi_{m}+T\ln a.$ (38)

###  2.4 Deriving the kinetic equation for hot trapped electrons (region II)

The kinetic equation in region II is given by

$\frac{\partial f_{\mathrm{II}}}{\partial t}+v_{\parallel}\frac{\partial
f_{\mathrm{II}}}{\partial z}-e\frac{\partial\phi}{\partial t}\frac{\partial
f_{\mathrm{II}}}{\partial\mathcal{E}_{\parallel}}=C(f_{\mathrm{II}},f_{\mathrm{I}})+C(f_{\mathrm{II}},f_{\mathrm{II}}+f_{\mathrm{III}})+\sum_{k}C_{e,ik}(f_{\mathrm{II}}).$
(39)

As the intermediate region, the ordering is most complex for the hot trapped
electrons. More care must also be taken when considering terms with time
derivatives as the collision frequency is lower in region II than region I.
The typical velocity in region I is of order $\sqrt{2e\phi_{m}/m_{e}}\sim
v_{T_{a}}$, so the collision frequency with the plasmoid in region II is of
the same order as in region III. In region II the bounce frequency is of the
same order as in region I (Eq. (25)). Hence we write

$C(f_{\mathrm{II}},f_{\mathrm{I}})\sim\sum_{k}C_{e,ik}(f_{\mathrm{I}})\sim\nu_{p}(v_{T_{a}})f_{\mathrm{II}},$
(40)
$C(f_{\mathrm{II}},f_{\mathrm{II}}+f_{\mathrm{III}})\sim\nu_{h}\mathrm{e}^{-\nu_{h}t}f_{\mathrm{II}},$
(41) $v_{\parallel}\frac{\partial f_{\mathrm{II}}}{\partial
z}\sim\nu_{B}f_{\mathrm{II}},$ (42)

noting that both the collision frequency with the plasmoid and the heating
rate are much smaller than the bounce frequency:

$\frac{\nu_{p}(v_{T_{a}})}{\nu_{B}}=\mu\ll 1,$ (43)
$\frac{\nu_{h}\mathrm{e}^{-\nu_{h}t}}{\nu_{B}}=\frac{4}{3\sqrt{\pi}}\frac{n_{a}L_{p}\mathrm{e}^{-\nu_{h}t}}{N_{p}}\mu\ll
1.$ (44)

As in region I, we assume that the terms containing time derivatives
correspond to a timescale much longer than the bounce period.

Then, the shortest timescale in the system is the bounce period, which leads
to the lowest-order equation

$v_{\parallel}\frac{\partial f_{\mathrm{II}}}{\partial z}=0,$ (45)

meaning that $f_{\mathrm{II}}$ (and hence $f$ as a whole) is independent of
$z$. The higher-order kinetic equation can then be bounce-averaged, yielding

$\frac{\partial f_{\mathrm{II}}}{\partial t}-\frac{1}{\tau}\frac{\partial
J}{\partial t}\frac{\partial
f_{\mathrm{II}}}{\partial\mathcal{E}_{\parallel}}=\left<C(f_{\mathrm{II}},f_{\mathrm{I}})+\sum_{k}C_{e,ik}(f_{\mathrm{II}})\right>+\left<C(f_{\mathrm{II}},f_{\mathrm{II}}+f_{\mathrm{III}})\right>$
(46)

where

$J(\mathcal{E}_{\parallel},t)=\oint
m_{e}v_{\parallel}\,\mathrm{d}z=\sqrt{2m_{e}}\oint\sqrt{\mathcal{E}_{\parallel}+e\phi}\,\mathrm{d}z$
(47)

is the second adiabatic invariant for an electron bouncing in the well.

Now we analyse the timescales on which the time-dependent terms act and
compare them to other timescales. Since $f_{\mathrm{II}}$ is equal to $f_{0}$
at $\mathcal{E}=\mathcal{E}_{\mathrm{I/II}}$ and equal to $f_{a}$ at
$\mathcal{E}_{\parallel}=0$ it serves to perform this analysis at each
boundary. At the trapped-passing separatrix, $\mathcal{E}_{\parallel}=0$,
$f_{\mathrm{II}}$ must be equal to $f_{a}$, which is constant in time; the
first term on the left hand side of Eq. (46) vanishes. The second term
represents the adiabatic change in electron energy as the well expansion,
which in Aleynikov et al., (2019) was shown to occur on the heating timescale:

$-\frac{1}{\tau}\frac{\partial J}{\partial t}\frac{\partial
f_{\mathrm{II}}}{\partial\mathcal{E}_{\parallel}}\sim\nu_{h}\mathrm{e}^{-\nu_{h}t}f_{\mathrm{II}}.$
(48)

At $\mathcal{E}=\mathcal{E}_{\mathrm{I/II}}$ Eq. (48) also holds. However,
since the Maxwellian in region I has a temperature that changes in time, the
time derivative of $f_{\mathrm{II}}$ is not zero here:

$\frac{\partial f_{\mathrm{II}}}{\partial
t}\bigg{|}_{\mathcal{E}=\mathcal{E}_{\mathrm{I/II}}}=\left[\frac{1}{\eta}\frac{\mathrm{d}\eta}{\mathrm{d}t}+\left(\frac{\mathcal{E}_{\mathrm{I/II}}}{T}-\frac{3}{2}\right)\frac{1}{T}\frac{\mathrm{d}T}{\mathrm{d}t}\right]f_{\mathrm{II}}\bigg{|}_{\mathcal{E}=\mathcal{E}_{\mathrm{I/II}}}.$
(49)

We can choose $|\mathcal{E}_{\mathrm{I/II}}|\sim T$, and we see that
$(\partial\eta/\partial t)/\eta\sim(\partial T/\partial t)/T$, so Eq. (49) can
be approximated by

$\frac{\partial f_{\mathrm{II}}}{\partial
t}\bigg{|}_{\mathcal{E}=\mathcal{E}_{\mathrm{I/II}}}\sim\frac{1}{T}\frac{\mathrm{d}T}{\mathrm{d}t}f_{\mathrm{II}}\bigg{|}_{\mathcal{E}=\mathcal{E}_{\mathrm{I/II}}}.$
(50)

That is, the timescale on which the term acts, which we define via the
frequency

$\nu_{t}=\frac{1}{T}\frac{\mathrm{d}T}{\mathrm{d}t}=\frac{T_{a}}{T}\nu_{h}\mathrm{e}^{-\nu_{h}t},$
(51)

is the time taken for $T$ to increase by a factor of $\mathrm{e}$. When $T$ is
small this can be a very short time, so $\nu_{t}$ cannot simply be assumed to
be small compared to the collision time.

The profiles from the modified self-similar solution give

$\frac{\nu_{t}}{\nu_{p}(v_{T_{a}})}\sim\frac{4}{3\sqrt{\pi}}\mathrm{e}^{-\nu_{h}t}(1-\mathrm{e}^{-\nu_{h}t})^{\frac{1}{2}}\left(1+\frac{2}{\pi\sqrt{3}}\sqrt{\frac{m_{i}}{m_{e}}}(\nu_{h}t)^{-\frac{3}{2}}\right)^{-1},$
(52)

which is always much less than unity for the plasma parameters given in the
Introduction; we write

$\frac{\nu_{t}}{\nu_{p}(v_{T_{a}})}\ll 1.$ (53)

So, the collision term with the with the plasmoid corresponds to the shortest
timescale in Eq. (46); the lowest-order kinetic equation is therefore

$\left<C(f_{\mathrm{II}},f_{\mathrm{I}})+\sum_{k}C_{e,ik}(f_{\mathrm{II}})\right>=0.$
(54)

The distribution function must continuous in collisional kinetic problems,
hence Eq. (54) must be solved with boundary conditions ensuring continuity:

$f_{\mathrm{II}}(\mathcal{E}=\mathcal{E}_{\mathrm{I/II}})=f_{0}(\mathcal{E}=\mathcal{E}_{\mathrm{I/II}}),$
(55)
$f_{\mathrm{II}}(\mathcal{E}_{\parallel}=0)=f_{a}(\mathcal{E}_{\parallel}=0).$
(56)

The higher-order equation in region II is

$\frac{\partial f_{\mathrm{II}}}{\partial t}-\frac{1}{\tau}\frac{\partial
J}{\partial t}\frac{\partial
f_{\mathrm{II}}}{\partial\mathcal{E}_{\parallel}}=\left<C(f_{\mathrm{II}},f_{\mathrm{II}}+f_{\mathrm{III}})\right>,$
(57)

which describes the heating and expansion timescales.

###  2.5 The quasi-equilibrium problem

The distribution function in region II is obtained by solving Eq. (54), which
must be be supplemented by boundary conditions (55),(56) enforcing continuity
of the distribution function into regions I and III. Conceptually, the kinetic
problem in region II describes a quasi-equilibrium (QE); hot trapped electrons
experience rapid collisions against a Maxwellian (and are isotropised by
collisions with ions), but the tail of the distribution is forced to meet a
Maxwellian of a different temperature at the trapped-passing separatrix.

When $T\ll T_{a}$ collisions with the distribution in region I are well-
approximated by collisions with the full Maxwellian:
$C(\cdot,f_{\mathrm{I}})\approx C(\cdot,f_{0})$. Since it only remains to
solve the kinetic problem in region II, it is unnecessary to have a subscript,
so we write $f=f_{\mathrm{II}}$ in region II. Hence the QE kinetic equation
can be written as

$\left<C_{\mathrm{QE}}(f)\right>=0$ (58)

for

$C_{\mathrm{QE}}(f)=C(f,f_{0})+\sum_{k}C_{e,ik}(f).$ (59)

Further, owing to the fact that collisions are linearised in terms of
collisions against a full Maxwellian, the lower boundary condition Eq. (55)
can actually be applied at $\mathcal{E}=-e\phi_{m}$ rather than
$\mathcal{E}_{\mathrm{I/II}}$.

###  2.6 Range of validity of the ordering

The ordering developed in this section is based upon the self-similar
expansion in Aleynikov et al., (2019), modified to provide plausible profiles
at the later stages of the expansion, given a line-integrated plasmoid density
$N_{p}=10^{22}\,\mathrm{m}^{-2}$ in an ambient plasma of density electron
$n_{a}=5\times 10^{19}$ at a temperature of $5\,\mathrm{keV}$. The
requirements of the ordering are that during most of the expansion the
potential height is of order the ambient temperature:

$e\phi_{m}\sim T_{a},$ (60)

that the plasmoid is transparent to passing and hot trapped electrons:

$\mu=\frac{\nu_{p}(v_{T_{a}})}{\nu_{T}}\sim\frac{\nu_{p}(v_{T_{a}})}{\nu_{B}}\ll
1,$ (61)

that the heating rate is much lower than the collision frequency with the
plasmoid:

$\frac{\nu_{h}\mathrm{e}^{-\nu_{h}t}}{\nu_{p}(v_{T_{a}})}\ll 1,$ (62)
$\frac{\nu_{h}\mathrm{e}^{-\nu_{h}t}}{\nu_{p}(v_{T})}\ll 1,$ (63)

and that the time taken for the plasmoid electron temperature to increase by a
factor of $\mathrm{e}$ is much larger than the collision time:

$\frac{1}{T}\frac{\mathrm{d}T}{\mathrm{d}t}\ll\nu_{p}(v_{T_{a}}),$ (64)
$\frac{1}{T}\frac{\mathrm{d}T}{\mathrm{d}t}\ll\nu_{p}(v_{T}).$ (65)

The transparency of the plasmoid is dependent upon the line-integrated
plasmoid density not being too large, and $e\phi_{m}$ being of order $T_{a}$
is dependent upon the line-integrated plasmoid density not being too small;
satisfying both conditions does formally somewhat constrain the values of the
line-integrated density.

However, the relative simplicity of the resulting kinetic problem provides
strong motivation for using the formalism: we immediately obtain the
distribution function in two out of three regions, and in the remaining region
we must solve a steady-state kinetic equation. Then, the macroscopic expansion
is described by a kinetic equation of which velocity moments can be taken, in
order to obtain much simpler time-dependent equations than one would by
including the time-dependent terms directly in the kinetic problem. In this
sense, the formalism is analogous to the Braginskii equations, but valid for
systems with a long rather than short hot electron mean free path.

We reiterate that the kinetic problem in region I and II was formally derived
assuming $T\ll T_{a}$, which permitted a choice of
$\mathcal{E}_{\mathrm{I/II}}$ that guaranteed a Maxwellian distribution in
region I, collisions with which are well-approximated by collisions with the
full Maxwellian. This ultimately allowed the linearisation of the kinetic
problem in region II in terms of collisions with the full Mawwellian. However,
we note that when $T=T_{a}$, the entire distribution function will be a single
Maxwellian, yet Eq. (58) would still be satisfied. This is because electrons
in region II would be colliding with Maxwellian electrons with the same
temperature in regions I and III. We therefore conclude that the whole
formalism is valid when $T\ll T_{a}$ and when $T\rightarrow T_{a}$.

Therefore the formalism can actually accurately model the expansion outside of
its formal ordering (which has $T\ll T_{a}$), which is characteristic of
robust simplifications of kinetic problems, such as the Braginskii equations,
which often achieve a level of qualitative correctness even when the mean free
path is long and the distribution function is not very close to a Maxwellian.
By the same argument one could expect that the formalism here is qualitatively
correct somewhat outside of the range of parameters that leads to the
ordering.

As mentioned in the previous subsection, the boundary condition (55) can be
applied at $\mathcal{E}=-e\phi_{m}$ in the QE problem. This solves the problem
of how to choose $\mathcal{E}_{\mathrm{I/II}}$ when $T$ approaches $T_{a}$ and
the well becomes too shallow to contain several $T$: when solving the QE
problem we can always choose $\mathcal{E}_{\mathrm{I/II}}=-e\phi_{m}$.

The formalism has been developed specifically with pellet plasmoids in mind,
and essentially models plasmoid expansion with ‘intermediate’ line-integrated
densities. An alternate approach, which is more suited to the abstract study
of plasmoid expansion, is to consider the limit as the line-integrated density
goes to zero or infinity. Then, the ratio $e\phi_{m}/T_{a}$ is a large or
small parameter on which an ordering may be based.

When the line-integrated density is very large, the plasmoid and ambient
temperatures will equilibrate before the plasmoid density is comparable to the
ambient density. Then, the expansion can be described simply with Maxwellian
electrons from an early stage. If instead it is very small, the densities
become comparable well before the temperatures have equilibrated. In our case,
with an intermediate line-integrated density, the two occur at approximately
the same time; certainly one cannot assume $T=T_{a}$ or $n_{pm}\approx n_{a}$
from an early stage.

When constructing the ordering, the opacity $\mu$ was given assuming that
intra-species collisions dominate; ambient ions collide most quickly with
plasmoid ions and ambient electrons collide most quickly with plasmoid
electrons. This is the case when the thermal velocities of ions and electrons
are disparate. However, if the thermal velocities of an ion population $k$ and
an electron population are comparable, then the friction of the ions on the
electrons is actually $m_{ik}/m_{e}$ times larger than the ion-ion collision
frequency (Helander & Sigmar,, 2002). The thermal velocities of the ambient
ions and plasmoid electrons are actually comparable for a brief window of time
where $T$ is extremely small. However, collisions of the ambient ions with
plasmoid electrons in this regime cause rapid heating of the plasmoid
electrons, therefore driving $T$ up and out of the regime where the plasmoid
electrons and ambient ions have a comparable thermal velocity. Hence these
collisions are negligible outside of the very early stages of the plasmoid
expansion, where the ordering is, anyway, not satisfied; we restrict our
attention to times later than this.

###  2.7 Expressing the quasi-equilibrium equation in the variables
$(\mathcal{E}_{\parallel},\mathcal{E}_{\perp})$

The collision operator against $f_{0}$ in the variables $(v,\theta,z,t)$ for
pitch-angle $\theta$ (assuming that $f$ is independent of the azimuthal angle
of the velocity $\varphi$) is given by

$\begin{split}C(f,f_{0})=\frac{e^{4}\ln\Lambda}{4\pi\varepsilon_{0}^{2}m_{e}^{2}}\Bigg{\\{}&\frac{1}{v^{3}}\frac{g^{\prime}(x)}{v_{T}}\frac{1}{2\sin\theta}\frac{\partial}{\partial\theta}\left(\sin\theta\frac{\partial
f}{\partial\theta}\right)+\\\ &\frac{1}{v^{2}}\frac{\partial}{\partial
v}\left[\frac{x^{3}g^{\prime\prime}(x)}{v_{T}}\left(f+\frac{T}{m_{e}v}\frac{\partial
f}{\partial v}\right)\right]\Bigg{\\}},\end{split}$ (66)

where $x=v/v_{T}$, $g(x)$ is the function

$g(x)=\int\left|\mathbf{v}-\mathbf{v}^{\prime}\right|f_{0}\,\mathrm{d}^{3}v=n_{0}v_{T}\left[\left(x+\frac{1}{2x}\right)\mathrm{erf}(x)+\frac{1}{\sqrt{\pi}}\mathrm{e}^{-x^{2}}\right],$
(67)

and

$n_{0}=\eta\mathrm{e}^{\frac{e\phi}{T}}$ (68)

is the density of core electrons (Helander & Sigmar,, 2002). We note that when
$x$ is large, i.e. we consider collisions of an electron with an energy much
larger than $T$ with the Maxwellian, then both $g^{\prime}(x)$ and
$x^{3}g^{\prime\prime}(x)$ are well-approximated by $n_{0}v_{T}$.

Similarly, assuming that collisions with ions are well-approximated by pitch-
angle scattering, we have

$C_{e,ik}(f)=\frac{e^{4}\ln\Lambda}{4\pi\varepsilon_{0}^{2}m_{e}^{2}}\left[\frac{1}{v^{3}}Z_{k}^{2}n_{ik}\frac{1}{2\sin\theta}\frac{\partial}{\partial\theta}\left(\sin\theta\frac{\partial
f}{\partial\theta}\right)\right]$ (69)

for ion charge $Z_{k}$ and density $n_{ik}$.

The collision operator is given by the divergence of the collisional flux
$\mathbf{F}$ in velocity space:

$C(f,f_{0})=\nabla_{\mathbf{v}}\cdot\mathbf{F},$ (70)

so it can always be expressed in the form

$C(f,f_{0})=|J|\nabla_{\mathbf{w}}\cdot\tilde{\mathbf{F}},$ (71)

where $|J|$ is the Jacobian of the transformation between coordinates
$\mathbf{v}$ and $\mathbf{w}$,

$J=\mathrm{det}\left(\frac{\partial w_{i}}{\partial v_{j}}\right),$ (72)

and $\tilde{\mathbf{F}}$ is the collisional flux in $\mathbf{w}$ phase-space.
Noting that

$\mathrm{det}\left(\frac{\partial(\mathcal{E}_{\parallel},\mathcal{E}_{\perp},\varphi)}{\partial(\mathbf{v})}\right)=m_{e}^{2}v_{\parallel}$
(73)

we seek to transform Eq. (66) into the form

$C(f,f_{0})=Av_{\parallel}\nabla_{(\mathcal{E}_{\parallel},\mathcal{E}_{\perp})}\cdot\tilde{\mathbf{F}}$
(74)

for some constant $A$. We find that

$\begin{split}C(f,f_{0})=&\frac{e^{4}\ln\Lambda}{4\pi\varepsilon_{0}^{2}m_{e}}v_{\parallel}\nabla_{(\mathcal{E}_{\parallel},\mathcal{E}_{\perp})}\cdot\\\
&\left[f_{0}\bigg{(}\mathbf{r}\mathbf{r}\mathcal{E}_{\perp}\frac{v_{\parallel}}{v^{3}}\frac{g^{\prime}(x)}{v_{T}}+T\frac{\mathbf{s}\mathbf{s}}{v^{5}v_{\parallel}}\frac{x^{3}g^{\prime\prime(x)}}{v_{T}}\bigg{)}\nabla_{(\mathcal{E}_{\parallel},\mathcal{E}_{\perp})}\left(\frac{f}{f_{0}}\right)\right],\end{split}$
(75)

where

$\mathbf{r}=(1,-1),\quad\mathbf{s}=\left(v_{\parallel}^{2},v_{\perp}^{2}\right)$
(76)

and $\mathbf{r}\mathbf{r}$, $\mathbf{s}\mathbf{s}$ represent dyadic products.
Similarly,

$C_{e,ik}(f)=\frac{e^{4}\ln\Lambda}{4\pi\varepsilon_{0}^{2}m_{e}}v_{\parallel}\nabla_{(\mathcal{E}_{\parallel},\mathcal{E}_{\perp})}\cdot\left[f_{0}\left(\mathbf{r}\mathbf{r}\mathcal{E}_{\perp}\frac{v_{\parallel}}{v^{3}}Z_{k}^{2}n_{ik}\right)\right].$
(77)

$\mathbf{r}\mathbf{r}$ is the tensor associated with pitch-angle scattering,
altering parallel and perpendicular energies such that their sum is unchanged.
$\mathbf{s}\mathbf{s}$ is associated with energy exchange, which affects
parallel and perpendicular energies equally.

We must bounce-average the collision operators (75) and (77). In order to
obtain a ‘useful’ expression, we must be able to commute divergence in
$(\mathcal{E}_{\parallel},\mathcal{E}_{\perp})$ and the orbit integral. The
following Lemma is useful with regards to commuting the divergence and orbit
integral: for a vector-valued function $\mathbf{F}$, we have

$\nabla_{(\mathcal{E}_{\parallel},\mathcal{E}_{\perp})}\cdot\oint\mathbf{F}\,\mathrm{d}z=\oint\nabla_{(\mathcal{E}_{\parallel},\mathcal{E}_{\perp})}\cdot\mathbf{F}\,\mathrm{d}z+4\nabla_{(\mathcal{E}_{\parallel},\mathcal{E}_{\perp})}z_{c}\cdot\mathbf{F}(z_{c}).$
(78)

In order for the divergence and orbit integral to commute, we must have
$\nabla_{(\mathcal{E}_{\parallel},\mathcal{E}_{\perp})}z_{c}\cdot\mathbf{F}(z_{c})=0$;
either $\mathbf{F}(z_{c})$ or
$\nabla_{(\mathcal{E}_{\parallel},\mathcal{E}_{\perp})}z_{c}$ vanishes, or
$\nabla_{(\mathcal{E}_{\parallel},\mathcal{E}_{\perp})}z_{c}$ is orthogonal to
$\mathbf{F}(z_{c})$.

We observe that the term associated with pitch-angle scattering is
proportional to $v_{\parallel}$, which (by definition) vanishes when
$z=z_{c}$. So, the divergence and orbit integral commute in the pitch-angle
scattering term.

As for the energy exchange term, we observe that

$\nabla_{(\mathcal{E}_{\parallel},\mathcal{E}_{\perp})}z_{c}=\frac{\partial
z_{c}}{\partial\mathcal{E}_{\parallel}}\left(\begin{matrix}1\\\
0\end{matrix}\right),$ (79)

and

$(\mathbf{s}\mathbf{s})(z=z_{c})=\left(\begin{matrix}0&0\\\ 0&v_{\perp}^{4}\\\
\end{matrix}\right),$ (80)

which means that with respect to this term,

$\nabla_{(\mathcal{E}_{\parallel},\mathcal{E}_{\perp})}z_{c}\cdot\mathbf{F}(z_{c})\propto\left(\begin{matrix}1\\\
0\end{matrix}\right)\cdot\left[\left(\begin{matrix}0&0\\\
0&v_{\perp}^{4}\end{matrix}\right)\nabla_{(\mathcal{E}_{\parallel},\mathcal{E}_{\perp})}\left(\frac{f}{f_{0}}\right)\right]=0,$
(81)

so the orbit integral and divergence commute for the energy exchange term.
Since both $f$ and $f_{0}$ are independent of $z$, they and their derivatives
may be brought outside the orbit integral, giving the bounce-averaged QE
collision operator

$\begin{split}\left<C_{\mathrm{QE}}(f)\right>=\frac{e^{4}\ln\Lambda}{4\pi\varepsilon_{0}^{2}m_{e}\tau}\nabla_{(\mathcal{E}_{\parallel},\mathcal{E}_{\perp})}\cdot\bigg{[}&f_{0}\bigg{(}\mathbf{r}\mathbf{r}\mathcal{E}_{\perp}\oint\frac{v_{\parallel}}{v^{3}}\left(\frac{g^{\prime}(x)}{v_{T}}+\sum_{k}Z_{k}^{2}n_{ik}\right)\,\mathrm{d}z\\\
&+T\oint\frac{\mathbf{s}\mathbf{s}}{v^{5}v_{\parallel}}\frac{x^{3}g^{\prime\prime(x)}}{v_{T}}\,\mathrm{d}z\bigg{)}\nabla_{(\mathcal{E}_{\parallel},\mathcal{E}_{\perp})}\left(\frac{f}{f_{0}}\right)\bigg{]}.\end{split}$
(82)

The quasi-equilibrium equation is given by setting the above to zero:

$\begin{split}\nabla_{(\mathcal{E}_{\parallel},\mathcal{E}_{\perp})}\cdot\bigg{[}&f_{0}\bigg{(}\mathbf{r}\mathbf{r}\mathcal{E}_{\perp}\oint\frac{v_{\parallel}}{v^{3}}\left(\frac{g^{\prime}(x)}{v_{T}}+\sum_{k}Z_{k}^{2}n_{ik}\right)\,\mathrm{d}z\\\
&+T\oint\frac{\mathbf{s}\mathbf{s}}{v^{5}v_{\parallel}}\frac{x^{3}g^{\prime\prime(x)}}{v_{T}}\,\mathrm{d}z\bigg{)}\nabla_{(\mathcal{E}_{\parallel},\mathcal{E}_{\perp})}\left(\frac{f}{f_{0}}\right)\bigg{]}=0,\end{split}$
(83)

which is in the form of an anisotropic steady-state diffusion problem in
$(\mathcal{E}_{\parallel},\mathcal{E}_{\perp})$ space:

$\nabla_{(\mathcal{E}_{\parallel},\mathcal{E}_{\perp})}\cdot\left[D_{\mathrm{QE}}\nabla_{(\mathcal{E}_{\parallel},\mathcal{E}_{\perp})}\left(\frac{f}{f_{0}}\right)\right]=0$
(84)

for

$D_{\mathrm{QE}}=D_{\mathrm{QE,S}}+D_{\mathrm{QE,F}},$ (85)

where the diffusion tensor associated with pitch-angle scattering is given by

$D_{\mathrm{QE,S}}=f_{0}\mathcal{E}_{\perp}\oint\frac{v_{\parallel}}{v^{3}}\left(\frac{g^{\prime}(x)}{v_{T}}+\sum_{k}Z_{k}^{2}n_{ik}\right)\,\mathrm{d}z\left(\begin{matrix}1&-1\\\
-1&1\end{matrix}\right),$ (86)

and that associated with energy exchange is given by

$D_{\mathrm{QE,F}}=f_{0}T\oint\frac{1}{v^{5}v_{\parallel}}\left(\begin{matrix}v_{\parallel}^{4}&v_{\parallel}^{2}v_{\perp}^{2}\\\
v_{\parallel}^{2}v_{\perp}^{2}&v_{\perp}^{4}\end{matrix}\right)\frac{x^{3}g^{\prime\prime}(x)}{v_{T}}\,\mathrm{d}z.$
(87)

Together with quasineutrality,

$\int f\,\mathrm{d}^{3}v=n_{e}=\sum_{k}Z_{k}n_{ik},$ (88)

Eq. (84) with the boundary conditions (55),(56) (noting that we write
$f=f_{\mathrm{II}}$) provides a unique solution for $f$ and $\phi$ in terms of
the parameters $\eta$ and $T$. However, these parameters are not know a
priori.

It should be noted that up to this point we have assumed that the potential is
monotonically decreasing, which is the case when the density profile is
monotonically decreasing. If this is not the case, i.e. the potential has more
than one peak, then there are actually multiple trapped electron populations
that must be treated independently. Some electrons can explore the region
encompassed by only one peak, and others, still trapped in the potential as a
whole, can explore more than one. This situation greatly complicates the
kinetic problem and is of secondary importance in this paper as we are
concerned with the expansion of a plasmoid with a potential that is initially
single-peaked; we expect, and observe, as will be shown in Section 3.4, the
profile to remain single-peaked when the high temperature of the ambient
plasma is accounted for. There is one exception to the inapplicability of the
foregoing model to multiply-peaked electric potential wells: when $T=T_{a}$.
In this case the solution to the QE problem is, anyway, the ambient Maxwellian
$f_{a}$, which is correct even when multiple peaks are present. It will be
seen that the cold-fluid model for ions does produce a multiply-peaked
electric potential during later stages of the expansion, but at this point $T$
is of order $T_{a}$, so the solution to the QE problem is close to a
Maxwellian and we expect the resulting expansion to be qualitatively correct.

###  2.8 The no-net-flux condition

Given some $\eta$ and $T$ we may solve the QE problem as specified in the
previous subsection. However, most of the combinations of $\eta$ and $T$ are
not physically meaningful, since they would not actually establish a steady-
state. Quasineutrality requires that there is no net charge, but most
combinations of $\eta$ and $T$ would cause a very large collisional flux of
electrons into or out of the trapped region of phase-space, causing the
plasmoid to ‘charge up’, quickly resulting in the violation of global
quasineutrality. Therefore a closer look at quasineutrality during the
establishment of the QE state is required.

The global quasineutrality condition is given by

$N_{t}+N_{p}=\sum_{k}Z_{k}N_{ik}$ (89)

for $N_{t}$ the line-integrated density of trapped electrons, $N_{p}$ the
line-integrated density of passing electrons, and $N_{ik}$ the line-integrated
density of ion species $k$. Formally, the magnetic field line we consider is
infinite. However, the entire plasmoid structure is localised, with the
possible exception of the plasmoid density approaching zero asymptotically (if
we use, for example, the self-similar ion density profile from Ref. Aleynikov
et al., (2019)). So, rather than the whole field line, we consider the global
quasineutrality condition on some interval $z\in[-L_{S}/2,L_{S}/2]$ for some
$L_{S}$ much larger than the plasmoid; large enough that the plasmoid density
at the endpoints is negligible compared to the ambient density, and the
electric potential at the endpoints is negligible compared to $T$, $T_{a}$, or
$e\phi_{m}$.

In order to maintain global quasineutrality we require

$\frac{\mathrm{d}N_{t}}{\mathrm{d}t}+\frac{\mathrm{d}N_{p}}{\mathrm{d}t}=\sum_{k}Z_{k}\frac{\mathrm{d}N_{ik}}{\mathrm{d}t}.$
(90)

Since the density is constant far the plasmoid, the terms in the above cease
to change as $L_{S}\rightarrow\infty$; we may take $L_{S}$ arbitrarily large
as we never directly evaluate $N_{p}$ or $N_{ik}$.

Using the results from Appendix A, we see that the time derivative of the
line-integrated density of trapped electrons is given by

$\frac{\mathrm{d}N_{t}}{\mathrm{d}t}=\frac{2\pi}{m_{e}^{2}}\frac{\mathrm{d}}{\mathrm{d}t}\int^{\infty}_{0}\int^{J_{m}}_{0}f\,\mathrm{d}J\,\mathrm{d}\mathcal{E}_{\perp},$
(91)

where $J_{m}=J(\mathcal{E}_{\parallel}=0)$ is the maximum value of the second
adiabatic invariant for a trapped electron, and, knowing that $f=f_{a}$ in the
$\mathcal{E}_{\parallel}>0$ region, the time derivative of the line-integrated
density of passing electrons is given by

$\frac{\mathrm{d}N_{p}}{\mathrm{d}t}=\frac{\mathrm{d}}{\mathrm{d}t}\int^{\frac{L_{S}}{2}}_{-\frac{L_{S}}{2}}n_{a}\mathrm{e}^{\frac{e\phi}{T_{a}}}\mathrm{erfc}\left(\sqrt{\frac{e\phi}{T_{a}}}\right)\,\mathrm{d}z.$
(92)

The full kinetic equation (6) can be bounce-averaged and the left hand side
changed to the independent variables $(J,\mathcal{E}_{\perp},t)$ to yield

$\frac{\partial f}{\partial
t}\bigg{|}_{J}=\left<C(f,f_{\mathrm{I}})\right>+\sum_{k}\left<C_{e,ik}(f)\right>+\left<C(f,f_{\mathrm{II}}+f_{\mathrm{III}})\right>,$
(93)

which, along with the approximation $C(f,f_{\mathrm{I}})\approx C(f,f_{0})$,
gives

$\begin{split}\frac{\mathrm{d}N_{t}}{\mathrm{d}t}=&\frac{2\pi}{m_{e}^{2}}\int^{\infty}_{0}\int^{J_{m}}_{0}\left<C_{\mathrm{QE}}(f)\right>\,\mathrm{d}J\,\mathrm{d}\mathcal{E}_{\perp}+\\\
&\frac{2\pi}{m_{e}^{2}}\int^{\infty}_{0}\int^{J_{m}}_{0}\left<C(f,f_{\mathrm{II}}+f_{\mathrm{III}})\right>\,\mathrm{d}J\,\mathrm{d}\mathcal{E}_{\perp}+\\\
&\frac{2\pi}{m_{e}^{2}}\int^{\infty}_{0}\frac{\partial J_{m}}{\partial
t}f_{a}(\mathcal{E}_{\perp})\,\mathcal{E}_{\perp},\end{split}$ (94)

where we have used the fact that the distribution function is continuous:
$f(J=J_{m})=f(\mathcal{E}_{\parallel}=0)=f_{a}(\mathcal{E}_{\parallel}=0)$. We
then see that Eq. (90) becomes

$\begin{split}&\frac{2\pi}{m_{e}^{2}}\int^{\infty}_{0}\int^{J_{m}}_{0}\left<C_{\mathrm{QE}}(f)\right>\,\mathrm{d}J\,\mathrm{d}\mathcal{E}_{\perp}+\\\
&\frac{2\pi}{m_{e}^{2}}\int^{\infty}_{0}\int^{J_{m}}_{0}\left<C(f,f_{\mathrm{II}}+f_{\mathrm{III}})\right>\,\mathrm{d}J\,\mathrm{d}\mathcal{E}_{\perp}+\\\
&\int^{\frac{L_{S}}{2}}_{-\frac{L_{S}}{2}}n_{a}\frac{e}{T_{a}}\frac{\partial\phi}{\partial
t}\mathrm{e}^{\frac{e\phi}{T_{a}}}\mathrm{erfc}\left(\sqrt{\frac{e\phi}{T_{a}}}\right)\,\mathrm{d}z=\sum_{k}Z_{k}\frac{\mathrm{d}N_{ik}}{\mathrm{d}t},\end{split}$
(95)

where the term associated with $\partial J_{m}/\partial t$ has cancelled out
between $\mathrm{d}N_{t}/\mathrm{d}t$ and $\mathrm{d}N_{p}/\mathrm{d}t$.

The first term on the left hand side is associated with fluxes due to
collisions with the plasmoid; we may write

$\frac{2\pi}{m_{e}^{2}}\int^{\infty}_{0}\int^{J_{m}}_{0}\left<C_{\mathrm{QE}}(f)\right>\,\mathrm{d}J\,\mathrm{d}\mathcal{E}_{\perp}\sim\nu_{p}(v_{T_{a}})N_{p}.$
(96)

The second term on the left hand side is due to heating; we may write

$\frac{2\pi}{m_{e}^{2}}\int^{\infty}_{0}\int^{J_{m}}_{0}\left<C(f,f_{\mathrm{II}}+f_{\mathrm{III}})\right>\,\mathrm{d}J\,\mathrm{d}\mathcal{E}_{\perp}\sim\nu_{h}\mathrm{e}^{-\nu_{h}t}N_{p}.$
(97)

The term on the right hand side is due to the plasma at infinity acting as a
source or sink of ions.

The third term on the left hand side is due to the constant replenishment of
the passing distribution by plasma at infinity, leading to it always having
the form $f_{a}$. We can estimate this term by noting that
$\exp(e\phi/T_{a})\mathrm{erfc}(\sqrt{e\phi/T_{a}})\leq 1$ for $\phi\geq 0$
and using the approximation

$\int^{\frac{L_{S}}{2}}_{-\frac{L_{S}}{2}}\frac{e\phi}{T_{a}}\,\mathrm{d}z\sim
L_{p}\frac{e\phi_{m}}{T_{a}}\sim L_{p}.$ (98)

Writing $L_{p}=N_{p}/n_{m}$ and using Eq. (3) then yields

$\int^{\frac{L_{S}}{2}}_{-\frac{L_{S}}{2}}n_{a}\frac{e}{T_{a}}\frac{\partial\phi}{\partial
t}\mathrm{e}^{\frac{e\phi}{T_{a}}}\mathrm{erfc}\left(\sqrt{\frac{e\phi}{T_{a}}}\right)\,\mathrm{d}z\sim\nu_{h}N_{p}\frac{\mu\sqrt{\frac{3m_{i}}{m_{e}}\nu_{h}t}}{\left((\nu_{h}t)^{\frac{3}{2}}+\frac{2}{\sqrt{3}}\sqrt{\frac{m_{i}}{m_{e}}}\mu\right)^{2}}.$
(99)

With the plasma parameters used in the Introduction, the above is at most of
order $\nu_{h}N_{p}$, and decreases as time advances; it is of the same order
as the heating term (the second on the left hand side of Eq. (95)).

The change in $N_{ik}$ depends upon the (yet unchosen) model for the ions. Of
course, the system for plasmoid ions is necessarily conservative (plasmoid
ions are localised), so the line-integrated plasmoid ion density is constant.
If the system also conserves ambient ions, i.e. the ambient plasma cannot act
as a source of ions, then the line-integrated ambient ion densities are
constant. On the other hand, if the plasma can act as a source for the ambient
ions then the terms associated with the change in their line-integrated
densities is at most on the same order as Eq. (99). This is because an ambient
ion density $n_{ik,a}$ is at most of order $n_{a}/Z_{k}$ and its change is
associated purely with the change in the electric potential.

Therefore the most significant term in Eq. (95) is the first on the left hand
side, which is associated with frequency with which an electron collides with
the plasmoid; the other terms are associated with the heating frequency. Hence
Eq. (95) is well-approximated with simply a vanishing flux associated with
collisions with the plasmoid, which can be expressed as

$\frac{e^{4}\ln\Lambda}{2\varepsilon_{0}^{2}m_{e}^{3}}\int^{\infty}_{0}\left[D_{\mathrm{QE}}\nabla_{(\mathcal{E}_{\parallel},\mathcal{E}_{\perp})}\left(\frac{f}{f_{0}}\right)\right]\bigg{|}_{\mathcal{E}_{\parallel}=0}\cdot\left(\begin{matrix}1\\\
0\end{matrix}\right)\,\mathrm{d}\mathcal{E}_{\perp}=0.$ (100)

Hence, to maintain global quasineutrality, there can (approximately) be no net
collisional flux due to collisions with the plasmoid into the trapped region
of phase-space; we call the above the no-net-flux condition. Intuitively this
makes sense; a steady-state due entirely to collisional fluxes cannot exist if
a consequence of the flux is an immediate violation of quasineutrality.

Since we had two free parameters, $\eta$ and $T$, the no-net-flux condition
fixes one in terms of the other. We choose to keep $T$ as the free parameter.
From Eq. (84) it is clear that although the QE problem implies a steady-state,
the collisional fluxes themselves do not vanish. This is the sense in which QE
is a dynamical steady-state rather than the static steady-state characteristic
of a thermal equilibrium.

###  2.9 Numerical solution to the QE problem

The QE problem (84) with boundary conditions (55), (56) (noting that we write
$f=f_{\mathrm{II}}$ in region II), a self-consistent potential given by
quasineutrality (88), and the zero-net-flux condition (100) was solved
numerically. The plasma was assumed to be hydrogenic: there was a single
species of singly-charged ion ($Z=1$) with the proton mass ($m_{i}$ =
$m_{p}$), following the profile

$n_{i}=n_{a}+N_{ic}\frac{1}{L_{p}\sqrt{\pi}}\mathrm{e}^{-\left(\frac{z}{L_{p}}\right)^{2}},$
(101)

where $L_{p}=2.8\,\mathrm{m}$, $N_{ic}=10^{22}\,\mathrm{m}^{-2}$,
$T=615\,\mathrm{eV}$, $n_{a}=5\times 10^{19}\,\mathrm{m}^{-3}$, and
$T_{a}=5\,\mathrm{keV}$. The Gaussian term in the above is consistent with the
profile in Aleynikov et al., (2019) at $t=20\,\mathrm{\mu s}$ given these
parameters.

Figure 3 shows the properties of the electron distribution function. The top
left plot shows the distribution function in velocity space at $z=0$.
$\mathcal{E}=0$ is indicated by the dashed circle and the trapped-passing
separatrix by the vertical dashed lines. The isotropic passing distribution is
clearly visible as concentric circles, as is the very isotropic core
($\mathcal{E}<0$). Significant flattening of the distribution in the high-
energy part region II is observed. The bottom left plot shows the phase-space
trajectories of electrons in $(\mathcal{E}_{\parallel},\mathcal{E}_{\perp})$
space, clearly indicating flow into and out of the trapped region. The dashed
line indicates $\mathcal{E}=0$. The right boundary is the trapped-passing
separatrix. The top right plot shows the effective phase-space flow velocity
$\mathbf{u}^{*}$, defined via

$\nabla_{(\mathcal{E}_{\parallel},\mathcal{E}_{\perp})}\cdot\left(\mathbf{u}^{*}f\right)=-\left<C_{\mathrm{QE}}(f)\right>,$
(102)

i.e.

$\mathbf{u}^{*}\propto-\frac{D_{\mathrm{QE}}\nabla_{(\mathcal{E}_{\parallel},\mathcal{E}_{\perp})}\left(\frac{f}{f_{0}}\right)}{f}.$
(103)

Electron phase-space flow for $\mathcal{E}<0$ is very weak, since this region
is highly isotropised and essentially conforms to a Maxwellian, which exhibits
no collisional phase-space flux.

The bottom right plot shows the flux through the trapped-passing separatrix,
where

$\Gamma_{S}=-\frac{e^{4}\ln\Lambda}{2\varepsilon_{0}^{2}m_{e}^{3}}\left[D_{\mathrm{QE,S}}\nabla_{(\mathcal{E}_{\parallel},\mathcal{E}_{\perp})}\left(\frac{f}{f_{0}}\right)\right]\bigg{|}_{\mathcal{E}_{\parallel}=0}\cdot\left(\begin{matrix}1\\\
0\end{matrix}\right)$ (104)
$\Gamma_{F}=-\frac{e^{4}\ln\Lambda}{2\varepsilon_{0}^{2}m_{e}^{3}}\left[D_{\mathrm{QE,F}}\nabla_{(\mathcal{E}_{\parallel},\mathcal{E}_{\perp})}\left(\frac{f}{f_{0}}\right)\right]\bigg{|}_{\mathcal{E}_{\parallel}=0}\cdot\left(\begin{matrix}1\\\
0\end{matrix}\right),$ (105)

and $\Gamma=\Gamma_{S}+\Gamma_{F}$. Collisions with the cold Maxwellian always
produces an inflow of electrons in region II due to friction (see $\Gamma_{F}$
in Fig. 3). Pitch-angle scattering may only cause a flow along lines of
constant $\mathcal{E}$, and is seen to eject electrons at low perpendicular
energies and causes an inflow at higher energies. The net flux through the
separatrix vanishes due to the no-net-flux condition.

Figure 3: Top left: Numerical distribution function at $z=0$ in SI units. Top
right: effective phase-space flow velocity $\mathbf{u}^{*}$ (Eq. (103)).
Bottom left: phase-space trajectories of electrons (i.e. the streamlines of
$\mathbf{u}^{*}$). Bottom right: collisional flux into the trapped region
(Eqs. (104),(105)). $v_{c}=\sqrt{2\phi_{m}/m_{e}}$ is the parallel escape
velocity at $z=0$.

###  2.10 Analytical solution to the QE problem

The main difficulty in the QE problem is the presence of bounce integrals,
which are difficult to evaluate analytically in a self-consistent potential.
In a square well, however, owing to the constancy of the electric potential
and plasmoid density, the bounce average operator is the identity,
significantly simplifying matters. Since the phase-space domain of the QE
problem depends only upon the potential height $\phi_{m}$, the distribution
function obtained as a solution to the QE problem in a square-well potential
of height $\phi_{m}$ may be used as an approximation of the solution to the QE
problem in a self-consistent potential.

Assuming that there is a single ion species of charge $Z$, that hot trapped
electrons have energies much larger than $v_{T}$, and that the plasmoid
density greatly exceeds the ambient density (hence quasineutrality is
approximately given by $n_{0}=Zn_{i}$), from Eqs. (66),(69) we see that the
bounce-averaged QE collision operator in a square well is given by

$\left<C_{\mathrm{QE}}(f)\right>=\frac{n_{0}e^{4}\ln\Lambda}{4\pi\varepsilon_{0}^{2}m_{e}^{2}}\left[\frac{1+Z}{v^{3}}\frac{1}{2\sin\theta}\frac{\partial}{\partial\theta}\left(\sin\theta\frac{\partial
f}{\partial\theta}\right)+\frac{1}{v^{2}}\frac{\partial}{\partial
v}\left(f+\frac{T}{m_{e}v}\frac{\partial f}{\partial v}\right)\right].$ (106)

We consider $\left<C_{\mathrm{QE}}(f)\right>=0$ separately in the
$\mathcal{E}<0$ (i.e. $v<v_{c}=\sqrt{2e\phi_{m}/m_{e}}$ in a square well) and
$\mathcal{E}>0$ ($v>v_{c}$) regions of phase-space, using the most convenient
coordinates in each case. It is convenient to write $f=f_{0}+f_{1}$ and solve
for $f_{1}$; in both regions we neglect $v$-diffusion for $f_{1}$ (i.e. the
term proportional to $T/(m_{e}v)$ in the above), leaving only $v$-friction.

In the $\mathcal{E}>0$ region we use the variables
$(v,v_{\parallel}=v\cos\theta)$, meaning we must solve $\mathcal{D}(f_{1})=0$
for

$\mathcal{D}(f_{1})=\frac{1+Z}{2}\frac{\partial}{\partial
v_{\parallel}}\left[(v^{2}-v_{\parallel}^{2})\frac{\partial f_{1}}{\partial
v_{\parallel}}\right]+v\frac{\partial f_{1}}{\partial
v}+v_{\parallel}\frac{\partial f_{1}}{\partial v_{\parallel}}.$ (107)

It is convenient to consider the limit where $v\gg v_{\parallel}$, which
represents the correct limit in the majority of $\mathcal{E}>0$,
$\mathcal{E}_{\parallel}<0$ phase-space:

$\mathcal{D}(f_{1})=\frac{1+Z}{2}v^{2}\frac{\partial^{2}f_{1}}{\partial
v_{\parallel}^{2}}+v\frac{\partial f_{1}}{\partial v},$ (108)

which has superposable solutions (that are even in $v_{\parallel}$) provided
by the separation of variables:

$f_{k}=C_{k}\mathrm{e}^{-\frac{1}{2}\lambda_{k}v^{2}}\cosh\left(v_{\parallel}\sqrt{\frac{2}{1+Z}\lambda_{k}}\right)$
(109)

for constants $\left\\{C_{k}\right\\}$ and $\left\\{\lambda_{k}\right\\}$. The
boundary condition (56) is satisfied by the sum of two solutions:

$f_{1}=f_{a}(\mathcal{E})\frac{\cosh\left(\frac{v_{\parallel}}{v_{T_{a}}}\sqrt{\frac{4}{1+Z}}\right)}{\cosh\left(\frac{v_{c}}{v_{T_{a}}}\sqrt{\frac{4}{1+Z}}\right)}-f_{0}(\mathcal{E})\frac{\cosh\left(\frac{v_{\parallel}}{v_{T}}\sqrt{\frac{4}{1+Z}}\right)}{\cosh\left(\frac{v_{c}}{v_{T}}\sqrt{\frac{4}{1+Z}}\right)}.$
(110)

In the $\mathcal{E}<0$ region we use the variables variables
$(v,\xi=\cos\theta)$. So, we must solve $\mathcal{G}(f_{1})=0$ where

$\mathcal{G}(f_{1})=\frac{1+Z}{2}\frac{\partial}{\partial\xi}\left[\left(1-\xi^{2}\right)\frac{\partial
f_{1}}{\partial\xi}\right]+v\frac{\partial f_{1}}{\partial v}.$ (111)

We note that Legendre polynomials are the eigenfunctions of the operator in
$\xi$, so we write $f$ in this basis:

$f_{1}=\sum_{n=0}^{\infty}a_{n}(v)P_{n}(\xi),$ (112)

which gives the following equations for $\\{a_{n}(v)\\}$:

$v\frac{\partial a_{n}}{\partial v}-\frac{1+Z}{2}n(n+1)a_{n}=0$ (113)

with the solutions

$a_{n}=c_{n}v^{\frac{1+Z}{2}n(n+1)}$ (114)

for $\\{c_{n}\\}$ constants. The continuity of the distribution function at
$v=v_{c}$ provides the expressions for $c_{n}$:

$\begin{split}c_{n}=\frac{2n+1}{2}v_{c}^{-\frac{1+Z}{2}n(n+1)}\Bigg{[}\frac{f_{a}(\mathcal{E}=0)}{\mathrm{cosh}\left(\frac{v_{c}}{v_{T_{a}}}\sqrt{\frac{4}{1+Z}}\right)}\int^{1}_{-1}\mathrm{cosh}\left(\frac{v_{c}}{v_{T_{a}}}\sqrt{\frac{4}{1+Z}}\xi\right)P_{n}(\xi)\,\mathrm{d}\xi-\\\
\frac{f_{0}(\mathcal{E}=0)}{\mathrm{cosh}\left(\frac{v_{c}}{v_{T}}\sqrt{\frac{4}{1+Z}}\right)}\int^{1}_{-1}\mathrm{cosh}\left(\frac{v_{c}}{v_{T}}\sqrt{\frac{4}{1+Z}}\xi\right)P_{n}(\xi)\,\mathrm{d}\xi\Bigg{]}.\end{split}$
(115)

To summarise, the analytical solution to the QE problem in a square well is
given by Eq. (110) in the $\mathcal{E}>0$, $\mathcal{E}_{\parallel}<0$ region,
Eq. (112) in the $\mathcal{E}\leq 0$ region, and $f=f_{a}$ in the
$\mathcal{E}_{\parallel}\geq 0$ region. The phase-space domain of the QE
problem is the same given any potential, so the substitution of
$\sqrt{(2/m_{e})(\mathcal{E}+e\phi_{m})}$ for $v$ and
$\sqrt{(2/m_{e})(\mathcal{E}_{\parallel}+e\phi_{m})}$ for $v_{\parallel}$
yields an (approximate) analytical solution to the QE problem valid in a self-
consistent potential. We refer to this analytical solution in figures as
$f_{\mathrm{an}}$.

Figure 4 shows the analytical distribution function for the same parameters as
those used to produce Fig. 3. The top left plot shows the distribution in
velocity space at $z=0$. The top right plot shows the percentage difference
from the numerical solution given in Fig. 3. The bottom left plot shows the
distribution function at $\mathcal{E}_{\perp}=0$. The bottom right plot shows
the distribution function at $\mathcal{E}_{\parallel}=-e\phi_{m}$. The
qualitative behaviour of the distribution is captured well, in particular the
‘flattening’ of the distribution function in the high-energy part of region
II. The observation that in here the contours of the distribution function are
horizontal lines in $(v_{\parallel},v_{\perp})$ can be explained by the fact
that for the analytical QE distribution function

$\frac{\partial}{\partial
v_{\parallel}}f\left(v_{\parallel},v_{\perp}\right)\propto\frac{2}{v_{T_{a}}^{2}}v_{\parallel}\left(\frac{2}{1+Z}-1\right)+O(v_{\parallel}^{3}),$
(116)

which vanishes to lowest order if $Z=1$.

We note that the simplification made by neglecting the $v$-diffusion term for
$f_{1}$ reduces a formerly second-order problem in $v$ to a first-order
problem. Solving the QE problem in the $\mathcal{E}>0$ part of region II with
the continuity boundary condition then fixes the value of $f$ at
$\mathcal{E}=0$. Similarly, in the $\mathcal{E}<0$ part of region II we only
have the opportunity to apply the continuity boundary condition at
$\mathcal{E}=0$, but not at $\mathcal{E}=\mathcal{E}_{\mathrm{I/II}}$. As a
consequence, as $\mathcal{E}\rightarrow\mathcal{E}_{\mathrm{I/II}}$,
$f_{1}\rightarrow c_{0}$, which does violate this continuity boundary
condition. On the other hand, $c_{0}$ is several orders of magnitude smaller
than $f_{0}$ in region $I$, so the discontinuity is small. This is purely an
artefact of the approximations made to obtain the analytical solution; no such
behaviour is seen in the numerical solution.

Figure 4: Top left: analytical electron distribution $f_{\mathrm{an}}$ in SI
units. Top right: percentage difference between the numerical $f$ (see Fig. 3)
and analytical $f_{\mathrm{an}}$. Bottom left: distributions at
$\mathcal{E}_{\perp}=0$. Bottom right: distributions at
$\mathcal{E}_{\parallel}=-e\phi_{m}$. $v_{c}=\sqrt{2\phi_{m}/m_{e}}$ is the
parallel escape velocity at $z=0$.

###  2.11 Expressions for $\eta$ and $\phi$ in terms of $T$

We have found an analytical solution to the QE problem in terms of the lowest-
order distribution function $f_{0}$, which is uniquely determined by the
parameters $\eta$ and $T$. The electric potential $\phi$ is as yet unknown,
but must be such that quasineutrality (88) is satisfied. Additionally, the no-
net-flux condition (100) must be satisfied. Therefore, we may reduce the
number of unknowns in the system from three ($\eta$, $T$, and $\phi$) to one
(which we choose to be $T$) by imposing quasineutrality and no-net-flux.

The analytical solution to QE, however, does not capture the ‘ejection
structure’ near $\mathcal{E}_{\parallel}=\mathcal{E}_{\perp}=0$ characterising
the flux of electrons out of trapped phase-space, so the no-net-flux condition
may not be used directly. Instead, analysis of the flow of electrons in phase-
space allows us to formulate a condition that is approximately equivalent to
Eq. (100). We note that, owing to the flattening of the distribution function
high-energy part of region II, pitch-angle scattering acts to scatter newly-
trapped electrons to lower parallel velocity while leaving their energy
unchanged. Simultaneously the scattered electrons lose energy via friction.
This can be seen in the streamline plot of Fig. 3. The net effect is that a
large fraction of electrons entering the trapped region of phase-space are
eventually drawn into the $\mathcal{E}<0$ region. In order for there to be no
net flux into the trapped region, the same number of electrons entering the
$\mathcal{E}<0$ region must escape from it, eventually being scattered and
ejected from trapped phase-space (forming the ejection structure).

Therefore, an approximation to the condition (100), which states that there is
no net flux into the trapped region of phase-space, is that there is no net
flux into the $\mathcal{E}<0$ region of phase-space. Since the flux into and
out of this region is characterised by integrals on the $\mathcal{E}=0$ line
it is not necessary to have a description of the ejection structure, whereas
directly applying Eq. (100) requires us to integrate along
$\mathcal{E}_{\parallel}=0$, where the ejection structure strongly affects the
flux. Since the analytical solution to QE potentially has a discontinuous
derivative at $\mathcal{E}=0$, the flux into the $\mathcal{E}<0$ region has
contributions from $\mathcal{E}\rightarrow 0^{+}$ and $\mathcal{E}\rightarrow
0^{-}$. We apply this approximation of the no-net-flux condition to the
analytical solution to QE in a square well.

Since the analytical solution to the QE problem is derived from that posed in
a square well, we will calculate the flux in terms of the variables
$(v,\xi,\varphi)$. Pitch-angle scattering may not change the energy of
electrons, so the pitch-angle scattering terms of the collision operator
cannot contribute to the flux across $\mathcal{E}=0$. In a square well, the
number of electrons entering the $\mathcal{E}<0$ region of phase-space is
equal to the number of electrons entering the $v<v_{c}$ region of phase space.
From Eq. (106) we see that the collisional flux in the $\hat{\mathbf{v}}$
direction at $v=v_{c}$ (assuming that $v_{c}\gg v_{T}$) is given by

$F=-\frac{n_{0}e^{4}\ln\Lambda}{4\pi\varepsilon_{0}^{2}m_{e}^{2}}\frac{1}{v_{c}^{2}}\left(f+\frac{T}{m_{e}v}\frac{\partial
f}{\partial v}\right)\Bigg{|}_{\mathcal{E}=0}.$ (117)

The phase-space area element on the $v=v_{c}$ sphere is given by
$v_{c}^{2}\mathrm{d}\xi\,\mathrm{d}\varphi$, so the net flux into the
$v<v_{c}$ region is given by

$G=\int^{2\pi}_{0}\int^{1}_{-1}v_{c}^{2}\left(F\big{|}_{\mathcal{E}=0^{+}}+F\big{|}_{\mathcal{E}=0^{-}}\right)\,\mathrm{d}\xi\,\mathrm{d}\varphi.$
(118)

Setting the above to zero will provide the relation $\eta$, $T$, and $\phi$
required for no net flux into the $\mathcal{E}<0$ region. It will be seen that
the expression obtained from $G=0$ in the limit $Z\rightarrow\infty$ is quite
accurate, even for $Z=1$. In this limit the analytical solution is given by
$f=f_{a}$ for $\mathcal{E}>0$ and $f=c_{0}+f_{0}$ for $\mathcal{E}<0$, where
$c_{0}=f_{a}(\mathcal{E}=0)-f_{0}(\mathcal{E}=0)$. Then, $G=0$ yields the
relation

$\eta=n_{a}\left(\frac{T}{T_{a}}\right)^{\frac{3}{2}}\left(2-\frac{T}{T_{a}}\right).$
(119)

Now, quasineutrality allows us to express both $\eta$ and $\phi$ in terms of
$T$ and the ion profile. In the variables
$(\mathcal{E}_{\parallel},\mathcal{E}_{\perp},z,t)$ the electron distribution
is independent of $z$, so the electron density is a function of $\phi$, time,
and the parameter $T$ ($\eta$ being already expressed in terms of $T$ via the
above). Therefore the quasineutrality condition (88) demands that

$n_{e}(\phi,t;T)=\sum_{k}Z_{k}n_{ik},$ (120)

which gives an implicit expression for $\phi$ at every point, given some $T$
and some ion density profiles $\\{n_{ik}\\}$. We can immediately obtain an
expression for the potential height when the plasmoid density is high; in this
case the electron density at $z=0$ is well-approximated by

$n_{e}(z=0)=\eta\mathrm{e}^{\frac{e\phi_{m}}{T}},$ (121)

so

$e\phi_{m}=T\left[\ln\left(\frac{n_{e}(z=0)}{n_{a}}\right)+\frac{3}{2}\ln\left(\frac{T_{a}}{T}\right)-\ln\left(2-\frac{T}{T_{a}}\right)\right].$
(122)

This is a modification of the Boltzmann relation for a plasma in thermal
equilibrium, the extra contributions accounting for the fact that in the
quasi-equilibrium state the distribution function is not necessarily
Maxwellian. As expected, the above reduces to the Boltzmann relation for
$T=T_{a}$ (when the distribution becomes Maxwellian), agreeing with the rough
estimate (4) which was used to justify the ordering leading to QE. As noted
earlier, $e\phi_{m}$ given by (122) exceeds that given by the Boltzmann
relation when $T<T_{a}$.

The height of the potential in the numerical solution to the QE problem (Fig.
3), which is calculated entirely self-consistently, and leads to no net flux
into the trapped region, is within 3% of the estimate (122). This is a
remarkable agreement considering that the estimate was derived in the limit
$Z\rightarrow\infty$, but the numerical solution is with $Z=1$. The estimate
also correctly predicts the height of the self-consistent potential for the
solution to the time-dependent kinetic problem with isotropic electrons given
in Arnold et al., (2023), which is strong evidence that the QE state was
established there too.

##  3 Plasmoid expansion

###  3.1 Treatment of ions

Two different models for the ions are considered: a collisionless (Vlasov)
system with hot ambient ions but cold plasmoid ions, and a cold-fluid
expansion. We restrict our attention to a single ion species of charge $Z$.
The collisionless and fluid models are in opposite regimes of collisionality,
which will provide the broadest range of qualitative results for the plasmoid
expansion. The collisionless system is, however, the most physically accurate
model for the plasmoid expansion since ambient ions have a long mean-free-path
relative to the plasmoid size. In fact, the ratio of plasmoid size to the
ambient ion mean free path is the same as ambient electrons, expressed by the
opacity (11). The collisionless system will still model cold plasmoid ions
accurately since they will be initialised with zero velocity spread.

For the collisionless ion expansion we solve the Vlasov equation for ions:

$\frac{\partial f_{i}}{\partial t}+v_{\parallel}\frac{\partial f}{\partial
z}-\frac{Ze}{m_{i}}\frac{\partial\phi}{\partial z}\frac{\partial
f_{i}}{\partial v_{\parallel}}=0.$ (123)

For the cold-fluid expansion we solve the equations

$\frac{\partial n_{i}}{\partial t}+\frac{\partial}{\partial z}(n_{i}u_{i})=0,$
(124) $\frac{\partial u_{i}}{\partial t}+u_{i}\frac{\partial u_{i}}{\partial
z}+\frac{Ze}{m_{i}}\frac{\partial\phi}{\partial z}=0.$ (125)

###  3.2 Treatment of electrons on expansion timescale

Equations (33) and (57) describe the electron dynamics on the expansion and
heating timescales. We note that because they take exactly the same form, Eq.
(57) and the bounce average of Eq. (33) may be combined to yield

$\frac{\partial f}{\partial t}-\frac{1}{\tau}\frac{\partial J}{\partial
t}\frac{\partial
f}{\partial\mathcal{E}_{\parallel}}=\left<C(f,f_{\mathrm{II}}+f_{\mathrm{III}})\right>$
(126)

to describe the both regions I and II (writing the electron distribution in
both regions as $f$). This equation describes the long-term expansion of the
plasmoid, which is driven by the heating term on the right hand side.

By solving the QE problem with the restrictions of quasineutrality and the no-
net-flux condition, the distribution is known up to the parameter $T$; the
evolution of the electron distribution function is completely characterised by
how $T$ changes in time. In analogy with the Braginskii equations we obtain an
evolution equation for $T$ by taking moments of the kinetic equation over
phase-space. When we take moments the quasi-equilibrium distribution serves
the same role as the near-Maxwellian distribution in the Braginskii equations.

In contrast to a local equilibrium distribution, where the parameters $\eta$
and $T$ defining the Maxwellian distribution are dependent upon $z$, in our
case (analogous to a global thermal equilibrium) the parameters of are
independent of $z$, so we also integrate the kinetic equation over $z$ as well
as the ‘momentum-like’ variables (we note that we have already integrated the
kinetic equation over $z$ when it is bounce-averaged). We also need only take
moments over the trapped region of phase-space since the passing electron
distribution is known.

As we seek an equation for $T$ it serves to take the $\mathcal{E}$-moment over
trapped phase-space. The line-integrated energy density of trapped electrons
is given by

$W_{t}=\int_{V_{t}}\mathcal{E}f\,\mathrm{d}^{3}v\,\mathrm{d}z,$ (127)

where $V_{t}$ is the trapped region of phase-space. Taking the
$\mathcal{E}$-moment of Eq. (126) yields

$\frac{\mathrm{d}W_{t}}{\mathrm{d}t}=\frac{\mathrm{d}W_{t}}{\mathrm{d}t}\bigg{|}_{\mathrm{adiabatic}}+\frac{\mathrm{d}W_{t}}{\mathrm{d}t}\bigg{|}_{\mathrm{separatrix}}+\frac{\mathrm{d}W_{t}}{\mathrm{d}t}\bigg{|}_{\mathrm{heating}}$
(128)

for

$\displaystyle\frac{\mathrm{d}W_{t}}{\mathrm{d}t}\bigg{|}_{\mathrm{adiabatic}}$
$\displaystyle=-\frac{2\pi}{m_{e}^{2}}\int^{\infty}_{0}\int^{0}_{-e\phi_{m}}\frac{\partial
J}{\partial
t}f\,\mathrm{d}\mathcal{E}_{\parallel}\,\mathrm{d}\mathcal{E}_{\perp},$ (129)
$\displaystyle\frac{\mathrm{d}W_{t}}{\mathrm{d}t}\bigg{|}_{\mathrm{separatrix}}$
$\displaystyle=\frac{2\pi}{m_{e}^{2}}\frac{\partial J_{m}}{\partial
t}\int^{\infty}_{0}\mathcal{E}_{\perp}f_{a}\big{|}_{\mathcal{E}_{\parallel}=0}\,\mathrm{d}\mathcal{E}_{\perp},$
(130)
$\displaystyle\frac{\mathrm{d}W_{t}}{\mathrm{d}t}\bigg{|}_{\mathrm{heating}}$
$\displaystyle=\int_{V_{t}}\mathcal{E}C(f,f_{\mathrm{II}}+f_{\mathrm{III}})\,\mathrm{d}^{3}v\,\mathrm{d}z.$
(131)

The details of the procedure are contained in Appendix A.

The terms in Eq. (128) are descriptive: ‘adiabatic’ corresponds to the
adiabatic change in the electron energy as the well changes shape,
‘separatrix’ corresponds to electrons crossing the trapped-passing separatrix
due to the changing depth of the potential well $e\phi_{m}$, and ‘heating’
corresponds to collisions of $f$ with the hot electrons.

Although $W_{t}$ is in principle expressible in terms of $T$, any tiny change
in $W_{t}$ (and hence $T$) results in a huge deviation from quasineutrality
due to the exponential dependency of $n_{0}$ on $e\phi/T$. Solving the
evolution equation for $W_{t}$ numerically, inverting the relation to find
$T$, and maintaining quasineutrality requires an impractically small timestep.
Instead, we derive an energy conservation law for electrons and ions which can
be used as an evolution equation for $T$.

The energy conservation law will contain contributions from both passing and
trapped electrons, so it is more convenient to consider the energy contained
within some interval $z\in[-L_{S}/2,L_{S}/2]$ for $L_{S}$ much larger than the
plasmoid. The ambient plasma then acts as an infinite source and sink of
electrons and energy for this section of the field line. The passing
distribution function on this interval being $f_{a}$ can be understood as the
ambient plasma instantly replenishing the passing distribution if it is
altered in any way by interaction with the plasmoid; the ambient plasma
essentially acts as a ‘thermostat’ for the passing distribution. The interval
$L_{S}$ must be large enough that the line-integrated trapped electron energy
density $W_{t}$ changes negligibly as $L_{S}$ increases (i.e. the trapped
electron density is negligible at $|z|=L_{S}/2$).

It is more convenient to work with the line-integrated kinetic energy

$K_{t}=W_{t}+\int^{\frac{L_{S}}{2}}_{-\frac{L_{S}}{2}}e\phi
n_{e,t}\,\mathrm{d}z$ (132)

(where $n_{e,t}$ is the trapped electron density) rather than $W_{t}$, since
we expect the sum of kinetic energies to exhibit a conservation law. Taking
the time derivative of $K_{t}$ yields

$\begin{split}\frac{\mathrm{d}K_{t}}{\mathrm{d}t}&=\frac{\mathrm{d}W_{t}}{\mathrm{d}t}+\int^{\frac{L_{S}}{2}}_{-\frac{L_{S}}{2}}e\frac{\partial\phi}{\partial
t}n_{e,t}\,\mathrm{d}z+\int^{\frac{L_{S}}{2}}_{-\frac{L_{S}}{2}}e\phi\frac{\partial
n_{e,t}}{\partial t}\,\mathrm{d}z\\\
&=\frac{\mathrm{d}W_{t}}{\mathrm{d}t}-\frac{\mathrm{d}W_{t}}{\mathrm{d}t}\bigg{|}_{\mathrm{adiabatic}}+\int^{\frac{L_{S}}{2}}_{-\frac{L_{S}}{2}}e\phi\left(Z\frac{\partial
n_{i}}{\partial t}-\frac{\partial n_{e,p}}{\partial
t}\right)\,\mathrm{d}z,\end{split}$ (133)

where we have used the quasineutrality condition (assuming that there is a
single species of ions with charge $Z$) and the fact that
$n_{e,t}+n_{e,p}=n_{e}$ for $n_{e,p}$ the passing electron density. The latter
two terms correspond to the changing kinetic energies of the ions and passing
electrons. The term involving ion density is given by

$\int^{\frac{L_{S}}{2}}_{-\frac{L_{S}}{2}}Ze\phi\frac{\partial n_{i}}{\partial
t}\,\mathrm{d}z=-\frac{\partial K_{i}}{\partial t}$ (134)

for $K_{i}$ the line-integrated kinetic energy of the ions. This can be
derived by taking velocity moments Eq. (123) or by constructing an energy
equation from Eqs. (124),(125). The same procedure cannot be carried out for
the passing electrons since these are restricted to a certain region of phase
space.

Extending the notion of an orbit integral to electrons with positive parallel
energy is straightforward on a finite $z$ interval:

$\oint
g(\mathcal{E}_{\parallel}>0)\,\mathrm{d}z=2\int^{\frac{L_{S}}{2}}_{-\frac{L_{S}}{2}}g(\mathcal{E}_{\parallel}>0)\,\mathrm{d}z,$
(135)

allowing us to define the second adiabatic invariant $J$ for positive parallel
energies. The line-integrated passing electron energy density on the interval
$z\in[-L_{S}/2,L_{S}/2]$ is given by

$W_{p}=\frac{2\pi}{m_{e}^{2}}\int^{\infty}_{0}\int^{\infty}_{J_{m}}\mathcal{E}f_{a}\,\mathrm{d}J\,\mathrm{d}\mathcal{E}_{\perp},$
(136)

hence

$\frac{\mathrm{d}W_{p}}{\mathrm{d}t}=\frac{2\pi}{m_{e}^{2}}\int^{\infty}_{0}\int^{\infty}_{0}\left(\frac{\mathcal{E}}{T_{a}}-1\right)\frac{\partial
J}{\partial
t}f_{a}\,\mathrm{d}\mathcal{E}_{\parallel}\,\mathrm{d}\mathcal{E}_{\perp}-\frac{2\pi}{m_{e}^{2}}\int^{\infty}_{0}\mathcal{E}_{\perp}\frac{\mathrm{d}J_{m}}{\mathrm{d}t}f_{a}(\mathcal{E}_{\perp})\,\mathrm{d}\mathcal{E}_{\perp}.$
(137)

We note that

$\frac{2\pi}{m_{e}^{2}}\int^{\infty}_{0}\int^{\infty}_{0}\frac{\partial
J}{\partial
t}f_{a}\,\mathrm{d}\mathcal{E}_{\parallel}\,\mathrm{d}\mathcal{E}_{\perp}=-\int^{\frac{L_{S}}{2}}_{-\frac{L_{S}}{2}}e\frac{\partial\phi}{\partial
t}n_{e,p}\,\mathrm{d}z$ (138)

and

$\frac{2\pi}{m_{e}^{2}}\int^{\infty}_{0}\mathcal{E}_{\perp}\frac{\mathrm{d}J_{m}}{\mathrm{d}t}f_{a}(\mathcal{E}_{\perp})\,\mathrm{d}\mathcal{E}_{\perp}=\frac{\mathrm{d}W_{t}}{\mathrm{d}t}\bigg{|}_{\mathrm{separatrix}},$
(139)

so, given that

$K_{p}=W_{p}+\int^{\frac{L_{S}}{2}}_{-\frac{L_{S}}{2}}e\phi
n_{e,p}\,\mathrm{d}z,$ (140)

we obtain the energy conservation law on the interval $z\in[-L_{S}/2,L_{S}/2]$

$\frac{\mathrm{d}}{\mathrm{d}t}\left(K_{t}+K_{p}+K_{i}\right)=\frac{\mathrm{d}W_{t}}{\mathrm{d}t}\bigg{|}_{\mathrm{heating}}+\frac{2\pi}{m_{e}^{2}}\int^{\infty}_{0}\int^{\infty}_{0}\frac{\mathcal{E}}{T_{a}}\frac{\partial
J}{\partial
t}f_{a}\,\mathrm{d}\mathcal{E}_{\parallel}\,\mathrm{d}\mathcal{E}_{\perp}.$
(141)

Evidently, the inclusion of the ion and passing electron energies in the
energy conservation law accounts for the absence of
$\mathrm{d}W_{t}/\mathrm{d}t\big{|}_{\mathrm{adiabatic}}$, since this
represents energy that is extracted from trapped electrons during the
expansion and given to other species. The absence of the separatrix term is
simply due to this contribution cancelling out between $K_{t}$ and $K_{p}$.

The second term on the right hand side of the above arises from the fact that
passing electrons have their energy altered when passing through the time-
varying potential well, and this energy gain (or loss) is absorbed by (or
suffered by) the ambient plasma, which continually supplies passing electrons
following the distribution $f_{a}$. The heating due to this effect is
essentially negligible when $T\ll T_{a}$, but provides a considerable fraction
of the heating power when $T\sim T_{a}$.

The first term on the right hand side represents the collisional heating of
trapped electrons. In Arnold et al., (2023), which treated a high-$Z$
plasmoid, it was shown that the heating rate for cold electrons in a plasmoid
was 3/4 that expected for cold electrons in a homogeneous plasma, since the
acceleration of passing electrons through the potential well decreases their
density and collisionality. That is, given that the per-electron heating rate
of a cold Maxwellian in a homogeneous plasma is $3\nu_{h}(T_{a}-T)$, the per-
electron heating rate for electrons trapped in the potential well was found to
be $(9/4)\nu_{h}(T_{a}-T)$. In our case the distribution function is also
highly isotropic for $\mathcal{E}<0$, so we approximate the collisional
heating term by

$\frac{\mathrm{d}W_{t}}{\mathrm{d}t}\bigg{|}_{\mathrm{heating}}=\frac{9}{4}\nu_{h}N_{\mathcal{E}<0}(T_{a}-T),$
(142)

where $N_{\mathcal{E}<0}$ is the line-integrated density of trapped electrons
with $\mathcal{E}<0$. As shown in Section 2.9, the distribution function is
somewhat less than $f_{a}$ in the $\mathcal{E}>0$ region for $Z=1$, so the
above is an upper bound for the heating rate.

Consolidating trapped and passing electron energies into $K_{e}=K_{t}+K_{p}$
we have the energy conservation law

$\frac{\mathrm{d}}{\mathrm{d}t}\left(K_{e}+K_{i}\right)=Q$ (143)

for heating power

$Q=\frac{9}{4}\nu_{h}N_{\mathcal{E}<0}(T_{a}-T)+\frac{2\pi}{m_{e}^{2}}\int^{\infty}_{0}\int^{\infty}_{0}\frac{\mathcal{E}}{T_{a}}\frac{\partial
J}{\partial
t}f_{a}\,\mathrm{d}\mathcal{E}_{\parallel}\,\mathrm{d}\mathcal{E}_{\perp}.$
(144)

###  3.3 Comparison of the system with earlier work

At this point a clear comparison can be drawn between this investigation and
Aleynikov et al., (2019); Arnold et al., (2021). In Aleynikov et al., (2019);
Arnold et al., (2021) a cold-fluid system for ions was coupled with an energy
conservation law for the plasmoid ions and electrons. The dynamics of the
passing electron distribution were neglected save for the assertion that the
presence of the ambient plasma resulted in a per-electron heating term
$3\nu_{h}T_{a}$ for plasmoid electrons (only the linear heating stage was
treated). The potential was given by the Boltzmann relation and was unbounded
as $|z|\rightarrow\infty$ due to the neglect of passing electrons.

Here, we use an energy conservation law with modified heating terms and
electron kinetic energy derived from a solution to the electron kinetic
problem. The electric potential is given by the quasineutrality condition
according to the electron density of this distribution, and vanishes as
$|z|\rightarrow\infty$. Taking the limit $n_{a}\rightarrow 0$ in Eq. (88) and
the energy conservation law (143) recovers the same system as Aleynikov et
al., (2019) with the exception of the 3/4 factor on the collisional heating
term.

###  3.4 Numerical solutions of the plasmoid expansion system

Numerical solutions to the system created by coupling the energy conservation
law (143), quasineutrality (120), and one of systems describing the ions (the
Vlasov equation (123) or the cold-ion system (124),(125)) were obtained using
the plasma parameters $n_{a}=5\times 10^{19}\,\mathrm{m}^{-3}$,
$T_{a}=5\,\mathrm{keV}$ and $N_{ic}=10^{22}\,\mathrm{m}^{-2}$, the same as in
Section 2.9. The plasma was once more assumed to be hydrogenic. The heating
timescale with these parameters is $\nu_{h}^{-1}=162\,\mathrm{\mu s}$ and the
expansions were run to $300\,\mathrm{\mu s}$, at which point the plasmoid and
ambient temperatures have nearly equalised and the plasmoid density has
dropped to nearly the ambient.

In the collisionless ion expansion the ambient ion distribution was
initialised to a Maxwellian of density $n_{a}$ and temperature $T_{a}$. The
plasmoid ions were initialised at a temperature of $50\,\mathrm{eV}$. In both
the cold-fluid and collisionless expansion the plasmoid was initialised in the
lowest-$z$ grid cell and the ambient uniformly across $z$

The analytical solution to the QE problem given in Section 2.10 was used to
calculate densities and energy densities, respectively appearing in the
quasineutrality condition (120) and the energy conservation law (143). We
neglect any deviation of $f$ from $f_{0}$ in the $\mathcal{E}<0$ region when
evaluating these densities as the difference is negligible. However, we do
take into account the fact that $f$ is flattened in the $\mathcal{E}>0$,
$\mathcal{E}_{\parallel}<0$ region, as this will have a relatively large
impact on the density and energy density. We use Eq. (119) as the expression
for $\eta$ in terms of $T$ when evaluating $f_{0}$. Additionally, we expand
the $\mathrm{cosh}$ functions found in the analytical expression of $f$ to
second order in order to obtain analytical expressions for the moments.

The electron density is then given by

$n_{e}=n_{\mathcal{E}<0}+n^{\mathcal{E}>0,c}_{\mathcal{E}_{\parallel}<0}+n^{\mathcal{E}>0,a}_{\mathcal{E}_{\parallel}<0}+n_{\mathcal{E}_{\parallel}>0}$
(145)

for

$n_{\mathcal{E}<0}=\eta\left(\mathrm{e}^{\frac{e\phi}{T}}\mathrm{erf}\left(\sqrt{\frac{e\phi}{T}}\right)-\frac{2}{\sqrt{\pi}}\sqrt{\frac{e\phi}{T}}\right),$
(146)
$n^{\mathcal{E}>0,c}_{\mathcal{E}_{\parallel}<0}=\frac{2}{\sqrt{\pi}}\eta\left[\sqrt{\frac{e\phi}{T}}-\frac{1}{1+\frac{2}{1+Z}\frac{e\phi_{m}}{T}}\left(\sqrt{\frac{e\phi}{T}}+\frac{1}{3}\frac{2}{1+Z}\left(\frac{e\phi}{T}\right)^{\frac{3}{2}}\right)\right],$
(147)
$n^{\mathcal{E}>0,a}_{\mathcal{E}_{\parallel}<0}=\frac{2}{\sqrt{\pi}}n_{a}\frac{1}{1+\frac{2}{1+Z}\frac{e\phi_{m}}{T_{a}}}\left(\sqrt{\frac{e\phi}{T_{a}}}+\frac{1}{3}\frac{2}{1+Z}\left(\frac{e\phi}{T_{a}}\right)^{\frac{3}{2}}\right),$
(148)

and

$n_{\mathcal{E}_{\parallel}>0}=n_{a}\mathrm{e}^{\frac{e\phi}{T_{a}}}\mathrm{erfc}\left(\sqrt{\frac{e\phi}{T_{a}}}\right).$
(149)

The line-integrated kinetic energies are given by

$K_{e}=K_{\mathcal{E}<0}+K^{\mathcal{E}>0,c}_{\mathcal{E}_{\parallel}<0}+K^{\mathcal{E}>0,a}_{\mathcal{E}_{\parallel}<0}+K_{\mathcal{E}_{\parallel}>0},$
(150)

where

$K_{\mathcal{E}<0}=\frac{3}{2}N_{\mathcal{E}<0}T-\frac{2}{\sqrt{\pi}}\eta\int^{\frac{L_{S}}{2}}_{-\frac{L_{S}}{2}}\left(\frac{e\phi}{T}\right)^{\frac{3}{2}}\,\mathrm{d}z,$
(151)
$K^{\mathcal{E}>0,c}_{\mathcal{E}_{\parallel}<0}=N^{\mathcal{E}>0,c}_{\mathcal{E}_{\parallel}<0}T+\int^{\frac{L_{S}}{2}}_{-\frac{L_{S}}{2}}e\phi
n^{\mathcal{E}>0,c}_{\mathcal{E}_{\parallel}<0}\,\mathrm{d}z,$ (152)
$K^{\mathcal{E}>0,a}_{\mathcal{E}_{\parallel}<0}=N^{\mathcal{E}>0,a}_{\mathcal{E}_{\parallel}<0}T_{a}+\int^{\frac{L_{S}}{2}}_{-\frac{L_{S}}{2}}e\phi
n^{\mathcal{E}>0,a}_{\mathcal{E}_{\parallel}<0}\,\mathrm{d}z,$ (153)

and

$K_{\mathcal{E}_{\parallel}>0}=\frac{3}{2}N_{\mathcal{E}_{\parallel}>0}T_{a}+\frac{1}{\sqrt{\pi}}n_{a}T_{a}\int^{\frac{L_{S}}{2}}_{-\frac{L_{S}}{2}}\sqrt{\frac{e\phi}{T_{a}}}\,\mathrm{d}z,$
(154)

where

$\displaystyle N_{\mathcal{E}<0}=$
$\displaystyle\int^{\frac{L_{S}}{2}}_{-\frac{L_{S}}{2}}n_{\mathcal{E}<0}\,\mathrm{d}z,$
(155) $\displaystyle N^{\mathcal{E}>0,c}_{\mathcal{E}_{\parallel}<0}=$
$\displaystyle\int^{\frac{L_{S}}{2}}_{-\frac{L_{S}}{2}}n^{\mathcal{E}>0,c}_{\mathcal{E}_{\parallel}<0}\,\mathrm{d}z,$
(156) $\displaystyle N^{\mathcal{E}>0,a}_{\mathcal{E}_{\parallel}<0}=$
$\displaystyle\int^{\frac{L_{S}}{2}}_{-\frac{L_{S}}{2}}n^{\mathcal{E}>0,a}_{\mathcal{E}_{\parallel}<0}\,\mathrm{d}z,$
(157) $\displaystyle N_{\mathcal{E}_{\parallel}>0}=$
$\displaystyle\int^{\frac{L_{S}}{2}}_{-\frac{L_{S}}{2}}n_{\mathcal{E}_{\parallel}>0}\,\mathrm{d}z.$
(158)

Figure 5 shows the ion distribution of the collisionless ion expansion at
various times. We see that the self-similar flow velocity matches in the
regions of highest ion density even at late times. Figures 6 and 7 show
quantities derived from the ion distribution or from the flow velocities and
densities of the cold-fluid expansion. In the fluid expansion an ion front
rapidly forms, resulting in very large density gradients. The front appears to
develop an oscillatory character at later times. Sound waves and perhaps even
solitons would develop near the ion front and propagate into the ambient
plasma, but the $z$ grid used was too course to resolve such waves. In
contrast the collisionless expansion exhibits no steep front owing to the fact
that the hot ambient plasma does not ‘pile-up’ on the moving plasmoid; the
majority of passing electrons have large enough parallel velocity to either
pass over the plasmoid or be reflected from the front much more quickly than
the expansion. In both cases the electron temperature approaches $T_{a}$
somewhat more quickly than the estimate (2), a consequence of the term in Eq.
(144) corresponding to the energy lost adiabatically by the passing electrons
but gained by trapped electrons and ions. This term constitutes the majority
of the heating when $T\sim T_{a}$.

Of particular note, and ultimately what is sought in this investigation, is
the energy balance between electrons and ions, plotted in the top left of
Figs. 6 and 7. The two lines represent the fraction of the heating energy
($\int^{t}_{0}Q\,\mathrm{d}t$ for $Q$ in Eq. (143)) deposited into each
species, and in all cases the energy balance tends to a near equal split
between electrons and ions, which is quite remarkable considering the opposite
collisionality regimes for ions in the fluid and Vlasov models. This energy
balance is similar to those predicted by the self-similar expansions in
Aleynikov et al., (2019) and Arnold et al., (2021).

Although not shown here, numerical simulations of the expansion with larger
line-integrated plasmoid densities exhibit an energy balance even closer to an
equal split. A larger line-integrated plasmoid density also results in the
plasmoid and ambient temperature equilibrating before the plasmoid density
drops to the ambient, which furnishes the possibility of modelling the
expansion (and wave generation) from that point with the much simplified
equations resulting from a Maxwellian electron distribution function.

Figure 5: Ion distribution function of the collisionless ion expansion at
various times. The dashed line is the self-similar flow velocity given in
Aleynikov et al., (2019). Figure 6: Derived quantities of the collisionless
ion expansion. Top left: the relative amounts of energy deposited into the
electrons and ions. Top right: the plasmoid electron temperature $T$ and the
estimated temperature evolution given a homogeneous plasma. Bottom: plots of
the electric potential and electron density at various times. Figure 7:
Derived quantities of the cold-fluid ion expansion. Top left: the relative
amounts of energy deposited into the electrons and ions. Top right: the
plasmoid electron temperature $T$ and the estimated temperature evolution
given a homogeneous plasma. Bottom: plots of the electric potential and
electron density at various times.

##  4 Discussion and Conclusions

A model for the expansion of a plasmoid produced by fuel pellet ablation has
been developed that, considering the high complexity of the problem, provides
a relatively simple steady-state electron kinetic problem that was rigorously
derived from the ordering of timescales. The long-term evolution of the
electrons is decsribed by an energy balance obtained by taking moments of the
electron kinetic equation. Ions were described with cold-fluid equations or
with the Vlasov equation; in the latter case the large temperature difference
between ambient and plasmoid ions can be accounted for.

[Discussion of shock; happens when $T_{a}$ is small]

The model presented in this paper is contrasted with earlier work that simply
assumed a Maxwellian electron distribution for plasmoid electrons and entirely
neglected ambient electrons except in a simple heating term. Earlier work also
only considered a cold-fluid model for the ions. The largest pitfall of
earlier work was the combination of an unbounded (and therefore unphysical)
electric potential and the lack of a proper treatment of passing electrons,
which called into question the conclusions about the qualitative character of
the expansion, most notably the electron to ion energy transfer. The
rectification of these issues was the main motivation for developing a model
in a more rigorous manner.

It has been shown that during the expansion electrons reach a ‘quasi-
equilibrium’ state: a dynamical steady-state on the fastest collisional
timescale which establishes an electron distribution that has no net flux
through the trapped-passing separatrix. An analytical solution to the QE
electron kinetic problem was obtained and compared to a numerical solution. An
estimate of the height of the self-consistent electric potential that supports
quasi-equilibrium has been derived. The estimate is consistent with the
Boltzmann relation when the temperatures of the plasmoid electrons and ambient
electrons are equal, and is in agreement with the height of the self-
consistent potential found in the solution to the time-dependent electron
kinetic problem in a high-$Z$ plasmoid (Arnold et al.,, 2023), providing
strong evidence for the establishment of the QE state. The QE kinetic equation
and energy balance can therefore be incorporated into established codes to
describe electron dynamics in a pellet plasmoid.

The quasineutrality condition, no-net-flux condition, and QE kinetic problem
allow a description of the QE distribution function and electric potential in
terms of the plasmoid electron temperature $T$ and the ion densities.
Analogous to the Braginskii equations, the evolution equation for $T$ was
obtained by taking the energy moment over the electron kinetic equation that
holds on the expansion timescale; this evolution equation takes the form of an
energy conservation law for both electrons and ions. The evolution of $T$ is
driven by the energy exchange between passing electrons, trapped electrons,
and ions; heating power initially deposited in the plasmoid electrons by
collisions with the hot ambient electrons can be redistributed between the
species.

When modelling the expansion collisionless and fluid models for ions were used
because of their opposite collisionality regimes; it is expected that the
shared qualitative properties of these expansions hold with a more accurate
model for the ions. The most important qualitative feature of the expansions
is the near-equal split of the heating power between electrons and ions. This
energy balance is in agreement with that of the self-similar solutions to the
expansion found in Aleynikov et al., (2019); Arnold et al., (2021). The
explanation for this result is that self-similar solutions tend to be
‘attractors’ for more complicated systems, with the particularly robust
predictions being those that do not contain reference to any parameters at
all, such as the energy balance (Barenblatt,, 1996). It is reasonable
therefore to suggest that the energy balance holds with a more accurate model
of the ions, perhaps also one including motion transverse to magnetic field
lines, which would allow a description of the assimilation of the plasmoid
into the core of the device.

Since this energy balance entails a considerable transfer of energy from
electrons to ions we conclude that the ambipolar expansion of a pellet
plasmoid is a potent mechanism for the heating of ions on a much faster
timescale than that on which electron-ion collisions occur; the expansion
happens on the hot electron-hot electron collision timescale and the resulting
ion flow energy is converted to thermal energy on the ion-ion collision
timescale, which is approximately $\sqrt{m_{i}/m_{e}}$ times smaller than the
electron-ion collision timescale. Hence fuel pellet injection should be
considered as not just a method for replenishing lost plasma, but also as a
technique for rapidly heating ions if their temperature is exceeded by that of
the electrons.

##  Appendix A Calculation of the $\mathcal{E}$-moment of the electron kinetic
equation on the expansion timescale

The volume element in the variables
$(\mathcal{E}_{\parallel},\mathcal{E}_{\perp},z)$ is given by

$\mathrm{d}^{3}v\,\mathrm{d}z=\frac{4\pi}{m_{e}^{2}}\frac{1}{v_{\parallel}}\,\mathrm{d}\mathcal{E}_{\parallel}\,\mathrm{d}\mathcal{E}_{\perp}\,\mathrm{d}z,$
(A.1)

so phase-space integrals over the trapped domain $V_{t}$ may be expressed as

$\begin{split}\int_{V_{t}}h\,\mathrm{d}v^{3}\,\mathrm{d}z&=\frac{4\pi}{m_{e}^{2}}\int^{\infty}_{-\infty}\int^{\infty}_{0}\int^{0}_{-e\phi}\frac{h}{v_{\parallel}}\,\mathrm{d}\mathcal{E}_{\parallel}\,\mathrm{d}\mathcal{E}_{\perp}\,\mathrm{d}z\\\
&=\frac{2\pi}{m_{e}^{2}}\int^{\infty}_{0}\int^{0}_{-e\phi_{m}}\oint\frac{h}{v_{\parallel}}\,\mathrm{d}z\,\mathrm{d}\mathcal{E}_{\parallel}\,\mathrm{d}\mathcal{E}_{\perp}\\\
&=\frac{2\pi}{m_{e}^{2}}\int^{\infty}_{0}\int^{0}_{-e\phi_{m}}\left<h\right>\tau\,\mathrm{d}\mathcal{E}_{\parallel}\,\mathrm{d}\mathcal{E}_{\perp},\\\
&=\frac{2\pi}{m_{e}^{2}}\int^{\infty}_{0}\int^{J_{m}}_{0}\left<h\right>\,\mathrm{d}J\,\mathrm{d}\mathcal{E}_{\perp}.\end{split}$
(A.2)

Hence, for a term $g$ that is already bounce-averaged, the integral over
phase-space is given by

$\frac{2\pi}{m_{e}^{2}}\int^{\infty}_{0}\int^{0}_{-e\phi_{m}}\tau
g\,\mathrm{d}\mathcal{E}_{\parallel}\,\mathcal{E}_{\perp}=\frac{2\pi}{m_{e}^{2}}\int^{\infty}_{0}\int^{J_{m}}_{0}g\,\mathrm{d}J\mathrm{d}\mathcal{E}_{\perp}.$
(A.3)

We note that Eq. (126) may be expressed as

$\frac{\partial f}{\partial
t}\bigg{|}_{J}=\left<C(f,f_{\mathrm{II}}+f_{\mathrm{III}})\right>,$ (A.4)

where $\cdot|_{J}$ indicates a derivative at constant $J$ rather than constant
$\mathcal{E}_{\parallel}$. Hence its integral over phase-space (after being
multiplied by $\mathcal{E}$) may be expressed as

$\frac{2\pi}{m_{e}^{2}}\int^{\infty}_{0}\int^{J_{m}}_{0}\mathcal{E}\frac{\partial
f}{\partial
t}\bigg{|}_{J}\,\mathrm{d}J\,\mathrm{d}\mathcal{E}_{\perp}=\int_{V_{t}}\mathcal{E}C(f,f_{\mathrm{II}}+f_{\mathrm{III}})\,\mathrm{d}^{3}v\,\mathrm{d}z,$
(A.5)

where we have used the fact

$\frac{2\pi}{m_{e}^{2}}\int^{\infty}_{0}\int^{J_{m}}_{0}\mathcal{E}\left<C(f,f_{\mathrm{II}}+f_{\mathrm{III}})\right>\,\mathrm{d}J\,\mathrm{d}\mathcal{E}_{\perp}=\int_{V_{t}}\mathcal{E}C(f,f_{\mathrm{II}}+f_{\mathrm{III}})\,\mathrm{d}^{3}v\,\mathrm{d}z.$
(A.6)

From Eq. (127) we find

$\begin{split}\frac{\mathrm{d}W_{t}}{\mathrm{d}t}=-\frac{2\pi}{m_{e}^{2}}&\int^{\infty}_{0}\int^{0}_{-e\phi_{m}}\frac{\partial
J}{\partial
t}f\,\mathrm{d}\mathcal{E}_{\parallel}\,\mathrm{d}\mathcal{E}_{\perp}+\frac{2\pi}{m_{e}^{2}}\frac{\partial
J_{m}}{\partial
t}\int^{\infty}_{0}\mathcal{E}_{\perp}f_{a}\big{|}_{\mathcal{E}_{\parallel}=0}\,\mathrm{d}\mathcal{E}_{\perp}\\\
&+\int_{V_{t}}\mathcal{E}C(f,f_{\mathrm{II}}+f_{\mathrm{III}})\,\mathrm{d}^{3}v\,\mathrm{d}z,\end{split}$
(A.7)

where we have used the fact that

$\frac{\partial\mathcal{E}}{\partial
t}\bigg{|}_{J}=-\frac{1}{\tau}\frac{\partial J}{\partial
t}\bigg{|}_{\mathcal{E}_{\parallel}}.$ (A.8)

## Funding

This work has been carried out within the framework of the EUROfusion
Consortium, funded by the European Union via the Euratom Research and Training
Programme (Grant Agreement No 101052200 – EUROfusion). Views and opinions
expressed are however those of the author(s) only and do not necessarily
reflect those of the European Union or the European Commission. Neither the
European Union nor the European Commission and be held responsible for them.
This work was supported by the U.S. Department of Energy under Contract Nos.
DEFG02-04ER54742 and DESC0016283.

## Declaration of interests

The authors report no conflict of interest.

## References

  * Aleynikov et al., (2019) Aleynikov, P., Breizman, B. N., Helander, P., & Turkin, Y. 2019.  Plasma ion heating by cryogenic pellet injection.  Journal of Plasma Physics, 85(905850105).
  * Arnold et al., (2021) Arnold, A. M., Aleynikov, P., & Helander, P. 2021.  Self-similar expansion of a plasmoid supplied by pellet ablation.  Plasma Physics and Controlled Fusion, 63(095008).
  * Arnold et al., (2023) Arnold, A. M., Aleynikov, P., & Breizman, B. N. 2023.  Electron kinetics in a high-$Z$ plasmoid.  Journal of Plasma Physics, 89(905890203).
  * Baldzuhn et al., (2019) Baldzuhn, J., et al. 2019\.  Pellet fueling experiments in Wendelstein 7-X.  Plasma Physics and Controlled Fusion, 61(095012).
  * Baldzuhn et al., (2020) Baldzuhn, J., et al. 2020\.  Enhanced energy confinement after series of pellets in Wendelstein 7-X.  Plasma Physics and Controlled Fusion, 62(055012).
  * Barenblatt, (1996) Barenblatt, G. I. 1996.  Scaling, self-similarity, and intermediate asymptotics. 1 edn.  Cambridge University Press.
  * Bozhenkov et al., (2020) Bozhenkov, S. A., et al. 2020\.  High-performance plasmas after pellet injections in Wendelstein 7-X.  Nuclear Fusion, 60(066011).
  * Helander & Sigmar, (2002) Helander, P., & Sigmar, D. J. 2002.  Collisional Transport in Magnetised Plasmas. 1st edn.  Cambridge University Press.
  * Parks et al., (1977) Parks, P. B., et al. 1977\.  A model for the ablation rate of a solid hydrogen pellet in plasma.  Nuclear Fusion, 17(3), 539–556.
  * Runov et al., (2021) Runov, A., Aleynikov, P., Arnold, A. M., Breizman, B. N., & Helander, P. 2021.  Modelling of parallel dynamics of a pellet-produced plasmoid.  Journal of Plasma Physics, 87(905870407).