A two-dimensional numerical study of ion-acoustic turbulence

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


Introduction
A defining property of weakly collisional plasmas is their ability to support resonant energy transfer between waves and particles.The best-known manifestation of this phenomenon is the famous Landau damping, whereby waves are resonantly damped via energy transfer to the plasma particles (Landau 1946).A contrasting possibility is the case where energy flows in the opposite direction, from the particles to the waves, thereby growing these to nonlinear amplitudes and triggering a wide variety of complex, and often poorly understood, nonlinear plasma behavior.Such kinetic instabilities are fundamental and ubiquitous plasma processes, thought to regulate, or at least significantly contribute to, particle behavior 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 magnetized collisionless plasmas (Roberg-Clark et al. 2018); and mirror and firehose instabilities are thought to regulate the ion distribution and, thus, limit the temperature anisotropy in the Earth's magnetotail (Zhang et al. 2018) and in solar wind (Alexandrova et al. 2013).
Kinetic instabilities may also be crucial in magnetic reconnection events in collisionless plasmas, because the intense wave-particle interactions that they trigger may create an anomalous resistivity (e.g., Sagdeev 1967;Galeev & Sagdeev 1984;LaBelle & Treumann 1988) that breaks the frozen flux condition and potentially sets the reconnection rate (e.g., Ji et al. 1998;Kulsrud 1998Kulsrud , 2001;;Treumann 2001;Uzdensky 2003).Indeed, there is observational (e.g., Deng & Matsumoto 2001;Farrell et al. 2002;Matsumoto et al. 2003) and numerical (e.g., Drake et al. 2003;Zhang et al. 2023) evidence that verifies the existence of turbulent phenomena and wave emissions around diffusion regions of reconnection sites.In addition, kinetic instabilities may influence the reconnection onset (i.e., the sudden transition from a relatively quiescent state to the reconnection stage proper) by disturbing the current sheet formation process and ultimately setting the properties of the reconnecting current sheet (e.g., Alt & Kunz 2019;Winarto & Kunz 2022).
As one of the most prominent examples of streaming-type kinetic instabilities, ionacoustic 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 & II spacecraft revealed the presence of ion-acoustic waves (IAWs) at heliocentric distances between 0.3 and 1 AU (Gurnett & Anderson 1977;Gurnett & Frank 1978).Recently, IAWs have been receiving progressively more attention thanks to ongoing space missions, namely NASA's Magnetospheric Multiscale Mission (MMS) (Burch et al. 2016) and Parker Solar Probe (PSP) (Fox et al. 2016); and European Space Agency's Solar Orbiter (Müller et al. 2013).For example, recent data from the Time Domain Sampler receiver on Solar Orbiter identifies the IAW as a dominant wave mode in the near-Sun solar wind below the local electron plasma frequency (Píša et al. 2021;Graham et al. 2021).Nonlinear structures in the particles' distribution functions associated with nonlinear IAWs, such as ion holes and electron holes, were also recently reported (Mozer et al. 2020(Mozer et al. , 2021a)).In addition, IAWs have also been identified in the separatrix and outflow region on the magnetospheric side of the reconnecting magnetopause (Uchino et al. 2017;Steinvall et al. 2021), which further suggests that IAI may indeed be one essential component of magnetic reconnection in collisionless environments.Understanding the nonlinear evolution of IAI, or 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 (Rudakov & Korablev 1966;Kovrizhnykh 1966;Zavoiski & Rudakov 1967); whereas the latter is mainly studied semiquantitatively based on the Kadomtsev-Petviashvili (KP) model considering wave-ion nonlinear interactions (Kadomtsev & Petviashvili 1962;Petviashvili 1963).In 1969, R. Z. Sagdeev derived an expression for IAT-induced anomalous resistivity by calculating the nonlinear wavenumber spectrum resulting (exclusively) from the nonlinear scattering of waves by ions (Sagdeev 1967;Sagdeev & Galeev 1969).Although Sagdeev's resistivity formula has been widely adopted (e.g., Uzdensky 2003), it is important to note that it is a partially qualitative result, because it fails to consider the angular distribution of IAT.Additionally, 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 V. Yu.Bychenkov and colleagues solved the kinetic equations analytically.Their work simultaneously accounted for the quasi-linear interaction of electrons with waves and the induced scattering of waves by ions, enabling the quantitative analysis of anomalous turbulent transport (Bychenkov & Silin 1981, 1982;Bychenkov et al. 1982Bychenkov et al. , 1984)).By establishing the spectral and angular distribution of ion-acoustic turbulent pulsations, they verified the correctness of Sagdeev's resistivity formula in cases where the scattering of waves by ions dominates the saturation mechanisms.
Despite these achievements, significant uncertainties remain.First, many assumptions or approximations are made in the above-mentioned theories, often excluding possibilities that may be crucial in determining the dynamics of the system.For example, nearly all analytical models assume that a nonlinear steady state will be reached, which is not necessarily true in a realistic situation (and will indeed not be 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 one.
Numerical studies of IAT started in the 1970s, with most simulations performed using the particle-in-cell (PIC) method (e.g., Biskamp & Chodura 1971;Biskamp et al. 1972;Boris et al. 1970;Ishihara & Hirose 1981, 1983;Dum & Chodura 1979;DeGroot et al. 1977).Most of these studies focused on identifying the saturation mechanism of IAT and quantifying the resulting turbulent heating.However, rather than gradually driving the system toward an unstable state -which is the scenario that better conforms to the theoretical approaches to this problem -these simulations start with super-critical initial conditions (their initial condition is such that the electron-ion drift velocity significantly exceeds the threshold for IAI); in some cases, the electron drift velocity is even kept constant throughout the simulation.As a result, they are inadequate to directly verify the analytical theories.Moreover, the relatively small numbers of macro-particles employed in those PIC simulations (due to computational constraints) produced unresolved results.Indeed, even with today's computational resources, simulating resonant kinetic instabilities with the PIC method can be challenging (see, e.g., Tavassoli et al. (2021) for an example concerning the Buneman instability).
The availability of efficient continuum Vlasov-Poisson solvers which arose at the turn of the century (e.g., Fijalkow 1999;Horne & Freeman 2001) enabled a resurgence of numerical studies of IAT.These new continuum Vlasov-based simulations systematically studied the relationship between anomalous resistivity and different simulation parameters such as different types of initial distribution functions (Watt et al. 2002;Petkaki et al. 2003), the initial drift velocity of electrons (Petkaki & Freeman 2008), and phase space structures (Büchner & Elkina 2006;Lesur et al. 2014).Some of these simulations adopted realistic electron-ion temperature and mass ratios (e.g., Petkaki et al. 2006;Hellinger et al. 2004).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 1984).
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 realization depends on the relative rates of cur-rent 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 (2D) effects has been identified in the early numerical work (Biskamp & Chodura 1971).Analytical studies (e.g.Zavoiski & Rudakov 1967;Sagdeev & Galeev 1969) also considered the influence of oblique modes on the instability dynamics.Perhaps as a result of this artificial restriction, the nonlinear evolution observed in these simulations is mostly highly stochastic.Consequently, none of these studies show a direct comparison against the theoretical spectra of IAT pulsations and Sagdeev's formula for anomalous resistivity.
These limitations mean that considerable uncertainty remains about the nonlinear properties of the IAI; in particular, fundamental questions such as its saturation mechanism, the resulting particle heating, and the magnitude of the anomalous resistivity that it drives remain unanswered.In this paper, we conduct a comprehensive numerical study of the nonlinear evolution of the current-driven IAI aimed at providing detailed answers to these questions.We numerically solve the Vlasov-Poisson equations in two dimensions both in real and in velocity space.To ensure the physical realizability of our system, our simulations start from a stable initial condition which is driven gradually towards instability via an imposed external electric field.This setup allows us to present the most complete and detailed understanding of IAT to date.
This paper is organized as follows.Numerical details are provided in section 2. The results from our main simulation are presented and discussed in section 3.In section 4, we analyze 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 section 5.

Numerical Setup
We investigate IAI in a uniform, spatially 2D plasma with periodic boundary conditions (thus effectively mimicking an infinitely large 2D 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 e0 /T i0 ≫ 1).An external electric field E ext is applied in the z (called parallel) direction throughout the simulation to drive the current and, thus, the instability.The reference value against which E ext must be compared is (Bychenkov et al. 1988) where c s0 = T e0 /m i is the ion sound speed, m e is the electron mass, ω pi is the ion plasma frequency, and λ De and λ Di are the electron and ion Debye lengths.In this work, we focus on the so-called weak field case, E ext ≪ E NL .This choice is guided by simplicity considerations: in this regime, wave-wave interactions are predicted to be sub-dominant, and quasi-linear theory is expected to apply (Bychenkov et al. 1988).This is, therefore, a necessary first step in the development of a complete understanding of IAT.We use Gkeyll (Shi et al. 2017(Shi et al. , 2019;;Mandell et al. 2020;Juno et al. 2018;Hakim & Juno 2020), a recently developed state-of-the-art code featuring energy-conserving, Table 1.Summary of the key parameters of the simulations.In the run name, the number following the letter M refers to the mass ratio.The number following the letter E refers to the multiple of 2.5 × 10 −4 for Ẽext.
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 normalized external electric field Ẽext ≡ E ext /(4πen 0 λ De ), where λ De = v T e0 /ω pe is the electron Debye length corresponding to the initial electron temperature.Table 1 lists all simulations and the corresponding values of Ẽext , E NL , mass ratio, and T e0 /T i0 .Run Main is our fiducial simulation, with a mass ratio m i /m e = 100, temperature ratio T e0 /T i0 = 50, and an external electric field of Ẽext = 2.5 × 10 −3 .Most other runs vary the mass ratio or electric field at fixed T 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 × 10 −4 for Ẽ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λDe in the z direction and 25.0λ De in the y direction, with grid resolution ∆z = ∆y = 0.5λ 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 impact 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 section 3, we use our fiducial simulation, Run Main, to demonstrate our main results.The other simulations are employed in section 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.

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 z direction and total wave energy.Based on this figure, we identify five distinct phases of evolution, which we analyze in detail in the subsections below.

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., dJ z /dt = E ext e 2 n 0 /m 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 dJ z /dt) is comparable to that of the most unstable IAW predicted by linear theory (Jackson 1960) (linear theory is briefly described in appendix A).We observe that this criterion predicts the onset of the instability very well.According to linear theory, the growth rate of an ion-acoustic mode, γ 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 γ e (t, k) ≈ γ e (t 0 , k)(t/t 0 ).Therefore, the wave amplitude should grow as We calculate γ 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 Eq. (3.1).This linear phase ends when the wave energy becomes large enough that waveparticle interaction becomes significant and the electron velocity distribution begins to be modified, at around ω pe t ≈ 500.

Phase II: saturation
Phase II captures the saturation stage of IAI.We first investigate the saturation mechanism.As mentioned in section 1, two processes are potentially responsible for saturation: quasi-linear interactions between electrons and waves, and nonlinear interactions with ions (ion trapping (Sagdeev & Galeev 1969;Biskamp & Chodura 1971)).The electrons' and ions' responses to the growth of the IAWs are represented in Figures 2 (a) and (b), respectively, by plotting the one-dimensional (1D) velocity distribution function at different moments of time, obtained by averaging over the spatial domain (z, y) and v y ; namely, F α (v z ) ∝ f α (v z , v y , z, y)dz dy dv 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 ω pe t ≈ 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 (same figure) that strong ion heating (which is the signature of ion trapping) only begins later in time, at around t ≈ 650 ω −1 pe , when the growth rate of wave energy is already significantly reduced 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 summarizing the maximum wave energy density we observe in simulations with different mass ratios and external electric fields W max in table 2 and comparing them with the quasi-linear estimate (Zavoiski & Rudakov 1967;Bychenkov et al. 1988), where W 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 ext and on the square root of mass ratio (m i /m e ) 1/2 .In addition, the absolute value of maximum wave energy density is only a little bit larger than Eq.(3.2).Therefore, we conclude that, under the weak external electric field condition, quasilinear relaxation of the electron distribution, rather than ion trapping, is the main IAI saturation mechanism, consistent with the theoretical prediction (Bychenkov et al. 1988).This conclusion is seemingly at odds with the numerical study by Biskamp & Chodura (1971), which argues instead that ion trapping is the mechanism responsible for saturation.We think that this discrepancy is due to the fact that their simulations start with a super-critical electron distribution.This causes the waves to grow rapidly to large amplitudes.Thus, nonlinear effects, rather than quasi-linear relaxation of electron distribution, dominate the saturation of the IAI in their simulations.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 Table 2.The saturation wave energy of simulations with different mass ratios and external electric fields.Wmax is the maximum wave energy measured from the simulation.Wsat is the quasilinear estimate of the saturated wave energy, i.e., Eq. (3.2).The wave energy exhibits approximate linear dependence on the external electric field and on the square root of mass ratio, consistent with the quasi-linear prediction.
(a) (b) emission (γ e (k); weakened due to plateau formation) and wave damping (γ 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 Eq. (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 (ω pe t = 750) and the cyan (ω pe t = 1000) curves in Figure 2(a).They display two noteworthy features: (i) the (1D) 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 2D feature.
The much more efficient scattering of electrons by the waves in the 2D 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 1D electron velocity distribution at the end of, and after, saturation can be well approximated by a double Maxwellian function, with one Maxwellian representing the bulk of the electron population (denoted with subscript 'b' in the following formula), and the other representing the tail (subscript 't'): . (3.4) 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 (ω pe t ≈ 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 Eq. ( 3.3) still holds at the end of the saturation even though the electron distribution no longer conforms to the plateau shape described by the standard quasi-linear theory (a similar method was previously employed in studies of Buneman instability in reconnection layers (Drake et al. 2003;Che et al. 2009Che et al. , 2010;;Jain et al. 2011)).Because ions do not deviate significantly from a Maxwellian distribution around saturation (see Figure 2(b)), the 1D dielectric function obtained from linear theory can be written as where λ Dα = v T α /ω pα is the Debye length of species α (ions, bulk electrons or tail electrons), ζ α ≡ (ω − ku dα )/kv T α , with ω = ω r + iγ e , u dα the drift velocity, and Z(ζ) the plasma dispersion function.Using the fitting model provided by Eq. (3.4) evaluated at ω pe t = 750, we obtain from this equation γ(k z , ω pe t = 750) ≈ 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 Eq. (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 Eq. (3.3).We note that, as mentioned earlier, Eq. (3.3) is not valid at later times, which voids the assumption made by the quasi-linear theory; see section 3.3.
To summarize, during the saturation process in the 2D 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 behavior of the new bulk, while electrons located on the right-hand side of the resonance region experience continuous acceleration.

Phase III: shutdown of IAT
Phase III is a post-saturation stage characterized by a significantly lower current growth rate (about a factor of 3 lower than the free acceleration rate); see Figure 1.As explained in section 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 α(t = 1000)e 2 E ext /m e , where α(t = 1000) represents the fraction of the fitted tail at t = 1000.This approximation yields a value close to 0.38e 2 E ext /m 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.

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 /(en) ≫ c s the ion-acoustic mode becomes stable for T e /T i ≲ 10 (Papadopoulos 1977; Benz 2012).We confirm this prediction by plotting the temperature ratio, in Figure 3.Even though both electrons and ions are strongly heated after saturation, we observe an abrupt decrease in the electron-to-ion temperature ratio, which happens mostly in phase II and early phase III.This plummeting of the temperature ratio can be directly responsible for a strong damping rate of IAWs for most wave modes.
However, linear theory assumes Maxwellian electron and ion velocity distributions, which is not the case in the nonlinear evolution of IAI at this stage.To prove this conjecture more rigorously, we return to the linear analysis of the modified distribution function to show that IAWs will indeed be Landau damped after saturation.By plugging the fitting parameters at ω pe t = 1000 into Eq.(3.5) and assuming Maxwellian ions, we observe wavenumber to satisfy this condition (3 × 10 −3 ω 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 ∼ 300ω −1 pe for modes with wavenumber satisfying k z ⩾ k c during phase III.By comparing the wavenumber (k) spectrum between ω pe t = 1000 (blue curve) and ω 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λ De /2π ≳ 0.2).However, these wave modes are only about an order of magnitude weaker than those in ω pe t = 1000, which are less damped than predicted.Using the fitting parameters at ω pe t = 1600 in Eq. (3.5), we obtain γ(k z , ω pe t = 1600) < 0 for nearly all values of k z , while the damping rates of wave modes with k z λ De /2π > 0.08 are strong (again, compared to −3 × 10 −3 ω pe ).Even though most wave modes are less damped than predicted, the fact that nearly all the wave modes at ω 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.

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 vortexes (also referred to as "phasespace holes" or "electrostatic solitary waves" (Hutchinson 2017)), which is the signature of electron trapping, can be identified right after saturation (a snapshot of electron phase space at ω pe t = 2000 is shown in Figure 6).This phenomenon is expected by nonlinear theories, and it is also reported in many simulations of two-stream instability (e.g., Morse & Nielson 1969;Berk et al. 1970), space observations where electrostatic streaming waves are present (e.g., Malaspina et al. 2013;Pickett et al. 2008;Mozer et al. 2021a), and laboratory experiments (e.g.Saeki et al. 1979).
One signature of electron trapping is that it can trigger positive frequency shifts of nonlinear IAWs, which may allow the sub-harmonic decay of IAWs (Berger et   As mentioned above, this frequency shift is expected to induce sub-harmonic decay of IAWs and, consequently, modulate the energy spectrum.In Chapman et al. (2014), the onset of IAT leads to ϕ k ∝ k −4 at relatively large wavenumbers.The corresponding wave energy spectrum is then . 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).

Kadomtsev-Petviashvili 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 Kadomtsev-Petviashvili (KP) spectrum, characterized by N (k) ∼ 1/k 3 ln(1/kλ De ), immediately after the saturation stage (ω pe t = 1000).While this spectrum has been observed in certain turbulent heating experiments (Hamberger & Jancarik 1972;Perepelkin et al. 1973;de Kluiver et al. 1991), we believe this to be its first direct numerical validation.
The agreement between our simulation results and the KP spectrum may initially seem puzzling, as the original derivation of the KP spectrum necessitates strong ion nonlinear effects and negligible ion Landau damping (Kadomtsev & Petviashvili 1962;Petviashvili 1963).Although we confirmed earlier in section 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 previously confirmed.In fact, even though quasi-linear effects of electron distribution dominate around saturation, ion nonlinear effects (induced scattering by ions) can still become noticeable afterward, as evidenced by the emerging plateau in the ion distribution function around ω pe t = 1000 (see Figure 2 (b)).Bychenkov & Silin (1982) also demonstrated that the KP spectrum can be obtained by simultaneously considering the quasi-linear relaxation of the electron distribution, leading to a weakening of γ e , and the induced scattering process, leading to a damping rate of γ NL .Specifically, they showed that the condition γ(k) = γ e (k) + γ i (k)+γ NL (k) = 0 can lead to the emergence of the KP spectrum, even when (Bychenkov & Silin 1982;Bychenkov et al. 1988).Our observation confirms the validity of their theory during the times around saturation and emphasizes the significance of nonlinear effects in the subsequent order of weak turbulence theory in shaping the behavior of IAT, even if they are of relatively minor importance in the saturation itself.

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 2D features of the system in Figure 6.Snapshots of the 2D electron velocity distribution are shown in the second column.At ω pe t = 2000, the bulk of the electron distribution is centered 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 ω pe t = 3500.By the end of phase IV, the electron velocity distribution is greatly broadened in the parallel direction with respect to what it was at the end of phase III.Furthermore, we observe that the tail of the electron distribution exhibits a triangular shape, consistent with early analytical conjectures (Sagdeev & Galeev 1969).This is because oblique unstable modes have a maximum propagation angle with respect to the parallel direction.The electron distribution population that lies outside of the resonance region created by the wave modes is more easily accelerated by the external electric field.

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 ω pe t = 2000 and can be clearly seen at ω pe t = 3500 (note that the colormap of the ion velocity distribution plots is in logarithmic scale to better emphasize 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 2D simulation (Dum et al. 1974) and can be explained by the quasi-linear theory for ions (Ishihara & Hirose 1981).Initially, as the ions are extremely cold, only a very small fraction of waves with low phase velocities (v ph < c s ) can directly heat the cold bulk ions through weak ion trapping.However, as time progresses, the quasilinear effect gradually diffuses the non-resonant portion of the ion distribution towards the phase-velocity range of the IAWs, allowing for stronger direct resonant interactions with the waves.This quasi-linear diffusion of non-resonant ions is the primary mechanism contributing to the formation of the high-energy ion tail (Ishihara & Hirose 1981, 1983).
It is interesting to note that similar ion distribution is also observed in a recent 2D2V Vlasov simulation of current-driven instabilities (Chan et al. 2022).
The ion tail formation can also help explain the fact that the most intense wave modes propagate with a non-zero angle with respect to the parallel direction during phase IV, as evidenced by the yellow curve in Figure 4 (b).According to linear theory, parallel modes are expected to dominate over oblique modes because of their higher linear growth rates and saturation levels (because parallel modes have higher effective drift velocity).This is confirmed by our simulations in the linear stage (phase I, the dark blue curve in Figure 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 ω pe t = 2000 and ω 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 ω 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.

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 section 3.3, we mentioned the emergence of electron holes right after saturation in section 3.3.1.During the later course of the nonlinear evolution, these phase-space vortexes 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 ω 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 favor the triggering of new instabilities because they asymmetrically reflect lowenergy 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 ω pe t = 3500, we see clearly larger and more asymmetric potential structures compared with the electric potential at ω pe t = 2000.The impacts on the electron velocity distribution function can be seen in Figure 2(a), where the depletion of the electrons at around v T e0 is already visible at ω pe t = 3500.Another possible consequence of particle trapping is modulational instability.In Rose (2005), it was demonstrated that negative nonlinear frequency shift and wave diffraction of Langmuir waves can trigger a self-focusing effect.It is possible that a similar effect could extend to the IAWs.In Figure 7(a), we indeed observe a stronger negative frequency shift during phase IV, which emphasizes the effects of nonlinear ion trapping at the late stage of the nonlinear evolution (Berger et al. 2013).The modulated wavefronts (the localized maxima in electric potential) observed in the last column of Figure 6 at ω 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 -clearly visible in the rightmost column of Figure 6 at ω pe t = 4300.The formation of double layers associated with IAI has been shown via experiments (e.g., Okuda & Ashour-Abdalla 1982;Chanteur et al. 1983), simulations (e.g., Chanteur 1985;Sato & Okuda 1980), and spatial measurements (e.g., Ergun et al. 2001Ergun et al. , 2009)).We here confirm the production of the double layer by IAT with a 2D2V simulation with a more self-consistent setup.

Phase V: onset of electron-acoustic waves
As evidenced by the wave energy time trace in Figure 1, another instability is triggered at this moment.We interpret these new wave modes as electron-acoustic waves (EAWs) (Holloway & Dorning 1991;Valentini et al. 2006Valentini et al. , 2012;;Anderegg et al. 2009) for the reasons described below.
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ω pe to around 0.6ω pe and a much faster phase velocity (∼ 1.2v T e0 ), which excludes the possibility of (slow) ion modes.Instead, the intermediate frequency range between ion and electron plasma frequencies is standard for EAWs (Holloway & Dorning 1991).Second, as seen from the electron phase structure at ω pe t = 4300 in the first column of Figure 6, the electron distribution has a trapped population around v T e0 with a trapping width extending to ∼ 4v T e0 .These trapped electrons around v T e0 effectively create a plateau shape around the phase velocity of the waves (i.e., ∼ v T e0 ) in the electron velocity distribution (which can be seen from

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

Anomalous resistivity intensity
We define the effective collision frequency, ν eff , from the equation dJ z /dt = e 2 n 0 /m e E ext − ν eff J z .A reference value is the quasi-linear result (Bychenkov et al. 1988): Shown in Figure 8(a) are time traces of ν eff for runs with the same external electric field but different mass ratios.It can be seen that the observed ν eff in phase III is only about 10% ∼ 20% of the value predicted by Eq. (4.1) for Run Main.The maximum ν eff , which is observed during saturation (in phase II), is also smaller than the quasi-linear value (about 38% of Eq. (4.1)).In terms of the mass ratio dependence, although the peak values depend weakly on the mass ratio, ν eff after saturation shows a roughly linear dependence on m i /m 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 1966;Bychenkov et al. 1988).This assumption is only valid for the resonant part of the electrons.However, both resonant and non-resonant parts of electrons are accounted for in the current.Due to the fact that most of the tail, with a significantly higher drift velocity, does not resonate with the waves and is freely accelerated by the external electric field, the ν 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 ν eff /ν QL eff for simulations with the same mass ratio (m i /m 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 ext ≪ E NL ), the linear dependence of ν eff on the external electric field holds.For example, ν eff in Run M25E10, Run M25E8, and Run M25E6 have nearly the same normalized 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 Eq. (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 section 3 still holds, including the saturation and shutdown mechanisms, and transition into bursts of EAWs.An estimate of the threshold electric field and more details on the extremely weak field regime can be found in appendix E.
Another well-known formula for anomalous resistivity, due to Sagdeev (Sagdeev 1967;Zavoiski & Rudakov 1967), predicts an effective collision frequency of ∼ 0.06ω 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 ext > E NL .Finally, to compare our simulations with previous numerical studies ( LaBelle & Treumann 1988;Watt et al. 2002;Hellinger et al. 2004), it is convenient to calculate a more qualitative estimate, namely, 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.Eq. ( 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 & IV).This conclusion is at odds with values of anomalous resistivity exceeding this estimate by at least one order of magnitude reported in many previous numerical studies (e.g., Watt et al. 2002;Hellinger et al. 2004;Petkaki et al. 2003Petkaki et al. , 2006)).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 Eq. ( 4.3).
To summarize, 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, Eq. (4.1), is a good approximation.Otherwise, the anomalous resistivity after saturation is about 10% ∼ 20% of Eq. (4.1) in our simulations.We also find that Eq. ( 4.3) provides a reasonable fit.However, even though we have validated Eq. (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.

Particle heating and anomalous resistivity duration
As explained in section 3.3, anomalous resistivity will eventually become negligible as a result of the effective shutdown of the IAT.Therefore, in addition to characterizing 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.

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 ∼ 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 i,f and T 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 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 i,f for simulations with different mass ratios and external electric fields, showing a linear dependence on m i /m e and E ext .This is the same as the dependence of W sat on E ext and m i /m e in Eq. (3.2) and confirmed in table 2. We thus fit the ion temperatures in most of our simulations as According to Eq. (4.4) and Eq.(3.2), the electron temperature at the end of phase III can then be fitted by Ẽext n 0 T e0 .(4.6) However, there are two simulations that do not conform to Eq. (4.5) in Figure 10 -Run M25E1 and Run M25E2 (the leftmost two blue rectangular data points).As discussed in section 4.1, these two simulations are under the extremely weak field regime.In fact, the final ion temperature T i,f given by Eq. (4.5) can potentially become insufficiently low to satisfy the shutdown criterion, as described in Eq. (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 NL -which we call the extremely weak field regime.Under this regime, as mentioned above, 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 Eq. (4.5) until the shutdown criterion Eq. (4.4) is finally met.

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/2 ω pe , which has the same dependence on external electric field and mass ratio as the Eq.(4.1).Then we can use τ res ≈ (T e,f − T e0 ) /(dT e /dt) and substitute T e,f with Eq. (4.6).We arrive at While the exact value of θ depends on parameters, our simulations suggest that θ ≈ 10 (it is slightly larger for cases with a larger mass ratio and a smaller external electric field, see Figure 9).This is consistent with linear theory (Papadopoulos 1977;Benz 2012), and is justified by the fact that IAI becomes stable when T e /T 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 switchoff criterion accurately.This may explain the deviations of θ from the value 10.Note that Eq. (4.7) does not apply to extremely weak field regime due to the underestimate of electron and ion heating.With Eq. (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, Eq. (4.7) yields τ res ∼ 2000ω −1 pe if we adopt θ = 10, which is in reasonable agreement with the numerical data shown in Figure 8(a).

Conclusion and discussion
In this study, we present a comprehensive first principle (i.e., continuum Vlasov-Poisson) numerical investigation of the current-driven IAI in unmagnetized plasmas.By starting from a stable plasma equilibrium and gradually driving it towards instability using a (weak) external electric field, we ensure the self-consistent evolution of the system towards instability and the physical realizability of the unstable states that are reached.
We identify 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 2D (in position space) system is much more efficient than its 1D 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, contrary 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.Additionally, 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 Kadomtsev-Petviashvili (KP) spectrum around saturation, which had not been reproduced numerically ever before and emphasizes 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 waveinduced 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 shift our focus to the investigation of anomalous resistivity resulting from IAT and its dependence on various simulation parameters.We study the dependence of electron and ion heating on those parameters and derive an estimate for the time duration of significant anomalous resistivity.We find that the quasilinear 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 1D numerical studies, where stronger anomalous resistivity is typically observed.We also observe, and derive analytically, a new electric field threshold below which the runaway electron population is small and the quasi-linear prediction of anomalous resistivity holds.
Our study of the classical version of the IAT can be directly applicable to some physical systems.For instance, IAI can exist in the separatrix region of reconnecting current sheets, where the magnetic field may be negligible (in the no guide-field case).Situations where electrons are heated to high temperatures while ions remain cold, such as stellar flares (Polito et al. 2018) and the low-speed solar wind near the Sun (Verscharen et al. 2022;Mozer et al. 2022), have been observed.The findings in this paper are therefore directly relevant to such environments.To be more specific, our simulations reveal features that match with direct observations in various plasma environments.For example, the plateaued-shaped ion velocity distribution with a high-energy tail which we find in our simulations is qualitatively similar to those measured in the solar wind (Mozer et al. 2020) and interplanetary shocks (Wilson et al. 2020), where IAWs are believed to be present.We also demonstrate the IAT's capability to significantly heat the electron core and potentially create the "strahl" in the electron distribution commonly observed in the solar wind (Verscharen et al. 2019).The development of anisotropy in particle velocity distribution is qualitatively similar to the measured electron distribution on the dayside of reconnecting magnetopause during IAW bursts (Uchino et al. 2017;Steinvall et al. 2021).Additionally, bursts of EAWs following IAT are experimentally observed in laser-driven reconnection events (Zhang et al. 2023).We also observe that oblique modes dominate over parallel modes in the later stages of nonlinear evolution.When the driving field is strong, the dominance of oblique modes is expected to occur even earlier (Bychenkov et al. 1988).These findings align with recent IAW observations in the solar wind (Mozer et al. 2021b), where the most intense waves propagate obliquely.The broadening of the spectrum during the nonlinear stage, as reported in other studies (Mozer et al. 2021a), is also consistent with the wavenumber spectrum observed in our simulations (see Figure 4 (a)).
Our results can also be adapted to study magnetic reconnection influenced by IAT.The electric field associated with reconnection is E rec ∼ γV A B 0 , where γ is the reconnection rate and V A is the Alfvén speed computed with the reconnecting field B 0 .Situations where this electric field is smaller than E NL (Eq.(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 e /T i ∼ 5 and T e ∼ 50eV (Mozer et al. 2021b).For these parameters, we find E rec /E NL ∼ 1.2π(1/β e )(v T e /c) m i /m e (T i /T e ) < 1.In this case, assuming a fast reconnection rate of γ ≈ 0.1, the reconnecting electric field is about E rec /(4πen 0 λ De ) ≈ 2γβ e m e /m i (v T e /c) ≈ 5 × 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 Eq. (4.5) (ion heating in Run M25E2 is about twice as intense as the prediction given by Eq. (4.5)).Then, according to Eq. (4.7), the corresponding τ res is at most on the order of 10 3 ω −1 pi .On the other hand, the typical reconnection time, τ rec = γ −1 L/V 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 km at distances of this order (Coles & Harmon 1988;Lotz et al. 2023;Huang et al. 2023).At such a scale, the corresponding V A can be estimated using the Alfvén velocity at the outer scale and the observed magnetic power spectrum; a value of V A ≈ 100km/s seems reasonable (Huang et al. 2023).Assuming a fast reconnection rate of γ ≈ 0.1, τ rec can easily extend beyond 10 5 ω −1 pi .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 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 λ 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 λ De ) is small compared to 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 m and, so, 10 4 λ De corresponds to a distance ∼ 100km, which is still smaller than the ion skin depth (∼ 140km) (Verscharen et al. 2019).On the other hand, the limitation imposed by the small size of the simulation box, especially in the parallel direction, excludes the effects of long-wavelength wave modes that become prominent in phase V.The box-size double layer that we observe can result from this artificial constraint.Despite this limitation, we believe that the simulations conducted up to phase V and the mechanism elucidated for triggering the second instability remain plausible.Expanding the simulation size in the parallel direction for a 2D2V simulation, while desirable, was deemed impractical due to the large computational resources it would require.In summary, our numerical investigation represents a significant stride toward a self-consistent understanding of the nonlinear dynamics of IAT.While we concentrate on the weak field regime in this step, which contains simpler (though nonetheless complex) nonlinear physics than the strong field case, our study yields substantial insights applicable to various physical scenarios where IAWs are present.Our study lays a robust foundation for future investigations into more complex IAT situations.These encompass strong external fields, the presence of magnetic fields, and less extreme initial temperature ratios.It is also important to acknowledge that our simulations do not encompass the hypothesized 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.Figure 12 shows plots of the growth rates (γ), growth rates divided by wavenumbers (γ/k), and phase space velocities (ω/k) as functions of the wavenumber, obtained by solving Eq. (A 3) for cases with different mass ratios.Both electrons and ions have a Maxwellian velocity distribution.The electron distribution has a drift velocity of 0.5v T e0 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 one-dimensional simulation
We run a 1D1V simulation with the same physical parameters as Run Main. Figure 13 shows the differences between the 1D and 2D simulations.The simulation with the label of "2D" is Run Main.The evolution of the 1D 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 1D case, which can be due to the different retained plateau shapes of the electron velocity distribution as explained below.
We plot 1D electron velocity distribution functions in Figure 14 right after saturation (ω pe t = 1000) and at the end of phase III for the 2D case (ω pe t = 1900).In the 1D 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 righthand 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.On the contrary, in the 2D 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 1D case, and scatters most of the electrons back to the resonance region (the region with lower v z ).Consequently, the 1D electron distribution exhibits in the 2D 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 1D and 2D runs result in very different nonlinear dynamics and demonstrate that the physics of IAT is not well captured by 1D models.where u(t) = eE ext t/m e is the bulk electron drift velocity (Schekochihin 2022).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 sat , can be estimated by quasi-linear theory; it is . (E 3) The time at which saturation occurs, i.e., t 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 ∆v trap is the velocity width (in the parallel direction) of an electron hole in phase space (i.e.trapping width).∆v trap can be estimated using eϕ ∼ m e ∆v 2 trap , (E 7) where ϕ is the typical electric potential at the moment of saturation in the system.We can estimate ϕ using the saturation wave energy density, i.e.
where we take k to be the wave number of the mode with the largest linear growth rate .

(E 9)
The corresponding trapping width is then . (E 10) Substituting t sat and Eq.(E 10) into Eq.(E 6), we obtain the inequality  which is consistent with our simulation results.Run 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) and (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 e /T i ∼ 10) is reached at around ω 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 (b) and (d) of Figure 17, are also quickly turned off.where k max is the largest unstable wavenumber of IAWs, which can be predicted from linear theory.Conveniently, we find that k 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 max λ De ≈ 2π × 0.36 for all mass-ratio cases.Plugging this k max into Eq.(F 1), we obtain

Appendix F. Determining the numerical resolution
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 ph,min , v ph,max ] i.e., This method was adopted in Petkaki et al. (2006).
The second method estimates the width of the electron resonance region due to the IAWs when the system enters the non-linear 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 γ k (Sagdeev 1967).If the maximum electric potential of the wave mode is ϕ, the bouncing frequency is ω B = ek 2 ϕ/2m e .Letting ω B ∼ γ k , we obtain ϕ ∼ m e γ 2 k /ek 2 .The resonance width in electron velocity space, ∆v, caused by this wave is ∆v ∼ eϕ/m e , and the corresponding resolution is To resolve most of the unstable waves, estimates given by Eq. (F 3) and Eq.(F 4) produce similar results according to Figure 12 (b) and (c) (we truncate the wave number at kλ De = 0.3 when estimating with Eq. (F 4)).From Figure 12 (b) and (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 summarize them here: ≈ 0.03, (F 5) ∆v v T e0 mi/me=100 ≈ 0.015, (F 6) ∆v v T e0 mi/me=400 ≈ 0.007.(F 7)

F.2. Linear benchmarks
We perform a set of benchmarks against linear theory to check both the simulation setup and the numerical resolution.In these simulations, no external electric field is applied.Electrons have an initial drift velocity (u e = 0.5v T e0 , which is a typical drift velocity during phase I of most of our simulations in the paper) such that the initial 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 ∆v e,z = 0.042v T e0 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 setup for these turbulent tests is the same as 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 Ẽext = 2.5 × 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 behavior.
The simulation that deviates most from the others is combination #1: it exhibits higher ion heating than the others and a somewhat shorter phase III and phase 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 behavior 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 is less impactful 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 summarize, as the real space numerical resolution does not depend on mass ratio, the grid size for all the simulations can be chosen as {∆z, ∆y} = {0.5, 0.5}λ De .
(F 8) The required velocity space numerical resolution linearly depends on the inverse square root of mass ratio (See Figure 12 and Eq.(F 5-F 7)), 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 {∆v e,z , ∆v e,y } = {0.021,0.042}v T e0 .
(F 9) For Run M200E10, we are forced by computational resource constraints to keep this numerical resolution (Eq.(F 9)).While insufficient to adequately capture the linearly unstable spectrum, these nonlinear tests demonstrate that phases III-IV of the evolution should still be reasonably resolved.

Figure 2 .
Figure 2. 1D electron (a) and ion (b) (in logarithmic scale) velocity distribution at different times.The shadowed region is the approximate electron resonance region at the moments around saturation (ωpet ≈ 750).The black dashed line represents the fitted bulk in the electron velocity distribution after saturation, and the dashed blue and green curves represent the fitted tail at different times.The fitting model is the double Maxwellian function, Eq. (3.4).See text for more details.

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

Figure 4 .
Figure 4. Wavenumber (a) and angular (b) spectra of IAWs at different times during the simulation.Angular spectra at different times are normalized to their value at θ = 0 • , where θ is the angle between the parallel direction and wave propagating direction.A fourth-order polynomial is used to fit the angular spectra to smooth the curves.The wavenumber spectrum at ωpet = 1000 agrees with the Kadomtsev-Petviashvili (KP) spectrum (the dashed curve).The wavenumber spectra at ωpet = 2000 and ωpet = 3500 with relatively larger wavenumbers agree with the spectrum observed in Chapman et al. (2014) (the dash-dotted curves).

Figure 5 .
Figure 5. Parallel wave energy spectrum |Ez(ω, kz)| 2 during the IAW burst phase (i.e., Phase I) (a), and the phase right after saturation (i.e., Phase II-III) (b) for Run M25E10.The dispersion relation of IAWs (Eq.(A 3)) is plotted with the red dashed curve in (a) and (b).The blue dash line in (b) indicates a set of waves with positive frequency shifts in the nonlinear stage.

Figure 6 .
Figure 6.Snapshots of electron phase-space (first column), electron velocity distribution (second column), ion velocity distribution (third column), and electric potential (last column) at the end of phase III, the end of phase IV, and the middle of phase V, respectively.The color bar for the ion velocity distribution is in logarithmic scale.The 2D distributions are obtained by averaging the 4-dimensional (2D2V) distribution over the other two dimensions.

Figure 7 .
Figure 7. Parallel wave energy spectrum |Ez(ω, kz)| 2 after effective shutdown of IAT (i.e., Phase IV) (a), and the EAW burst phase (i.e., Phase V) (b) for Run M25E10.The dispersion relations of IAWs and EAWs are plotted with the red dashed curve in (a) and (b), respectively.The blue dash-dotted line in (a) indicates a set of waves with negative frequency shifts.The blue dash line in (b) indicates a set of waves with a phase velocity of 1.2vT e0.

Figure 8 .
Figure 8.(a) Time traces of effective collision frequency for simulations with different mass ratios and the same external electric field ( Ẽext = 2.5 × 10 −3 ).For Run Main we also plot, as a reference, the curve corresponding to 0.5W/n0Te, where W is the wave energy density (dotted curve), and 10% of the quasi-linear value (Eq.(4.1)) (dashed horizontal line).(b) Time traces of effective collision frequency normalized to the quasi-linear prediction, Eq. (4.1), for simulations with different external electric fields and the same mass ratio (mi/me = 25).

Figure 9 .
Figure 9.Time traces of temperature ratios for simulations with different mass ratios, electric fields, and initial temperature ratios.The black dashed line marks Te/Ti = 10.

Figure 10 .
Figure 10.Final ion temperature (normalized to its initial value) as a function of the square root of mass ratio at fixed Ẽext = 2.5 × 10 −3 (top horizontal axis, red) and as a function of the external electric field at fixed mi/me = 25 (bottom horizontal axis, blue).

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

Figure 12 .
Figure 12.Growth rates (a), growth rates divided by wavenumbers (b), and phase velocities (c) of IAWs for different mass ratio cases obtained from solving Eq. (A 3).

Figure 13 .
Figure 13.Time traces of wave energy (a), ion temperature (b), current in the parallel direction (c), and effective collision frequency (d) from 1D and 2D simulations with the same physical parameters.

Figure 14 .
Figure 14.1D electron velocity distribution function, Fe(vz) obtained from 1D (dashed curves) and 2D (solid curves) simulations with the same physical parameters.The first time moment plotted (red curves) is right after saturation, and the second time moment (blue curves) is at the end of phase III for the 2D simulation.

Figure 15 .Figure 16 .
Figure 15.Growth rates of different wave modes obtained by solving Eq. (3.5) at different times.The electron distribution function is represented by the double Maxwellian model, Eq. (3.4), with the fitted parameters summarized in table 3. Ion distribution is a single Maxwellian function with the ion temperature and the drift velocity at the corresponding time.The blue dashed curve is obtained by artificially increasing the ion temperature to that at ωpet = 1900.The black dashed marks γ = 0.The black dotted line marks the relatively strong damping rate during phase III (γ = −0.003ωpe).

Figure 17 .
Figure 17.Time traces of wave energy (a), ion temperature (b), current in the parallel direction (c), and effective collision frequency (d) from Run M25E2.A clear shutdown moment is visible from these time traces.

F. 1 .
Ideal resolutionThe 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

Figure 19 .
Figure 19.Time traces of wave energy (a), ion temperature (b), current in the parallel direction (c), and effective collision frequency (d) from nonlinear evolution tests with different numerical resolution combinations (indicated in table 4).
Time trace of total wave energy (red curve, right axis) and parallel current (blue curve, left axis) of Run Main.We split the evolution of the system into five different stages; see text for details.The red dotted vertical line indicates the moment when the current growth rate is comparable to that of the most unstable IAW predicted by linear theory.The black dashed line is W (t) predicted by linear theory, see Eq. (3.1).The green dashed line indicates the current growth rate during phase III (which is 30% of the free acceleration rate), and the orange dashed line indicates the current growth rate during phase IV (which is 85% of the free acceleration rate).

Table 3 .
Fitting parameters of the double Maxwellian model, Eq. (3.4), for Run Main at different times.