Nonlinear kinetic simulation study of the ion–ion streaming instability in single- and multi-ion species plasmas

The electrostatic instability driven in plasma by ions that counterstream at approximately sonic speeds is of interest in several active research areas due to its ability to excite a large electrostatic potential and heat ions. Laboratory experiments^1–3 using opposed laser-irradiated foils have demonstrated the efficacy of the ion–ion streaming instability^4,5 (IISI) to increase ion temperatures at rates that cannot be explained by collisional (fluid-like) processes in the counterstreaming plasma. An IISI has also been shown to enhance the collisional friction between plasma flows, modifying flow speeds.⁶ Collisonless shocks, where ions are reflected and counterstream at the shock interface, are thought to be present in many astrophysical contexts with scales ranging from the Earth's bow shock to galactic clusters, and have long been investigated as an explanation for observed cosmic ray spectra.^7–11 Similarly, collisonless shocks have recently been utilized in proof-of-principle experiments^12–14 to generate high-quality $≥10$ MeV proton beams. The IISI can lead to the breakup of such shocks^15,16 and has been suggested as a mechanism of enhanced ion acceleration.¹⁷ 

The IISI can be saturated by both ion heating^4,5,18–20 and ion trapping.^5,21 The two mechanisms are not independent: ion trapping can convert the electrostatic energy of the ion waves into particle kinetic energy (which, after phase space filamentation and subsequent diffusion in velocity, can be viewed as a change in temperature that is generally anisotropic), while a change in ion temperature will alter the strength of all wave-particle interaction processes in addition to directly stabilizing the IISI.

In the following, we demonstrate various nonlinear behaviors of the IISI, covering saturation, sensitivity to initial flow velocity, and the unique consequence of ion trapping when ion species of different charge-to-mass ratios are present. Our results are obtained using a fully kinetic continuum method^22–24 that is sixth-order accurate in all phase space directions and time to describe the evolution of distribution functions. This approach is novel in the study of this instability and enables accurate diagnosis of growth rates and instability thresholds. In our simulations, we measure growth rates that agree with linear theory and show phase space structures that can be understood using simple analytic estimates of the ion trapping regions in phase space. These structures persist long after the ions are detrapped. We also show that ion trapping in a multi-ion species plasma can preferentially heat the heavier ion species and is able to dramatically alter the streaming velocity of the lighter species, offering a new species velocity separation mechanism.

We study interpenetrating streams of hydrogen (H) and carbon (C). Pure hydrogen plasmas (H–H interpenetration) are commonly occurring in a variety of applications, such as the solar wind or other astrophysical contexts, making a natural case of interest. To exhibit the novel phenomena that arise when ions of differing mass (and, importantly, charge-to-mass ratio) are present, we study interpenetration where one stream is pure H and the other is pure C. To demonstrate velocity separation via the IISI, we also present results from two interpenetrating streams of CH₂. Such C/H mixtures are common in experiments featuring laser-driven plasma ablation and expansion, such as in Refs. 1 and 2.

The layout of this article is as follows: first, the plasma parameters of focus here are introduced. Linear quantities, such as growth rates and the wave vector at which the maximum growth rate occurs, are calculated. The simulation method is then discussed. Next, results from hydrogen streams of differing initial flow velocity are presented, followed by hydrogen–carbon mixtures. Finally, we summarize and discuss the implications of our results.

We study the IISI using (i) two counterstreaming flows of hydrogen with a range of initial relative flow speeds u that change the qualitative behavior of the IISI, and (ii) counterstreaming flows of hydrogen and carbon. To remove ambiguity, electromagnetic (Weibel) instabilities are precluded by choosing to solve a Vlasov–Poisson system. The Weibel instability is typically slow compared to the IISI,²⁵ meaning that the IISI will generally saturate before the subsequent onset of the Weibel instability.¹⁵ 

Due to the timescales considered as well as for simplicity, the plasma is treated as collisionless. The system is discretized in two Cartesian space and velocity dimensions (2D2V) and solved using the LOKI^22–24 code. Continuum methods are not subject to numerical noise above machine precision and our simulations are therefore free of spurious heating, although we introduce thermal fluctuations to provide a broad-spectrum initial density perturbation (these density perturbations are $O(10−8)$ in Fourier space). The sixth-order scheme permits larger grid point intervals in space and time for a given numerical accuracy than lower-order methods. The latter property is particularly desirable due to the computationally demanding fully kinetic description of both the electrons and the ions (with a physical mass ratio) in our simulations. Simpler quasilinear descriptions lack the necessary physics of ion trapping.⁵ In addition to the 2D2V Vlasov simulations presented here, we were also able to reproduce the nonlinear states using the particle-in-cell code OSIRIS.²⁶ 

The plasma is initially homogeneous apart from the initial low density fluctuation level. The spatial boundaries are chosen to be periodic, facilitating a clear view of the instability (the methodology is similar to Ref. 27). The velocity grids are truncated at the maximum and minimum velocities for each species (see Table II) that are of sufficient magnitude such that the details of the boundaries are unimportant to the solution; further information is provided in Ref. 24. In all cases, the distribution functions of the electrons and ions must be evolved on separate grids since they have different masses and charges. For H/H interpenetration (but not C/H or CH₂/CH₂), one may choose to evolve the two ions streams as a single combined distribution function or as two separate distribution functions; we confirmed that this choice has no bearing on the solution, although for diagnostic purposes, evolving the two ion streams separately is convenient and instructive.

The plasma properties are summarized in Table I. In all cases, the plasma is fully ionized. The plasma is initialized in a frame where there are two (or four, as in case 5) shifted isotropic Maxwellian ion distributions of number density n_i and temperature T_i, counterstreaming with velocities $ui=(±u,0)$⁠. For all cases, the total n_i is equal for the opposing streams. The electrons are initialized with a single charge-neutralizing and current-canceling Maxwellian distribution of temperature T_e, and are hot relative to the ions. Table I lists the wave numbers $k=(k||,k⊥)$ and growth rates γ_k of the fastest-growing (superscript “max”) Fourier modes excited by the IISI in each case. The components of vectors are defined as parallel (‖) or perpendicular (⁠$⊥$⁠) to the flow direction. Throughout, velocities are explicitly normalized to the cold long-wavelength hydrogen plasma sound speed $cs=Te/mp$⁠, frequencies to the hydrogen plasma frequency $ωpi=4πnee2/mp$⁠, where n_e is the average electron density in the system, and distances to the electron Debye length λ_De. Note that for consistency, these definitions of the reference quantities c_s and ω_pi remain fixed throughout, even when considering C/H interactions. Note also that the actual sound speed of a given mode is generally lower than c_s and will depend on its wavelength; this correction is not necessarily small since one has $kmaxλDe∼1$⁠. “Hat” notation denotes a unit-vector, e.g., $k̂=k/k$⁠.

The kinetic dispersion relation of the system,

may be solved to find the unstable modes. Here, $Ωkj=ωk−k·uj+ ιγk$⁠; $ωk,k,γk∈ℝ$⁠; $ι=−1$⁠; and the species index j includes both electrons (e) and ions (i). $fj=fj(t,r,v)$ and ω_pj are, respectively, the species distribution function and plasma frequency. t, r, and v are the time, space, and velocity variables. In all cases shown here, the susceptibilities $χj=−[1/(2k2λDj2)]Z′[Ωkj/(2kvtj)]$ are evaluated using the plasma dispersion function²⁸ Z′, where $λDj=vtj/ωpj$ is the species Debye length for which $vtj=Tj/mj$⁠, and T_j and m_j are the species temperature and mass, respectively. We take m_H = m_p, $mC=12mp$⁠, and $mp=1836me$⁠. Standard fluid approximations to χ_i for the parameters used in cases 1–5 are only in qualitative agreement with exact numerical solutions, while for the electrons, one has $vte≫|Re[Ωke/k]|$ and the susceptibility is well-approximated by their adiabatic response, $χe≃1/(k2λDe2)$⁠. In the chosen (flow symmetric) frame, one has for cases 1–3 and 5 the real frequency $ωk=0$ for unstable modes (growth rate $γk>0$⁠), as noted in Ref. 4, where the subscript k denotes a given Fourier mode with wave vector k.

The numerical resolutions and system sizes are summarized in Table II. The spatial dimensions of the system in each case are tailored to describe the unstable modes in the system.

The linear theory of the IISI has been studied extensively previously.^4,5,18–20 For identical ion streams with equal and opposite flows, the IISI is unstable when $u≳1.3vti$ and $Ti/(ZTe)$ is less than a critical value. For $u≳2.1vti$⁠, instability requires $Ti/(ZTe) ≲ 0.28$⁠, and this criterion becomes more stringent²⁰ for lower u. For values of u close to the stability boundary, the fastest growing mode has $k⊥max=0$ and the IISI is approximately 1D in nature. However, when the IISI is unstable, there are always unstable modes with $k⊥>0$ and therefore, after the system has become nonlinear, there will also be some amount of transverse ion heating (this is discussed later). For $|k̂max·ui|≳cs$⁠, one has $k⊥max>0$ and the instability is necessarily 2D even in the linear stage of evolution. For sufficiently large u, the IISI is aligned nearly perpendicular to the flow and there are no unstable modes with $k⊥=0$⁠. In this section, we explore the nonlinear states of ion streams as u is varied between these limits.

The evolution of the simulated system is qualitatively similar in all cases shown here: From a quiescent plasma, the IISI excites a spectrum of modes with an associated electrostatic potential. The fastest-growing modes are saturated by ions trapped in this potential (a form of “wave breaking”²⁹). Since the IISI excites a spectrum of modes, a velocity phase-space diffusion then takes place with a partitioning of heating in the directions parallel and transverse to the flow that depends on $kmax$ (and therefore on the initial flow velocity) and the ion streams undergo some amount of slowing. The electrostatic field energy associated with the excited potential is transferred back to the particles and the IISI is fully stabilized by a skewed heating of the ion distributions. The time evolution of the electrostatic field energy, kinetic energy, flow velocity, and ion temperature is shown in Fig. 1; a detailed discussion of these quantities follows.

In Fig. 2, numerical solutions to the dispersion relation of Eq. (1) are compared with values of γ_k measured in simulations of cases 1–3 for the linear stage of growth using fits within the interval $ωpit=[0,∼80]$⁠. Measuring γ_k for the exponentially growing unstable IISI modes in our simulations is straightforward due to the numerical method employed, where Fourier modes grow from thermal fluctuations through $∼8$ orders of magnitude before reaching saturation. An electron–ion streaming instability (EISI) is also present in the system.⁴ However, the EISI is here one or more orders of magnitude slower than the IISI and does not play a role in our conclusions. Our results show no change if instead two electron populations are initialized, where each cancels the current of its associated ion stream.

The nonlinear evolution of the IISI is strongly affected by u. In Fig. 3, the hydrogen ion distributions $⟨fi⟩$ are shown for the three qualitatively different cases, where for the most unstable mode the ion streams are marginally subsonic (top row); marginally supersonic (middle row); and strongly supersonic (bottom row), corresponding to cases 1–3 in Table I, respectively. $⟨…⟩$ denotes averaging over all space.

The differing phase space structures can be understood by considering the electrostatic potential $Φ$ excited by the IISI, which in turn is dictated by the most unstable wave vectors (see Table I), and the associated ion trapping. In Fig. 4(a), $Φ$ at the maximum in time of the electrostatic energy E_es [the time history of E_es is shown in Fig. 1(a)] is shown for case 1. In Fig. 4(b), f_i sampled at $r⊥=0$ and integrated over $v⊥$ is shown, exhibiting characteristic ion trapping structures. Later in time, $Φ$ decays to a turbulent spectrum [see Fig. 4(c)], the ions detrap [see Fig. 4(d)], and E_es is converted via velocity space diffusion to particle kinetic energy, shown in Fig. 1(a). In Figs. 4(e) and 4(f), $Φ$ at the maximum in time of E_es is shown for cases 2 and 3. The structures in $Φ$ shown here correspond to large density modulations of $∼10%$ in case 1 and $∼30%$ in cases 2 and 3. The eventual turbulent states of cases 2 and 3 are similar to that shown for case 1.

Electrons and ions with velocities near the phase velocity of the wave are resonant with the unstable (⁠$γk>0$⁠) modes. In symmetric systems such as cases 1–3 and 5, the phase velocity $vk=(ωk/k)k̂$ is zero because the real frequency of the unstable modes is zero in the simulation frame (see Table I). Such particles can become trapped by the electrostatic potential and follow trapped orbits, generally with excursions in both $r||$ and $r⊥$⁠. γ_k is symmetric under changes of signs of $k||$ and $k⊥$ and as a result the resonant region of velocity for species j is enclosed by a diamond (except for the special case $k⊥=0$⁠) given by $k̂·(v−vk)=±vtr,j$ where the trapping width $vtr,j$ is given by $vtr,j/vtj=2e|Zj|ϕ/Tj$⁠. For a system dominated by the fastest-growing linear mode, $k=(±k||max,±k⊥max)$⁠. Z_j is the charge state of species j and $ϕ$ is a potential that characterizes the half-depth of the potential well experienced by the particles.

In Fig. 3, regions in velocity space resonant with $⟨fi⟩$ are shown for the modes $k=kmax$ (see Table I). Regions bound by both $vtr,imax=vtr,i(0.5[max(Φ)−min(Φ)])$⁠, and $vtr,iRMS=vtr,i[RMS(Φ)]$ are indicated by dotted and dashed lines, respectively. RMS denotes a root-mean-squared value, calculated by averaging over the entire system. $vtr,imax$ defines the furthest extent of resonant interaction, while $vtr,iRMS$ is an approximate average over the 2D structures in $Φ$⁠. The arrows indicate the direction of the trapped ion oscillations. Not shown, the electron distribution function f_e is also modified, exhibiting a trapping-induced flattening of a few percent centered about v_k. This modification of f_e does not appear to impact the evolution of the IISI in the cases studied here.

By examining $vtr,i$⁠, it is apparent that the amplitude of a mode driven by the IISI is limited to the intersection of the resonant region with the bulk of the ion distribution in cases 1–3 (known as wave breaking²⁹), i.e., $Φ$ will grow until $vtr,iRMS≃k̂max·ui$⁠, which is accurate to within 10% for cases 1–3. After saturation of the $k∼kmax$ modes, the spectrum becomes sufficiently broad such that $Φ$ causes a diffusion in velocity space, fully stabilizing the IISI.

In cases 1, 2, and 3, the average flow velocity defined as $ui(t)=⟨∫vv||fi dv⟩/ni$ approaches a reduction in magnitude by 30%, 38%, and just 4%, respectively, shown in Fig. 1(b). However, the maximum of $⟨fi⟩$ in velocity space moves little if at all, as is apparent in Fig. 3. In fact, for case 2, the position of the maximum increases in magnitude by 14% relative to its initial value due to particles being pulled away by trapping predominantly from the side of f closest to the resonance. We define the time-varying spatially averaged directional temperature as $Tj=(mj/nj)⟨∫v(v−uj)2fj dv⟩$⁠. This is the temperature that a Maxwellian distribution of species j would need in order to have a kinetic energy equal to that of f_j in the frame moving at u_j, calculated independently for the parallel and perpendicular directions. T_j is shown for cases 1–3 in Fig. 1(c). Despite f_e exhibiting flattening near v = 0 due to trapping, T_e deviates from its initial value by only $∼1%$⁠.

For supersonic flows, the scalings of the final flow velocities and temperatures with the initial flow velocities follow from the behavior of $kmax$⁠: at higher flow velocities, the most unstable mode moves to larger angles relative to the flow (so that $kmax·u$ approximately satisfies the acoustic dispersion relation), and the action of trapping on the distribution is likewise directed at larger angles relative to the flow. In case 3, the ion bounce motion is almost perpendicular to the flow, resulting in almost exclusively transverse ion heating and little slowing down. However, the final value of $Ti||$ is not simply a monotonic function of $u/cs$⁠. As u decreases to subsonic values and the IISI becomes effectively one-dimensional (as in case 1), one has at saturation $vtr,iRMS=u$ and the final value of $Ti||$ falls accordingly until, for a sufficiently low u, the instability is simply below threshold and does not occur.

A quasilinear description, such as that explored by Forslund,⁵ is in qualitative agreement with the results shown for cases 1–3, correctly recovering that perpendicular ion heating dominates parallel heating at high flow velocities and that the slowing down due to the IISI for symmetric single-ion species flows is modest. However, any feature of a kinetic simulation that is a direct consequence of trapping will be missed by a quasilinear treatment. For example, a quasilinear model will not include the excursion of trapped ions to a distance $vtr,i$ beyond the point in velocity space at which they are resonant, leading to a significant underestimate of the final perpendicular ion temperature for supersonic flows: in case 3, one has $Ti⊥/Te≈0.5$⁠, whereas the prediction from a quasilinear model⁵ is $Ti⊥/Te≈0.25$⁠.

We now consider case 4 (see Table I), where one stream is entirely comprised of hydrogen (H) ions and the other is fully ionized carbon (C). The initial ion temperatures and number densities are equal. Shown in Figs. 5(a) and 5(b), the theoretical linear kinetic growth rate and the measured numeric growth rate are again in quantitative agreement. In Fig. 5(c), $Φ$ taken at the maximum in time of E_es is shown. The structure of $Φ$ is qualitatively similar to case 2 due to the similar range of unstable modes. However, in cases 1–3, the electron and ion density modulations are in phase. In contrast, the electrons and C ions are in phase in case 4, while the H ions (notably with $O(1)$ density modulations) are approximately anti-phased. In fact, the phase relationship between each component [shown at the peak of E_es in Figs. 5(d)–5(f)] evolves as T_j increases in time, similar to the “slow mode” described in Refs. 30 and 31.

The trapping widths of H and C ions are given by $vtr,i/cs=2Zi/Aieϕ/Te$ and therefore $vtr,H/vtr,C=2$⁠, where A_i is the ion mass number. As a result, the conventional wave-breaking limit²⁹ of H ions (when the trapping width intersects the bulk) will be reached before that of the C ions. The evolution of $⟨fi⟩$ for case 4 is shown in Fig. 6, and the associated energy partitioning, flow velocity, and temperature evolution in time are shown in Fig. 7.

As before, the most unstable linear IISI mode and ion trapping together dictate the structure in velocity space. However, in this case, the H ions alone are not able to saturate the growth of the IISI by trapping: $vtr,j$ grows until the resonant region fully encompasses the entire H ion distribution [Fig. 6(c)], and is instead saturated by trapping of the C ions (at the C ion wave breaking limit) where $vtr,CRMS≃k̂max·uC$⁠. Note that for the phase velocity of the most unstable mode, one has $vk=ωkmax/kmax≪k̂max·uj$ (see Table I) and therefore the center of the resonance is far from the bulk of either ion distribution. This greater likelihood of H ions to become trapped due to their larger trapping width relative to C ions is then reinforced by the time required for ions to complete a trapped orbit: for the bounce frequency of a trapped particle, one has $ωb,j=ke|Zj|ϕ/mj$ and the time taken to complete a trapped orbit is $τb,j=2π/ωb,j$⁠. As $ϕ$ falls during the latter stages of the evolution of the system, $τb,j$ becomes large and ions may not fill out phase space as fully as if they were to experience an adiabatically decreasing amplitude. This effect is more pronounced for the C ions since $τb,C/τb,H=2$⁠. The longer bounce period of the C ions is readily apparent in the kinetic energy oscillations shown in Fig. 7(a) (compare K_H and K_C).

As a result of the total trapping of the H ions, the average flow speed of the H ions becomes that of the dominant IISI mode, which initially has $vk≪k̂max·uj$⁠, i.e., the H ions are effectively stopped by the C ions due to the IISI. Conserving momentum, the C ions are also slowed, but to a lesser degree (⁠$∼8%$⁠) due to their larger mass [see Fig. 7(b)]. The electrons cancel out the current of the ions and so are slowed relative to their initial speed by an amount similar to the slowing of the ions, although this slowing is small (⁠$∼1$%) compared to v_te. Similar to the cases discussed in Sec. III A, the electron distribution is flattened by a few percent but does not deviate from its initial temperature by more than $∼1$%. In Fig. 7(c), the evolution of T_j is shown. The heating of the C ions in the direction parallel to the flow is striking (⁠$TC||/Te→∼0.5$⁠) and exceeds the perpendicular C ion heating or the H ion heating in either direction by a factor of 2, while as before T_e is effectively unchanged. Analysis of $⟨fC⟩$ reveals a ring-like formation [Fig. 6(f)] that is qualitatively similar to case 2 [Fig. 3(f)].

Finally, we discuss case 5, where we simulate counterstreaming flows of CH₂. Here, there are two symmetric counterstreaming flows of C ions and two symmetric counterstreaming flows of H ions, where there are twice as many H ions as C ions and initially one has $uH,C=(±u,0)$⁠. As before, growth rates for the IISI from theory and our simulations (not shown) agree in the linear stage, and this case is qualitatively similar to case 4. However, due to the symmetry of case 5, the flow of H ions is now separated from the flow of C ions in each stream: the H ions are trapped more readily than the C ions, and are pulled toward the center of the frame by the action of the IISI. This phenomenon of velocity separation is readily apparent in Fig. 8(a), showing the C ions slowing by $∼26%$ and H ions slowing by $∼50%$⁠. The C ion heating shown in Fig. 8(b) is even stronger than that in case 4, most notably in the transverse direction, resulting in $TC⊥/Te→ ∼1.5$ and $TC||/Te→ ∼0.8$ (⁠$TH||,⊥/Te→ ∼0.25$⁠, unchanged from case 4).

We have studied the action of the ion–ion streaming instability (IISI) for single- and multi-ion species flows. The linear stage of evolution, and the subsequent initial saturation by trapping, is uncomplicated: in high-fidelity 2D2V Vlasov simulations, the most unstable modes are those predicted by linear theory, and the amplitude of the unstable modes is limited by ion trapping (the intersection of the trapping width with the bulk of the ion distribution, or “wave breaking”). The degree of ion heating parallel or perpendicular to the flow is dictated by the wave vector of the most unstable mode, with perpendicular heating dominating at higher flow velocities.

In plasma flows containing multiple ion species of differing charge-to-mass ratios, the trapping widths of each species will differ. Trapping of the lighter ion species alone may be insufficient to saturate the IISI; saturation instead occurs when the trapping width of the heavier ion species intersects the bulk of its distribution. The heating and slowing of each species in the flow may be qualitatively different. This provides a mechanism for species (including isotope) separation in velocity and, in a nonperiodic system, space.

In collisionless shocks, the anisotropic heating of ions by the IISI has been observed in single-ion species simulations¹⁵ to give rise to a complicated picture, with simultaneous distinct Weibel instabilities arising from ions reflected at the shock front and from the temperature anisotropy driven by the IISI. Interpreting structures excited by Weibel instabilities may require careful consideration of the role of the IISI in modifying the plasma conditions from which other slower instabilities may grow, particularly if multiple ion species are present. A significant interplay between the IISI and the electron cyclotron drift instability is also known to occur.³² In plasmas with some degree of collisionality, the IISI will also modify the collisional drag and heating by reducing relative flow speeds, as in cases 1–4, or indeed by introducing a relative flow as in case 5 (counterstreaming flows of CH₂), where the hydrogen ions are slowed more than the carbon ions within each flow. The IISI may also modify the collisional drag more directly, such as in Ref. 6.

We are pleased to acknowledge the valuable discussions with Wojciech Rozmus. This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract No. DE-AC52–07NA27344 and funded by the Laboratory Research and Development Program at LLNL under project tracking code 17-ERD-081. Computing support for this work came from the Lawrence Livermore National Laboratory (LLNL) Institutional Computing Grand Challenge program. B. J. Winjum acknowledges support from DOE under Grant Nos. DE-SC0019010 and DE-NA0003842 and UCOP under Grant No. LFR-17–449059. This document was prepared as an account of work sponsored by an agency of the United States government. Neither the United States government nor Lawrence Livermore National Security, LLC, nor any of their employees makes any warranty, expressed or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States government or Lawrence Livermore National Security, LLC. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States government or Lawrence Livermore National Security, LLC, and shall not be used for advertising or product endorsement purposes.

The data that support the findings of this study are available from the corresponding author upon reasonable request.