Parallel expansion of a fuel pellet plasmoid

The problem of the assimilation of a cryogenic fuel pellet injected into a hot plasma is considered. Due to the transparency to ambient particles of the plasmoid, the localised region of high-density plasma created by ionisation of the ablated pellet material, electrons reach a ‘quasiequilibrium’ (QE) state which is characterised by a steady-state on the fastest collisional time scale. The simplified electron kinetic equation of the QE state is solved. Taking a velocity moment of the higher-order electron kinetic equation, which is valid on the expansion time scale, permits a fluid closure, yielding an evolution equation for the macroscopic 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 with the size of density perturbations, and permits an anisotropic and highly non-Maxwellian distribution function. As the QE 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 nearly equal amounts 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 time scale shorter than that of ion–electron thermalisation.


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(Baldzuhn et al. , 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 electronion 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 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×10 19 m −3 at a temperature T a = 5 keV.A typical line-integrated density along the field line of a fuel pellet plasmoid in W7-X is N p = 10 22 m −2 (Arnold et al. 2021).
The fuel pellets in W7-X contain approximately 10 20 electrons and penetrate roughly 0.1 m into the plasma, resulting in the average line-density in the radial direction of 10 21 m −1 (Baldzuhn et al. 2019).In W7-X the flux surface located at minor radius r = 0.3 m, given the major radius R = 5 m, has a flux-surface integrated density of 3 × 10 21 m −1 if the density of this flux surface is 5×10 19 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πR)(2πr)/r I , since this is the length after which the flux tube of diameter r I selfintersects.For W7-X pellets r I ≈ 0.1 m shortly after injection (Baldzuhn et al. 2019;Arnold et al. 2021), giving a maximum connection length of 600 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 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 eV, but are immersed in an ambient plasma with a temperature on the order of several 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 , where v ∥ is the velocity parallel to the field line, perpendicular energy E ⊥ = m e v 2 ⊥ /2, where v ⊥ 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 electronelectron 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 ϕ 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 = ν h T a t for t ≪ ν −1 h , where is the inverse heating time of the cold plasmoid electrons by a hot population of density n a and temperature T a .The expression 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 → 0 and vanishes for t → ∞, neglecting the fact that the electron density approaches n a as time advances.Therefore the expression 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) et al. 2023).Example trapped (with turning points ±zc) 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 nm and ϕm at z = 0.
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 ≪ T a .Equation (1.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.

Electron kinetics
The kinetic equation for the electron distribution function f is given by in the variables (v ∥ , v ⊥ , z, t), where v ∥ is the velocity parallel to the field line, v ⊥ 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 showing that collisionless change in electron energy is associated with time-variation of the electric potential.We note that with a stationary potential both E ∥ and E ⊥ are constants of motion in the absence of collisions.Passing electrons have E ∥ ⩾ 0 and trapped E ∥ < 0. The minimum parallel energy of an electron is E ∥ = −eϕ m , which corresponds to an electron with v ∥ = 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 I , f II , and f III .That is, when we are in region I, for example, we solve the kinetic equation for f I .No matter where we are in phase-space, collisions are experienced with every other region of phase-space via C(•, f I + f II + f III ).In this sense it is understood that f I (or f II or f III ) is the entire distribution function in region I (or II or III), but is zero outside of its own region.

Solving the kinetic equation for passing electrons (region III)
The kinetic equation for electrons in region III is given by The collision frequency of an electron inside the plasmoid with the plasmoid electrons and ions is approximately given by where v is the electron speed, n p is the plasmoid electron density, and is the effective charge of the ions.Equation (2.4) is derived from the frequency with which an electron experiences pitch-angle scattering, assuming that quasineutrality 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 Ta = 2T a /m e .The inverse of the time taken for a passing electron to transit the plasmoid is ν T = v Ta /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 (2.7) With the plasma and plasmoid parameters given in the Introduction this ratio is much less than unity: (2.8) where for the purpose of calculation the Coulomb logarithm ln Λ 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 µ as the opacity of the plasmoid.µ is independent of any parameters that change during the expansion, hence it is small throughout the entire expansion.
The terms in the kinetic equation (2.3) 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  (2.11)where z c (E ∥ , t) > 0 is the turning point such that (2.12) (cf.Fig. 1).The bounce average of g is given by for bounce period (2.14) We note that on an infinitely long magnetic field line the bounce-average of any function g for E ∥ > 0 is given by (2.15) Since the potential vanishes at infinity and we expect the distribution function to be constant at infinity, the bounce average of Eq (2.10) is solved by the Maxwellian defining the ambient plasma: e − E Ta . (2.16) We note that, as mentioned earlier, we have used the expression for peak plasmoid density to deduce the ordering (2.8) and obtain Eq. (2.16), despite the fact that plasmoid electrons spend little time inside the plasmoid.The same result can be obtained by a slightly different approach: Eq. (2.7) 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| → ∞, resulting in the distribution identical to that at infnity; hence Eq. (2.16).

Solving the kinetic equation for deeply trapped electrons (region I)
The kinetic equation for electrons in region I is given by (2.17) 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 ϕ m and width L p , the bounce frequency is given by and is associated with the convective term, so we write Substituting temperature (1.2) and density (1.3) from the modified self-similar solution into the Boltzmann relation (1.4) yields (2.20) Given that (2/(π √ 3)) m i /m e µ is at least order unity, the above is order unity for several heating times ν −1 h , 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ϕ m ∼ T a for smaller values of µ.
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: (2.21) Since we expect f I to correspond to a dense population of cold electrons, we associate the collision terms against f I with the frequency of collisions with the plasmoid; we write where v T = 2T /m e the typical electron speed within region I. Since f II and f 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 and write Comparing the frequency of collisions with the plasmoid with the heating using the modified self-similar temperature (1.2) and density (1.3) gives 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 Since region I represents the cold electrons, we expect the time-dependent terms in Eq. (2.17) 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. (2.21)), 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 (2.27) When T ≪ T a , T ≪ eϕ m , which allows us to define E 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 for parameters η(t), T (t) (Aleynikov et al. 2019).We see now that in the earlier work which captures the collisionless change in electron energy due to the expanding well and the heating of the cold Maxwellian by the hot electrons.

Choosing E I/II
Owing to the ordering, when T ≪ T a we understand that the potential well is deep enough for a Maxwellian to reside in region I, provided E I/II is chosen close enough to zero.We now decide upon an explicit definition for E 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 (1.4) and the density of the full Maxwellian distribution (2.28) being n p = η exp(eϕ/T ) we find that (2.30) Considering an electron with with parallel energy E I/II , it has a turning point z c such that eϕ(z c ) + E I/II = 0. Therefore, at this turning point, (2.31) During its orbit, this electron collides with a plasmoid density that strictly larger than the above.Therefore, if we choose E I/II such that n p (z c ) > an a for a ≫ 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 E I/II is expressed as (2.32) 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 for a ≫ 1.Then, the lower bound for E I/II is given by (2.34) 2.4.Deriving the kinetic equation for hot trapped electrons (region II) The kinetic equation in region II is given by 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 2eϕ m /m e ∼ v Ta , 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.(2.21)).Hence we write noting that both the collision frequency with the plasmoid and the heating rate are much smaller than the bounce frequency: 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 meaning that f II (and hence f as a whole) is independent of z.The higher-order kinetic equation can then be bounce-averaged, yielding where 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 II is equal to f 0 at E = E I/II and equal to f a at E ∥ = 0 it serves to perform this analysis at each boundary.At the trapped-passing separatrix, E ∥ = 0, f II must be equal to f a , which is constant in time; the first term on the left hand side of Eq. (2.42) 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: (2.44) At E = E I/II Eq. (2.44) also holds.However, since the Maxwellian in region I has a temperature that changes in time, the time derivative of f II is not zero here: That is, the timescale on which the term acts, which we define via the frequency is the time taken for T to increase by a factor of e.When T is small this can be a very short time, so ν t cannot simply be assumed to be small compared to the collision time.
The profiles from the modified self-similar solution give which is always much less than unity for the plasma parameters given in the Introduction; we write So, the collision term with the with the plasmoid corresponds to the shortest timescale in Eq. (2.42); the lowest-order kinetic equation is therefore (2.50) The distribution function must continuous in collisional kinetic problems, hence Eq. (2.50) must be solved with boundary conditions ensuring continuity: (2.52) The higher-order equation in region II is which describes the heating and expansion timescales.

The quasi-equilibrium problem
The distribution function in region II is obtained by solving Eq. (2.50), which must be be supplemented by boundary conditions (2.51),(2.52)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 trappedpassing separatrix.
When T ≪ T a collisions with the distribution in region I are well-approximated by collisions with the full Maxwellian: C(•, f I ) ≈ C(•, 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 II in region II.Hence the QE kinetic equation can be written as (2.55) Further, owing to the fact that collisions are linearised in terms of collisions against a full Maxwellian, the lower boundary condition Eq. (2.51) can actually be applied at E = −eϕ m rather than E I/II .

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 m −2 in an ambient plasma of density electron n a = 5 × 10 19 at a temperature of 5 keV.The requirements of the ordering are that during most of the expansion the potential height is of order the ambient temperature: (2.56) that the plasmoid is transparent to passing and hot trapped electrons: that the heating rate is much lower than the collision frequency with the plasmoid: and that the time taken for the plasmoid electron temperature to increase by a factor of e is much larger than the collision time: (2.61) The transparency of the plasmoid is dependent upon the line-integrated plasmoid density not being too large, and eϕ m being of order T a is dependent upon the lineintegrated 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 ≪ T a , which permitted a choice of E 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. (2.54) 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 ≪ T a and when T → T a .
Therefore the formalism can actually accurately model the expansion outside of its formal ordering (which has T ≪ 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 (2.51) can be applied at E = −eϕ m in the QE problem.This solves the problem of how to choose E 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 E I/II = −eϕ 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ϕ 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 ≈ n a from an early stage.
When constructing the ordering, the opacity µ 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 (E ∥ , E ⊥ ) The collision operator against f 0 in the variables (v, θ, z, t) for pitch-angle θ (assuming that f is independent of the azimuthal angle of the velocity φ) is given by (2.62) is the function and 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 ′ (x) and x 3 g ′′ (x) are well-approximated by n 0 v T .
Similarly, assuming that collisions with ions are well-approximated by pitch-angle scattering, we have for ion charge Z k and density n ik .The collision operator is given by the divergence of the collisional flux F in velocity space: (2.66) so it can always be expressed in the form where |J| is the Jacobian of the transformation between coordinates v and w, and F is the collisional flux in w phase-space.Noting that det we seek to transform Eq. (2.62) into the form for some constant A. We find that and rr, ss represent dyadic products.Similarly, rr is the tensor associated with pitch-angle scattering, altering parallel and perpendicular energies such that their sum is unchanged.ss is associated with energy exchange, which affects parallel and perpendicular energies equally.
We must bounce-average the collision operators (2.71) and (2.73).In order to obtain a 'useful' expression, we must be able to commute divergence in (E ∥ , E ⊥ ) and the orbit integral.The following Lemma is useful with regards to commuting the divergence and orbit integral: for a vector-valued function F, we have (2.74) In order for the divergence and orbit integral to commute, we must have We observe that the term associated with pitch-angle scattering is proportional to v ∥ , 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 which means that with respect to this term, 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 (2.78) The quasi-equilibrium equation is given by setting the above to zero: which is in the form of an anisotropic steady-state diffusion problem in (E ∥ , E ⊥ ) space: for where the diffusion tensor associated with pitch-angle scattering is given by and that associated with energy exchange is given by (2.83) Together with quasineutrality, Eq. (2.80) with the boundary conditions (2.51),(2.52)(noting that we write f = f II ) provides a unique solution for f and ϕ in terms of the parameters η 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 singlepeaked; 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.

The no-net-flux condition
Given some η and T we may solve the QE problem as specified in the previous subsection.However, most of the combinations of η 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 η 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 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 ∈ [−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ϕ m .In order to maintain global quasineutrality we require Since the density is constant far the plasmoid, the terms in the above cease to change as L S → ∞; 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 lineintegrated density of trapped electrons is given by where J m = J(E ∥ = 0) is the maximum value of the second adiabatic invariant for a trapped electron, and, knowing that f = f a in the E ∥ > 0 region, the time derivative of the line-integrated density of passing electrons is given by n a e eϕ Ta erfc eϕ T a dz. (2.88) The full kinetic equation (2.2) can be bounce-averaged and the left hand side changed to the independent variables (J, E ⊥ , t) to yield (2.90) where we have used the fact that the distribution function is continuous: where the term associated with ∂J m /∂t has cancelled out between dN t /dt and dN p /dt.The first term on the left hand side is associated with fluxes due to collisions with the plasmoid; we may write (2.92) The second term on the left hand side is due to heating; we may write (2.93) 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ϕ/T a )erfc( eϕ/T a ) ⩽ 1 for ϕ ⩾ 0 and using the approximation (2.94) Writing L p = N p /n m and using Eq.(1.3) then yields (2.95) With the plasma parameters used in the Introduction, the above is at most of order ν 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. (2.91)).
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.(2.95).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. (2.91) 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. (2.91) is well-approximated with simply a vanishing flux associated with collisions with the plasmoid, which can be expressed as 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, η 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. (2.80) 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.

Numerical solution to the QE problem
The QE problem (2.80) with boundary conditions (2.51), (2.52) (noting that we write f = f II in region II), a self-consistent potential given by quasineutrality (2.84), and the zero-net-flux condition (2.96) 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 where L p = 2.8 m, N ic = 10 22 m −2 , T = 615 eV, n a = 5 × 10 19 m −3 , and T a = 5 keV.The Gaussian term in the above is consistent with the profile in Aleynikov et al. (2019) at t = 20 µ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. 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 (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 (E ∥ , E ⊥ ) space, clearly indicating flow into and out of the trapped region.The dashed line indicates E = 0.The right boundary is the trapped-passing separatrix.The top right plot shows the effective phase-space flow velocity u * , defined via (2.99) Electron phase-space flow for 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 and Γ = Γ S + Γ F .Collisions with the cold Maxwellian always produces an inflow of electrons in region II due to friction (see Γ F in Fig. 3).Pitch-angle scattering may only cause a flow along lines of constant 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.

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 ϕ m , the distribution function obtained as a solution to the QE problem in a square-well potential of height ϕ 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. (2.62),(2.65)we see that the bounce-averaged QE collision operator in a square well is given by We consider ⟨C QE (f )⟩ = 0 separately in the E < 0 (i.e.v < v c = 2eϕ m /m e in a square well) and 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 E > 0 region we use the variables (v, v ∥ = v cos θ), meaning we must solve D(f 1 ) = 0 for It is convenient to consider the limit where v ≫ v ∥ , which represents the correct limit in the majority of E > 0, E ∥ < 0 phase-space: which has superposable solutions (that are even in v ∥ ) provided by the separation of variables: for constants {C k } and {λ k }.The boundary condition (2.52) is satisfied by the sum of two solutions: . (2.106) In the E < 0 region we use the variables variables (v, ξ = cos θ).So, we must solve G(f 1 ) = 0 where (2.107) We note that Legendre polynomials are the eigenfunctions of the operator in ξ, so we write f in this basis: which gives the following equations for {a n (v)}: with the solutions (2.110) for {c n } constants.The continuity of the distribution function at v = v c provides the expressions for c n : (2.111) To summarise, the analytical solution to the QE problem in a square well is given by Eq. ( 2.106) in the E > 0, E ∥ < 0 region, Eq. (2.108) in the E ⩽ 0 region, and f = f a in the E ∥ ⩾ 0 region.The phase-space domain of the QE problem is the same given any potential, so the substitution of (2/m e )(E + eϕ m ) for v and (2/m e )(E ∥ + eϕ m ) for v ∥ 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 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 E ⊥ = 0.The bottom right plot shows the distribution function at E ∥ = −eϕ 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 ∥ , v ⊥ ) can be explained by the fact that for the analytical QE distribution function (2.112) 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 E > 0 part of region II with the continuity boundary condition then fixes the value of f at E = 0. Similarly, in the E < 0 part of region II we only have the opportunity to apply the continuity boundary condition at E = 0, but not at E = E I/II .As a consequence, as E → E I/II , f 1 → 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.
2.11.Expressions for η and ϕ 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 η and T .The electric potential ϕ is as yet unknown, but must be such that quasineutrality (2.84) is satisfied.Additionally, the no-net-flux condition (2.96) must be satisfied.Therefore, we may reduce the number of unknowns in the system from three (η, T , and ϕ) 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 E ∥ = E ⊥ = 0 characterising the flux of electrons out of trapped phase-space, so the no-netflux condition may not be used directly.Instead, analysis of the flow of electrons in phasespace allows us to formulate a condition that is approximately equivalent to Eq. (2.96).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 E < 0 region.In order for there to be no net flux into the trapped region, the same number of electrons entering the 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 (2.96), which states that there is no net flux into the trapped region of phase-space, is that there is no net flux into the  E < 0 region of phase-space.Since the flux into and out of this region is characterised by integrals on the E = 0 line it is not necessary to have a description of the ejection structure, whereas directly applying Eq. (2.96) requires us to integrate along E ∥ = 0, where the ejection structure strongly affects the flux.Since the analytical solution to QE potentially has a discontinuous derivative at E = 0, the flux into the E < 0 region has contributions from E → 0 + and E → 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, ξ, φ).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 E = 0.In a square well, the number of electrons entering the E < 0 region of phase-space is equal to the number of electrons entering the v < v c region of phase space.From Eq. (2.102) we see that the collisional flux in the v direction at v = v c (assuming that v c ≫ v T ) is given by (2.113) The phase-space area element on the v = v c sphere is given by v 2 c dξ dφ, so the net flux into the v < v c region is given by (2.114) Setting the above to zero will provide the relation η, T , and ϕ required for no net flux into the E < 0 region.It will be seen that the expression obtained from G = 0 in the limit Z → ∞ is quite accurate, even for Z = 1.In this limit the analytical solution is given by f = f a for E > 0 and f = c 0 + f 0 for E < 0, where c 0 = f a (E = 0) − f 0 (E = 0).Then, G = 0 yields the relation (2.115) Now, quasineutrality allows us to express both η and ϕ in terms of T and the ion profile.In the variables (E ∥ , E ⊥ , z, t) the electron distribution is independent of z, so the electron density is a function of ϕ, time, and the parameter T (η being already expressed in terms of T via the above).Therefore the quasineutrality condition (2.84) demands that which gives an implicit expression for ϕ 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 so 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 (1.4) which was used to justify the ordering leading to QE.As noted earlier, eϕ m given by (2.118) 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 (2.118).This is a remarkable agreement considering that the estimate was derived in the limit Z → ∞, 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.

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 (2.7).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: For the cold-fluid expansion we solve the equations

Treatment of electrons on expansion timescale
Equations (2.29) and (2.53) describe the electron dynamics on the expansion and heating timescales.We note that because they take exactly the same form, Eq. (2.53) and the bounce average of Eq. (2.29) may be combined to yield 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 η 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 E-moment over trapped phase-space.The line-integrated energy density of trapped electrons is given by where V t is the trapped region of phase-space. (3.9) The details of the procedure are contained in Appendix A.
The terms in Eq. (3.6) 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ϕ 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ϕ/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 ∈ [−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  (3.11)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 for K i the line-integrated kinetic energy of the ions.This can be derived by taking velocity moments Eq. (3.1) or by constructing an energy equation from Eqs. (3.2),(3.3).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: (3.13) allowing us to define the second adiabatic invariant J for positive parallel energies.The line-integrated passing electron energy density on the interval z ∈ [−L S /2, L S /2] is given by Evidently, the inclusion of the ion and passing electron energies in the energy conservation law accounts for the absence of dW t /dt 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 ≪ T a , but provides a considerable fraction of the heating power when T ∼ 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ν h (T a − T ), the per-electron heating rate for electrons trapped in the potential well was found to be (9/4)ν h (T a − T ).In our case the distribution function is also highly isotropic for E < 0, so we approximate the collisional heating term by where N E<0 is the line-integrated density of trapped electrons with E < 0. As shown in Section 2.9, the distribution function is somewhat less than f a in the 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 for heating power  2021) a coldfluid 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ν 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| → ∞ 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| → ∞.Taking the limit n a → 0 in Eq. (2.84) and the energy conservation law (3.21)recovers the same system as Aleynikov et al. (2019) with the exception of the 3/4 factor on the collisional heating term.

Numerical solutions of the plasmoid expansion system
Numerical solutions to the system created by coupling the energy conservation law (3.21),quasineutrality (2.116), and one of systems describing the ions (the Vlasov equation (3.1) or the cold-ion system (3.2),(3.3))were obtained using the plasma parameters n a = 5 × 10 19 m −3 , T a = 5 keV and N ic = 10 22 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 ν −1 h = 162 µs and the expansions were run to 300 µ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 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 (2.116) and the energy conservation law (3.21).We neglect any deviation of f from f 0 in the 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 E > 0, E ∥ < 0 region, as this will have a relatively large impact on the density and energy density.We use Eq.(2.115) as the expression for η in terms of T when evaluating f 0 .Additionally, we expand the 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 The line-integrated kinetic energies are given by where and where (3.36) 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 (1.2), a consequence of the term in Eq. (3.22) 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 ∼ 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 ( t 0 Q dt for Q in Eq. (3.21)) 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 lineintegrated 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.

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 coldfluid 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 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.where we have used the fact that .8)

Figure 1 .
Figure1.Schematic of the electric potential induced by the presence of the plasmoid(Arnold et al. 2023).Example trapped (with turning points ±zc) 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 nm and ϕm at z = 0.

Figure 2 .
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 E = E ∥ + E ⊥ = 0.The diagonal dashed line indicates E = E I/II , which separates regions I and II.The vertical dashed line indicates E ∥ = 0, the trapped-passing separatrix.

Figure 4 .
Figure 4. Top left: analytical electron distribution fan in SI units.Top right: percentage difference between the numerical f (see Fig. 3) and analytical fan.Bottom left: distributions at E ⊥ = 0. Bottom right: distributions at E ∥ = −eϕm.vc = 2ϕm/me is the parallel escape velocity at z = 0.
.22) 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. (

Figure 5 .Figure 6 .Figure 7 .
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).

E
Appendix A. Calculation of the E-moment of the electron kinetic equation on the expansion timescaleThe volume element in the variables (E ∥ , E ⊥ , z) is given byd 3 v dz = 4π m 2 e 1 v ∥ dE ∥ dE ⊥ dz, (A.1)so phase-space integrals over the trapped domain V t may be expressed asVt a term g that is already bounce-averaged, the integral over phase-space is given f,f II + f III )⟩ , (A.4)where •| J indicates a derivative at constant J rather than constant E ∥ .Hence its integral over phase-space (after being multiplied by E) may be expresseddE ⊥ = Vt EC(f, f II + f III ) d 3 v dz, (A.5)where we have used the fact2π f, f II + f III )⟩ dJ dE ⊥ = Vt EC(f, f II + f III ) d 3 v dz.⊥ f a E ∥ =0 dE ⊥ + Vt EC(f, f II + f III ) d 3 v dz,(A.7) Arnold et al. (2021)19);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.
(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