3.1. Phase I: linear stage

At the start of the simulation, both electrons and ions are at rest. They are freely accelerated by the external electric field, thus ramping up the current and driving the system toward instability. During this stage, the current growth rate is the free acceleration rate, i.e. $\textrm {d} J_z/\textrm {d} t = E_\textrm {ext} e^2 n_0 / m_\textrm {e}$ (which neglects the small ion contribution to the current), where $J_z$ is the parallel current in the system.

In figure 1, the red dotted vertical line indicates the moment of time, $t_0$, when the current growth rate ($1/J_z\,\textrm {d} J_z/\textrm {d} t$) is comparable to that of the most unstable IAW predicted by linear theory (Jackson Reference Jackson1960) (linear theory is briefly described in Appendix A). We observe that this criterion predicts the onset of the instability very well. According to linear theory, the growth rate of an ion-acoustic mode, $\gamma _\textrm {e}(k)$, is roughly proportional to the drift velocity of electrons (as long as the drift velocity of the electrons is much larger than the ion sound speed). Because the electrons are freely accelerated by the external electric field, the linear growth rate of an ion-acoustic mode at time $t$ is then approximately $\gamma _\textrm {e}(t, k) \approx \gamma _\textrm {e}(t_0, k)(t/t_0)$. Therefore, the wave amplitude should grow as

(3.1)\begin{equation} W(t) = W(t_0) \exp({2\gamma_{\rm e} (t,k) t}) \approx W(t_0)\exp\left({\frac{2\gamma_{\rm e}(t_0,k)}{t_0} t^2}\right). \end{equation}

We calculate $\gamma _\textrm {e}(t_0, k)$ of the most unstable mode using linear theory, and the resulting $W(t)$ is plotted with the black dashed line in figure 1. The observed wave energy agrees reasonably well with (3.1).

This linear phase ends when the wave energy becomes large enough that wave–particle interaction becomes significant and the electron velocity distribution begins to be modified, at around $\omega _{\textrm {pe}}t \approx 500$.

3.2. Phase II: saturation

Phase II captures the saturation stage of IAI. We first investigate the saturation mechanism. As mentioned in § 1, two processes are potentially responsible for saturation: quasi-linear interactions between electrons and waves, and nonlinear interactions with ions (ion trapping Sagdeev & Galeev Reference Sagdeev and Galeev1969; Biskamp & Chodura Reference Biskamp and Chodura1971). The electron and ion responses to the growth of the IAWs are represented in figure 2(a,b), respectively, by plotting the one-dimensional (1-D) velocity distribution function at different moments of time, obtained by averaging over the spatial domain ($z,y$) and $v_y$; namely, $F_{\alpha }(v_z) \propto \int f_{\alpha }(v_z,v_y, z,y)\,\textrm {d} z\,{\textrm {d}y}\,\textrm {d} v_y$. As the dark blue curve in figure 2(a) shows, a quasi-linear plateau starts forming in the resonance region in the electron velocity distribution at around $\omega _{\textrm {pe}}t\approx 500$, which weakens the growth of the IAI. This moment of time coincides with the onset of strong electron heating, as evidenced by the electron temperature time trace shown in the inset plot of figure 3. Importantly, we observe (in the same figure) that strong ion heating (which is the signature of ion trapping) only begins later in time, at around $t\approx 650 \omega _{\textrm {pe}}^{-1}$, when the growth rate of wave energy is already reduced significantly and the wave energy is about to arrive at its peak. This suggests that ion trapping is not the main saturation mechanism. We further verify this conjecture by summarising the maximum wave energy density we observe in simulations with different mass ratios and external electric fields $W_\textrm {max}$ in table 2 and comparing them with the quasi-linear estimate (Zavoiski & Rudakov Reference Zavoiski and Rudakov1967; Bychenkov et al. Reference Bychenkov, Silin and Uryupin1988),

(3.2)\begin{equation} W_{\rm{sat}} \approx n_0T_{\rm e0} \tilde{E}_{\rm{ext}}\left(\frac{m_{\rm i}}{2m_{\rm e}}\right)^{1/2}, \end{equation}

where $W_\textrm {sat}$ is the theoretical wave energy density at saturation. Our numerical results indeed show an approximately linear dependence of the peak wave energy on the external electric field $E_\textrm {ext}$ and on the square root of mass ratio $(m_\textrm {i}/m_\textrm {e})^{1/2}$. In addition, the absolute value of maximum wave energy density is only a little bit larger than (3.2). Therefore, we conclude that, under the weak external electric field condition, quasi-linear relaxation of the electron distribution, rather than ion trapping, is the main IAI saturation mechanism, consistent with the theoretical prediction (Bychenkov et al. Reference Bychenkov, Silin and Uryupin1988). This conclusion is seemingly at odds with the numerical study by Biskamp & Chodura (Reference Biskamp and Chodura1971), which argues instead that ion trapping is the mechanism responsible for saturation. We think that this discrepancy is due to the fact that their simulations start with a super-critical electron distribution. This causes the waves to grow rapidly to large amplitudes. Thus, nonlinear effects, rather than quasi-linear relaxation of electron distribution, dominate the saturation of the IAI in their simulations.

Figure 2. 1-D electron (a) and ion (b) (in logarithmic scale) velocity distribution at different times. The shadowed region is the approximate electron resonance region at the moments around saturation ($\omega _{\textrm {pe}} t \approx 750$). The black dashed line represents the fitted bulk in the electron velocity distribution after saturation, and the dashed blue and green curves represent the fitted tail at different times. The fitting model is the double Maxwellian function (3.4). See the text for more details.

Figure 3. Time traces of electron and ion temperatures, normalised to their respective initial temperature (red dashed and red dotted lines, left axis), and time trace of electron-to-ion temperature ratio (blue line, right axis). The temperatures are defined as the second moment of their own velocity distribution functions.

Table 2. Saturation wave energy of simulations with different mass ratios and external electric fields. Here $W_\textrm {max}$ is the maximum wave energy measured from the simulation and $W_\textrm {sat}$ is the quasi-linear estimate of the saturated wave energy, i.e. (3.2). The wave energy exhibits approximate linear dependence on the external electric field and on the square root of mass ratio, consistent with the quasi-linear prediction.

Despite a good agreement between theory and simulations on wave energy at the moment of saturation, the quasi-linear theory assumes that IAT will reach a steady state (i.e. approximately constant wave energy) determined by the balance between wave emission ($\gamma _\textrm {e}(k)$; weakened due to plateau formation) and wave damping ($\gamma _i(k)$), namely,

(3.3)\begin{equation} \gamma(\boldsymbol{k}) = \gamma_{\rm e}(\boldsymbol{k}) + \gamma_i(\boldsymbol{k}) = 0. \end{equation}

However, it can be seen from figure 1 that a steady state is not achieved: the wave energy quickly starts to drop after saturation. We show later in this section that (3.3) is indeed achieved at the moment of saturation, but does not remain valid afterward.

To understand the nonlinear evolution after saturation, observe the blue ($\omega _{\textrm {pe}}t=750$) and the cyan ($\omega _{\textrm {pe}}t=1000$) curves in figure 2(a). They display two noteworthy features: (i) the (1-D) electron velocity distribution does not retain the quasi-linear plateau at the end of the saturation process; (ii) the non-resonant part is mostly represented by a tail (see the blue dashed curve in figure 2a), which appears due to significant electron heating (see figure 3) and acceleration by the external electric field during saturation. The disappearance of the plateau shape in the electron velocity distribution is, in fact, attributed to the saturation of oblique modes, which is a notable 2-D feature. The much more efficient scattering of electrons by the waves in the 2-D situation brings the majority of electrons back to the resonance regime. More evidence of this can be found in Appendix B.

Based on these observations, we find that the 1-D electron velocity distribution at the end of, and after, saturation can be well approximated by a double Maxwellian function, with one Maxwellian representing the bulk of the electron population (denoted with subscript ‘b’ in the following formula) and the other representing the tail (subscript ‘t’):

(3.4)\begin{align} F_{\rm e}(v_z) & \approx (1-\alpha) \frac{1}{\sqrt{2{\rm \pi}}v_{{\rm Te,b}}} \exp \left(-\frac{(v_z-u_{d, {\rm b}})^2}{2v_{{\rm Te,b}}^2}\right) \nonumber\\ & \quad +\alpha \frac{1}{\sqrt{2{\rm \pi}}v_{{\rm Te,t}}} \exp \left(-\frac{(v_z-u_{d,{\rm t}})^2}{2v_{{\rm Te,t}}^2}\right). \end{align}

This bulk–tail model is predicated on the notion that, as the wave energy saturates, a resonant and a non-resonant part should be manifest in the electron velocity distribution. The acceleration of the resonant part by the external electric field should be countered by the anomalous resistivity created by IAT, while the non-resonant part can still be freely accelerated. The fitted distribution at the moment of saturation ($\omega _{\textrm {pe}}t \approx 750$) is plotted with a black dashed curve (the bulk) and a dashed blue curve (the tail) in figure 2(a). The exact fitting parameters can be found in Appendix C.

Using the empirical fit, we now demonstrate that (3.3) still holds at the end of the saturation even though the electron distribution no longer conforms to the plateau shape described by the standard quasi-linear theory (a similar method was previously employed in studies of Buneman instability in reconnection layers Drake et al. Reference Drake, Swisdak, Cattell, Shay, Rogers and Zeiler2003; Che et al. Reference Che, Drake, Swisdak and Yoon2009, Reference Che, Drake, Swisdak and Yoon2010; Jain, Umeda & Yoon Reference Jain, Umeda and Yoon2011). Because ions do not deviate significantly from a Maxwellian distribution around saturation (see figure 2b), the 1-D dielectric function obtained from linear theory can be written as

(3.5)\begin{equation} 1+ \frac{1+\zeta_{i} Z(\zeta_{i})}{k_z^{2} \lambda_{D i}^{2}} +(1-\alpha)\frac{1+\zeta_{\rm e,b} Z(\zeta_{\rm e,b})}{k_z^{2} \lambda_{\rm D e,b}^{2}} + \alpha \frac{1+\zeta_{\rm e,t} Z(\zeta_{\rm e,t})}{k_z^{2} \lambda_{\rm D e,t}^{2}} = 0, \end{equation}

where $\lambda _{\textrm {D}\alpha } = v_{T\alpha }/\omega _{\textrm {p}\alpha }$ is the Debye length of species $\alpha$ (ions, bulk electrons or tail electrons), $\zeta _{\alpha } \equiv (\omega - ku_{\textrm {d}\alpha }) / kv_{T\alpha }$, with $\omega = \omega _{r} + i \gamma _\textrm {e}$, $u_{\textrm {d}\alpha }$ the drift velocity and $Z(\zeta )$ the plasma dispersion function. Using the fitting model provided by (3.4) evaluated at $\omega _{\textrm {pe}} t = 750$, we obtain from this equation $\gamma (k_z,\omega _{\textrm {pe}}t = 750) \approx 0$ for nearly all values of $k_z$ that are originally unstable. While not strictly rigorous, because we only consider waves in the parallel direction in (3.5), this result strongly suggests that the saturation process of IAI indeed tends to bring the system to a marginally stable state as described by (3.3). We note that, as mentioned earlier, (3.3) is not valid at later times, which voids the assumption made by the quasi-linear theory; see § 3.3.

To summarise, during the saturation process in the 2-D case, electrons are efficiently scattered back to lower velocities by both parallel and oblique wave modes within the resonance region. Consequently, the electrons within the resonance region become the dominant population, forming a new electron bulk. After saturation, a balance is established between the effective friction provided by the waves and the acceleration from the external electric field, leading to the stationary behaviour of the new bulk, while electrons located on the right-hand side of the resonance region experience continuous acceleration.

3.3. Phase III: shutdown of IAT

Phase III is a post-saturation stage characterised by a significantly lower current growth rate (about a factor of three lower than the free acceleration rate); see figure 1. As explained in § 3.2, the growth of the parallel current after saturation is primarily contributed by this fitted tail. We can estimate the current growth rate in the early stages of phase III to be approximately $\alpha (t = 1000) e^2E_{\textrm {ext}}/m_\textrm {e}$, where $\alpha (t = 1000)$ represents the fraction of the fitted tail at $t = 1000$. This approximation yields a value close to $0.38 e^2E_{\textrm {ext}}/m_\textrm {e}$, which is only slightly higher than the current growth rate indicated by the green dashed line in figure 1. During phase III, the wave energy gradually decreases and reaches a level approximately one order of magnitude smaller than its peak value. Eventually, the intensity of the IAT becomes insufficient to provide enough friction to the bulk electrons.

3.3.1. Landau damping induced by strong ion heating

It has long been conjectured that the IAI will be switched off due to the reduction in the electron-to-ion temperature ratio. This reasoning is based on the observation, from linear theory, that when $J_z/(e n)\gg c_s$ the ion-acoustic mode becomes stable for $T_\textrm {e}/T_\textrm {i}\lesssim 10$ (Papadopoulos Reference Papadopoulos1977; Benz Reference Benz2012). We confirm this prediction by plotting the temperature ratio, in figure 3. Even though both electrons and ions are strongly heated after saturation, we observe an abrupt decrease in the electron-to-ion temperature ratio, which happens mostly in phase II and early phase III. This plummeting of the temperature ratio can be directly responsible for a strong damping rate of IAWs for most wave modes.

However, linear theory assumes Maxwellian electron and ion velocity distributions, which is not the case in the nonlinear evolution of IAI at this stage. To prove this conjecture more rigorously, we return to the linear analysis of the modified distribution function to show that IAWs will indeed be Landau damped after saturation. By plugging the fitting parameters at $\omega _{\textrm {pe}}t=1000$ into (3.5) and assuming Maxwellian ions, we observe $\gamma (k_z \geq k_c,\omega _{\textrm {pe}}t=1000) < - 3 \times 10^{-3} \omega _{\textrm {pe}}$, where $k_c \lambda _{\textrm {De}}/2{\rm \pi} \approx 0.15$ is the smallest wavenumber to satisfy this condition ($3 \times 10^{-3} \omega _{\textrm {pe}}$ is the inverse of about one-third time duration of phase III). This result implies that we should see about an order of magnitude decrease in wave energy in $\sim 300\omega _{\textrm {pe}}^{-1}$ for modes with wavenumber satisfying $k_z \geq k_c$ during phase III. By comparing the wavenumber ($k$) spectrum between $\omega _{\textrm {pe}}t=1000$ (blue curve) and $\omega _{\textrm {pe}}t=1600$ (cyan curve) plotted in figure 4(a), we indeed see a significant decrease in amplitude for the wave modes with larger wavenumbers ($k\lambda _{\textrm {De}}/2{\rm \pi} \gtrsim 0.2$). However, these wave modes are only about an order of magnitude weaker than those in $\omega _{\textrm {pe}}t = 1000$, which are less damped than predicted. Using the fitting parameters at $\omega _{\textrm {pe}}t=1600$ in (3.5), we obtain $\gamma (k_z,\omega _{\textrm {pe}}t=1600) < 0$ for nearly all values of $k_z$, while the damping rates of wave modes with $k_z \lambda _{\textrm {De}}/2{\rm \pi} > 0.08$ are strong (again, compared to $- 3\times 10^{-3}\omega _{\textrm {pe}}$). Even though most wave modes are less damped than predicted, the fact that nearly all the wave modes at $\omega _{\textrm {pe}}t=2000$ (the green curve in figure 4a) show significantly lower energy confirms our analysis. The exact fitting parameters of the bulk–tail double Maxwellian model and more details on this linear analysis can be found in Appendix C.

Figure 4. Wavenumber (a) and angular (b) spectra of IAWs at different times during the simulation. Angular spectra at different times are normalised to their value at $\theta = 0^{\circ }$, where $\theta$ is the angle between the parallel direction and wave propagating direction. A fourth-order polynomial is used to fit the angular spectra to smooth the curves. The wavenumber spectrum at $\omega _{\textrm {pe}}t=1000$ agrees with the KP spectrum (the dashed curve). The wavenumber spectra at $\omega _{\textrm {pe}}t=2000$ and $\omega _{\textrm {pe}}t=3500$ with relatively larger wavenumbers agree with the spectrum observed in Chapman et al. (Reference Chapman, Brunner, Banks, Berger, Cohen and Williams2014) (the dash-dotted curves).

3.3.2. Effects of particle trapping

An important mechanism affecting the wave energy spectrum not considered by linear or weak-turbulence theory is particle trapping, which can strongly weaken Landau damping. The formation of electron phase-space vortices (also referred to as ‘phase-space holes’ or ‘electrostatic solitary waves’ Hutchinson Reference Hutchinson2017), which is the signature of electron trapping, can be identified right after saturation (a snapshot of electron phase space at $\omega _{\textrm {pe}}t = 2000$ is shown in figure 6). This phenomenon is expected by nonlinear theories, and it is also reported in many simulations of two-stream instability (e.g. Morse & Nielson Reference Morse and Nielson1969; Berk, Nielsen & Roberts Reference Berk, Nielsen and Roberts1970), space observations where electrostatic streaming waves are present (e.g. Pickett et al. Reference Pickett, Chen, Mutel, Christopher, Santolı, Lakhina, Singh, Reddy, Gurnett and Tsurutani2008; Malaspina et al. Reference Malaspina, Newman, Willson III, Goetz, Kellogg and Kerstin2013; Mozer et al. Reference Mozer, Bonnell, Hanson, Gasque and Vasko2021a) and laboratory experiments (e.g. Saeki et al. Reference Saeki, Michelsen, Pécseli and Rasmussen1979).

One signature of electron trapping is that it can trigger positive frequency shifts of nonlinear IAWs, which may allow the sub-harmonic decay of IAWs (Berger et al. Reference Berger, Brunner, Chapman, Divol, Still and Valeo2013; Chapman et al. Reference Chapman, Berger, Brunner and Williams2013, Reference Chapman, Brunner, Banks, Berger, Cohen and Williams2014). Modes with positive frequency shifts can spread across the region in $\omega$–$k$ space approximately bounded by $\omega /k = c_{s0} + \Delta v_\textrm {i}$, where $\Delta v_\textrm {i}$ corresponds to a region where ion Landau damping is reduced by ion trapping. To obtain evidence of this, we examine the wave energy spectrum by performing 2-D Fourier transforms (in the parallel direction ($z$) and in time) of the parallel electric field $E_z$ between different time ranges. Note that due to limited storage and computational power, we are not able to save the files often enough to have high resolution in frequency for run Main. Instead, the wave energy spectra $|E_z(k,\omega )|^2$ are calculated using run M25E10 at the corresponding evolution phases. In the linear phase, the $\omega$–$k$ diagram is consistent with the IAI linear dispersion relation, as shown in figure 5(a). In phases II and III, as shown in figure 5(b), we indeed see the positive frequency shift of IAWs in a wide range of wavenumbers. The maximum shift of the phase speed, $\Delta v_\textrm {i}$, can be estimated using $\Delta v_\textrm {e}$, which is the half-width of the plateau of the electron velocity distribution before saturation (the plateau for run Main can be seen in figure 2a), i.e. $\Delta v_\textrm {i} \approx \sqrt {m_\textrm {e}/m_\textrm {i}} \Delta v_\textrm {e} \approx 0.2 v_{\textrm {Te}0}$. This estimate roughly agrees with the observed spectrum in figure 5(b).

Figure 5. Parallel wave energy spectrum $|E_z(\omega, k_z)|^2$ during the IAW burst phase (i.e. phase I) (a), and the phase right after saturation (i.e. phases II and III) (b) for run M25E10. The dispersion relation of IAWs (A3) is plotted with the red dashed curve in (a,b). The blue dashed line in (b) indicates a set of waves with positive frequency shifts in the nonlinear stage.

As mentioned above, this frequency shift is expected to induce sub-harmonic decay of IAWs and, consequently, modulate the energy spectrum. In Chapman et al. (Reference Chapman, Brunner, Banks, Berger, Cohen and Williams2014), the onset of IAT leads to $\phi _k \propto k^{-4}$ at relatively large wavenumbers. The corresponding wave energy spectrum is then $N(k) \propto |E_k|^2 \propto k^2 \phi _k^2 \propto k^{-6}$. This spectrum is observed in our simulation at the late stage of phase III and in phase IV, as evidenced by the dash-dotted curves in figure 4(a).

3.4. Phase IV: post-shutdown phase

At the end of phase III, the current growth rate returns to a value that is close to the free acceleration rate, as shown by the orange dashed line in figure 1. The ion heating rate also relaxes to a relatively low level (only about 25 % ion heating happens after phase III, see figure 3). After IAT is shut down, the system still undergoes a long evolution process before entering a new regime. We refer to this nonlinear evolution stage as phase IV.

To gain a better understanding of this stage, we plot the 2-D features of the system in figure 6. Snapshots of the 2-D electron velocity distribution are shown in the second column. At $\omega _{\textrm {pe}} t = 2000$, the bulk of the electron distribution is centred around the resonance region. The subsequent shutdown of the IAI decreases the effective friction on the electrons, and allows the bulk to migrate to the tail region by $\omega _{\textrm {pe}} t = 3500$. By the end of phase IV, the electron velocity distribution is greatly broadened in the parallel direction with respect to what it was at the end of phase III. Furthermore, we observe that the tail of the electron distribution exhibits a triangular shape, consistent with early analytical conjectures (Sagdeev & Galeev Reference Sagdeev and Galeev1969). This is because oblique unstable modes have a maximum propagation angle with respect to the parallel direction. The electron distribution population that lies outside of the resonance region created by the wave modes is more easily accelerated by the external electric field.

Figure 6. Snapshots of electron phase-space (a,e,i), electron velocity distribution (b,f,j), ion velocity distribution (c,g,k) and electric potential (d,h,l) at the end of phase III, the end of phase IV and the middle of phase V, respectively. The colour bar for the ion velocity distribution is in logarithmic scale. The 2-D distributions are obtained by averaging the four-dimensional (2D2V) distribution over the other two dimensions.

3.4.2. Formation of double layer

Finally, we look at the electron phase space structures in the first column of figure 6 and the corresponding electric potential in the last column. In § 3.3, we mentioned the emergence of electron holes right after saturation in § 3.3.1. During the later course of the nonlinear evolution, these phase-space vortices gradually merge and finally become one single large electron hole at the end of phase IV (which can be seen from the electron phase space plot at $\omega _{\textrm {pe}}t = 3500$). In the meantime, the potential wells created by the waves become asymmetric (which means there is a non-vanishing electric potential change across a single potential well), with a steeper edge. The steep asymmetric potential wells favour the triggering of new instabilities because they asymmetrically reflect low-energy electrons and accelerate high-energy electrons, which enhances the two-stream structure in the electron distribution and depletes the particles in the middle of the well. In the electric potential depicted in the second row of figure 6 at $\omega _{\textrm {pe}}t=3500$, we see clearly larger and more asymmetric potential structures compared with the electric potential at $\omega _{\textrm {pe}}t=2000$. The effects on the electron velocity distribution function can be seen in figure 2(a), where the depletion of the electrons at around $v_{\textrm {Te}0}$ is already visible at $\omega _{\textrm {pe}}t =3500$. Another possible consequence of particle trapping is modulational instability. In Rose (Reference Rose2005), it was demonstrated that negative nonlinear frequency shift and wave diffraction of Langmuir waves can trigger a self-focusing effect. It is possible that a similar effect could extend to the IAWs. In figure 7(a), we indeed observe a stronger negative frequency shift during phase IV, which emphasises the effects of nonlinear ion trapping at the late stage of the nonlinear evolution (Berger et al. Reference Berger, Brunner, Chapman, Divol, Still and Valeo2013). The modulated wavefronts (the localised maxima in electric potential) observed in the last column of figure 6 at $\omega _{\textrm {pe}}t=3500$, may result from this self-focusing effect.

Figure 7. Parallel wave energy spectrum $|E_z(\omega, k_z)|^2$ after effective shutdown of IAT (i.e. phase IV) (a), and the EAW burst phase (i.e. phase V) (b) for run M25E10. The dispersion relations of IAWs and EAWs are plotted with the red dashed curve in (a,b), respectively. The blue dash-dotted line in (a) indicates a set of waves with negative frequency shifts. The blue dashed line in (b) indicates a set of waves with a phase velocity of $1.2v_{\textrm {Te}0}$.

Extending from phase IV to phase V, the potential wells coalesce and steepen further and, finally, a double layer forms, as is clearly visible in the rightmost column of figure 6 at $\omega _{\textrm {pe}}t = 4300$. The formation of double layers associated with IAI has been shown via experiments (e.g. Okuda & Ashour-Abdalla Reference Okuda and Ashour-Abdalla1982; Chanteur et al. Reference Chanteur, Adam, Pellat and Volokhitin1983), simulations (e.g. Sato & Okuda Reference Sato and Okuda1980; Chanteur Reference Chanteur, Matsumoto and Sato1985) and spatial measurements (e.g. Ergun et al. Reference Ergun, Su, Andersson, Carlson, McFadden, Mozer, Newman, Goldman and Strangeway2001, Reference Ergun, Andersson, Tao, Angelopoulos, Bonnell, McFadden, Larson, Eriksson, Johansson and Cully2009). We here confirm the production of the double layer by IAT with a 2D2V simulation with a more self-consistent set-up.

3.5. Phase V: onset of electron-acoustic waves

As evidenced by the wave energy time trace in figure 1, another instability is triggered at this moment. We interpret these new wave modes as electron-acoustic waves (EAWs) (Holloway & Dorning Reference Holloway and Dorning1991; Valentini, O'Neil & Dubin Reference Valentini, O'Neil and Dubin2006; Anderegg et al. Reference Anderegg, Driscoll, Dubin, O'Neil and Valentini2009; Valentini et al. Reference Valentini, Perrone, Califano, Pegoraro, Veltri, Morrison and O'Neil2012) for the reasons given in the following.

First, we examine the wave energy spectrum with run M25E10 again, which is shown in figure 7(b). We can see that the new wave modes in phase V exhibit much higher frequencies extending from $0.1 \omega _{\textrm {pe}}$ to around $0.6 \omega _{\textrm {pe}}$ and a much faster phase velocity (${\sim }1.2v_{\textrm {Te}0}$), which excludes the possibility of (slow) ion modes. Instead, the intermediate frequency range between ion and electron plasma frequencies is standard for EAWs (Holloway & Dorning Reference Holloway and Dorning1991). Second, as seen from the electron phase structure at $\omega _{\textrm {pe}} t = 4300$ in the first column of figure 6, the electron distribution has a trapped population around $v_{\textrm {Te}0}$ with a trapping width extending to $\sim 4 v_{\textrm {Te}0}$. These trapped electrons around $v_{\textrm {Te}0}$ effectively create a plateau shape around the phase velocity of the waves (i.e. ${\sim }v_{\textrm {Te}0}$) in the electron velocity distribution (which can be seen from the second column in figure 6 at $\omega _{\textrm {pe}}t=4300$) and allow the EAWs to stay immune to Landau damping (Holloway & Dorning Reference Holloway and Dorning1991). Third, we plot the theoretically obtained dispersion relation of EAWs by approximating the electron velocity distribution with two Maxwellian functions and a plateau in figure 7(b), which agrees reasonably well with the wavenumbers and frequencies we observe in the simulation. These EAWs have lower wavenumbers compared to IAWs, which is consistent with the domain-size electron-hole observed in the last row in figure 6. In Appendix D, we describe in detail how we obtain the dispersion relation of EAWs. The formation of double layers and bursts of EAWs after bursts of IAWs are also reported in recent experimental (Zhang et al. Reference Zhang, Chien, Gao, Ji, Blackman, Follett, Froula, Katz, Li and Birkel2023) and numerical (Hara & Treece Reference Hara and Treece2019; Vazsonyi, Hara & Boyd Reference Vazsonyi, Hara and Boyd2020; Chen et al. Reference Chen, Khrabrov, Kaganovich and Li2022; Zhang et al. Reference Zhang, Chien, Gao, Ji, Blackman, Follett, Froula, Katz, Li and Birkel2023) studies.

4.1. Anomalous resistivity intensity

We define the effective collision frequency, $\nu _\textrm {eff}$, from the equation $\textrm {d} J_z/\textrm {d} t = e^2n_0/m_\textrm {e} E_\textrm {ext} - \nu _\textrm {eff}J_z$. A reference value is the quasi-linear result (Bychenkov et al. Reference Bychenkov, Silin and Uryupin1988):

(4.1)\begin{equation} \nu_{\rm{eff}}^{\rm{QL}} \approx 0.47 \omega_{{\rm pe}} \tilde{E}_{\rm{ext}}\left(\frac{m_{\rm i}}{m_{\rm e}}\right)^{1/2}. \end{equation}

Time traces of $\nu _\textrm {eff}$ for runs with the same external electric field but different mass ratios are shown in figure 8(a). It can be seen that the observed $\nu _\textrm {eff}$ in phase III is only about $10\,\% \sim 20\,\%$ of the value predicted by (4.1) for run Main. The maximum $\nu _\textrm {eff}$, which is observed during saturation (in phase II), is also smaller than the quasi-linear value (about 38 % of (4.1)). In terms of the mass ratio dependence, although the peak values depend weakly on the mass ratio, $\nu _\textrm {eff}$ after saturation shows a roughly linear dependence on $\sqrt {m_\textrm {i}/m_\textrm {e}}$. For example, run Main (green curve) exhibits approximately twice the amplitude of run M25E10 (blue curve) after saturation.

Figure 8. (a) Time traces of effective collision frequency for simulations with different mass ratios and the same external electric field ($\tilde {E}_\textrm {ext} = 2.5\times 10^{-3}$). For run Main we also plot, as a reference, the curve corresponding to $0.5 W/n_0T_\textrm {e}$, where $W$ is the wave energy density (dotted curve), and 10 % of the quasi-linear value (4.1) (dashed horizontal line). (b) Time traces of effective collision frequency normalised to the quasi-linear prediction (4.1), for simulations with different external electric fields and the same mass ratio ($m_\textrm {i}/m_\textrm {e} = 25$).

The main disagreement between the simulation results and the quasi-linear theory is that the latter assumes the electron–ion drift velocity will be fixed at around the ion sound speed after IAI saturates (Rudakov & Korablev Reference Rudakov and Korablev1966; Bychenkov et al. Reference Bychenkov, Silin and Uryupin1988). This assumption is only valid for the resonant part of the electrons. However, both resonant and non-resonant parts of electrons are accounted for in the current. Due to the fact that most of the tail, with a significantly higher drift velocity, does not resonate with the waves and is freely accelerated by the external electric field, the $\nu _\textrm {eff}$ we obtain in our simulations is much smaller than the theoretical one and continues to decrease as the current grows, even if the friction force on the bulk electrons provided by the scattering of waves is unchanged.

Figure 8(b) shows the time traces of $\nu _\textrm {eff}/\nu _\textrm {ff}^\textrm {QL}$ for simulations with the same mass ratio ($m_\textrm {i}/m_\textrm {e}=25$) but different external electric fields. As long as the electric field is relatively strong (but still in the weak field regime, i.e. $E_\textrm {ext} \ll E_\textrm {NL}$), the linear dependence of $\nu _\textrm {eff}$ on the external electric field holds. For example, $\nu _\textrm {eff}$ in runs M25E10, M25E8 and M25E6 have nearly the same normalised amplitudes. As the electric field becomes weaker, the effective collision frequency starts to approach the quasi-linear value during saturation gradually. For run M25E1, the effective collision frequency roughly agrees with the quasi-linear prediction around saturation.

Weak enough external electric fields such that (4.1) approximately holds will be referred to as extremely weak. In this extremely weak field regime, electrons have a drift velocity very close to the ion sound speed (i.e. the phase velocity of the IAW) at the moment of saturation and, therefore, most of the electron population, including the tail, can be trapped in the resonance region after saturation, which validates the assumptions underlying quasi-linear theory. Because the tail portion cannot effectively run away in this case, the electron distribution barely changes during phase III. The wave energy also remains approximately steady, pinned at the saturation level. Otherwise, most physics discussed in § 3 still holds, including the saturation and shutdown mechanisms, and transition into bursts of EAWs. An estimate of the threshold electric field and more details on the extremely weak field regime can be found in Appendix E.

Another well-known formula for anomalous resistivity, due to Sagdeev (Sagdeev Reference Sagdeev and Grad1967; Zavoiski & Rudakov Reference Zavoiski and Rudakov1967),

(4.2)\begin{equation} \nu_{\rm{eff}}^{\rm{Sagdeev}} \approx 0.2\omega_{{\rm pe}}\left(\frac{T_{\rm e}}{T_{\rm i}}\right)^{1/2} \left( \frac{T_{\rm e0}}{2T_{\rm e}} \tilde{E}_{\rm{ext}}^2\right )^{1/4}, \end{equation}

predicts an effective collision frequency of $\sim 0.06\omega _{\textrm {pe}}$ for run Main, which is over an order of magnitude larger than observed. Such discrepancy is unsurprising because the simulations are performed in the weak electric field regime, whereas Sagdeev's formula is only expected to apply when $E_\textrm {ext}>E_\textrm {NL}$.

Finally, to compare our simulations with previous numerical studies (LaBelle & Treumann Reference LaBelle and Treumann1988; Watt et al. Reference Watt, Horne and Freeman2002; Hellinger et al. Reference Hellinger, Trávníček and Menietti2004), it is convenient to calculate a more qualitative estimate, namely,

(4.3)\begin{equation} \nu_{\rm{eff}} \sim \omega_{{\rm pe}} \frac{W}{n_0T_{\rm e}}, \end{equation}

with $W$ being the wave energy density. This estimate is obtained by balancing the electron momentum loss due to the emission of waves and momentum gain from the external electric field. In fact, Sagdeev's formula also originates from this estimate, but it severely overestimates the wave energy by assuming that the nonlinear scattering of ions is the only saturation mechanism. Equation (4.3) for run Main is plotted as the black dotted curve in figure 8(a), and turns out to be a reasonable match, especially after saturation (phases III and IV). This conclusion is at odds with values of anomalous resistivity exceeding this estimate by at least one order of magnitude reported in many previous numerical studies (e.g. Watt et al. Reference Watt, Horne and Freeman2002; Petkaki et al. Reference Petkaki, Watt, Horne and Freeman2003, Reference Petkaki, Freeman, Kirk, Watt and Horne2006; Hellinger et al. Reference Hellinger, Trávníček and Menietti2004). Again, the fact that most previous simulations start with a super-critical configuration can lead to artificially high wave energy, which might be responsible for the large discrepancy between previous numerical experiments and (4.3).

To summarise, no nonlinear theory that we are aware of gives a reliable estimate of the anomalous resistivity created by IAT in the weak electric field regime. When the electric field is extremely weak, the quasi-linear value (4.1), is a good approximation. Otherwise, the anomalous resistivity after saturation is about 10–$20\,\%$ of (4.1) in our simulations. We also find that (4.3) provides a reasonable fit. However, even though we have validated (3.2) for the wave energy at saturation, the fact that $W$ does not remain at that level, but rather immediately starts to decay, precludes a priori usage of this estimate.

4.2. Particle heating and anomalous resistivity duration

As explained in § 3.3, anomalous resistivity will eventually become negligible as a result of the effective shutdown of the IAT. Therefore, in addition to characterising the intensity of the anomalous resistivity, it is of interest to understand how long it is sustained at appreciable levels as a function of mass ratio, external electric field and initial temperature ratio. Because the shutdown of IAT is mainly associated with temperatures, we investigate how electron and ion heating depends on these parameters.

Before we proceed with our analysis, we stress again that we define temperature as the second moment of the distribution function. This deviates from temperature in the normal sense, because neither the electron's nor the ion's distribution functions are Maxwellian when the system enters the nonlinear stage. This makes it challenging to derive an analytical expression that precisely describes the temporal evolution of electron and ion heating and the shutdown criterion of IAT.

4.2.1. Ion and electron heating

The time traces of electron-to-ion temperature ratio for simulations with different mass ratios, external electric fields, and initial temperature ratios are shown in figure 9. It can be seen that all our simulations reach a final temperature ratio of ${\sim }10$. Moreover, as shown by T20 (pink curve) and T100 (silver curve) in figure 9, the fact that the final temperature ratio (compared with the blue curve) is independent of the initial temperature ratio confirms the conjecture that IAT is indeed shut down by reaching a specific temperature ratio.

Figure 9. Time traces of temperature ratios for simulations with different mass ratios, electric fields and initial temperature ratios. The black dashed line marks $T_\textrm {e}/T_\textrm {i} = 10$.

For use in what follows, let us define the final temperature ratio as

(4.4)\begin{equation} \frac{T_{\rm e,f}}{T_{\rm i,f}} = \theta, \end{equation}

where $T_\textrm {i,f}$ and $T_\textrm {e,f}$ denote the final ion and electron temperature, i.e. the temperature when the ion heating is effectively shut down. The effective shutdown is defined as the moment when the ion heating rate reduces to 10 % of its peak value. Note that $T_\textrm {i,f}$ is a little bit larger than the ion temperature at the end of phase III because there is still some ion heating happening during phase IV. Figure 10 plots $T_\textrm {i,f}$ for simulations with different mass ratios and external electric fields, showing a linear dependence on $\sqrt {m_\textrm {i}/m_\textrm {e}}$ and $E_\textrm {ext}$. This is the same as the dependence of $W_\textrm {sat}$ on $E_\textrm {ext}$ and $\sqrt {m_\textrm {i}/m_\textrm {e}}$ in (3.2) and confirmed in table 2. We thus fit the ion temperatures in most of our simulations as

(4.5)\begin{equation} n_0 T_{\rm i,f} \approx 28.6 W_{\rm sat}. \end{equation}

According to (4.4) and (3.2), the electron temperature at the end of phase III can then be fitted by

(4.6)\begin{align} n_0T_{\rm e,f} & \approx 28.6 \theta W_{\rm sat} \nonumber\\ & \approx 20.2 \theta\left(\frac{m_{\rm i}}{m_{\rm e}}\right)^{1/2} \tilde{E}_{\rm{ext}} n_0T_{\rm e0}. \end{align}

However, there are two simulations that do not conform to (4.5) in figure 10: runs M25E1 and M25E2 (the leftmost two blue rectangular data points). As discussed in § 4.1, these two simulations are under the extremely weak field regime. In fact, the final ion temperature $T_\textrm {i,f}$ given by (4.5) can potentially become insufficiently low to satisfy the shutdown criterion, as described in (4.4). This scenario arises particularly when the external electric field is extremely weak. This implies the existence of another threshold for the external electric field below $E_\textrm {NL}$, which we call the extremely weak field regime. Under this regime, as mentioned previously, wave energy will be sustained at a high level during phase III. As the ion heating rate is proportional to the wave energy, being in this regime results in more intense ion heating than the prediction given by (4.5) until the shutdown criterion (4.4) is finally met.

Figure 10. Final ion temperature (normalised to its initial value) as a function of the square root of mass ratio at fixed $\tilde {E}_\textrm {ext}=2.5 \times 10^{-3}$ (top horizontal axis, red) and as a function of the external electric field at fixed $m_\textrm {i}/m_\textrm {e}= 25$ (bottom horizontal axis, blue).

4.2.2. Duration of significant anomalous resistivity

Finally, to estimate the duration of phase III, we plot time traces of electron heating rates for simulations with different mass ratios and external electric fields in figure 11. Both the time-dependent patterns and amplitudes of the electron heating rates are very similar to the effective collision frequencies shown in figure 8(a), especially during phase II and the early stage of phase III, which is expected because electron heating is caused by anomalous resistivity. However, unlike the anomalous resistivity, the electron heating rate does not die away but remains flat at later times in phase III. This is due to the definition of temperature that we have adopted: the increase in electron temperature is, in fact, mainly caused by the acceleration of the electron tail during phase III (so the overall electron distribution becomes broader). To provide a rough quantitative estimate of the duration of significance anomalous resistivity, we may approximate the electron heating rate with a constant value according to figure 11 as $(1/T_\textrm {e0})\,\textrm {d} T_\textrm {e}/\textrm {d} t \approx 0.08 \tilde {E}_\textrm {ext} (m_\textrm {i}/m_\textrm {e})^{1/2} \omega _{\textrm {pe}}$, which has the same dependence on external electric field and mass ratio as the (4.1). Then we can use $\tau _\textrm {res} \approx (T_\textrm {e,f} - T_\textrm {e0}) / (\textrm {d} T_\textrm {e}/\textrm {d} t)$ and substitute $T_\textrm {e,f}$ with (4.6). We arrive at

(4.7)\begin{equation} \tau_{\rm{res}} \approx \frac{20.2 \theta (m_{\rm i}/m_{\rm e})^{1/2} \tilde{E}_{\rm{ext}} -1 }{0.08 (m_{\rm i}/m_{\rm e})^{1/2}\tilde{E}_{\rm{ext}}} \omega_{{\rm pe}}^{{-}1}. \end{equation}

While the exact value of $\theta$ depends on parameters, our simulations suggest that $\theta \approx 10$ (it is slightly larger for cases with a larger mass ratio and a smaller external electric field, see figure 9). This is consistent with linear theory (Papadopoulos Reference Papadopoulos1977; Benz Reference Benz2012), and is justified by the fact that IAI becomes stable when $T_\textrm {e}/T_\textrm {i}<10$ (as long as the drift velocity between electrons and ions is much smaller than electron thermal velocity). However, it is important to note that the linear theory is based on Maxwellian distribution functions, while in our simulations both the electron and ion distributions are rather better approximated by double Maxwellians. Therefore, using temperature as the sole parameter cannot capture the exact dynamics and the corresponding switch-off criterion accurately. This may explain the deviations of $\theta$ from the value $10$. Note that (4.7) does not apply to extremely weak field regime due to the underestimate of electron and ion heating.

Figure 11. Time traces of electron heating rate for some example simulations with different mass ratios and electric fields.

With (4.7), we are able to estimate the extent of the time interval over which the IAT-induced anomalous resistivity will last in a particular system. For our run Main, (4.7) yields $\tau _{res} \sim 2000 \omega _{\textrm {pe}}^{-1}$ if we adopt $\theta = 10$, which is in reasonable agreement with the numerical data shown in figure 8(a).