3.1 Frequency response of the momentum rate of change

To analyse this instability, we first consider the steady-state solution of an EH. The steady-state EH potential $\unicode[STIX]{x1D719}(x)$ is considered to extend from $x_{a}$ to $x_{b}$ in the hole frame, $x_{a}$ and $x_{b}$ are taken to be far away from the centre of the hole so that both $\unicode[STIX]{x1D719}(x)$ and its derivatives vanish at these limits: $\unicode[STIX]{x1D719}(x_{a})=\unicode[STIX]{x1D719}(x_{b})=\unicode[STIX]{x1D719}^{\prime }(x_{a})=\unicode[STIX]{x1D719}^{\prime }(x_{b})=0$ . The ions and the bulk electrons are assumed to be Maxwellian and at rest in the laboratory frame with their background density being denoted by $n_{\infty }$ . The EH moves at a velocity $U$ in the laboratory frame. The sign convention is such that $U<0$ . We first adopt a cold beam approximation for ions and the finite ion temperature effect will be treated later on in this paper. The distance is normalized to $\unicode[STIX]{x1D706}_{\text{De}}$ , $\unicode[STIX]{x1D719}$ is measured in $T_{e}/e$ , the velocity is in units of $c_{s}=\sqrt{T_{e}/m_{i}}$ and the time is normalized to $1/\unicode[STIX]{x1D714}_{pi}$ . The schematic of a steady-state EH is shown in figure 3. The steady-state velocity and density of the ions in the hole frame can be derived from conservation of energy and the continuity equation:

(3.1) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}v_{0}(x)=-\displaystyle \frac{U}{|U|}\sqrt{U^{2}-2\unicode[STIX]{x1D719}(x)}\\ n_{0}(x)=n_{\infty }\displaystyle \frac{|U|}{\sqrt{U^{2}-2\unicode[STIX]{x1D719}(x)}}.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

Figure 3. Schematic of a steady-state EH with the associated phase-space structure and the ion response. (a) EH potential, (b) electron phase-space orbits, the trapped orbits are shaded, (c) the steady-state ion velocity $v_{0}$ and density $n_{0}$ in the hole frame.

In our previous paper (Hutchinson & Zhou Reference Hutchinson and Zhou2016), we have established that the motion of an EH is governed by the momentum conservation when acceleration or growth are steady. The parallel momentum contained in the electromagnetic field can be ignored for a field-aligned solitary electrostatic structure. The momentum balance is between the two components of the plasma: ${\dot{P}}_{i}+{\dot{P}}_{e}=0$ , where ${\dot{P}}_{i/e}$ represents the total inertial frame momentum rate of change for the two species. Here we extend the analysis to oscillatory acceleration or the frequency domain. To a first approximation, we are going to assume that the EH potential $\unicode[STIX]{x1D719}(x)$ in the hole frame does not change and there is a small perturbation in the hole velocity represented by an ansatz $\dot{U}\sim \text{Re}(\exp [-\text{i}\unicode[STIX]{x1D714}t])$ . Our ansatz corresponds to an eigenmode which results in a displacement of the steady-state equilibrium. It is sometimes referred to as the Goldstone mode of a soliton solution (Purwins et al. Reference Purwins, Bödeker and Liehr2005, pp. 267–308). The frequency $\unicode[STIX]{x1D714}$ we consider is much lower than the average passing electron transit frequency. The passing electrons feel a nearly constant hole acceleration during their transit (this amounts to a short-transit-time approximation for electrons), so that we can use the previous expression of ${\dot{P}}_{e}$ for $U\ll v_{\text{th,e}}$ with a full ion response correction (Zhou & Hutchinson Reference Zhou and Hutchinson2016)

(3.2) $$\begin{eqnarray}\displaystyle {\dot{P}}_{e}=-m_{e}\dot{U}n_{\infty }\int _{x_{a}}^{x_{b}}h(\sqrt{\unicode[STIX]{x1D719}(x)})+1-\frac{|U|}{\sqrt{U^{2}-2\unicode[STIX]{x1D719}(x)}}\,\text{d}x, & & \displaystyle\end{eqnarray}$$

where $h(\unicode[STIX]{x1D712})=-(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$ . The ion response correction $1-|U|/\sqrt{U^{2}-2\unicode[STIX]{x1D719}(x)}$ accounts for the ion density accumulation inside the EH. The exact shape of the trapped electron distribution does not appear directly in our approach. The total number of trapped electrons inside the EH is deduced from global charge neutrality of the solitary structure. Our ansatz treats the trapped electron phase-space structure as a holistic object.

The ions however feel an oscillating potential when they transit the hole region. The ion momentum change can be decomposed into two different terms, a momentum outflow term ${\dot{P}}_{io}$ at the boundaries and a contained momentum term ${\dot{P}}_{ic}$ . The conservation of momentum needs to be evaluated at a fixed time. Let subscripts $s$ and $f$ refer to the starting time and the final time, $a$ and $b$ refer to the positions $x_{a}$ and $x_{b}$ in the hole frame and bar denote velocities in the inertial frame (the unbarred velocities are evaluated in the hole frame). An ion particle enters the control volume at $x_{a}$ when $t=t_{s}$ exits at $x_{b}$ when $t=t_{f}$ . At $t=t_{f}$ , we have therefore

(3.3) $$\begin{eqnarray}\displaystyle \frac{{\dot{P}}_{io}}{m_{i}} & = & \displaystyle n_{bf}v_{bf}\bar{v}_{bf}-n_{af}v_{af}\bar{v}_{af}\nonumber\\ \displaystyle & = & \displaystyle n_{bf}v_{bf}\bar{v}_{bf}-n_{bf}v_{bf}\bar{v}_{af}+n_{bf}v_{bf}\bar{v}_{af}-n_{af}v_{af}\bar{v}_{af}\nonumber\\ \displaystyle & = & \displaystyle n_{bf}v_{bf}(\bar{v}_{bf}-\bar{v}_{af})+(n_{bf}v_{bf}-n_{af}v_{af})\bar{v}_{af}.\end{eqnarray}$$

The term $\bar{v}_{bf}-\bar{v}_{af}$ represents the ‘jetting’ effect due to the hole acceleration. In the comoving frame of the EH, the equation of motion of a single ion particle admits a first integral

(3.4) $$\begin{eqnarray}\displaystyle \frac{1}{2}v^{2}+\unicode[STIX]{x1D719}(x)+\int \dot{U}v\,\text{d}t=\text{Const}. & & \displaystyle\end{eqnarray}$$

Applying this conservation law between the time $t_{s}$ and $t_{f}$ , we have

(3.5) $$\begin{eqnarray}\displaystyle v_{bf}^{2}-v_{as}^{2}+2\int _{t_{s}}^{t_{f}}\dot{U}v\,\text{d}t=0. & & \displaystyle\end{eqnarray}$$

The equilibrium velocity for $\dot{U}=0$ is $v_{0}(x)=-(U/|U|)\sqrt{U^{2}-2\unicode[STIX]{x1D719}(x)}$ . The idea is to do an expansion around the equilibrium orbit. Perturbation expansion gives $v=v_{0}+v_{1}$ and $|v_{1}|/|v_{0}|\sim |t_{ab}\dot{U}/U|\ll 1$ where $t_{ab}=t_{f}-t_{s}$ is the single ion transit time. In principle, the amplitude of the hole acceleration can be made arbitrarily small to satisfy this ordering. To the leading order in the small parameter $t_{ab}\dot{U}/U$ , we expand the difference between the ion velocity exiting and entering the hole region $v_{bf}-v_{as}$

(3.6) $$\begin{eqnarray}\displaystyle v_{bf}-v_{as} & = & \displaystyle \frac{-2}{v_{bf}+v_{as}}\int _{t_{s}}^{t_{f}}\dot{U}v\,\text{d}t\nonumber\\ \displaystyle & \simeq & \displaystyle \int _{x_{a}}^{x_{b}}\frac{\dot{U}(t(x,x_{b}))}{U}\,\text{d}x.\end{eqnarray}$$

where $t(x,x_{b})=t_{f}-\int _{x}^{x_{b}}(\text{d}u/v(u))$ is an intermediate time. To keep notations simple, we will omit its explicit form while keeping in mind that unless stated otherwise, $\dot{U}$ and $v$ are evaluated when the considered ion particle is at the position indicated by the dummy variable of the integration. Taking into account the change in the hole velocity between $t_{s}$ and $t_{f}$ , we have

(3.7) $$\begin{eqnarray}\displaystyle \bar{v}_{bf}-\bar{v}_{af} & = & \displaystyle v_{bf}-v_{as}+\int _{t_{s}}^{t_{f}}\dot{U}\,\text{d}t\nonumber\\ \displaystyle & \simeq & \displaystyle \int _{x_{a}}^{x_{b}}\left(\frac{1}{U}+\frac{1}{v_{0}}\right)\dot{U}\,\text{d}x.\end{eqnarray}$$

Thus we have obtained the ‘jetting’ effect due to the acceleration of the EH to the relevant order and the first term in (3.3) can be evaluated as

(3.8) $$\begin{eqnarray}\displaystyle n_{bf}v_{bf}(\bar{v}_{bf}-\bar{v}_{af})\simeq n_{\infty }\int _{x_{a}}^{x_{b}}\left(-1-\frac{U}{v_{0}}\right)\dot{U}\,\text{d}x. & & \displaystyle\end{eqnarray}$$

To calculate the second term in (3.3), we need to know how the ion flux changes with $\dot{U}$ . We apply the continuity of an ion fluid element from $x_{a}$ to $x_{b}$

(3.9) $$\begin{eqnarray}\displaystyle n_{as}v_{as}\unicode[STIX]{x1D6FF}t_{as}=n_{bf}v_{bf}\unicode[STIX]{x1D6FF}t_{bf}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6FF}t_{as}$ and $\unicode[STIX]{x1D6FF}t_{bf}$ are two infinitesimal time duration for the same ion fluid element to enter and exit the control volume. They are related by the derivative of transit time $t_{ab}$ with respect to the starting time $t_{s}$ : $\unicode[STIX]{x1D6FF}t_{bf}\simeq \unicode[STIX]{x1D6FF}t_{as}(1+(\text{d}t_{ab}/\text{d}t_{s}))$ . Thus, to the leading order,

(3.10) $$\begin{eqnarray}\displaystyle n_{bf}v_{bf}-n_{af}v_{af} & = & \displaystyle n_{bf}v_{bf}-n_{bf}v_{bf}\frac{\unicode[STIX]{x1D6FF}t_{bf}}{\unicode[STIX]{x1D6FF}t_{as}}\frac{v_{af}}{v_{as}}\nonumber\\ \displaystyle & = & \displaystyle n_{bf}v_{bf}\left[1-\frac{\unicode[STIX]{x1D6FF}t_{bf}}{\unicode[STIX]{x1D6FF}t_{as}}\left(1+\frac{v_{af}-v_{as}}{v_{as}}\right)\right]\nonumber\\ \displaystyle & \simeq & \displaystyle n_{bf}v_{bf}\left(-\frac{\text{d}t_{ab}}{\text{d}t_{s}}-\frac{v_{af}-v_{as}}{v_{as}}\right),\end{eqnarray}$$

where we used the constancy of inflow density $n_{as}=n_{af}=n_{\infty }$ and that $\text{d}t_{ab}/\text{d}t_{s}$ is of the same order as $(v_{af}-v_{as})/v_{as}$ . The derivative $\text{d}t_{ab}/\text{d}t_{s}$ describes the non-constancy of the ion transit time due to the hole acceleration. It can be evaluated to the first order using $t_{ab}=\int _{x_{a}}^{x_{b}}(\text{d}x/v)$ and the conservation law of (3.4)

(3.11) $$\begin{eqnarray}\displaystyle \frac{\text{d}t_{ab}}{\text{d}t_{s}} & = & \displaystyle \frac{\text{d}}{\text{d}t_{s}}\int _{x_{a}}^{x_{b}}\frac{1}{v}\,\text{d}x\nonumber\\ \displaystyle & = & \displaystyle \left.\int _{x_{a}}^{x_{b}}-\frac{1}{v^{3}}\frac{\unicode[STIX]{x2202}(v^{2}/2)}{\unicode[STIX]{x2202}t_{s}}\right|_{x}\,\text{d}x\nonumber\\ \displaystyle & = & \displaystyle \left.\int _{x_{a}}^{x_{b}}-\frac{1}{v^{3}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t_{s}}\right|_{x}\left(v_{as}^{2}/2-\unicode[STIX]{x1D719}(x)-\int _{x_{a}}^{x}\dot{U}\,\text{d}x_{1}\right)\,\text{d}x\nonumber\\ \displaystyle & \simeq & \displaystyle \int _{x_{a}}^{x_{b}}\frac{1}{v^{3}}\left[-U\dot{U}(t_{s})+\int _{x_{a}}^{x}\ddot{U} \,\text{d}x_{1}\right]\,\text{d}x.\end{eqnarray}$$

$\ddot{U}$ is the rate of change of hole acceleration or the jerk evaluated when the ion particle is at $x_{1}$ . We used interchangeably $\unicode[STIX]{x2202}/\unicode[STIX]{x2202}t_{s}$ and $\unicode[STIX]{x2202}/\unicode[STIX]{x2202}t$ as $\text{d}t=\text{d}t_{s}(1+O(t_{ab}\dot{U}/U))$ for $x_{1}$ fixed. We can further remove the $\dot{U}(t_{s})$ term in (3.11) performing integration by parts

(3.12) $$\begin{eqnarray}\displaystyle \dot{U}(t_{s})-\frac{1}{U}\int _{x_{a}}^{x}\ddot{U} \,\text{d}x_{1} & = & \displaystyle \dot{U}(t_{s})-\frac{1}{U}\int _{x_{a}}^{x}v\,\text{d}\dot{U}\nonumber\\ \displaystyle & \simeq & \displaystyle \dot{U}(t_{s})-\left[\frac{\dot{U}(t)v_{0}(x)}{U}-\frac{\dot{U}(t_{s})v_{0}(x_{a})}{U}\right]+\frac{1}{U}\int _{x_{a}}^{x}\dot{U}\frac{\text{d}v_{0}}{\text{d}x_{1}}\,\text{d}x_{1}\nonumber\\ \displaystyle & \simeq & \displaystyle -\frac{\dot{U}(t)v_{0}(x)}{U}-\frac{1}{U}\int _{x_{a}}^{x}\frac{\dot{U}\unicode[STIX]{x1D719}^{\prime }(x_{1})}{v_{0}(x_{1})}\,\text{d}x_{1}.\end{eqnarray}$$

We used here $\text{d}v_{0}/\text{d}x_{1}=-\unicode[STIX]{x1D719}^{\prime }/v_{0}$ , where $\unicode[STIX]{x1D719}^{\prime }$ is the spatial derivative of $\unicode[STIX]{x1D719}$ . Combining (3.11) and (3.12), we can evaluate the right-hand side of (3.10) as

(3.13) $$\begin{eqnarray}\displaystyle n_{bf}v_{bf}\left(-\frac{\text{d}t_{ab}}{\text{d}t_{s}}-\frac{v_{af}-v_{as}}{v_{as}}\right) & \simeq & \displaystyle n_{\infty }(-U)\int _{x_{a}}^{x_{b}}-\frac{1}{v_{0}^{3}}\left[-U\dot{U}(t_{s})+\int _{x_{a}}^{x}\ddot{U} \,\text{d}x_{1}\right]-\frac{1}{U}\frac{\dot{U}}{v_{0}}\,\text{d}x\nonumber\\ \displaystyle & \simeq & \displaystyle n_{\infty }\int _{x_{a}}^{x_{b}}\frac{U^{2}}{v_{0}^{3}}\left[\frac{\dot{U}v_{0}}{U}+\frac{1}{U}\int _{x_{a}}^{x}\frac{\dot{U}\unicode[STIX]{x1D719}^{\prime }}{v_{0}}\,\text{d}x_{1}\right]+\frac{\dot{U}}{v_{0}}\,\text{d}x\nonumber\\ \displaystyle & = & \displaystyle n_{\infty }\int _{x_{a}}^{x_{b}}\frac{U}{v_{0}^{2}}\dot{U}+\frac{U}{v_{0}^{3}}\left(\int _{x_{a}}^{x}\frac{\dot{U}\unicode[STIX]{x1D719}^{\prime }}{v_{0}}\,\text{d}x_{1}\right)+\frac{\dot{U}}{v_{0}}\,\text{d}x.\end{eqnarray}$$

We choose the inertial frame as the instantaneous rest frame of the EH so that $\bar{v}_{af}=v_{af}\simeq -U$ , we have a final expression for the total rate of momentum outflow by using (3.8) and (3.13). It is first order in $t_{ab}\dot{U}/U$ , thus linear in $\dot{U}$

(3.14) $$\begin{eqnarray}\displaystyle {\dot{P}}_{io}=m_{i}n_{\infty }\int _{x_{a}}^{x_{b}}\left(-1-2\frac{U}{v_{0}}-\frac{U^{2}}{v_{0}^{2}}\right)\dot{U}-\frac{U^{2}}{v_{0}^{3}}\left(\int _{x_{a}}^{x}\frac{\dot{U}\unicode[STIX]{x1D719}^{\prime }}{v_{0}}\,\text{d}x_{1}\right)\,\text{d}x. & & \displaystyle\end{eqnarray}$$

Now we proceed to calculate the rate of change of ion momentum contained inside the control volume between $x_{a}$ and $x_{b}$ in the same inertial frame at $t=t_{f}$

(3.15) $$\begin{eqnarray}\displaystyle {\dot{P}}_{ic}=m_{i}\int _{x_{a}}^{x_{b}}\frac{\unicode[STIX]{x2202}(nv)}{\unicode[STIX]{x2202}t}\,\text{d}x+m_{i}\dot{U}\int _{x_{a}}^{x_{b}}n\,\text{d}x. & & \displaystyle\end{eqnarray}$$

The derivation of ${\dot{P}}_{ic}$ is similar in spirit to what we have shown for ${\dot{P}}_{io}$ but involves heavier algebra, we leave it to the appendix A. The final result, which is of the same order as ${\dot{P}}_{io}$ in (3.14), can be expressed as

(3.16) $$\begin{eqnarray}\displaystyle {\dot{P}}_{ic}=-m_{i}n_{\infty }\int _{x_{a}}^{x_{b}}\int _{x_{a}}^{x}\frac{U}{v_{0}^{3}}\int _{x_{a}}^{x_{1}}\dot{U}(t(x_{2},x))\unicode[STIX]{x1D719}^{\prime \prime }(x_{2})\,\text{d}x_{2}\,\text{d}x_{1}\,\text{d}x. & & \displaystyle\end{eqnarray}$$

We combine (3.14) and (3.16) to give a full expression for ${\dot{P}}_{i}={\dot{P}}_{io}+{\dot{P}}_{ic}$ . The conservation of total momentum gives an eigenmode equation for $\unicode[STIX]{x1D714}$ :

(3.17) $$\begin{eqnarray}\displaystyle {\dot{P}}_{i}(\unicode[STIX]{x1D714})+{\dot{P}}_{e}(\unicode[STIX]{x1D714})=0. & & \displaystyle\end{eqnarray}$$

The imaginary part of $\unicode[STIX]{x1D714}$ determines the stability of the corresponding eigenmode. We apply Nyquist stability analysis (Nyquist Reference Nyquist1932) to determine the stability.

The equation can be rewritten as ${\dot{P}}_{i}/{\dot{P}}_{e}+1=0$ , where ${\dot{P}}_{i}/{\dot{P}}_{e}$ is given by a long integral expression

(3.18) $$\begin{eqnarray}\displaystyle \frac{{\dot{P}}_{i}}{{\dot{P}}_{e}}(\unicode[STIX]{x1D714},U,\unicode[STIX]{x1D719}) & = & \displaystyle -\frac{m_{i}}{m_{e}}\left[\int _{x_{a}}^{x_{b}}\left(-1-2\frac{U}{v_{0}(x)}-\frac{U^{2}}{v_{0}^{2}(x)}\right)\exp \left(\text{i}\unicode[STIX]{x1D714}\int _{x}^{x_{b}}\frac{\text{d}x_{3}}{v_{0}(x_{3})}\right)\,\text{d}x\right.\nonumber\\ \displaystyle & & \displaystyle -\,\int _{x_{a}}^{x_{b}}\frac{U^{2}}{v_{0}^{3}(x)}\int _{x_{a}}^{x}\frac{\unicode[STIX]{x1D719}^{\prime }(x_{1})}{v_{0}(x_{1})}\exp \left(\text{i}\unicode[STIX]{x1D714}\int _{x_{1}}^{x_{b}}\frac{\text{d}x_{3}}{v_{0}(x_{3})}\right)\,\text{d}x_{1}\,\text{d}x\nonumber\\ \displaystyle & & \displaystyle \left.-\,\int _{x_{a}}^{x_{b}}\int _{x_{a}}^{x}\frac{U}{v_{0}^{3}(x_{1})}\int _{x_{a}}^{x_{1}}\exp \left(\text{i}\unicode[STIX]{x1D714}\int _{x_{2}}^{x}\frac{\text{d}x_{3}}{v_{0}(x_{3})}\right)\unicode[STIX]{x1D719}^{\prime \prime }(x_{2})\,\text{d}x_{2}\,\text{d}x_{1}\,\text{d}x\right]\nonumber\\ \displaystyle & & \displaystyle \bigg/\left(\int _{x_{a}}^{x_{b}}h(\sqrt{\unicode[STIX]{x1D719}(x)})+1-\frac{|U|}{\sqrt{U^{2}-2\unicode[STIX]{x1D719}(x)}}\,\text{d}x\right).\end{eqnarray}$$

${\dot{P}}_{i}/{\dot{P}}_{e}$ can be expanded to give a much simpler form in the limit where the ion kinetic energy in the hole frame is much greater than their electrostatic potential energy $U^{2}\gg 2\unicode[STIX]{x1D713}$ with $\unicode[STIX]{x1D713}$ being the maximum of $\unicode[STIX]{x1D719}$ . This approximation is very well satisfied at the onset of instability observed in our simulation. The leading term of the expanded form is

(3.19) $$\begin{eqnarray}\displaystyle \frac{{\dot{P}}_{i}}{{\dot{P}}_{e}}(\unicode[STIX]{x1D714},U,\unicode[STIX]{x1D719})\simeq -\frac{m_{i}}{m_{e}}\frac{\unicode[STIX]{x1D713}^{2}}{U^{4}}\frac{4\text{i}\displaystyle \frac{\unicode[STIX]{x1D714}}{U}I\left(\displaystyle \frac{\unicode[STIX]{x1D714}}{U}\right)+\text{i}\displaystyle \frac{\unicode[STIX]{x1D714}^{2}}{U^{2}}I^{\prime }\left(\displaystyle \frac{\unicode[STIX]{x1D714}}{U}\right)-3I_{0}}{\displaystyle \int _{x_{a}}^{x_{b}}h(\sqrt{\unicode[STIX]{x1D719}(x)})-\displaystyle \frac{\unicode[STIX]{x1D719}(x)}{U^{2}}\,\text{d}x}, & & \displaystyle\end{eqnarray}$$

where

(3.20) $$\begin{eqnarray}\displaystyle & \displaystyle I_{0}=\int _{x_{a}}^{x_{b}}\tilde{\unicode[STIX]{x1D719}}(x)^{2}\,\text{d}x, & \displaystyle\end{eqnarray}$$

(3.21) $$\begin{eqnarray}\displaystyle & \displaystyle I\left(\frac{\unicode[STIX]{x1D714}}{U}\right)=\int _{x_{a}}^{x_{b}}\int _{x_{a}}^{y}\tilde{\unicode[STIX]{x1D719}}(x)\tilde{\unicode[STIX]{x1D719}}(y)\exp \left(\text{i}\frac{\unicode[STIX]{x1D714}(x-y)}{U}\right)\,\text{d}x\,\text{d}y, & \displaystyle\end{eqnarray}$$

with $\tilde{\unicode[STIX]{x1D719}}(x)=\unicode[STIX]{x1D719}(x)/\unicode[STIX]{x1D713}$ being the normalized potential profile function. This leading term is second order in the small expansion parameter $2\unicode[STIX]{x1D713}/U^{2}$ . The details of this expansion are given in appendix B. Both the full expression and the leading-order expansion of ${\dot{P}}_{i}/{\dot{P}}_{e}(\unicode[STIX]{x1D714})$ can be evaluated numerically for real frequencies $\unicode[STIX]{x1D714}$ using for example the widely cited Schamel type (Schamel Reference Schamel1972) of EH potential $\unicode[STIX]{x1D719}(x)=\unicode[STIX]{x1D713}\,\text{sech}^{4}(x/4)$ . In the evaluation, we use the sign convention that $x_{a}$ is $-\infty$ and $x_{b}$ is $+\infty$ in the hole frame. The resulting contours are plotted in figure 4(a). The number of encirclements of the point $-1$ in the complex plane by the ${\dot{P}}_{i}/{\dot{P}}_{e}(\unicode[STIX]{x1D714})$ contour gives the number of unstable $\unicode[STIX]{x1D714}$ solutions to the eigenmode equation (3.17). There is a critical speed $U_{c}$ for $|U|$ below which the system is unstable. The leading-order term is within a few per cent of the full expression at the onset of instability. From now on, we will work with the leading-order term instead of the full expression for the purpose of studying this instability. This approximation makes the mathematics much more tractable. Our analysis is general and can be applied to any type of equilibrium EH potential, including but not limited to Schamel type of EHs.

Figure 4. (a) ${\dot{P}}_{i}/{\dot{P}}_{e}$ evaluated on the real axis for $\unicode[STIX]{x1D719}=0.23\,\text{sech}^{4}(x/4),m_{i}/m_{e}=1836$ and three different hole speeds. ${\dot{P}}_{i}/{\dot{P}}_{e}(\unicode[STIX]{x1D714})+1=0$ has two unstable zeros when $|U|<U_{c}=4.6c_{s}$ here. (b) $F(\unicode[STIX]{x1D714}/U)$ function defined in (3.22) evaluated for $\unicode[STIX]{x1D714}$ on the real axis using $\tilde{\unicode[STIX]{x1D719}}(x)=\text{sech}^{4}(x/4)$ . $F$ contour is invariant for different hole velocity $U$ .

The ${\dot{P}}_{i}/{\dot{P}}_{e}(\unicode[STIX]{x1D714})$ contour is essential to the study of this instability. We are going to take advantage of its scaling property to solve for the critical speed $U_{c}$ . To simplify the notations, we introduce two auxiliary functions $F$ and $G$ defined as

(3.22) $$\begin{eqnarray}\displaystyle & \displaystyle F\left(\frac{\unicode[STIX]{x1D714}}{U}\right)=4\text{i}\frac{\unicode[STIX]{x1D714}}{U}I\left(\frac{\unicode[STIX]{x1D714}}{U}\right)+\text{i}\frac{\unicode[STIX]{x1D714}^{2}}{U^{2}}I^{\prime }\left(\frac{\unicode[STIX]{x1D714}}{U}\right)-3I_{0}, & \displaystyle\end{eqnarray}$$

(3.23) $$\begin{eqnarray}\displaystyle & \displaystyle G(U)=\frac{m_{e}}{m_{i}}\frac{1}{\unicode[STIX]{x1D713}^{2}}\left[U^{4}\int _{x_{a}}^{x_{b}}h(\sqrt{\unicode[STIX]{x1D719}(x)})\,\text{d}x-U^{2}\int _{x_{a}}^{x_{b}}\unicode[STIX]{x1D719}(x)\,\text{d}x\right]. & \displaystyle\end{eqnarray}$$

We have therefore

(3.24) $$\begin{eqnarray}\displaystyle \frac{{\dot{P}}_{i}}{{\dot{P}}_{e}}(\unicode[STIX]{x1D714},U)=-\frac{F\left(\displaystyle \frac{\unicode[STIX]{x1D714}}{U}\right)}{G(U)}. & & \displaystyle\end{eqnarray}$$

We look for the critical speed $U_{c}$ for a given equilibrium hole potential $\unicode[STIX]{x1D719}$ such that the ${\dot{P}}_{i}/{\dot{P}}_{e}(\unicode[STIX]{x1D714})$ contour crosses the point $-1$ in the complex plane. $F$ depends on $U$ through $\unicode[STIX]{x1D714}/U$ , it gives the same contour for different $U$ values when $\unicode[STIX]{x1D714}$ is evaluated on the real axis, although at different $\unicode[STIX]{x1D714}$ values. While $F$ is a complex-valued function, function $G$ only takes real values. $U$ scales the size of the ${\dot{P}}_{i}/{\dot{P}}_{e}(\unicode[STIX]{x1D714})$ contour through $G$ . This property is demonstrated in figure 4. The identical $F(\unicode[STIX]{x1D714}/U)$ contour for different values of $U$ using a Schamel type of EH potential is shown in figure 4(b). We denote its intersection with the positive real axis by $C(\tilde{\unicode[STIX]{x1D719}})$ . The existence of this intersection $C(\tilde{\unicode[STIX]{x1D719}})$ is guaranteed for a general class of admissible hole potential $\unicode[STIX]{x1D719}(x)$ , which we will show later in this paper. The critical speed $U_{c}$ , below which the system is unstable, satisfies an equation

(3.25) $$\begin{eqnarray}\displaystyle -\frac{C(\tilde{\unicode[STIX]{x1D719}})}{G(-U_{c})}=-1. & & \displaystyle\end{eqnarray}$$

The above equation gives a quadratic equation in $U_{c}^{2}$

(3.26) $$\begin{eqnarray}\displaystyle U_{c}^{4}\int _{x_{a}}^{x_{b}}h(\sqrt{\unicode[STIX]{x1D719}(x)})\,\text{d}x-U_{c}^{2}\int _{x_{a}}^{x_{b}}\unicode[STIX]{x1D719}(x)\,\text{d}x-\unicode[STIX]{x1D713}^{2}\frac{m_{i}}{m_{e}}C(\tilde{\unicode[STIX]{x1D719}})=0. & & \displaystyle\end{eqnarray}$$

The unique real and positive solution of $U_{c}^{2}$ is

(3.27) $$\begin{eqnarray}\displaystyle U_{c}^{2}=\frac{\displaystyle \int _{x_{a}}^{x_{b}}\unicode[STIX]{x1D719}(x)\,\text{d}x+\sqrt{\left(\displaystyle \int _{x_{a}}^{x_{b}}\unicode[STIX]{x1D719}(x)\,\text{d}x\right)^{2}+4\unicode[STIX]{x1D713}^{2}\displaystyle \frac{m_{i}}{m_{e}}C(\tilde{\unicode[STIX]{x1D719}})\displaystyle \int _{x_{a}}^{x_{b}}h(\sqrt{\unicode[STIX]{x1D719}(x)})\,\text{d}x}}{2\displaystyle \int _{x_{a}}^{x_{b}}h(\sqrt{\unicode[STIX]{x1D719}(x)})\,\text{d}x}. & & \displaystyle\end{eqnarray}$$

Now we calculate the oscillation frequency of the unstable eigenmode. At marginal instability, the ${\dot{P}}_{i}/{\dot{P}}_{e}(\unicode[STIX]{x1D714})$ contour crosses $-1$ . We need to find the frequency for which this crossing happens. The imaginary part of ${\dot{P}}_{i}/{\dot{P}}_{e}$ is

(3.28) $$\begin{eqnarray}\displaystyle \text{Im}\left(\frac{{\dot{P}}_{i}}{{\dot{P}}_{e}}(\unicode[STIX]{x1D714})\right)=-\frac{\text{Im}\left(F\left(\displaystyle \frac{\unicode[STIX]{x1D714}}{U}\right)\right)}{G(U)}, & & \displaystyle\end{eqnarray}$$

where

(3.29) $$\begin{eqnarray}\displaystyle \text{Im}\left(F\left(\frac{\unicode[STIX]{x1D714}}{U}\right)\right)=4\frac{\unicode[STIX]{x1D714}}{U}\text{Re}\left(I\left(\frac{\unicode[STIX]{x1D714}}{U}\right)\right)+\frac{\unicode[STIX]{x1D714}^{2}}{U^{2}}\text{Re}\left(I^{\prime }\left(\frac{\unicode[STIX]{x1D714}}{U}\right)\right). & & \displaystyle\end{eqnarray}$$

To calculate the imaginary part of ${\dot{P}}_{i}/{\dot{P}}_{e}$ , we need to evaluate $\text{Re}(I(\unicode[STIX]{x1D714}/U))$ . It can be shown by taking $x_{a}$ and $x_{b}$ to infinity that

(3.30) $$\begin{eqnarray}\displaystyle \text{Re}\left(I\left(\frac{\unicode[STIX]{x1D714}}{U}\right)\right) & = & \displaystyle \frac{1}{2}\left(\int _{-\infty }^{\infty }\tilde{\unicode[STIX]{x1D719}}(x)\exp \left(\text{i}\frac{\unicode[STIX]{x1D714}}{U}x\right)\,\text{d}x\right)\left(\int _{-\infty }^{\infty }\tilde{\unicode[STIX]{x1D719}}(y)\exp \left(-\text{i}\frac{\unicode[STIX]{x1D714}}{U}y\right)\,\text{d}y\right)\nonumber\\ \displaystyle & = & \displaystyle \frac{1}{2}\tilde{\unicode[STIX]{x1D6F7}}\left(\frac{\unicode[STIX]{x1D714}}{U}\right)^{2},\end{eqnarray}$$

where $\tilde{\unicode[STIX]{x1D6F7}}$ is the modulus of the Fourier transform of $\tilde{\unicode[STIX]{x1D719}}$ . Thus $\text{Im}({\dot{P}}_{i}/{\dot{P}}_{e})(\unicode[STIX]{x1D714})=0$ gives

(3.31) $$\begin{eqnarray}\displaystyle \text{Im}\left(F\left(\displaystyle \frac{\unicode[STIX]{x1D714}}{U}\right)\right)\equiv \frac{\unicode[STIX]{x1D714}}{U}\tilde{\unicode[STIX]{x1D6F7}}\left(\frac{\unicode[STIX]{x1D714}}{U}\right)\left(2\tilde{\unicode[STIX]{x1D6F7}}\left(\frac{\unicode[STIX]{x1D714}}{U}\right)+\frac{\unicode[STIX]{x1D714}}{U}\tilde{\unicode[STIX]{x1D6F7}}^{\prime }\left(\frac{\unicode[STIX]{x1D714}}{U}\right)\right)=0. & & \displaystyle\end{eqnarray}$$

Equation (3.31) admits three real solutions for $\unicode[STIX]{x1D714}$ , one is the trivial $\unicode[STIX]{x1D714}=0$ , the other two solutions are given by the equation

(3.32) $$\begin{eqnarray}\displaystyle -\frac{\tilde{\unicode[STIX]{x1D6F7}}^{\prime }\left(\displaystyle \frac{\unicode[STIX]{x1D714}}{U}\right)}{2\tilde{\unicode[STIX]{x1D6F7}}\left(\displaystyle \frac{\unicode[STIX]{x1D714}}{U}\right)}=\frac{1}{\unicode[STIX]{x1D714}/U}. & & \displaystyle\end{eqnarray}$$

Since $\tilde{\unicode[STIX]{x1D719}}$ is real, we have $\tilde{\unicode[STIX]{x1D6F7}}$ is an even function and (3.32) gives two solutions that have the opposite sign. We define the positive $\unicode[STIX]{x1D714}$ solution of (3.32) as $\unicode[STIX]{x1D714}_{0}$ and a length scale $L$

(3.33) $$\begin{eqnarray}\displaystyle L\equiv \frac{|U|}{\unicode[STIX]{x1D714}_{0}}. & & \displaystyle\end{eqnarray}$$

The ${\dot{P}}_{i}/{\dot{P}}_{e}$ contour crosses the real axis at frequency $\unicode[STIX]{x1D714}_{0}$ for a given EH potential $\unicode[STIX]{x1D719}$ and $U$ . $L$ is the characteristic length of the EH and it is entirely determined by the EH potential shape $\tilde{\unicode[STIX]{x1D719}}$ . At the onset of instability, we have ${\dot{P}}_{i}/{\dot{P}}_{e}(\unicode[STIX]{x1D714}_{0}(U_{c}),-U_{c})=-1$ , $\unicode[STIX]{x1D714}_{0}$ evaluated for the critical speed $U_{c}$ is therefore the angular frequency of the initially growing unstable eigenmode. We define this frequency as

(3.34) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714}_{c}\equiv \unicode[STIX]{x1D714}_{0}(U_{c})\equiv \frac{U_{c}}{L}. & & \displaystyle\end{eqnarray}$$

It is the critical ion transit frequency through the hole potential. The frequency of the growing oscillation corresponds to a physical frequency of the system.

The existence of this critical frequency $\unicode[STIX]{x1D714}_{c}$ and the crossing point $C(\tilde{\unicode[STIX]{x1D719}})$ are guaranteed by the continuous differentiability of the EH potential $\unicode[STIX]{x1D719}(x)$ . A physical EH potential $\unicode[STIX]{x1D719}(x)$ should possess a second derivative as it satisfies Poisson’s equation $\unicode[STIX]{x1D719}^{\prime \prime }(x)+\unicode[STIX]{x1D70C}(x)/\unicode[STIX]{x1D716}_{0}=0$ , and a physical $\unicode[STIX]{x1D70C}(x)$ should have bounded variation. This smoothness requirement constrains the asymptotic behaviour of its Fourier transform. Function $\tilde{\unicode[STIX]{x1D6F7}}(p)$ decays at least as fast as $p^{-3}$ when $p\rightarrow \infty$ (see Katznelson Reference Katznelson2002, pp. 188–200). Thus we have $-\tilde{\unicode[STIX]{x1D6F7}}^{\prime }(p)/2\tilde{\unicode[STIX]{x1D6F7}}(p)=-(1/2)\,\text{d}\ln (\tilde{\unicode[STIX]{x1D6F7}}(p))/\text{d}p\geqslant 3/2p>1/p$ as $p\rightarrow \infty$ . While as $p\rightarrow 0$ , we have $-\tilde{\unicode[STIX]{x1D6F7}}^{\prime }(p)/2\tilde{\unicode[STIX]{x1D6F7}}(p)\rightarrow -\tilde{\unicode[STIX]{x1D6F7}}^{\prime }(0)/2\tilde{\unicode[STIX]{x1D6F7}}(0)\ll 1/p$ . Solutions are guaranteed for (3.32). The behaviour of the $F$ contour is as follows. As $\unicode[STIX]{x1D714}/U\rightarrow 0^{+}$ , we have $\text{Re}(F(\unicode[STIX]{x1D714}/U))\rightarrow -3I_{0}<0$ and $\text{Im}(F(\unicode[STIX]{x1D714}/U))\rightarrow 0^{+}$ . Asymptotic analysis as $\unicode[STIX]{x1D714}/U\rightarrow +\infty$ gives $\text{Re}(F(\unicode[STIX]{x1D714}/U))\rightarrow 0^{+}$ and $\text{Im}(F(\unicode[STIX]{x1D714}/U))\rightarrow 0^{-}$ . In other words, with $\unicode[STIX]{x1D714}/U$ increasing from $0$ to infinity, the $F$ contour starts from a point on the negative real axis, goes into the upper half-plane, crosses the positive real axis at $1/L$ and returns to zero. The crossing point $C(\tilde{\unicode[STIX]{x1D719}})$ shown in figure 4(b), which is crucial to this instability, is an universal feature for all physically admissible hole potentials $\unicode[STIX]{x1D719}(x)$ and we have $C(\tilde{\unicode[STIX]{x1D719}})=F(1/L)$ . Contours without a crossing can be obtained only from unphysical hole shapes. For example, $\tilde{\unicode[STIX]{x1D719}}=\exp (-|x|/\unicode[STIX]{x1D706})$ does not give a crossing point and is therefore stable; but it is unphysical, as the electric field is undefined at $x=0$ .

Figure 5. The critical values of hole speed in the ion frame below which the instability occurs for different sized EHs and two different mass ratios. The theoretical stability boundaries ( $\unicode[STIX]{x1D6FE}=0$ ) and the $\unicode[STIX]{x1D6FE}=0.1$ growth rate boundaries for Schamel type of EHs $\unicode[STIX]{x1D719}(x)=\unicode[STIX]{x1D713}\,\text{sech}^{4}(x/4)$ are plotted as reference lines. The observational data point and the numerical calculation of the same $\unicode[STIX]{x1D713}$ correspond to the same run. The ion reflection limit is much lower than the instability threshold, hence our approximation $U^{2}\gg 2\unicode[STIX]{x1D713}$ is well satisfied. All the PIC runs have $T_{e}/T_{i}=20$ .

Figure 6. The oscillations seen in our simulation are Fourier analysed to extract the main frequency for the first few periods of unstable oscillations. The uncertainty in the theoretically predicted frequency due to the uncertainty of $U_{c}$ used in (3.34) is shown by the grey uncertainty bands. Notice that the unstable oscillation frequency is in general a few times the ion plasma frequency.

Having obtained the analytic solution for the critical speed and the unstable oscillation frequency, we now compare these results with our PIC simulation observations. The hole pushing technique enables us to explore continuously the EH velocity in the ion frame. We performed a series of runs with different initialization to create EHs of different sizes. Then we determined the critical speed $U_{c}$ by inspecting the onset of unstable velocity oscillations as the hole speed is decreased. It is compared with the threshold speeds obtained by solving the Nyquist stability problem numerically using the $\unicode[STIX]{x1D719}(x)$ right before the instability onset from the same run. The electrostatic potential output $\unicode[STIX]{x1D719}(x)$ from our PIC simulation is used to construct the $F$ contour numerically from (3.22) and find its crossing point $C(\tilde{\unicode[STIX]{x1D719}})$ . We use the formula for $U_{c}$ in (3.27) to calculate its predicted value. This method takes into account the exact potential shape of the EH in our PIC simulation which is different from one run to anotherFootnote ³ . The results are presented in figure 5. The solid lines are obtained using (3.27) assuming a Schamel type of EH potential. This solution’s asymptotic behaviour comes from the special function $h(\unicode[STIX]{x1D712})$ (Hutchinson & Zhou Reference Hutchinson and Zhou2016): $h(\unicode[STIX]{x1D712})\rightarrow \unicode[STIX]{x1D712}^{2}-(8/3\sqrt{\unicode[STIX]{x03C0}})\unicode[STIX]{x1D712}^{3}$ as $\unicode[STIX]{x1D712}\rightarrow 0$ and $h(\unicode[STIX]{x1D712})\rightarrow 1-(2/\sqrt{\unicode[STIX]{x03C0}}\unicode[STIX]{x1D712})$ as $\unicode[STIX]{x1D712}\sim 1$ . For shallow holes $\unicode[STIX]{x1D713}\ll 1$ , we have $U_{c}\simeq 1+O(\unicode[STIX]{x1D713})$ and for deep holes $\unicode[STIX]{x1D713}\sim 1$ , $U_{c}\sim (m_{i}/m_{e})^{1/4}\unicode[STIX]{x1D713}^{1/2}$ . The slight deviation of our data points from the solid curves represents the deviation of the hole potential in our PIC simulation from the Schamel type. The full calculation using the exact hole potential yields a good agreement with the observation. We have runs with two different mass ratios and our results show the $(m_{i}/m_{e})^{1/4}$ scaling of $U_{c}$ with the mass ratio as predicted by the theory. The linear growth rate $\unicode[STIX]{x1D6FE}=\text{Im}(\unicode[STIX]{x1D714})$ of the instability when $|U|<U_{c}$ can be evaluated by solving the equation ${\dot{P}}_{i}/{\dot{P}}_{e}(\unicode[STIX]{x1D714})+1=0$ numerically for given $\unicode[STIX]{x1D719}(x)$ , $U$ and the mass ratio. In figure 5, we show in dashed lines the EH velocity calculated as a function of $\unicode[STIX]{x1D713}$ for $\unicode[STIX]{x1D6FE}=0.1$ and Schamel type of EH potential as a useful reference line. We shall give a detailed analysis of the growth rate in § 3.3.

In figure 6, we show the frequencies of the unstable velocity oscillations seen in our simulation. The EH position and hence velocity $U$ is obtained from the hole-tracking module for each PIC simulation time step of $0.3/\unicode[STIX]{x1D714}_{pe}$ . A discrete Fourier transform of $U$ during the first few unstable periods gives a sharp peak centred at the oscillation frequency. The error bar is given by the frequency range above half-peak power. They are plotted against the characteristic hole width $L$ defined in (3.33) calculated using the potential $\unicode[STIX]{x1D719}(x)$ right before the oscillation onset from our PIC simulations. $\unicode[STIX]{x1D719}(x)$ gives numerically the $F$ contour and it crosses the positive real axis at $F(1/L)$ . EH initialization is adjusted in our PIC code to give EHs of different width $L$ and a narrow range of potential height $\unicode[STIX]{x1D713}\sim 0.45$ , hence $U_{c}$ . The oscillation frequency predicted by our theory is inversely proportional to the hole width $L$ : $\unicode[STIX]{x1D714}_{c}=U_{c}/L$ . Good quantitative agreement is achieved between the observations and our theory. Our analysis captures the correct scaling with the mass ratio, which highlights the importance of ion dynamics for this instability.

The instability threshold $U_{c}$ is scale invariant. If we apply a change of scale $x\rightarrow \unicode[STIX]{x1D706}x$ , the equation (3.26) giving the critical speed $U_{c}$ is invariant under this change of scale as each term is multiplied by the same factor $1/\unicode[STIX]{x1D706}$ . This property is obvious for the first two terms in (3.26). For the third term, it can be shown that

(3.35) $$\begin{eqnarray}\displaystyle F_{\unicode[STIX]{x1D706}}\left(\frac{\unicode[STIX]{x1D714}}{U}\right)=\frac{1}{\unicode[STIX]{x1D706}}F\left(\frac{\unicode[STIX]{x1D714}}{U\unicode[STIX]{x1D706}}\right). & & \displaystyle\end{eqnarray}$$

Hence the crossing point satisfies the same scaling relation $C_{\unicode[STIX]{x1D706}}(\tilde{\unicode[STIX]{x1D719}})=(1/\unicode[STIX]{x1D706})C(\tilde{\unicode[STIX]{x1D719}})$ and $U_{c}$ remains invariant. However, the oscillation frequency scales linearly with $\unicode[STIX]{x1D706}$ : $\unicode[STIX]{x1D714}_{c,\unicode[STIX]{x1D706}}=\unicode[STIX]{x1D706}\unicode[STIX]{x1D714}_{c}$ . For example, two different EH potentials $\unicode[STIX]{x1D719}(x)=\unicode[STIX]{x1D713}\,\text{sech}^{4}(x)$ and $\unicode[STIX]{x1D719}(x)=\unicode[STIX]{x1D713}\,\text{sech}^{4}(x/4)$ have the same threshold $U_{c}$ , while the unstable oscillation frequency for the first potential profile is four times as high. This argument explains why the runs in figure 6 have a similar $U_{c}$ but different oscillation frequencies.

3.2 Counter-streaming ions

We have shown an example of the instability observed in a plasma with counter-propagating ions in figure 2. If the EH potential $\unicode[STIX]{x1D719}(x)$ is symmetric, then the counter-streaming situation with an EH at rest $U=0$ and two ion streams travelling at $\pm v_{i}$ is equivalent to having one single ion stream at rest and the EH travelling at $U=v_{i}$ for the described instability mechanism. The sign convention is our analysis is such that the ions enter from $-\infty$ and exit at $+\infty$ . The change of hole velocity from $U$ to $-U$ in the ion frame results in flipping the sign convention thus $x_{a}$ and $x_{b}$ . When the potential $\unicode[STIX]{x1D719}(x)$ is symmetric so that $\unicode[STIX]{x1D719}^{\prime }(-x)=-\unicode[STIX]{x1D719}^{\prime }(x)$ and $\unicode[STIX]{x1D719}^{\prime \prime }(-x)=\unicode[STIX]{x1D719}^{\prime \prime }(x)$ , equations (3.14) and (3.16) show that the resulting ${\dot{P}}_{i}$ is exactly the opposite as $\dot{U}$ has an opposite sign under the two opposite sign conventions. Therefore, the contribution to the total ${\dot{P}}_{i}$ from the two ion streams, evaluated with the same sign convention, should be exactly equal and add up. More concretely, an ion particle arriving from the left sees the same potential as an ion particle arriving from the right. However, the same EH acceleration $\dot{U}$ works in an opposite way for them. This argument explains why the instability threshold observed in the counter-streaming ion plasma is identical to the threshold value for the single ion stream case. The two situations are equivalent in terms of linear stability. Once the instability has fully grown, the nonlinear stage of the instability can be different for the two cases.

3.3 Linear growth rate

The linear growth rate of the instability is obtained by solving the eigenmode equation ${\dot{P}}_{i}/{\dot{P}}_{e}(\unicode[STIX]{x1D714},U)+1=0$ . The growth rate $\unicode[STIX]{x1D6FE}$ is the imaginary part of the solution $\unicode[STIX]{x1D714}$ : $\unicode[STIX]{x1D6FE}=\text{Im}(\unicode[STIX]{x1D714})$ . Although an analytic solution for arbitrary $U$ is too difficult to obtain, we can obtain $\unicode[STIX]{x1D6FE}$ by an expansion near marginal instability. Recall that at marginal instability, we have

(3.36) $$\begin{eqnarray}\displaystyle \frac{{\dot{P}}_{i}}{{\dot{P}}_{e}}\left(\unicode[STIX]{x1D714}_{c}=\frac{U_{c}}{L},-U_{c}\right)=-1. & & \displaystyle\end{eqnarray}$$

If the hole velocity is $U=-U_{c}+\unicode[STIX]{x0394}U$ such that $|\unicode[STIX]{x0394}U/U_{c}|\ll 1$ . We need to find $\unicode[STIX]{x1D714}=\unicode[STIX]{x1D714}_{c}+\unicode[STIX]{x0394}\unicode[STIX]{x1D714}$ with $|\unicode[STIX]{x0394}\unicode[STIX]{x1D714}/\unicode[STIX]{x1D714}_{c}|\ll 1$ that satisfies the eigenmode equation ${\dot{P}}_{i}/{\dot{P}}_{e}(\unicode[STIX]{x1D714},U)=-1$ . A linear expansion gives

(3.37) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x0394}\unicode[STIX]{x1D714}\left.\frac{\unicode[STIX]{x2202}({\dot{P}}_{i}/{\dot{P}}_{e})}{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}}\right|_{\unicode[STIX]{x1D714}_{c},-U_{c}}+\unicode[STIX]{x0394}U\left.\frac{\unicode[STIX]{x2202}({\dot{P}}_{i}/{\dot{P}}_{e})}{\unicode[STIX]{x2202}U}\right|_{\unicode[STIX]{x1D714}_{c},-U_{c}}=0. & & \displaystyle\end{eqnarray}$$

Substituting ${\dot{P}}_{i}/{\dot{P}}_{e}(\unicode[STIX]{x1D714},U)=-F(\unicode[STIX]{x1D714}/U)/G(U)$ , the two partial derivatives in (3.37) can be evaluated with functions $F$ and $G$

(3.38) $$\begin{eqnarray}\displaystyle & \displaystyle \left.\frac{\unicode[STIX]{x2202}({\dot{P}}_{i}/{\dot{P}}_{e})}{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}}\right|_{\unicode[STIX]{x1D714}_{c},-U_{c}}\,=\,\frac{F^{\prime }(-\unicode[STIX]{x1D714}_{c}/U_{c})}{U_{c}G(-U_{c})}, & \displaystyle\end{eqnarray}$$

(3.39) $$\begin{eqnarray}\displaystyle & \displaystyle \left.\frac{\unicode[STIX]{x2202}({\dot{P}}_{i}/{\dot{P}}_{e})}{\unicode[STIX]{x2202}U_{c}}\right|_{\unicode[STIX]{x1D714}_{c},-U_{c}}\,=\,\frac{F^{\prime }(-\unicode[STIX]{x1D714}_{c}/U_{c})G(-U_{c})(\unicode[STIX]{x1D714}_{c}/U_{c}^{2})+F(-\unicode[STIX]{x1D714}_{c}/U_{c})G^{\prime }(-U_{c})}{G(-U_{c})^{2}}. & \displaystyle\end{eqnarray}$$

Hence

(3.40) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x0394}\unicode[STIX]{x1D714} & = & \displaystyle \unicode[STIX]{x0394}U\left\{-\frac{\unicode[STIX]{x1D714}_{c}}{U_{c}}-U_{c}\frac{F(-\unicode[STIX]{x1D714}_{c}/U_{c})}{G(-U_{c})}\frac{G^{\prime }(-U_{c})}{F^{\prime }(-\unicode[STIX]{x1D714}_{c}/U_{c})}\right\}\nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x0394}U\left\{-\frac{\unicode[STIX]{x1D714}_{c}}{U_{c}}-U_{c}\frac{G^{\prime }(-U_{c})}{F^{\prime }(-\unicode[STIX]{x1D714}_{c}/U_{c})}\right\},\end{eqnarray}$$

where we used $F(-\unicode[STIX]{x1D714}_{c}/U_{c})/G(-U_{c})=1$ . $G$ is an even polynomial function defined in (3.23) and $-U_{c}G^{\prime }(-U_{c})\equiv U_{c}G^{\prime }(U_{c})$ can be evaluated as

(3.41) $$\begin{eqnarray}\displaystyle U_{c}G^{\prime }(U_{c})=4F\left(\frac{\unicode[STIX]{x1D714}_{c}}{U_{c}}\right)+2U_{c}^{2}\frac{m_{e}}{m_{i}}\frac{1}{\unicode[STIX]{x1D713}^{2}}\int _{x_{a}}^{x_{b}}\unicode[STIX]{x1D719}(x)\,\text{d}x. & & \displaystyle\end{eqnarray}$$

Hence the final expression for $\unicode[STIX]{x0394}\unicode[STIX]{x1D714}$ is

(3.42) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x0394}\unicode[STIX]{x1D714}=\unicode[STIX]{x0394}U\left\{-\frac{\unicode[STIX]{x1D714}_{c}}{U_{c}}+\frac{4F(\unicode[STIX]{x1D714}_{c}/U_{c})}{F^{\prime }(-\unicode[STIX]{x1D714}_{c}/U_{c})}+2U_{c}^{2}\frac{m_{e}}{m_{i}}\frac{1}{F^{\prime }(-\unicode[STIX]{x1D714}_{c}/U_{c})}\displaystyle \frac{1}{\unicode[STIX]{x1D713}^{2}}\displaystyle \int _{x_{a}}^{x_{b}}\unicode[STIX]{x1D719}(x)\,\text{d}x\right\}. & & \displaystyle\end{eqnarray}$$

The growth rate $\unicode[STIX]{x1D6FE}$ is the imaginary part of $\unicode[STIX]{x0394}\unicode[STIX]{x1D714}$ . The real part of $\unicode[STIX]{x0394}\unicode[STIX]{x1D714}$ gives a small correction to the oscillation frequency $\unicode[STIX]{x1D714}_{c}$ when $U$ is different from $U_{c}$ . We also have $\unicode[STIX]{x1D714}_{c}/U_{c}=1/L$ , a constant only depending on the hole shape, and $F(\unicode[STIX]{x1D714}_{c}/U_{c})=F(1/L)=C(\tilde{\unicode[STIX]{x1D719}})$ . While $F(\unicode[STIX]{x1D714}_{c}/U_{c})$ is a real number, the derivative $F^{\prime }(-\unicode[STIX]{x1D714}_{c}/U_{c})$ is complex and $\unicode[STIX]{x1D6FE}$ is given by

(3.43) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FE}=\unicode[STIX]{x0394}U\,\text{Im}\left\{\frac{4F(1/L)}{F^{\prime }(-1/L)}+2U_{c}^{2}\frac{m_{e}}{m_{i}}\frac{1}{F^{\prime }(-1/L)}\displaystyle \frac{1}{\unicode[STIX]{x1D713}^{2}}\int _{x_{a}}^{x_{b}}\unicode[STIX]{x1D719}(x)\,\text{d}x\right\}. & & \displaystyle\end{eqnarray}$$

The first term is only a function of the hole shape $\tilde{\unicode[STIX]{x1D719}}$ while the second term depends on hole size $\unicode[STIX]{x1D713}$ and the mass ratio $m_{i}/m_{e}$ . However, this second term is not important except for extremely shallow EHs such that $\unicode[STIX]{x1D713}\ll 1$ . For example, for a Schamel type of EH of size $\unicode[STIX]{x1D713}=0.1$ and $m_{i}/m_{e}=1836$ , the magnitude of the second term is approximately 4 % of the first one. It is thus a good approximation that for not too shallow EHs we have

(3.44) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FE}\simeq \unicode[STIX]{x0394}U\,\text{Im}\left\{\frac{4F(1/L)}{F^{\prime }(-1/L)}\right\}=\unicode[STIX]{x0394}U\,\text{Im}\left\{\frac{4F(1/L)}{F^{\prime }(1/L)}\right\}=-\unicode[STIX]{x0394}|U|\,\text{Im}\left\{\frac{4F(1/L)}{F^{\prime }(1/L)}\right\}. & & \displaystyle\end{eqnarray}$$

We define $\unicode[STIX]{x0394}|U|=|U|-U_{c}$ and the imaginary part of $F$ is odd so its derivative is even: $\text{Im}(F^{\prime }(-1/L))=\text{Im}(F^{\prime }(1/L))$ . This growth rate scales linearly with $\unicode[STIX]{x0394}|U|$ . If the hole potential shape is of Schamel type, a numerical evaluation of the constants gives $\unicode[STIX]{x1D6FE}\simeq -\unicode[STIX]{x0394}|U|/1.74$ for $|U|$ evaluated in $c_{s}$ and $\unicode[STIX]{x1D6FE}$ evaluated in $\unicode[STIX]{x1D714}_{pi}$ . The instability grows fast once $|U|$ is slower than $U_{c}$ .

Figure 7. Instability growth rate $\unicode[STIX]{x1D6FE}$ as a function of $\unicode[STIX]{x0394}U$ . The line represents (3.44) for fixed hole shape. Its uncertainty bands represent the small variation of shape from one run to another, giving uncertainty in the comparison. The triangles are obtained from solving numerically the full eigenmode equation ${\dot{P}}_{i}/{\dot{P}}_{e}+1=0$ using the PIC potential output. Circles are the growth rate observed in PIC runs.

In the PIC simulations, we measured the growth rate of unstable velocity oscillations by fitting an exponential growth model to it. We used the unstable runs with different counter-streaming ion velocity and $\unicode[STIX]{x1D713}\sim 0.8$ , in which $\unicode[STIX]{x0394}|U|$ can be precisely measured. The growth rates and the error bars are obtained from the regression. They are compared with the expanded linear solution in (3.44) and numerical solutions of the full eigenmode equation. The results are shown in figure 7. Our analytic theory agrees with the observed instability growth rate for the weakly unstable cases up to $\unicode[STIX]{x1D6FE}\sim 0.1$ and the linear expansion gives a very good approximation to the solution of the full eigenmode equation. The growth rates for strongly unstable cases are less than the predicted values. We attribute this discrepancy to our inability to observe the growth rate in a truly linear stage when $\unicode[STIX]{x1D6FE}$ is large. We discuss these nonlinear effects in § 4.

3.4 Finite ion temperature

So far, we have treated the ions as a cold beam. In reality, they have a thermal velocity spread. For an ion of velocity $v$ , the EH has a velocity $U-v$ in its frame. Let us consider the ion thermal speed to be small compared to $|U|$ such that $v_{\text{th,i}}\ll |U|$ and there are no reflected ions by the hole potential. As $U_{c}$ is several times $c_{s}$ , this assumption holds approximately even when $T_{i}\geqslant T_{e}$ . We can integrate the contributions from ions of different velocities to get the total ${\dot{P}}_{i}$

(3.45) $$\begin{eqnarray}\displaystyle {\dot{P}}_{i}(\unicode[STIX]{x1D714},U,\unicode[STIX]{x1D719},v_{\text{th,i}})=m_{i}\dot{U}\int _{-\infty }^{\infty }f_{\infty ,i}(v)\frac{\unicode[STIX]{x1D713}^{2}}{(U-v)^{4}}F\left(\frac{\unicode[STIX]{x1D714}}{U-v}\right)\,\text{d}v. & & \displaystyle\end{eqnarray}$$

We apply a Taylor series expansion to (3.45) assuming $|v|\ll |U|$ . Consider $f_{\infty ,i}(v)$ to be a Maxwellian and only the even-order moments of $v$ survive after the integration over velocity. This expansion gives

(3.46) $$\begin{eqnarray}\displaystyle {\dot{P}}_{i}(\unicode[STIX]{x1D714},U,\unicode[STIX]{x1D719},v_{\text{th,i}}) & = & \displaystyle n_{\infty }m_{i}\dot{U}\frac{\unicode[STIX]{x1D713}^{2}}{U^{4}}\left\{F\left(\frac{\unicode[STIX]{x1D714}}{U}\right)+F_{2}\left(\frac{\unicode[STIX]{x1D714}}{U}\right)\left(\frac{v_{\text{th,i}}}{U}\right)^{2}+O\left(\left(\frac{v_{\text{th,i}}}{U}\right)^{4}\right)\right\}\nonumber\\ \displaystyle & = & \displaystyle n_{\infty }m_{i}\dot{U}\frac{\unicode[STIX]{x1D713}^{2}}{U^{4}}\left\{F_{\text{th,i}}\left(\frac{\unicode[STIX]{x1D714}}{U}\right)+O\left(\left(\frac{v_{\text{th,i}}}{U}\right)^{4}\right)\right\},\end{eqnarray}$$

where

(3.47) $$\begin{eqnarray}\displaystyle F_{2}\left(\frac{\unicode[STIX]{x1D714}}{U}\right)=10F\left(\frac{\unicode[STIX]{x1D714}}{U}\right)+5F^{\prime }\left(\frac{\unicode[STIX]{x1D714}}{U}\right)\left(\frac{\unicode[STIX]{x1D714}}{U}\right)+\frac{1}{2}F^{\prime \prime }\left(\frac{\unicode[STIX]{x1D714}}{U}\right)\left(\frac{\unicode[STIX]{x1D714}}{U}\right)^{2}. & & \displaystyle\end{eqnarray}$$

We have right now ${\dot{P}}_{i}/{\dot{P}}_{e}(\unicode[STIX]{x1D714},U,\unicode[STIX]{x1D719},v_{\text{th,i}})=-F_{\text{th,i}}(\unicode[STIX]{x1D714}/U)/G(U)$ to the leading order in $|v_{\text{th,i}}/U|$ . It suffices to substitute $F$ with $F_{\text{th,i}}$ in our previous analysis and everything follows as before.

The leading-order term of the finite ion temperature correction is second order in $|v_{\text{th,i}}/U|$ . The effect of finite ion temperature can be visualized through the $F_{\text{th,i}}$ contours. In figure 8, we show the $F_{\text{th,i}}$ contours evaluated on the real axis for different values of $|v_{\text{th,i}}/U|$ using Schamel type of EH potential. The most salient effect is that the finite ion temperature moves the crossing point $C(\tilde{\unicode[STIX]{x1D719}})$ outwards, resulting in a higher value for $U_{c}$ . Because of the leading $U^{4}$ term in $G(U)$ , the resulting change in $U_{c}$ is actually relatively small. In terms of $U_{c}$ , this correction is ${\sim}5\,\%$ when $|U|=5v_{\text{th,i}}$ , it grows to ${\sim}10\,\%$ for $|U|=4v_{\text{th,i}}$ and ${\sim}20\,\%$ when $|U|=3v_{\text{th,i}}$ . This property holds similarly for other EH potential models such as the Gaussian. The same trend is noticed in our PIC simulation. A higher ion temperature $T_{i}$ leads to a slightly higher threshold velocity $U_{c}$ .

When $|U|<3v_{\text{th,i}}$ , the ion reflection from the EH potential becomes important and can no longer be neglected in the global momentum balance. The unbalanced scattering of ions tends to accelerate the EH to a higher velocity in the ion frame (Hutchinson & Zhou Reference Hutchinson and Zhou2016). With hole pushing technique, we were able to explore the situation with the presence of mild ion reflection from the hole potential, the instability is still observed in these cases. After the onset of instability, we stopped hole pushing and observed the hole velocity to oscillate while its mean velocity accelerates due to the ion reflection. Resonant ion effects such as ion Landau damping or reflection only become important when $|U|\sim v_{\text{th,i}}$ , which requires the ions to be extremely hot. Our analysis holds well for the usual range of $T_{e}/T_{i}$ in space plasmas.

Figure 8. Finite ion temperature effect on the $F$ contour for a Schamel type of EH. The contour shape is approximately preserved while its size grows with a larger $T_{i}$ .