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. (Reference Arnold, Aleynikov and Helander2021)).

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. Reference Arnold, Aleynikov and Breizman2023). 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. Reference Aleynikov, Breizman, Helander and Turkin2019). We see now that in the earlier work Aleynikov et al. (Reference Aleynikov, Breizman, Helander and Turkin2019), Runov et al. (Reference Runov, Aleynikov, Arnold, Breizman and Helander2021), Arnold et al. (Reference Arnold, Aleynikov and Helander2021) and Aleynikov et al. (Reference Aleynikov, Arnold, Breizman, Helander and Runov2023), 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. (Reference Aleynikov, Breizman, Helander and Turkin2019) 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. (Reference Aleynikov, Breizman, Helander and Turkin2019), 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 Reference Chapman and Cowling1970)), 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 Reference Helander and Sigmar2002). 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 Reference Helander and Sigmar2002). 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. (Reference Aleynikov, Breizman, Helander and Turkin2019)). 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. (Reference Aleynikov, Breizman, Helander and Turkin2019) 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. (Reference Arnold, Aleynikov and Breizman2023), which is strong evidence that the QE state was established there.

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. (Reference Arnold, Aleynikov and Breizman2023), 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. Reference Aleynikov, Breizman, Helander and Turkin2019), 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. (Reference Aleynikov, Breizman, Helander and Turkin2019), Arnold et al. (Reference Arnold, Aleynikov and Helander2021) and Aleynikov et al. (Reference Aleynikov, Arnold, Breizman, Helander and Runov2023). 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 (Reference Rozhanskij and Veselova1994) 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 Reference Luciani, Mora and Virmont1983). 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. (Reference Aleynikov, Breizman, Helander and Turkin2019).

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 (Reference Rozhanskij and Veselova1994) and may well have been observed experimentally in the TFTR (Tokamak Fusion Test Reactor) tokamak (Mansfield et al. Reference Mansfield, Janos, Owens, Schmidt, Bell, Cavallo, Frederickson, Ramsey and Taylor1991); 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 (Reference Rozhanskij and Veselova1994) with cold ions and in Mansfield et al. (Reference Mansfield, Janos, Owens, Schmidt, Bell, Cavallo, Frederickson, Ramsey and Taylor1991) 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. (Reference Aleynikov, Breizman, Helander and Turkin2019), Arnold et al. (Reference Arnold, Aleynikov and Helander2021) and Aleynikov et al. (Reference Aleynikov, Arnold, Breizman, Helander and Runov2023), 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. (Reference Aleynikov, Breizman, Helander and Turkin2019), Arnold et al. (Reference Arnold, Aleynikov and Helander2021) and Aleynikov et al. (Reference Aleynikov, Arnold, Breizman, Helander and Runov2023).

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