Nonlinear current-driven ion-acoustic instability driven by phase-space structures

The model describes the collisionless evolution of a two-species, 1D electrostatic plasma. In addition to academic interest, a 1D model is relevant for plasma immersed in a strong, relatively homogeneous magnetic field [36]. The evolution of each particle distribution, f_s(x, v, t), where s = i, e, is given by the Vlasov equation,

$$\begin{equation} {\frac{\partial {f_s}}{\partial {t}}} \,+\, v \, {\frac{\partial {f_s}}{\partial {x}}} \,+\, \frac{q_s E}{m_s} \, {\frac{\partial {f_s}}{\partial {v}}} \;=\; 0, \end{equation} \tag{ 1 }$$

where q_s and m_s are the particle charge and mass, respectively. The evolution of the electric field E satisfies a current equation,

$$\begin{equation} \frac{\partial E}{\partial t} \;=\; - \sum _s \frac{m_s \omega _{p_s}^2}{n_0 q_s} \, \int v f_s(x,v,t) \, {\rm d}v, \end{equation} \tag{ 2 }$$

where ω_p,s is the plasma frequency and n₀ is the spatially averaged plasma density. The initial electric field is given by solving Poisson's equation. We denote δf_s ≡ f_s − f_0,s and $\tilde{f_s}\equiv f_s-\bar{f_s}$, where $f_{0,s}(v)=\bar{f_s}(v,0)$ is the initial velocity distribution, and $\bar{f_s}(v,t)$ is the spatial average of f_s. In a 1D periodic system, a spatially uniform current drives a uniform electric field, which oscillates at a frequency ω_u = ω_p,e(1 + m_e/m_i)^1/2. This rapid oscillation of both the uniform electric field and the uniform current is of little interest here [37]. Numerically, the average part of E is set to zero, following common practice [10, 38, 39].

On the one hand, in [6, 40], we described, verified, validated and benchmarked a semi-Lagrangian kinetic code COBBLES, capable of long-time simulations of 1D plasmas. A particular feature of this code is that the surface elements of phase-space density are locally conserved, up to the machine precision. Hereafter, we refer to these simulations as Vlasov simulations. On the other hand, to reproduce the results of [3, 10], we developed a simple particle-in-cell code, PICKLES (Particle-In-Cell Kinetic Lazy Electrostatic Solver), which can be switched between full-f and δf treatments. Hereafter, we refer to these simulations as full-f PIC simulations and δf PIC simulations, respectively. We denote the number of marker-particles per species as N_p. The PIC code is used only in section 3.2. The rest of this paper is based on Vlasov simulations.

Hereafter, we adopt the physical parameters of [3, 10]. The mass ratio is m_i/m_e = 4 (small mass ratio improves numerical tractability and the readability of phase-space contour plots). The system size is L = 2π/k₁, where $k_1 = 0.2 \lambda_D ^{-1}$. The initial velocity distribution for each species is a Gaussian, $f_{0,s}(v)=n_0/[(2\pi)^{1/2}v_{T,s}]\exp [-(v-v_{0,s})^2/2 v_{T,s}^2]$, with v_0,i = 0. The ion and electron temperatures are equal. Boundary conditions are periodic in real space. In COBBLES, we ensure zero-particle flux at the velocity cut-offs v_cut,s. We choose v_cut,s = v_0,s ± 7v_{T, s}.

All Vlasov simulations are performed with at least N_x = 768 and N_v = 1024 grid points in configuration-space and velocity-space, respectively, and with a time-step width at most $\Delta t = 0.1 \omega_{\rm p,e}^{-1}$. The grid cell size in real space is Δx = 0.04λ_D. Although the length-scales of interest are larger than the Debye length, such a small cell size is necessary to reduce numerical artifacts.

We run a series of Vlasov simulations for different values of initial drift v_d ≡ v_0,e − v_0,i. The linear stability threshold of the ion-acoustic mode is v_d = v_d,cr, with v_d,cr/v_T,i = 3.92. The initial velocity distributions are shown in figure 1 for v_d/v_T,i = 3.8. The initial perturbation is an ensemble of waves in the electron distribution,

$$\begin{equation} \left. f_{\rm e}\right|_{t=0} \;=\; \left[1+\sum_{m=1}^{m_{\rm max}} \frac{k_m}{k_1} \epsilon \cos (k_m x + \phi _m) \right] f_{0,{\rm e}}(v), \end{equation} \tag{ 3 }$$

where k_m = m k₁, m_max = 20 $(k_{\rm max}=4.0 \lambda _D^{-1})$, φ_m are random phases, and controls the initial electric field amplitude. We measure the electric field energy by the mean square potential, φ = 〈φ²〉^1/2. We define the potential energy as eφ.

Zoom In Zoom Out Reset image size

Figure 2 shows the time-evolution of potential energy, normalized to thermal energy, eφ/T, for three linearly stable cases, and two linearly unstable cases. For all stable cases, we hide data for times longer than half a numerical recurrence time, T_R/2 = π N_v/(10k_maxv_T,e). The oscillations for v_d/v_T,i < 4 are due to the beating of waves with opposite phase velocities. We observe that the solutions are consistent with linear theory. Note that for drift velocities of 4.2 and 4.5, the growth rate increases before saturation. This is in contrast with the conventional wave saturation, where the growth rate decreases in time. This phenomenon was predicted and confirmed in [7]. It appears for barely unstable cases, above a threshold amplitude, for which the nonlinear growth rate overcomes the linear growth rate.

Zoom In Zoom Out Reset image size

The important conclusion here, is that we observe no subcritical instability. This is in contradiction with [3], where subcritical instabilities are reported for the same parameters and v_d/v_T,i ⩾ 1.5, even for low-amplitude initial fluctuations with eφ/T ≪ 10⁻². We scanned the parameter space of velocity drift and initial amplitude, and concluded that subcritical instabilities emerge only when the drift is very close to linear marginal stability, and the initial perturbation is relatively large. Figure 3 shows the time-evolution of normalized potential energy, for v_d/v_T,i = 3.8, which is only 3% below marginal stability, and for large initial amplitudes. We do observe a subcritical instability when the initial amplitude is eφ/T ≈ 0.2, but not below. For v_d/v_T,i = 0.5, 1.5, 2.5 or 3.0, we did not obtain any subcritical instability even for initial amplitudes as high as eφ/T ≈ 2. These results are summarized in table 1, which shows whether an incoherent ensemble of waves is nonlinearly stable (W↓) or unstable (W↑). From this table, we conclude that the nonlinear stability threshold is approximately v_d/v_T,i = 3.5 (v_d/v_d,cr = 0.89).

Zoom In Zoom Out Reset image size

Subcritical instabilities are in many aspects qualitatively similar to linearly unstable cases. The saturation level of potential energy is similar, we observe wide particle redistribution in phase-space, especially of the electrons, turbulent heating, and significant anomalous resistivity. Figure 1 includes snapshots of ion and electron velocity distributions after the nonlinear growth for the subcritical simulation (v_d/v_T,i = 3.8) with high initial amplitude (eφ/T ≈ 0.2). At saturation, the electron distribution is flattened over a large range, −1 < v/v_T,i < 6. The ion distribution develops a plateau around v/v_T,i = 4, which is due to accumulating ion phase-space vortices.

Phase-space redistribution is associated with anomalous resistivity. We define the anomalous resistivity η as

$$\begin{eqnarray} &&\fl n_0 q_{\rm i}^2 \left( p_{\rm i}-p_{\rm e} \right) \eta = q_{\rm i} \left\langle E \right\rangle - \left(\frac{1}{m_{\rm i}}+\frac{1}{m_{\rm e}}\right)^{-1} \frac{{\rmd} \left( p_{\rm i}-p_{\rm e}\right)}{{\rmd} t} ,\nonumber\\ \end{eqnarray} \tag{ 4 }$$

where 〈E〉 is the spatial average of the electric field (here it is zero), and p_s ≡ ∫vf_s dx dv is the momentum of species s. For typical tokamak plasma conditions, n₀ = 10¹⁹ m⁻³ and T_i = T_e = 1 keV, the ion-electron collision frequency is of the order of ν_ei ∼ 10⁻⁷ω_p,e (with real mass ratio). Figure 4(a) shows the moving average (over $\delta t=16 \omega_{\rm p,e}^{-1})$ of η/η_coll where $\eta_{\rm coll}= m_{\rm e} \nu_{\rm ei}/(n_0 q_{\rm i}^2)$ is a typical value of the collisional resistivity. The maximum anomalous resistivity is 4 orders of magnitude higher than typical collisional resistivity, for both subcritical and supercritical cases. In addition, ion-acoustic waves cause both ion and electron heating. We define the mean thermal energy T_s ≡ (m_s/2 n₀)∫(v − p_s)²f_s dx dv. It reduces to the temperature for spatially uniform, Boltzmann distributions. Figure 4(b) shows the moving average of the mean thermal energy perturbation δT_s = T_s(t) − T_s(0). Both ion and electron thermal energies roughly double, for both subcritical and supercritical cases. These results indicate that the saturated level of turbulence, the anomalous resistivity, the turbulent heating, etc. are not directly affected by linear stability. In other words, essential nonlinear phenomena do not undergo any bifurcation at the linear stability threshold.

Zoom In Zoom Out Reset image size

The simulations of [3], which are in disagreement with the above Vlasov simulations, were obtained with a full-f PIC code, with a number of particles limited by the computing power of the time. We reproduced those simulations with N_p = 102 400, and indeed recovered the results of that reference, that is, subcritical growth for v_d/v_T,i ⩾ 1.5. Figure 5 includes one example of subcritical instability, for v_d/v_T,i = 2.5 and initial amplitude eφ/T ≈ 0.03. Figure 6 is a snapshot of both ion and electron phase-spaces, at ω_p,et = 400. We observe several electron holes, most noticeably at (k₁x, v/v_T,i) = (0.9π, 3.0). We checked that particles follow trapped orbits within the latter hole.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We ran a series of such full-f PIC simulations for different values of initial drift. We define the nonlinear growth rate as $\gamma _{\rm NL} = \partial \ln \tilde{f}_{\rm h,e} /\partial t$, where $\tilde{f}_{\rm h,e}(t) \equiv\max_{x,v}{|\tilde{f_{\rm e}}|}$ is the depth of the deepest negative perturbation. Taking an alternative definition of the nonlinear growth rate as the growth rate of the potential energy, yields similar growth rates. This is because the potential is dominated by the largest hole, owing to the relation $\varphi \sim \Delta v ^2 \sim \tilde{f}_{\rm h,e} ^2$ (see equation (7)). Since the growth rate depends on the amplitude, we must measure it at some fixed value of $\tilde{f}_{\rm h,e}$. We choose the same value as the reference, $\tilde{f}_{\rm h,e}\approx 0.1 n_0/v_{T,{\rm e}}$. Figure 7 shows the measured nonlinear growth rate as a function of initial drift. The data points are in agreement with simulation results of [3] (see figure 4). Besides, they are in qualitative agreement with the nonlinear instability theory of [41]. The theoretical growth rate was obtained by solving equation (99) (with the plus sign) in [41], while assuming a structure velocity v₊ = v_T,i (the holes have a wide range of different speeds, but v₊ = v_T,i is where the hole growth is expected to be the strongest), a self-binding factor b = 3 [42], a ratio $c=\tau_{\rm cl}^s/\tau^s=2.9$ [41], and generalized electron and ion structure lifetime as $\tau_{\rm e}=40 \omega_{\rm p,e}^{-1}$ and τ_i = 2τ_e [3]. This theory is only meant to give an estimate of the qualitative behavior. Here we only note that the simulation data happen to agree qualitatively with the theory for this particular set of mass ratio, temperature ratio and perturbation amplitude. A careful validation of the theory is out of the scope of this paper.

Zoom In Zoom Out Reset image size

Also shown is the linear growth rate for the most unstable wave-number, taking into account that the system size allows wave-numbers that are multiples of $0.2 \lambda _D^{-1}$ only. Given the agreement with theory, one is tempted to conclude that the observed instability is physical, rather than a numerical phenomenon, and that it is indeed the nonlinear ion-acoustic instability. However, these subcritical instabilities disappear (except for large initial amplitude and close to linear marginality) when the number of particles is increased or when the δf approach is adopted. Figure 5 includes time-series of potential energy in a full-f PIC simulation with N_p = 2²⁴, and in a δf PIC simulation with N_p = 102 400 with a similar initial amplitude (actually thrice larger at t = 0, but similar at ω_p,et = 30). Both cases are stable, indicating that growth in the simulations of [3] is due to numerical noise. Indeed, we found that the initial, unperturbed velocity distribution is so noisy that it is linearly unstable.

With a number of particles as small as 102 400, the velocity distributions in full-f PIC simulations are very noisy, and only approximately Gaussian. Figure 8 shows, for v_d/v_T,i = 2.5, the initial ion and electron distribution functions, which were obtained by distributing the particles into 128 boxes. Substituting the distribution at N_p = 102 400 into the linearized model equations, and solving the corresponding eigenvalue problem, yields a positive linear growth rate, γ_L/ω_p,e = 0.008. Varying the number of boxes from 32 to 512, or modeling the distribution by a spline and increasing the resolution of the discretization yielded similar values, γ_L/ω_p,e = 0.006–0.009, positive in any case. We therefore conclude that the instabilities observed in [3] were not subcritical, but linearly unstable, even when the drift was below the linear threshold for Gaussian distribution.

Zoom In Zoom Out Reset image size

The agreement between noisy simulations and theory can be explained as follows. The noisy distributions are such that, in small regions of velocity space, ∂_vf_0,e/m_e + ∂_vf_0,i/m_i > 0. This enables waves to grow linearly from small amplitudes eφ/T ≪ 10⁻², which leads to the formation of phase-space structures by particle trapping. Subsequently, these phase-space structures grow nonlinearly, in agreement with theory, but the initial state is not subcritical.

To be clear, our conclusion is not that the theory is wrong. The theory applies to the growth of a hole that is already present, not to the growth of a hole from random perturbations. To assess the nonlinear stability, we must understand what kind of initial conditions are unstable. In the following section, we address the stability of an artificial seed structure.

Since phase-space density is conserved along particle trajectories, the center of a hole, where particles are deeply trapped, and which therefore follows particle orbits, must conserve f. Therefore, an isolated hole can grow (decay) by climbing (descending) a velocity gradient. When the gradient is positive (negative), it must accelerate (decelerate) to grow and decelerate (accelerate) to decay. In figure 11, we observe that the artificial electron seed initially grows, from ω_p,et = 0 to 400, by climbing the positive velocity gradient, thereby accelerating from v_h/v_T,i = 0.8 to 3.0. The growth stops when the structure reaches the top of the electron distribution (v_h = v_d). Then, from ω_p,et = 0 to 400, it decays by descending the velocity gradient, while still accelerating. It decays until it reaches a velocity v such that $f_{0,{\rm e}}(v)\approx f_{0,{\rm e}}(v_{\rm h}(0))-\tilde{f}_{\rm h}(0)$. Before this final velocity is reached, however, the diagnostic switches to other holes, which are then deeper than the initial seed. It is these new holes that drive the nonlinear instability at later times. This process is described in the next section, 4.2.

Whether an artificial hole-like seed initially grows or not depends on both its characteristics and the plasma drift velocity. Figure 12 shows the evolution of an electron hole-like seed, initially located in the region of strong overlap between ion and electron distributions, for various initial conditions. The ratio Δx_h/Δv_h = 20/ω_p,e is arbitrary. For v_d = 0 (figure 12(a)), we observe that all seeds (for the shapes and sizes we tested, as listed in the legend of figure 12) are damped. This is expected since in this configuration, there is no free-energy. Note that in this configuration, an isolated hole must accelerate (descend the velocity gradient) to decay. Trapped particles accelerate with the hole. Thus, a phase-space structure can drive transient velocity-space particle transport, even as it decays.

Zoom In Zoom Out Reset image size

For v_d/v_T,i = 3.0 (figure 12(c)), which is relatively close to linear threshold (v_d = 0.76 v_d,cr), all seeds (for the shapes and sizes we tested) initially grow. The evolution is similar to the reference case S. For v_d/v_T,i = 2.5, results are qualitatively similar and are not shown here.

For an intermediate value of drift, v_d/v_T,i = 1.5 (figure 12(b)), which is far from the linear threshold (v_d = 0.38v_d,cr), the hole growth and subsequent decay depends on its size and location. When they are located in the velocity-region of strong overlap between ion and electrons, holes seem to grow more easily. This is consistent with the underlying growth mechanism (local momentum exchange between ions, electrons and waves) and with the predicted theoretical growth rate of an isolated hole. The hole growth is expected to be of the order of $v_{T,{\rm e}}^2 v_{T,{\rm i}}^2 \partial _v f_{\rme,0}(v_{\rm h}) \partial _v f_{\rmi,0}(v_{\rm h}) \omega_{\rm b}$, where ω_b is the bounce frequency of a particle deeply trapped in the hole [32]. However, in all cases, the hole eventually decays and no other phase-space structure is generated. Whether a hole with a different shape can or cannot trigger phase-space turbulence remains an open question, since we have studied only six cases.

In (figure 12(b)), we observe that holes initially decay for a while before starting to grow. The reason may be that the shape of the artificial initial seed is not optimal for growth. This hypothesis is supported by the following analysis. Figure 13 compares two cases, where Δv_h(0) is fixed to 0.2 v_T,i, but Δx_h(0) differs. If Δx_h(0)/λ_D = 2, which corresponds to a ratio $\Delta x_{\rm h} (0) / \Delta v_{\rm h}(0) = 20 \omega _{\rm p,e}^{-1}$, the ratio Δx_h(t)/Δv_h(t) executes damped oscillations around a value of $200\omega _{\rm p,e}^{-1}$, until it stabilizes. Then the hole starts its nonlinear growth. If Δx_h(0)/λ_D = 5, which corresponds to a ratio $\Delta x_{\rm h} (0) / \Delta v_{\rm h}(0) \approx 50 \omega _{\rm p,e}^{-1}$, the hole starts to grow almost straightaway, suggesting that its shape is closer to the optimal shape for nonlinear growth.

Zoom In Zoom Out Reset image size

For a velocity drift v_d/v_T,i = 1.5, we mentioned that, although the hole initially grows, it eventually decays and no other phase-space structure is generated. In such cases, we observe that the electric potential does not increase. The velocity redistribution is negligible (〈δf_s〉/f_s,0 < 10%), the mean thermal energies are constants (δT_s/T_s < 0.1%), and the anomalous resistivity is small (η/η_coll < 500). Thus, the nonlinear growth of one isolated hole is observed without system-wide instability. We refer to this situation as local hole growth.

This local hole growth contrasts with the global subcritical instability observed e.g. for the reference case S. Figure 14 shows snapshots of the velocity distributions in the latter case. We observe wide particle redistribution between $t=1000 \omega_{\rm p,e}^{-1}$ and $t=2000 \omega_{\rm p,e}^{-1}$, which is after the initial hole has decayed, and during the growth of the new holes. We can check that particles are indeed trapped into the self-emerging structures. Figure 15 shows the velocity v_p of a test electron, which is deeply trapped by a hole formed around $t=800 \omega_{\rm p,e}^{-1}$. The central velocity of the hole v_h is estimated by tracking the local minimum of f. The velocity of the test electron in the framework of the hole, v_p − v_h, is shown to oscillate around zero, which indicate that particles follow trapped orbits in phase-space, in the reference-frame of the structure. The bounce frequency is measured from the time-trace of v_p − v_h as ω_b ≈ 0.15ω_p,e. Another way to estimate the bounce frequency is to measure the local maximum φ₀ of the electric potential, and the spatial extent Δx_h of the negative perturbation. The bounce frequency is then given by $\omega_{\rm b}^2=|q_{\rm e}|k_{\rm h} ^2\varphi_0 /m_{\rm e}$, with k_h = 2π/Δx_h. This method also yields ω_b ≈ 0.15ω_p,e for the same hole. The oscillation is quasi-periodic from $t_1=850 \omega_{\rm p,e}^{-1}$ to $t_2=1280 \omega_{\rm p,e}^{-1}$. At t₂, the hole appears to be sheared by the tidal forces exerted by a neighboring, larger hole. The lifetime of the hole is thus estimated as τ_c = t₂ − t₁ = 430. The Kubo number in the region of phase-space within this hole is K = ω_bτ_c/(2π) ≈ 10. Repeating this analysis for other holes among the largest ones, we found Kubo numbers in the range K ≈ 3–20. This confirms that the simulation is in a regime of large Kubo number. The merging of holes reduces their lifetime. However, merging is rare enough that particles bounce many times during the life of most large holes.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 16 shows the moving average (over $\delta t=16 \omega_{\rm p,e}^{-1})$ of the potential energy time-series in three cases, including the reference case S, for v_h/v_T,i = 0.8. We observe that the potential energy grows to eφ/T ∼ 0.3, which is of the same order as the saturated potential in linearly unstable cases. In the reference case S, we see a clear difference between the first phase (t = 0 − 700) and the second phase (t = 700 − 2000). In the first phase, a single hole develops (as seen in figures 9(c), (e) and (g)), and although the field energy grows, it is only transiently. The field energy then decays back to a value close to the initial perturbation. In the second phase, where multiple holes develop (as seen in figure 9(i) and (b)), the field energy grows to eφ/T ∼ 0.3 before relaxing to eφ/T ∼ 0.1, thereby driving a global subcritical instability. Therefore, for a given initial field energy, phase-space turbulence (multiple interacting holes) is more efficient than a single hole to drive the instability.

Zoom In Zoom Out Reset image size

Figure 4 includes the time-evolution of anomalous resistivity and perturbed mean thermal energies for the reference case S. We observe large anomalous resistivity and turbulent heating, qualitatively similar to the cases with an initial ensemble of waves. These subcritical instabilities occur relatively far from marginal stability, even when the initial potential energy is as low as eφ/T ∼ 10⁻². This is in sharp contrast with the case where the initial perturbation is a collection of waves. In other words, phase-space structures, even with non-optimal shapes, are much more efficient than coherent waves for driving nonlinear instabilities.

To summarize, we observe that many holes, even small, but when scattered in phase-space, can drive global subcritical instabilities. In contrast, one single hole, even a large one, can evolve while leaving most of the phase-space intact, without system-wide instability (without significant potential energy growth, redistribution, heating or anomalous resistivity). One single hole can drive instabilities indirectly though, by triggering the formation of many smaller holes in its wake. This process is detailed as follows. As can be seen in figure 9, as the initial, artificial hole accelerates within the region v_h < v_d, its depth increases (along with its width in velocity). This increase in depth is due to the trapping of additional particles. This leaves a trail of negative phase-space density perturbations in the region v_h < v_d (see figures 9(c) or (e)). Then, as the hole enters the region v_h > v_d, its width in velocity decreases, and de-trapping occurs. This, in turn, leaves a trail of positive phase-space density perturbation in the region v_h > v_d (see figure 9(g)). In analogy to self-gravitating matter organizing into hierarchical structures via the mechanism of Jeans collapse, negative phase-space density perturbations have a natural propensity to coalesce [8, 43]. The negative trail bunches into a collection of small holes, scattered in phase-space (see figures 9(g) and (i)). The turbulent interaction of these many holes (phase-space turbulence), is shown in figure 10. From our analysis, we conclude that phase-space turbulence is a very efficient source of particle-transport in velocity space (or redistribution), turbulent heating and anomalous resistivity.

The above nonlinear stability analysis of phase-space holes is summarized in table 1, which shows whether a hole decays (H↓), grows but then decays without triggering system-wide electron redistribution (local growth, H∼) or grows and trigger such redistribution (global subcritical instability, H↑). From this table, we conclude that the nonlinear stability threshold with an initial hole in terms of v_d/v_T,i is between 1.5 and 2.5.

To study the effect of larger wave-lengths, we ran many of the simulations above with a quadrupled system size, allowing wave-numbers as small as $0.05 \lambda_D^{-1}$. We did not find any qualitative difference in the results.

Recently, we have discovered a new kind of self-organized structure, called a phase-space jet [35]. In figure 9(i), we superimposed a dashed curve of constant f_e (not constant $\tilde{f}_\rme)$, which spans a velocity range 3v_T,i = 1.5v_T,e. This structure is highly anisotropic. It has a lifetime $\tau _{\rm jet} \approx 20 \omega_{\rm p,e}^{-1}$, which is comparable to the average time it takes a particle to change its velocity by v_T,e, $\tau _{\rm travel} \approx 25\omega_{\rm p,e}^{-1}$. Thus it is a phase-space jet, which can facilitate electron redistribution. We conclude that phase-space jets are spontaneously created in subcritical (as well as linearly unstable) conditions.

Our numerical analysis clarify the process of phase-space structures formation. If the system is linearly unstable, a turbulent state can be reached, in which particles are randomly scattered, leading to fine-grain structures that act as seeds. If the system is linearly stable, we can speculate that at least four routes to instability are possible. The first route was demonstrated in figure 3. It corresponds to the growth from random fluctuations. This route is limited to an initial barely stable equilibrium and requires large amplitude perturbations. The initial perturbation (e.g. thermal noise) can be seen as an ensemble of waves, which will trap particles and form seed structures that can grow nonlinearly. The second route was demonstrated in figure 12(c). It corresponds to the growth of a single hole. Such holes may be externally driven by the experimental setup or by physical processes that are not included in this model. As we have seen, a single hole may or may not lead to phase-space turbulence. We speculate the existence of a third route, which would be a transition from supercritical to subcritical instability on a fluid time-scale. Fluid parameters of the background plasma such as fluid velocities or temperatures may evolve, on a slow time-scale, due to processes that are not included in this model, such as an applied electric field, external heating, or other magneto-hydrodynamic instabilities. This may lead to a transition from a linearly unstable system to a linearly stable state. Phase-space turbulence that originates from linearly-driven seeds should be able to survive this transition. A fourth route was demonstrated by NGuyen [44, 45]. It is a transition from supercritical to metastable (subcritical steady-state) on a kinetic time-scale. As the wave grows linearly, the ion resonance width increases. Eventually, trapped ions absorb wave energy at a rate for which the total nonlinear growth rate vanishes. Such process can result in a metastable steady-state, where phase-space structures may be continuously created and dissipated. These scenarios are summarized in table 2.

This work is concerned with collisionless plasma, but even small collisions can have qualitative effects on the nonlinear evolution of wave-particle interactions [12, 46]. If a collision operator that tries to recover a Gaussian distribution is introduced, we expect to find regimes of intermittent, rather than transient, turbulence. This is a speculation, based on the known effects of collisions on phase-space structures in the bump-on-tail instability.

Phase-space holes in pair plasmas with small mass ratios are ubiquitous in semiconductors and space [47]. However, the small mass ratio m_i/m_e = 4 adopted in this work brings the question of applicability of our findings to the most common hydrogen-ion plasma. In the opposite limit of an electron–oxygen-ion plasma with mass ratio m_i/m_e = 29 500, a single electron hole remains stationary for a hundred $\omega_{\rm p,e}^{-1}$, until an ion density cavity is formed [48]. Whether phase-space turbulence and subcritical instabilities develop or not on a ion kinetic time-scale $(\omega_{\rm p,i}^{-1})$ in large mass-ratio plasma remains an open question. Moreover for larger mass ratio, the phase-space turbulence will probably not affect the ion distribution and its mean thermal energy.

While in this work we focused on current-driven ion-acoustic instability in unmagnetized plasmas, similar scenarios can be developed for magnetized plasmas. Formation of a single phase-space structure is reported in [25]. In that paper, it was shown that a drift-hole extracts free energy more efficiently than linear waves do. The turbulent case with many structures is discussed in [15]. In this work it was shown that in trapped-ion resonance driven, ion temperature gradient instability, transport is determined, not by quasi-linear turbulent diffusion, but rather by dynamical frictions exerted on turbulent trapped-ion granulations. More recently, it was shown that both a coherent drift-hole and an ensemble of granulations can interact with zonal flows [20, 49, 50]. The impact of zonal flows on transport driven by trapped-ion granulations was formulated as a part of dynamical friction [51]. In magnetized space plasmas, phase-space turbulence and jets are promising candidates mediators of magnetic reconnection via anomalous resistivity.