A two-dimensional numerical study of ion-acoustic turbulence

Abstract

We investigate the linear and nonlinear evolution of the current-driven ion-acoustic instability in a collisionless plasma via two-dimensional (2-D) Vlasov–Poisson numerical simulations. We initialise the system in a stable state and gradually drive it towards instability with an imposed, weak external electric field, thus avoiding physically unrealisable super-critical initial conditions. A comprehensive analysis of the nonlinear evolution of ion-acoustic turbulence (IAT) is presented, including the detailed characteristics of the evolution of the particles’ distribution functions, (2-D) wave spectrum and the resulting anomalous resistivity. Our findings reveal the dominance of 2-D quasi-linear effects around saturation, with nonlinear effects, such as particle trapping and nonlinear frequency shifts, becoming pronounced during the later stages of the system's nonlinear evolution. Remarkably, the Kadomtsev–Petviashvili (KP) spectrum is observed immediately after the saturation of the instability. Another crucial and noteworthy result is that no steady saturated nonlinear state is ever reached: strong ion heating suppresses the instability, which implies that the anomalous resistivity associated with IAT is transient and short-lived, challenging earlier theoretical results. Towards the conclusion of the simulation, electron-acoustic waves are triggered by the formation of a double layer and strong modifications to the particle distribution induced by IAT.

1. Introduction

A defining property of weakly collisional plasmas is their ability to support resonant energy transfer between waves and particles. The best-known manifestation of this phenomenon is the famous Landau damping, whereby waves are resonantly damped via energy transfer to the plasma particles (Landau 1946). A contrasting possibility is the case where energy flows in the opposite direction, from the particles to the waves, thereby growing these to nonlinear amplitudes and triggering a wide variety of complex, and often poorly understood, nonlinear plasma behaviour. Such kinetic instabilities are fundamental and ubiquitous plasma processes, thought to regulate, or at least significantly contribute to, particle behaviour in weakly collisional plasmas, thereby determining their large-scale properties. For example, it is conjectured that the whistler instability can greatly reduce the heat flux in weakly magnetised collisionless plasmas (Roberg-Clark et al. 2018); and mirror and firehose instabilities are thought to regulate the ion distribution and, thus, limit the temperature anisotropy in the Earth's magnetotail (Zhang et al. 2018) and in solar wind (Alexandrova et al. 2013).

Kinetic instabilities may also be crucial in magnetic reconnection events in collisionless plasmas, because the intense wave–particle interactions that they trigger may create an anomalous resistivity (e.g. Sagdeev 1967; Galeev & Sagdeev 1984b; LaBelle & Treumann 1988) that breaks the frozen flux condition and potentially sets the reconnection rate (e.g. Ji et al. 1998; Kulsrud 1998, 2001; Treumann 2001; Uzdensky 2003). Indeed, there is observational (e.g. Deng & Matsumoto 2001; Farrell et al. 2002; Matsumoto et al. 2003) and numerical (e.g. Drake et al. 2003; Zhang et al. 2023) evidence that verifies the existence of turbulent phenomena and wave emissions around diffusion regions of reconnection sites. In addition, kinetic instabilities may influence the reconnection onset (i.e. the sudden transition from a relatively quiescent state to the reconnection stage proper) by disturbing the current sheet formation process and ultimately setting the properties of the reconnecting current sheet (e.g. Alt & Kunz 2019; Winarto & Kunz 2022).

As one of the most prominent examples of streaming-type kinetic instabilities, ion-acoustic instability (IAI), which arises when the drift velocity between electrons and (relatively cold) ions exceeds a threshold of the order of the ion sound speed, is a strong candidate for explaining many observations in various plasma environments. Dating from as far back as the 1970s, observations by the Helios I and II spacecraft revealed the presence of ion-acoustic waves (IAWs) at heliocentric distances between 0.3 and 1 AU (Gurnett & Anderson 1977; Gurnett & Frank 1978). Recently, IAWs have been receiving progressively more attention thanks to ongoing space missions, namely NASA's Magnetospheric Multiscale Mission (MMS) (Burch et al. 2016) and Parker Solar Probe (PSP) (Fox et al. 2016), and the European Space Agency's Solar Orbiter (Müller et al. 2013). For example, recent data from the Time Domain Sampler receiver on Solar Orbiter identifies the IAW as a dominant wave mode in the near-Sun solar wind below the local electron plasma frequency (Graham et al. 2021; Píša et al. 2021). Nonlinear structures in the particles’ distribution functions associated with nonlinear IAWs, such as ion holes and electron holes, were also recently reported (Mozer et al. 2020, 2021a). In addition, IAWs have also been identified in the separatrix and outflow region on the magnetospheric side of the reconnecting magnetopause (Uchino et al. 2017; Steinvall et al. 2021), which further suggests that IAI may indeed be one essential component of magnetic reconnection in collisionless environments. Understanding the nonlinear evolution of IAI, or ion-acoustic turbulence (IAT), is thus central to the interpretation of these observations.

Attempts to understand IAT began as early as the 1950s. Early analytical works either considered electron scattering by turbulence pulsations or induced scattering of waves by ions as the mechanisms underlying the saturation of the instability, and driving anomalous resistivity. The former is commonly described using a quasi-linear approach that assumes a weak turbulence level and unperturbed ions (Kovrizhnykh 1966; Rudakov & Korablev 1966; Zavoiski & Rudakov 1967); whereas the latter is mainly studied semi-quantitatively based on the Kadomtsev–Petviashvili (KP) model considering wave–ion nonlinear interactions (Kadomtsev & Petviashvili 1962; Petviashvili 1963). In 1969, Sagdeev derived an expression for IAT-induced anomalous resistivity by calculating the nonlinear wavenumber spectrum resulting (exclusively) from the nonlinear scattering of waves by ions (Sagdeev 1967; Sagdeev & Galeev 1969). Although Sagdeev's resistivity formula has been widely adopted (e.g. Uzdensky 2003), it is important to note that it is a partially qualitative result, because it fails to consider the angular distribution of IAT. In addition, its general validity is unclear because the scattering of waves by ions is not the only mechanism determining the nonlinear evolution of IAI: changes in the particles’ velocity distribution as the waves grow and saturate can, in principle, be just as important and are not considered in Sagdeev's derivation. A comprehensive quantitative nonlinear theory of IAT was not developed until the 1980s when Bychenkov and colleagues solved the kinetic equations analytically. Their work simultaneously accounted for the quasi-linear interaction of electrons with waves and the induced scattering of waves by ions, enabling the quantitative analysis of anomalous turbulent transport (Bychenkov & Silin 1981, 1982; Bychenkov, Gradov & Silin 1982; Bychenkov, Gradov & Silin 1984). By establishing the spectral and angular distribution of ion-acoustic turbulent pulsations, they verified the correctness of Sagdeev's resistivity formula in cases where the scattering of waves by ions dominates the saturation mechanisms.

Despite these achievements, significant uncertainties remain. First, many assumptions or approximations are made in the above-mentioned theories, often excluding possibilities that may be crucial in determining the dynamics of the system. For example, nearly all analytical models assume that a nonlinear steady state will be reached, which is not necessarily true in a realistic situation (and indeed is not the case in this study). Second, no theory properly considers how the modification of the ion velocity distribution due to the heating of resonant ions changes the IAT. Thus, numerical simulations are needed to identify the key nonlinear mechanisms and either validate the current theoretical understanding or, instead, guide the development of a new, improved understanding.

Numerical studies of IAT started in the 1970s, with most simulations performed using the particle-in-cell (PIC) method (e.g. Boris et al. 1970; Biskamp & Chodura 1971; Biskamp, Von Hagenow & Welter 1972; DeGroot et al. 1977; Dum & Chodura 1979; Ishihara & Hirose 1981, 1983). Most of these studies focused on identifying the saturation mechanism of IAT and quantifying the resulting turbulent heating. However, rather than gradually driving the system toward an unstable state, which is the scenario that better conforms to the theoretical approaches to this problem, these simulations start with super-critical initial conditions (their initial condition is such that the electron–ion drift velocity significantly exceeds the threshold for IAI); in some cases, the electron drift velocity is even kept constant throughout the simulation. As a result, they are inadequate to directly verify the analytical theories. Moreover, the relatively small numbers of macro-particles employed in those PIC simulations (due to computational constraints) produced unresolved results. Indeed, even with today's computational resources, simulating resonant kinetic instabilities with the PIC method can be challenging (see, e.g., Tavassoli et al. (2021) for an example concerning the Buneman instability).

The availability of efficient continuum Vlasov–Poisson solvers which arose at the turn of the century (e.g. Fijalkow 1999; Horne & Freeman 2001) enabled a resurgence of numerical studies of IAT. These new continuum Vlasov-based simulations systematically studied the relationship between anomalous resistivity and different simulation parameters such as different types of initial distribution functions (Watt, Horne & Freeman 2002; Petkaki et al. 2003), the initial drift velocity of electrons (Petkaki & Freeman 2008) and phase space structures (Büchner & Elkina 2006; Lesur, Diamond & Kosuga 2014). Some of these simulations adopted realistic electron–ion temperature and mass ratios (e.g. Hellinger, Trávníček & Menietti 2004; Petkaki et al. 2006). One conclusion from most of these numerical studies was that the anomalous resistivity created by IAT was at least an order of magnitude larger and more strongly dependent on the drift velocity of electrons than the theoretical estimation (Galeev & Sagdeev 1984a).

However, despite the valuable insights provided by these simulations, some potential problems persist. One is related to the starting configurations used in the simulations, which are often super-critical. It remains unclear whether such initial configurations are achievable in realistic scenarios, as their realisation depends on the relative rates of current development and instability growth. Another concern is that, due to computational constraints, these simulations are restricted to one dimension in both real and velocity space (1D1V), whereas the importance of two-dimensional (2-D) effects has been identified in the early numerical work (Biskamp & Chodura 1971). Analytical studies (e.g. Zavoiski & Rudakov 1967; Sagdeev & Galeev 1969) also considered the influence of oblique modes on the instability dynamics. Perhaps as a result of this artificial restriction, the nonlinear evolution observed in these simulations is mostly highly stochastic. Consequently, none of these studies show a direct comparison against the theoretical spectra of IAT pulsations and Sagdeev's formula for anomalous resistivity.

These limitations mean that considerable uncertainty remains about the nonlinear properties of the IAI; in particular, fundamental questions such as its saturation mechanism, the resulting particle heating, and the magnitude of the anomalous resistivity that it drives remain unanswered. In this paper, we conduct a comprehensive numerical study of the nonlinear evolution of the current-driven IAI aimed at providing detailed answers to these questions. We numerically solve the Vlasov–Poisson equations in two dimensions both in real and in velocity space. To ensure the physical realisability of our system, our simulations start from a stable initial condition which is driven gradually towards instability via an imposed external electric field. This set-up allows us to present the most complete and detailed understanding of IAT to date.

This paper is organised as follows. Numerical details are provided in § 2. The results from our main simulation are presented and discussed in § 3. In § 4, we analyse the sensitivity of our results to the amplitude of the external electric field, the electron–ion mass ratio, and the initial temperature ratio. Lastly, conclusions from our work and a discussion of the implications of our results in broader contexts are presented in § 5.

2. Numerical set-up

We investigate IAI in a uniform, spatially 2-D plasma with periodic boundary conditions (thus, effectively mimicking an infinitely large 2-D plasma). All simulations start with uniform Maxwellian electrons and ions, with zero drift velocity (i.e. no net current). We focus on the classic IAI and, therefore, consider only the case when electrons are initially much hotter than the ions ($T_{\rm e0}/T_{i0} \gg 1$). An external electric field $E_{\rm {ext}}$ is applied in the $\textit{z}$ (called parallel) direction throughout the simulation to drive the current and, thus, the instability. The reference value against which $E_{\rm {ext}}$ must be compared is (Bychenkov, Silin & Uryupin 1988)

(2.1)\begin{equation} E_{\rm{NL}} = \frac{m_{\rm e} c_{s0} \omega_{{\rm pi}} T_{\rm e0}}{6{\rm \pi} e T_{i0}}, \end{equation}

where $c_{s0} = \sqrt {T_{\rm e0}/m_{\rm i}}$ is the ion sound speed, $m_{\rm e}$ is the electron mass, $\omega _{{\rm pi}}$ is the ion plasma frequency, and $\lambda _{{\rm De}}$ and $\lambda _{{\rm Di}}$ are the electron and ion Debye lengths. In this work, we focus on the so-called weak field case, $E_{\rm {ext}} \ll E_{\rm {NL}}$. This choice is guided by simplicity considerations: in this regime, wave–wave interactions are predicted to be sub-dominant, and quasi-linear theory is expected to apply (Bychenkov et al. 1988). This is, therefore, a necessary first step in the development of a complete understanding of IAT.

We use Gkeyll (Shi et al. 2017, 2019; Juno et al. 2018; Hakim & Juno 2020; Mandell et al. 2020), a recently developed state-of-the-art code featuring energy-conserving, high-order discretisation methods to solve the 2-D Vlasov–Poisson equations:

(2.2)\begin{equation} \left.\begin{gathered} \frac{\partial f_{\alpha}}{\partial t}+\boldsymbol{v} \boldsymbol{\cdot}\boldsymbol{\nabla} f_{\alpha}+ \left(\boldsymbol{E_{\rm{ext}}} - \frac{q_{\alpha}}{m_{\alpha}} \boldsymbol{\nabla} \varphi\right) \boldsymbol{\cdot}\frac{\partial f_{\alpha}}{\partial \boldsymbol{v}}=0, \\ \nabla^{2} \varphi={-}4 {\rm \pi}\sum_{\alpha} q_{\alpha} \int \boldsymbol{d}^{3} v f_{\alpha}. \end{gathered}\right\} \end{equation}

We perform a set of simulations varying the magnitude of the external electric field (restricted to the weak electric field regime), the mass ratio and the initial temperature ratio. For convenience, we introduce the normalised external electric field $\tilde {E}_\textrm {ext} \equiv E_\textrm {ext}/ (4{\rm \pi} e n_0 \lambda _{\textrm {De}})$, where $\lambda _{\textrm {De}} = v_{\textrm {Te}0}/\omega _{\textrm {pe}}$ is the electron Debye length corresponding to the initial electron temperature. Table 1 lists all simulations and the corresponding values of $\tilde {E}_\textrm {ext}$, $E_\textrm {NL}$, mass ratio and $T_\textrm {e0}/T_{i0}$. Run Main is our fiducial simulation, with a mass ratio $m_\textrm {i}/m_\textrm {e}=100$, temperature ratio $T_\textrm {e0}/T_{i0}=50$ and an external electric field of $\tilde {E}_\textrm {ext} = 2.5\times 10^{-3}$. Most other runs vary the mass ratio or electric field at fixed $T_\textrm {e0}/T_{i0}=50$. In the run name, the number that follows the letter M refers to the mass ratio, whereas the number following E refers to the multiple of $2.5\times 10^{-4}$ for $\tilde {E}_\textrm {ext}$. The last two runs listed in table 1 are meant to investigate the effect of the initial temperature ratio.

Table 1. Summary of the key parameters of the simulations. In the run name, the number following the letter M refers to the mass ratio. The number following the letter E refers to the multiple of $2.5\times 10^{-4}$ for $\tilde {E}_\textrm {ext}$.

For all simulations, the spatial domain size is $50.0 \lambda _{\textrm {De}}$ in the $\textit {z}$ direction and $25.0 \lambda _{\textrm {De}}$ in the $\textit {y}$ direction, with grid resolution $\Delta z = \Delta y = 0.5\lambda _{\textrm {De}}$. The resolution of electron and ion velocity space grids are jointly determined by linear theory, linear benchmarking and a series of nonlinear evolution tests and, therefore, vary among the different simulations. We note, however, that simulation M200E10 is performed at numerical resolutions that cannot resolve the linear stage of ion-acoustic modes with small wavelengths (but it can resolve the most unstable linear mode), due to limited computational resources. However, according to our nonlinear evolution tests, this shortcoming has only a very limited effect on the evolution of the current and heating in the nonlinear stage. A detailed description of how we determine spatial and velocity space resolution is provided in Appendix F.

In § 3, we use our fiducial simulation, run Main, to demonstrate our main results. The other simulations are employed in § 4 to discuss how the results depend on the parameters listed in table 1 and to extrapolate our results to the realistic mass-ratio case.

3. IAT driven by a weak electric field

In this section, we report on run Main, for which we plot in figure 1 the time traces of electron current in the $\textit {z}$ direction and total wave energy. Based on this figure, we identify five distinct phases of evolution, which we analyse in detail in the following subsections.

Figure 1. Time trace of total wave energy (red curve, right axis) and parallel current (blue curve, left axis) of run Main. We split the evolution of the system into five different stages; see the text for details. The red dotted vertical line indicates the moment when the current growth rate is comparable to that of the most unstable IAW predicted by linear theory. The black dashed line is $W(t)$ predicted by linear theory, see (3.1). The green dashed line indicates the current growth rate during phase III (which is 30 % of the free acceleration rate), and the orange dashed line indicates the current growth rate during phase IV (which is 85 % of the free acceleration rate).

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 1960) (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 1969; Biskamp & Chodura 1971). 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 1967; Bychenkov et al. 1988),

(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. 1988). This conclusion is seemingly at odds with the numerical study by Biskamp & Chodura (1971), 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. 2003; Che et al. 2009, 2010; Jain, Umeda & Yoon 2011). 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 1977; Benz 2012). 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. (2014) (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 2017), 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 1969; Berk, Nielsen & Roberts 1970), space observations where electrostatic streaming waves are present (e.g. Pickett et al. 2008; Malaspina et al. 2013; Mozer et al. 2021a) and laboratory experiments (e.g. Saeki et al. 1979).

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. 2013; Chapman et al. 2013, 2014). 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. (2014), 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.3.3. KP spectrum

Another noteworthy nonlinear feature observed in the wavenumber spectrum, as shown by the dashed curve in figure 4(a), is its agreement with the KP spectrum, characterised by $N(k) \sim 1/k^3 \ln (1/k\lambda _{\textrm {De}})$, immediately after the saturation stage ($\omega _{\textrm {pe}}t=1000$). While this spectrum has been observed in certain turbulent heating experiments (Hamberger & Jancarik 1972; Perepelkin et al. 1973; de Kluiver, Perepelkin & Hirose 1991), we believe this to be its first direct numerical validation.

The agreement between our simulation results and the KP spectrum may initially seem puzzling, as the original derivation of the KP spectrum necessitates strong ion nonlinear effects and negligible ion Landau damping (Kadomtsev & Petviashvili 1962; Petviashvili 1963). Although we confirmed earlier in § 3.3.1 that Landau damping is indeed negligible right after saturation, the requirement for strong ion nonlinear effects still seems to contradict the quasi-linear saturation mechanism that we confirmed previously. In fact, even though quasi-linear effects of electron distribution dominate around saturation, ion nonlinear effects (induced scattering by ions) can still become noticeable afterwards, as evidenced by the emerging plateau in the ion distribution function around $\omega _{\textrm {pe}}t = 1000$ (see figure 2b). Bychenkov & Silin (1982) also demonstrated that the KP spectrum can be obtained by simultaneously considering the quasi-linear relaxation of the electron distribution, leading to a weakening of $\gamma _\textrm {e}$, and the induced scattering process, leading to a damping rate of $\gamma _\textrm {NL}$. Specifically, they showed that the condition $\gamma (\boldsymbol {k}) = \gamma _\textrm {e} (\boldsymbol {k}) + \gamma _i (\boldsymbol {k}) + \gamma _\textrm {NL} (\boldsymbol {k}) = 0$ can lead to the emergence of the KP spectrum, even when $|\gamma _\textrm {NL}(\boldsymbol {k})|\ll \min \{|\gamma _{e}(\boldsymbol {k})|,|\gamma _{i}(\boldsymbol {k})|\}$ (Bychenkov & Silin 1982; Bychenkov et al. 1988). Our observation confirms the validity of their theory during the times around saturation and emphasises the significance of nonlinear effects in the subsequent order of weak turbulence theory in shaping the behaviour of IAT, even if they are of relatively minor importance in the saturation itself.

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 1969). 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.1. Formation of high-energy ion tail and dominance of oblique modes

Another interesting feature is that a much more obvious high energy ion-tail forms, which becomes visible in the third column of figure 6 at $\omega _{\textrm {pe}}t=2000$ and can be clearly seen at $\omega _{\textrm {pe}}t=3500$ (note that the colourmap of the ion velocity distribution plots is in logarithmic scale to better emphasise its features). By comparing the cyan, green and yellow curves in figure 2(b), we can observe that the cold bulk ions experience some heating during phase II and the early stages of phase III. However, the formation of the high-temperature tail primarily occurs during phase III and phase IV. This phenomenon is consistent with a previous 2-D simulation (Dum, Chodura & Biskamp 1974) and can be explained by the quasi-linear theory for ions (Ishihara & Hirose 1981). Initially, as the ions are extremely cold, only a very small fraction of waves with low phase velocities ($v_{\textrm {ph}}< c_s$) can directly heat the cold bulk ions through weak ion trapping. However, as time progresses, the quasi-linear effect gradually diffuses the non-resonant portion of the ion distribution towards the phase-velocity range of the IAWs, allowing for stronger direct resonant interactions with the waves. This quasi-linear diffusion of non-resonant ions is the primary mechanism contributing to the formation of the high-energy ion tail (Ishihara & Hirose 1981, 1983). It is interesting to note that similar ion distribution is also observed in a recent 2D2V Vlasov simulation of current-driven instabilities (Chan et al. 2022).

The ion tail formation can also help explain the fact that the most intense wave modes propagate with a non-zero angle with respect to the parallel direction during phase IV, as evidenced by the yellow curve in figure 4(b). According to linear theory, parallel modes are expected to dominate over oblique modes because of their higher linear growth rates and saturation levels (because parallel modes have higher effective drift velocity). This is confirmed by our simulations in the linear stage (phase I, the dark blue curve in figure 4b) and early times during the nonlinear evolution (phase III, the blue curve). However, as the quasi-linear diffusion takes place, more ions are diffused in the parallel direction from the bulk (non-resonant portion) to the resonant region, which is confirmed by the third column of figure 6 at $\omega _{\textrm {pe}}t=2000$ and $\omega _{\textrm {pe}}t=3500$. As a result, modes with small propagation angles are more heavily damped because there are more ions participating in the resonant interactions with the waves. Indeed, as shown by the yellow curve in figure 4(b), we see oblique modes dominating over parallel modes at $\omega _{\textrm {pe}}t=3500$. This observation also clearly demonstrates the relevance of oblique modes in the nonlinear evolution of IAI, even in the weak electric field regime that we consider here. At the later time of phase IV, the damping of the waves is gradually reduced by the formation of the plateau in the high-energy ion tail visible in figure 2(b), which is confirmed by the red curve in figure 1.

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 (2005), 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. 2013). 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 1982; Chanteur et al. 1983), simulations (e.g. Sato & Okuda 1980; Chanteur 1985) and spatial measurements (e.g. Ergun et al. 2001, 2009). 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 1991; Valentini, O'Neil & Dubin 2006; Anderegg et al. 2009; Valentini et al. 2012) 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 1991). 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 1991). 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. 2023) and numerical (Hara & Treece 2019; Vazsonyi, Hara & Boyd 2020; Chen et al. 2022; Zhang et al. 2023) studies.

4. Dependence of results on simulation parameters

In this section, we investigate the dependence of the results discussed previously on mass ratio, external electric field and initial temperature ratio.

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. 1988):

(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 1966; Bychenkov et al. 1988). 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 1967; Zavoiski & Rudakov 1967),

(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 1988; Watt et al. 2002; Hellinger et al. 2004), 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. 2002; Petkaki et al. 2003, 2006; Hellinger et al. 2004). 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 1977; Benz 2012), 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).

5. Conclusion and discussion

In this study, we have presented a comprehensive first principle (i.e. continuum Vlasov–Poisson) numerical investigation of the current-driven IAI in unmagnetised plasmas. By starting from a stable plasma equilibrium and gradually driving it towards instability using a (weak) external electric field, we have ensured the self-consistent evolution of the system towards instability and the physical realisability of the unstable states that are reached.

We have identified five distinct phases in the evolution of IAI and focus on the nonlinear phases. After the linear stage (phase I), the onset of which occurs when the linear growth rate of the instability approximately matches the current growth rate, we observe (phase II) that the saturation of IAI primarily occurs through the quasi-linear relaxation of the electron velocity distribution. We also observe that the 2-D (in position space) system is much more efficient than its 1-D counterpart in scattering electrons to the resonance region with lower drift velocity, which helps to form a resonant bulk and a tail in the electron velocity distribution. Moving into phase III, in contrast to the assumption made in many nonlinear or quasi-linear theories that a steady state is reached, we find that the current continues to grow, albeit at a significantly reduced rate. This ongoing growth is attributed to the non-resonant electrons, which can still be accelerated by the external electric field. In addition, the wave energy decreases due to the intense ion heating, ultimately leading to the suppression of IAT and related phenomena, including anomalous resistivity and ion heating. We also observe the KP spectrum around saturation, which had not been reproduced numerically ever before and emphasises the role of the nonlinear effect in setting the wavenumber spectrum, even though the effect itself is not the main driver of the saturation. The effect of electron trapping and sub-harmonic decay of IAWs is also identified to modulate the wave spectrum. In phase IV, we observe that the current growth rate approximately returns to the rate of free electron acceleration, due to the absence of significant wave-induced friction. Furthermore, we note the emergence of a high-energy ion tail and the prevalence of oblique modes during this phase. Finally, at the conclusion of our simulation, we observe the formation of a double layer in the electric potential, which triggers the generation of EAWs.

In the second part of our study, we have shifted our focus to the investigation of anomalous resistivity resulting from IAT and its dependence on various simulation parameters. We have studied the dependence of electron and ion heating on those parameters and derive an estimate for the time duration of significant anomalous resistivity. We have found that the quasi-linear approach yields the correct dependence on the mass ratio and external electric field; however, it generally predicts a larger absolute value of resistivity compared to our simulations. This discrepancy contradicts the findings of many previous 1-D numerical studies, where stronger anomalous resistivity is typically observed. We have also observed, and derived analytically, a new electric field threshold below which the runaway electron population is small and the quasi-linear prediction of anomalous resistivity holds.

Our study of the classical version of the IAT can be directly applicable to some physical systems. For instance, IAI can exist in the separatrix region of reconnecting current sheets, where the magnetic field may be negligible (in the no-guide-field case). Situations where electrons are heated to high temperatures while ions remain cold, such as stellar flares (Polito et al. 2018) and the low-speed solar wind near the Sun (Mozer et al. 2022; Verscharen et al. 2022), have been observed. The findings in this paper are therefore directly relevant to such environments. To be more specific, our simulations reveal features that match with direct observations in various plasma environments. For example, the plateau-shaped ion velocity distribution with a high-energy tail which we find in our simulations is qualitatively similar to those measured in the solar wind (Mozer et al. 2020) and interplanetary shocks (Wilson et al. 2020), where IAWs are believed to be present. We also demonstrate the IAT's capability to significantly heat the electron core and potentially create the ‘strahl’ in the electron distribution commonly observed in the solar wind (Verscharen, Klein & Maruca 2019). The development of anisotropy in particle velocity distribution is qualitatively similar to the measured electron distribution on the dayside of reconnecting magnetopause during IAW bursts (Uchino et al. 2017; Steinvall et al. 2021). In addition, bursts of EAWs following IAT are experimentally observed in laser-driven reconnection events (Zhang et al. 2023). We also observe that oblique modes dominate over parallel modes in the later stages of nonlinear evolution. When the driving field is strong, the dominance of oblique modes is expected to occur even earlier (Bychenkov et al. 1988). These findings align with recent IAW observations in the solar wind (Mozer, Vasko & Verniero 2021b), where the most intense waves propagate obliquely. The broadening of the spectrum during the nonlinear stage, as reported in other studies (Mozer et al. 2021a), is also consistent with the wavenumber spectrum observed in our simulations (see figure 4a).

Our results can also be adapted to study magnetic reconnection influenced by IAT. The electric field associated with reconnection is $E_\textrm {rec}\sim \gamma V_\textrm {A} B_0$, where $\gamma$ is the reconnection rate and $V_\textrm {A}$ is the Alfvén speed computed with the reconnecting field $B_0$. Situations where this electric field is smaller than $E_\textrm {NL}$ (2.1) can exist in space plasmas, such as locations of the solar wind near the Sun (at distances around 20 solar radii) where $T_\textrm {e}/T_\textrm {i} \sim 5$ and $T_\textrm {e} \sim 50\,\textrm {eV}$ (Mozer et al. 2021b). For these parameters, we find $E_\textrm {rec}/E_\textrm {NL} \sim 1.2{\rm \pi} (1/\beta _\textrm {e}) (v_\textrm {Te}/c) \sqrt {m_\textrm {i}/m_\textrm {e}} (T_\textrm {i}/T_\textrm {e}) < 1$. In this case, assuming a fast reconnection rate of $\gamma \approx 0.1$, the reconnecting electric field is about $E_\textrm {rec}/(4{\rm \pi} e n_0\lambda _{\textrm {De}}) \approx 2\gamma \beta _\textrm {e} \sqrt {m_\textrm {e}/m_\textrm {i}} (v_\textrm {Te}/c) \approx 5\times 10^{-5}$, which should be in the extremely weak field regime that we have identified in this paper. We may estimate the upper limit of the duration of significant resistivity by assuming that ion heating caused by IAT is an order of magnitude larger than the value given by (4.5) (ion heating in run M25E2 is about twice as intense as the prediction given by (4.5)). Then, according to (4.7), the corresponding $\tau _\textrm {res}$ is at most of the order of $10^3\omega _{\textrm {pi}}^{-1}$. On the other hand, the typical reconnection time, $\tau _\textrm {rec}=\gamma ^{-1} L/V_\textrm {A}$, where $L$ is the characteristic length scale of the reconnecting field, can also be estimated in this context. There is a range of potential values of $L$ that one may consider; a reasonable lower bound appears to be the scale at which the turbulent inertial range starts: around $10^3\,\textrm {km}$ at distances of this order (Coles & Harmon 1988; Huang et al. 2023; Lotz et al. 2023). At such a scale, the corresponding $V_\textrm {A}$ can be estimated using the Alfvén velocity at the outer scale and the observed magnetic power spectrum; a value of $V_\textrm {A}\approx 100\,\textrm {km}\,\textrm {s}^{-1}$ seems reasonable (Huang et al. 2023). Assuming a fast reconnection rate of $\gamma \approx 0.1$, $\tau _\textrm {rec}$ can easily extend beyond $10^{5} \omega _{\textrm {pi}}^{-1}$. Therefore, we conclude that there is ample time for IAT to induce particle heating and change the particle distribution during a reconnection event in this system (and note that this conclusion should remain valid even if our estimates of $L$ and $V_\textrm {A}$ are somewhat off).

Finally, we discuss how our choice of simulation domain size and boundary conditions may limit the validity of our results. The employment of periodic boundary conditions, which effectively imitates an infinitely extended plasma, raises the question of whether the simulated system size corresponds to realistic plasma scales. In run Main, we can estimate the travel distance of an electron by integrating the area under the blue curve in figure 1, which is approximately $10^{4} \lambda _{\textrm {De}}$. Hence, our conclusions should remain applicable in a real plasma as long as the system size is larger than this estimate. In fact, this length scale ($10^4\lambda _{\textrm {De}}$) is small compared with the system size in various plasmas. Considering the solar wind as an example, the Debye length at 1 AU is on the order of $10\,\textrm {m}$ and, so, $10^{4} \lambda _{\textrm {De}}$ corresponds to a distance ${\sim }100\,\textrm {km}$, which is still smaller than the ion skin depth (${\sim }140\,\textrm {km}$) (Verscharen et al. 2019). On the other hand, the limitation imposed by the small size of the simulation box, especially in the parallel direction, excludes the effects of long-wavelength wave modes that become prominent in phase V. The box-size double layer that we observe can result from this artificial constraint. Despite this limitation, we believe that the simulations conducted up to phase V and the mechanism elucidated for triggering the second instability remain plausible. Expanding the simulation size in the parallel direction for a 2D2V simulation, while desirable, was deemed impractical due to the large computational resources it would require.

In summary, our numerical investigation represents a significant stride toward a self-consistent understanding of the nonlinear dynamics of IAT. While we concentrate on the weak field regime in this step, which contains simpler (though nonetheless complex) nonlinear physics than the strong field case, our study yields substantial insights applicable to various physical scenarios where IAWs are present. Our study lays a robust foundation for future investigations into more complex IAT situations. These encompass strong external fields, the presence of magnetic fields, and less extreme initial temperature ratios. It is also important to acknowledge that our simulations do not encompass the hypothesised ultimate state of the system, wherein the current finally ceases to increase. The mechanism responsible for stopping particle acceleration when a collisionless plasma is exposed to an external electric field thus remains unknown. To better understand this problem, longer and more computationally demanding simulations would be required, which are currently beyond our available resources.