Plasma electron hole oscillatory velocity instability

TL;DR

This paper introduces a new instability in electron holes caused by their interaction with passing ions, leading to oscillations, wave emissions, and density perturbations, supported by theoretical analysis and PIC simulations.

Contribution

It presents a novel theoretical framework for electron hole-ion interactions and validates it with high-fidelity PIC simulation results.

Findings

Electron hole velocity oscillations become unstable under certain conditions.

Unstable oscillations emit ion-acoustic waves and cause hole decay.

Theoretical predictions match PIC simulation observations.

Abstract

In this paper, we report a new type of instability of electron holes (EHs) interacting with passing ions. The nonlinear interaction of EHs and ions is investigated using a new theory of hole kinematics. It is shown that the oscillation in the velocity of the EH parallel to the magnetic field direction becomes unstable when the hole velocity in the ion frame is slower than a few times the cold ion sound speed. This instability leads to the emission of ion-acoustic waves from the solitary hole and decay in its magnitude. The instability mechanism can drive significant perturbations in the ion density. The instability threshold, oscillation frequency and instability growth rate derived from the theory yield quantitative agreement with the observations from a novel high-fidelity hole-tracking Particle-In-Cell (PIC) code.

Figures

Equations

$$\{\begin{array}{rl}& v_0(x)=-\frac{U}{|U|}\sqrt{U^2-2\varphi (x)}\\ & n_0(x)=n_{\infty }\frac{|U|}{\sqrt{U^2-2\varphi (x)}}\end{array}$$

$${\dot{P}}_e=-m_e\dot{U}{n}_{\infty }{\int }_{x_a}^{x_b}h(\sqrt{\varphi (x)})+1-\frac{|U|}{\sqrt{U^2-2\varphi (x)}}dx,$$

$$\begin{array}{rl}\frac{{\dot{P}}_{io}}{m_i}& =n_{bf}{v}_{bf}{\bar{v}}_{bf}-n_{af}{v}_{af}{\bar{v}}_{af}\\ & =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}\\ & =n_{bf}{v}_{bf}({\bar{v}}_{bf}-{\bar{v}}_{af})+(n_{bf}{v}_{bf}-n_{af}{v}_{af}){\bar{v}}_{af}.\end{array}$$

$$\frac{1}{2}{v}^2+\varphi (x)+\int \dot{U}vdt=\mathrm{Constant}.$$

$$v_{bf}^2-v_{as}^2+2{\int }_{t_s}^{t_f}\dot{U}vdt=0.$$

$$v_{bf}-v_{as}$$

$$\simeq {\int }_{x_a}^{x_b}\frac{\dot{U}(t(x,x_b))}{U}dx.$$

$${\bar{v}}_{bf}-{\bar{v}}_{af}$$

$$\simeq {\int }_{x_a}^{x_b}(\frac{1}{U}+\frac{1}{v_0})\dot{U}dx.$$

$$n_{bf}{v}_{bf}({\bar{v}}_{bf}-{\bar{v}}_{af})\simeq n_{\infty }{\int }_{x_a}^{x_b}(-1-\frac{U}{v_0})\dot{U}dx.$$

$$n_{as}{v}_{as}\delta t_{as}=n_{bf}{v}_{bf}\delta t_{bf},$$

$$\displaystyle\begin{aligned} n_{bf}v_{bf}-n_{af}v_{af}&\,=\,n_{bf}v_{bf}-n_{bf}v_{bf}\frac{\delta t_{bf}}{\delta t_{as}}\frac{v_{af}}{v_{as}}\\ &\,=\,n_{bf}v_{bf}\bigg{[}1-\frac{\delta t_{bf}}{\delta t_{as}}(1+\frac{v_{af}-v_{as}}{v_{as}})\bigg{]}\\ &\,\simeq n_{bf}v_{bf}(-\frac{dt_{ab}}{dt_{s}}-\frac{v_{af}-v_{as}}{v_{as}}),\end{aligned}$$

$$\frac{d{t}_{ab}}{d{t}_s}$$

$$\displaystyle=\,\int_{x_{a}}^{x_{b}}-\frac{1}{v^{3}}\frac{\partial(v^{2}/2)}{\partial t_{s}}\bigg{|}_{x}\,dx$$

$$\displaystyle=\,\int_{x_{a}}^{x_{b}}-\frac{1}{v^{3}}\frac{\partial}{\partial t_{s}}\bigg{|}_{x}\left(v_{as}^{2}/2-\phi(x)-\int_{x_{a}}^{x}\dot{U}\,dx_{1}\right)\,dx$$

$$\simeq {\int }_{x_a}^{x_b}\frac{1}{v^3}[-U\dot{U}(t_s)+{\int }_{x_a}^x\ddot{U}d{x}_1]dx.$$

$$\dot{U}(t_s)-\frac{1}{U}{\int }_{x_a}^x\ddot{U}d{x}_1$$

$$\displaystyle\simeq\,\dot{U}(t_{s})-\bigg{[}\frac{\dot{U}(t)v_{0}(x)}{U}-\frac{\dot{U}(t_{s})v_{0}(x_{a})}{U}\bigg{]}+\frac{1}{U}\int_{x_{a}}^{x}\dot{U}\frac{dv_{0}}{dx_{1}}\,dx_{1}$$

$$\simeq -\frac{\dot{U}(t)v_0(x)}{U}-\frac{1}{U}{\int }_{x_a}^x\frac{\dot{U}{\varphi }^{\prime }(x_1)}{v_0(x_1)}d{x}_1.$$

$$n_{bf}{v}_{bf}(-\frac{d{t}_{ab}}{d{t}_s}-\frac{v_{af}-v_{as}}{v_{as}})$$

$$\displaystyle\simeq\,n_{\infty}\int_{x_{a}}^{x_{b}}\frac{U^{2}}{v_{0}^{3}}\bigg{[}\frac{\dot{U}v_{0}}{U}+\frac{1}{U}\int_{x_{a}}^{x}\frac{\dot{U}\phi^{\prime}}{v_{0}}\,dx_{1}\bigg{]}+\frac{\dot{U}}{v_{0}}\,dx$$

$$\displaystyle=\,n_{\infty}\int_{x_{a}}^{x_{b}}\frac{U}{v_{0}^{2}}\dot{U}+\frac{U}{v_{0}^{3}}\bigg{(}\int_{x_{a}}^{x}\frac{\dot{U}\phi^{\prime}}{v_{0}}\,dx_{1}\bigg{)}+\frac{\dot{U}}{v_{0}}\,dx.$$

$$\dot{P}_{io}=m_{i}n_{\infty}\int_{x_{a}}^{x_{b}}(-1-2\frac{U}{v_{0}}-\frac{U^{2}}{v_{0}^{2}})\dot{U}-\frac{U^{2}}{v_{0}^{3}}\bigg{(}\int_{x_{a}}^{x}\frac{\dot{U}\phi^{\prime}}{v_{0}}\,dx_{1}\bigg{)}\,dx.$$

$${\dot{P}}_{ic}=m_i{\int }_{x_a}^{x_b}\frac{\partial (nv)}{\partial t}dx+m_i\dot{U}{\int }_{x_a}^{x_b}ndx.$$

$${\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)){\varphi }^{\prime \prime }(x_2)d{x}_{2}d{x}_{1}dx.$$

$${\dot{P}}_i(\omega )+{\dot{P}}_e(\omega )=0.$$

$$\displaystyle\begin{aligned} \displaystyle\frac{\dot{P}_{i}}{\dot{P}_{e}}(\omega,U,\phi)=&-\displaystyle\frac{m_{i}}{m_{e}}\bigg{[}\displaystyle\int_{x_{a}}^{x_{b}}\left(-1-2\displaystyle\frac{U}{v_{0}(x)}-\displaystyle\frac{U^{2}}{v_{0}^{2}(x)}\right)\exp\left(i\omega\displaystyle\int_{x}^{x_{b}}\displaystyle\frac{dx_{3}}{v_{0}(x_{3})}\right)\,dx\\ &\hskip 28.45274pt-\displaystyle\int_{x_{a}}^{x_{b}}\displaystyle\frac{U^{2}}{v_{0}^{3}(x)}\displaystyle\int_{x_{a}}^{x}\displaystyle\frac{\phi^{\prime}(x_{1})}{v_{0}(x_{1})}\exp\left(i\omega\displaystyle\int_{x_{1}}^{x_{b}}\displaystyle\frac{dx_{3}}{v_{0}(x_{3})}\right)\,dx_{1}\,dx\\ &\hskip 28.45274pt-\displaystyle\int_{x_{a}}^{x_{b}}\displaystyle\int_{x_{a}}^{x}\displaystyle\frac{U}{v_{0}^{3}(x_{1})}\displaystyle\int_{x_{a}}^{x_{1}}\exp\left(i\omega\displaystyle\int_{x_{2}}^{x}\displaystyle\frac{dx_{3}}{v_{0}(x_{3})}\right)\phi^{\prime\prime}(x_{2})\,dx_{2}\,dx_{1}\,dx\bigg{]}\\ &\bigg{/}\left(\displaystyle\int_{x_{a}}^{x_{b}}h(\sqrt{\phi(x)})+1-\frac{\lvert U\rvert}{\sqrt{U^{2}-2\phi(x)}}\,dx\right).\end{aligned}$$

$$\frac{{\dot{P}}_i}{{\dot{P}}_e}(\omega ,U,\varphi )\simeq -\frac{m_i}{m_e}\frac{{\psi }^2}{U^4}\frac{4i\frac{\omega }{U}I(\frac{\omega }{U})+i\frac{{\omega }^2}{U^2}{I}^{\prime }(\frac{\omega }{U})-3I_0}{{\int }_{x_a}^{x_b}h(\sqrt{\varphi (x)})-\frac{\varphi (x)}{U^2}dx},$$

$$I_0={\int }_{x_a}^{x_b}\tilde{\varphi }(x{)}^{2}dx,$$

$$I(\frac{\omega }{U})={\int }_{x_a}^{x_b}{\int }_{x_a}^y\tilde{\varphi }(x)\tilde{\varphi }(y)\exp (i\frac{\omega (x-y)}{U})dxdy,$$

$$F(\frac{\omega }{U})=4i\frac{\omega }{U}I(\frac{\omega }{U})+i\frac{{\omega }^2}{U^2}{I}^{\prime }(\frac{\omega }{U})-3I_0,$$

$$\displaystyle G(U)\,=\,\displaystyle\frac{m_{e}}{m_{i}}\displaystyle\frac{1}{\psi^{2}}\bigg{[}U^{4}\displaystyle\int_{x_{a}}^{x_{b}}h(\sqrt{\phi(x)})\,dx-U^{2}\displaystyle\int_{x_{a}}^{x_{b}}\phi(x)\,dx\bigg{]}.$$

$$\frac{{\dot{P}}_i}{{\dot{P}}_e}(\omega ,U)=-\frac{F(\frac{\omega }{U})}{G(U)}.$$

$$-\frac{C(\tilde{\varphi })}{G(-U_c)}=-1.$$

$$U_c^4{\int }_{x_a}^{x_b}h(\sqrt{\varphi (x)})dx-U_c^2{\int }_{x_a}^{x_b}\varphi (x)dx-{\psi }^2\frac{m_i}{m_e}C(\tilde{\varphi })=0.$$

$$U_c^2=\frac{{\int }_{x_a}^{x_b}\varphi (x)dx+\sqrt{{\left({\int }_{x_a}^{x_b}\varphi (x)dx\right)}^2+4{\psi }^2\frac{m_i}{m_e}C(\tilde{\varphi }){\int }_{x_a}^{x_b}h(\sqrt{\varphi (x)})dx}}{2{\int }_{x_a}^{x_b}h(\sqrt{\varphi (x)})dx}.$$

Full text

I Introduction

A plasma electron hole is a solitary vortex structure of trapped electrons in the electron phase-space. The deficit of electrons on the trapped orbits gives rise to a coherent self-trapping structure. We have performed 1-D Particle-In-Cell simulations of EHs with different velocities relative to the ions. Oscillations in the velocity of a slowly moving EH are observed to be unstable. Using a novel theory of hole kinematics that treats the EH as a holistic object, we give a physical explanation for this newly observed instability. The analytical results agree quantitatively with our simulation observations. EHs occur widely in nature and they are important to many space physics phenomena.

The EH was first discovered as a nonlinear equilibrium of a collisionless Vlasov-Poisson plasma under the name of Bernstein-Greene-Kruskal mode [1]. The study of these nonlinear structures in plasma gained increasing interest after space probe measurements confirmed their wide-spread existence in the Earth’s auroral zone [2], magnetosphere [3] and in the solar wind [4]. They are implicated in magnetic reconnection [5], electron acceleration [6], collisionless shocks [7] and other important plasma dynamics in space. EHs are also detected in the laboratory plasma during magnetic reconnection [8] and beam injection [9].

The early theoretical research of EHs neglected the ion dynamics for simplicity and considered them as a uniform neutralizing background [10]. Later, Saeki et al. [11] showed using PIC simulations that an EH can be disrupted by ions when its velocity is slower than the ion sound speed. Eliasson et al. [12] reported that a standing EH can be ejected by the ion density cavity it created and is attracted to ion density maxima. The recent observations of “slow” EHs also suggest a more important role for the ions. EHs traveling with the ion sound speed $c_{s}$ have been recently reported at a magnetic reconnection site [13] and the magnetopause [14] measured by the Cluster spacecrafts, a velocity much slower than what was frequently observed before ($\sim v_{th,e}$). The authors suggested that Buneman instability resulting from dynamic reconnecting current sheets generated these slow EHs. Schamel gave an upper limit [15] for the speed of the EHs by a structural existence argument. How slow can an EH travel? Saeki et al. [11] briefly touched upon this question by deriving the nonlinear dispersion relation using the Sagdeev pseudo-potential method. However, the nonlinear dispersion relation is based on the existence of a stationary solution of which the stability is not guaranteed. An EH can experience different kinds of instabilities in higher dimensions, e.g. the whistler instability [16] in the strongly magnetized case and the transverse instability [17] in the weakly magnetized case. These instability mechanisms don’t involve ion dynamics. Ion-acoustic wave radiation from a solitary structure in plasma has been studied in the case of Langmuir soliton [18] and ion hole [19]. Dyrud et al. [20] reported the observation of ion-acoustic waves emitted from a chain of electron holes in PIC simulations. Dokgo et al. [21] reported the generation of coherent ion-acoustic solitary waves from an EH as it propagates from a lower plasma density region to a higher plasma density region. In this paper, we show that the ion dynamics is important for the stability of an EH even in the 1-D equilibrium, causing an oscillatory velocity instability for slow electron holes to decay into ion-acoustic waves. This discovery suggests that the solitary solutions constructed using the Sagdeev pseudo-potential method can be unstable to small perturbations in its velocity when the ion dynamics becomes important.

The instability mechanism discussed in this paper is closely related to the velocity thus the kinematics of an EH. The EH kinematics has been studied by the authors [22] treating the EH as a composite object. When an EH accelerates, ion and electron plasma momentum changes both inside and outside the EH. The EH should move in a way to conserve the total plasma momentum. This quasi-particle approach using the global momentum conservation can be found in the early EH theory developed by Dupree et al. [23] Our theory is a major improvement over Dupree’s and was successful in explaining quantitatively the dynamics of EHs observed in PIC simulations [24], such as the transient self-acceleration and the “hole pushing/pulling” effect due to steady-state hole momentum coupling to the ions. In this paper, we extend our theory to the frequency domain and use multiple-scale analysis to give a mathematically rigorous treatment of the instability.

The paper is organized as follows: in Section II, we report the observational details of this instability from our PIC simulation. In Section III, a first principle analytic calculation using hole kinematics theory is presented, the instability boundary, unstable mode frequency and growth rate are analytically derived and compared with the PIC simulation observation. In Section IV, we are going to discuss the nonlinear stage of the instability and its potential implication in space plasma. Section V is the summary.

II PIC observation of the instability

The simulations are performed using a 1-D electrostatic PIC code with fully kinetic ions, which is designed to study highly-resolved EH dynamics. A solitary EH is created in our simulation using an electron phase-space density deficit as the initial seed [24]. The thermal noise in our PIC simulation is controlled by using more than $10^{6}$ particles per cell. There are $\sim 10$ cells per Debye length to resolve the detail of an EH. We performed box simulation of a solitary EH with a computation domain that can self-consistently follow the EH motion [24]. The size of the computation domain is only $\sim 50\lambda_{\mathrm{De}}$ by virtue of using hole tracking. Our hole-tracking PIC allows us to study the detail of hole motion with reasonable computational cost. Once a steady-state EH is obtained in our simulation, we apply an artificial slow ion acceleration to slowly “push” the EH with ions so that it slows down in the ion frame. This “pushing” process has been demonstrated to be quasi steady-state and reversible [24]. We discovered when performing these “pushing” runs that there is a limit velocity of the EH in the ion frame, below which the EH becomes unstable. This threshold velocity is well above other physical limit of the system such as the ion reflection velocity limit.

An example of the instability observed in our PIC simulation is presented in Figure 1. On the top left, we show the characteristic solitary potential structure of a stable EH that extends over several Debye lengths. In a steady-state EH, the ions are slowed down by the hole potential and their density is slightly higher inside the EH. This ion density compressional pulse is the ion-acoustic soliton attached to the phase-space EH described by Saeki et al.[11] It is clearly visible inside the stable EH shown on the left. When the EH slows down in the ion frame, the oscillation amplitude in its velocity begins to grow once its speed is slower than a threshold value. The EH potential and the ion density at a later time step ($\omega_{pe}t=2925$) after the instability has grown are presented on the right for comparison. The EH keeps its potential shape while its velocity oscillates. The down stream ion density becomes unsteady as the velocity oscillation amplitude grows. Unsteady ion density perturbations are emitted from the EH after the instability onset. The perturbations propagate in the ion frame with the ion sound speed, mainly in the opposite direction to the EH velocity. The EH velocity was obtained from the hole-tracking module in our code [24] and a low pass filter has been applied to it to filter out the statistical noise. The “hole pushing” can be turned off at any moment before the instability onset and the EH will enter a stable steady-state with the same velocity, but not after the EH is slower than the threshold velocity.

A similar phenomenon also happens in a plasma with counter-streaming ions. We initialize an EH at rest ($U=0$ initially) in the lab frame on top of the electron distribution with counter-streaming ions traveling at $\pm v_{i}$ in the lab frame. This time we don’t need to apply any special technique such as the hole pushing and hole tracking. The initialization will naturally favor the formed EH to stay at $U=0$ and we can do a regular box simulation of the EH with a static domain. We observe that there is a minimum value of $v_{i}$ below which the system is unstable. A case of the observed instability is shown in Figure 2. The self-consistently formed EH is unstable. Perturbations in its velocity grow exponentially. Ion-acoustic perturbations grow and are emitted from both sides of the EH because the ions are counter-streaming. The simulation was performed with warm ions $T_{i}=T_{e}$, corresponding to a case where ion-acoustic waves are strongly Landau damped, ion-ion type of instability and Buneman instability are ruled out by the simulation parameter setting. For this particular case shown in Figure 2, using a slightly higher $v_{i}=7c_{s}$ can stabilize the instability. We shall see later in the paper that from the stability point of view, the counter-streaming ion case is equivalent to the single ion stream case.

We have repeated the numerical experiments with different parameter settings. It is observed that for the same EH, the instability only depends significantly on the relative velocity between the EH and the ions $U-v_{i}$. The EH velocity with respect to the bulk electrons $U$ and the relative velocity drift between the ions and the bulk electrons $v_{i}$ alone have no significant influence on the instability onset and its growth. We have performed simulations with an electron-to-ion temperature ratio $T_{e}/T_{i}$ from $20$ down to $0.5$. We observe that this velocity instability clearly persists in the regime $T_{e}\leq T_{i}$, where ion-acoustic type of instability is unexpected. A hotter ion population leads to a higher threshold value of $U-v_{i}$ and damped ion-acoustic wings. We shall discuss the finite ion temperature effect on the instability in detail later in our paper.

Our PIC simulation contains no imprecision beyond the statistical noise which was largely mitigated by using a large number of particles. It is clear that a self-consistent solitary solution with a complete ion response can be constructed in the case of instability using the Sagadeev pseudo potential or the BGK approach. Our PIC code actually does this by solving Poisson’s equation numerically. However, once this steady-state solitary solution is allowed to evolve in time, it becomes unstable. The characteristics of this instability does not fit into any existing linear plasma stability theory. The core of this problem is very nonlinear because of the strong particle trapping nonlinearity in the EH. We need to adopt a new approach to analyze its mechanism.

III Hole velocity stability deduced from kinematics

III.1 Frequency response of the momentum rate of change

To analyze this instability, we first consider the steady-state solution of an EH. The steady-state EH potential $\phi(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 center of the hole so that both $\phi(x)$ and its derivatives vanish at these limits: $\phi(x_{a})=\phi(x_{b})=\phi^{\prime}(x_{a})=\phi^{\prime}(x_{b})=0$. The ions and the bulk electrons are assumed to be Maxwellian at rest in the lab frame with their background density being denoted by $n_{\infty}$. The EH moves at a velocity $U$ in the lab 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 $\lambda_{\mathrm{De}}$, $\phi$ 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/\omega_{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:

[TABLE]

In our previous paper [22], 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 in a control volume that extends from $x_{a}$ to $x_{b}$. 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 $\phi(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\mathrm{Re}(\exp[-i\omega t])$. The frequency $\omega$ 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 [22] of $\dot{P}_{e}$ for $U\ll v_{th,e}$ with a full ion response correction [24]

[TABLE]

where $h(\chi)\,=\,-\frac{2}{\sqrt{\pi}}\chi+(2\chi^{2}-1)e^{\chi^{2}}\operatorname{erfc}(\chi)+1$. The ion response correction accounts for the ion density accumulation inside the EH. The exact shape of the trapped electron distribution doesn’t appear 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

[TABLE]

The term $\bar{v}_{bf}-\bar{v}_{af}$ represents the “jetting” effect [22] 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

[TABLE]

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

[TABLE]

The equilibrium velocity for $\dot{U}=0$ is $v_{0}(x)=-\frac{U}{\lvert U\rvert}\sqrt{U^{2}-2\phi(x)}$. The idea is to do an expansion around the equilibrium orbit. Perturbation expansion gives $v=v_{0}+v_{1}$ and $\lvert v_{1}\rvert/\lvert v_{0}\rvert\sim\lvert t_{ab}\dot{U}/U\rvert\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}$

[TABLE]

where $t(x,x_{b})=t_{f}-\int_{x}^{x_{b}}\frac{du}{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

[TABLE]

Thus we have obtained the “jetting” effect due to the acceleration of the EH to the relevant order and the first term in Eqn. (3) can be evaluated as

[TABLE]

To calculate the second term in Eqn. (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}$

[TABLE]

where $\delta t_{as}$ and $\delta t_{bf}$ are two infinitesimal time durations for the same ion fluid element to enter and exit the control volume. They are related by the derivative of transit time $t_{ab}$ w.r.t. the starting time $t_{s}$: $\delta t_{bf}\simeq\delta t_{as}(1+\frac{dt_{ab}}{dt_{s}})$. Thus, to the leading order,

[TABLE]

where we used the constancy of inflow density $n_{as}=n_{af}=n_{\infty}$ and that $dt_{ab}/dt_{s}$ is of the same order as $(v_{af}-v_{as})/v_{as}$. The derivative $dt_{ab}/dt_{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}}\frac{dx}{v}$ and the conservation law of Eqn. (4)

[TABLE]

$\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 $\partial/\partial t_{s}$ and $\partial/\partial t$ as $dt=dt_{s}(1+\mathcal{O}(t_{ab}\dot{U}/U))$ for $x_{1}$ fixed. We can further get rid of the $\dot{U}(t_{s})$ term in Eqn. (11) performing integration by parts

[TABLE]

We used here $dv_{0}/dx_{1}=-\phi^{\prime}/v_{0}$, where $\phi^{\prime}$ is the spatial derivative of $\phi$. Combining Eqns. (11) and (12), we can evaluate the right hand side of Eqn. (10) as

[TABLE]

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 Eqns. (8) and (13). It is first order in $t_{ab}\dot{U}/U$, thus linear in $\dot{U}$

[TABLE]

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

[TABLE]

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 the same order as $\dot{P}_{io}$ in Eqn. (14) can be expressed as

[TABLE]

We combine Eqs. (14) and (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 $\omega$:

[TABLE]

The imaginary part of $\omega$ determines the stability of the corresponding eigenmode.We apply Nyquist stability analysis [25] 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

[TABLE]

$\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\psi$ with $\psi$ being the maximum of $\phi$. This approximation is very well satisfied at the onset of instability observed in our simulation. The leading term of the expanded form is

[TABLE]

where

[TABLE]

with $\tilde{\phi}(x)=\phi(x)/\psi$ being the normalized potential profile function. This leading term is second order in the small expansion parameter $2\psi/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}(\omega)$ can be evaluated numerically for real frequencies $\omega$ using for example the widely cited Schamel type of EH potential $\phi(x)=\psi\operatorname{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 the panel (a) of Figure 4. The number of encirclements of the point $-1$ in the complex plane by the $\dot{P}_{i}/\dot{P}_{e}(\omega)$ contour gives the number of unstable $\omega$ solutions to the eigenmode Eqn. (17). There is a critical speed $U_{c}$ for $|U|$ below which the system is unstable. The leading order term is within a few percent 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.

The $\dot{P}_{i}/\dot{P}_{e}(\omega)$ 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

[TABLE]

We have therefore

[TABLE]

We look for the critical speed $U_{c}$ for a given equilibrium hole potential $\phi$ such that the $\dot{P}_{i}/\dot{P}_{e}(\omega)$ contour crosses the point $-1$ in the complex plane. $F$ depends on $U$ through $\omega/U$, it gives the same contour for different $U$ values when $\omega$ is evaluated on the real axis, although at different $\omega$ 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}(\omega)$ contour through $G$. This property is demonstrated in Figure 4. The identical $F(\omega/U)$ contour for different values of $U$ using a Schamel type of EH potential is shown in the panel (b) of Figure 4. We denote its intersection with the positive real axis by $C(\tilde{\phi})$. The existence of this intersection $C(\tilde{\phi})$ is guaranteed for a general class of admissible hole potential $\phi(x)$, which we will show later in this paper. The critical speed $U_{c}$, below which the system is unstable, satisfies an equation

[TABLE]

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

[TABLE]

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

[TABLE]

Now we calculate the oscillation frequency of the unstable eigenmode. At marginal instability, the $\dot{P}_{i}/\dot{P}_{e}(\omega)$ 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

[TABLE]

where

[TABLE]

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

[TABLE]

where $\tilde{\Phi}$ is the modulus of the Fourier transform of $\tilde{\phi}$. Thus $\mathrm{Im}(\dot{P}_{i}/\dot{P}_{e})(\omega)=0$ gives

[TABLE]

Eqn. (31) admits three real solutions for $\omega$, one is the trivial $\omega=0$, the other two solutions are given by the equation

[TABLE]

Since $\tilde{\phi}$ is real, we have $\tilde{\Phi}$ is an even function and Eqn. (32) gives two solutions that have the opposite sign. We define the positive $\omega$ solution of Eqn. (32) as $\omega_{0}$ and a length scale $L$

[TABLE]

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

[TABLE]

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 $\omega_{c}$ and the crossing point $C(\tilde{\phi})$ are guaranteed by the continuous differentiablility of the EH potential $\phi(x)$. A physical EH potential $\phi(x)$ should possess a second derivative as it satisfies Poisson’s equation $\phi^{\prime\prime}(x)+\rho(x)/\epsilon_{0}=0$, and a physical $\rho(x)$ should have bounded variation. This smoothness requirement constrains the asymptotic behavior of its Fourier transform. Function $\tilde{\Phi}(p)$ decays at least as fast as $p^{-3}$ when $p\rightarrow\infty$ [26]. Thus we have $-\tilde{\Phi}^{\prime}(p)/2\tilde{\Phi}(p)=-\frac{1}{2}d\ln(\tilde{\Phi}(p))/dp\geq 3/2p>1/p$ as $p\rightarrow\infty$. While as $p\rightarrow 0$, we have $-\tilde{\Phi}^{\prime}(p)/2\tilde{\Phi}(p)\rightarrow-\tilde{\Phi}^{\prime}(0)/2\tilde{\Phi}(0)\ll 1/p$. Solutions are guaranteed for Eqn. (32). The behavior of the $F$ contour is as follows. As $\omega/U\rightarrow 0^{+}$, we have $\mathrm{Re}(F(\omega/U))\rightarrow-3I_{0}<0$ and $\mathrm{Im}(F(\omega/U))\rightarrow 0^{+}$. Asymptotic analysis as $\omega/U\rightarrow+\infty$ gives $\mathrm{Re}(F(\omega/U))\rightarrow 0^{+}$ and $\mathrm{Im}(F(\omega/U))\rightarrow 0^{-}$. In other words, with $\omega/U$ increasing from [math] 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{\phi})$ shown in Figure fig. 4, which is crucial to this instability, is an universal feature for all physically admissible hole potentials $\phi(x)$ and we have $C(\tilde{\phi})=F(1/L)$. Contours without a crossing can be obtained only from unphysical hole shapes. For example, $\tilde{\phi}=\exp(-\lvert x\rvert/\lambda)$ does not give a crossing point and is therefore stable; but it is unphysical, as the electric field is undefined at $x=0$.

Having obtained the analytic solution for the critical speed and the unstable oscillation frequency, we now compare these results with our PIC 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 $\phi(x)$ right before the instability onset from the same run. The electrostatic potential output $\phi(x)$ from our PIC simulation is used to construct the $F$ contour numerically from Eqn. (22) and find its crossing point $C(\tilde{\phi})$. We use the formula for $U_{c}$ in Eqn. (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 another 111Our hole-tracking PIC simulation produces relatively low-noise and highly resolved EH potential. We applied some post processing to the PIC potential output to make the numerical calculation more accurate. In our analysis, $\phi(x)$ is considered to fall to zero far away from the hole center. However, there is always some non-zero intrinsic statistical noise in the PIC simulation. In post processing, we find the positions where the electric field first becomes zero outside the hole center and consider them to be the limits of the hole spatial extent. The values of $\phi(x)$ beyond these limits are forced to decay to zero by multiplying a Debye decaying exponential to them. We use this slightly smoothed $\phi(x)$ in our numerical calculation of $U_{c}$.. The results are presented in Figure 5. The solid lines are obtained using Eqn. (27) assuming a Schamel type of EH potential. This solution’s asymptotic behavior comes from the special function $h(\chi)$ [22]: $h(\chi)\rightarrow\chi^{2}-\frac{8}{3\sqrt{\pi}}\chi^{3}$ as $\chi\rightarrow 0$ and $h(\chi)\rightarrow 1-\frac{2}{\sqrt{\pi}\chi}$ as $\chi\sim 1$. For shallow holes $\psi\ll 1$, we have $U_{c}\simeq 1+\mathcal{O}(\psi)$ and for deep holes $\psi\sim 1$, $U_{c}\sim(m_{i}/m_{e})^{1/4}\psi^{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 $\gamma=\mathrm{Im}(\omega)$ of the instability when $\lvert U\rvert<U_{c}$ can be evaluated by solving the equation $\dot{P}_{i}/\dot{P}_{e}(\omega)+1=0$ numerically for given $\phi(x)$, $U$ and the mass ratio. In Figure 5, we show in dashed lines the EH velocity calculated as a function of $\psi$ for $\gamma=0.1$ and Schamel type of EH potential as a useful reference line. We shall give a detailed analysis of the growth rate in III.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 time step of $0.3/\omega_{pe}$. A discrete Fourier transform of $U$ during the first few unstable periods gives a sharp peak centered 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 Eqn. (33) calculated using the potential $\phi(x)$ right before the oscillation onset from our PIC simulations. $\phi(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 $\psi\sim 0.45$, hence $U_{c}$. The oscillation frequency predicted by our theory is inversely proportional to the hole width $L$: $\omega_{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\lambda x$, the Eqn. (26) giving the critical speed $U_{c}$ is invariant under this change of scale as each term is multiplied by the same factor $1/\lambda$. This property is obvious for the first two terms in Eqn. (26). For the third term, it can be shown that

[TABLE]

Hence the crossing point satisfies the same scaling relation $C_{\lambda}(\tilde{\phi})=(1/\lambda)C(\tilde{\phi})$ and $U_{c}$ remains invariant. However, the oscillation frequency scales linearly with $\lambda$: $\omega_{c,\lambda}=\lambda\omega_{c}$. For example, two different EH potentials $\phi(x)=\psi\operatorname{sech}^{4}(x)$ and $\phi(x)=\psi\operatorname{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.

III.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 $\phi(x)$ is symmetric, then the counter-streaming situation with an EH at rest $U=0$ and two ion streams traveling at $\pm v_{i}$ is equivalent to having one single ion stream at rest and the EH traveling 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 $\phi(x)$ is symmetric so that $\phi^{\prime}(-x)=-\phi^{\prime}(x)$ and $\phi^{\prime\prime}(-x)=\phi^{\prime\prime}(x)$, Eqns. (14) (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.

III.3 Linear growth rate

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

[TABLE]

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

[TABLE]

Substituting $\dot{P}_{i}/\dot{P}_{e}(\omega,U)=-F(\omega/U)/G(U)$, the two partial derivatives in Eqn. (37) can be evaluated with functions $F$ and $G$

[TABLE]

Hence

[TABLE]

where we used $F(-\omega_{c}/U_{c})/G(-U_{c})=1$. $G$ is an even polynomial function defined in Eqn. (23) and $-U_{c}G^{\prime}(-U_{c})\equiv U_{c}G^{\prime}(U_{c})$ can be evaluated as

[TABLE]

Hence the final expression for $\Delta\omega$ is

[TABLE]

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

[TABLE]

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

[TABLE]

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

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 $\psi\sim 0.8$, in which $\Delta\lvert U\rvert$ 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 Eqn. (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 $\gamma\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 $\gamma$ is large. We are going to discuss these nonlinear effects in Section IV.

III.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’s consider the ion thermal speed to be small compared to $\lvert U\rvert$ such that $v_{\mathrm{th,i}}\ll\lvert U\rvert$ 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}\geq T_{e}$. We can integrate the contributions from ions of different velocities to get the total $\dot{P}_{i}$

[TABLE]

We apply a Taylor series expansion to Eqn. (45) assuming $\lvert v\rvert\ll\lvert U\rvert$. 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

[TABLE]

where

[TABLE]

We have right now $\dot{P}_{i}/\dot{P}_{e}(\omega,U,\phi,v_{\mathrm{th,i}})=-F_{\mathrm{th,i}}(\omega/U)/G(U)$ to the leading order in $|v_{\mathrm{th,i}}/U|$. It suffices to substitute $F$ with $F_{\mathrm{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_{\mathrm{th,i}}/U|$. The effect of finite ion temperature can be visualized through the $F_{\mathrm{th,i}}$ contours. In Figure 8, we show the $F_{\mathrm{th,i}}$ contours evaluated on the real axis for different values of $|v_{\mathrm{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{\phi})$ 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_{\mathrm{th,i}}$, it grows to $\sim 10\%$ for $|U|=4v_{\mathrm{th,i}}$ and $\sim 20\%$ when $|U|=3v_{\mathrm{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_{\mathrm{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 [22]. 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_{\mathrm{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.

IV Discussion

The deformation of hole potential during oscillation has been neglected in our analysis. It is a next order correction for the “jetting” effect we’ve calculated. We now show that the ion momentum change due to the hole potential variation can be ignored in the parameter regime we are interested in. When an ion transits through the EH potential, its momentum changes when there is a temporal variation of the hole potential height. We call this the hole growth “jetting” effect and a formula is given in our previous paper [22]

[TABLE]

where $\partial\phi/\partial t$ represents the temporal variation of the hole potential in its rest frame. It is related to the change in charge density by Poisson’s equation. To a first approximation, we assume that the frequency of velocity oscillation is low so that the electron density remains the same in the hole frame. Therefore, the only density variation comes from the ions and we have approximately $\delta\phi/\phi\sim\delta n_{i}/n_{i}$. Using $n_{i}\partial\phi/\partial t\sim\phi\partial n_{i}/\partial t$, we get

[TABLE]

The $(v_{bf}-v_{as})_{\mathrm{accel}}$ is the “jetting” effect due to hole acceleration. It was shown in Eqn. (6). The ion momentum change due to self-consistent EH potential variation is on the order of a factor $\psi/U^{2}\ll 1$ smaller than the momentum change due to EH acceleration. Thus we can ignore it in the momentum balance.

In our analysis, the trapped electrons are assumed to move with the hole potential while remaining on their trapped orbits in an oscillating hole. However, there are always shallowly trapped electrons whose slow orbits will resonate with the oscillation frequency. The fraction of these resonant particles is small and they do not much affect the linear stability analysis. Once the instability has fully developed, some resonant particles become detrapped, causing the EH to shrink in size. This is the nonlinear stage of the instability. Trapped electrons in a steady-state EH are tagged in our simulations to follow their motion [24]. The detrapping of trapped electrons by the instability is shown in Figure 9. It is observed that the instability can be nonlinearly saturated by the shrinking of the EH. In some other cases, the oscillation amplitude in the hole velocity grows until the EH velocity is in the close vicinity of the ion velocity. The EH is then disrupted by the ions through the mechanism described by Saeki et al. [11] The presence of a parallel electric field can slow down the EHs in the ion frame and lead to the instability. Happening in space, this instability can cause ion heating by ion Landau damping and drive anomalous resistivity [20]. The EHs are considered to stem from phase-space instability and they are reservoirs of wave energy. The described instability provides a mechanism to couple the stored wave energy in the EH to the ion and electron plasma energy.

Our solitary wave velocity stability theory is generic and appears to apply to ion-acoustic solitons. Why are ion-acoustic solitary waves, propagating at a velocity slightly higher than $c_{s}$ in the ion frame, stable? $\dot{P}_{e}$ goes to [math] when $|U|$ approaches $c_{s}$ for small wave amplitude $\psi$. Hence, for a solitary wave propagating around this velocity, solutions of our dispersion relation, $\dot{P}_{i}/\dot{P}_{e}+1=0$, exist only at very high freuqency where the short-transit-time approximation for electrons break down. So it is possible to have a stable propagation region for the solitary wave at a speed around $c_{s}$. For a solitary wave in this regime, $|\dot{P}_{i}/\dot{P}_{e}(\omega)|\gg 1$ for all physical frequencies of the system. The ion dynamics therefore dominates over the electron dynamics inside the solitary wave and it propagates like a Korteweg-de Vries-type ion-acoustic soliton even though an electron phase-space structure might be attached to it.

In contrast, a solitary wave propagating much faster than $U_{c}$ is dominated by electron dynamics: $|\dot{P}_{i}/\dot{P}_{e}(\omega)|\ll 1$ for all physical $\omega$, and can be considered a pure electron hole. When the ion dynamics and electron dynamics are comparable inside a propagating solitary wave, our stability theory predicts that the propagation can be driven unstable by positive dynamical feedback between the two species. The velocity instability reported in this paper naturally separates these two major types of plasma electrostatic solitary waves. We will address this point in greater detail in follow up work.

V Conclusion

In this paper, we’ve reported a new kind of instability for an EH propagating with the presence of heavier ions and we have presented the theoretical understanding behind it. An EH at low speed in the ion frame experiences unstable velocity oscillations that can be understood treating it as a holistic object. Our analytic treatment is in full agreement with the PIC simulation observations. The “slow” EHs that are observed in space might be susceptible to this instability. We demonstrated that the velocity oscillations initially take place at the critical ion transit frequency through the hole potential and it grows quickly to a noticeable level once the speed is below the threshold. The instability happens when the electron and ion dynamics are $180$ degrees out of phase. Our discovery is a type of solitary wave instability driven by the different inertial scales of the two constituent species of the plasma. Analogies may also be found in other physical systems.

Appendix A: rate of change of contained ion momentum

We calculate the rate of change of inertial-frame ion momentum contained in the control volume from $x_{a}$ to $x_{b}$ to the leading order in $t_{ab}\dot{U}/U$

[TABLE]

$U_{\mathrm{ref}}$ is the hole velocity in the reference frame. $U_{\mathrm{ref}}$ is the characteristic of the reference frame and does not depend on $x$. We have $\dot{U}_{\mathrm{ref}}=\dot{U}$. Choose the inertial frame such that $U_{\mathrm{ref}}=0$ (the instantaneous hole rest frame). Then at $t=t_{f}$

[TABLE]

We apply the continuity of an ion fluid element from $x_{a}$ to a position $x$ between $x_{a}$ and $x_{b}$

[TABLE]

where the subscript $s_{x}$ refers to the time $t_{s_{x}}$ for an ion particle at $x$ when $t=t_{f}$ to first enter the control volume at $x_{a}$, so $t_{s_{x}}=t_{f}-\int_{x_{a}}^{x}\frac{du}{v}$. Using the constancy of the inflow density, we express $nv$ as

[TABLE]

with

[TABLE]

$t_{ax}=\int_{x_{a}}^{x}\frac{dx_{1}}{v}$ is the ion transit time between $x_{a}$ and $x$. The ion velocity is governed by the conservation of its first integral of motion

[TABLE]

We can then give an expression for $\frac{\partial t_{ax}}{\partial t_{s_{x}}}|_{x}$ to its leading order

[TABLE]

We used interchangeably $dt$ and $dt_{sx}$ as $t(x_{2},x)=t_{f}-\int_{x_{2}}^{x}\frac{du}{v}$ and $t(x)=t_{s_{x}}+\int_{x_{a}}^{x}\frac{dx_{1}}{v_{0}(x_{1})+v_{1}(x_{1},t_{s_{x}})}$ so that $dt=dt_{s_{x}}(1+\mathcal{O}(t_{ab}\dot{U}/U))$ for $x$ fixed. The first term of $\dot{P}_{ic}$ in Eq. (A-2) can be expressed as

[TABLE]

Combining Eqn. (Appendix A: rate of change of contained ion momentum) and Eqn. (A-5), we have an expression for $\frac{\delta t_{as_{x}}}{\delta t_{xf}}$

[TABLE]

Therefore, to the lowest order, the term I in Eq. (Appendix A: rate of change of contained ion momentum) is

[TABLE]

To calculate the term II in Eq. (Appendix A: rate of change of contained ion momentum), we need

[TABLE]

We use the results from Eq. (Appendix A: rate of change of contained ion momentum) to get

[TABLE]

Thus, the term II in Eq. (Appendix A: rate of change of contained ion momentum) can be written as

[TABLE]

Combining everything above, we have to the leading order

[TABLE]

Now we apply integration by parts to the last term in Eq. (Appendix A: rate of change of contained ion momentum) using $\frac{\partial t(x_{2},x)}{\partial x_{2}}|_{x}=\frac{1}{v(t(x_{2},x))}$

[TABLE]

Plug this expression in Eq. (Appendix A: rate of change of contained ion momentum) and remember that $v(t_{s_{x}})\simeq-U$, we have

[TABLE]

where

[TABLE]

and

[TABLE]

Simplify Eq. (Appendix A: rate of change of contained ion momentum) by replacing the $\ddot{U}$ terms using Eqs. (Appendix A: rate of change of contained ion momentum) (Appendix A: rate of change of contained ion momentum). Most terms cancel out to the relevant order and we are left with only one leading order term

[TABLE]

where we used the equation of motion for a single ion particle

[TABLE]

We combine $\dot{P}_{io}$ and $\dot{P}_{ic}$ to get $\dot{P}_{i}$

[TABLE]

Appendix B: Leading order expansion of $\dot{P}_{i}/\dot{P}_{e}$ in $\psi/U^{2}$

First order in $\psi/U^{2}$

The full expression of $\dot{P}_{i}$ is

[TABLE]

We will begin by expanding the expression of $\dot{P}_{i}$ to the first order in $\psi/U^{2}$. Recall that the equilibrium ion velocity in the hole frame $v_{0}$ is

[TABLE]

And the phase term inside the integral sign of Eqn. (B-1) is

[TABLE]

The term inside the first integral can be rewritten taking $v_{0}\simeq-U$ as

[TABLE]

which is already a second order term. The only first order contributions come from the last two integrals. To the first order in $\psi/U^{2}$, $\dot{P}_{i}$ can be written as

[TABLE]

We first apply integration by parts to the triple integral term in Eqn. (B-5)

[TABLE]

We interchange the order of integration to integrate the triple integral. The integration domain of this triple integral is

[TABLE]

We define an indicator function $\mathds{1}_{\mathcal{D}}$ associated with this measurable subset of $\mathds{R}^{3}$ and we have

[TABLE]

Interchanging the order of integration is the mathematical technique we use throughout this derivation to deal with multiple integrals. Since the first term in Eqn. (B-8) cancels the first term in Eqn. (B-6) and

[TABLE]

Using this result in Eqn. (B-5) and we get the exact cancellation of the two first-order terms. Therefore, we need to push the expansion to the next order, which is $\psi^{2}/U^{4}$.

Second order in $\psi/U^{2}$

First, we identify all the terms that are of order $\psi^{2}/U^{4}$ in the expansion of $\dot{P}_{i}$. We examine each term in Eqn. (B-1). We have for the first term

[TABLE]

The second term will give contributions both from $v_{0}$ and the phase term

[TABLE]

The third term can be expanded as

[TABLE]

We will deal with these integrals one by one, starting from the simpler double integrals and moving our way to the quadruple integral. The idea is to use integration by parts and interchanging the order of integration to simplify them as much as possible as we showed previously. In the end, we should have simpler integral expressions involving only $\tilde{\phi}$ and not its derivatives. Let’s start with the double integral that is multiplying the second order term in Eqn. (B-11).

[TABLE]

The triple integral in Eqn. (B-11) is

[TABLE]

Now we simplify the triple integral in Eqn. (B-12) using twice integration by parts and interchanging the order of integration

[TABLE]

Now we deal with the only term left, which is the quadruple integral in the Eqn. (B-12). We will try to relate it to the integrals we’ve already calculated. Observing the integrand, we notice that it is easy to integrate with respect to the variable $x_{1}$. We interchange the order of integration as we did before

[TABLE]

Notice that the above integral can be expressed with the derivative w.r.t. $(i\omega/U)$ of another integral that we’ve calculated previously in the beginning of Eqn. (B-15)

[TABLE]

Therefore, the quadruple integral can be simplified using the final result of Eqn. (B-15)

[TABLE]

Now we add up all the second order contributions using Eqns. (B-10) (B-13) (B-14) (B-15) (B-18). The majority of terms cancel out and we are left with

[TABLE]

The expansion of $\dot{P}_{e}$ in $\psi/U^{2}$ is trivial

[TABLE]

We introduce the constant $I_{0}$ and a new function $I$

[TABLE]

We thus get the final leading order expansion of $\dot{P}_{i}/\dot{P}_{e}$ as it appears in Eqn. (19)

[TABLE]