Influence of neutral pressure on instability enhanced friction and ion velocities at the sheath edge of two-ion-species plasmas

The speed at which ions exit a plasma is a critical parameter for the interpretation of plasma diagnostics,¹ the design of plasma-based manufacturing devices,² and plasma-material interactions in fusion experiments.³ It has recently been appreciated that plasmas consisting of multiple ion species are susceptible to ion-ion two-stream instabilities in the presheath and that these instabilities influence the flow speed of ions at the sheath edge.⁴ Ion-neutral collisions can disrupt this instability and further modify the transport properties of ions in the presheath. This work presents theory and simulations describing how the flow speed of each ion species at the sheath edge changes with neutral pressure.

Langmuir⁵ and Bohm⁶ predicted in single ion species plasmas ions reach the sheath edge at a speed exceeding the sound speed. In this paper, the sheath-edge is defined as the region where the ion density exceeds the electron density by 1% in the direction from the bulk to the sheath.⁷ For plasmas with two ion species, the Bohm criterion was generalized to be^8,9

where V_j is the flow speed at the sheath edge, n_j is the density of ion species j at the sheath edge, and n_e is the electron density. The individual sound speed of ions is $csi=kBTe/mi$⁠, where k_B is the Boltzmann constant, and T_e and m_i are the electron temperature and ion mass, respectively. The convention used in this paper is that the lighter ion species is labeled 1, and the heavier labeled 2. The plasma is expected to be quasi-neutral in the presheath so that n₁ + n₂ = n_e.

Equation (1) provides only one constraint and therefore cannot uniquely determine both V₁ and V₂. One possible solution, motivated by theoretical and numerical studies of Franklin,¹⁰ suggested that each species obtains its individual sound speed at the sheath edge. This Individual Sound Speed (ISS) limit is obtained if the ion-ion drag is weak in the presheath, which is expected based on the predicted Coulomb collision rate. However, the ISS model was challenged by experimental measurements made using Laser Induced Fluorescence^11–15 (LIF). Initial experiments seemed to indicate that ion-ion friction was strong and observed that the ion flow speeds were close to a common system sound speed [the System Sound Speed (SSS) model]

Later work¹⁶ proposed that an Instability-Enhanced Frictional Force (IEF) associated with ion-ion two-stream instabilities may be the physical mechanism responsible for the merging of the flow speeds observed in experiments.^11–15 Ions of different mass will accelerate at different rates in the presheath, creating a flow separation between species. Ion-ion two-stream instabilities can be excited if the flow separation, V₁ – V₂, is larger than a critical value ΔV_c. When two-stream instabilities occur in the presheath, V₁ − V₂ cannot significantly exceed ΔV_c because the two-stream instabilities generate an instability-enhanced friction force preventing further separation. Thus, the instability creates a “stiff transport” causing a restoring force that holds the flow difference to be ΔV_c. If no ion-ion streaming instabilities are present, then the ions are sufficiently decoupled that the ISS model is expected to be valid. The final condition for predicting the differential flow speed at the sheath edge is then

Equation (1) (assuming equality in the inequality¹⁸) and Eq. (3) constitute a closed system that uniquely determines V₁ and V₂. The threshold ΔV_c is calculated from the dielectric response function given the local plasma parameters (such as species temperature, density, and mass) in the presheath. This theory predicts that V₁ and V₂ have a strong dependence on the ion species concentrations, generally predicting that the speeds range from the ISS solution in the dilute (impurity) limit, merging toward the SSS condition at similar concentrations. Experiments that varied the ion species concentrations at low neutral pressure (<1 mTorr) were found to support the IEF predictions.^17,18 More recent work using particle-in-cell Monte-Carlo-collision (PIC-MCC) simulations of Helium-Xenon plasmas have also been shown to agree with the predictions when ΔV_c is accurately calculated numerically.¹⁹ This work also cleared up a discrepancy between theoretical predictions based on an approximate model for ΔV_c and previous PIC-MCC simulations.²⁰ 

In this paper, the IEF theory is extended to account for the influence of ion-neutral collisions. Neutral damping of the ion-ion two-stream instability changes the instability threshold flow difference, ΔV_c. The pressure dependence of ΔV_c is accounted for in the theory by including a Bhatnagar–Gross–Krook (BGK) type collision operator²¹ in the computation of the linear response function. It is shown that for modest neutral pressures, the neutral damping increases the requisite flow separation for instability, ΔV_c. Therefore, the flow speeds at the sheath edge of both ions, V₁ and V₂, become dependent upon the neutral pressure. At sufficiently high neutral pressure (10s of mTorr), the instability is completely suppressed and the plasma is well described by the ISS solution. Recent experiments were conducted that demonstrated increasing the neutral pressure of the system caused the ion flows to return to their individual sound speeds at high pressures.²² Our work presents simulations that agree with these experimental measurements and provides a theoretical framework by which the sheath edge speeds can be modeled.

This paper is divided as follows: Sec. II presents a theoretical model for ion-neutral collisions and the resulting modification to the dielectric response function. Section III uses the dielectric response function presented in Sec. II to predict the speed of each ion species at the sheath edge as neutral pressure is increased. Sections IV and V present new PIC-MCC simulations of Helium-Xenon plasma over a wide range of neutral pressures. Agreement is shown between the predicted flows at the sheath edge and the observations in the simulations. Further, a direct calculation is made connecting the instabilities with an enhancement of the ion-ion friction force which has not been done in previous works. Section VI presents simulations that suggest deviations from Maxwellian ion velocity distribution functions (IVDF) at the sheath edge resulting from ion holes in phase space associated with the two-stream instability. Section VII presents further theoretical predictions relevant to experimental conditions.

In this section, three different models for calculating ΔV_c from the dielectric linear response function are presented. The three models considered are no ion-neutral collisions (Sec. II A), the BGK collision operator (Sec. II B), and the Krook collision operator (Sec. II C). A linear response description is justified because the scale over which the plasma potential varies in the presheath for the experimental conditions of interest (typically on the order of centimeters) is much longer than the Debye-length-scale wavelengths of interest (typically a fraction of a millimeter).

In each case, the ion distribution functions are modeled as drifting Maxwellians. This basic assumption can break down at high neutral pressures where charge exchange collisions can alter the IVDFs. However, comparison with the PIC-MCC simulations will show that the predictions based on drifting Maxwellians still agree well with the simulated ion flow speeds based on moments of the IVDFs. Ion-ion collisions are negligible in the presheath at low temperature plasma conditions (⁠$ne∼108−1010 cm−3,Te∼1−10 eV$⁠) if it is stable and do not significantly modify the thresholds for instability. Ion-ion collisions only become dominant after the plasma surpasses the two-stream instability threshold.

For a collisionless one-dimensional plasma, the linearized distribution function evolves in accordance with the Vlasov equation

where $fi(0)$ and $fi(1)$ are the ambient and perturbed distribution functions, and q_i is the charge of the species. Solving Eq. (4) for the density perturbation and using Gauss's law to relate this to the potential perturbation provides the linear dielectric response function.²³ For a plasma with two ion species described by drifting Maxwellian velocity distribution functions, and adiabatic electrons, this can be written as

where k is the wavenumber, $ξ1=ΔV(Ω−1/2)/vT1, ξ2=ΔV(Ω+1/2)/vT2, ΔV=V1−V2$ is the difference between the ion flow speeds, $Z′(ξ)=dZ/dξ$⁠, Z is the plasma dispersion function, $vTi=2kBTi/mi$ is the ion thermal speed, and λ_De is the electron Debye length. Here, the parameter Ω is defined by the substitution

in which ω  =  ω_R + iγ is the complex wave frequency. Here, γ = kΔVΩ_I is the growth rate, $ωR=12k(V1+V2)+kΔVΩR$ is the real frequency, and Ω_R and Ω_I are the real and imaginary parts of Ω, respectively.

Neutral damping of the ion-ion two-stream instability can be modeled using the BGK kinetic equation²¹ 

where ν_i–n is the ion-neutral collision frequency, $ni(1)$ is the first order density perturbation, $Φi(v)=mi/2πkBTi exp (−v2/vTi2)$ is a Maxwellian velocity distribution associated with stationary ions with respect to the neutral gas.

Again obtaining the linear dielectric response and assuming drifting Maxwellian distributions provides²⁴ 

where $ζ1=ΔV(Ω−1/2)/vT1+iν1−n/kvT1, ζ2=ΔV(Ω+1/2)/vT2+iν2−n/kvT2$⁠, and Ω is defined as in Eq. (6).

The conversion between the ion-neutral collision frequency, ν_i–n, and the neutral pressure used was $νi−n=ng(P) σ(vR) vR$⁠, where σ(v_R) is the velocity dependent cross-section, n_g(P) is the neutral gas density which is dependent upon the neutral pressure, P, and $vR=|vi−vn|$ with v_i the projectile ion speed and v_n the target neutral speed. The ion flow speed, V_i, is much greater than both the ion and neutral thermal speeds in the presheath, V_i ≫ v_Ti ≫ v_Tn. Therefore, the relative velocity between ions and neutrals is accurately approximated as v_R ≃ V_i.

Ion-neutral cross sections which did not have readily available data were modeled using the Langevin cross section²⁵ 

where α_p is the polarizability of the neutral, and m_R is the reduced mass, and ϵ_o is the permittivity of free space. Helium has a polarizability of $1.38 ao3$ and Xenon $27.06 ao3$⁠, where a_o is the Bohr radius. The Langevin cross section was used to model all the elastic scattering collisions between ions and neutrals. Two cross sections from LXcat for He⁺-He and Xe⁺-Xe charge exchange were used in the theoretical calculations.^26,27 The He⁺-Xe and Xe⁺-He charge exchange collisions were modeled to be half that of the corresponding elastic cross sections (Ref. 2, pg. 76). Figure 1 illustrates all the energy dependent cross sections used in simulations (to be presented in Sec. IV) and used to model the ion-neutral collision frequency in the theory (to be presented in Sec. III).

We also include the Krook collision model²¹ 

because it provides a commonly used estimate for neutral damping of waves. Applying the same calculation for the dielectric response function as in Sec. II B provides

The Krook collision operator predicts that the growth rate of instabilities is damped by subtracting the ion-neutral collision frequency: ω_K = ω_C – iν_i–n, where ω_K and ω_C are the dispersions calculated by the Krook and collisionless theory.

Direct numerical solutions for $ε̂BGK=0$ and $ε̂K=0$ were calculated for a Helium-Xenon plasma with T_e = 4.5 eV, T_i = 0.3 eV, n_He/n_e = 0.07, n_e = 7 × 10⁹ cm⁻³, and $ΔV=0.55|csHe−csXe|$⁠. These parameters are characteristic of the simulations presented in Sec. IV. Figures 2(a) and 2(d) illustrate the calculation technique. Contours identifying where the real part of the dielectric function is zero are shown in black and where the imaginary part of the dielectric function is zero are shown in magenta. The intersection of these contours represents the wave modes. The mode corresponding to the two-stream instability was chosen by identifying the root with the maximum growth rate and a real frequency between the Helium plasma frequency and the Xenon plasma frequency. The dispersion relation as kλ_De varies from 0.03 to 3.0 is shown as a multiple-colored line in Figs. 2(a) and 2(d).

Figure 2(b) shows the growth rate in the case of no ion-neutral collisions and the growth rate for ν_i–n = 0.032 c_s/λ_De, which is the smallest ion-neutral collision frequency required to completely damp the two stream wave predicted by the Krook model. Figure 2(e) shows the BGK predicted growth rate for an ion-neutral collision frequency of 0.032 and 0.35 c_s/λ_De. Figure 2(e) shows that for the BGK model a collision frequency of $νi−n=0.35 cs/λDe$ is required to damp the two-stream instability. This illustrates the general trend that the Krook model predicts a lower ion-neutral collision frequency cutoff for instability than the BGK. Figures 2(c) and 2(f) present the predicted real frequency dispersion for the two stream instability. No significant alterations to the real frequency of the wave occur in either model that could be testable in a PIC-MCC simulation.

The neutral collision frequency cutoffs for the Krook and BGK models convert to a neutral pressure of 6 mTorr and 64 mTorr, respectively. These cutoffs were calculated using $P=νi−n kB Tn/(σi−n vi)$ with σ_i–n = 5 × 10⁻¹⁵ cm², v_i = 0.55 c_s1, and T_n is the neutral gas temperature set to 300 K. The cross section σ_i–n is the total elastic cross section calculated from Eq. (9) for He⁺ at a speed v_i. This representative calculation shows that the Krook and BGK models differ by an order of magnitude when calculating the neutral pressure cutoff for the ion-ion two stream instability.

It is expected that the BGK dielectric should perform better than the Krook as the Krook collision operator is derived in the limit that the projectile and target particles have disparate masses. We have included the Krook predictions to highlight the deviations from BGK results.

The original IEF theory calculated V₁ and V₂ using Eqs. (1) and (3). The critical differential flow for the onset of the two stream instability, ΔV_c, in Eq. (3) was calculated from the collisionless dielectric response function, Eq. (5). Here, we extend the theory to include neutral collisions by calculating ΔV_c from the Krook dielectric, Eq. (11), and the BGK dielectric, Eq. (8).

We focus on a plasma with a concentration mix of 7% He⁺ and 93% Xe⁺ because this mix was shown in the previous work to produce strong instabilities.¹⁹ The simulations presented in Sec. IV were set to attempt to achieve the same mixture.

Numerical calculations of ΔV_c are presented in Fig. 3. These were created by iteratively calculating the growth rate of the two-stream wave in order to find the differential flow and pressures where the maximum growth rate is zero. The dotted line represents the difference between individual sound speeds in units of the He⁺ thermal velocity, which is the maximum energetically obtainable differential flow speed in the presheath. The threshold differential velocities were calculated for parameters relevant to the simulations presented in Sec. IV. Two solid black lines labeled (A) and (B) represent the differential flow boundary for instability predicted by the BGK model. These two contours represent the extremes of the plasma parameters observed in our simulations. Contour (A) was calculated using parameters T_i = 0.40 eV, T_e = 4.0 eV, n₁/n_e = 0.07, and n_e = 7 × 10⁹ cm⁻³ and contour (B): T_i = 0.46 eV, T_e = 3.7 eV, n₁/n_e = 0.035, and n_e = 1.5× 10¹⁰ cm⁻³. The average energy for ions in the presheath observed in our simulations was 1 eV. The total elastic cross sections for He⁺ and Xe⁺ at 1 eV were 4 × 10⁻¹⁵ cm² for both species, calculated using Eq. (9). This cross section approximated the average effect of ion-neutral collisions in the presheath. The dashed line is the differential velocity boundary calculated from the Krook model using the same parameters as contour (A).

The main result illustrated in Fig. 3 is that ion-neutral collisions act to increase the flow separation required to excite the instability. To calculate the individual flow speeds of each ion at the sheath edge, ΔV_c was taken from Fig. 3 at every neutral pressure and input into Eq. (1) and Eq. (3) to solve for V₁ and V₂. The predictions are displayed in Fig. 4 which shows that as neutral pressure increases ions transition from a speed influenced by IEF toward their individual sound speeds gradually over a decade in the neutral pressure.

Simulations were carried out using the PHOENIX1D code, which was also used in previous work to examine the two ion species sheath problem.¹⁹ PHOENIX1D is an electrostatic single spatial dimensional, three velocity dimensional PIC-MCC code. The domain was 10 cm long and sectioned into 1000 grid segments. Particles were generated throughout the domain by ionization reactions determined by a Monte Carlo electron-neutral collision routine.²⁸ Electrons were heated by a sinusoidal electric field perpendicular to the domain with a frequency of 13.56 MHz²⁹ (modeling an inductive radio-frequency discharge). The time step of the simulation was 1/2000 of the 13.56 MHz wave period or 37 ps. Particles were lost in the simulation upon traversing across grounded completely absorbing boundaries with no secondary electron emission. The gas constituents were Helium and Xenon. The code was run until a steady state was reached, where the rate of particle loss to the walls balanced the rate of particles gained by ionization. This equilibration time was approximately 200 μs. Data addressing the plasma dynamics were collected after this equilibration stage.

A Monte Carlo ion-neutral collision scheme was introduced that followed the method outlined in Vahedi and Surendra.²⁸ The ion-neutral elastic cross sections used in the simulations are displayed in Fig. 1. We made one alteration to the numerical scheme. A running total of the fractional number of ions undergoing a collision was kept. When the total was greater than one, a neutral collision was modeled and the counter lost one collision. This routine was made to capture the number of ion-neutral collisions when the neutral pressure was small and when the simulations were in the equilibration process.

Generating particles from ionization collisions caused the species temperature and density profiles to be dependent upon the neutral pressure (i.e., the electron temperature depends on pressure in the type of inductive discharge being simulated²). A singular energy dependent collision probability for electron-neutral collisions was established to make the sourcing rate of charged particles consistent for all simulations.

Figure 4 displays the results of our simulations overlaid with theoretical predictions for the ISS (dashed), the IEF-Krook (dashed-dot), and the IEF-BGK (solid-shaded) flow speeds of each ion species at the sheath edge as neutral pressure is varied over four orders of magnitude. Data points represent simulated flow speeds of He⁺ (blue) and Xe⁺ (red) at the sheath edge where charge neutrality was first violated by 1%, (n_He + n_Xe)/n_e − 1 = 0.01, when traversing from the bulk to the sheath. Error bars correspond to velocities measured when the quasi-neutrality condition was violated by 0.5% (low error bar) to 1.5% (high error bar).

The IEF-Krook and IEF-BGK predictions for the flow speed of both ions at the sheath edge in Fig. 4 were calculated from the thresholds shown in Fig. 3. Good agreement is observed between the IEF-BGK theory region and the simulated flow speeds, but not for the Krook model. This supports the expectation that the more complete BGK collision model leads to more accurate predictions for the flow speeds, and it highlights that substantial errors can arise when applying the simple, but crude, Krook model.

Figure 4 shows three regimes of neutral pressure when predicting the flow speed labeled IEF, IEF-BGK, and ISS. The IEF regime is the low pressure regime where the predicted flow speeds do not depend significantly on the neutral pressure. For the simulated plasma parameters this was for pressures less than 1 mTorr. The IEF-BGK regime is an intermediate pressure regime where the predicted flow speeds are well modeled using the BGK ion-neutral collision operator. This regime of pressures was recently tested experimentally.²² The Individual Sound Speeds (ISS) regime is the high neutral pressure regime where the two stream instability is predicted to shut off and the flow speed at the sheath edge is the individual sound speed of each ion species. For the simulation parameters, this occurred above 40 mTorr.

Electron density fluctuations were analyzed in the presheath for simulations representative of the IEF, IEF-BGK, and ISS pressure regimes. Figures 5(a)–5(c) display the Fourier transforms of electron density fluctuations, δn_e(k, ω), normalized to the steady state density, $neo(k)$⁠, over a segment of the domain ranging from 0.3 cm to 1.6 cm from the left wall for simulations corresponding to 0.02 mTorr, 20 mTorr, and 55 mTorr. The largest amplitude fluctuations were observed in the low pressure IEF regime, while virtually no fluctuations were observed in the high pressure ISS regime. The amplitude of fluctuations was observed to diminish continuously over the intermediate IEF-BGK pressure regime. Overlaid on each FFT spectrum is the real frequency and growth-rate calculated from Eq. (8) using the ions' temperature and density at 0.3 and 1.6 cm from the left wall. Good agreement is found between the fluctuations observed in the simulations and the dispersion calculated using local plasma parameters. Note that the differences in the predicted real frequency dispersion across the pressure regimes results from the local plasma conditions changing. The density fluctuation spectra provide evidence that support the ability of the BGK dielectric function to accurately quantify the damping of the two-stream wave.

In the sections above, we have shown that if the flow difference is above a critical threshold, ion-ion two-stream instabilities are excited. The ion flow speeds were also observed to merge when the instabilities were present. In this section, we now explicitly demonstrate that this merging occurs because of an enhanced ion-ion friction force that acts to speed up the Xe⁺ and slow down the He⁺.

Previous simulations¹⁹ had approximated the friction between species as the difference between the electric force density supplied by the potential gradient in the presheath and the kinetic force density of the particles. This section advances the simulation tests of the theory by removing this assumption and showing directly that the momentum exchange between ion species is due entirely to the friction force produced by the ion-ion two-stream instability.

Quantification of the friction force can be obtained from the momentum moment of the kinetic equation for the He⁺ species (species 1)

Here, C_1–n represents the ion-neutral collision operator associated with the MCC model described above. $C1−2=−(q1/m1)∂/∂vx·〈δf1δE〉$ represents the ion-ion collision operator associated with particles interacting with electric field fluctuations resolved by the grid (in particular, wave-particle interactions associated with the instability, but not standard Coulomb collisions, which are not resolved by the grid because they occur at spatial scales much less than the Debye length). Taking the momentum moment of Eq. (12) yields a time averaged momentum balance equation

where $R1−2=∫d3v m1vxC1−2(f1)$ is the instability enhanced frictional force density, and $R1−n=∫d3v m1vxC1−n(f1)$ is the ion neutral drag density. In Eq. (13), the kinetic force density of He⁺ has been label K, the electric force density E, and the difference has been labeled R for residual.

The instability enhanced friction force density, R_1–2, was calculated directly by writing the ion-ion collision operator in terms of its definition as the time correlation of the perturbations of the electric field and velocity gradient of the distribution function³⁰ 

which can then be integrated by parts to yield

where E(t) and n₁(t) are the electric field and density profiles at a specific time, t, $n1o$ and E^o are the steady state density and electric field profiles, and $〈…〉$ represents a time average.

Each of the four terms in Eq. (13) were calculated directly. Simulations were run until a steady-state was reached. The steady-state density n^o, potential ϕ^o, and second moment of the distribution function were obtained by averaging over 14 μs (200 rf periods). The steady state density and potential gradient yielded the electric term, E, while the gradient of the second moment of the distribution function yielded the kinetic term, K. R_1–2 was calculated using Eq. (15) at every grid cell then time averaging over 14 μs. The friction force density due to ion-neutral drag was directly calculated each time step by calculating the momentum of each particle before and after the ion-neutral collision routine. The particles' change in momentum due to neutral collisions was weighted on the grid and the drag density was determined at each grid point.

Figures 5(d)–5(f) displays the magnitude of each term, K (black), E (blue), R_1–n (green), R_1–2 (red), from Eq. (13) obtained from the simulations in a domain spanning from 0.3 cm to 1.6 cm from the left wall. At low pressure, the instability-enhanced friction force (red) is observed to be the dominant drag term balancing the electric field and kinetic force densities. However, as neutral pressure increases the neutral drag term (green) becomes important. At 20 mTorr each is approximately the same magnitude. At high pressure, the neutral drag becomes the dominant term balancing the electric field force.

The magnitude of the IEF term is observed to be highly correlated with the relative magnitude of instabilities observed in Figs. 5(a)–5(c). This further supports the IEF theory because the friction observed in the simulations is correlated with the growth rate of the instabilities. Furthermore, the observed decrease in the ion-ion friction force with the decrease in the growth rate as neutral pressure increases supports the prediction that neutral damping of two stream instabilities influences the ion flow speeds.

In the previous work,¹⁹ the magnitude of the instability enhanced friction term, R_1–i, in Eq. (13) was approximated by taking the difference of the kinetic and electric terms. The difference was denoted by the residual, R. In each plot in Figs. 5(d)–5(f), the cyan line represents residual term of the electric and kinetic terms. Figure 5(d) shows that the direct calculation of the instability enhanced friction term shown in red is consistent with the residual of the electric and kinetic terms shown in cyan, therefore showing that the momentum exchange between ions is entirely due to the instability enhanced friction at low pressure. The differences between the residual and the direct calculation can be attributed to statistical noise in the PIC-MCC simulations. At higher pressures, the residual consists of contributions from both IEF and ion-neutral drag. These results have verified that the momentum exchange between particles is exclusively explained by the instability-enhanced friction force from the ion-ion two stream instability at low pressure.

Simulations were conducted over a range of neutral pressures from 1 mTorr to 110 mTorr including both charge exchange and elastic ion-neutral collisions. In the simulations, the charge exchange cross sections for Xe⁺-Xe and He⁺-He were taken from LXcat, while those for Xe⁺-He and He⁺-Xe were taken to be half of the Langevin cross section. Again, the average ion energy was 1 eV. The total charge exchange cross section for He⁺ and Xe⁺ at 1 eV was 3.8 × 10⁻¹⁵ cm² and 9.0 × 10⁻¹⁵ cm². The total cross section (elastic plus charge exchange) was taken to be 1.3 × 10⁻¹⁴ cm² in the IEF-BGK theory.

Figure 6 presents the measured flow speed for He⁺ (blue) and Xe⁺ (red) ions at the sheath edge in the simulations. This shows that even with charge exchange collisions the IEF-BGK theory still accurately predicts the flow speeds below the neutral pressure threshold for instability. The plot displays three theoretical models. First, the ISS model for c_sHe and c_sXe is shown as dashed lines. Second, the IEF-BGK theory is shown in solid lines bounding a shaded region indicating the predicted range in flow velocities. Third, a collisional Bohm criterion is shown in a dashed-dot line for both species.

Examining Fig. 6, pressure regimes can again be assigned based upon the observed data as before. For our simulation parameters, the same three pressure regimes IEF (0–3 mTorr), IEF-BGK (3–13 mTorr), and ISS (13–20 mTorr) were observed. However, an additional pressure regime was observed at neutral pressures over 40 mTorr, where the speed of both ion species were observed to decrease below their respective sound speeds. This behavior violates the generalized Bohm criterion from Eq. (1). Results from a collisional Bohm criterion³¹ 

are displayed in Fig. 6. Here, $λi−n=(ngσi−n)−1$ is the ion-neutral mean free path and $λDes$ is the electron Debye length at the sheath edge. The data in Fig. 6 agree with the collisional Bohm criterion for pressures > 40 mTorr, but the onset of the collisional Bohm speed is not well modeled for lower pressures between 13 and 20 mTorr.

The IEF-BGK theory was observed to accurately predict the ion flow speeds up until the neutral pressure cutoff for instability. Above the neutral pressure cutoff for instabilities, the solution transitions to the traditional ISS, or the pressure becomes sufficiently high that the Bohm criterion becomes suspect, and modifications due to strong ion-neutral drag must be considered.

Helium and xenon ion velocity distribution functions (IVDFs) time averaged for 37 μs at the sheath edge are presented in Fig. 7 from simulations that included charge exchange collisions. Figure 7 shows that the He⁺ IVDF was observed to be double peaked when the two stream instability was present. This was an unexpected result as the two humped distribution had a sub sonic (relative to the helium sound speed) population that was predicted by the IEF-BGK theory, as well as a supersonic population. However, the flow moment of the distribution function produced flow speeds that agreed with the theoretical predictions.

Figure 8 shows the phase space distribution of helium ions in the 0.02 mTorr simulation from the midpoint of the domain to the right end wall that has been time averaged for 3 μs. At 3.5 cm from the midpoint of the simulation, a depleted region in velocity space forms in the plasma that expands all the way to the sheath edge located around 4.7 cm from the center. This depleted region is thought to be associated with an ion phase space hole generated by the two-stream instability. Figure 8 shows that helium ion holes propagate in the presheath and cause the helium IVDF to have two distinct populations corresponding to trapped and untrapped particles.

Holes in velocity space are an extensively studied³² phenomena and have been shown to be the steady state solution for two-stream unstable plasmas. Our simulations provide evidence that a large amplitude two-stream instability may generate these ion holes in the presheath and effectively create two distinct populations of ions. Figure 5(a) shows density fluctuations existing outside the wavenumber range predicted by the linear theory indicating that the instability has non-linear components that may be associated with ion velocity-space holes.

However, we note that this effect was not seen in IVDFs measured using LIF.¹⁴ The formation of ion holes in the simulations is possibly an artifact associated with having a single spatial dimension. It has been noted in the literature^33,34 that the generation of holes in velocity space in three spatial dimensions requires an anisotropy in the distribution function usually supplied by a magnetic field. Simulations of streaming instabilities in one dimension, however, have no such criterion and can generate phase-space holes readily. It is possible that the ion holes observed in our simulations are a result of our one dimensional setup and that the experiment in three dimensions is sufficiently isotropic that it precludes the generation of velocity-space holes. Further analysis will be needed to diagnose if ion holes can play a role in the 3D dynamics of the presheath of multiple ion species plasmas.

Since the simulations support the theoretical model, we applied the theory to predict the stability boundaries at other experimentally relevant conditions. Argon-Xenon is a common mixture used in experiments because both species velocity distribution functions can be measured using LIF. Figure 9 presents the differential flow velocity threshold of instability as a function of neutral pressure for an Ar-Xe plasma with plasma parameters: T_e = 1.0 eV, T_i = 0.04 eV, and σ_i–n = 1 × 10⁻¹⁴ cm². Figures 9(a) and 9(b) correspond to an Argon/Xenon mix of 50–50 and 10–90, respectively. The 50–50 mix has been traditionally examined experimentally in the two ion species sheath problem.^11–15,22

Both figures show how increasing the electron density increases the neutral pressure cutoff for the ion-ion two stream instability. The main result is that the two stream instability persists into the tens of mTorr range at low plasma densities (10⁸ cm⁻³) up until one hundred mTorr at a high plasma density (10¹⁰ cm⁻³). This suggests that ion-ion two stream instabilities could exist in higher pressure discharges than previous measurements (10s of mTorr). One caveat is that at such high pressures, some of the assumptions made start to break down, such as Maxwellian ion distributions, and the validity of the Bohm criterion, as shown in Fig. 6.

The IEF theory was extended to account for ion-neutral collisions. Tests using PIC-MCC simulations were found to agree well with the predicted flow speed of each ion species over a wide range of neutral pressures. Neutral damping of the two-stream instability was identified as the mechanism that causes the ion flow speeds to separate with increasing pressure, eventually reaching their individual sound speeds at high neutral pressure. The PIC-MCC simulations suggest this theory is accurate when charge exchange collisions are included, and therefore could be used in the interpretation of experimental diagnostics, the design of plasma-based materials processing tools, and in global plasma models,³⁵ all of which depend on the accurate modeling of ion flow speeds at the sheath edge. These results may also contribute to progress in modeling plasma-boundary interactions in plasmas with more than two ion species.³⁶ 

The authors thank the Iowa Center for Research by Undergraduates for sponsoring this work. This work was also supported in part by DOE Grant Award No. DE-SC0016473.