## 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 Reference Landau1946). 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.* Reference Roberg-Clark, Drake, Reynolds and Swisdak2018); 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.* Reference Zhang, Angelopoulos, Artemyev and Liu2018) and in solar wind (Alexandrova *et al.* Reference Alexandrova, Chen, Sorriso-Valvo, Horbury and Bale2013).

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 Reference Sagdeev and Grad1967; Galeev & Sagdeev Reference Galeev, Sagdeev, Galeev and Sudan1984*b*; LaBelle & Treumann Reference LaBelle and Treumann1988) that breaks the frozen flux condition and potentially sets the reconnection rate (e.g. Ji *et al.* Reference Ji, Yamada, Hsu and Kulsrud1998; Kulsrud Reference Kulsrud1998, Reference Kulsrud2001; Treumann Reference Treumann2001; Uzdensky Reference Uzdensky2003). Indeed, there is observational (e.g. Deng & Matsumoto Reference Deng and Matsumoto2001; Farrell *et al.* Reference Farrell, Desch, Kaiser and Goetz2002; Matsumoto *et al.* Reference Matsumoto, Deng, Kojima and Anderson2003) and numerical (e.g. Drake *et al.* Reference Drake, Swisdak, Cattell, Shay, Rogers and Zeiler2003; Zhang *et al.* Reference Zhang, Chien, Gao, Ji, Blackman, Follett, Froula, Katz, Li and Birkel2023) 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 Reference Alt and Kunz2019; Winarto & Kunz Reference Winarto and Kunz2022).

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 Reference Gurnett and Anderson1977; Gurnett & Frank Reference Gurnett and Frank1978). Recently, IAWs have been receiving progressively more attention thanks to ongoing space missions, namely NASA's Magnetospheric Multiscale Mission (MMS) (Burch *et al.* Reference Burch, Moore, Torbert and Giles2016) and Parker Solar Probe (PSP) (Fox *et al.* Reference Fox, Velli, Bale, Decker, Driesman, Howard, Kasper, Kinnison, Kusterer and Lario2016), and the European Space Agency's Solar Orbiter (Müller *et al.* Reference Müller, Marsden, St. Cyr and Gilbert2013). 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.* Reference Graham, Khotyaintsev, Vaivads, Edberg, Eriksson, Johansson, Sorriso-Valvo, Maksimovic, Souček and Píša2021; Píša *et al.* Reference Píša, Souček, Santolík, Hanzelka, Nicolaou, Maksimovic, Bale, Chust, Khotyaintsev and Krasnoselskikh2021). 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.* Reference Mozer, Bonnell, Bowen, Schumm and Vasko2020, Reference Mozer, Bonnell, Hanson, Gasque and Vasko2021*a*). In addition, IAWs have also been identified in the separatrix and outflow region on the magnetospheric side of the reconnecting magnetopause (Uchino *et al.* Reference Uchino, Kurita, Harada, Machida and Angelopoulos2017; Steinvall *et al.* Reference Steinvall, Khotyaintsev, Graham, Vaivads, André and Russell2021), 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 Reference Kovrizhnykh1966; Rudakov & Korablev Reference Rudakov and Korablev1966; Zavoiski & Rudakov Reference Zavoiski and Rudakov1967); whereas the latter is mainly studied semi-quantitatively based on the Kadomtsev–Petviashvili (KP) model considering wave–ion nonlinear interactions (Kadomtsev & Petviashvili Reference Kadomtsev and Petviashvili1962; Petviashvili Reference Petviashvili1963). 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 Reference Sagdeev and Grad1967; Sagdeev & Galeev Reference Sagdeev and Galeev1969). Although Sagdeev's resistivity formula has been widely adopted (e.g. Uzdensky Reference Uzdensky2003), 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 Reference Bychenkov and Silin1981, Reference Bychenkov and Silin1982; Bychenkov, Gradov & Silin Reference Bychenkov, Gradov and Silin1982; Bychenkov, Gradov & Silin Reference Bychenkov, Gradov and Silin1984). 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.* Reference Boris, Dawson, Orens and Roberts1970; Biskamp & Chodura Reference Biskamp and Chodura1971; Biskamp, Von Hagenow & Welter Reference Biskamp, Von Hagenow and Welter1972; DeGroot *et al.* Reference DeGroot, Barnes, Walstead and Buneman1977; Dum & Chodura Reference Dum, Chodura, Palmadesso and Papadopoulos1979; Ishihara & Hirose Reference Ishihara and Hirose1981, Reference Ishihara and Hirose1983). 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.* (Reference Tavassoli, Chapurin, Jimenez, Papahn Zadeh, Zintel, Sengupta, Couëdel, Spiteri, Shoucri and Smolyakov2021) 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 Reference Fijalkow1999; Horne & Freeman Reference Horne and Freeman2001) 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 Reference Watt, Horne and Freeman2002; Petkaki *et al.* Reference Petkaki, Watt, Horne and Freeman2003), the initial drift velocity of electrons (Petkaki & Freeman Reference Petkaki and Freeman2008) and phase space structures (Büchner & Elkina Reference Büchner and Elkina2006; Lesur, Diamond & Kosuga Reference Lesur, Diamond and Kosuga2014). Some of these simulations adopted realistic electron–ion temperature and mass ratios (e.g. Hellinger, Trávníček & Menietti Reference Hellinger, Trávníček and Menietti2004; Petkaki *et al.* Reference Petkaki, Freeman, Kirk, Watt and Horne2006). 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 Reference Galeev, Sagdeev, Galeev and Sudan1984*a*).

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 Reference Biskamp and Chodura1971). Analytical studies (e.g. Zavoiski & Rudakov Reference Zavoiski and Rudakov1967; Sagdeev & Galeev Reference Sagdeev and Galeev1969) 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 Reference Bychenkov, Silin and Uryupin1988)

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.* Reference Bychenkov, Silin and Uryupin1988). This is, therefore, a necessary first step in the development of a complete understanding of IAT.

We use Gkeyll (Shi *et al.* Reference Shi, Hammett, Stoltzfus-Dueck and Hakim2017, Reference Shi, Hammett, Stoltzfus-Dueck and Hakim2019; Juno *et al.* Reference Juno, Hakim, TenBarge, Shi and Dorland2018; Hakim & Juno Reference Hakim and Juno2020; Mandell *et al.* Reference Mandell, Hakim, Hammett and Francisquez2020), a recently developed state-of-the-art code featuring energy-conserving, high-order discretisation methods to solve the 2-D Vlasov–Poisson equations:

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.

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.

### 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

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),

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.

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,

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 2*a*), 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’):

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 2*b*), the 1-D dielectric function obtained from linear theory can be written as

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 4*a*) 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.

#### 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 Vasko2021*a*) 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 2*a*), 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*).

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.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 Reference Hamberger and Jancarik1972; Perepelkin *et al.* Reference Perepelkin, Suprunenko, Vasil'ev and Diky1973; de Kluiver, Perepelkin & Hirose Reference de Kluiver, Perepelkin and Hirose1991), 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 Reference Kadomtsev and Petviashvili1962; Petviashvili Reference Petviashvili1963). 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 2*b*). Bychenkov & Silin (Reference Bychenkov and Silin1982) 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 Reference Bychenkov and Silin1982; Bychenkov *et al.* Reference Bychenkov, Silin and Uryupin1988). 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 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.

#### 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 Reference Dum, Chodura and Biskamp1974) and can be explained by the quasi-linear theory for ions (Ishihara & Hirose Reference Ishihara and Hirose1981). 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 Reference Ishihara and Hirose1981, Reference Ishihara and Hirose1983). 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.* Reference Chan, Hara, Wang, Jain, Mirjalili and Boyd2022).

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 4*b*) 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 (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.

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. 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.* Reference Bychenkov, Silin and Uryupin1988):

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.

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),

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,

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.

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

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

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

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.

#### 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

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.

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.* Reference Polito, Dudík, Kašparová, Dzifčáková, Reeves, Testa and Chen2018) and the low-speed solar wind near the Sun (Mozer *et al.* Reference Mozer, Bale, Cattell, Halekas, Vasko, Verniero and Kellogg2022; Verscharen *et al.* Reference Verscharen, Chandran, Boella, Halekas, Innocenti, Jagarlamudi, Micera, Pierrard, Štverák and Vasko2022), 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.* Reference Mozer, Bonnell, Bowen, Schumm and Vasko2020) and interplanetary shocks (Wilson *et al.* Reference Wilson, Chen, Wang, Schwartz, Turner, Stevens, Kasper, Osmane, Caprioli and Bale2020), 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 Reference Verscharen, Klein and Maruca2019). 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.* Reference Uchino, Kurita, Harada, Machida and Angelopoulos2017; Steinvall *et al.* Reference Steinvall, Khotyaintsev, Graham, Vaivads, André and Russell2021). In addition, bursts of EAWs following IAT are experimentally observed in laser-driven reconnection events (Zhang *et al.* Reference Zhang, Chien, Gao, Ji, Blackman, Follett, Froula, Katz, Li and Birkel2023). 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.* Reference Bychenkov, Silin and Uryupin1988). These findings align with recent IAW observations in the solar wind (Mozer, Vasko & Verniero Reference Mozer, Vasko and Verniero2021*b*), where the most intense waves propagate obliquely. The broadening of the spectrum during the nonlinear stage, as reported in other studies (Mozer *et al.* Reference Mozer, Bonnell, Hanson, Gasque and Vasko2021*a*), is also consistent with the wavenumber spectrum observed in our simulations (see figure 4*a*).

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.* Reference Mozer, Vasko and Verniero2021*b*). 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 Reference Coles and Harmon1988; Huang *et al.* Reference Huang, Sioulas, Shi, Velli, Bowen, Davis, Chandran, Matteini, Kang and Shi2023; Lotz *et al.* Reference Lotz, Nel, Wicks, Roberts, Engelbrecht, Strauss, Botha, Kontar, Pitňa and Bale2023). 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.* Reference Huang, Sioulas, Shi, Velli, Bowen, Davis, Chandran, Matteini, Kang and Shi2023). 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.* Reference Verscharen, Klein and Maruca2019). 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.

## Acknowledgements

The authors thank I. Hutchinson, and N. Mandell for insightful discussions. The authors thank the Gkeyll team for their help with the simulations.

*Editor Alex Schekochihin thanks the referees for their advice in evaluating this article.*

## Funding

This work was supported by the NSF-DOE Partnership in Basic Plasma Science and Engineering Award No. PHY-2010136 (ZL and NFL) and the Partnership for Multiscale Gyrokinetic Turbulence (MGK), part of the US Department of Energy (DOE) Scientific Discovery Through Advanced Computing (SciDAC) program, via DOE contract DE-AC02-09CH11466 for the Princeton Plasma Physics Laboratory and subaward No. UTA18-000276 to MIT under US DOE Contract DE-SC0018429 (MF). This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the US Department of Energy under Contract No. DE-AC02-05CH11231 using NERSC award FES-ERCAP0026577. This research also used resources of the MIT-PSFC partition of the Engaging cluster at the MGHPCC facility, funded by DOE award No. DE-FG02-91-ER54109 and the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the US Department of Energy under Contract No. DE-AC02-05CH11231 using NERSC award FES-ERCAP0020063.

## Declaration of interests

The authors report no conflict of interest.

## Data availability statement

Readers may reproduce our results and also use Gkeyll for their applications. The code used in this study is available online. Full installation instructions for Gkeyll are provided on the Gkeyll website (https://gkeyll.readthedocs.io). The input files used in the study are under version control and can be obtained from the repository at https://github.com/ammarhakim/gkyl-paper-inp/tree/master/2023_JPP_iat.

## Appendix A. Linear theory of IAI

Starting from the linearised Vlasov–Poisson equations (2.2), the plasma dielectric function can be written as (Schekochihin Reference Schekochihin2022)

where $p = -\textrm {i}\omega + \gamma$ and $\alpha$ refers to different species. Given a wavevector $\boldsymbol {k}$, one can always choose the $z$-axis to be along $\boldsymbol {k}$. Therefore, the dielectric function can be rewritten as

where $F_\alpha (v_z)$ is the 1-D velocity distribution function, and prime denotes differentiation with respect to $v_z$. The analytical dispersion relation can be obtained by solving

Figure 12 shows plots of the growth rates ($\gamma$), growth rates divided by wavenumbers ($\gamma /k$) and phase space velocities ($\omega /k$) as functions of the wavenumber, obtained by solving (A3) for cases with different mass ratios. Both electrons and ions have a Maxwellian velocity distribution. The electron distribution has a drift velocity of $0.5 v_{\textrm {Te}0}$ with respect to ions. The electron-to-ion temperature ratio is 50, as in most of our runs. Figures 12(*b*) and 12(*c*) are useful for determining the numerical resolution necessary; see Appendix F.

## Appendix B. Comparison with 1-D simulation

We run a 1D1V simulation with the same physical parameters as run Main. Figure 13 shows the differences between the 1-D and 2-D simulations. The simulation with the label of ‘2D’ is run Main. The evolution of the 1-D system is close to the scenario described by the quasi-linear theory in that it exhibits relatively steady wave energy after saturation, as shown by the top left panel in the figure. As a result, ions can constantly gain energy from the waves at an approximately constant rate via nonlinear trapping. We indeed see ion temperature growing linearly on the top right panel. Moreover, as shown by the bottom left panel, although a sharp drop in the current happens during saturation, its growth rate quickly returns to the free-acceleration level. The corresponding anomalous resistivity, plotted in the bottom right panel, shows only a sharp peak around saturation. No shutdown of IAT nor transition into EAWs is observed in the 1-D case, which can be due to the different retained plateau shapes of the electron velocity distribution as explained in the following.

We plot 1-D electron velocity distribution functions in figure 14 right after saturation ($\omega _{\textrm {pe}}t=1000$) and at the end of phase III for the 2-D case ($\omega _{\textrm {pe}}t=1900$). In the 1-D case, after saturation, the electron distribution remains flat at the phase velocities of the waves. This plateau greatly reduces the wave emissions by electrons and keeps the wave energy approximately constant afterward (see the red and blue dashed curves in figure 14). On the other hand, the bulk of the electrons is still located on the right-hand side of the plateau and are still freely accelerated by the external electric field. Therefore, after the quick drop in current due to plateau formation, the current growth rate quickly returns to a value that is close to the free acceleration rate. In contrast, in the 2-D case, plateau formation in the parallel direction cannot shut down the growth of IAWs due to the presence of oblique modes. The saturation of these oblique modes makes the scattering much more efficient than the 1-D case, and scatters most of the electrons back to the resonance region (the region with lower $v_z$). Consequently, the 1-D electron distribution exhibits in the 2-D case a bulk–tail shape, where both bulk and tail can be reasonably well approximated with Maxwellian functions (see the red and blue solid curves in figure 14). As the bulk remains stationary after saturation, the current growth rate is significantly reduced. In sum, the qualitative and quantitative differences between the electron distribution functions obtained in the 1-D and 2-D runs result in very different nonlinear dynamics and demonstrate that the physics of IAT is not well captured by 1-D models.

## Appendix C. Electron velocity distribution fitting parameters

We summarise the fitting parameters at different moments for run Main in table 3.

The growth rates of different wave modes at different times, obtained by solving the linear dispersion relation, are plotted in figure 15. It can be seen that all the wave modes have growth rates close to zero at the moment of saturation ($\omega _{\textrm {pe}}t=750$). However, as the electron-to-ion temperature decreases, wave modes begin to be more strongly damped, especially those with larger wavenumbers. The blue dashed curve is obtained by using the fitted electron distribution at $\omega _{\textrm {pe}}t=1000$ but artificially increasing the ion temperature to its value at $\omega _{\textrm {pe}}t=1900$; this yields a similar curve to that obtained at $\omega _{\textrm {pe}}t=1900$. Therefore, we conclude that the strong damping of IAWs is indeed mainly caused by ion heating.

## Appendix D. EAW dispersion relation

To identify the wave modes in phase V, we calculate the theoretical dispersion relation of EAWs and compare it with the observed wavenumbers and frequencies in our simulation. For electron waves, the Landau dispersion relation takes the form (Landau Reference Landau1946)

where $\omega = \omega _{r} + i\gamma$ is the complex frequency, $k$ the wavenumber and $F_0$ the electron velocity distribution in the direction of wave propagation. The EAW dispersion relation is obtained by retaining only the principal part of the velocity integral in (D1), i.e.

This is because EAW theory assumes that there is a narrow plateau in the electron velocity distribution around the phase velocities of the EAWs, thereby rendering the waves immune to Landau damping (Holloway & Dorning Reference Holloway and Dorning1991). The assumption can be valid because the trapped particle distribution can flatten the velocity distribution at the wave phase velocity.

To solve (D2) analytically, we consider parallel wave modes and use the double Maxwellian model to approximate the 1-D electron distribution during phase V of run M25E10. The fitted result is shown in figure 16(*a*), and the dispersion relation obtained by solving (D2) is plotted in figure 16(*b*). These theoretically predicted wavenumbers, frequencies and phase velocities agree well with the observed properties of waves during phase V (see the $\omega$–$k$ diagram in figure 7).

## Appendix E. Extremely weak field regime

The threshold external electric field of the *extreme weak field regime* can be estimated as follows. The growth rate of IAI predicted by the linear theory is approximately

where $u(t) = eE_\textrm {ext}t/m_\textrm {e}$ is the bulk electron drift velocity (Schekochihin Reference Schekochihin2022). The wave energy density, $W$, as a function of time is then

where $W_0$ is the initial wave energy density (noise). The wave energy density at the moment of saturation, $W_\textrm {sat}$, can be estimated by quasi-linear theory; it is

The time at which saturation occurs, i.e. $t_\textrm {sat}$, can thus be estimated by solving

We can estimate the drift velocity of electrons in the tail at the moment of saturation

To be in the *extremely weak field regime*, the difference between the ion sound speed (phase velocity of the waves) and the drift velocity of the electron tail at the moment of saturation should be smaller than the resonance width so that most of the tail can be effectively trapped by the waves, i.e.

where $\Delta v_\textrm {trap}$ is the velocity width (in the parallel direction) of an electron hole in phase space (i.e. trapping width). Here $\Delta v_\textrm {trap}$ can be estimated using

where $\phi$ is the typical electric potential at the moment of saturation in the system. We can estimate $\phi$ using the saturation wave energy density, i.e.

where we take $k$ to be the wavenumber of the mode with the largest linear growth rate ($k\approx 2{\rm \pi} \times 0.12 \lambda _{\textrm {De}}^{-1}$). Therefore,

The corresponding trapping width is then

Substituting $t_\textrm {sat}$ and (E10) into (E6), we obtain the inequality

We now can find a critical $E_\textrm {ext}$, i.e. $E_\textrm {ext}^\textrm {crit}$, that makes the approximate sign in inequality (E11) hold. Numerically, we find the $E_\textrm {ext}^\textrm {crit}$ under our simulation set-up for the mass ratio of 25 is (assuming $W_0/n_0 T_\textrm {e0} \sim 10^{-8}$)

which is consistent with our simulation results. Runs M25E2 and M25E1 should indeed be in the *extremely weak field regime*.

In the *extremely weak field regime*, the dynamics in phase III exhibit some distinct features. For example, in run M25E2, as shown in figure 17(*a*,*c*), both the wave energy and current remain approximately constant after saturation. The key difference between this regime and run Main is that a large fraction of electrons, including those in the tails, become resonant with IAWs. As the tail portion cannot efficiently run away, the shape of electron velocity distribution barely changes during phase III. Consequently, wave energy will not decrease until a critical temperature ratio in reached. In this simulation, after the critical temperature ratio ($T_\textrm {e}/T_\textrm {i} \sim 10$) is reached at around $\omega _{\textrm {pe}}t = 3000$, the IAWs are quickly damped, leading to a rapid decrease in wave energy. The corresponding ion heating and anomalous resistivity, plotted in figure 17(*b*,*d*), are also quickly turned off.

## Appendix F. Determining the numerical resolution

#### F.1 Ideal resolution

The numerical resolution in real space is exclusively determined by the unstable wavenumbers of IAI. Assuming $n$ grid points are required to resolve a single wavelength, the corresponding grid size is approximately

where $k_\textrm {max}$ is the largest unstable wavenumber of IAWs, which can be predicted from linear theory. Conveniently, we find that $k_\textrm {max}$ depends weakly on the drift velocity between electrons and ions according to linear theory, which implies that the ideal resolution in real space is unchanged throughout the course of the simulations. According to the numerical scheme we are using, $n$ is around 10. We can see in figure 12(*a*) that $k_\textrm {max}\lambda _{\textrm {De}} \approx 2{\rm \pi} \times 0.36$ for all mass-ratio cases. Plugging this $k_\textrm {max}$ into (F1), we obtain

With regard to velocity space, we use two methods to guide the choice of numerical resolution. The first method considers the phase velocities of linearly unstable IAWs. We assume that at least three velocity grid points are required to be located in the interval of phase velocities of unstable IAWs $[v_{\textrm {ph},\textrm {min}}, v_{\textrm {ph},\textrm {max}}]$, i.e.

This method was adopted in Petkaki *et al.* (Reference Petkaki, Freeman, Kirk, Watt and Horne2006).

The second method estimates the width of the electron resonance region due to the IAWs when the system enters the nonlinear regime. We may assume that a mode will enter the nonlinear regime when the bouncing frequency of an electron in the potential well caused by this wave is comparable to the wave's linear growth rate $\gamma _k$ (Sagdeev Reference Sagdeev and Grad1967). If the maximum electric potential of the wave mode is $\phi$, the bouncing frequency is $\omega _B = \sqrt {ek^2\phi /2m_\textrm {e}}$. Letting $\omega _B \sim \gamma _{k}$, we obtain $\phi \sim m_\textrm {e} \gamma _{k}^2/ek^2$. The resonance width in electron velocity space, $\Delta v$, caused by this wave is $\Delta v \sim \sqrt {e\phi /m_\textrm {e}}$, and the corresponding resolution is

To resolve most of the unstable waves, estimates given by (F3) and (F4) produce similar results according to figure 12(*b*,*c*) (we truncate the wavenumber at $k\lambda _{\textrm {De}} = 0.3$ when estimating with (F4)). From figure 12(*b*,*c*), we also see that the ideal numerical resolution in velocity space depends linearly on the inverse of the square root of mass ratio. We summarise them here:

#### F.2 Linear benchmarks

We perform a set of benchmarks against linear theory to check both the simulation set-up and the numerical resolution. In these simulations, no external electric field is applied. Electrons have an initial drift velocity ($u_\textrm {e} = 0.5 v_{\textrm {Te}0}$, which is a typical drift velocity during phase I of most of our simulations in the paper) such that the initial condition is ion-acoustic unstable. An ion-to-electron mass ratio of 25 and an electron-to-ion initial temperature of 50 are used. We only perturb one wave mode in each simulation. All the simulations are run to the end of the linear stage.

The numerical resolution we use in our simulations differs from (F2) and (F5). This is because adopting (F2) demands unaffordable computational power. However, due to the fact that modes with very small or large wavenumbers are weak in amplitude and, thus, unimportant to the nonlinear evolution of the system, we only need to resolve modes within a certain range of wavenumbers that have relatively strong growth rates. Several different resolution combinations are proposed in table 4 based on the estimates given by (F2) and (F5). We refer to grid sizes in real space as ($\Delta z, \Delta y$) and grid sizes in electron and ion velocity spaces as ($\Delta v_{\textrm {e},z}, \Delta v_{\textrm {e},y}$, $\Delta v_{\textrm {i},z}$, $\Delta v_{\textrm {i},y}$). The combination #3–1 has the same resolution in the $z$ (parallel) direction with #3, and it is only used for the nonlinear test described in § F.3. The ion velocity domain is small and thus demands much fewer computational resources. We set $\Delta v_{\textrm {i},z} = \Delta v_{\textrm {i},y} = v_{Ti0}/3$ for all simulations.

The benchmark results are concluded in figure 18. Combination #1 can correctly resolve the growth rates with wavenumbers smaller than about $2{\rm \pi} \times 0.14 \lambda _{\textrm {De}}^{-1}$, while combinations #3 and #4 can resolve higher wavenumber modes up to $2{\rm \pi} \times 0.22 \lambda _{\textrm {De}}^{-1}$. Different from combination #1, combination #2 cannot resolve wave modes with higher wavenumbers because of insufficient numerical resolution in electron velocity space. The fact that combinations #3 and #4 can resolve a similar range of wave modes implies that $\Delta v_{\textrm {e},z} = 0.042 v_{\textrm {Te}0}$ is enough to resolve most of the linearly unstable wave modes.

#### F.3 Nonlinear tests

To understand how these numerical resolutions influence the nonlinear evolution of our system, we perform turbulence simulations with different numerical resolution combinations in table 4.

The set-up for these turbulent tests is the same as for the results presented in the main paper. Both electrons and ions have zero drift velocity at the start, and an external electric field is applied to drive IAI. The external electric field is chosen to be $\tilde {E}_\textrm {ext} =2.5\times 10^{-3}$, which corresponds to E10 in run names. All simulations are long enough to identify the bursts of EAWs.

The time traces for perturbed electric field energy, ion temperature, current in the parallel direction and anomalous resistivity are plotted in figure 19. It can be seen that all of these resolution combinations exhibit similar linear and nonlinear evolution behaviour.

The simulation that deviates most from the others is combination #1: it exhibits higher ion heating than the others and somewhat shorter phases III and IV. On the other hand, combination #2 behaves much better than resolution combination #1, even though it is only able to resolve a similar wavenumber range as combination #1 in the linear benchmarks. As expected, not being able to resolve the linear growth rates of some modes correctly does not imply that the nonlinear behaviour will also be substantially influenced because wave amplitude is small during the linear phase and thus requires finer grids to resolve it. This observation implies that lower resolution in velocity space has less of an effect in the nonlinear stage. In addition, the highly overlapped time traces between combination #3 and #3–1 imply that lower resolution in the $y$ direction of electron velocity space is harmless. Based on these observations, we finally conclude that resolution combination #3–1 can capture nonlinear physics accurately enough for the mass ratio of 25 cases.

To summarise, as the real space numerical resolution does not depend on mass ratio, the grid size for all the simulations can be chosen as

The required velocity space numerical resolution linearly depends on the inverse square root of mass ratio (see figure 12 and (F5)–(F7)), so the resolution of electron velocity space for the other cases can be calculated based on combination #3–1. For example, the numerical resolution of electron velocity space of run Main should be

For run M200E10, we are forced by computational resource constraints to keep this numerical resolution (F9). While insufficient to adequately capture the linearly unstable spectrum, these nonlinear tests demonstrate that phases III and IV of the evolution should still be reasonably resolved.