Effects of multi-dimensionality and energy exchange on electrostatic current-driven plasma instabilities and turbulence

Large-amplitude current-driven plasma instabilities, which can transition to the Buneman instability, were observed in one-dimensional simulations to generate high-energy back-streaming ions. We investigate the saturation of multi-dimensional plasma instabilities and its effects on energetic ion formation. Such ions directly impact spacecraft thruster lifetimes and are associated with magnetic reconnection and cosmic ray inception. An Eulerian Vlasov–Poisson solver employing the grid-based direct kinetic method is used to study the growth and saturation of 2D2V collisionless, electrostatic current-driven instabilities spanning two dimensions each in the configuration (D) and velocity (V) spaces supporting ion and electron phase-space transport. Four stages characterise the electric potential evolution in such instabilities: linear modal growth, harmonic growth, accelerated growth via quasi-linear mechanisms alongside nonlinear fill-in and saturated turbulence. Its transition and isotropisation process bears considerable similarities to the development of hydrodynamic turbulence. While a tendency to isotropy is observed in the plasma waves, followed by electron and then ion phase spaces after several ion-acoustic periods, the formation of energetic back-streaming ions is more limited in the 2D2V than in the 1D1V simulations. Plasma waves formed by two-dimensional electrostatic kinetic instabilities can propagate in the direction perpendicular to the net electron drift. Thus, large-amplitude multi-dimensional waves generate high-energy transverse-streaming ions and eventually limit energetic backward-streaming ions along the longitudinal direction. The multi-dimensional study sheds light on interactions between longitudinal and transverse electrostatic plasma instabilities, as well as fundamental characteristics of the inception and sustenance of unmagnetised plasma turbulence.


Introduction
Two categories of electrostatic current-carrying plasma instabilities are typically considered (Omura et al. 2003;Mikellides et al. 2005).The current-carrying ion-acoustic instability manifests when the net electron drift speed Ũd = Ũe − Ũi exceeds the ionacoustic speed k B Te / mi and Te ≫ Ti , where Ũe , Ũi , k B , Te , Ti , and mi are the electron bulk velocity, ion bulk velocity, Boltzmann constant, electron temperature, ion temperature, and ion mass, respectively, and the tilde notation over problem parameters denotes dimensional quantities (Stringer 1964 and references therein).Over a broad range of Te / Ti , the Buneman instability arises when the net electron drift speed exceeds the electron thermal speed.In a 1D configuration, the threshold is Ũd ≳ 1.3c e , where ce = k B Te / me and me is the electron mass (Buneman 1959).Depending on Ũd , both instabilities can be physically relevant in the same plasma configuration.In such instabilities, a sufficiently fast net electron drift triggers plasma wave growth.At a significant wave amplitude, a transition to multi-scale electric potential structures occurs, which in turn induces localised ion acceleration with broadband phase-space structures of an ostensibly turbulent nature.
An example in which multi-dimensional current-carrying instabilities play an important role is the plume of a hollow cathode, e.g., for electric spacecraft thrusters, surface processing, and plasma discharges (Oks & Schanin 1999;Becker et al. 2006;Goebel et al. 2021;and references therein).The erosion of cathode structures can limit the lifetime of years-long outer-space missions by an order of magnitude to less than O(10 4 ) hours (e.g., Friedly & Wilbur 1992;Kameyama & Wilbur 2000;Williams Jr. et al. 2000;Mikellides et al. 2005Mikellides et al. , 2007Mikellides et al. , 2008Mikellides et al. , 2015;;Goebel et al. 2007;Jorns et al. 2014;Lev et al. 2019).A key cause of cathode sputtering is surface collisions of energetic ions at its orifice with energies corresponding to up to 100 eV (Friedly & Wilbur 1992;Williams & Wilbur 1992;Kameyama & Wilbur 2000;Williams Jr. et al. 2000;Boyd & Crofton 2004;Goebel et al. 2007;Mikellides et al. 2008;Farnell et al. 2011).Such highenergy ions are currently postulated to arise in part from the aforementioned collisionless, electrostatic current-carrying instabilities (Williams Jr. et al. 2000;Mikellides et al. 2005Mikellides et al. , 2007Mikellides et al. , 2008Mikellides et al. , 2015;;Goebel et al. 2007;Jorns et al. 2014;Lopez Ortega & Mikellides 2016;Jorns et al. 2017;Sary et al. 2017a,b;Hara & Hanquist 2018;Lopez Ortega et al. 2018;Hara 2019).In particular, the presence of transversely energetic ions has been observed experimentally (Boyd & Crofton 2004;Goebel et al. 2007;Farnell et al. 2011;Hall et al. 2019), and recent laser Thomson scattering measurements show evidence of bi-directional plasma waves at oblique orientations to the local magnetic field (and uni-directional along the parallel direction) in the plume of a hollow cathode (Tsikata et al. 2021).Despite the increasing experimental evidence of multi-dimensional electrostatic plasma waves amplified due to kinetic instabilities, there have not been many detailed numerical studies simulating and analysing the long-term behaviour of such waves to date.Characterising these precursors of the erosion process can improve predictions of thruster lifetime and complement accelerated life testing.Additionally, high-energy ions that are accelerated in the longitudinal and transverse directions can be relevant to the inception of magnetic reconnection and its resulting auroral emissions, where high-frequency turbulence and subsequently intense currents are formed during changes in the magnetic topology (Quon & Wong 1976;Mozer et al. 1980;Temerin et al. 1982;Lysak 1990;Cairns et al. 1995;Ergun et al. 1998;Drake et al. 2003;and references therein).
The generation of axially energetic ions has been numerically investigated through Eulerian Vlasov-Poisson simulations employing the grid-based direct kinetic (DK) method in a single spatial dimension (1D) (Hara & Treece 2019;Vazsonyi et al. 2020).In this paper, we extend the analysis to two spatial dimensions (2D) to investigate how ions are energised in the transverse (radial) direction.More broadly, we seek a fundamental understanding of the primary mechanisms driving the growth and saturation of multidimensional plasma waves driven by current-carrying instabilities, as well as the resulting ion and electron velocity distributions after the system reaches a state of saturated turbulence.Such a study involves simultaneous analysis of the transition pathway to turbulence of the accompanying electric field and then the ion phase-space distribution.We build on previous numerical studies of the 2D Buneman instability (Amano & Hoshino 2009) but focus on late ion acceleration instead of early electron acceleration and use a physical mass ratio close to that of a hydrogen plasma ( mi / me = 1.8229×10 3 ).Following the pioneering work by Forslund & Shonk (1970) on electrostatic counterstreaming instabilities, Chapman et al. (2021) investigated the non-linear evolution of the ionion streaming instability for warm electrons and counter-streaming ions with streaming speeds on the order of the ion-acoustic speed.We analyse the growth and saturation characteristics of 2D current-carrying instabilities for an ion-electron system with larger net drift speeds on the order of the electron thermal speed, focusing on plasma wave growth due to the presence of the net current.In addition, we take a deeper dive into the key stages underlying the genesis of multi-dimensional electrostatic plasma turbulence, and examine their implications on energetic ion formation and hollow cathode sputtering.Several Ũd /c e ratios are considered to span possible cathode operating conditions.A number of 1D and 2D instability characteristics are also compared to establish connections with previous 1D work and elucidate the role of multi-dimensionality.In this work, we consider Te / Ti = 10, which is representative of cathode operating conditions (Goebel et al. 2005;Mikellides et al. 2005Mikellides et al. , 2007Mikellides et al. , 2008;;Farnell et al. 2011;Jorns et al. 2017).
The objective of this work is to analyse the stages of collisionless, electrostatic currentcarrying instability growth, as well as the inception of multi-dimensional electrostatic plasma turbulence and its spectral characteristics.In § 2, we describe our computational and physical set-up.In § 3, we discuss results of the 2D instability.These are compared with key results from the 1D instability in § 4. Conclusions are provided in § 5.

Methodology
While the particle-in-cell method is the predominant kinetic simulation approach, such particle-based kinetic methods are susceptible to statistical noise.This is because a sufficiently large number of super-particles is required in each computational cell, together with sufficiently long time averaging for steady and quasi-steady configurations, to adequately sample the particle distribution function in the configuration and velocity spaces (Chen & Boyd 1996;Kannenberg & Boyd 2000;Farbar & Boyd 2010).This issue is particularly exacerbated in current-carrying instabilities whose phase velocity coincides with the tails of the ion velocity distribution function (VDF), where the distribution magnitude is small but the associated ion energy is large.The corresponding ion trapping needs to be resolved with a large number of super-particles per cell.Moreover, these instabilities are transient and preclude time averaging as a viable technique for statistical convergence.Such statistical noise is eliminated in the grid-based direct kinetic method, where the distribution function is directly simulated in a discretised phase space without the use of super-particles (Filbet & Sonnendrücker 2003;Kolobov et al. 2007;Thomas et al. 2012;Dimarco & Pareschi 2014;Hara & Hanquist 2018;Palmroth et al. 2018;and references therein), keeping in mind that an accurate solution necessitates sufficient resolution in both the physical and velocity space dimensions, the latter particularly to resolve particle trapping (e.g., Schamel 2000;Dieckmann et al. 2004a).We leverage this superiority of the direct kinetic method to gain more precise insights into instability development and transition to turbulence.The direct kinetic method has been shown to yield more accurate linear instability growth rates than the particle-in-cell method due to the statistical noise inherent in the latter, which is seen to vary with the number of particles per cell and the wavenumber of interest (Dieckmann et al. 2004b;Tavassoli et al. 2021).We verify these growth rates in § 3.1 and also demonstrate the utility of the large dynamic range of direct kinetic solvers in resolving turbulence structures in § 3.2.

Direct kinetic solver
The direct kinetic solver employed here was originally developed at the University of Michigan with verification and validation against canonical and complex plasma problems, such as waves, electron-emitting sheaths, and Hall thruster discharges (Hara et al. 2012(Hara et al. , 2014(Hara et al. , 2015(Hara et al. , 2017;;Hara & Hanquist 2018;Raisanen et al. 2019;Vazsonyi et al. 2020).In contrast to state-of-the-art particle-in-cell solvers, direct kinetic solvers eliminate statistical noise and are suitable for precise investigations of instability growth and turbulence inception as discussed above.For collisionless, unmagnetised plasmas, the solver computes the time evolution of the probability density function, f * , for some particle type * = i, e (ions, electrons) according to the Vlasov equation where x, ṽ, and Ẽ respectively denote the position, velocity, and electric field vectors, q * denotes the charge of the simulated particle type, and t denotes the time.This is consistent with the conservation of charge in phase space.The computational domain is discretised in x-ṽ space with a parallelised second-order finite-volume method, which is described by Chan & Boyd (2022a,b).Assuming small induced magnetic fields and their rates of change, we adopt the electrostatic approximation, where Gauss's law is expressed using Ẽ = −∇ x φ as a Poisson equation for the electric potential φ of the form where ε 0 and e respectively denote the vacuum permittivity and elementary charge, and ñi and ñe respectively denote the ion and electron number densities ñi x; t = ṽ fi x, ṽ′ ; t dṽ ′ ; ñe x; t = ṽ fe x, ṽ′ ; t dṽ ′ .
(2.3) Periodic and no-flux boundary conditions are employed for x and ṽ, respectively.

Problem set-up and linear stability analysis
We consider 1D1V and 2D2V current-carrying instabilities, where D and V denote the configuration and velocity spaces, respectively.The corresponding simulations are respectively two-and four-dimensional.Since the initial ion-acoustic wave excites longwavelength modes, we choose domain lengths sufficiently larger than the Debye length, ).The initial species temperatures are Te = 2 eV and Ti = 0.2 eV in consistency with the temperature ratio described in § 1.In addition, all species are initialised with Maxwellian velocity distributions.The electrons have a initial net axial drift described by the initial electron Mach number M e = Ũe /c e , which corresponds to a shifted Maxwellian initial condition.The ions have zero initial mean speed, i.e., M i = 0, which corresponds to a stationary Maxwellian initial condition.This configuration is relevant, for example, to hollow cathode plumes, which are biased by direct-current voltages.As such, the initial conditions for this study assume a stable current.This is to be distinguished from counterstreaming studies (e.g., Forslund & Shonk 1970;Bret & Dieckmann 2008;Chapman et al. 2021).Hereinafter, we define x as the axial direction and y as the transverse direction in 2D.Correspondingly, v x and v y denote the axial and transverse velocity dimensions, respectively.Linear growth rates for the current-carrying instability can be analytically predicted via solution of the following linear dispersion relation, which is derived for one-and multidimensional electrostatic waves travelling at an arbitrary angle θ to the axial direction across the background ion and electron populations defined above (stationary and shifted Maxwellians, respectively) with initial species Mach numbers M * as follows where k = k and ω are, respectively, the modal wavenumber magnitude and angular frequency of the electrostatic wave in question, ω * = (ñ * e 2 ) / ( m * ε 0 ) is the plasma frequency, c * = k B T * / m * is the species thermal speed, and Z is the plasma dispersion function corresponding to a Maxwellian velocity distribution (2.5) Unless otherwise stated, lengths, velocities, and time are hereinafter respectively nondimensionalised by λD , c * , and the inverse electron frequency 1/ω e = λD /c e .Specifically, k = kλ D , x = x/ λD , f * = f * c * , v = ṽ/c * , and t = ωe t.Note that the characteristic ion time scale, 1/ω i , and correspondingly the ion-acoustic period, exceed the characteristic electron time scale, 1/ω e , by a factor of mi / me ≈ 40.In this work, we consider the evolution of the following energy modes (all listed in their dimensional and volume-integrated forms): the electrostatic potential energy , and the random kinetic energy ñ * k B T * d Ṽ .In particular, we consider the kinetic energies in detail in § 4.

2D current-carrying instability
Building on the preliminary work of Vazsonyi (2021), the 2D2V instability is simulated with domain extents x, y ∈ [0, L = 80], v x,i , v y,i ∈ [−45, 45], v x,e , v y,e ∈ [−8.5, 8.5] and resolutions ∆x = ∆y = 0.4, ∆v x,i = ∆v y,i = 0.7, ∆v x,e = ∆v y,e = 0.1 for the spatial, ion velocity, and electron velocity dimensions, respectively.These are in line with the gridpoint recommendations of Chan & Boyd (2022b) for adequate resolution of macroscopic quantities and gradients.The corresponding number of grid points are N x = N y = 200, N vx,i = N vy,i = 128, N vx,e = N vy,e = 192.The simulations were performed on 96 Intel Xeon Gold CPUs over a total wall-clock duration of about 4 months.The baseline simulation is converged with respect to the spatial grid and domain extent (not shown here), as well as the velocity grid (see appendix A), for the macroscopic quantities of interest.For the baseline simulation, ∆t = 0.06 and M e = 2.3 using the initial and boundary conditions detailed in § 2. Hereinafter, electric potentials and field strengths are normalised by the thermal potential, φth = k B Te /e, and the corresponding characteristic field strength, φth / λD , respectively.Specifically, ϕ = φ/ φth and E = Ẽλ D / φth .

Stages of instability growth: electric potential and energy spectra
Figures 1 and 2 plot the modal potential growth rates over three time intervals for the baseline simulation introduced above but with twice the velocity resolution.These are based on the potential spectrum E ϕϕ = φ φ * , where the hat notation denotes a Fourier transform and the superscripted asterisk denotes the complex conjugate operation.In the linear stage depicted in figure 1, the simulated growth rates in panel (a) agree with the theoretically predicted rates in panel (b), which were obtained through numerical  solution of the dispersion relation in (2.4) using standard iterative techniques for nonlinear equations, as well as those reported by Amano & Hoshino (2009).The maximum growth rate occurs at k x M e ≈ 1.2 and matches the corresponding 1D growth rate, since the 2D dispersion relation in (2.4) is equivalent to its 1D counterpart when θ = 0.By symmetry of the dispersion relation, the growth rates are reflectionally symmetric about the k x axis, i.e., the wave solution is identical for positive and negative k y (not shown here).Note that the analysis method adopted here effectively sums the contributions of positive and negative k x by virtue of the complex conjugate product in the definition of E ϕϕ .Conversely, the linear instability is a one-way instability where positive k x are most unstable.As k x increases, the growth rates decrease for axial x harmonics of the fastest-growing fundamental (k x ̸ = 0, k y = 0) along the k x axis, as well as their neighbouring modes.The subsequent stages of the instability development process are depicted in figure 2. Growth of these harmonics at a comparable rate to the fastest-growing fundamental is first observed in figure 2(a).At later times, accelerated growth in the forms of quasi-linear harmonic growth (see Rajawat & Sengupta 2017, and references therein) and a strong non-linear fill-in process at intermediate wavenumbers via a catch-up mechanism is then observed in figure 2(b).Here, primarily non-harmonic wavenumbers experience accelerated growth beyond the linear growth rate of the fastestgrowing fundamental mode.Eventually, no further modal growth occurs on average and the potential reaches a saturated turbulence state.
To further illustrate the transition between the linear, harmonic growth, and accelerated growth regimes of the instability growth, figures 3 and 4 plot the electrostatic potential energy modal coefficients, Êx Ê * x + Êy Ê * y , as functions of time spanning the time intervals discussed in figures 1 and 2, as well as beyond.Here, E x and E y respectively denote the axial and transverse electric field strengths.In each panel of figure 3, k y is held fixed while k x is varied, i.e., the growth of waves of different axial wavenumbers is directly compared.In each panel of figure 4, k x is held fixed while k y is varied, i.e., the growth of waves of different transverse wavenumbers is directly compared.
Modes first grow linearly at their analytically predicted modal growth rates, with high-k modes exhibiting slow growth and even decay in agreement with linear stability theory and in contrast to their faster-growing low-k counterparts, as marked by the left shaded region in each panel.Particularly, the growth rate of the fastest-growing mode, as indicated by the darkest bolded curve in figure 3(a), matches the theoretical maximum growth rate, as indicated by the red dashed line in the same panel.The curve corresponding to the same axial wavenumber is also bolded in figure 3(b), corroborating the observation in figure 1 that the axial wavenumber of the fastest-growing fundamental is insensitive to the transverse wavenumber.In further corroboration with the growth rates of figure 1, modes of various transverse wavenumbers exhibit larger linear growth at k x = 0.3 than at k x = 0, as shown in figure 4. The duration of the linear growth phase depends on the wavenumber magnitude.This is approximately visualised by the different time extents of the left shaded regions in figures 3 and 4, although the shaded regions are only intended as visual guides since the phase duration differs for each individual mode.
Subsequently, harmonics of the fastest-growing fundamental start to grow at a comparable (but slightly slower) rate to the fundamental, as marked by the middle shaded region in each panel.This process begins with the axial x harmonics [k x ̸ = 0, k y = 0, figure 3(a)]: as also evidenced by figure 2(a), purely axial modes whose x wavenumbers are approximately integer-valued multiples of the most rapidly growing fundamental lead the harmonic growth process.Similar dynamics occur for weakly oblique x harmonics [figure 3(b)], where the x wavenumbers are identically almost integer-valued multiples of the fastest-growing fundamental while the y wavenumbers are slightly non-zero.Figure 4 indicates that this is then followed by the purely transverse (k x = 0, k y ̸ = 0) and eventually strongly oblique (k x , k y ̸ = 0 harmonics where both wavenumber components have large magnitudes) resonant modes.The observed delay in this onset is about 100-200 dimensionless times in figure 4(a) and 300-500 dimensionless times in figure 4(b).Note that the k x = 0 purely transverse modes in figure 4(a), which experienced limited growth in the linear phase, are eventually excited to significant amplitudes in the harmonic growth phase.
The final accelerated growth phase is indicated approximately by the right shaded regions in figures 3 and 4. Here, the remaining non-harmonic modes, which include both axis-parallel and oblique wavevector orientations, experience accelerated super-linear growth, i.e., beyond that predicted for any mode by linear theory.This abrupt growth fills in the intermediate wavenumbers to form a saturated and persistent broadband spectrum.While the non-harmonic modes grow significantly in this phase to catch up to the harmonic and fundamental wave amplitudes in the broadband spectrum, the growth is much less pronounced for these latter wave categories.The time interval between the onset of harmonic growth and the development of saturated turbulence in the electric potential and field is about O(10 3 ) dimensionless times, or O(10) ion-acoustic periods.
The aforementioned growth of harmonics (locking) and subsequent fill-in (non-linear regime) are features of turbulence also seen in, e.g., hydrodynamic simulations of turbulent boundary layers, where energy cascades from long-wavelength modes to their short-wavelength harmonics through quasi-linear interactions, while non-linear triadic and polyadic interactions induce eventual fill-in of the intermediate-wavenumber spectral content (e.g., Sayadi et al. 2011).Parallel analysis of an ensemble of comparable 1D simulations (not shown here) suggests that the growth rates and order of growth of modes are insensitive to the details of the initial system noise keeping its magnitude constrained.As with hydrodynamic turbulence, this suggests that the system has limited memory of the initial conditions, as self-similarity eventually percolates across the system scales to achieve some form of statistical equilibrium.Note that longer modes (k ⩽ 1) mostly lead and exceed shorter modes (k ⩾ 1) in the modal growth process illustrated in figure 3 and particularly figure 4, possibly casting in doubt the existence of an inverse energy cascade [see, e.g., the discussion of causality in modal development by Lozano-Durán & Arranz (2022)].In other words, strong parallels with hydrodynamic turbulence are exhibited, where the dominant dynamics are associated with a forward energy cascade predominantly carrying energy from large to small scales.As such, electrostatic plasma turbulence of a sufficiently large separation of scales should also reasonably satisfy local isotropy, where small-and intermediate-scale statistics are mostly independent of largescale forcing details, and are thus weak functions of spatial direction in small spatial neighbourhoods.This is largely supported by the time-averaged electrostatic potential energy spectrum in its saturated turbulence state, as depicted in figure 5.The wave energy is maximum at the lower-k modes (i.e., long-wavelength modes) and decays toward highk modes.There exists a minor tendency for waves to travel at an angle of 35 • to the axial direction due to ion-induced wave scattering.When the electric field is oriented in this direction, its axial component is slightly larger than its transverse component, which is in corroboration with the initial instability being driven in the axial direction.In addition, this tendency is consistent with the analytical prediction in the review by Bychenkov et al. (1988), which discusses the role of ion-induced scattering of ion-acoustic waves on turbulent fluctuations in more detail.Apart from this minor asymmetry, the energy spectrum exhibits isotropy at intermediate scales on the whole.Such local isotropy is consistent with the principles of scale separation and self similarity underlying a turbulent field with largely uni-directional inter-scale transport from large to small scales, as introduced in the seminal works of Richardson (1922), Kolmogorov (1941), andOnsager (1945) in the context of the classical energy cascade in hydrodynamic turbulence [for generalisation beyond kinetic energy transport, see, e.g., Chan et al. (2021)].It should be emphasised that the current work considers unmagnetised instabilities where there is no bulk magnetic field introducing further anisotropy into the system.This study focuses on possible connections and similarities between hydrodynamic and electrostatic plasma turbulence as a starting point to invoke parallels in any cascading dynamics, and the effects of magnetic anisotropy are to be a subject of further study.

Turbulent saturation: ion velocity and energy distributions
Figure 6 plots the ion velocity distribution, weighted by the squared velocity magnitude, i.e., |v| 2 f i (v), to highlight the phase-space distribution of energy.The results for the baseline simulation introduced in § 3 are shown at three spatial locations for two time instances soon after the potential reaches its saturated turbulence state as observed in § 3.1.Shortly after potential saturation, figures 6(a)-(c) demonstrate highenergy ion formation in the direction of the net electron drift at t = 2.3 × 10 3 , which is past the accelerated growth stage identified in figures 3 and 4, with similarities to the 1D instability (Hara & Treece 2019).The energetic ions adopt speeds that exceed the initial ion-acoustic speed, which is approximately 3c i for the chosen Te / Ti , by under an order of magnitude.While these elevated speeds are suggestive of electron heating (see also figure 8), the ions remain strongly resonant with dominant ion-acoustic modes and exhibit coherence in physical space at these early times when not many ion-acoustic periods have elapsed since potential saturation.Conversely, at t = 3.8 × 10 3 , i.e., about 30 ion-acoustic periods later, figures 6(d)-(f) indicate that more multi-scale ion phase- accumulated response to the sustained broadband spectrum of ion-acoustic modes.Once again, the corresponding energetic ion speeds exceed the initial ion-acoustic speed and may be associated with concomitant electron heating.Such elevated energies can be of interest to downstream physical and engineering phenomena, such as the onset of hollow cathode sputtering in spacecraft thrusters.However, at the selected initial M e , more backward-streaming ions were observed in the 1D instability (Hara & Treece 2019) and appear to be less abundant in the 2D case.
In order to depict vortical structures in phase space more clearly, figure 7 plots the time evolution of the ion distribution function on spatial-velocity axes in the longitudinal f y,vy i (x, v x ) and transverse f x,vx i (y, v y ) directions in each row of panels, respectively, with averaging (•) performed over the remaining two phase-space axes (in superscript in the compact notation).The snapshots are extracted from the non-linear saturation phase across a duration of about 30 ion-acoustic periods.Figure 7(a) suggests that plasma waves initially propagate with net motion in the forward axial direction in agreement with 1D theory and observations (Hara & Treece 2019).Conversely, figure 7(b) shows that the wave energy along the transverse direction is initially weak, implying insignificant transverse ion trapping.As the saturated plasma waves evolve, the increasing emergence of overturning phase-space structures is indicative of gradually intensifying longitudinal particle trapping as observed in figure 7(c).The ensuing heating broadens the distribution tails in velocity space.At this time, ions travelling along the transverse direction also demonstrate gradually intensifying trapping in both the forward and backward directions, as indicated in figure 7(d), with comparable maximum energies along both the transverse and longitudinal directions.The observed structures correspond to waves of a travelling nature along the axial direction and of a standing nature along the transverse direction.As shown in figures 7(e) and (f), the increasing homogeneity in physical space is a signature of particle trajectories spanning an increasingly isotropic phase-space distribution.It should be noted that there remains a strong signature of ions travelling in the forward axial direction even at these late times, as shown in figure 7(e), while fast backward-moving ions are observable but not in excessive amounts.This is in contrast with corresponding 1D simulations, where a significant population of high-energy ions eventually travels in the direction opposite to the initial net electron drift (Hara & Treece 2019).The differences between 1D and 2D simulation results will be discussed in more detail in § 4 with reference to exchanges between various energy modes.
As the fidelity of the phase-space distributions can be further improved with increased velocity resolution due to the filamentous nature of the Vlasov equation, we focus on qualitative trends here and defer a more detailed analysis to future work.In line with this approach, figure 8 plots the time evolution of the spatially averaged ion f x,y i (v x , v y ) and electron f x,y e (v x , v y ) velocity distributions, as well as the electric potential (ϕ), all for the baseline simulation as well.Electron trapping and isotropisation occur more rapidly than their corresponding ion processes owing to the faster thermal speed and shorter response time of electrons.Particle trapping occurs at the phase speed of the excited plasma waves, which is on the order of the ion-acoustic speed k B Te / mi and becomes respectively O(10 −1 ) and O(10) on the non-dimensional electron and ion velocity axes for the subsequently heightened Te posited in the discussion of figure 6.The isotropisation of both velocity distributions is centred around these characteristic speeds in the positive axial direction, continuing the particle trapping process observed at early times in figure 7.However, note that while plasma waves are observed along both the transverse and longitudinal directions in figures 8(c,f,i), the wave amplitude is dampened at later times due to energy exchange with particles, and the plasma waves become less coherent.
Shortly after the potential reaches a saturated turbulence state as indicated in figure 8(c), the ion phase-space distribution remains largely Maxwellian with the nascent generation of high-energy forward-streaming ions in figure 8(a).The net electron drift can still be observed in figure 8(b) in corroboration with the nascent axial ion acceleration in figure 7(a).Then, the isotropisation of the spatially averaged electron phasespace distribution illustrated in figure 8(e) is accompanied by an increase in transversestreaming high-energy ions in figure 8(d) in corroboration with the heightened ion velocities in figure 7(e).As observed in figures 8(g,h) and later times (not shown here), this eventually results in the appearance of backward-streaming ions, as well as higherenergy forward-streaming ions, while the phase-space distributions themselves approach a quasi-steady saturated turbulence state.In summary, the development of multi-scale ion phase-space structures and accompanying cascades, if any, occurs after the potential reaches a saturated turbulence state and the electron phase-space distribution obeys isotropic statistics in the absence of any magnetic fields.It may then be reasonable to assume that multi-scale ion phase-space transport occurs in the background of a turbulent potential field with statistically isotropic electrons.The generation of the corresponding broadband ion phase-space structures occurs over a duration of approximately O(10 4 ) dimensionless times, or O(100) ion-acoustic periods-an order of magnitude larger than the time taken for the electric potential and field to reach a saturated turbulence state.Figure 9 plots the time evolution of the axial ion velocity distribution function, averaged over physical space and integrated over all transverse velocities, i.e., f x,y,vy i (v x ), for two M e values (M e = 2.3 and 2.8).The ion distributions with a larger initial M e deviate more rapidly from a Maxwellian and exhibit broader tails, particularly in the backward direction.These tails are evident five to six orders of magnitude lower than the peak distribution values, which necessitate similar variations in the particle-to-cell ratio for adequate resolution by particle solvers and showcase the relative superiority of the direct kinetic solver in this respect (Hara 2019).For both M e , the distributions are seen to approach an asymptotic state by the final plotted time instant, indicating statistical saturation in the phase-space distribution itself.An interesting observation is the transient appearance of a forward-streaming distribution plateau at about t = 4×10 3 that steepens at late times.The impermanence of the plateau is likely a consequence of the multi-dimensional nature of the system and can be explained as follows.Strong trapping primarily occurs in the forward-streaming direction at early times, resulting in a discernible distribution plateau similar to that observed in 1D simulations (Hara & Treece 2019, see also figure 13(a) of this manuscript).A subsequent tendency to statistical isotropy at late times then promotes trapping in more off-axis velocity directions at the expense of axial trapping.Additionally trapped ions can take oblique velocities with axial components smaller than the mean speed corresponding to the original plateau, thereby adding contributions to the left of the plateau and steepening it.The weakening of the plateau also implies continued damping of the plasma waves and is congruent with their gradually decreasing amplitude at late times [see also figure 8(i) and the time evolution of the electrostatic potential energy in figure 16 in appendix A].
It is interesting to note that the formation of backward-streaming ions is also present in the 2D case but with much smaller amplitude than was observed in its 1D counterpart (Hara & Treece 2019).In the latter, it was predicted and shown that M e ⩾ 1.3 results in a transition from the current-carrying ion-acoustic instability (leading to a unidirectional plasma wave) to the Buneman instability (leading to a bi-directional plasma wave in 1D).A key hypothesis of this study was that a similar transition to the Buneman instability would be observed in 2D.However, figure 9 indicates that the plasma wave amplitude is not large enough to excite enough ions with negative axial velocities.This is in large part due to the nature of energy transfer in multi-dimensional instabilities.

Comparisons between 1D and 2D instabilities
In order to relate key conclusions from previous 1D work to their multi-dimensional counterparts, it is instructive to compare several metrics, such as exchanges between different energy modes and macroscopic quantity variations, between the current 2D simulations and their 1D analogues.These metrics enable more insights into the temporal evolution of the ion and electron distribution functions, as well as their corresponding electric potential.Comparable 1D simulations of the same domain extents and resolutions as the baseline 2D simulations introduced in § 3 were performed to this end.We focus on the M e = 2.8 case in this section for brevity.

Energy exchange: transfer between different modes
We will analyse the transfer of energy between the modes defined in § 2.2 and reiterated here: the electrostatic potential energy 1 2 ε 0 | Ẽ| 2 d Ṽ , the bulk kinetic energy Figure 10 plots the time evolution of the exchange between different modes of energy for the larger initial electron Mach number considered in § 3.2, in a manner analogous to the energy decomposition considered by Chan & Boyd (2022a).For both dimensionalities, energy is initially transferred from electron bulk kinetic energy (due to the initial non-zero M e ) to the plasma waves, leading to an exponential growth of the electrostatic potential energy.The plasma waves then perturb the ion motion, increasing the ion bulk kinetic energy (particularly through ion trapping at positive axial velocities).Ion and electron heating are simultaneously observed as the plasma waves increase the ion and electron random kinetic energies in tandem.The primary net energy exchange transaction is from electron bulk kinetic energy to electron random kinetic energy with some minor ion heating.
The ion bulk kinetic energies evolve in a similar fashion in 1D and 2D.However, in 2D, the ion bulk kinetic energy consistently exceeds the electrostatic potential energy following plasma wave saturation, as shown in figure 10(a), while it only does so at late times in 1D, as shown in figure 10(b).The saturated electrostatic potential energy is much larger in 1D, indicating that the 1D instability results in plasma waves with larger amplitudes that eventually trap ions in both axial directions.The more oscillatory electrostatic potential energy curve in 1D may be attributed to this standing nature of the plasma waves, in contrast to the travelling nature of the axial 2D plasma waves noted earlier in figure 7.With fewer accessible degrees of freedom in 1D, the ion and electron random kinetic energies reach larger multiples of their initial values.This larger increase in effective ion and electron temperatures in the 1D instability is also driven by stronger electric fields and potential amplitudes as noted above.In the 2D case, the saturated electrostatic potential energy and electron bulk kinetic energy are concurrently lower, as the initial energy in the system is transferred to more degrees of freedom in the ion and electron random kinetic energies, and the magnitude of energy transfer is larger in absolute terms in the 2D instability.
Figure 11 plots the decomposition of the ion and electron kinetic energies into their axial and transverse components for the 2D instability.These may be written dimensionally as 1 2 ñ * m * Ũ 2 * ,x or y d Ṽ for the bulk energy components and for the random energy components.For the ion and electron bulk kinetic energies shown in figure 11(a), the axial energy dominates the transverse energy throughout the instability growth and saturation.Note from the transverse energy curves that anisotropy appears early in the system in the process of turbulent saturation, where specific modes dominate prior to the non-linear fill-in stage and the potential spectrum is not uniformly broadband.Some isotropisation in species motion then occurs at later times amid the isotropic turbulent potential field to reduce the mean transverse velocities.For the ion and electron random kinetic energies shown in figure 11(b), axial heating occurs before and more intensely than transverse heating.The two energy components quickly equilibrate in the case of the electrons, marking a rapid return to statistical isotropy.On the other hand, the ions clearly exhibit a reduced tendency to isotropy, which is indicative of lower trapping efficiency away from the forward axial direction.

Distribution isotropy and mean drifts
Figure 12(a) plots a comparison of the averaged axial electron velocity distribution functions, f x,y,vy e (v x ), for the 1D and 2D instabilities.The 1D distribution extends to larger speeds as an indication of a larger effective electron temperature, in corroboration with the discussion of figure 10.At the late times in the non-linear saturation regime considered here, the 1D distribution retains a peak around in the positive axial direction exhibiting the effects of the initial electron drift.This is absent in the 2D distribution, indicating that the initial electron drift energy is almost entirely transferred to the multidimensional plasma waves alongside electron heating.It should be emphasised that there is an additional degree of freedom for the deposition of random kinetic energy in the 2D case that constitutes strong heating in absolute terms.Owing to the multi-dimensional nature of the 2D instability, the phase-space distribution plateau is allowed to extend in multiple dimensions, resulting in the lower net electron axial drift seen in figure 12(b).
Figure 13(a) plots a similar comparison of the averaged axial ion velocity distribution functions, f x,y,vy i (v x ).Likewise, the 1D distribution extends to larger speeds as an indication of a larger effective ion temperature.In 1D, plateaus indicative of ion-acoustic wave behaviour are present at positive and negative axial velocities.A plateau is only observed in the positive axial direction in 2D.While the late-time ion velocity distribution functions in figure 8 suggested the presence of energetic backward-streaming ions in the 2D instability, these represent a smaller proportion of the ion population compared to those generated from a 1D instability.In other words, the multi-dimensional nature of the 2D instability appears to promote the earlier formation of transverse-streaming ions at the expense of backward-streaming ion formation due to significant early growth of transverse plasma waves (see also figures 3 and 4).Correspondingly, the mean axial ion velocity in the 1D simulations, where there are almost as many backward-streaming ions as forward-streaming ions, is smaller than that in the 2D simulations, where forwardstreaming ions outnumber backward-streaming ions, as seen in figure 13(b).The greater symmetry of the ion distribution along the axial velocity direction in the 1D instability in figure 13(a), which corresponds to a less pronounced axial ion drift in figure 13(b), implies that the 1D plasma waves take on more of a standing than a travelling character in contrast to 2D axial waves as observed earlier, such as in the vortical structures of the left column of figure 7.This accounts for the more oscillatory nature of some of the 1D curves remarked earlier in the discussion of figure 10.

Potential amplitude and trapping frequency
We take a closer look at the effects of the electrostatic potential energy magnitude on ion and electron trapping.Trapping intensity may be quantified in more detail through the bounce frequency g, which is proportional to the square root of the potential amplitude, ϕ rms , keeping all else constant.Comparisons of the root-mean-squared potential are plotted in figure 14(a).The 1D potential amplitude consistently exceeds its 2D counterpart by a factor of 4 or more, suggesting that trapping occurs at least twice as rapidly in the 1D instability than in the 2D one given that g ∼ √ ϕ rms .This may account for the lower degree of 2D energisation and the smaller tendency of the 2D instability to trap ions in the backward-streaming direction.As seen on a logarithmic scale in figure 14(b), the difference in electrostatic potential energies between the 1D and 2D instabilities is of a couple of orders of magnitude, and not simply because energy is partitioned between the axial and transverse degrees of freedom.The smaller potential amplitudes and field strengths in the 2D system account for smaller effective ion and electron temperatures.

Conclusions
The investigation of electrostatic current-driven plasma instabilities is crucial to determine the origin and fluxes of high-energy ions.Such ions have practical implications including hollow cathode erosion in electric spacecraft thrusters, as well as the inception of magnetic reconnection and cosmic rays.A direct kinetic solver is used to study the evolution and directional dependence of the ion and electron velocity distribution functions, as well as their moments and accompanying electric potential, without contamination from statistical noise inherent in state-of-the-art particle methods.
The electric potential field underlying such current-driven instabilities exhibits four developmental stages: linear modal growth, harmonic growth, accelerated growth via quasilinear mechanisms alongside non-linear fill-in, and saturated turbulence.The maximum linear modal growth rate matches analytical predictions from the linear plasma dispersion relation.Harmonics of the fastest-growing fundamental, followed by intermediatewavenumber modes, grow in a process that resembles the development of hydrodynamic turbulence.Ion and electron phase-space structures further exhibit a statistical tendency to isotropy also reminiscent of classical fluid behaviour.However, unlike hydrodynamic turbulence, which only fully emerges in 3D, such plasma turbulence occurs even in lowerdimensional instabilities as postulated by Buneman (1959) since ions and electrons are allowed to interpenetrate unlike fluids.
While the electrostatic potential energy quickly saturates and reaches a turbulent state, transverse-streaming and then backward-streaming ions are only formed after several ion trapping cycles.In other words, high-energy ions appear to be generated amid a turbulent potential field with statistically isotropic electrons.This formation process is accelerated and energised with a larger initial electron Mach number.A schematic summarising the instability growth and non-linear saturation process is shown in figure 15.
The additional transverse degree of freedom in the 2D instability has several physical implications.A larger forward axial ion drift is observed due to a preference for transversestreaming ions over backward-streaming ions, and the corresponding axial plasma waves are of a travelling rather than a standing nature.Steepening of a transient plateau occurs in the ion velocity distribution at the trapping speed as an indication of continued longtime plasma wave damping.The axial electron drift is smaller due to effective isotropisation of the electron distribution in the multi-dimensional velocity space.Despite the larger energy source resulting from this, smaller effective ion and electron temperatures are exhibited due to multi-dimensional energy redistribution, alongside smaller potential amplitudes and field strengths, which result in slower and weaker trapping as the system approaches a saturated turbulence state.On the whole, the degree of energisation in 2D appears to be lower than that in 1D, and the return to statistical isotropy of ion motion in phase space is correspondingly arrested, causing a less substantial presence of backward-streaming high-energy ions.
Under the warm-electron conditions studied in this work, it is shown that forwardstreaming ions of up to 50 eV, as well as transverse-streaming and backward-streaming ions of up to 20 eV, are readily generated by electrostatic current-carrying instabilities.A detailed characterisation of such high-energy ion formation and directional distribution, which constitute precursors of the hollow cathode erosion process, can improve predictions of spacecraft thruster lifetime and complement accelerated life testing.A complete characterisation, however, necessitates the inclusion of additional physics such as threedimensionality, collisions, and magnetisation, which will be addressed through further development and deployment of the direct kinetic method.

Figure 1 .
Figure1.Numerical growth rates of spectral modes {kxMe, kyMe} for the first quadrant of the potential spectrum E ϕϕ , which is obtained via a Fourier transform of ϕ in physical space into its Fourier coefficients φ, and then computation of φ φ * , which is the modal coefficient multiplied by its complex conjugate.The growth rates are plotted for t = [3.4× 10 1 , 2.7 × 10 2 ] in panel (a) and are obtained via linear regression of the spectral magnitudes in time.The maximum growth rate predicted by the dispersion relation in (2.4) is 0.017, and the analytical growth rates corresponding to other modes of interest are obtained through numerical solution of the dispersion relation in (2.4) and plotted in panel (b).Growth rates less than 5 × 10 −3 are excluded to remove cases where oscillations confound the regression.Hereinafter, potentials and wavenumbers are normalised by φth and 1/ λD, respectively, and times by 1/ωe.

Figure 3 .
Figure 3.Time evolution of the electrostatic potential energy spectrum, which is obtained via a Fourier transform of the axial and transverse field strengths in physical space, Ex and Ey, into their Fourier coefficients Êx and Êy, and then computation of the quantity Êx Ê * x + Êy Ê * y .The spectral coefficients are plotted for (a) ky = 0 and (b) ky = 0.3.The sloped dashed lines denote the maximum growth rate from (2.4).Every second mode is plotted (∆k plot = 2∆k = 2kmin = 4π/80 ≈ 0.16) and the curves are coloured from blue to yellow (dark to light in greyscale) in increasing k (kmax = Nxkmin/2 ≈ 7.9).The curves corresponding to the fundamental x wavenumber are bolded in panels (a) and (b), while those corresponding to the first two harmonics are also bolded in panel (a).The shaded regions denoting the instability stage are intended only as visual guides, as the linear, harmonic growth, and accelerated growth stages differ for each mode.Hereinafter, field strengths are normalised by φth / λD.

Figure 4 .
Figure 4. Time evolution of the electrostatic potential energy spectrum for (a) kx = 0 and (b) kx = 0.3.For a description of the dashed lines and colours, refer to the caption of figure 3.

Figure 8 .
Figure 8.Time evolution of spatially averaged (a,d,g) ion and (b,e,h) electron velocity distribution functions (VDFs), as well as (c,f,i) the electric potential field.The provided snapshots are (a-c) near the time instant where the potential reaches its saturated turbulence state, as well as about (d-f) 35 and (g-i) 120 ion-acoustic times after.

Figure 9 .
Figure 9.Time evolution of the axial ion velocity distribution function in the non-linear saturation phase, averaged over physical space and integrated over all transverse velocities, for (a) Me = 2.3 and (b) Me = 2.8.The first three time snapshots correspond to those in figure 8, while the final time instant is about 210 ion-acoustic times after the onset of potential turbulence in the first snapshot.

Figure 10 .
Figure 10.Time evolution of potential (PE) and various kinetic (KE) energy modes in the (a,c) 2D and (b,d) 1D instabilities for Me = 2.8.The energy densities are plotted on a logarithmic axis with a zoomed-in vertical axis for panels (c,d).Axial and transverse energy densities are summed for the 2D case.All energy densities are computed by dividing the relevant dimensional energies by the domain measure (L λD in 1D and L λD2 , and the random kinetic energy ñ * k B T * d Ṽ .

Figure 11 .
Figure 11.Time evolution of the axial and transverse kinetic energy (KE) modes for the (a) bulk and (b) random energies in the 2D2V Me = 2.8 case.Note the difference in vertical axis ranges in the two panels.For the energy normalisation, refer to the caption of figure 10.

Figure 12 .
Figure 12.(a) Comparison of 1D and 2D axial electron velocity distribution functions (VDFs), averaged over all points in space and t ∈ [7 × 10 3 , 9 × 10 3 ] in the non-linear saturation regime, and integrated over all transverse velocities for the 2D VDF.(b) Comparison of time evolution of 1D and 2D axial electron drifts.Both plots are for Me = 2.8.

Figure 13 .Figure 14 .
Figure 13.(a) Comparison of 1D and 2D axial ion velocity distribution functions (VDFs), averaged over all points in space and t ∈ [7 × 10 3 , 9 × 10 3 ] in the non-linear saturation regime, and integrated over all transverse velocities for the 2D VDF.(b) Comparison of time evolution of 1D and 2D axial ion drifts.Both plots are for Me = 2.8.

Figure 15 .
Figure15.Summary of instability growth and non-linear saturation process of various macroscopic and phase-space quantities as a function of the dimensionless time, t.Ion phase space isotropisation is arrested in the 2D instability.

Figure 16 .
Figure 16.Time evolution of the axial and transverse electrostatic potential energies for the (a) baseline and (b) velocity-refined simulations introduced in § 3 with Me = 2.3.The analytical lines denote the growth rate from linear stability theory via solution of (2.4).