2.1. Direct kinetic solver

The DK 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, Boyd & Kolobov Reference Hara, Boyd and Kolobov2012; Hara et al. Reference Hara, Sekerak, Boyd and Gallimore2014, Reference Hara, Chapman, Banks, Brunner, Joseph, Berger and Boyd2015, Reference Hara, Barth, Kaminski, Dodin and Fisch2017; Hara & Hanquist Reference Hara and Hanquist2018; Raisanen, Hara & Boyd Reference Raisanen, Hara and Boyd2019; Vazsonyi et al. Reference Vazsonyi, Hara and Boyd2020). In contrast to state-of-the-art particle-in-cell solvers, DK 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, $\tilde {f}_*$, for some particle type $*=i,e$ (ions, electrons) according to the Vlasov equation

(2.1)\begin{equation} \frac{\partial \tilde{f}_* \left(\boldsymbol{\tilde{x}},\boldsymbol{\tilde{v}}; \tilde{t}\right)}{\partial \tilde{t}} + \boldsymbol{\tilde{v}}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{ \tilde{x}}} \tilde{f}_*\left(\boldsymbol{\tilde{x}},\boldsymbol{\tilde{v}};\tilde{t}\right) + \frac{\tilde{q}_* \boldsymbol{\tilde{E}}}{\tilde{m}_*}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{\tilde{v}}} \tilde{f}_*\left(\boldsymbol{\tilde{x}},\boldsymbol{\tilde{v}};\tilde{t}\right) = 0, \end{equation}

where $\boldsymbol {\tilde {x}}$, $\boldsymbol {\tilde {v}}$ and $\boldsymbol {\tilde {E}}$ respectively denote the position, velocity and electric field vectors, $\tilde {q}_*$ denotes the charge of the simulated particle type and $\tilde {t}$ denotes the time. This is consistent with the conservation of charge in phase space. The computational domain is discretised in $\boldsymbol {\tilde {x}}\unicode{x2013}\boldsymbol {\tilde {v}}$ space with a parallelised second-order finite-volume method, which is described by Chan & Boyd (Reference Chan and Boyd2022a,Reference Chan and Boydb). Assuming small induced magnetic fields and their rates of change, we adopt the electrostatic approximation, where Gauss's law is expressed using $\boldsymbol {\tilde {E}}=-\boldsymbol {\nabla }_{\boldsymbol {\tilde {x}}}\tilde {\phi }$ as a Poisson equation for the electric potential $\tilde {\phi }$ of the form

(2.2)\begin{equation} \boldsymbol{\nabla}_{\boldsymbol{\tilde{x}}}^2 \tilde{\phi} ={-}\frac{e\left(\tilde{n}_i-\tilde{n}_e\right)}{\varepsilon_0}, \end{equation}

where $\varepsilon _0$ and $e$ respectively denote the vacuum permittivity and elementary charge, and $\tilde {n}_i$ and $\tilde {n}_e$ respectively denote the (singly charged) ion and electron number densities

(2.3a,b)\begin{equation} \tilde{n}_i\left(\boldsymbol{\tilde{x}};\tilde{t}\right) = \int_{\boldsymbol{\tilde{v}}} \tilde{f}_i\left(\boldsymbol{\tilde{x}},\boldsymbol{\tilde{v}'};\tilde{t}\right) \mathrm{d} \boldsymbol{\tilde{v}'}; \quad \tilde{n}_e\left(\boldsymbol{\tilde{x}};\tilde{t}\right) = \int_{\boldsymbol{\tilde{v}}} \tilde{f}_e\left(\boldsymbol{\tilde{x}},\boldsymbol{\tilde{v}'};\tilde{t}\right) \mathrm{d} \boldsymbol{\tilde{v}'}. \end{equation}

Periodic and no-flux boundary conditions are employed for $\boldsymbol {\tilde {x}}$ and $\boldsymbol {\tilde {v}}$, respectively.

2.2. 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 long-wavelength modes, we choose domain lengths sufficiently larger than the Debye length, $\tilde {\lambda }_{D} = \sqrt {(\varepsilon _0 k_{B} \tilde {T}_e)/(\tilde {n}_e e^2)}$. The initial species temperatures are $\tilde {T}_e = 2\ \text {eV}$ and $\tilde {T}_i = 0.2\ \text {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 = \tilde {U}_e/\tilde {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 Reference Forslund and Shonk1970; Bret & Dieckmann Reference Bret and Dieckmann2008; Chapman et al. Reference Chapman, Winjum, Berger, Dimits, Banks, Brunner, Joseph and Ghosh2021). Hereinafter, we define $x$ as the axial direction and $y$ as the transverse direction in two dimensions. 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 multi-dimensional electrostatic waves travelling at an arbitrary angle $\theta$ 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:

(2.4)\begin{equation} 1 + \sum_* \frac{\tilde{\omega}_*^2}{\tilde{c}_*^2\tilde{k}^2}\left[1 + \left(\frac{\tilde{\omega}}{\sqrt{2}\tilde{c}_* \tilde{k}} - \frac{ M_*\cos\theta}{\sqrt{2}}\right) Z\left(\frac{\tilde{\omega}}{\sqrt{2}\tilde{c}_* \tilde{k}} - \frac{ M_*\cos\theta}{\sqrt{2}}\right) \right] = 0,\end{equation}

where $\tilde {k} = |\boldsymbol {\tilde {k}}|$ and $\tilde {\omega }$ are, respectively, the modal wavenumber magnitude and angular frequency of the electrostatic wave in question, $\tilde {\omega }_* = \sqrt {(\tilde {n}_* e^2)/(\tilde {m}_* \varepsilon _0)}$ is the plasma frequency, $\tilde {c}_* = \sqrt {k_{B} \tilde {T}_*/\tilde {m}_*}$ is the species thermal speed and $Z$ is the plasma dispersion function corresponding to a Maxwellian velocity distribution

(2.5a,b)\begin{equation} Z(\xi) = \frac{1}{\sqrt{\rm \pi}}\int_{-\infty}^{\infty} \frac{{\rm e}^{{-}z^2}}{z - \xi} \, \mathrm{d} z = \mathrm{i} \sqrt{\rm \pi} {\rm e}^{-\xi^2} \,\text{erfc}(-\mathrm{i} \xi); \quad \frac{\mathrm{d} Z(\xi)}{\mathrm{d} \xi} ={-}2\left[1 + \xi Z(\xi)\right]. \end{equation}

Unless otherwise stated, lengths, velocities and time are hereinafter respectively non-dimensionalised by $\tilde {\lambda }_{D}$, $\tilde {c}_*$ and the inverse electron frequency $1/\tilde {\omega }_e = \tilde {\lambda }_{D}/\tilde {c}_e$. Specifically, $k = \tilde {k}\tilde {\lambda }_{D}$, $\boldsymbol {x} = \boldsymbol {\tilde {x}}/\tilde {\lambda }_{D}$, $f_* = \tilde {f}_*\tilde {c}_*$, $\boldsymbol {v} = \boldsymbol {\tilde {v}}/\tilde {c}_*$ and $t = \tilde {\omega }_e\tilde {t}$. Note that the characteristic ion time scale, $1/\tilde {\omega }_i$, and correspondingly the ion-acoustic period, exceed the characteristic electron time scale, $1/\tilde {\omega }_e$, by a factor of $\sqrt {\tilde {m}_i/\tilde {m}_e} \approx 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 $\int \frac {1}{2} \varepsilon _0 |\boldsymbol {\tilde {E}}|^2 \, \mathrm {d} \tilde {V}$, the bulk kinetic energy $\int \frac {1}{2} \tilde {n}_* \tilde {m}_* \tilde {U}_*^2 \, \mathrm {d} \tilde {V}$ and the random kinetic energy $\int \tilde {n}_* k_{B} \tilde {T}_* \, \mathrm {d} \tilde {V}$. In particular, we consider the kinetic energies in detail in § 4.

3.1. 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_{\phi \phi } = \hat {\phi }\hat {\phi }^*$, 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 (Reference Amano and Hoshino2009). The maximum growth rate occurs at $k_x M_e \approx 1.2$ and matches the corresponding 1-D growth rate, since the 2-D dispersion relation in (2.4) is equivalent to its 1-D counterpart when $\theta = 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_{\phi \phi }$. 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 \neq 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 (Reference Rajawat and Sengupta2017), and references therein) and a strong nonlinear 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 fastest-growing fundamental mode. Eventually, no further modal growth occurs on average and the potential reaches a saturated turbulence state.

Figure 1. Numerical growth rates of spectral modes $\{k_x M_e,k_y M_e\}$ for the first quadrant of the potential spectrum $E_{\phi \phi }$, which is obtained via a Fourier transform of $\phi$ in physical space into its Fourier coefficients $\hat {\phi }$, and then computation of $\hat {\phi }\hat {\phi }^*$, which is the modal coefficient multiplied by its complex conjugate. The growth rates are plotted for $t=[3.4\times 10^1,2.7\times 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\times 10^{-3}$ are excluded to remove cases where oscillations confound the regression. Hereinafter, potentials and wavenumbers are normalised by $\tilde {\phi }_{\text {th}}$ and $1/\tilde {\lambda }_{D}$, respectively, and times by $1/\tilde {\omega }_e$.

Figure 2. The same growth rates plotted in figure 1 for (a) $t=[4.6\times 10^2,1.1\times 10^3]$ and (b) $t=[1.1\times 10^3,1.6\times 10^3]$. Note the different colour bar maximum in panel (b). The extent of plotted wavenumbers approaches $k_x, k_y \approx 2{\rm \pi}$, which corresponds to the Debye length itself.

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, $\hat {E}_x\hat {E}_x^* + \hat {E}_y\hat {E}_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.

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, $E_x$ and $E_y$, into their Fourier coefficients $\hat {E}_x$ and $\hat {E}_y$, and then computation of the quantity $\hat {E}_x\hat {E}_x^* + \hat {E}_y\hat {E}_y^*$. The spectral coefficients are plotted for (a) $k_y = 0$ and (b) $k_y = 0.3$. The sloped dashed lines denote the maximum growth rate from (2.4). Every second mode is plotted $(\Delta k_\text {plot} = 2\Delta k = 2 k_\text {min} = 4{\rm \pi} /80 \approx 0.16)$ and the curves are coloured from blue to yellow (dark to light in greyscale) in increasing $k$ $(k_\text {max} = N_x k_\text {min}/2 \approx 7.9)$. The curves corresponding to the fundamental $x$ wavenumber are bolded in panels (a,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 $\tilde {\phi }_{\text {th}}/\tilde {\lambda }_{D}$.

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

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 \neq 0, k_y = 0$, figure 3a); 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 occurs for weakly oblique $x$ harmonics (figure 3b), 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 \neq 0)$ and eventually strongly oblique $(k_x, k_y \neq 0\text { 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 approximately $O(10^3)$ dimensionless times, or $O(10)$ ion-acoustic periods.

The aforementioned growth of harmonics (locking) and subsequent fill-in (nonlinear 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 nonlinear triadic and polyadic interactions induce eventual fill-in of the intermediate-wavenumber spectral content (e.g. Sayadi, Hamman & Moin Reference Sayadi, Hamman and Moin2011). Parallel analysis of an ensemble of comparable 1-D 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 \leqslant 1$) mostly lead and exceed shorter modes ($k \geqslant 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 Reference Lozano-Durán and Arranz2022). In other words, strong parallels with hydrodynamic turbulence are exhibited, where the dominant dynamics is 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 large-scale 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 high-$k$ modes. There exists a minor tendency for waves to travel at an angle of $35^\circ$ 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, Silin & Uryupin (Reference Bychenkov, Silin and Uryupin1988), 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 (Reference Richardson1922), Kolmogorov (Reference Kolmogorov1941) and Onsager (Reference Onsager1945) in the context of the classical energy cascade in hydrodynamic turbulence (for generalisation beyond kinetic energy transport, see, e.g. Chan, Johnson & Moin Reference Chan, Johnson and Moin2021). 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.

Figure 5. Electrostatic potential energy spectrum in the nonlinear saturation phase averaged over the time interval $t \in [2.9 \times 10^3,3.7 \times 10^3]$. The dashed lines mark $35^\circ$ angles from the horizontal.

3.2. Turbulent saturation: ion velocity and energy distributions

Figure 6 plots the ion velocity distribution, weighted by the squared velocity magnitude, i.e. $|\boldsymbol {v}|^2 f_i(\boldsymbol {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)–6(c) demonstrate high-energy ion formation in the direction of the net electron drift at $t=2.3\times 10^3$, which is past the accelerated growth stage identified in figures 3 and 4, with similarities to the 1-D instability (Hara & Treece Reference Hara and Treece2019). The energetic ions adopt speeds that exceed the initial ion-acoustic speed, which is approximately $3\tilde {c}_i$ for the chosen $\tilde {T}_e/\tilde {T}_i$, 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\times 10^3$, i.e. approximately 30 ion-acoustic periods later, figures 6(d)–6( f) indicate that more multi-scale ion phase-space structures, i.e. with a broad range of characteristic scales, are observed. High-energy ions are generated in significant amounts, first in the forward ($+x$) direction, and subsequently in the transverse ($\pm y$) and backward ($-x$) directions due to the multi-dimensional plasma wave formation, at equivalent temperatures of between 20 and 50 eV. The velocity distributions now demonstrate greater incoherence in physical space due to 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 1-D instability (Hara & Treece Reference Hara and Treece2019) and appear to be less abundant in the 2-D case.

Figure 6. Local ion VDFs, weighted by the squared velocity magnitude to highlight the phase-space distribution of energy, at three different locations at (a–c) $t = 2.3\times 10^3$ and (d–f) $t = 3.8\times 10^3$. The contours are logarithmically spaced with their corresponding exponents labelled in the colour bar and the two concentric circles represent 20 and 45 eV contours. Hereinafter, velocities are normalised by the species thermal speed $\tilde {c}_*$.

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 $[\bar {f}_i^{y,v_y}(x,v_x)]$ and transverse $[\bar {f}_i^{x,v_x}(\kern0.7pt y,v_y)]$ directions in each row of panels, respectively, with averaging $(\bar {\cdot })$ performed over the remaining two phase-space axes (in superscript in the compact notation). The snapshots are extracted from the nonlinear saturation phase across a duration of approximately 30 ion-acoustic periods. Figure 7(a) suggests that plasma waves initially propagate with net motion in the forward axial direction in agreement with 1-D theory and observations (Hara & Treece Reference Hara and Treece2019). 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 7( 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 1-D simulations, where a significant population of high-energy ions eventually travels in the direction opposite to the initial net electron drift (Hara & Treece Reference Hara and Treece2019). The differences between 1-D and 2-D simulation results will be discussed in more detail in § 4 with reference to exchanges between various energy modes.

Figure 7. Ion distribution functions on spatial–velocity axes in the (a,c,e) longitudinal ($x, v_x$) and (b,d, f) transverse ($y, v_y$) directions, after averaging over the remaining two phase-space axes, at (a,b) $t = 2.3\times 10^3$, (c,d) $t = 3.0\times 10^3$ and (e, f) $t = 3.8\times 10^3$.

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 $[\bar {f}_i^{x,y}(v_x,v_y)]$ and electron $[\bar {f}_e^{x,y}(v_x,v_y)]$ velocity distributions, as well as the electric potential ($\phi$), 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 of the order of the ion-acoustic speed $\sqrt {k_{B}\tilde {T}_e/\tilde {m}_i}$ and becomes respectively $O(10^{-1})$ and $O(10)$ on the non-dimensional electron and ion velocity axes for the subsequently heightened $\tilde {T}_e$ 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 figure 8(c, f,i), the wave amplitude is damped at later times due to energy exchange with particles, and the plasma waves become less coherent.

Figure 8. Time evolution of spatially averaged (a,d,g) ion and (b,e,h) electron 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 approximately (d–f) 35 and (g–i) 120 ion-acoustic times after.

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 phase-space distribution illustrated in figure 8(e) is accompanied by an increase in transverse-streaming high-energy ions in figure 8(d) in corroboration with the heightened ion velocities in figure 7(e). As observed in figure 8(g,h) and later times (not shown here), this eventually results in the appearance of backward-streaming ions, as well as higher-energy 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. $\bar {f}_i^{x,y,v_y}(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 DK solver in this respect (Hara Reference Hara2019). 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 approximately $t = 4\times 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 1-D simulations (Hara & Treece Reference Hara and Treece2019, 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 the Appendix).

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

It is interesting to note that the formation of backward-streaming ions is also present in the 2-D case but with much smaller amplitude than was observed in its 1-D counterpart (Hara & Treece Reference Hara and Treece2019). In the latter, it was predicted and shown that $M_e \geqslant 1.3$ results in a transition from the current-carrying ion-acoustic instability (leading to a uni-directional plasma wave) to the Buneman instability (leading to a bi-directional plasma wave in one dimension). A key hypothesis of this study was that a similar transition to the Buneman instability would be observed in two dimensions. 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.

4.1. 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 $\int \frac {1}{2} \varepsilon _0 |\boldsymbol {\tilde {E}}|^2 \, \mathrm {d} \tilde {V}$, the bulk kinetic energy $\int \frac {1}{2} \tilde {n}_* \tilde {m}_* \tilde {U}_*^2 \, \mathrm {d} \tilde {V}$ and the random kinetic energy $\int \tilde {n}_* k_{B} \tilde {T}_* \, \mathrm {d} \tilde {V}$.

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 (Reference Chan and Boyd2022a). 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.

Figure 10. Time evolution of potential (PE) and various kinetic (KE) energy modes in the (a,c) 2-D and (b,d) 1-D instabilities for $M_e = 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 2-D case. All energy densities are computed by dividing the relevant dimensional energies by the domain measure ($L\tilde {\lambda }_{D}$ in one dimension and $[L\tilde {\lambda }_{D}]^2$ in two dimensions) and then normalising by the initial electron random kinetic energy density for a single spatial degree of freedom, $\tilde {n}_e k_{B} \tilde {T}_e/2$.

The ion bulk kinetic energies evolve in a similar fashion in one dimension and two dimensions. However, in two dimensions, 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 one dimension, as shown in figure 10(b). The saturated electrostatic potential energy is much larger in one dimension, indicating that the 1-D instability results in plasma waves with larger amplitudes that eventually trap ions in both axial directions. The more oscillatory electrostatic potential energy curve in one dimension may be attributed to this standing nature of the plasma waves, in contrast to the travelling nature of the axial 2-D plasma waves noted earlier in figure 7. With fewer accessible degrees of freedom in one dimension, 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 1-D instability is also driven by stronger electric fields and potential amplitudes as noted above. In the 2-D 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 2-D instability.

Figure 11 plots the decomposition of the ion and electron kinetic energies into their axial and transverse components for the 2-D instability. These may be written dimensionally as

(4.1)\begin{equation} \int \tfrac{1}{2} \tilde{n}_* \tilde{m}_* \tilde{U}_{*,x\text{ or }y}^2 \, \mathrm{d} \tilde{V}, \end{equation}

for the bulk energy components and

(4.2)\begin{equation} \int \tfrac{1}{2} \tilde{n}_* k_{B} \tilde{T}_{*,x\text{ or }y} \, \mathrm{d} \tilde{V}, \end{equation}

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 nonlinear 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.

Figure 11. Time evolution of the axial and transverse kinetic energy (KE) modes for the (a) bulk and (b) random energies in the 2D2V $M_e = 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.

4.2. Distribution isotropy and mean drifts

Figure 12(a) plots a comparison of the averaged axial electron velocity distribution functions, $\bar {f}_e^{x,y,v_y}(v_x)$, for the 1-D and 2-D instabilities. The 1-D 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 nonlinear saturation regime considered here, the 1-D distribution retains a peak around in the positive axial direction exhibiting the effects of the initial electron drift. This is absent in the 2-D distribution, indicating that the initial electron drift energy is almost entirely transferred to the multi-dimensional 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 2-D case that constitutes strong heating in absolute terms. Owing to the multi-dimensional nature of the 2-D 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 12. (a) Comparison of 1-D and 2-D axial electron VDFs, averaged over all points in space and $t \in [7\times 10^3,9\times 10^3]$ in the nonlinear saturation regime, and integrated over all transverse velocities for the 2-D VDF. (b) Comparison of time evolution of 1-D and 2-D axial electron drifts. Both plots are for $M_e = 2.8$.

Figure 13(a) plots a similar comparison of the averaged axial ion velocity distribution functions, $\bar {f}_i^{x,y,v_y}(v_x)$. Likewise, the 1-D distribution extends to larger speeds as an indication of a larger effective ion temperature. In one dimension, 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 two dimensions. While the late-time ion velocity distribution functions in figure 8 suggested the presence of energetic backward-streaming ions in the 2-D instability, these represent a smaller proportion of the ion population compared with those generated from a 1-D instability. In other words, the multi-dimensional nature of the 2-D 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 1-D simulations, where there are almost as many backward-streaming ions as forward-streaming ions, is smaller than that in the 2-D simulations, where forward-streaming 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 1-D instability in figure 13(a), which corresponds to a less pronounced axial ion drift in figure 13(b), implies that the 1-D plasma waves take on more of a standing than a travelling character in contrast to 2-D 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 1-D curves remarked earlier in the discussion of figure 10.

Figure 13. (a) Comparison of 1-D and 2-D axial ion VDFs, averaged over all points in space and $t \in [7\times 10^3,9\times 10^3]$ in the nonlinear saturation regime, and integrated over all transverse velocities for the 2-D VDF. (b) Comparison of time evolution of 1-D and 2-D axial ion drifts. Both plots are for $M_e = 2.8$.

4.3. 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, $\phi _\text {rms}$, keeping all else constant. Comparisons of the root-mean-squared potential are plotted in figure 14(a). The 1-D potential amplitude consistently exceeds its 2-D counterpart by a factor of 4 or more, suggesting that trapping occurs at least twice as rapidly in the 1-D instability than in the 2-D one given that $g \sim \sqrt {\phi _\text {rms}}$. This may account for the lower degree of 2-D energisation and the smaller tendency of the 2-D 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 1-D and 2-D 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 2-D system account for smaller effective ion and electron temperatures.

Figure 14. (a) Comparison of time evolution of 1-D and 2-D root-mean-squared potentials. (b) Comparison of time evolution of 1-D and 2-D electrostatic potential energies. Both plots are for $M_e = 2.8$.