Parallel expansion of a fuel pellet plasmoid

Abstract

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.

1. Introduction

During a recent experimental campaign of the W7-X (Wendelstein 7-X) stellarator, fuel pellet injection was found to be associated with a subsequent increase in the ion–electron temperature ratio (Baldzuhn et al. 2019, 2020; Bozhenkov et al. 2020). Such an increase is generally desirable, as a higher ion temperature results in a larger fusion cross-section. It has been conjectured that this effect is due to stabilisation of the ITG modes, resulting in better ion confinement (Xanthopoulos et al. 2020). However, the rapid parallel expansion of the ionised pellet material, the pellet plasmoid, along magnetic field lines has been proposed as at least partially contributing to the increased temperature ratio by directly transferring energy from electrons to ions in the form of significant ion flow velocity, which is eventually converted to heat (Aleynikov et al. 2019; Arnold, Aleynikov & Helander 2021; Runov et al. 2021; Aleynikov et al. 2023).

The aim of this paper is to provide a rigorous model of the parallel plasmoid expansion that resolves the inconsistencies of earlier models. 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 the new approach, which does not suffer from these pitfalls, is provided before being realised mathematically.

When a fuel pellet is injected into a magnetic confinement fusion (MCF) device, the incoming energy flux from the ambient plasma ablates the surface of the pellet and produces a gas cloud (Parks, Turnbull & Foster 1977). The pellet and gas are composed of electrically neutral molecules, but plasma is continually generated within the gas cloud by the collisions of the high-energy ions and electrons composing the multikiloelectronvolt 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 happens 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.

As the plasmoid is a localised increase in density, 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. As 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 electric potential drives parallel dynamics much more quickly than transverse dynamics, the latter being due to drifts and collisional diffusion. Therefore, as with previous investigations, we consider the parallel expansion of the pellet plasmoid on a given field line.

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 in $z$ and monotonically decreasing in $|z|$, with the electron density and potential reaching their maxima $n_m$ and $\phi _m$ at $z = 0$.

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 those produced by pellet injection in a state-of-the-art MCF device. We take as a reference point the W7-X stellarator, as the success of its pellet injection experiments 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 electron density (along a field line) of a fuel pellet plasmoid in W7-X is $N = 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 approximately $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 will 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{\rm \pi} R)(2{\rm \pi} r)/r_I$, as this is the length after which the flux tube of diameter $r_I$ self-intersects.

The transverse size of the pellet plasmoid in W7-X is difficult to estimate from the camera observations (Baldzuhn et al. 2019), and the question of the transverse size of the pellet plasmoid is complicated by the fact that the ablated and ionised material tends to split up into smaller, closely spaced (often overlapping) discrete plasmoids that are seen as ‘striations’ in the density profile, as observed in ASDEX (Axially Symmetric Divertor EXperiment) Upgrade experiments (Mueller et al. 2002). In those experiments, the overall plasmoid split up into tens of plasmoids of size ${\sim }1\ \mathrm {cm}$. Therefore, we consider $r_I \sim 1\unicode{x2013}10\ \mathrm {cm}$, which can be interpreted as the size range of one small plasmoid to the overall plasmoid. This range gives a connection length of $600\unicode{x2013}6000\ \mathrm {m}$. As the parallel size of the plasmoid is expected to reach this after its density has dropped to at most the value of the ambient plasma for $N = 10^{22}\ \mathrm {m}^{-2}$ (Aleynikov et al. 2019), we formally take the connection length to be infinite. This formulation is naturally consistent with the ambient plasma temperature remaining constant. Arnold, Aleynikov & Breizman (2023), in contrast, treated electron kinetics in a high-Z plasmoid on a field line of finite connection length, accounting for the quenching of the ambient plasma temperature.

As plasmoid electrons are ‘born’ at energies comparable to the ionisation energy, of order tens of electronvolts, but are immersed in an ambient plasma with a temperature of the order of several kiloelectronvolts, 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 considered 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, only 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 a treatment of 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 electron dynamics is required to resolve these issues. We will consider electron kinetics in the variables of (i) 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, (ii) perpendicular energy $\mathcal {E}_\perp = m_e v_\perp ^2/2$, where $v_\perp$ is the speed perpendicular to the field line, (iii) $z$, the position along the field line and (iv) $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 (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 problem. 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.

In each region we employ the separation of time scales appropriate for the plasmoid and plasma parameters mentioned earlier in this section in order to obtain a simplified kinetic equation for the electrons. We find that trapped electrons collide with the cold, dense plasmoid electrons (and the plasmoid ions) much more frequently 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 time scales leads to electrons reaching a ‘quasiequilibrium’ (QE) state which is characterised by the electron distribution exhibiting a steady-state on the time scale 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, in order to prevent the ‘charging up’ of the plasmoid (violation of quasineutrality) on this time scale. The QE 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 time. 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 time scale on which QE is established. This allows us to express one parameter of the lowest-order distribution function in terms of the other. Once the electron distribution is known in terms of the remaining parameter, which we choose to be $T$, 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 $T$. The velocity moment corresponding to line-integrated energy density is then taken over the higher-order electron kinetic equation, which is valid on the expansion time scale, in order 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 are 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 much larger line-integrated density $N$ the temperature equilibration would occur well before the densities are comparable. With a much smaller line-integrated density, the densities would 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 plasmoid rapidly expanding into the ambient plasma, causing a shock near the extremities of the plasmoid. This shock would generate sound waves and solitons that propagate into the ambient plasma; wave generation would then sap the kinetic energy of the plasmoid expansion and possibly alter 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, as a deviation from quasineutrality would occur near the shock and within the waves. 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 the hot ambient ions either traverse the entire plasmoid or are rapidly reflected from the potential, hence do not ‘pile up’ at the moving 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 of 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. We develop these estimates assuming that plasmoid ions are singly charged and have mass $m_i$. However, as we wish to model the 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. give 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

(1.1)\begin{equation} \nu_h = \frac{n_a e^4 \ln\varLambda}{6\sqrt{2}{\rm \pi}^{3/2} \varepsilon_0^2 m_e^{1/2}T_a^{3/2}} \end{equation}

is the inverse heating time of the cold plasmoid electrons by a hot population of density $n_a$ and temperature $T_a$. The expression

(1.2)\begin{equation} T \sim T_a(1 - \exp({-\nu_h t})) \end{equation}

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 for both the short and long-term evolution of the plasmoid.

In the self-similar solution, the density becomes infinite at $z = 0$ 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

(1.3)\begin{equation} n_m \sim N \nu_h\sqrt{\frac{3 m_i}{8{\rm \pi} (\nu_h t)^3 T_a}} + n_a \end{equation}

for the peak electron 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

(1.4)\begin{equation} e\phi_m \sim T \ln\left(\frac{n_m}{n_a}\right) \end{equation}

provides an adequate estimate for the height of the potential. This estimate is also supported by the 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 < T_a$.

Equation (1.3) expresses the peak electron density, which we will subsequently use to obtain an ordering. Of course, electrons move within the plasmoid, so the most rigorous approach would be to consider ‘average’ quantities throughout an orbit. This would be very complicated, and relies on detailed knowledge of the shape of the potential which we have not yet obtained. Therefore, we obtain the ordering using 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 average over the trapped orbit of any quantity will be heavily weighted by the value at the peak. When considering passing electrons, we note that the expressions we obtain are the same as those obtained by integrating over the $z$-dependent plasmoid density.

2. Electron kinetics

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

(2.1)\begin{equation} \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) \end{equation}

in the variables $(v_\parallel, v_\perp, z, t)$, where $C(\,f,f)$ is the electron self-collision operator and $C_{e,ik}(\,f)$ is the collision operator for electrons colliding with ion species $k$.

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

(2.2)\begin{equation} \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), \end{equation}

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 \ge 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 equation in regions I, II and III after splitting 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

(2.3)\begin{equation} \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}).\end{equation}

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

(2.4)\begin{equation} \nu_p(v) = \frac{n_m(1 + Z_\mathrm{eff})e^4 \ln\varLambda}{8{\rm \pi}\varepsilon_0^2 m_e^2 v^3}, \end{equation}

where $v$ is the electron speed and

(2.5)\begin{equation} Z_\mathrm{eff} = \frac{\sum_k Z_k^2 n_{ik}(z = 0)}{\sum_k Z_k n_{ik}(z = 0)} \end{equation}

is the effective charge of the ions at $z = 0$, where $Z_k$ and $n_{ik}$ are the charge (in multiples of $e$) and density of ion species $k$, respectively. Equation (2.4) is the frequency at which an electron experiences pitch-angle scattering, assuming that quasineutrality

(2.6)\begin{equation} n_e(z = 0) = n_m = \sum_k Z_k n_{ik}(z = 0), \end{equation}

where $n_e$ is the electron density, holds. 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/n_m$. 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)\begin{equation} \frac{\nu_p(v_{T_a})}{\nu_T} := \mu = \frac{N e^4 \ln\varLambda}{16{\rm \pi}\varepsilon_0^2 T_a^2}.\end{equation}

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

(2.8)\begin{equation} \mu \ll 1, \end{equation}

where for the purpose of calculation the Coulomb logarithm $\ln \varLambda$ 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. Here $\mu$ is independent of any parameters that change during the expansion, hence it is small throughout the entire expansion. The expression for opacity obtained by taking into account the changing plasmoid density is the same as (2.7) (see (6) in Arnold et al. (2021)).

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 time scale than the transit time. Therefore, the shortest time scale in 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.9)\begin{equation} v_\parallel \frac{\partial f_\mathrm{III}}{\partial z} = 0, \end{equation}

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

(2.10)\begin{align} \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}), \end{align}

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

(2.11)\begin{equation} \oint g\,\mathrm{d}z := 2 \int^{z_c}_{-z_c}g\,\mathrm{d}z, \end{equation}

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

(2.12)\begin{equation} \mathcal{E}_\parallel+ e\phi(z_c,t) = 0 \end{equation}

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

(2.13)\begin{equation} \langle g\rangle = \frac{1}{\tau}\oint g\,\mathrm{d}z \end{equation}

for bounce period

(2.14)\begin{equation} \tau = \oint \frac{\mathrm{d}z}{v_\parallel}, \end{equation}

where $v_\parallel = \sqrt {(2/m_e)(\mathcal {E}_\parallel + e\phi (z))}$. We note that on an infinitely long magnetic field line the bounce-average of any function $g$ for $\mathcal {E}_\parallel \ge 0$ is given by

(2.15)\begin{equation} \langle g\rangle|_{\mathcal{E}_\parallel \ge 0} = \lim_{|z| \rightarrow \infty} g. \end{equation}

As the potential vanishes at infinity, and we expect the distribution function to be constant at infinity, the bounce-average of (2.10) is solved by the Maxwellian defining the ambient plasma,

(2.16)\begin{equation} f_\mathrm{III} = f_a = n_a \left(\frac{m_e}{2{\rm \pi} T_a}\right)^{3/2} \exp\left({-\frac{\mathcal{E}}{T_a}}\right). \end{equation}

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

The kinetic equation for electrons in region I is given by

(2.17)\begin{equation} \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}). \end{equation}

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 of height $\phi _m$ and width $L_p$, the bounce frequency is given by

(2.18)\begin{equation} \nu_B \sim \frac{L_p}{\sqrt{2e\phi_m/m_e}} \end{equation}

and is associated with the convective term, so we write

(2.19)\begin{equation} v_\parallel \frac{\partial f_\mathrm{I}}{\partial z} \sim \nu_B\kern0.08em f_\mathrm{I}. \end{equation}

Substituting expressions for the temperature (1.2) and density (1.3) from the modified self-similar solution into the Boltzmann relation (1.4) yields

(2.20)\begin{equation} \frac{e\phi_m}{T_a} \sim (1-\exp({-\nu_h t}))\ln\left(1 + \frac{2}{{\rm \pi}\sqrt{3}} \sqrt{\frac{m_i}{m_e}} \mu (\nu_h t)^{-({3}/{2})}\right). \end{equation}

The above is always order unity given $\mu \lesssim 1$, except at very early or very late times.

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

(2.21)\begin{equation} \nu_B \sim \frac{L_p}{\sqrt{2T_a/m_e}}.\end{equation}

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

(2.22)\begin{equation} C(\,f_\mathrm{I},f_\mathrm{I}) \sim \sum_k C_{e,ik}(\,f_\mathrm{I}) \sim\nu_p(v_T) f_\mathrm{I}, \end{equation}

where $v_T = \sqrt {2T/m_e}$ is the typical electron speed in region I.

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

(2.23)\begin{equation} \frac{1}{T_a}\frac{\mathrm{d} T}{\mathrm{d} t} = \nu_h \exp({-\nu_h t}) \end{equation}

and write

(2.24)\begin{equation} C(\,f_\mathrm{I},f_\mathrm{II} + f_\mathrm{III}) \sim \nu_h \exp({-\nu_h t}) f_\mathrm{I}. \end{equation}

The ratio of the heating rate and the frequency of collisions with the plasmoid (using the modified self-similar temperature (1.2) and density (1.3)) gives

(2.25)\begin{align} \frac{\nu_h \exp({-\nu_h t})} {\nu_p(v_T)} \sim \frac{4}{3\sqrt{\rm \pi}}\exp({-\nu_h t})\left(1 - \exp({-\nu_h t})\right)^{3/2}\left(1 + \frac{2}{{\rm \pi}\sqrt{3}}\sqrt{\frac{m_i}{m_e}}\mu (\nu_h t)^{-({3}/{2})}\right)^{-1}. \end{align}

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

(2.26)\begin{equation} \frac{\nu_h \exp({-\nu_h t})} {\nu_p(v_T)} \ll 1. \end{equation}

As region I represents the cold electrons, we expect the time-dependent terms in (2.17) to act on a much longer time scale 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. (2.21)), remains the same as the collision frequency increases. Therefore, collisions with the plasmoid and bounce motion are associated with the two shortest time scales in region I. Accordingly, the lowest-order kinetic equation in region I is

(2.27)\begin{equation} v_\parallel \frac{\partial f_\mathrm{I}}{\partial z} = C(\,f_\mathrm{I},f_\mathrm{I}) + \sum_k C_{e,ik}(\,f_\mathrm{I}). \end{equation}

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 $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

(2.28)\begin{equation} f_\mathrm{I} = f_0 = \eta \left(\frac{m_e}{2{\rm \pi} T}\right)^{3/2} \exp\left({-\frac{\mathcal{E}}{T}}\right) \end{equation}

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) and Aleynikov et al. (2023), 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 higher-order kinetic equation in region I, corresponding to the heating time scale, is given by

(2.29)\begin{equation} \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}), \end{equation}

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 of temperature $T$ 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 full Maxwellian. However, we must not artificially extend region I past the point where the distribution function would actually be Maxwellian.

From the Boltzmann relation (1.4) and the density of the full Maxwellian distribution (2.28) being $n_0 = \eta \exp (e\phi /T)$ we find that

(2.30)\begin{equation} n_0 \sim n_a \exp\left(\frac{e\phi}{T}\right). \end{equation}

Consider an electron with parallel energy $\mathcal {E}_\mathrm {I/II}$, which has a turning point $z_c$ such that $\mathcal {E}_\mathrm {I/II} + e\phi (z_c) = 0$. At this turning point,

(2.31)\begin{equation} n_0(z_c) \sim n_a \exp\left({-\frac{\mathcal{E}_\mathrm{I/II}}{T}}\right). \end{equation}

During its orbit this electron collides with a plasmoid electron density that is at least as large as the above. Therefore, if we choose $\mathcal {E}_\mathrm {I/II}$ such that $n_0(z_c) > a n_a$ for $a \gg 1$, then the collisions the electron experiences are dominated by collisions with the plasmoid 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. The upper bound for $\mathcal {E}_\mathrm {I/II}$ is correspondingly expressed as

(2.32)\begin{equation} \mathcal{E}_\mathrm{I/II} <-T\ln a. \end{equation}

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

(2.33)\begin{equation} \frac{f_0(\mathcal{E}_\mathrm{I/II})}{f_0(-e\phi_m)} < \frac{1}{a}. \end{equation}

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

(2.34)\begin{equation} \mathcal{E}_\mathrm{I/II} >-e\phi_m + T\ln a. \end{equation}

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

The kinetic equation in region II is given by

(2.35)\begin{equation} \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}). \end{equation}

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 in region I. The typical velocity in region II 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 (2.21). Hence, we write

(2.36)\begin{gather} C(\,f_\mathrm{II},f_\mathrm{I}) \sim \sum_k C_{e,ik}(\,f_\mathrm{II}) \sim \nu_p(v_{T_a})f_\mathrm{II}, \end{gather}

(2.37)\begin{gather}C(\,f_\mathrm{II},f_\mathrm{II} + f_\mathrm{III}) \sim \nu_h \exp({-\nu_h t})f_\mathrm{II}, \end{gather}

(2.38)\begin{gather}v_\parallel \frac{\partial f_\mathrm{II}}{\partial z} \sim \nu_B f_\mathrm{II}, \end{gather}

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

(2.39)\begin{gather} \frac{\nu_p(v_{T_a})}{\nu_B} = \mu \ll 1, \end{gather}

(2.40)\begin{gather}\frac{\nu_h \exp({-\nu_h t})}{\nu_B} = \frac{4}{3\sqrt{\rm \pi}} \frac{n_a L_p \exp({-\nu_h t})}{N} \mu \ll 1. \end{gather}

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

Then, the shortest time scale in region II is the bounce period, which leads to the lowest-order equation

(2.41)\begin{equation} v_\parallel \frac{\partial f_\mathrm{II}}{\partial z} = 0, \end{equation}

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

(2.42)\begin{equation} \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\langle C(\,f_\mathrm{II},f_\mathrm{I}) + \sum_k C_{e,ik}(\,f_\mathrm{II})\right\rangle + \langle C(\,f_\mathrm{II},f_\mathrm{II} + f_\mathrm{III})\rangle \end{equation}

where

(2.43)\begin{equation} 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 \end{equation}

is the second adiabatic invariant for a trapped electron.

Now we analyse the time scales on which the time-dependent terms act and compare them with other time scales. Collisional kinetic problems require continuity of the distribution, so $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 the analysis of the time scales at each boundary. We note that $f_a$ is constant in time, hence at $\mathcal {E}_\parallel = 0$ the first term on the left-hand side of (2.42) vanishes. The second term represents the adiabatic change in electron energy as the well expands, which in Aleynikov et al. (2019) was shown to occur on the heating time scale,

(2.44)\begin{equation} \left|\frac{1}{\tau}\frac{\partial J}{\partial t}\frac{\partial f_\mathrm{II}}{\partial \mathcal{E}_\parallel}\right| \sim \nu_h \exp({-\nu_h t}) f_\mathrm{II}.\end{equation}

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

(2.45)\begin{equation} \left.\frac{\partial f_\mathrm{II}}{\partial t}\right|_{\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}}. \end{equation}

Owing to the logarithmic dependence on the large number $a$ in (2.32), we can always choose $\mathcal {E}_\mathrm {I/II}$ such that $|\mathcal {E}_\mathrm {I/II}| \sim T$. We expect that $(\partial \eta /\partial t)/\eta \sim (\partial T/\partial t)/T$, so (2.45) can be approximated by

(2.46)\begin{equation} \left.\frac{\partial f_\mathrm{II}}{\partial t}\right|_{\mathcal{E} = \mathcal{E}_\mathrm{I/II}} \sim \frac{1}{T} \left.\frac{\mathrm{d} T}{\mathrm{d} t}f_\mathrm{II}\right|_{\mathcal{E} = \mathcal{E}_\mathrm{I/II}}. \end{equation}

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

(2.47)\begin{equation} \nu_t = \frac{1}{T}\frac{\mathrm{d} T}{\mathrm{d} t} = \frac{T_a}{T} \nu_h \exp({-\nu_h t}), \end{equation}

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 with the collision time.

However, the expressions from the modified self-similar solution give

(2.48)\begin{equation} \frac{\nu_t} {\nu_p(v_{T_a})} \sim \frac{4}{3\sqrt{\rm \pi}}\exp({-\nu_h t}) \sqrt{1 - \exp({-\nu_h t})}\left(1 +\frac{2}{{\rm \pi}\sqrt{3}}\sqrt{\frac{m_i}{m_e}} (\nu_h t)^{-({3}/{2})}\right)^{-1}, \end{equation}

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

(2.49)\begin{equation} \frac{\nu_t} {\nu_p(v_{T_a})} \ll 1. \end{equation}

We also find that the ratio of the heating rate to the collision frequency with the plasmoid is much less than unity at all times,

(2.50)\begin{equation} \frac{\nu_h \exp({-\nu_h t})} {\nu_p(v_{T_a})} \sim \frac{4}{3\sqrt{\rm \pi}}\exp({-\nu_h t})\left(1 +\frac{2}{{\rm \pi}\sqrt{3}} \sqrt{\frac{m_i}{m_e}}\mu (\nu_h t)^{-({3}/{2})}\right)^{-1} \ll 1. \end{equation}

So, the collision terms with the plasmoid corresponds to the shortest time scale in (2.42); the lowest-order kinetic equation with respect to (2.42) is therefore

(2.51)\begin{equation} \left\langle C(\,f_\mathrm{II},f_\mathrm{I}) + \sum_k C_{e,ik}(\,f_\mathrm{II})\right\rangle= 0.\end{equation}

This must be solved with boundary conditions ensuring continuity,

(2.52)\begin{gather} f_\mathrm{II}(\mathcal{E} = \mathcal{E}_\mathrm{I/II}) = f_0(\mathcal{E} = \mathcal{E}_\mathrm{I/II}), \end{gather}

(2.53)\begin{gather}f_\mathrm{II}(\mathcal{E}_\parallel= 0) = f_a(\mathcal{E}_\parallel= 0). \end{gather}

The higher-order equation in region II is

(2.54)\begin{equation} \frac{\partial f_\mathrm{II}}{\partial t} - \frac{1}{\tau} \frac{\partial J}{\partial t}\frac{\partial f_\mathrm{II}}{\partial \mathcal{E}_\parallel} = \langle C(\,f_\mathrm{II},f_\mathrm{II} + f_\mathrm{III})\rangle, \end{equation}

which again describes the heating and expansion time scales.

2.5. The QE problem

The distribution function in region II is obtained by solving (2.51), which must be be supplemented by boundary conditions (2.52), (2.53) enforcing continuity of the distribution function into regions I and III. Conceptually, the kinetic problem in region II describes a quasiequilibrium: hot trapped electrons experience rapid collisions against a Maxwellian (and experience rapid pitch-angle scattering against 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)$. As it only remains to solve the kinetic problem in region II, it is unnecessary to have a subscript; we write $f = f_\mathrm {II}$ in region II. Hence, the QE kinetic equation can be written as

(2.55)\begin{equation} \langle C_\mathrm{QE}(\,f)\rangle = 0 \end{equation}

for

(2.56)\begin{equation} C_\mathrm{QE}(\,f) = C(\,f,f_0) + \sum_k C_{e,ik}(\,f). \end{equation}

We note that, owing to the fact that collisions are now linearised in terms of collisions against a full Maxwellian, the lower boundary condition (2.52) 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 on 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 electron density $N = 10^{22}\ \mathrm {m}^{-2}$ in an ambient plasma of electron density $n_a = 5\times 10^{19}$ at a temperature $T_a = 5\ \mathrm {keV}$. The resulting QE formalism demands that during most of the expansion the potential height is of order the ambient temperature,

(2.57)\begin{equation} e\phi_m \sim T_a, \end{equation}

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

(2.58)\begin{equation} \mu = \frac{\nu_p(v_{T_a})}{\nu_T} \sim \frac{\nu_p(v_{T_a})}{\nu_B} \ll 1, \end{equation}

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

(2.59)\begin{gather} \frac{\nu_h \exp({-\nu_h t})}{\nu_p(v_{T_a})} \ll 1, \end{gather}

(2.60)\begin{gather}\frac{\nu_h \exp({-\nu_h t})}{\nu_p(v_T)} \ll 1, \end{gather}

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,

(2.61)\begin{gather} \frac{1}{T}\frac{\mathrm{d} T}{\mathrm{d} t} \ll \nu_p(v_{T_a}), \end{gather}

(2.62)\begin{gather}\frac{1}{T}\frac{\mathrm{d} T}{\mathrm{d} t} \ll \nu_p(v_T). \end{gather}

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 therefore does formally somewhat constrain the value of the line-integrated density, but the potential growing logarithmically with the plasmoid density means that even for quite small line-integrated plasmoid densities the latter condition is met.

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 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 $\partial /\partial t$ term directly in the kinetic equation. In this sense, the formalism is analogous to the Braginskii equations (or any Chapman–Enskog expansion (Chapman & Cowling 1970)), but valid for systems with a long rather than short hot electron mean free path.

We reiterate that the kinetic problem in regions 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 Maxwellian. However, we note that when $T = T_a$, which occurs in the late stage of the expansion, the potential well will be too shallow to contain several $T$. However, at this point the entire electron distribution is a single Maxwellian, so (2.55) will still be satisfied. We therefore conclude that the 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 (2.52) can be applied at $\mathcal {E} = -e\phi _m$ as we approximate collisions with $f_I$ by collisions with the full Maxwellian $f_0$. 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 kinetic equation 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 plasmoid density goes to zero or infinity. Then, the ratio $e\phi _m/T_a$ is a small or large 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 for most of the assimilation process. 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 \approx T_a$ or $n_m \approx n_a$ for most of the assimilation process.

When constructing the ordering, the opacity $\mu$ was given assuming that intraspecies 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 ions slow very rapidly on the electrons; the collision frequency of the ions with 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 causes rapid heating of the plasmoid electrons, 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 QE 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 $\varphi$ of the velocity) is given by

(2.63)\begin{align} C(\,f, f_0) & = \frac{e^4 \ln\varLambda}{4{\rm \pi}\varepsilon_0^2 m_e^2} \left\{\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) \right. \nonumber\\ & \quad \left.+\,\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]\right\}, \end{align}

where $x = v/v_T$, $g(x)$ is the function

(2.64)\begin{equation} g(x) = \int |\boldsymbol{v} - \boldsymbol{v}^\prime| f_0^\prime\,\mathrm{d}^3v^\prime = n_0 v_T\left[\left(x +\frac{1}{2x}\right)\mathrm{erf}(x) + \frac{1}{\sqrt{\rm \pi}}\exp\left(-x^2\right)\right], \end{equation}

and

(2.65)\begin{equation} n_0 = \eta \exp\left({\frac{e\phi}{T}}\right) \end{equation}

is the density of deeply trapped electrons (Helander & Sigmar 2002). We note that when $x$ is large, i.e. we consider collisions of electrons with 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

(2.66)\begin{equation} C_{e,ik}(\,f) = \frac{e^4 \ln\varLambda}{4{\rm \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]. \end{equation}

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

(2.67)\begin{equation} C(\,f,f_0) = \boldsymbol{\nabla}_{\boldsymbol{v}}\boldsymbol{\cdot} \boldsymbol{F}, \end{equation}

so it can always be expressed in the form

(2.68)\begin{equation} C(\,f,f_0) = |J| \boldsymbol{\nabla}_{\boldsymbol{w}}\boldsymbol{\cdot} \tilde{\boldsymbol{F}}, \end{equation}

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

(2.69)\begin{equation} J = \mathrm{det}\left(\frac{\partial w_i}{\partial v_j}\right), \end{equation}

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

(2.70)\begin{equation} \mathrm{det}\left(\frac{\partial (\mathcal{E}_\parallel, \mathcal{E}_\perp,\varphi)}{\partial (\boldsymbol{v})}\right) = m_e^2 v_\parallel \end{equation}

we seek to transform (2.63) into the form

(2.71)\begin{equation} C(\,f,f_0) = A v_\parallel \boldsymbol{\nabla}_{(\mathcal{E}_\parallel, \mathcal{E}_\perp)} \boldsymbol{\cdot} \tilde{\boldsymbol{F}} \end{equation}

for some constant $A$. We find that

(2.72)\begin{align} C(\,f,f_0) & = \frac{e^4 \ln \varLambda}{4{\rm \pi}\varepsilon_0^2 m_e}v_\parallel \boldsymbol{\nabla}_{(\mathcal{E}_\parallel,\mathcal{E}_\perp)} \nonumber\\ & \quad \boldsymbol{\cdot}\left[\,f_0\left(\boldsymbol{r}\boldsymbol{r} \mathcal{E}_\perp \frac{v_\parallel}{v^3}\frac{g^\prime(x)}{v_T} + T \frac{\boldsymbol{s}\boldsymbol{s}}{v^5 v_\parallel} \frac{x^3 g^{\prime\prime}(x)}{v_T}\right) \boldsymbol{\nabla}_{(\mathcal{E}_\parallel,\mathcal{E}_\perp)}\left(\frac{f}{f_0}\right)\right], \end{align}

where

(2.73a,b)\begin{equation} \boldsymbol{r} = (1,-1),\quad \boldsymbol{s} =(v_\parallel^2,v_\perp^2) \end{equation}

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

(2.74)\begin{equation} C_{e,ik}(\,f) = \frac{e^4 \ln \varLambda}{4{\rm \pi}\varepsilon_0^2 m_e}v_\parallel \boldsymbol{\nabla}_{(\mathcal{E}_\parallel,\mathcal{E}_\perp)}\boldsymbol{\cdot} \left[\,f_0\left(\boldsymbol{r}\boldsymbol{r} \mathcal{E}_\perp \frac{v_\parallel}{v^3}Z_k^2 n_{ik}\right)\right]. \end{equation}

Here $\boldsymbol {r}\boldsymbol {r}$ is the tensor associated with pitch-angle scattering, altering parallel and perpendicular energies such that their sum is unchanged. Here $\boldsymbol {s}\boldsymbol {s}$ is associated with collisions that alter energy, affecting parallel and perpendicular energies equally.

We must bounce-average the collision operators (2.72) and (2.74). 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 $\boldsymbol {F}$, we have

(2.75)\begin{equation} \boldsymbol{\nabla}_{(\mathcal{E}_\parallel, \mathcal{E}_\perp)}\boldsymbol{\cdot}\oint \boldsymbol{F}\,\mathrm{d}z = \oint \boldsymbol{\nabla}_{(\mathcal{E}_\parallel, \mathcal{E}_\perp)}\boldsymbol{\cdot} \boldsymbol{F}\,\mathrm{d}z + 4\boldsymbol{\nabla}_{(\mathcal{E}_\parallel, \mathcal{E}_\perp)}z_c \boldsymbol{\cdot} \boldsymbol{F}(z_c). \end{equation}

In order for the divergence and orbit integral to commute, we must have $\boldsymbol {\nabla }_{(\mathcal {E}_\parallel,\mathcal {E}_\perp )}z_c \boldsymbol {\cdot } \boldsymbol {F}(z_c) = 0$; either $\boldsymbol {F}(z_c)$ or $\boldsymbol {\nabla }_{(\mathcal {E}_\parallel,\mathcal {E}_\perp )}z_c$ vanishes, or $\boldsymbol {\nabla }_{(\mathcal {E}_\parallel,\mathcal {E}_\perp )}z_c$ is orthogonal to $\boldsymbol {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 with respect to the pitch-angle scattering term.

As for the energy-altering term, we observe that

(2.76)\begin{equation} \boldsymbol{\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), \end{equation}

and

(2.77)\begin{equation} (\boldsymbol{s}\boldsymbol{s})(z=z_c) =\left(\begin{matrix} 0 & 0\\ 0 & v_\perp^4\\ \end{matrix}\right), \end{equation}

which means that with respect to this term,

(2.78)\begin{equation} \boldsymbol{\nabla}_{(\mathcal{E}_\parallel, \mathcal{E}_\perp)}z_c \boldsymbol{\cdot} \boldsymbol{F}(z_c) \propto\left(\begin{matrix} 1\\ 0 \end{matrix}\right)\boldsymbol{\cdot} \left[\left(\begin{matrix} 0 & 0\\ 0 & v_\perp^4 \end{matrix}\right) \boldsymbol{\nabla}_{(\mathcal{E}_\parallel, \mathcal{E}_\perp)}\left(\frac{f}{f_0}\right)\right] = 0, \end{equation}

so the orbit integral and divergence commute for the energy-altering term. As 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.79)\begin{align} \langle C_\mathrm{QE}(\,f)\rangle & = \frac{e^4 \ln \varLambda}{4{\rm \pi}\varepsilon_0^2 m_e \tau} \boldsymbol{\nabla}_{(\mathcal{E}_\parallel,\mathcal{E}_\perp)} \boldsymbol{\cdot}\left[\,f_0\left(\boldsymbol{r}\boldsymbol{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 \right.\right. \nonumber\\ & \quad \left.\left.\vphantom{\left(\frac{g^\prime(x)}{v_T} + \sum_k Z_k^2 n_{ik}\right)} +T \oint\frac{\boldsymbol{s}\boldsymbol{s}}{v^5 v_\parallel} \frac{x^3 g^{\prime\prime}(x)}{v_T}\mathrm{d}z\right) \boldsymbol{\nabla}_{(\mathcal{E}_\parallel,\mathcal{E}_\perp)}\left(\frac{f}{f_0}\right)\right]. \end{align}

The QE kinetic equation is given by setting the above to zero,

(2.80)\begin{align} & \boldsymbol{\nabla}_{(\mathcal{E}_\parallel,\mathcal{E}_\perp)} \boldsymbol{\cdot}\left[\,f_0\left(\boldsymbol{r}\boldsymbol{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 \right.\right. \nonumber\\ & \quad \left.\left.\vphantom{\left(\frac{g^\prime(x)}{v_T} + \sum_k Z_k^2 n_{ik}\right)} + T \oint\frac{\boldsymbol{s}\boldsymbol{s}}{v^5 v_\parallel} \frac{x^3 g^{\prime\prime}(x)}{v_T}\mathrm{d}z\right) \boldsymbol{\nabla}_{(\mathcal{E}_\parallel,\mathcal{E}_\perp)}\left(\frac{f}{f_0}\right)\right]= 0, \end{align}

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

(2.81)\begin{equation} \boldsymbol{\nabla}_{(\mathcal{E}_\parallel,\mathcal{E}_\perp)}\boldsymbol{\cdot} \left[D_\mathrm{QE}\boldsymbol{\nabla}_{(\mathcal{E}_\parallel,\mathcal{E}_\perp)} \left(\frac{f}{f_0}\right)\right] = 0 \end{equation}

for

(2.82)\begin{equation} D_\mathrm{QE} = D_\mathrm{QE,S} + D_\mathrm{QE,F}, \end{equation}

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

(2.83)\begin{equation} 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), \end{equation}

and that associated with collisions that alter energy is given by

(2.84)\begin{equation} 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. \end{equation}

Together with quasineutrality,

(2.85)\begin{equation} \int f\,\mathrm{d}^3v = n_e = \sum_k Z_k n_{ik}, \end{equation}

(2.81) with the boundary conditions (2.52), (2.53) (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 known 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. In that case, 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 § 3.4, the profile to remain single-peaked when the high temperature of the ambient ions 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 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 combinations of $\eta$ and $T$ are not physically meaningful, as 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 quasineutrality. Therefore, a closer look at quasineutrality during the establishment of the QE state is required.

The global quasineutrality condition is given by

(2.86)\begin{equation} N_t + N_p = \sum_k Z_k N_{ik} \end{equation}

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 Gaussian ion density profile from 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 with the ambient density, and the electric potential at the endpoints is negligible compared with $T$, $T_a$ or $e\phi _m$.

In order to maintain global quasineutrality we require

(2.87)\begin{equation} \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}. \end{equation}

As the density is constant far from 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

(2.88)\begin{equation} \frac{\mathrm{d} N_t}{\mathrm{d} t} = \frac{2{\rm \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, \end{equation}

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

(2.89)\begin{equation} \frac{\mathrm{d} N_p}{\mathrm{d} t} = \frac{\mathrm{d}}{\mathrm{d} t} \int^{L_S/2}_{-({L_S}/{2})} n_a \exp\left({\frac{e\phi}{T_a}}\right) \mathrm{erfc}\left(\sqrt{\frac{e\phi}{T_a}}\right)\mathrm{d}z. \end{equation}

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

(2.90)\begin{equation} \left.\frac{\partial f}{\partial t}\right|_J = \langle C(\,f,f_\mathrm{I})\rangle + \sum_k \langle C_{e,ik}(\,f)\rangle + \langle C(\,f,f_\mathrm{II} + f_\mathrm{III})\rangle, \end{equation}

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

(2.91)\begin{align} \frac{\mathrm{d} N_t}{\mathrm{d} t} & = \frac{2{\rm \pi}}{m_e^2} \int^\infty_0\int^{J_m}_0 \langle C_\mathrm{QE}(\,f)\rangle\,\mathrm{d}J\,\mathrm{d}\mathcal{E}_\perp \nonumber\\ & \quad +\frac{2{\rm \pi}}{m_e^2} \int^\infty_0 \int^{J_m}_0 \langle C(\,f,f_\mathrm{II}+f_\mathrm{III})\rangle\,\mathrm{d}J\,\mathrm{d}\mathcal{E}_\perp \nonumber\\ & \quad +\frac{2{\rm \pi}}{m_e^2} \int^\infty_0 \frac{\mathrm{d} J_m}{\mathrm{d} t} f_a(\mathcal{E} = \mathcal{E}_\perp)\,\mathrm{d}\mathcal{E}_\perp, \end{align}

where we have used the fact that $f(J = J_m) = f(\mathcal {E}_\parallel = 0) = f_a(\mathcal {E}_\parallel = 0)$. We then see that (2.87) becomes

(2.92)\begin{align} & \frac{2{\rm \pi}}{m_e^2} \int^\infty_0 \int^{J_m}_0 \langle C_\mathrm{QE}(\,f)\rangle\,\mathrm{d}J\,\mathrm{d}\mathcal{E}_\perp \nonumber\\ & \quad +\frac{2{\rm \pi}}{m_e^2} \int^\infty_0 \int^{J_m}_0 \langle C(\,f,f_\mathrm{II}+f_\mathrm{III})\rangle\,\mathrm{d}J\,\mathrm{d}\mathcal{E}_\perp \nonumber\\ & \quad +\int^{{L_S}/{2}}_{-({L_S}/{2})}n_a \frac{e}{T_a} \frac{\partial \phi}{\partial t}\exp\left({\frac{e\phi}{T_a}}\right) \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{align}

where the term associated with $\mathrm {d} J_m/\mathrm {d} 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

(2.93)\begin{equation} \frac{2{\rm \pi}}{m_e^2} \int^\infty_0 \int^{J_m}_0 \langle C_\mathrm{QE}(\,f)\rangle\,\mathrm{d}J\,\mathrm{d}\mathcal{E}_\perp\sim \nu_p(v_{T_a})N, \end{equation}

where we have used the fact that $N \approx N_t$. The second term on the left-hand side is due to heating; we may write

(2.94)\begin{equation} \frac{2{\rm \pi}}{m_e^2} \int^\infty_0 \int^{J_m}_0 \langle C(\,f,f_\mathrm{II}+f_\mathrm{III})\rangle\,\mathrm{d}J\,\mathrm{d}\mathcal{E}_\perp\sim \nu_h \exp({-\nu_h t}) N. \end{equation}

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 $f$ always being equal to $f_a$ in the passing region. We can estimate this term by noting that $\exp (e\phi /T_a)\mathrm {erfc}(\sqrt {e\phi /T_a}) \sim 1$ for $\phi \lesssim T_a$ and using the approximation

(2.95)\begin{equation} \int^{{L_S}/{2}}_{-({L_S}/{2})} \frac{e\phi}{T_a}\mathrm{d}z \sim L_p \frac{e\phi_m}{T_a} \sim L_p. \end{equation}

Writing $L_p = N/n_m$ and using (1.3) then yields

(2.96)\begin{equation} \int^{{L_S}/{2}}_{-({L_S}/{2})} n_a \frac{e}{T_a} \frac{\partial \phi}{\partial t}\exp\left({\frac{e\phi}{T_a}}\right) \mathrm{erfc}\left(\sqrt{\frac{e\phi}{T_a}}\right)\mathrm{d}z \sim \nu_h N \frac{\mu \sqrt{\dfrac{3m_i}{m_e} \nu_h t}}{\left((\nu_h t)^{3/2} + \dfrac{2}{{\rm \pi}\sqrt{3}}\sqrt{\dfrac{m_i}{m_e}}\mu\right)^2}. \end{equation}

With the plasma parameters used in the Introduction, the above is at most of order $\nu _h N$, and decreases as time advances; it is of the same order as the heating term (the second on the left-hand side of (2.92)).

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 are at most of the same order as (2.96). 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 (2.92) is the first on the left-hand side, which is associated with the frequency at which an electron collides with the plasmoid; the other terms are associated with the heating frequency. Hence (2.92) is well-approximated by a vanishing net electron flux associated with collisions with the plasmoid, which can be expressed as

(2.97)\begin{equation} \varGamma = \frac{e^4 \ln \varLambda}{2\varepsilon_0^2 m_e^3} \int^\infty_0 \left.\left[D_\mathrm{QE}\boldsymbol{\nabla}_{(\mathcal{E}_\parallel,\mathcal{E}_\perp)} \left(\frac{f}{f_0}\right)\right]\right|_{\mathcal{E}_\parallel= 0}\boldsymbol{\cdot} \left(\begin{matrix} 1 \\ 0 \end{matrix}\right)\mathrm{d}\mathcal{E}_\perp= 0.\end{equation}

Hence, to maintain global quasineutrality, there can (approximately) be no net electron 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.

As we had two free parameters in the QE problem, $\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 (2.81) 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 (2.81) with boundary conditions (2.52), (2.53) (noting that we write $f = f_\mathrm {II}$ in region II) was solved numerically with a self-consistent potential given by quasineutrality (2.85) and the no-net-flux condition (2.97). The plasma was assumed to be hydrogenic: there was a single species of singly charged ion following the profile

(2.98)\begin{equation} n_i = n_a + N\frac{1}{L_p\sqrt{\rm \pi}}\exp\left({-\left(\frac{z}{L_p}\right)^2}\right), \end{equation}

where $L_p = 2.8\ \mathrm {m}$, $N = 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 self-similar profile in Aleynikov et al. (2019) at $t = 20\,\mathrm {\mu }\mathrm {s}$ given these parameters.

Figure 3 shows the properties of the electron distribution function. Figure 3(a) shows the distribution function in velocity space at $z = 0$. Here $\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-$v_\perp$ region of trapped phase space is observed. This is a result of the friction experienced by hot electrons causing them to fall to lower energies after being scattered into this region. Figure 3(c) 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-hand boundary is the trapped–passing separatrix. Figure 3(b) shows the effective phase-space flow velocity $\boldsymbol {u}^*$, defined via

(2.99)\begin{equation} \boldsymbol{\nabla}_{(\mathcal{E}_\parallel,\mathcal{E}_\perp)} \boldsymbol{\cdot} (\boldsymbol{u}^*f) =-\langle C_\mathrm{QE}(\,f)\rangle, \end{equation}

i.e.

(2.100)\begin{equation} \boldsymbol{u}^* \propto-\frac{D_\mathrm{QE} \boldsymbol{\nabla}_{(\mathcal{E}_\parallel,\mathcal{E}_\perp)} \left(\frac{f}{f_0}\right)}{f}. \end{equation}

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

Figure 3. (a) Numerical distribution function at $z = 0$ in SI units. (b) Effective phase-space flow velocity $\boldsymbol {u}^*$ (see (2.100)). (c) Phase-space trajectories of electrons (the streamlines of $\boldsymbol {u}^*$). (d) Collisional flux into the trapped region (see (2.101), (2.102)). Here $v_c = \sqrt {2e\phi _m/m_e}$ is the parallel escape velocity at $z = 0$.

Figure 3(d) shows the flux through the trapped–passing separatrix, where

(2.101)\begin{gather} \varGamma_S =-\frac{e^4 \ln \varLambda}{2\varepsilon_0^2 m_e^3}\left.\left[D_\mathrm{QE, S} \boldsymbol{\nabla}_{(\mathcal{E}_\parallel,\mathcal{E}_\perp)} \left(\frac{f}{f_0}\right)\right]\right|_{\mathcal{E}_\parallel= 0}\boldsymbol{\cdot} \left(\begin{matrix} 1 \\ 0 \end{matrix}\right), \end{gather}

(2.102)\begin{gather}\varGamma_F =-\frac{e^4 \ln \varLambda}{2\varepsilon_0^2 m_e^3}\left.\left[D_\mathrm{QE, F} \boldsymbol{\nabla}_{(\mathcal{E}_\parallel,\mathcal{E}_\perp)} \left(\frac{f}{f_0}\right)\right]\right|_{\mathcal{E}_\parallel= 0}\boldsymbol{\cdot}\left(\begin{matrix} 1 \\ 0 \end{matrix}\right), \end{gather}

and $\varGamma = \varGamma _S + \varGamma _F$. Collisions with the cold Maxwellian always produce an inflow of electrons in the $\mathcal {E} > 0$ region due to friction (see $\varGamma _F$ in figure 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 cause an inflow at higher energies. The net flux through the separatrix vanishes due to the no-net-flux condition.

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. As 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 to the solution of 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 = Z n_i$), from (2.63), (2.66) we see that the bounce-averaged QE collision operator in a square well is given by

(2.103)\begin{equation} \langle C_\mathrm{QE}(\,f)\rangle = \frac{n_0 e^4 \ln\varLambda}{4{\rm \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]. \end{equation}

We consider $\langle C_\mathrm {QE}(\,f)\rangle = 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

(2.104)\begin{equation} \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}. \end{equation}

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,

(2.105)\begin{equation} \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}, \end{equation}

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

(2.106)\begin{equation} f_k = C_k \exp\left({-\frac{1}{2}\lambda_k v^2}\right) \cosh\left(v_\parallel\sqrt{\frac{2}{1+Z}\lambda_k}\right) \end{equation}

for constants $\{C_k\}$ and $\{\lambda _k\}$. The boundary condition (2.53) is satisfied by the sum of two solutions

(2.107)\begin{equation} f_1 = f_a(\mathcal{E}) \frac{\cosh\left(\dfrac{v_\parallel}{v_{T_a}} \sqrt{\dfrac{4}{1+Z}}\right)}{\cosh\left(\dfrac{v_c}{v_{T_a}} \sqrt{\dfrac{4}{1+Z}}\right)} - f_0(\mathcal{E}) \frac{\cosh\left(\dfrac{v_\parallel}{v_{T}}\sqrt{\dfrac{4}{1+Z}}\right)}{\cosh \left(\dfrac{v_c}{v_{T}}\sqrt{\dfrac{4}{1+Z}}\right)}. \end{equation}

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

(2.108)\begin{equation} \mathcal{G}(\,f_1) = \frac{1 + Z}{2}\frac{\partial}{\partial \xi} \left[(1 - \xi^2)\frac{\partial f_1}{\partial \xi}\right] + v \frac{\partial f_1}{\partial v}. \end{equation}

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

(2.109)\begin{equation} f_1 = \sum_{n=0}^\infty a_n(v) P_n(\xi), \end{equation}

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

(2.110)\begin{equation} v \frac{\partial a_n}{\partial v} - \frac{1 + Z}{2}n(n+1) a_n = 0, \end{equation}

which has solutions

(2.111)\begin{equation} a_n = c_n v^{(({1 + Z})/{2})n(n+1)} \end{equation}

for $\{c_n\}$ constants. The continuity of the distribution function at $v = v_c$ provides expressions for $c_n$,

(2.112)\begin{align} c_n & = \frac{2n + 1}{2}v_c^{-({(1+Z)}/{2})n(n+1)} \left[\frac{f_a(\mathcal{E}=0)}{\mathrm{cosh} \left(\dfrac{v_c}{v_{T_a}}\sqrt{\dfrac{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 \right.\nonumber\\ & \quad \left.-\frac{f_0(\mathcal{E}=0)}{\mathrm{cosh}\left(\dfrac{v_c}{v_{T}} \sqrt{\dfrac{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\right]. \end{align}

To summarise, the analytical solution to the QE problem in a square well is given by (2.107) in the $\mathcal {E} > 0$, $\mathcal {E}_\parallel < 0$ region, (2.109) in the $\mathcal {E} \le 0$ region and $f = f_a$ in the $\mathcal {E}_\parallel \ge 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 figure 3. Figure 4(a) shows the distribution in velocity space at $z = 0$. Figure 4(b) shows the percentage difference from the numerical solution given in figure 3. Figure 4(c) shows the distribution function at $\mathcal {E}_\perp = 0$. Figure 4(d) 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-$v_\perp$ region of trapped phase space. The observation that in this region the contours of the distribution function are nearly horizontal lines in $(v_\parallel,v_\perp )$ can be explained by the fact that for the analytical QE distribution function

(2.113)\begin{equation} \frac{\partial}{\partial v_\parallel}f(v_\parallel, v_\perp) \propto \frac{2}{v_{T_a}^2} v_\parallel\left(\frac{2}{1 + Z} - 1\right) + O(v_\parallel^3), \end{equation}

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

Figure 4. (a) Analytical electron distribution $f_\mathrm {an}$ in SI units at $z = 0$. (b) Percentage difference between the numerical $f$ (see figure 3) and analytical $f_\mathrm {an}$. (c) Distributions at $\mathcal {E}_\perp = 0$. (d) Distributions at $\mathcal {E}_\parallel = -e\phi _m$. Here $v_c = \sqrt {2e\phi _m/m_e}$ is the parallel escape velocity at $z = 0$.

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$ region with the continuity boundary condition (2.53) then fixes the value of $f$ at $\mathcal {E} = 0$. Similarly, in the $\mathcal {E} < 0$ region we only have the opportunity to apply the continuity boundary condition at $\mathcal {E} = 0$, but not at $\mathcal {E} = -e\phi _m$. As a consequence, as $\mathcal {E} \rightarrow -e\phi _m$, $f_1 \rightarrow c_0$, which violates (2.52). On the other hand, $c_0$ is several orders of magnitude smaller than $f_0(-e\phi _m)$, so the discrepancy 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 $\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 (2.85) is satisfied. Additionally, the no-net-flux condition (2.97) 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 the QE problem, 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 (2.97). Pitch-angle scattering acts to diffuse newly trapped electrons to lower parallel velocities while leaving their energy unchanged. Simultaneously the scattered electrons lose energy via friction. These effects can be seen in the streamline plot of figure 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 (2.97), 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. As 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 (2.97) requires us to integrate along $\mathcal {E}_\parallel = 0$, where the ejection structure strongly affects the flux. As the analytical solution to the QE problem 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 the QE problem in a square well.

As 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 term 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 (2.103) we see that the collisional flux in the $\hat {\boldsymbol {v}}$ direction at $v = v_c$ is given by

(2.114)\begin{equation} F =-\frac{n_0 e^4 \ln\varLambda}{4{\rm \pi}\varepsilon_0^2 m_e^2} \left.\frac{1}{v_c^2}\left(\,f + \frac{T}{m_e v}\frac{\partial f}{\partial v}\right)\right|_{\mathcal{E} = 0}. \end{equation}

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

(2.115)\begin{equation} G =-\int^{2{\rm \pi}}_0 \int^1_{-1} v_c^2(F|_{\mathcal{E} = 0^+} + F|_{\mathcal{E} = 0^-})\,\mathrm{d}\xi\,\mathrm{d}\varphi. \end{equation}

Setting the above to zero will provide the relation between $\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

(2.116)\begin{equation} \eta = n_a \left(\frac{T}{T_a}\right)^{3/2}\left(2 - \frac{T}{T_a}\right). \end{equation}

Now, quasineutrality allows us to express both $\eta$ and $\phi$ in terms of $T$ and the ion profiles. 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 (2.85) demands that

(2.117)\begin{equation} n_e(\phi,t;T) = \sum_k Z_k n_{ik}, \end{equation}

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

(2.118)\begin{equation} n_e(z = 0) = \eta \exp\left(\frac{e\phi_m}{T}\right), \end{equation}

so

(2.119)\begin{equation} 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]. \end{equation}

This is a modification of the Boltzmann relation, the extra contributions accounting for the fact that in the QE state the distribution function resembles a cold Maxwellian of temperature $T$ for $\mathcal {E} < 0$ and a hot Maxwellian of temperature $T_a$ for $\mathcal {E} > 0$. As expected, the above reduces to the Boltzmann relation for $T = T_a$ (when the distribution becomes a single Maxwellian), agreeing with the rough estimate (1.4) which was used to justify the ordering leading to the QE problem. As noted earlier, $e\phi _m$ given by (2.119) exceeds that given by the Boltzmann relation when $T < T_a$.

The height of the potential in the numerical solution to the QE problem (figure 3), which is calculated entirely self-consistently, and leads to no net flux into the trapped region, is within 3 % of the estimate (2.119). 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.

3. Plasmoid expansion

3.1. Treatment of ions

We consider two different models for the ions: a collisionless (Vlasov) system with hot ambient ions, but cold plasmoid ions; and a cold-fluid system. We restrict our attention to plasmoids with 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 as 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 for ambient electrons: the opacity (2.7). The collisionless system still models cold plasmoid ions accurately as they are initialised with zero velocity spread.

For the collisionless ion expansion we solve the Vlasov equation,

(3.1)\begin{equation} \frac{\partial f_i}{\partial t} + v_\parallel \frac{\partial f_i}{\partial z} - \frac{Z e}{m_i}\frac{\partial \phi}{\partial z}\frac{\partial f_i}{\partial v_\parallel} = 0. \end{equation}

For the cold-fluid expansion we solve the equations

(3.2)\begin{gather} \frac{\partial n_i}{\partial t} + \frac{\partial}{\partial z}(n_i u_i) = 0, \end{gather}

(3.3)\begin{gather}\frac{\partial u_i}{\partial t} + u_i \frac{\partial u_i}{\partial z} + \frac{Z e}{m_i}\frac{\partial \phi}{\partial z} = 0. \end{gather}

3.2. Treatment of electrons on the expansion time scale

Equations (2.29) and (2.54) describe electron dynamics on the expansion and heating time scales. We note that because they take exactly the same form after bounce-averaging, (2.54) and the bounce-average of (2.29) may be combined to yield

(3.4)\begin{equation} \frac{\partial f}{\partial t} - \frac{1}{\tau}\frac{\partial J}{\partial t} \frac{\partial f}{\partial \mathcal{E}_\parallel} = \langle C(\,f,f_\mathrm{II} + f_\mathrm{III})\rangle \end{equation}

to describe both regions I and II.

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 a moment of the kinetic equation (3.4) over phase space. When we take moments, the QE 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 as 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

(3.5)\begin{equation} W_t = \int_{V_t} \mathcal{E} f\,\mathrm{d}^3v\,\mathrm{d}z, \end{equation}

where $V_t$ is the trapped region of phase space. Taking the $\mathcal {E}$-moment of (3.4) yields

(3.6)\begin{equation} \frac{\mathrm{d} W_t}{\mathrm{d} t} = \left.\frac{\mathrm{d} W_t}{\mathrm{d} t}\right|_\mathrm{adiabatic} + \left.\frac{\mathrm{d} W_t}{\mathrm{d} t}\right|_\mathrm{separatrix} + \left.\frac{\mathrm{d} W_t}{\mathrm{d} t}\right|_\mathrm{heating} \end{equation}

for

(3.7)\begin{gather} \left.\frac{\mathrm{d} W_{t}}{\mathrm{d} t}\right|_\mathrm{adiabatic} =-\frac{2{\rm \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, \end{gather}

(3.8)\begin{gather}\left.\frac{\mathrm{d} W_{t}}{\mathrm{d} t}\right|_\mathrm{separatrix} = \frac{2{\rm \pi}}{m_e^2}\frac{\mathrm{d} J_m}{\mathrm{d} t} \left.\int^\infty_0 \mathcal{E}_\perp f_a\right|_{\mathcal{E}_\parallel= 0}\,\mathrm{d}\mathcal{E}_\perp, \end{gather}

(3.9)\begin{gather}\frac{\mathrm{d} W_{t}}{\mathrm{d} t}\bigg|_\mathrm{heating} = \int_{V_t} \mathcal{E} C(\,f,f_\mathrm{II} + f_\mathrm{III})\,\mathrm{d}^3v\,\mathrm{d}z. \end{gather}

The details of the procedure are contained in Appendix A.

The terms in (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\phi _m$; ‘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 time step. 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. Here $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.10)\begin{equation} K_t = W_t + \int^{L_S/2}_{-({L_S}/{2})} e \phi n_{e,t}\,\mathrm{d}z \end{equation}

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

(3.11)\begin{align} \frac{\mathrm{d} K_t}{\mathrm{d} t} & = \frac{\mathrm{d} W_t}{\mathrm{d}t} + \int^{L_S/2}_{-({L_S}/{2})} e \frac{\partial \phi}{\partial t}n_{e,t}\,\mathrm{d}z + \int^{L_S/2}_{-({L_S}/{2})} e\phi \frac{\partial n_{e,t}}{\partial t}\,\mathrm{d}z \nonumber\\ & = \left.\frac{\mathrm{d} W_t}{\mathrm{d}t} - \frac{\mathrm{d} W_t}{\mathrm{d} t}\right|_\mathrm{adiabatic} + \int^{{L_S}/{2}}_{-({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{align}

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

(3.12)\begin{equation} \int^{L_S/2}_{-({L_S}/{2})} Z e \phi \frac{\partial n_i}{\partial t}\mathrm{d}z =-\frac{\partial K_i}{\partial t} \end{equation}

for $K_i$ the line-integrated kinetic energy of the ions. This can be derived by taking a velocity moment of (3.1) or by constructing an energy equation from (3.2), (3.3). The same procedure cannot be carried out for the passing electrons as 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)\begin{equation} \oint g(\mathcal{E}_\parallel> 0)\,\mathrm{d}z = 2 \int^{L_S/2}_{-({L_S}/{2})} g(\mathcal{E}_\parallel> 0)\,\mathrm{d}z, \end{equation}

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

(3.14)\begin{equation} W_p = \frac{2{\rm \pi}}{m_e^2}\int^\infty_0\int^\infty_{J_m} \mathcal{E} f_a\,\mathrm{d}J\,\mathrm{d} \mathcal{E}_\perp, \end{equation}

hence

(3.15)\begin{equation} \frac{\mathrm{d} W_p}{\mathrm{d} t} = \frac{2{\rm \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{\rm \pi}}{m_e^2}\int^\infty_0 \mathcal{E}_\perp \frac{\mathrm{d} J_m}{\mathrm{d} t}f_a(\mathcal{E} = \mathcal{E}_\perp)\,\mathrm{d}\mathcal{E}_\perp. \end{equation}

We note that

(3.16)\begin{equation} \frac{2{\rm \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^{L_S/2}_{-({L_S}/{2})} e \frac{\partial \phi}{\partial t}n_{e,p}\,\mathrm{d}z \end{equation}

and

(3.17)\begin{equation} \frac{2{\rm \pi}}{m_e^2}\int^\infty_0 \mathcal{E}_\perp \frac{\mathrm{d} J_m}{\mathrm{d} t}f_a(\mathcal{E} = \mathcal{E}_\perp) \,\mathrm{d}\mathcal{E}_\perp= \left.\frac{\mathrm{d} W_t}{\mathrm{d} t}\right|_\mathrm{separatrix}, \end{equation}

so, given that

(3.18)\begin{equation} K_p = W_p + \int^{L_S/2}_{-({L_S}/{2})} e\phi n_{e,p}\,\mathrm{d}z, \end{equation}

we obtain the energy conservation law on the interval $z \in [-L_S/2,L_S/2]$

(3.19)\begin{equation} \frac{\mathrm{d}}{\mathrm{d} t}(K_t + K_p + K_i) = \left.\frac{\mathrm{d} W_t}{\mathrm{d} t}\right|_\mathrm{heating} + \frac{2{\rm \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. \end{equation}

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|_\mathrm {adiabatic}$, as 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 of that expected for cold electrons in a homogeneous plasma, as 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 approximately $3\nu _h(T_a-T)$ (Aleynikov et al. 2019), the per-electron heating rate for electrons trapped in the potential well was found to be $(9/4) \nu _h(T_a-T)$. The QE distribution function is highly isotropic for $\mathcal {E} < 0$, so we approximate the collisional heating term by

(3.20)\begin{equation} \left.\frac{\mathrm{d} W_t}{\mathrm{d} t}\right|_\mathrm{heating} = \frac{9}{4} \nu_h N_{\mathcal{E}<0}(T_a - T), \end{equation}

where $N_{\mathcal {E}<0}$ is the line-integrated density of trapped electrons with $\mathcal {E} < 0$. As shown in § 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

(3.21)\begin{equation} \frac{\mathrm{d}}{\mathrm{d}t}(K_e + K_i) = Q \end{equation}

for heating power

(3.22)\begin{equation} Q = \frac{9}{4}\nu_hN_{\mathcal{E}<0}(T_a - T) + \frac{2{\rm \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. \end{equation}

3.3. Comparison of the system with earlier models

At this point a clear comparison can be drawn between this investigation and Aleynikov et al. (2019), Arnold et al. (2021) and Aleynikov et al. (2023). In those publications 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 rate $3\nu _h T_a$ for plasmoid electrons (only the $T \ll T_a$ stage of the expansion was treated). The potential was given by the Boltzmann relation and was unbounded as $|z| \rightarrow \infty$ due to the neglecting of passing electrons.

In this paper, on the other hand, we have an energy conservation law with heating terms and an electron kinetic energy derived from the solution to a rigorously derived electron kinetic problem. The electric potential is given by the quasineutrality condition, which now accounts for passing electrons, hence the potential vanishes as $|z| \rightarrow \infty$.

As an alternative to the more recent work on plasmoid expansion, Rozhanskij & Veselova (1994) used cold-fluid equations for ions and an energy transport equation for electrons. The energy transport equation considers a single electron temperature which may depend on position, accounts for decompressive cooling of the electrons (which corresponds to the adiabatic extraction of energy from electrons as the well expands), and attempts to account for the long mean free path of electrons with a non-local heat conductivity derived from an approximate electron kinetic problem (Luciani, Mora & Virmont 1983). However, this formulation is still fundamentally limited by the electron mean free path and does not account for the trapped and passing electrons populations with distinct dynamics: for a plasmoid almost completely transparent to hot electrons it is not appropriate. However, the model still has utility in describing the late-stage expansion and the formation of an ion front as it does account for the ambient plasma density.

In summary, a full kinetic treatment such as the QE formalism would seem necessary to simultaneously account for the long electron mean free path, the strikingly non-Maxwellian nature of the electron distribution function (facilitated by the separation of phase space into trapped and passing regions), and the non-negligible effect of the ambient plasma density associated with the passing electrons.

3.4. Numerical solutions to the plasmoid expansion system

Numerical solutions to the system created by coupling the energy conservation law (3.21), quasineutrality (2.117) and one of systems describing the ions (the Vlasov equation (3.1) or the cold-ion system (3.2), (3.3)) were obtained with the plasma parameters $n_a = 5\times 10^{19}\ \mathrm {m}^{-3}$, $T_a = 5\ \mathrm {keV}$ and $N = 10^{22}\ \mathrm {m}^{-2}$, the same as in § 2.9. The plasma was once more assumed to be hydrogenic. The heating time scale with these parameters is $\nu _h^{-1} = 162\,\mathrm {\mu } \mathrm {s}$ and the expansions were run to $300\,\mathrm {\mu } \mathrm {s}$, at which point $T$ is nearly equal to $T_a$ and the peak 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 § 2.10 was used to calculate densities and kinetic energy densities, respectively, appearing in the quasineutrality condition (2.117) and the energy conservation law (3.21). We neglected any deviation of $f$ from $f_0$ in the $\mathcal {E}<0$ region when evaluating these densities as the difference is negligible. However, we did take into account the fact that $f$ is flattened in the $\mathcal {E} > 0$, $\mathcal {E}_\parallel < 0$ region, as this has a relatively large impact on the density and energy density. We used (2.116) as the expression for $\eta$ in terms of $T$ when calculating densities and energy densities. Additionally, we expanded the $\mathrm {cosh}$ functions found in the analytical expression for $f$ to second order to obtain analytical expressions for the densities and energy densities.

The electron density is then given by

(3.23)\begin{equation} 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} \end{equation}

for

(3.24)\begin{gather} n_{\mathcal{E}<0} = \eta \left(\exp\left(\frac{e\phi}{T}\right) \mathrm{erf}\left(\sqrt{\frac{e\phi}{T}}\right) - \frac{2}{\sqrt{\rm \pi}} \sqrt{\frac{e\phi}{T}}\right), \end{gather}

(3.25)\begin{gather}n^{\mathcal{E}>0,c}_{\mathcal{E}_\parallel< 0} = \frac{2}{\sqrt{\rm \pi}}\eta\left[\sqrt{\frac{e\phi}{T}} - \frac{1}{1 + \dfrac{2}{1+Z}\dfrac{e\phi_m}{T}}\left(\sqrt{\frac{e\phi}{T}} + \frac{1}{3}\frac{2}{1+Z}\left(\frac{e\phi}{T}\right)^{3/2}\right)\right], \end{gather}

(3.26)\begin{gather}n^{\mathcal{E}>0,a}_{\mathcal{E}_\parallel< 0} = \frac{2}{\sqrt{\rm \pi}}n_a \frac{1}{1 + \dfrac{2}{1+Z}\dfrac{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)^{3/2}\right) \end{gather}

and

(3.27)\begin{equation} n_{\mathcal{E}_\parallel> 0} = n_a \exp\left(\frac{e\phi}{T_a}\right) \mathrm{erfc}\left(\sqrt{\frac{e\phi}{T_a}}\right). \end{equation}

The subdensities on the right-hand side of (3.23) are associated with the regions of phase space indicated in their subscripts and superscripts, with $a$ (for ‘ambient’) and $c$ (for ‘cold’), respectively, indicating that they are associated with the first or second term on the right-hand side of (2.107)

The line-integrated electron kinetic energy is given by

(3.28)\begin{equation} 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}, \end{equation}

where

(3.29)\begin{gather} K_{\mathcal{E}<0} = \frac{3}{2}N_{\mathcal{E}<0}T -\frac{2}{\sqrt{\rm \pi}} \eta\int^{L_S/2}_{-({L_S}/{2})}\left(\frac{e\phi}{T}\right)^{3/2}\mathrm{d}z, \end{gather}

(3.30)\begin{gather}K^{\mathcal{E}>0,c}_{\mathcal{E}_\parallel< 0} = N^{\mathcal{E}>0,c}_{\mathcal{E}_\parallel< 0}T + \int^{L_S/2}_{-({L_S}/{2})}e\phi n^{\mathcal{E}>0,c}_{\mathcal{E}_\parallel< 0}\,\mathrm{d}z, \end{gather}

(3.31)\begin{gather}K^{\mathcal{E}>0,a}_{\mathcal{E}_\parallel< 0} = N^{\mathcal{E}>0,a}_{\mathcal{E}_\parallel< 0}T_a + \int^{L_S/2}_{-({L_S}/{2})}e\phi n^{\mathcal{E}>0,a}_{\mathcal{E}_\parallel< 0}\,\mathrm{d}z \end{gather}

and

(3.32)\begin{equation} K_{\mathcal{E}_\parallel> 0} = \frac{3}{2}N_{\mathcal{E}_\parallel> 0}T_a + \frac{1}{\sqrt{\rm \pi}}n_a T_a \int^{L_S/2}_{-({L_S}/{2})} \sqrt{\frac{e\phi}{T_a}}\mathrm{d}z, \end{equation}

where

(3.33)\begin{gather} N_{\mathcal{E}<0} = \int^{L_S/2}_{-({L_S}/{2})}n_{\mathcal{E}<0}\,\mathrm{d}z, \end{gather}

(3.34)\begin{gather}N^{\mathcal{E}>0,c}_{\mathcal{E}_\parallel< 0} = \int^{L_S/2}_{-({L_S}/{2})} n^{\mathcal{E}>0,c}_{\mathcal{E}_\parallel< 0}\,\mathrm{d}z, \end{gather}

(3.35)\begin{gather}N^{\mathcal{E}>0,a}_{\mathcal{E}_\parallel< 0} = \int^{L_S/2}_{-({L_S}/{2})} n^{\mathcal{E}>0,a}_{\mathcal{E}_\parallel< 0}\,\mathrm{d}z, \end{gather}

(3.36)\begin{gather}N_{\mathcal{E}_\parallel> 0} = \int^{L_S/2}_{-({L_S}/{2})}n_{\mathcal{E}_\parallel> 0}\,\mathrm{d}z. \end{gather}

The subscripts and superscripts carry the same meaning as for the densities.

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 collisionless ion expansion and the cold-fluid expansion, respectively. 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 solitons would likely develop near the ion front and propagate into the ambient plasma if Poisson's equation were used to account for the expected deviation from quasineutrality near the front.

Figure 5. Ion distribution function (in SI units) 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. (a) The relative amounts of energy deposited into the electrons and ions. (b) The plasmoid electron temperature $T$ and the estimated temperature evolution given a homogeneous plasma. (c,d) Plots of the electric potential and electron density at various times.

Figure 7. Derived quantities of the cold-fluid ion expansion. (a) The relative amounts of energy deposited into the electrons and ions. (b) The plasmoid electron temperature $T$ and the estimated temperature evolution given a homogeneous plasma. (c,d) Plots of the electric potential and electron density at various times.

In contrast, the collisionless expansion exhibits no steep front owing to the fact that the hot ambient ions do not ‘pile-up’ on the moving plasmoid; the majority of ambient ions have a large enough parallel velocity to either pass over the plasmoid or be reflected from the front before experiencing a significant change in energy. In both cases the electron temperature approaches $T_a$ somewhat more quickly than the estimate (1.2), a consequence of the term in (3.22) corresponding to the energy lost by passing electrons but gained by trapped electrons and ions. This term constitutes the majority of the heating when $T \sim T_a$.

The ion front in the cold-fluid expansion is also found in Rozhanskij & Veselova (1994) and may well have been observed experimentally in the TFTR (Tokamak Fusion Test Reactor) tokamak (Mansfield et al. 1991); in that experiment two symmetric density ‘pulses’ were observed to move along magnetic field lines in opposite directions shortly after pellet injection. As discussed in the previous paragraph and in the Introduction, the ion front is more likely to develop with a lower ambient plasma temperature owing to the ‘piling up’ on the edge of the expanding plasmoid; the front was observed in Rozhanskij & Veselova (1994) with cold ions and in Mansfield et al. (1991) with an ambient ion temperature which is likely to be somewhat lower than the electron temperature of $2\unicode{x2013}2.5\ \mathrm {keV}$.

Naturally, the derivative of the electric potential decreases as the expansion proceeds, owing to the plasmoid density decreasing relative to the ambient density; the expansion velocity will not increase indefinitely, reaching instead some maximum as the electric potential offers increasingly weak acceleration. Correspondingly, the plasmoid electrons do not gain energy indefinitely, instead their temperature is tending to the ambient $T_a$. These limited electron and ion energies are in contrast to the models in Aleynikov et al. (2019), Arnold et al. (2021) and Aleynikov et al. (2023), which have an unbounded electric potential as $|z| \rightarrow \infty$; the potential offers unlimited acceleration to the extremities of the plasmoid and, owing to the form of the heating term, the plasmoid electrons do not reach a maximum temperature.

Of particular note, and ultimately what is sought in this investigation, is the energy balance between electrons and ions, plotted in figures 6(a) and 7(a). The two lines represent the fraction of the heating energy ($\int ^t_0Q\,\mathrm {d}t$ for $Q$ in (3.22)) 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 cold-fluid and Vlasov models. This energy balance is similar to those predicted by the self-similar expansions in Aleynikov et al. (2019), Arnold et al. (2021) and Aleynikov et al. (2023).

We stress that only plasmoid electrons are directly heated by collisions in these expansions; plasmoid ions acquire energy in the form of flow velocity due to acceleration by the electric potential, hence the physical interpretation of this energy balance is that energy extracted adiabatically from trapped plasmoid electrons (which themselves gain energy through heating by ambient electrons) is transferred to the ions.

This energy transfer mechanism occurs on the expansion time scale, i.e. the ambient electron–ambient electron collision time $\nu _h^{-1}$ ($162\,\mathrm {\mu } \mathrm {s}$ with the parameters used here). However, the ion flow energy will eventually be converted to heat on the ambient ion–ambient ion collision time $\nu _h^{-1} \sqrt {m_i/m_e}$ ($7\ \mathrm {ms}$). Ordinarily, ions acquire energy from electrons on the ambient electron–ambient ion collision time, $\nu _h^{-1}(m_i/m_e)$ ($300\ \mathrm {ms}$), hence the process of parallel expansion speeds up energy transfer from electrons to ions by a factor of $\sqrt {m_i/m_e}$ (${\approx }40$).

4. Discussion and conclusions

A model for the parallel expansion of the plasmoid produced by fuel pellet ablation has been developed with a particular focus on rigorously deriving a simplified kinetic problem for the electrons. The result is a relatively simple steady-state electron kinetic problem on the shortest collisional time scale. The long-term evolution of the electrons is described by an energy conservation law obtained by taking a velocity moment of the electron kinetic equation which is valid on the expansion time scale. Ions are described with cold-fluid equations or with the Vlasov equation.

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 for 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 ‘first principles’ manner.

It has been shown that during the expansion, electrons reach a ‘quasiequilibrium’ state: a dynamical steady-state on the fastest collisional time scale 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 with a numerical solution. A robust estimate of the height of the self-consistent electric potential that supports QE 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 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 time scale; with some manipulation this evolution equation can take 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 cold-fluid models were used for ions due to 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. Ions, in contrast to electrons, gain their energy in the form of flow velocity rather than heat owing to their rapid acceleration by the electric potential. This energy balance is in agreement with that of the self-similar solutions to the expansion found in Aleynikov et al. (2019) and 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 reference 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 plasma.

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

Section 6 of Aleynikov et al. (2019) contains a detailed discussion of the net effect of this energy transfer, particularly its competition with collisional thermalisation between ions and electrons. It also includes transport simulations of W7-X that incorporate the energy transfer resulting from exactly 50 % of the heating energy deposited in plasmoid electrons going to the ions. It is, however, difficult to estimate the true magnitude of the energy transfer with only a model of the parallel expansion; different field lines in the transverse profile of the plasmoid, with different line-integrated plasmoid densities, will exhibit different energy transfers (a fact discussed in more detail later in this section), so detailed knowledge of the transverse profile of the plasmoid is required to make an accurate prediction.

From a more general point of view, the QE formalism is essentially the analogue of the Braginskii equations which is appropriate for electrons with a long mean free path in the presence of a symmetric, positively oriented electric potential with a single peak. The QE distribution is the analogue of the near-Maxwellian distribution of the Braginskii equations, which is similarly derived from a steady-state kinetic equation. Therefore, the use of the QE formalism for electrons rather than the Braginskii equations should be considered in codes that attempt to model pellet assimilation in MCF devices.

The QE formalism is more complicated than earlier models of plasmoid expansion in Aleynikov et al. (2019), Arnold et al. (2021), Aleynikov et al. (2023) and Rozhanskij & Veselova (1994), but it has significant advantages. Firstly, the long mean free path of electrons and non-Maxwellian nature of the electron distribution function (with its trapped–passing dynamics) call into question the validity of the model given in Rozhanskij & Veselova (1994). Secondly, the QE formalism is valid over the entire assimilation process, an advantage over the remaining models which are fundamentally limited by the requirements $T \ll T_a$ and $n_m \gg n_a$, which are bound to be violated at late times.

Although not presented in detail in this paper, other simulations of the expansion with the QE formalism indicate that the magnitude of the transfer of energy to ions can change with the line-integrated plasmoid density, as opposed to it remaining constant at approximately 50 % in the models in Aleynikov et al. (2019) and Arnold et al. (2021). Indeed, even with $N = 10^{22}\ \mathrm {m}^{-2}$, which was used in the numerical simulations in § 3.4, the energy balance is closer to 60/40 electrons/ions rather than an equal split; it appears that the energy balance tends to an exact 50/50 split for larger line-integrated densities and a less favourable split for ions at lower line-integrated densities. This can be easily understood as the electric potential providing weaker acceleration to ions when the plasmoid density approaches the ambient density in the QE formalism, whereas in the other models the potential is unbounded and no such reduction in acceleration is seen. The natural conclusion is that the QE formalism is necessary for correctly describing the energy transfer at lower line-integrated densities, such as those found in the extremities of the transverse profile of the pellet plasmoid.