The nonlinear evolution of the ion–ion streaming instability (IISI) is studied using numerical techniques novel to this problem that afford direct insight into the evolution of the particle distributions of each species. During the linear phase of the instability, we demonstrate quantitative agreement with linear kinetic theory. Subsequently, the electrostatic field generated by the IISI causes ring-like velocity distributions of ions to form that are both heated and slowed to varying degrees relative to their initial flows. Due to variation in the trapping conditions for ion species of differing charge-to-mass ratio, when flows of multiple species interact, the nonlinear evolution of each species can be starkly different: we show a case where a lighter ion species is completely stopped by a heavier ion species via the IISI alone (i.e., without collisions) and, for the first time, demonstrate how the IISI can introduce a relative flow between ion species that initially have the same flow velocities, thereby separating them.

## I. INTRODUCTION

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.^{6} 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 $\u226510$ 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.^{17}

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_{2}. 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.

## II. PLASMA PARAMETERS AND SIMULATION METHOD

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,^{25} meaning that the IISI will generally saturate before the subsequent onset of the Weibel instability.^{15}

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\u22128)$ 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.^{5} 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.^{26}

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_{2}/CH_{2}), 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*, counterstreaming with velocities $ui=(\xb1u,0)$. For all cases, the total

_{i}*n*is equal for the opposing streams. The electrons are initialized with a single charge-neutralizing and current-canceling Maxwellian distribution of temperature

_{i}*T*, and are hot relative to the ions. Table I lists the wave numbers $k=(k||,k\u22a5)$ and growth rates

_{e}*γ*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 ($\u22a5$) 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 $\omega pi=4\pi nee2/mp$, where

_{k}*n*is the average electron density in the system, and distances to the electron Debye length

_{e}*λ*. Note that for consistency, these definitions of the reference quantities

_{De}*c*and

_{s}*ω*remain fixed throughout, even when considering C/H interactions. Note also that the actual sound speed of a given mode is generally lower than

_{pi}*c*and will depend on its wavelength; this correction is not necessarily small since one has $kmax\lambda De\u223c1$. “Hat” notation denotes a unit-vector, e.g., $k\u0302=k/k$.

_{s}Case . | Mats. . | $u/cs$ . | $k||max\lambda De$ . | $k\u22a5max\lambda De$ . | $\gamma kmax/\omega pi$ . | $\omega kmax/\omega pi$ . |
---|---|---|---|---|---|---|

1 | H/H | 0.51 | 1.08 | 0 | 0.14 | 0 |

2 | H/H | 0.88 | 0.68 | 0.68 | 0.16 | 0 |

3 | H/H | 2.57 | 0.23 | 0.93 | 0.16 | 0 |

4 | C/H | 0.64 | 0.68 | 0.82 | 0.11 | 0.031 |

5 | CH_{2}/CH_{2} | 0.64 | 0.57 | 1.24 | 0.13 | 0 |

Case . | Mats. . | $u/cs$ . | $k||max\lambda De$ . | $k\u22a5max\lambda De$ . | $\gamma kmax/\omega pi$ . | $\omega kmax/\omega pi$ . |
---|---|---|---|---|---|---|

1 | H/H | 0.51 | 1.08 | 0 | 0.14 | 0 |

2 | H/H | 0.88 | 0.68 | 0.68 | 0.16 | 0 |

3 | H/H | 2.57 | 0.23 | 0.93 | 0.16 | 0 |

4 | C/H | 0.64 | 0.68 | 0.82 | 0.11 | 0.031 |

5 | CH_{2}/CH_{2} | 0.64 | 0.57 | 1.24 | 0.13 | 0 |

The kinetic dispersion relation of the system,

may be solved to find the unstable modes. Here, $\Omega kj=\omega k\u2212k\xb7uj+\u2009\iota \gamma k$; $\omega k,k,\gamma k\u2208\mathbb{R}$; $\iota =\u22121$; 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 $\chi j=\u2212[1/(2k2\lambda Dj2)]Z\u2032[\Omega kj/(2kvtj)]$ are evaluated using the plasma dispersion function

^{28}

*Z′*, where $\lambda Dj=vtj/\omega pj$ is the species Debye length for which $vtj=Tj/mj$, and

*T*and

_{j}*m*are the species temperature and mass, respectively. We take

_{j}*m*=

_{H}*m*, $mC=12mp$, and $mp=1836me$. Standard fluid approximations to

_{p}*χ*for the parameters used in cases 1–5 are only in qualitative agreement with exact numerical solutions, while for the electrons, one has $vte\u226b|Re[\Omega ke/k]|$ and the susceptibility is well-approximated by their adiabatic response, $\chi e\u22431/(k2\lambda De2)$. In the chosen (flow symmetric) frame, one has for cases 1–3 and 5 the real frequency $\omega k=0$ for unstable modes (growth rate $\gamma k>0$), as noted in Ref. 4, where the subscript

_{i}*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.

Case . | $L\u0303||$ . | $L\u0303\u22a5$ . | $v\u0303H,||b\xb1$ . | $v\u0303H,\u22a5b\xb1$ . | $NvH,d$ . | $v\u0303C,||b\xb1$ . | $v\u0303C,\u22a5b\xb1$ . | $NvC,d$ . |
---|---|---|---|---|---|---|---|---|

1 | 29.2 | 62.8 | ±2.1 | ±1.7 | 64 | N/A | N/A | N/A |

2 | 46.2 | 46.2 | ±2.5 | ±2.1 | 72 | N/A | N/A | N/A |

3 | 137 | 33.1 | ±4.3 | ±3.9 | 128 | N/A | N/A | N/A |

4 | 31.4 | 31.4 | (2.1, −1.7) | ±1.7 | 64 | (1.0, −1.2) | ±1.0 | 96 |

5 | 55.1 | 25.3 | ±2.1 | ±1.7 | 64 | ±1.2 | ±1.0 | 96 |

Case . | $L\u0303||$ . | $L\u0303\u22a5$ . | $v\u0303H,||b\xb1$ . | $v\u0303H,\u22a5b\xb1$ . | $NvH,d$ . | $v\u0303C,||b\xb1$ . | $v\u0303C,\u22a5b\xb1$ . | $NvC,d$ . |
---|---|---|---|---|---|---|---|---|

1 | 29.2 | 62.8 | ±2.1 | ±1.7 | 64 | N/A | N/A | N/A |

2 | 46.2 | 46.2 | ±2.5 | ±2.1 | 72 | N/A | N/A | N/A |

3 | 137 | 33.1 | ±4.3 | ±3.9 | 128 | N/A | N/A | N/A |

4 | 31.4 | 31.4 | (2.1, −1.7) | ±1.7 | 64 | (1.0, −1.2) | ±1.0 | 96 |

5 | 55.1 | 25.3 | ±2.1 | ±1.7 | 64 | ±1.2 | ±1.0 | 96 |

## III. RESULTS

### A. Single-ion species flows

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\u22731.3vti$ and $Ti/(ZTe)$ is less than a critical value. For $u\u22732.1vti$, instability requires $Ti/(ZTe)\u2009\u2272\u20090.28$, and this criterion becomes more stringent^{20} for lower *u*. For values of *u* close to the stability boundary, the fastest growing mode has $k\u22a5max=0$ and the IISI is approximately 1D in nature. However, when the IISI is unstable, there are always unstable modes with $k\u22a5>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\u0302max\xb7ui|\u2273cs$, one has $k\u22a5max>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\u22a5=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”^{29}). 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 $\omega pit=[0,\u223c80]$. Measuring

*γ*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 $\u223c8$ orders of magnitude before reaching saturation. An electron–ion streaming instability (EISI) is also present in the system.

_{k}^{4}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 $\u27e8fi\u27e9$ 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. $\u27e8\u2026\u27e9$ denotes averaging over all space.

The differing phase space structures can be understood by considering the electrostatic potential $\Phi $ 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), $\Phi $ at the maximum in time of the electrostatic energy *E _{es}* [the time history of

*E*is shown in Fig. 1(a)] is shown for case 1. In Fig. 4(b),

_{es}*f*sampled at $r\u22a5=0$ and integrated over $v\u22a5$ is shown, exhibiting characteristic ion trapping structures. Later in time, $\Phi $ decays to a turbulent spectrum [see Fig. 4(c)], the ions detrap [see Fig. 4(d)], and

_{i}*E*is converted via velocity space diffusion to particle kinetic energy, shown in Fig. 1(a). In Figs. 4(e) and 4(f), $\Phi $ at the maximum in time of

_{es}*E*is shown for cases 2 and 3. The structures in $\Phi $ shown here correspond to large density modulations of $\u223c10%$ in case 1 and $\u223c30%$ in cases 2 and 3. The eventual turbulent states of cases 2 and 3 are similar to that shown for case 1.

_{es}Electrons and ions with velocities near the phase velocity of the wave are resonant with the unstable ($\gamma k>0$) modes. In symmetric systems such as cases 1–3 and 5, the phase velocity $vk=(\omega k/k)k\u0302$ 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\u22a5$. *γ _{k}* is symmetric under changes of signs of $k||$ and $k\u22a5$ and as a result the resonant region of velocity for species

*j*is enclosed by a diamond (except for the special case $k\u22a5=0$) given by $k\u0302\xb7(v\u2212vk)=\xb1vtr,j$ where the trapping width $vtr,j$ is given by $vtr,j/vtj=2e|Zj|\varphi /Tj$. For a system dominated by the fastest-growing linear mode, $k=(\xb1k||max,\xb1k\u22a5max)$.

*Z*is the charge state of species

_{j}*j*and $\varphi $ 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 $\u27e8fi\u27e9$ are shown for the modes $k=kmax$ (see Table I). Regions bound by both $vtr,imax=vtr,i(0.5[max(\Phi )\u2212min(\Phi )])$, and $vtr,iRMS=vtr,i[RMS(\Phi )]$ 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 $\Phi $. 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*does not appear to impact the evolution of the IISI in the cases studied here.

_{e}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^{29}), i.e., $\Phi $ will grow until $vtr,iRMS\u2243k\u0302max\xb7ui$, which is accurate to within 10% for cases 1–3. After saturation of the $k\u223ckmax$ modes, the spectrum becomes sufficiently broad such that $\Phi $ causes a diffusion in velocity space, fully stabilizing the IISI.

In cases 1, 2, and 3, the average flow velocity defined as $ui(t)=\u27e8\u222bvv||fi\u2009dv\u27e9/ni$ approaches a reduction in magnitude by 30%, 38%, and just 4%, respectively, shown in Fig. 1(b). However, the maximum of $\u27e8fi\u27e9$ 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)\u27e8\u222bv(v\u2212uj)2fj\u2009dv\u27e9$. 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*exhibiting flattening near

_{e}*v*=

*0 due to trapping,*

*T*deviates from its initial value by only $\u223c1%$.

_{e}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\xb7u$ 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,^{5} 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\u22a5/Te\u22480.5$, whereas the prediction from a quasilinear model^{5} is $Ti\u22a5/Te\u22480.25$.

### B. Multi-ion species flows

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), $\Phi $ taken at the maximum in time of *E _{es}* is shown. The structure of $\Phi $ 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*in Figs. 5(d)–5(f)] evolves as

_{es}*T*increases in time, similar to the “slow mode” described in Refs. 30 and 31.

_{j}The trapping widths of H and C ions are given by $vtr,i/cs=2Zi/Aie\varphi /Te$ and therefore $vtr,H/vtr,C=2$, where *A _{i}* is the ion mass number. As a result, the conventional wave-breaking limit

^{29}of H ions (when the trapping width intersects the bulk) will be reached before that of the C ions. The evolution of $\u27e8fi\u27e9$ 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\u2243k\u0302max\xb7uC$. Note that for the phase velocity of the most unstable mode, one has $vk=\omega kmax/kmax\u226ak\u0302max\xb7uj$ (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 $\omega b,j=ke|Zj|\varphi /mj$ and the time taken to complete a trapped orbit is $\tau b,j=2\pi /\omega b,j$. As $\varphi $ falls during the latter stages of the evolution of the system, $\tau 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 $\tau b,C/\tau 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\u226ak\u0302max\xb7uj$, 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 ($\u223c8%$) 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 ($\u223c1$%) 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 $\u223c1$%. 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\u2192\u223c0.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*is effectively unchanged. Analysis of $\u27e8fC\u27e9$ reveals a ring-like formation [Fig. 6(f)] that is qualitatively similar to case 2 [Fig. 3(f)].

_{e}Finally, we discuss case 5, where we simulate counterstreaming flows of CH_{2}. 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=(\xb1u,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 $\u223c26%$ and H ions slowing by $\u223c50%$. 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\u22a5/Te\u2192\u2009\u2009\u223c1.5$ and $TC||/Te\u2192\u2009\u2009\u223c0.8$ ($TH||,\u22a5/Te\u2192\u2009\u2009\u223c0.25$, unchanged from case 4).

## IV. CONCLUSIONS

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^{15} 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.^{32} 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_{2}), 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.

## ACKNOWLEDGMENTS

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.

## DATA AVAILABILITY

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