3.1 Eigen-analysis of the Poisson–Vlasov system

The perturbed Poisson equation (leaving the factor $\text{e}^{\text{i}(ky-\unicode[STIX]{x1D714}t)}$ implicit) is

(3.1) $$\begin{eqnarray}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D719}_{1}=\frac{\text{d}^{2}\hat{\unicode[STIX]{x1D719}}}{\text{d}z^{2}}-k^{2}\hat{\unicode[STIX]{x1D719}}=-\frac{q_{e}}{\unicode[STIX]{x1D716}_{0}}\int f_{\Vert 1}\,\text{d}v_{z}=-\frac{q_{e}}{\unicode[STIX]{x1D716}_{0}}\left(\hat{\unicode[STIX]{x1D719}}\frac{\text{d}n_{0}}{\text{d}\unicode[STIX]{x1D719}_{0}}+\int \tilde{f}_{\Vert }\,\text{d}v_{z}\right).\end{eqnarray}$$

The structure of the linear stability problem is then that we can regard the quantity ${\tilde{n}}=\int \tilde{f}_{\Vert }\,\text{d}v_{z}$ as a linear functional of $\hat{\unicode[STIX]{x1D719}}$ determined by the solution of the linearized Vlasov equation. The Poisson equation for given $k$ is then a generalized eigen-problem (Lewis & Symon Reference Lewis and Symon1979), of which the complex $\unicode[STIX]{x1D714}$ is effectively the eigenvalue and the spatial shape of $\hat{\unicode[STIX]{x1D719}}$ is the eigenmode. The non-uniformity of the equilibrium means that the eigenmode structure is of course not a Fourier mode. Determining all the eigenmodes and eigenvalues exactly is difficult, in part because ${\tilde{n}}$ is an integral, rather than differential, functional.

However, we are not attempting here to identify every mode; only the mechanism and growth rate of a specific type of perturbation observed to occur in simulations, namely the transverse instability. Therefore we concentrate on obtaining a good estimate of the eigenvalue, assuming we have a fair approximation to what the eigenmode shape is. For this purpose, we presuppose what characteristics of the eigenmode the simulations indicate, namely that it is present for small finite values of $k$ (long transverse wavelength) and the magnitude of its complex frequency is small. To lowest order, then, the eigenmode equation can omit the terms $k^{2}\unicode[STIX]{x1D719}$ and ${\tilde{n}}$ , in which case it must approximately satisfy

(3.2) $$\begin{eqnarray}\frac{\text{d}^{2}\hat{\unicode[STIX]{x1D719}}}{\text{d}z^{2}}=-\frac{q_{e}}{\unicode[STIX]{x1D716}_{0}}\frac{\text{d}n_{0}}{\text{d}\unicode[STIX]{x1D719}_{0}}\hat{\unicode[STIX]{x1D719}}=\hat{\unicode[STIX]{x1D719}}\frac{\text{d}}{\text{d}\unicode[STIX]{x1D719}_{0}}\left(\frac{\text{d}^{2}\unicode[STIX]{x1D719}_{0}}{\text{d}z^{2}}\right).\end{eqnarray}$$

It is easy to verify that this linear equation is satisfied by $\hat{\unicode[STIX]{x1D719}}=-\unicode[STIX]{x0394}(\text{d}\unicode[STIX]{x1D719}_{0}/\text{d}z)$ , a shift mode, which is sometimes called the Goldston mode in soliton and other literature (see e.g. Skryabin (Reference Skryabin2002)). It consists of simply a displacement of the entire hole structure in the $z$ -direction, which is obviously a neighbouring equilibrium, requiring no time derivative terms. In so far as there is a unique solution (modulo a position shift) for a BGK mode (electron hole) with specified distribution function, the shift appears to be the unique solution of this adiabatic equation. In any event, the simulations show that the mode does in fact consist predominantly of a shift, with some other minor distortions.

We then take advantage of the fact that an eigen-problem ${\mathcal{L}}\hat{\unicode[STIX]{x1D719}}=\unicode[STIX]{x1D706}{\mathcal{M}}\hat{\unicode[STIX]{x1D719}}$ where ${\mathcal{L}}$ and ${\mathcal{M}}$ are self-adjoint operatorsFootnote ¹ , is equivalent to a variational problem that finds the extremum of the quotient $Q=\langle \hat{\unicode[STIX]{x1D719}}{\mathcal{L}}\hat{\unicode[STIX]{x1D719}}\rangle /\langle \hat{\unicode[STIX]{x1D719}}{\mathcal{M}}\hat{\unicode[STIX]{x1D719}}\rangle$ , and that the eigenvalue therefore deviates from the quotient only by terms second order in any deviation of the $\hat{\unicode[STIX]{x1D719}}$ used in evaluating it from the exact eigenmode. Applying this approach we substitute the shift mode into the full Poisson equation, multiply by $\hat{\unicode[STIX]{x1D719}}$ and integrate over space. This process annihilates the adiabatic terms leaving

(3.3) $$\begin{eqnarray}\unicode[STIX]{x0394}k^{2}\int \unicode[STIX]{x1D716}_{0}\left(\frac{\text{d}\unicode[STIX]{x1D719}_{0}}{\text{d}z}\right)^{2}\,\text{d}z=-\int \frac{\text{d}\unicode[STIX]{x1D719}_{0}}{\text{d}z}q_{e}\int \tilde{f}_{\Vert }\,\text{d}v_{z}\,\text{d}z.\end{eqnarray}$$

The value of $\unicode[STIX]{x1D714}$ that satisfies this equation is then anticipated to be a good approximation to the eigen-frequency of the mode.

3.2 Force balance

A simple physical interpretation of (3.3) is that it is the conservation of momentum: balancing the electric field tension (left-hand side) with the electric force on the non-adiabatic electron perturbation (right-hand side). That makes it an overall approximate ‘kinematic’ constraint. Indeed, this interpretation is one heuristic explanation of the transverse instability mechanism. It is a balance between the effects of ‘jetting’ – the momentum transfer to particles because of the acceleration of the hole (Hutchinson & Zhou Reference Hutchinson and Zhou2016) – (right-hand side), and the stabilizing effects of electric field tension (left-hand side). We shall see shortly that in most circumstances the jetting is the predominant term and the stability boundary is approximately where the jetting is reduced to zero.

We can obtain this equation directly through force considerations without any immediate assumptions about the perturbation eigenmode shape, as follows. Force balance requires the total first-order ( $z$ -direction) force exerted by the electric field to be zero, because the electric field of a solitary potential structure is incapable of sustaining any net momentum transfer. Denoting the charge density by $\unicode[STIX]{x1D70C}$ , this is $E_{0}\unicode[STIX]{x1D70C}_{1}$ plus $E_{1}\unicode[STIX]{x1D70C}_{0}$ spatially integrated over the hole. These two forces (summed over species) are equal and opposite for a straight ( $k=0$ ) isolated hole because,

(3.4) $$\begin{eqnarray}\displaystyle \int \frac{\text{d}\unicode[STIX]{x1D719}_{0}}{\text{d}z}\frac{\unicode[STIX]{x1D70C}_{1}}{\unicode[STIX]{x1D716}_{0}}\,\text{d}z=-\int \frac{\text{d}\unicode[STIX]{x1D719}_{0}}{\text{d}z}\frac{\text{d}^{2}\unicode[STIX]{x1D719}_{1}}{\text{d}z^{2}}\,\text{d}z=-\left[\frac{\text{d}\unicode[STIX]{x1D719}_{0}}{\text{d}z}\frac{\text{d}\unicode[STIX]{x1D719}_{1}}{\text{d}z}\right]+\int \frac{\text{d}^{2}\unicode[STIX]{x1D719}_{0}}{\text{d}z^{2}}\frac{\text{d}\unicode[STIX]{x1D719}_{1}}{\text{d}z}\,\text{d}z=-\int \frac{\unicode[STIX]{x1D70C}_{0}}{\unicode[STIX]{x1D716}_{0}}\frac{\text{d}\unicode[STIX]{x1D719}_{1}}{\text{d}z}\,\text{d}z. & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

Integrating by parts the final form of (3.4) noting that $\unicode[STIX]{x1D70C}_{0}$ is a function of $\unicode[STIX]{x1D719}_{0}$ , we find

(3.5) $$\begin{eqnarray}\int \unicode[STIX]{x1D70C}_{0}\frac{\text{d}\unicode[STIX]{x1D719}_{1}}{\text{d}z}\,\text{d}z=-\int \frac{\text{d}\unicode[STIX]{x1D70C}_{0}}{\text{d}\unicode[STIX]{x1D719}_{0}}\frac{\text{d}\unicode[STIX]{x1D719}_{0}}{\text{d}z}\unicode[STIX]{x1D719}_{1}\,\text{d}z.\end{eqnarray}$$

A hole with $y$ -variation ( $k\neq 0$ ) has an additional term in $\unicode[STIX]{x1D70C}_{1}$ , arising from $\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D719}_{1}/\unicode[STIX]{x2202}y^{2}=-k^{2}\unicode[STIX]{x1D719}_{1}$ . Including it we get

(3.6) $$\begin{eqnarray}\tilde{F}\equiv -\int \frac{\text{d}\unicode[STIX]{x1D719}_{0}}{\text{d}z}\left[\unicode[STIX]{x1D70C}_{1}-\frac{\text{d}\unicode[STIX]{x1D70C}_{0}}{\text{d}\unicode[STIX]{x1D719}_{0}}\unicode[STIX]{x1D719}_{1}\right]\,\text{d}z=-\unicode[STIX]{x1D716}_{0}k^{2}\int \frac{\text{d}\unicode[STIX]{x1D719}_{0}}{\text{d}z}\unicode[STIX]{x1D719}_{1}\,\text{d}z\equiv F_{E}.\end{eqnarray}$$

The term in the square brackets is manifestly the non-adiabatic charge density $q_{e}{\tilde{n}}=q_{e}\int \tilde{f}\,\text{d}v$ . Notice that because $\unicode[STIX]{x1D719}_{0}$ is symmetric, this equation selects the antisymmetric part of $\unicode[STIX]{x1D719}_{1}$ . For a shift mode it is equivalent to the form $\int \unicode[STIX]{x1D719}_{1}^{\ast }\tilde{f}\,\text{d}^{3}v\,\text{d}^{3}x$ (whose variational Euler equation is the eigen-problem). The right-hand side becomes $F_{E}=\unicode[STIX]{x0394}k^{2}\int \unicode[STIX]{x1D716}_{0}(\text{d}\unicode[STIX]{x1D719}_{0}/\text{d}z)^{2}\,\text{d}z$ , which can be thought of as arising from the transfer of $z$ -momentum in the $y$ -direction via the electric field off-diagonal components of the Maxwell stress tensor ( $(\unicode[STIX]{x2202}/\unicode[STIX]{x2202}y)E_{y}E_{z}$ ). It causes the hole to resist kinking as if it had a tension in the $y$ -direction.

4.1 Shift mode perturbed distribution

Now specialize to the case when the eigenmode is the shift mode $\hat{\unicode[STIX]{x1D719}}(z(\unicode[STIX]{x1D70F}))=-\unicode[STIX]{x0394}(\text{d}\unicode[STIX]{x1D719}_{0}/\text{d}z)$ . Using

(4.3) $$\begin{eqnarray}\frac{\text{d}\unicode[STIX]{x1D719}_{0}}{\text{d}z}=-\frac{m_{e}}{q_{e}}v_{z}\frac{\text{d}v_{z}}{\text{d}z}=-\frac{m_{e}}{q_{e}}\frac{\text{d}v_{z}}{\text{d}t},\end{eqnarray}$$

we have

(4.4) $$\begin{eqnarray}\unicode[STIX]{x1D6F7}=\unicode[STIX]{x1D6E5}\frac{m_{e}}{q_{e}}\int _{-\infty }^{t}\frac{\text{d}v_{z}}{\text{d}\unicode[STIX]{x1D70F}}\text{e}^{-\text{i}\unicode[STIX]{x1D714}^{\prime }(\unicode[STIX]{x1D70F}-t)}\,\text{d}\unicode[STIX]{x1D70F}=\unicode[STIX]{x1D6E5}\frac{m_{e}}{q_{e}}\left\{[v_{z}\text{e}^{-\text{i}\unicode[STIX]{x1D714}^{\prime }(\unicode[STIX]{x1D70F}-t)}]_{-\infty }^{\unicode[STIX]{x1D70F}=t}+\int _{-\infty }^{t}v_{z}\text{i}\unicode[STIX]{x1D714}^{\prime }\text{e}^{-\text{i}\unicode[STIX]{x1D714}^{\prime }(\unicode[STIX]{x1D70F}-t)}\,\text{d}\unicode[STIX]{x1D70F}\right\}.\end{eqnarray}$$

For a passing particle it is better to write $v_{\infty }$ for the velocity at $\unicode[STIX]{x1D70F}\rightarrow \pm \infty$ on the orbit under consideration and then $\text{d}v_{z}/\text{d}\unicode[STIX]{x1D70F}=(\text{d}/\text{d}\unicode[STIX]{x1D70F})(v_{z}-v_{\infty })$ so that

(4.5) $$\begin{eqnarray}\unicode[STIX]{x1D6F7}=\unicode[STIX]{x1D6E5}\frac{m_{e}}{q_{e}}\left\{(v_{z}(t)-v_{\infty })+\int _{-\infty }^{t}(v_{z}-v_{\infty })\text{i}\unicode[STIX]{x1D714}^{\prime }\text{e}^{-\text{i}\unicode[STIX]{x1D714}^{\prime }(\unicode[STIX]{x1D70F}-t)}\,\text{d}\unicode[STIX]{x1D70F}\right\}.\end{eqnarray}$$

The advantage of this form is that $(v_{z}-v_{\infty })$ is zero outside the hole so the integral has a finite domain. For trapped particles this same form applies provided that $\unicode[STIX]{x1D714}^{\prime }$ has a positive imaginary part that ensures convergence and eliminates the lower limit contribution. Taking $v_{\infty }=0$ is more convenient for trapped particles.

In any case, the first (integrated) term $(v_{z}-v_{\infty })$ expresses the effect of an essentially rigid velocity shift of the trapped distribution function, because it gives rise to a term in $\tilde{f}$

(4.6) $$\begin{eqnarray}m_{e}\text{i}\unicode[STIX]{x1D714}^{\prime }\unicode[STIX]{x0394}v_{z}\frac{\unicode[STIX]{x2202}f_{0}}{\unicode[STIX]{x2202}W_{\Vert }}=-m_{e}v_{\text{hole}}v_{z}\frac{\unicode[STIX]{x2202}f_{0}}{\unicode[STIX]{x2202}W_{\Vert }}=-v_{\text{hole}}\frac{\unicode[STIX]{x2202}f_{0}}{\unicode[STIX]{x2202}v_{z}},\end{eqnarray}$$

where $v_{\text{hole}}=-\text{i}\unicode[STIX]{x1D714}^{\prime }\unicode[STIX]{x1D6E5}=\dot{\unicode[STIX]{x1D6E5}}$ is the hole incremental parallel velocity (observed in the $v_{y}$ frame). The term involving $v_{y}(\unicode[STIX]{x2202}f_{0}/\unicode[STIX]{x2202}W_{\bot })=\unicode[STIX]{x2202}f_{0}/m_{e}\unicode[STIX]{x2202}v_{y}$ is zero when integrated $\text{d}v_{y}$ .

A natural notation for the other term, because it has dimensions of velocity, is to denote it $\unicode[STIX]{x1D714}^{\prime }\tilde{L}$ with

(4.7) $$\begin{eqnarray}\tilde{L}(\unicode[STIX]{x1D714}^{\prime })\equiv \int _{-\infty }^{t}(v_{z}-v_{\infty })\text{i}\text{e}^{-\text{i}\unicode[STIX]{x1D714}^{\prime }(\unicode[STIX]{x1D70F}-t)}\,\text{d}\unicode[STIX]{x1D70F}.\end{eqnarray}$$

Although $\unicode[STIX]{x1D714}^{\prime }\tilde{L}$ is locally of higher order in the transit time, it cannot be simply ignored at positions outside the hole, because all the other terms in $\unicode[STIX]{x1D6F7}$ there are zero. It in fact gives rise to the net perturbation $\tilde{f}$ that crosses the outer boundary of the hole region. In other words, it is the source of the jetting perturbation. Besides which, when we are dealing with situations where $\unicode[STIX]{x1D714}^{\prime }$ times the transit time is not small, the $\tilde{L}$ term contains all the resulting new behaviour, and it is the thing that is hard to calculate.

4.2 Shift mode rigid velocity shift contributes no force

It is helpful to eliminate $\unicode[STIX]{x1D719}_{0}$ from the shift mode expressions by using the orbit equilibrium relation equation (4.3), allowing us to express spatial variations in terms of $v_{z}$ . It is also helpful for passing particles to transform distribution $v_{z}$ -integrals into $v_{\infty }$ -integrals. Then at constant $z$ (and $\unicode[STIX]{x1D719}$ )

(4.8) $$\begin{eqnarray}v_{\infty }\,\text{d}v_{\infty }=v_{z}\,\text{d}v_{z};\quad \text{so }\frac{\unicode[STIX]{x2202}v_{z}}{\unicode[STIX]{x2202}v_{\infty }}=\frac{v_{\infty }}{v_{z}}.\end{eqnarray}$$

Moreover, Vlasov’s equation means $f_{0}(v_{z})=f_{\infty }(v_{\infty })$ .

The non-adiabatic density perturbation ${\tilde{n}}=\int \tilde{f}\,\text{d}^{3}v$ acquires a contribution ${\tilde{n}}_{v}$ from the first term $(v_{z}-v_{\infty })$ of (4.5). However, its contribution to the total force (integrated over the hole) is zero because

(4.9) $$\begin{eqnarray}{\tilde{n}}_{v}\equiv \int \text{i}q_{e}\unicode[STIX]{x1D714}^{\prime }\frac{\unicode[STIX]{x2202}f_{0}}{\unicode[STIX]{x2202}W_{\Vert }}\unicode[STIX]{x1D6E5}\frac{m_{e}}{q_{e}}(v_{z}-v_{\infty })\,\text{d}^{3}v=\text{i}\unicode[STIX]{x1D714}^{\prime }\unicode[STIX]{x1D6E5}\int \frac{\unicode[STIX]{x2202}f_{\infty }}{\unicode[STIX]{x2202}v_{\infty }}\left(1-\frac{v_{\infty }}{v_{z}}\right)\,\text{d}v_{\infty };\end{eqnarray}$$

and so

(4.10) $$\begin{eqnarray}\displaystyle \tilde{F}_{v}\equiv \int \!-q_{e}\frac{\text{d}\unicode[STIX]{x1D719}_{0}}{\text{d}z}{\tilde{n}}_{v}\,\text{d}z & = & \displaystyle m_{e}\text{i}\unicode[STIX]{x1D714}^{\prime }\unicode[STIX]{x1D6E5}\int \int v_{z}\frac{\text{d}v_{z}}{\text{d}z}\frac{\unicode[STIX]{x2202}f_{\infty }}{\unicode[STIX]{x2202}v_{\infty }}\left(1-\frac{v_{\infty }}{v_{z}}\right)\text{d}v_{\infty }\,\text{d}z\nonumber\\ \displaystyle & = & \displaystyle m_{e}\text{i}\unicode[STIX]{x1D714}^{\prime }\unicode[STIX]{x1D6E5}\int \frac{\unicode[STIX]{x2202}f_{\infty }}{\unicode[STIX]{x2202}v_{\infty }}\left[\int v_{z}\left(1-\frac{v_{\infty }}{v_{z}}\right)\text{d}v_{z}\right]\text{d}v_{\infty }=0.\end{eqnarray}$$

The same annihilation will occur for any quantity that can be rendered into the form of a spatial integral $(\text{d}v_{z}/\text{d}z)\text{d}z\rightarrow \text{d}v_{z}$ (the square bracket here, which is zero) when all the $z$ -dependence of the integrand is in $v_{z}$ . This proof is equally valid for trapped particles as for passing. Therefore for momentum balance purposes we can ignore $\tilde{F}_{v}$ and ${\tilde{n}}_{v}$ .

The second term of equation (4.5), i.e.  $\unicode[STIX]{x1D714}^{\prime }\tilde{L}$ , is not of this form and gives non-zero total force. From now on, we shall consider only shift modes and use $\unicode[STIX]{x1D6F7}$ to denote just $\unicode[STIX]{x1D6E5}(m_{e}/q_{e})\unicode[STIX]{x1D714}^{\prime }\tilde{L}$ .

4.3 Low-frequency jetting

Major simplification of $\tilde{L}$ occurs for passing particles when $\unicode[STIX]{x1D714}^{\prime }$ is much smaller than the inverse transit time of the particles through the hole, because then the $\text{e}^{\text{i}\unicode[STIX]{x1D714}^{\prime }(\unicode[STIX]{x1D70F}-t)}$ term can be taken to be unity. We defer integration $\text{d}^{2}v_{\bot }$ over the perpendicular velocities and regard $\unicode[STIX]{x1D714}^{\prime }$ as fixed for now. The resulting contributionFootnote ² to the perturbed density from the $\unicode[STIX]{x2202}f_{0}/\unicode[STIX]{x2202}W_{\Vert }$ term is

(4.11) $$\begin{eqnarray}\displaystyle \frac{\text{d}{\tilde{n}}_{p}}{\text{d}^{2}v_{\bot }} & = & \displaystyle \int m_{e}\text{i}\unicode[STIX]{x1D714}^{\prime }\frac{\unicode[STIX]{x2202}f_{0}}{\unicode[STIX]{x2202}W_{\Vert }}\unicode[STIX]{x1D6E5}\text{i}\unicode[STIX]{x1D714}^{\prime }\int _{z_{s}}^{z}(v_{z^{\prime }}-v_{\infty })\frac{\text{d}z^{\prime }}{v_{z^{\prime }}}\,\text{d}v_{z}\nonumber\\ \displaystyle & = & \displaystyle -\unicode[STIX]{x1D714}^{\prime 2}\unicode[STIX]{x1D6E5}\int \frac{\unicode[STIX]{x2202}f_{\infty }}{\unicode[STIX]{x2202}v_{\infty }}\int _{z_{s}}^{z}\left(1-\frac{v_{\infty }}{v_{z^{\prime }}}\right)\text{d}z^{\prime }\,\frac{\text{d}v_{\infty }}{v_{z}}.\end{eqnarray}$$

Here the lower limit of the $z^{\prime }$ -integral $z_{s}$ is the start of the orbit, which depends on the sign of $v$ , which we will denote $\unicode[STIX]{x1D70E}_{v}$ .Footnote ³ Contributions come only from places with $v_{z}\neq v_{\infty }$ i.e. $\unicode[STIX]{x1D719}\neq 0$ , so $|z_{s}|$ must exceed the extent of the hole but need not actually be $\infty$ . We denote the other end position of an entire spatial integral as $z_{f}$ (finish of orbit).

The passing particle force is then

(4.12) $$\begin{eqnarray}\displaystyle \frac{\text{d}\tilde{F}_{p}}{\text{d}^{2}v_{\bot }} & = & \displaystyle \unicode[STIX]{x1D70E}_{v}\int _{z_{s}}^{z_{f}}-q_{e}\frac{\text{d}\unicode[STIX]{x1D719}_{0}}{\text{d}z}\frac{\text{d}{\tilde{n}}_{p}}{\text{d}^{2}v_{\bot }}\,\text{d}z\nonumber\\ \displaystyle & = & \displaystyle -m_{e}\unicode[STIX]{x1D714}^{\prime 2}\unicode[STIX]{x1D6E5}\unicode[STIX]{x1D70E}_{v}\int _{z_{s}}^{z_{f}}v_{z}\frac{\text{d}v_{z}}{\text{d}z}\int \frac{\unicode[STIX]{x2202}f_{\infty }}{\unicode[STIX]{x2202}v_{\infty }}\int _{z_{s}}^{z}\left(1-\frac{v_{\infty }}{v_{z^{\prime }}}\right)\text{d}z^{\prime }\frac{\text{d}v_{\infty }}{v_{z}}\,\text{d}z\nonumber\\ \displaystyle & = & \displaystyle -m_{e}\unicode[STIX]{x1D714}^{\prime 2}\unicode[STIX]{x1D6E5}\int \frac{\unicode[STIX]{x2202}f_{\infty }}{\unicode[STIX]{x2202}v_{\infty }}\unicode[STIX]{x1D70E}_{v}\int _{z_{s}}^{z_{f}}\frac{\text{d}v_{z}}{\text{d}z}\int _{z_{s}}^{z}\left(1-\frac{v_{\infty }}{v_{z^{\prime }}}\right)\text{d}z^{\prime }\,\text{d}z\,\text{d}v_{\infty }.\end{eqnarray}$$

Now we do the $\text{d}z$ integral by parts, making it

(4.13) $$\begin{eqnarray}v_{\infty }\int _{z_{s}}^{z_{f}}\left(1-\frac{v_{\infty }}{v_{z^{\prime }}}\right)\text{d}z^{\prime }-\int _{z_{s}}^{z_{f}}v_{z}\left(1-\frac{v_{\infty }}{v_{z}}\right)\text{d}z=-\int _{z_{s}}^{z_{f}}v_{z}\left(1-\frac{v_{\infty }}{v_{z}}\right)^{2}\,\text{d}z.\end{eqnarray}$$

Reversing again the order of integration, we obtain a $\text{d}v_{\infty }$ integral that can be done by parts

(4.14) $$\begin{eqnarray}\int -\frac{\unicode[STIX]{x2202}f_{\infty }}{\unicode[STIX]{x2202}v_{\infty }}v_{z}\left(1-\frac{v_{\infty }}{v_{z}}\right)^{2}\,\text{d}v_{\infty }=2f_{s}v_{z0}+\int f_{\infty }\frac{d}{\text{d}v_{\infty }}\left[v_{z}\left(1-\frac{v_{\infty }}{v_{z}}\right)^{2}\right]\text{d}v_{\infty },\end{eqnarray}$$

where $v_{z0}$ is the absolute value of $v_{z}$ at the separatrix, and $f_{s}$ is the distribution function on the separatrix. The leading term arises from the sign discontinuity of $v_{z}$ at $v_{\infty }=0$ . Now

(4.15) $$\begin{eqnarray}\displaystyle \frac{\text{d}}{\text{d}v_{\infty }}\left[v_{z}\left(1-\frac{v_{\infty }}{v_{z}}\right)^{2}\right] & = & \displaystyle \left[\frac{v_{\infty }}{v_{z}}\left(1-\frac{v_{\infty }}{v_{z}}\right)-2\left(1-\frac{v_{\infty }^{2}}{v_{z}^{2}}\right)\right]\left(1-\frac{v_{\infty }}{v_{z}}\right)\nonumber\\ \displaystyle & = & \displaystyle -2+3\frac{v_{\infty }}{v_{z}}-\frac{v_{\infty }^{3}}{v_{z}^{3}}.\end{eqnarray}$$

Hence

(4.16) $$\begin{eqnarray}\frac{\text{d}\tilde{F}_{p}}{\text{d}^{2}v_{\bot }}=-m_{e}\unicode[STIX]{x1D714}^{\prime 2}\unicode[STIX]{x1D6E5}\left\{\unicode[STIX]{x1D70E}_{v}\int _{z_{s}}^{z_{f}}2f_{s}v_{z0}+\int \left[-2+3\frac{v_{\infty }}{v_{z}}-\left(\frac{v_{\infty }}{v_{z}}\right)^{3}\right]f_{\infty }\,\text{d}v_{\infty }\,\text{d}z\right\},\end{eqnarray}$$

where $\unicode[STIX]{x1D70E}_{v}\int _{z_{s}}^{z_{f}}\text{d}z$ is simply $\int \text{d}z$ . This agreesFootnote ⁴ with the prior kinematic calculation (Hutchinson & Zhou Reference Hutchinson and Zhou2016), except for the presence of the term $2f_{s}v_{z0}$ .

The trapped particle contribution associated with $\unicode[STIX]{x2202}f/\unicode[STIX]{x2202}W_{\Vert }$ can be treated using the same sequence of partial integrations except that the $v_{\infty }=0$ choice means we cannot adopt $v_{\infty }$ as the velocity integration variable. The past orbit integral now extends to $\unicode[STIX]{x1D70F}=-\infty$ , because the orbit never escapes the hole, but assuming $\unicode[STIX]{x1D714}$ to have a positive imaginary part, the integral converges and the lower limit can be ignoredFootnote ⁵ . The phase-space element $(\unicode[STIX]{x2202}f/\unicode[STIX]{x2202}W_{\Vert })\text{d}W_{\Vert }=(\unicode[STIX]{x2202}f/\unicode[STIX]{x2202}v_{z})\,\text{d}v_{z}$ commutes with the $\text{d}z$ integral. We interpret $\text{d}W_{\Vert }$ as implying summing over positive and negative velocities, so we can consider the $\text{d}z$ integral to be over the entire relevant orbit range in the positive direction (avoiding the need for $\unicode[STIX]{x1D70E}_{v}$ ).

(4.17) $$\begin{eqnarray}\displaystyle \frac{\text{d}\tilde{F}_{t}}{\text{d}^{2}v_{\bot }} & = & \displaystyle m_{e}\unicode[STIX]{x1D6E5}\int v_{z}\frac{\text{d}v_{z}}{\text{d}z}\int \text{i}\unicode[STIX]{x1D714}^{\prime }\frac{\unicode[STIX]{x2202}f_{0}}{\unicode[STIX]{x2202}v_{z}}\text{i}\unicode[STIX]{x1D714}^{\prime }\int _{z_{s}}^{z}\,\text{d}z^{\prime }\frac{\text{d}v_{z}}{v_{z}}\,\text{d}z\nonumber\\ \displaystyle & = & \displaystyle -m_{e}\unicode[STIX]{x1D714}^{\prime 2}\unicode[STIX]{x1D6E5}\int \frac{\unicode[STIX]{x2202}f_{0}}{\unicode[STIX]{x2202}v_{z}}\int \frac{\text{d}v_{z}}{\text{d}z}\int _{z_{s}}^{z}\text{d}z^{\prime }\,\text{d}z\,\text{d}v_{z}.\nonumber\\ \displaystyle & = & \displaystyle m_{e}\unicode[STIX]{x1D714}^{\prime 2}\unicode[STIX]{x1D6E5}\int \frac{\unicode[STIX]{x2202}f_{0}}{\unicode[STIX]{x2202}v_{z}}\int v_{z}\,\text{d}z\,\text{d}v_{z}=m_{e}\unicode[STIX]{x1D714}^{\prime 2}\unicode[STIX]{x1D6E5}\int \int \frac{\unicode[STIX]{x2202}f_{0}}{\unicode[STIX]{x2202}v_{z}}v_{z}\,\text{d}v_{z}\,\text{d}z\nonumber\\ \displaystyle & = & \displaystyle m_{e}\unicode[STIX]{x1D714}^{\prime 2}\unicode[STIX]{x1D6E5}\int \left[2f_{s}v_{z0}-\int _{-v_{z0}}^{v_{z0}}f_{0}\,\text{d}v_{z}\right]\text{d}z.\end{eqnarray}$$

The term $2f_{s}v_{z0}$ , which arises from the limits of integration $\text{d}v_{z}$ at the separatrix, cancels the similar term in the passing particle force expression (4.16). Without that term, the trapped force can be considered to be simply the inertia of the trapped particles.

When $k=0$ , $\unicode[STIX]{x1D714}^{\prime }=\unicode[STIX]{x1D714}$ , we can integrate $\text{d}^{2}v_{\bot }$ , and no $\unicode[STIX]{x2202}f/\unicode[STIX]{x2202}W_{\bot }$ term need be considered. Then we find

(4.18) $$\begin{eqnarray}\tilde{F}=\tilde{F}_{p}+\tilde{F}_{t}=m_{e}\ddot{\unicode[STIX]{x1D6E5}}\int _{z_{s}}^{z_{f}}\left\{\int \left[-2+3\frac{v_{\infty }}{v_{z}}-\left(\frac{v_{\infty }}{v_{z}}\right)^{3}\right]f_{\infty }\,\text{d}v_{\infty }+\int _{-v_{z0}}^{v_{z0}}f\,\text{d}v_{z}\right\}\text{d}z,\end{eqnarray}$$

where $f$ here is the (unperturbed) one-dimensional ( $v_{z}$ )-distribution function. This force expression is in full agreement with the prior one-dimensional kinematic calculation (Hutchinson & Zhou Reference Hutchinson and Zhou2016).

We observe that since an electron hole has no net electric charge, for immobile ions the sum of the trapped electron charge and the integrated difference of the passing electron charge density from its external value must be zero. This allows us to deduce an alternative expression for the trapped particle number:

(4.19) $$\begin{eqnarray}\int _{z_{s}}^{z_{f}}\int _{-v_{z0}}^{v_{z0}}f\,\text{d}v_{z}\,\text{d}z=\int _{z_{s}}^{z_{f}}\int \left(1-\left|\frac{v_{\infty }}{v_{z}}\right|\right)f_{\infty }\,\text{d}v_{\infty }\,\text{d}z.\end{eqnarray}$$

4.4 Low-frequency momentum balance including $\unicode[STIX]{x2202}f/\unicode[STIX]{x2202}W_{\bot }$ terms

Let us suppose that the external distribution is isotropic Maxwellian (and the hole is stationary in the Maxwellian frame). Then the $kv_{y}$ terms in $\tilde{f}$ cancel each other for passing particles which has the effect of making one of the $\unicode[STIX]{x1D714}^{\prime }$ terms in equation (4.16) just $\unicode[STIX]{x1D714}$ . No such cancellation occurs for the trapped particles because the sign of $\text{d}f/\text{d}W_{\Vert }$ is reversed: the electron distribution is smaller at smaller $v_{z}$ . For shallow holes for trapped particles $|\text{d}f/\text{d}W_{\Vert }|\gg |\text{d}f/\text{d}W_{\bot }|$ so the $\text{d}f/\text{d}W_{\bot }$ term hardly contributes. A simple way to account for it is to suppose that the distribution is of the Schamel type (Schamel Reference Schamel1979) having a parallel Maxwellian of negative temperature in the trapped region. In that case,

(4.20) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}f}{\unicode[STIX]{x2202}W_{\Vert }}=\unicode[STIX]{x1D6FD}\frac{\unicode[STIX]{x2202}f}{\unicode[STIX]{x2202}W_{\bot }},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FD}$ is the inverse of the ratio (a negative quantity) of trapped parallel temperature to perpendicular. This ansatz is very convenient because it allows us simply to multiply $kv_{y}(\unicode[STIX]{x2202}f/\unicode[STIX]{x2202}W_{\Vert })$ by $(1+1/|\unicode[STIX]{x1D6FD}|)$ to account for the $\unicode[STIX]{x2202}f/\unicode[STIX]{x2202}W_{\bot }$ term. This extra factor accompanies all $k^{2}$ terms. For a slow-moving hole, $-\unicode[STIX]{x1D6FD}\approx 1+(15/16)\sqrt{\unicode[STIX]{x03C0}T_{e}/e\unicode[STIX]{x1D713}}$ , which is large for a shallow hole (small $\unicode[STIX]{x1D713}$ ). But the quantity $1/|\unicode[STIX]{x1D6FD}|\approx 0.346$ is still only a moderate correction for a deep hole $\unicode[STIX]{x1D713}=T_{e}/e$ .

We still need to integrate over $v_{y}$ to arrive at the total force. When we do so, first-order $kv_{y}$ terms coming from the cross-products $\unicode[STIX]{x1D714}kv_{y}$ integrate to zero for a symmetric $f(v_{y})$ distribution. The passing particle force $\tilde{F}_{p}$ then is unchanged in form except that only $\unicode[STIX]{x1D714}$ appears in it, not $\unicode[STIX]{x1D714}^{\prime }$ . The trapped force $\tilde{F}_{t}$ has an $\unicode[STIX]{x1D714}^{2}$ term that adds to the passing as before, plus a $k^{2}\langle v_{y}^{2}\rangle$ term that is otherwise of the same form as (4.17).

Consequently, using $\langle v_{y}^{2}\rangle =T_{y}/m_{e}$ , and considering $f$ to be the parallel distribution function, the full particle force can be taken as

(4.21) $$\begin{eqnarray}\displaystyle \tilde{F} & = & \displaystyle -m_{e}\unicode[STIX]{x1D6E5}\int \unicode[STIX]{x1D714}^{2}\left\{\int \left[-2+3\frac{v_{\infty }}{v_{z}}-\left(\frac{v_{\infty }}{v_{z}}\right)^{3}\right]f_{\infty }\,\text{d}v_{\infty }+\int _{-v_{z0}}^{v_{z0}}f\text{d}v_{z}\right\}\nonumber\\ \displaystyle & & \displaystyle \quad +\left(1+\frac{1}{|\unicode[STIX]{x1D6FD}|}\right)k^{2}\frac{T_{y}}{m_{e}}\left\{-2v_{s}f_{s}+\int _{-v_{z0}}^{v_{z0}}f\text{d}v_{z}\right\}\text{d}z.\end{eqnarray}$$

(4.22) $$\begin{eqnarray}\displaystyle & = & \displaystyle -m_{e}\unicode[STIX]{x1D6E5}\int \unicode[STIX]{x1D714}^{2}\left\{\int \left[-1+2\frac{v_{\infty }}{v_{z}}-\left(\frac{v_{\infty }}{v_{z}}\right)^{3}\right]f_{\infty }\,\text{d}v_{\infty }\right\}\nonumber\\ \displaystyle & & \displaystyle \quad +\left(1+\frac{1}{|\unicode[STIX]{x1D6FD}|}\right)k^{2}\frac{T_{y}}{m_{e}}\left\{\int _{-v_{z0}}^{v_{z0}}(f-f_{s})\,\text{d}v_{z}\right\}\text{d}z\nonumber\\ \displaystyle & = & \displaystyle m_{e}\unicode[STIX]{x1D6E5}n_{\infty }\int \left[\unicode[STIX]{x1D714}^{2}h(\unicode[STIX]{x1D712})+\left(1+\frac{1}{|\unicode[STIX]{x1D6FD}|}\right)k^{2}\frac{T_{y}}{m_{e}}g(\unicode[STIX]{x1D712})\right]\text{d}z.\end{eqnarray}$$

The curly brace expressions have been denoted by the dimensionless functions $h$ and $g$ , which depend upon potential $\unicode[STIX]{x1D719}$ expressed in terms of $\unicode[STIX]{x1D712}^{2}=-q_{e}\unicode[STIX]{x1D719}/T_{e}=m_{e}v_{z0}^{2}/2T_{e}$ . They are both positive. In so far as the electric field stress is negligible, the dispersion relation is $\tilde{F}=0$ which immediately shows that $\unicode[STIX]{x1D714}^{2}=-(\langle g\rangle /\langle h\rangle )(1+1/|\unicode[STIX]{x1D6FD}|)k^{2}T_{y}/m_{e}$ , where $\langle g\rangle$ and $\langle h\rangle$ are the spatial averages of $g$ and $h$ . The frequency is therefore pure imaginary, one root being positive, which is a growing unstable perturbation.

For an unshifted Maxwellian $f_{\infty }$ , the velocity integrals needed for (4.22) can be carried out. They are

(4.23) $$\begin{eqnarray}\displaystyle h(\unicode[STIX]{x1D712})=-2\int _{0}^{\infty }\!\left[-1+2\frac{v_{\infty }}{v_{z}}-\left(\frac{v_{\infty }}{v_{z}}\right)^{3}\right]\frac{f_{\infty }}{n_{\infty }}\,\text{d}v_{\infty }=-\frac{2}{\sqrt{\unicode[STIX]{x03C0}}}\unicode[STIX]{x1D712}+(2\unicode[STIX]{x1D712}^{2}-1)\text{e}^{\unicode[STIX]{x1D712}^{2}}\,\text{erfc}(\unicode[STIX]{x1D712})+1,\hspace{-18.0pt} & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

and (using (4.19))

(4.24) $$\begin{eqnarray}g(\unicode[STIX]{x1D712})=2\left\{\frac{f_{s}}{n_{\infty }}v_{z0}-\int _{0}^{\infty }\left(\frac{1-v_{\infty }}{v_{z}}\right)\frac{f_{\infty }}{n_{\infty }}\,\text{d}v_{\infty }\right\}=\frac{2}{\sqrt{\unicode[STIX]{x03C0}}}\unicode[STIX]{x1D712}-[1-\text{e}^{\unicode[STIX]{x1D712}^{2}}\,\text{erfc}(\unicode[STIX]{x1D712})].\end{eqnarray}$$

For small $\unicode[STIX]{x1D712}$ the power series is $\text{erfc}(\unicode[STIX]{x1D712})=1-(2/\sqrt{\unicode[STIX]{x03C0}})[\unicode[STIX]{x1D712}-\unicode[STIX]{x1D712}^{3}/3+O(\unicode[STIX]{x1D712}^{5})]$ , so $h(\unicode[STIX]{x1D712})=\unicode[STIX]{x1D712}^{2}-(2/\sqrt{\unicode[STIX]{x03C0}})(4/3)\unicode[STIX]{x1D712}^{3}+O(\unicode[STIX]{x1D712}^{4})$ and $g(\unicode[STIX]{x1D712})=\unicode[STIX]{x1D712}^{2}-(2/\sqrt{\unicode[STIX]{x03C0}})(2/3)\unicode[STIX]{x1D712}^{3}+O(\unicode[STIX]{x1D712}^{4})$ . The ratio $g/h\rightarrow 1$ as $\unicode[STIX]{x1D712}\rightarrow 0$ ; at $\unicode[STIX]{x1D712}=1$ , $g/h=1.86$ and $\langle g\rangle /\langle h\rangle =1.63$ (for a $\text{sech}^{4}z$ shape potential).

The electric tension force when $\unicode[STIX]{x1D719}_{0}=\unicode[STIX]{x1D713}\,\text{sech}^{4}(z/4\unicode[STIX]{x1D706}_{D})$ can be evaluated as

(4.25) $$\begin{eqnarray}F_{E}=\unicode[STIX]{x1D6E5}\unicode[STIX]{x1D716}_{0}k^{2}\int \left(\frac{\text{d}\unicode[STIX]{x1D719}_{0}}{\text{d}z}\right)^{2}\,\text{d}z=\unicode[STIX]{x1D6E5}\unicode[STIX]{x1D716}_{0}k^{2}\frac{\unicode[STIX]{x1D713}^{2}}{\unicode[STIX]{x1D706}_{D}}\frac{128}{315}=\unicode[STIX]{x1D6E5}k^{2}\unicode[STIX]{x1D706}_{D}\frac{q_{e}^{2}\unicode[STIX]{x1D713}^{2}n_{\infty }}{T_{e}}\frac{128}{315}.\end{eqnarray}$$

Therefore the full dispersion relation is

(4.26) $$\begin{eqnarray}m_{e}\unicode[STIX]{x1D714}^{2}\int h(\unicode[STIX]{x1D712})\,\text{d}z=-k^{2}T_{y}\left(\left(1+\frac{1}{|\unicode[STIX]{x1D6FD}|}\right)\int g(\unicode[STIX]{x1D712})\,\text{d}z-\unicode[STIX]{x1D706}_{D}\frac{q_{e}^{2}\unicode[STIX]{x1D713}^{2}}{T_{y}T_{e}}\frac{128}{315}\right).\end{eqnarray}$$

A crucial observation is that as $\unicode[STIX]{x1D713}\rightarrow 0$ , $g,h\rightarrow \unicode[STIX]{x1D712}^{2}\sim \unicode[STIX]{x1D713}$ , whereas the electric tension scales like $\unicode[STIX]{x1D713}^{2}$ . Therefore for shallow holes ( $\unicode[STIX]{x1D713}\ll eT_{e}$ ) the $F_{E}$ term is ignorable. Even for a deep hole such as $\unicode[STIX]{x1D713}=eT_{e}$ , numerical evaluation shows that $\int g(\unicode[STIX]{x1D712})\,\text{d}z=0.55\times 4\unicode[STIX]{x1D706}_{D}$ , which makes the $g$ term 5.4 times larger than $F_{E}$ . Therefore, for all but exceptionally deep holes, $F_{E}$ can be ignored for low- $k$ modes.

In summary then, ignoring $F_{E}$ , the predicted instability at low $k$ , where the transit, or bounce, time of the electrons is short compared with $1/\unicode[STIX]{x1D714}^{\prime }$ , is that the imaginary part of the frequency is

(4.27) $$\begin{eqnarray}\unicode[STIX]{x1D714}_{i}=k\sqrt{\frac{T_{y}}{m_{e}}}\sqrt{\frac{\langle g\rangle }{\langle h\rangle }\left(1+\frac{1}{|\unicode[STIX]{x1D6FD}|}\right)}\end{eqnarray}$$

where the second square root factor is unity for shallow holes and rises only to 1.48 for $\unicode[STIX]{x1D713}\simeq T_{e}/e$ . (The real part of the frequency is zero.)

In view of the proportionality of $\unicode[STIX]{x1D714}_{i}$ and $k$ , one expects that the fastest growing mode has large $k$ . However, the calculation so far is for low- $k$ . Therefore the observed growth rate in a simulation or in nature is anticipated to be at the upper end of the $k$ -range for which the low- $k$ approximations apply. We therefore need to analyse the breakdown of the low- $k$ approximations, to find the behaviour of the fastest growing modes.

4.5 Full dispersion relation including finite transit time effects

In order to determine the behaviour at high- $k$ near the instability threshold, we must consider situations where the transit time is comparable to $1/\unicode[STIX]{x1D714}^{\prime }$ . We must therefore abandon the low- $k$ approximation and fully account for $\unicode[STIX]{x1D714}^{\prime }$ in the integral $\tilde{L}$ in (4.7). Rather than pursue further analytical approximation, numerical integration is adopted. Since this requires a quadruple integration over $\unicode[STIX]{x1D70F}$ , $v_{y}$ , $v_{z}$ and $z$ , which becomes computationally expensive if not done efficiently, it is helpful to recognize that one can actually combine the evaluations corresponding to all positions $z$ into a single orbit integral $\text{d}\unicode[STIX]{x1D70F}$ . In other words, one does not have to do a different orbit integral for every position $z$ for a certain parallel energy $v_{z}^{2}$ . That integral can be done once, accumulating values for all positions $z$ and can be scaled to provide $\tilde{L}$ , $\unicode[STIX]{x1D6F7}$ and hence $\tilde{f}$ and the force contribution.

4.5.1 Evaluation of the past orbit integral for passing particles

For passing particles of specified $v_{\infty }$ , we perform spatial integrals on a uniform $z$ -grid, starting at a negative $z$ position $z_{s}$ far enough outside the hole to have its integrand negligible. Call the past orbit time there $\unicode[STIX]{x1D70F}_{s}=0$ . For each succeeding position $z_{j}=z_{s}+j\unicode[STIX]{x1D6FF}z$ , we find its corresponding orbit time as $\unicode[STIX]{x1D70F}_{j+1}-\unicode[STIX]{x1D70F}_{j}=\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70F}=\unicode[STIX]{x1D6FF}z(1/v_{j}+1/v_{j+1})/2$ which is an appropriate trapezoid increment for $\unicode[STIX]{x1D70F}=\int \text{d}z/v$ . The $v_{j}$ depend only on the known potential $\unicode[STIX]{x1D719}_{j}$ at $z_{j}$ , and on $v_{\infty }$ . The next value of $\tilde{L}\equiv \int _{\unicode[STIX]{x1D70F}_{s}}^{t}(v-v_{\infty })\text{i}\text{e}^{-\text{i}\unicode[STIX]{x1D714}^{\prime }(\unicode[STIX]{x1D70F}-t)}\,\text{d}\unicode[STIX]{x1D70F}$ is calculated writing $v_{j+1/2}=(v_{j+1}+v_{j})/2$ and using

(4.28) $$\begin{eqnarray}\displaystyle \tilde{L}_{j+1}-\text{e}^{\text{i}\unicode[STIX]{x1D714}^{\prime }\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70F}}\tilde{L}_{j} & {\approx} & \displaystyle (v_{j+1/2}-v_{\infty })\int _{\unicode[STIX]{x1D70F}_{j}}^{\unicode[STIX]{x1D70F}_{j+1}}\text{i}\text{e}^{-\text{i}\unicode[STIX]{x1D714}^{\prime }(\unicode[STIX]{x1D70F}-\unicode[STIX]{x1D70F}_{j+1})}\,\text{d}\unicode[STIX]{x1D70F}\nonumber\\ \displaystyle & = & \displaystyle -(v_{j+1/2}-v_{\infty })\frac{1}{\unicode[STIX]{x1D714}^{\prime }}[1-\text{e}^{\text{i}\unicode[STIX]{x1D714}^{\prime }\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70F}}].\end{eqnarray}$$

This approach preserves accuracy since $v$ is slowly varying even if $\text{e}^{-\text{i}\unicode[STIX]{x1D714}^{\prime }\unicode[STIX]{x1D70F}}$ is not. Thus all values of $\tilde{L}_{j}$ and hence of $\unicode[STIX]{x1D6F7}$ (4.5) and the contribution to the hole force on the entire uniform $z$ -mesh are obtained from one cumulative integral procedure, that gives

(4.29) $$\begin{eqnarray}\frac{\text{d}F_{p}}{\text{d}^{2}v_{\bot }\text{d}v_{\infty }}=\text{i}m_{e}\unicode[STIX]{x1D6E5}\left(\unicode[STIX]{x1D714}^{\prime }\frac{\unicode[STIX]{x2202}f_{0}}{\unicode[STIX]{x2202}W_{\Vert }}+(\unicode[STIX]{x1D714}-\unicode[STIX]{x1D714}^{\prime })\frac{\unicode[STIX]{x2202}f_{0}}{\unicode[STIX]{x2202}W_{\bot }}\right)\int -q_{e}\frac{\text{d}\unicode[STIX]{x1D719}}{\text{d}z}\unicode[STIX]{x1D714}^{\prime }\tilde{L}(z,\unicode[STIX]{x1D714}^{\prime })\frac{\text{d}z}{v_{z}}v_{\infty },\end{eqnarray}$$

which can then be integrated over the (passing) velocity $v_{\infty }$ and $v_{y}$ (and trivially $v_{x}$ ).

4.5.2 Evaluation of the past orbit integral for trapped particles

For trapped particles take $v_{\infty }=0$ and then denote orbits instead by the quantity $v_{\unicode[STIX]{x1D713}}=\sqrt{2(-q_{e}\unicode[STIX]{x1D713}+W_{\Vert })/m_{e}}>0$ , which is the orbit speed at $z=0$ , $\unicode[STIX]{x1D719}=\unicode[STIX]{x1D713}$ , when the parallel energy is the negative quantity $W_{\Vert }({\geqslant}q_{e}\unicode[STIX]{x1D713})$ . To obtain high resolution near the separatrix, it is best to space orbits by equal intervals of $\sqrt{-W_{\Vert }}$ . For a given $v_{\unicode[STIX]{x1D713}}$ the orbit has a finite $z$ -extent, and turning points ( $v_{z}=0$ ) at its ends. Equally spaced $z$ -positions are used, but spanning just the $z$ -extent of the orbit (different for each energy $W_{\Vert }$ ). Near the ends of the orbit, integration interpolation is optimized by determining $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70F}$ via $\unicode[STIX]{x1D6FF}t=\unicode[STIX]{x1D6FF}v/\dot{v}$ ; but near $z=0$ using $\unicode[STIX]{x1D6FF}t=\unicode[STIX]{x1D6FF}z/v$ is better. Again $\unicode[STIX]{x1D70F}$ integrals are combined with $z$ -integrals over the hole’s spatial extent, by making the $v_{y}$ integral the outermost and converting the combined $\text{d}v_{\unicode[STIX]{x1D713}}\text{d}z$ phase-space integral of the orbit contribution over the entire hole into an equivalent $\text{d}\unicode[STIX]{x1D70F}$ integral round the closed phase-space orbit. Schematically $F_{t}=\int \ldots \text{d}\unicode[STIX]{x1D70F}^{\prime }\,\text{d}z\,\text{d}v_{z}\text{d}v_{y}=\int \ldots \text{d}\unicode[STIX]{x1D70F}^{\prime }\,\text{d}\unicode[STIX]{x1D70F}v_{\unicode[STIX]{x1D713}}\,\text{d}v_{\unicode[STIX]{x1D713}}\,\text{d}v_{y}$ .Footnote ⁶

As previously noted, one must assume that $\unicode[STIX]{x1D714}^{\prime }$ has a positive imaginary part $\unicode[STIX]{x1D714}_{i}$ which ensures the backward $\unicode[STIX]{x1D70F}$ integral giving $\tilde{L}$ converges. But it is obvious that the contribution from each preceding orbit period is the same except that they are multiplied by successive factors $\exp (\text{i}\unicode[STIX]{x1D714}^{\prime }t_{b})$ , where $t_{b}$ is the orbit (bounce) period. This factor accounts for the attenuation $\exp (\unicode[STIX]{x1D714}_{i}\unicode[STIX]{x1D70F})$ and the phase change relative to the succeeding period. Therefore it is necessary to perform the integral numerically only around a single period of the orbit; the total including all the prior periods can then be synthesized by multiplying by the infinite sumFootnote ⁷ :

(4.30) $$\begin{eqnarray}\mathop{\sum }_{\ell =0}^{\infty }[\exp (\text{i}\unicode[STIX]{x1D714}^{\prime }t_{b})]^{\ell }=1/[1-\exp (\text{i}\unicode[STIX]{x1D714}^{\prime }t_{b})].\end{eqnarray}$$

It does not matter where one chooses to start and stop the single orbit numerically evaluated. For definiteness, choose to start the integral at $\unicode[STIX]{x1D70F}=0$ saving the cumulative integral $L(t)=\int _{0}^{t}v_{z}\text{i}\text{e}^{-\text{i}\unicode[STIX]{x1D714}^{\prime }(\unicode[STIX]{x1D70F}-t)}\,\text{d}\unicode[STIX]{x1D70F}$ . Then take $\unicode[STIX]{x1D70F}=0$ to be the left end of the orbit (where $v_{z}$ changes sign from negative to positive). To construct the first prior orbit for any other position $z,v_{z}$ on it, take $t$ to be the value of $\unicode[STIX]{x1D70F}$ during the first orbit at which the numerically calculated orbit passes through that position.

Since the hole is reflectionally symmetric, it suffices to calculate only half of the orbit, where $0<t\leqslant t_{b}/2$ . Starting at the left-hand end and integrating up to the right-hand end is the same as starting at the right and integrating along the negative- $v$ part of the orbit to the left. The only difference is that the velocity is opposite in sign. The resulting reversal of $\unicode[STIX]{x1D6F7}$ is cancelled by the reflection of the $z$ position and the consequent reversal of $\text{d}\unicode[STIX]{x1D719}_{0}/\text{d}z$ . Hence positive and negative velocity parts of the integral give equal contribution to the force. The complete integral for a full prior period ending in the $0<t\leqslant t_{b}/2$ segment is obtained as the sum of four parts of the orbit: $(0\rightarrow t),(t-t_{b}\rightarrow -t_{b}/2),(-t_{b}/2\rightarrow -t_{b}/2+t),(-t_{b}/2+t\rightarrow 0)$ . Noting that $\int _{-t_{b}}^{t-t_{b}}v_{z}\text{i}\text{e}^{-\text{i}\unicode[STIX]{x1D714}^{\prime }(\unicode[STIX]{x1D70F}-t)}\,\text{d}\unicode[STIX]{x1D70F}=\text{e}^{\text{i}\unicode[STIX]{x1D714}^{\prime }t_{b}}L(t)$ , $\int _{-t_{b}/2}^{t-t_{b}/2}v_{z}\text{i}\text{e}^{-\text{i}\unicode[STIX]{x1D714}^{\prime }(\unicode[STIX]{x1D70F}-t)}\,\text{d}\unicode[STIX]{x1D70F}=-\text{e}^{\text{i}\unicode[STIX]{x1D714}^{\prime }t_{b}/2}L(t)$ and $\int _{t}^{t_{b}/2}v_{z}\text{i}\text{e}^{-\text{i}\unicode[STIX]{x1D714}^{\prime }(\unicode[STIX]{x1D70F}-t)}\,\text{d}\unicode[STIX]{x1D70F}=\text{e}^{\text{i}\unicode[STIX]{x1D714}^{\prime }(t-t_{b}/2)}L(t_{b}/2)-L(t)$ one finds

(4.31) $$\begin{eqnarray}\displaystyle \int _{t-t_{b}}^{t}v_{z}(\unicode[STIX]{x1D70F})\text{i}\text{e}^{-\text{i}\unicode[STIX]{x1D714}^{\prime }(\unicode[STIX]{x1D70F}-t)}\,\text{d}\unicode[STIX]{x1D70F} & = & \displaystyle \int _{0}^{t}+\int _{t-t_{b}}^{-t_{b}/2}+\int _{-t_{b}/2}^{-t_{b}/2+t}+\int _{-t_{b}/2+t}^{0}v_{z}(\unicode[STIX]{x1D70F})\text{i}\text{e}^{-\text{i}\unicode[STIX]{x1D714}^{\prime }(\unicode[STIX]{x1D70F}-t)}\,\text{d}\unicode[STIX]{x1D70F}\nonumber\\ \displaystyle & = & \displaystyle L(t)+\text{e}^{\text{i}\unicode[STIX]{x1D714}^{\prime }t_{b}}[\text{e}^{\text{i}\unicode[STIX]{x1D714}^{\prime }(t-t_{b}/2)}L(t_{b}/2)-L(t)]\nonumber\\ \displaystyle & & \displaystyle -\,\text{e}^{\text{i}\unicode[STIX]{x1D714}^{\prime }t_{b}/2}L(t)-\text{e}^{\text{i}\unicode[STIX]{x1D714}^{\prime }t_{b}/2}[\text{e}^{\text{i}\unicode[STIX]{x1D714}^{\prime }(t-t_{b}/2)}L(t_{b}/2)-L(t)]\nonumber\\ \displaystyle & = & \displaystyle (1-\text{e}^{\text{i}\unicode[STIX]{x1D714}^{\prime }t_{b}})L(t)+\text{e}^{\text{i}\unicode[STIX]{x1D714}^{\prime }t}(\text{e}^{\text{i}\unicode[STIX]{x1D714}^{\prime }t_{b}/2}-1)L(t_{b}/2).\end{eqnarray}$$

And the full prior time integral from $-\infty$ to give $\tilde{L}$ is obtained by multiplying by (4.30).

Then

(4.32) $$\begin{eqnarray}\frac{\text{d}F_{t}}{\text{d}^{2}v_{\bot }\text{d}v_{\unicode[STIX]{x1D713}}}=\text{i}m_{e}\unicode[STIX]{x1D6E5}\left(\unicode[STIX]{x1D714}^{\prime }\frac{\unicode[STIX]{x2202}f_{0}}{\unicode[STIX]{x2202}W_{\Vert }}+(\unicode[STIX]{x1D714}-\unicode[STIX]{x1D714}^{\prime })\frac{\unicode[STIX]{x2202}f_{0}}{\unicode[STIX]{x2202}W_{\bot }}\right)\int -q_{e}\frac{\text{d}\unicode[STIX]{x1D719}}{\text{d}z}\unicode[STIX]{x1D714}^{\prime }\tilde{L}(z,\unicode[STIX]{x1D714}^{\prime })\frac{\text{d}z}{v_{z}}v_{\unicode[STIX]{x1D713}},\end{eqnarray}$$

where the integral $\text{d}z/v_{z}=\text{d}t$ is over the half-orbit with positive $v_{\unicode[STIX]{x1D713}}$ , and an equal contribution also comes from the negative- $v_{\unicode[STIX]{x1D713}}$ half.

4.5.3 Results

Implementing these algorithms gives the forces $F_{t}$ and $F_{p}$ for specified general values of $\unicode[STIX]{x1D713}$ , $k$ and $\unicode[STIX]{x1D714}$ for some specified potential form $\unicode[STIX]{x1D719}_{0}(z)$ and distribution function $f$ . A verification of the numerics is obtained by comparing the values found at $k=0$ , $\text{Re}(\unicode[STIX]{x1D714})=0$ , and $\text{Im}(\unicode[STIX]{x1D714})$ small, with the analytic expressions leading up to equation (4.18), and employing the forms (4.23) and (4.24).

Figure 2.

Comparison of numerical and analytic calculations of trapped and passing particle forces, for small value of $\unicode[STIX]{x1D714}_{i}$ and zero $k$ and $\unicode[STIX]{x1D714}_{r}$ . The forces are real for symmetric $v_{y}$ -distributions. ( $\unicode[STIX]{x1D713}$ is measured in units of $T_{e}/e$ .)

As figure 2 indicates, we find essentially exact agreement for the passing particle force $F_{p}$ , using Maxwellian distributions and a potential $\unicode[STIX]{x1D719}=\unicode[STIX]{x1D713}\,\text{sech}^{4}(z/4)$ at all values of $\unicode[STIX]{x1D713}$ . However, the trapped particle force requires a self-consistent parallel velocity distribution, which is not readily available analytically except for shallow holes. Using the approximation $f(W_{\Vert })=f_{\infty }(0)\exp (-\unicode[STIX]{x1D6FD}W_{\Vert })$ with $-\unicode[STIX]{x1D6FD}=1+(15/16)\sqrt{\unicode[STIX]{x03C0}T_{e}/e\unicode[STIX]{x1D713}}$ , better than 2 % agreement in $F_{t}$ is obtained for $\unicode[STIX]{x1D713}<0.02$ , which verifies the integration coding.

Figure 3.

Total particle force $F_{t}+F_{p}$ , and electric force $F_{E}$ , normalized to $\unicode[STIX]{x1D714}^{2}$ versus wavenumber $k$ . The dispersion relation is $F_{t}+F_{p}=F_{E}$ .

As $k$ is increased from zero the main effect is observed to occur as a reduction of $F_{t}$ , with only small changes to $F_{p}$ . In figure 3 is shown the total particle force and the electric tension force as a function of $k$ for two (fixed) imaginary frequencies $\unicode[STIX]{x1D714}=i\unicode[STIX]{x1D714}_{i}$ . It illustrates the fact that the tension force $F_{E}$ (drawn for the corresponding $\unicode[STIX]{x1D713}$ values and varying $\propto \unicode[STIX]{x1D713}^{2}$ ) is a minor correction, probably less important in practice than the approximation arising from taking the trapped distribution to be $\propto \exp (-\unicode[STIX]{x1D6FD}v_{z}^{2})$ . For small $\unicode[STIX]{x1D714}_{i}$ (figure 3 a), this plot is essentially universal. It gives the small- $k$ dispersion root at the intersection of corresponding curves, where $kv_{t}/\unicode[STIX]{x1D714}_{i}$ is slightly below 1 in accordance with (4.27). However for larger $\unicode[STIX]{x1D714}_{i}$ and correspondingly $kv_{t}$ (figure 3 b), there is no intersection for shallow holes (small $\unicode[STIX]{x1D713}$ ) because the curves of $F_{t}+F_{p}$ have a minimum above zero. Figure 3(b) also illustrates the existence of a second solution at higher $k$ for higher $\unicode[STIX]{x1D713}$ values. We shall shortly see that this behaviour arises because of changes of resonance on trapped electrons.

Figure 4.

Dispersion relation showing scaled transverse instability growth rate versus scaled transverse wavenumber, for a wide range of hole potential $\unicode[STIX]{x1D713}$ . Stationary hole, $\unicode[STIX]{x1D719}=\unicode[STIX]{x1D713}\text{sech}^{4}(z/4)$ , Maxwellian background, immobile ions. Units of $\unicode[STIX]{x1D6FE}$ , $z$ (and $k^{-1}$ ) and $\unicode[STIX]{x1D713}$ are respectively $\unicode[STIX]{x1D714}_{p}$ , $\unicode[STIX]{x1D706}_{De}$ and $T_{e}/e$ .

The dispersion relation that solves $F_{t}+F_{p}=F_{E}$ can usefully be displayed in terms of $\unicode[STIX]{x1D6FE}\equiv \unicode[STIX]{x1D714}_{i}$ as a function of given $k$ . This can be accomplished for a stationary hole in a Maxwellian distribution by bisection real root finding because the real part of $\unicode[STIX]{x1D714}$ is zero by symmetry. Both $k$ and $\unicode[STIX]{x1D6FE}$ scale like $\sqrt{\unicode[STIX]{x1D713}}$ , so plotting the resulting $\unicode[STIX]{x1D6FE}/\sqrt{\unicode[STIX]{x1D713}}$ versus $k/\sqrt{\unicode[STIX]{x1D713}}$ , as shown in figure 4, gives an almost universal curve, at least for $\unicode[STIX]{x1D713}\lesssim 0.5$ . (Here and subsequently in presenting numerical results we use non-dimensionalized parameters where units of time are $1/\unicode[STIX]{x1D714}_{p}$ , of space $\unicode[STIX]{x1D706}_{D}$ and of energy $T_{e}/e$ ; that is equivalent to setting $|q_{e}|=m_{e}=1$ .) The iterative numerical solution breaks down at the upper end of the $k$ -range ( $k_{c}$ ) where $\unicode[STIX]{x1D714}_{i}$ returns to zero, because the integrals no longer converge. Numerically it breaks down just before $k_{c}$ , because of imprecision; that is why the curves do not extend all the way back to the $k$ axis. The forces are evaluated using the same first-order approximation for $\unicode[STIX]{x1D6FD}$ which is inaccurate at $\unicode[STIX]{x1D713}\sim 1$ . That is one reason why the curves only become truly invariant at low $\unicode[STIX]{x1D713}$ . Another is that $F_{E}$ begins to become significant for $\unicode[STIX]{x1D713}\sim 1$ . Observe that the slope at the origin also somewhat exceeds the limit value of 1 (indicated by the short solid line) except when $\unicode[STIX]{x1D713}\ll 1$ . That is consistent with (4.27). The peak growth rate occurs very close to the previously published (Hutchinson Reference Hutchinson2018) estimate $k=\sqrt{\unicode[STIX]{x1D713}}/8$ . The maximum growth rate, even at low $\unicode[STIX]{x1D713}$ somewhat exceeds the prior estimate $\unicode[STIX]{x1D6FE}/\sqrt{\unicode[STIX]{x1D713}}=1/16$ . This section has thus confirmed that the transverse instability of electron holes (in Maxwellian background plasma) is a kinematic effect arising from the jetting force balance in an accelerating hole; and it has given a detailed mathematical derivation and precise numerical evaluation of coefficients.