Ion dynamics in capacitively coupled argon-xenon discharges

An argon-xenon (Ar/Xe) plasma is used as a model system for complex plasmas. Based on this system, symmetric low-pressure capacitively coupled radio-frequency discharges are examined utilizing Particle-In-Cell/Monte Carlo Collisions (PIC/MCC) simulations. In addition to the simulation, an analytical energy balance model fed with the simulation data is applied to analyze the findings further. This work focuses on investigating the ion dynamics in a plasma with two ion species and a gas mixture as background. By varying the gas composition and driving voltage of the single-frequency discharge, fundamental mechanics of the discharge, such as the evolution of the plasma density and the energy dispersion, are discussed. Thereby, close attention is paid to these measures' influence on the ion energy distribution functions at the electrode surfaces. The results show that both the gas composition and the driving voltage can significantly impact the ion dynamics. The mixing ratio of argon to xenon allows for shifting the distribution function for one ion species from collisionless to collision dominated. The mixing ratio serves as a control parameter for the ion flux and the impingement energy of ions at the surfaces. Additionally, a synergy effect between the ionization of argon and the ionization of xenon is found and discussed.

Techniques using multiple driving frequencies, such as voltage waveform tailoring [12], succeed to independently control the plasma generation and the ion bombardment energy [17,18]. The plasmas investigated in these studies are predominantly argon plasmas [13,16,19]. However, industrially relevant etching plasmas consist of rather complex gas mixtures like CF 4 /H 2 [20][21][22] or SF 6 /O 2 [23,24]. For these plasmas, the interplay of several charged and neutral heavy species impacts the ion dynamics. The ion dynamics in the plasma eventually determine how ions reach the walls. Here, both the quantitative (e.g., how many ions reach the target/substrate?) and the qualitative perspective (e.g., how are the ions affected by collisions?) need to be considered.
Researching complex plasma chemistry in RF-CCPs is a tedious task. Experimental studies show that the ion energy distribution functions (IEDFs) at the electrodes become rather complicated [25][26][27][28][29]. Commonly used tools such as the retarding field analyzer filter the incident ions by energy and do not differentiate between the ion species [25,29,30]. There is recent and ongoing work to utilize ion mass spectrometry to overcome these issues [31]. Nevertheless, this technique is currently not widely applied as a diagnostic tool to analyze plasmas. Therefore, theoretical studies and simulations are necessary to help to interpret and to understand the measured data. However, the inherently complex chemistry renders a complete simulation cumbersome. The commonly used kinetic Particle-In-Cell/Monte Carlo collisions (PIC/MCC) method typically avoids complex chemistry mainly due to lack of cross section data (although conceptually possible). The reasoning is to keep the number of species and superparticles traceable [32]. Otherwise, the computational load of PIC/MCC simulation would not be feasible.
Combining complex discharge chemistry with the multi-frequency approaches mentioned above makes a detailed assessment of ion dynamics' features too cumbersome to conduct collectively. Hence, we decide to investigate the fundamental principles of a discharge with two ion species for this study. The mixture of the noble gases argon and xenon has some history of being an adequate model for complex chemistry. In lowpressure plasmas, the plasma chemistry of noble gases becomes relatively simple [35,36]. Therefore, studies on ion acoustic waves [33,34] and a generalized Bohm criterion [35][36][37] depicted this mixture as a simple example for a multi-ion discharge. Recent studies by Kim et al. [38] and Adrian et al. [39] contributed to those discussions using or referring to Ar/Xe plasmas.
Apart from being a model system, there are some academic applications of Ar/Xe plasmas (e.g., as trace gas for mass spectrometry [40,41], for the diagnostics of the electron temperature [42], or in halide lamp simulations [43]). Furthermore, the mixture has had great success as the illuminant [44] or as part of the illuminant mixture [32,45,46] of plasma display panels (PDPs). This historical background causes both gases to be relatively well researched. This fact entails many valuable data for theory and simulation.
This work aims to add to the existing studies conducted for various gas mixtures [47][48][49][50], investigating their intrinsic mixture dynamics. This knowledge will eventually enable the adaptation of the known means of plasma control to the complex discharges of industrial relevance. In contrast to the existing studies, our work is focused on the impact of the gas mixture composition on the ion dynamics. We will show that the gas composition is a suitable control parameter for the ion dynamics (e.g., the impingement energy of ions at the surface).
This manuscript presents our findings as follows: in section 2, we introduce our simulation framework and a model for the missing cross section data. Moreover, we introduce an energy balance model for CCPs with multiple ion species. The findings of these models are interpreted in section 3. We first discuss the influence of a variation of the gas composition on the ion dynamics. Then, we validate the energy balance model with our simulation data. Afterward, we apply this energy balance model to support and to analyze our gas composition variation findings. We conclude section 3 by discussing and examining the influence of a variation of the gas composition combined with a variation of the driving voltage on the ion dynamics. Finally, in section 4, we summarise our findings, draw a conclusion, and set this work into the context of industrial applications.

Particle-In-Cell Simulation
The first particle simulations were introduced in the 1940s [51], and the PIC/MCC scheme was developed in the 1960s [52]. Since then, the PIC/MCC method became a commonly used tool to self-consistently simulate low-pressure plasmas [32,52,53]. Despite having the disadvantage of a substantial computational load, its most significant advantage is the statistical representation of distribution functions in phase-space, allowing the method to capture non-local dynamics [53,54].
For this work, a benchmarked PIC/MCC implementation called yapic1D [56] is used to generate the results. The original code is modified to include two background gases and multiple ion species. Aside from that, diagnostics for the energy balance model mentioned above is added to the original code.
This simulation setup is taken to be fully geometrically symmetric (compare Wilczek et al. for details [54]). 1d3v electrostatic simulations are executed using a Cartesian grid with 800 grid cells representing an electrode gap of 25 mm. The resulting cell size ∆x meets the requirement to resolve the Debye length λ D [52,54,56]. Similarly, the single harmonic driving frequency f RF " 13.56 MHz is sampled with 3000 points per RF period. The time step ∆t is sufficiently small to fulfill the requirement regarding the electron plasma frequency ω pe [52,54,56]. Several other studies mention the influence of the number of superparticles on the statistics and the plasma density [56][57][58]. For this work, we did not include individual weighting for different particle species. To have an acceptable resolution for each ion species at all values of the xenon fraction x Xe , we simulated about 800.000 super-electrons for each case. The advantage of this choice is that an average of 3000 converged RF-cycles provides satisfactory results.
The ideal gas law defines the neutral species' total density, and the neutral fraction x i is varied. Thereby, the gas pressure p gas is kept constant at 3 Pa, and the gas temperature T gas at 300 K. First, we choose the amplitude of the RF voltage V RF to be  Table 1. Plasma chemistry and collision processes considered in the simulation.
Meaning of the data sources: "Phelps" refers to the cross section data found initially in the JILA database [59] and now distributed by the LXCat project [60][61][62]. "LJ pot" refers to a cross section obtained based on a phenomenological Lennard-Jones potential as described by Laricchiuta et al. [63]. "Viehland" marks a cross section calculated from an interaction potential given by Viehland et al. [67]. Details of the calculations can be found in section 2.3.
100 V. Later in section 3.3, we discuss the implications of a voltage variation between 100 V and 1000 V on the ion dynamics. All the parameters presented in this section are typical for baseline studies of RF-CCPs [14-16, 54, 56].

Discharge chemistry
For PIC simulations to provide a realistic representation of the particle distribution functions and physics in a low-pressure discharge, collisions need to be considered. The method of choice is the Monte Carlo collision technique [52][53][54] that is combined with a so-called null collision scheme [54][55][56]. Both techniques require the knowledge of momentum transfer cross sections. The chemistry set for argon and xenon is in line with the work of Gumundsson et al. [35,36]. All reactions can be seen in detail in table 1. In contrast to Gumundsson et al., we decide to take advantage of the commonly used and acknowledged [19,54,56,68] cross section data obtained by Phelps. The data was initially distributed via the JILA database [59] and is now available at the LXCat project website [60][61][62]. Phelps combines the cross sections for all electronically excited states into one "effective excitation" cross section. This effective excitation reduces the total number of reactions and the numerical load.
The second difference compared to Gumundsson et al. is our treatment of the missing cross section data for Ar`/Xe and Xe`/Ar. Both conclude to neglect charge transfer collisions between argon and xenon due to it being a non-resonant process that requires a third particle to ensure momentum and energy conservation. The disparity lies in our treatment of the remaining scattering process. They assume the cross sections for processes 7 and 9 and the cross sections for processes 10 and 12, respectively, to be equal. We adopt a physical model to procure the necessary cross sections from interaction potentials. In this way, we create individual cross sections for processes 9 and 12. The details of how to calculate these cross sections will be presented in the following section.
The cross section data used for this work are depicted in figure 1. It is noticeable that the cross sections of processes involving xenon species generally have higher values than corresponding processes that involve argon species. In terms of a hard-sphere model [1], this deviation is explained by the different covalent atom radii of argon and xenon [69]. Xenon is, compared to argon, simply the bigger target. In terms of a more sophisticated collision model [1], one, for example, needs to consider the atomic polarisability of the neutral particle. Nevertheless, such a view leads to the same insight. Xenon has a higher atomic polarisability than argon [69] and stronger interaction with charged particles. Correspondingly, the cross section for charged particles interacting with xenon has to be larger than the cross section for the interaction of charged particles and argon.

The calculation of the cross sections
On an elementary level of theory, all cross sections are based on an interaction potential between the colliding particles. If the literature does not provide a cross section, a possible solution is to make it the modeler's task to develop an interaction potential by making several assumptions. A classic example of this is the Langevin capture cross section [72] used in studies to make up for unknown cross sections [73]. Despite the Langevin cross section's advantages, a complete implementation is numerically extensive and leads to anisotropic scattering [74,75]. The cross sections given by Phelps are a kind of momentum transfer cross sections [59]. There, the scattering angles are found in an isotropic manner [52]. Hence, it is questionable to apply anisotropic scattering for two collisions while the other collisions are treated isotropically. We perceive another approach to be more suitable for this work.
The approach used in this work is based on Laricchiuta et al. [63], who use a phenomenological potential to describe a two-body interaction given by where the standard exponents of the Lennard Jones potential, 12 and 6, are replaced by npx ij q and m. Depending on the type of interaction, m is either 4 for neutral-ion interactions or 6 for neutral-neutral interactions. In this work, the potential is applied to neutral-ion interactions only. Hence, m is always equal to 4. The dimensionless coordinate x " r{r m,ij depends on the parameterized position of the potential well The data source and a detailed description of each process are found in table 1. Abbreviations used in the legend: ela = elastic collision electron/neutral, exc = electronic excitation electron/neutral, ion = electron impact ionization, iso = isotropic scattering ion/neutral as defined by [59], back = backscattering ion/neutral as defined by [59].
r p,ij . The potential itself is scaled by the parameterized potential well depth p,ij . Both parameterizations are empirical approximations that depend on atomic properties like the polarizability. More details related to the exact empirical formulas can be found in Laricchiuta et al. [63], Cambi et al. [64], Cappelletti et al. [66], and Aquilanti et al. [65]. Two additional steps are required to obtain the cross section. The first step is calculating the scattering angle, χ ij according to with ij the kinetic energy in the center of mass frame, b the impact parameter, r the distance between the particles, and r 0 the distance of closest approach. The scattering angles are calculated using a program that is based on Colonna et al. [71]. The second step is calculating the cross section σ ij σ plq ij p ij q " 2π with l an integer that indicates which type of cross section is calculated. In this work, we used l " 1, which corresponds to the momentum transfer cross section. The cross section is integrated based on an algorithm developed by Viehland [70]. Finally, the scattering angle corresponding to the obtained momentum transfer cross sections is consistently taken to be isotropic in our simulations.

Energy balance model
The conservation of energy is one of the central continuity equations of physics and so knowing how the energy disperses into a system is key to understanding the process. In terms of low-temperature plasma physics, a frequently used model as given by Lieberman and Lichtenberg [1] for a geometrically symmetric situation reads: S abs denotes the total energy flux into the system, n s the plasma density at the sheath edge, u B denotes the Bohm velocity, Γ B is the ion flux at the Bohm point, and ε tot is the total energy loss in eV. The last transformation in equation (4) shows that the energy loss per electron-ion pair created may be split into an energy loss due to electrons hitting the bounding surface (ε e ), an energy loss due to collisions (ε c ), and an energy loss due to ions impinging at the bounding surface (ε i ). The loss terms ε e and ε i describe an averaged energy loss of the system per lost particle (neglecting particle reflections). The third term ε c is treated differently. It represents the collisional losses per newly created electron/ion pair. Previous work [76] has shown that an adaptation of equation (4) gives insight into the system's electron dynamics by calculating all necessary terms from a PIC/MCC simulation. An essential insight is that, due to flux conservation, the Bohm flux Γ B can be exchanged by the electron flux Γ e,el or the ion flux Γ i,el at the electrode.
In detail, the energy conversion through collisions ε c consists of an electron ε c,e and an ion contribution ε c,i . For low-pressure plasmas, it is argued that the energy loss due to ion collisions ε c,i is often negligible [1]. However, a PIC/MCC study by Jiang et al. [77] showed that ε c,i can significantly impact the energy balance of low-pressure plasmas.
For this and more complex systems, it is useful to split the total energy flux S abs,tot into a separate term for each species. This separation is done in equation (5). Besides, we split the collisional losses to the background gas for ions, ε c,i , into two terms. One represents the losses due to charge exchange collisions for Ar + ions (ε cx,Ar`) and Xe + ions (ε cx,Xe`) . The other term gives the losses caused by the remaining isotropic scattering. It separates the isotropic losses for Ar + ions (ε is,Ar`) , and Xe + ions (ε is,Xe`) . This distinction is based on the nomenclature of Phelps [59] and will prove useful for understanding the ion dynamics. The terms for the electron flux (Γ e ), the Ar + ion flux (Γ Ar`) and the Xe + ion flux (Γ Xe`) in this model are obtained from the PIC/MCC simulation at the surface of the electrodes.

Influence of the neutral gas composition on the discharge
Initially, the influence of the gas composition on the discharge is investigated by varying the Ar/Xe density ratio. Figure 2 a) shows the total ion density n i,tot as a function of the xenon gas fraction x Xe or the argon gas fraction x Ar , respectively. Here, the total ion density n i,tot is defined as the sum of the spatially and temporally averaged number densities of Ar + and Xe + ions. The gas fractions of argon and xenon are defined by the ratio of the respective species density and the total gas density. In analogy to this definition, we define an ion fraction, e.g., the fraction of Xe + ions x Xe`, as the ratio of the number density of Xe + ions and the total ion density n i,tot . Figure 2 b) depicts this Xe + ion ratio as a function of the xenon gas fraction x Xe or argon gas fraction x Ar , respectively.
When varying the gas mixture from pure argon to pure xenon by successively increasing the xenon fraction x Xe , the plasma density rises significantly over about one order of magnitude ( fig. 2 a)). The ratio of Xe + ions ( fig. 2 b)) reveals that even small admixtures of xenon to an argon gas produce a high amount of Xe + ions. A xenon fraction of x Xe « 0.15 is already sufficient for Xe + ions to become the dominant ion species. Xenon admixtures of about 30 percent (x Xe " 0.3) produce a strongly Xe + dominated discharge (x Xe`Á 0.8). Both the development of the plasma density and the fraction of Xe + ions in the discharge as a function of the gas composition show non-linear relations. Hereby, the trend of figure 2 a) approximates a compressed parabola whilst the trend of 2 b) resembles the function of the square root. In the following, the overall dominance of Xe + ions will be examined and explained in more detail. The difference in the ionization energies gives a basic explanation of the observed behavior. The ionization threshold for xenon (ε thr,i,Xe " 12.12 eV) is much smaller than the threshold for argon (ε thr,i,Ar " 15.8 eV) (comp. tab. 1). This disparity allows lower energetic electrons to ionize xenon in contrast to argon. Additionally, the ionization cross section of xenon σ i,Xe is about one order of magnitude bigger than the corresponding cross section σ i,Ar for argon (comp. fig. 1 a) and b)). As a result, Xe + ions are prevalent, even for low xenon admixtures, and dominate the discharge for a wide mixture range. This result agrees with previous works [48,50] that, for different mixtures, have shown at least one dominant ion species for a wide range of admixtures.
The influence of the gas composition also directly manifests in a variation of the IEDFs for both ion species. The plots of figure 3 show IEDFs for both Ar + and Xe + ions at the electrode surface. The energy, plotted on the abscissa, is given in eV. The ordinate shows the IEDF normed on the respective ion flux Γ i,s at the electrode. Each row of figure 3 represents results for both ion species and the same case. The cases are distinguished by the xenon fraction x Xe as indicated. Here, the plots in the right column show IEDFs of Ar + ions, and the results for Xe + ions are shown in the right column.
In section 2.2, we argue that the charge exchange between Ar + /Xe and Xe + /Ar, respectively, is a non-resonant process and a three-body collision. We conclude that this process is negligible. As a result, a variation of the gas composition changes the ions' probability to perform charge exchange collisions. Therefore, Ar + ions, for high argon fractions x Ar , show an IEDF clearly dominated by collisions. This IEDF becomes a collisionless distribution for small argon admixtures to a xenon background ( fig. 3 left). The IEDF of Xe + ions shows a similar trend except that Xe + ions have a less distinct  bimodal behaviour for the cases with high argon fraction. This difference is explained by the scaling of the width of the bimodal peak being proportional to a m´1 i [2,28]. Besides, an argon fraction x Ar of 0.2 and 0.3 or, vice versa, a xenon fraction x Xe of 0.2 and 0.3 creates an intermediate or hybrid regime. A significant number of ions experiences the discharge as being collision dominated, while the remaining ions cross the sheath collisionlessly. In figure 3, the described regime is visible for argon at x Xe " 0.2 and xenon at x Xe " 0.8. Several distinct peaks are visible at low energies that stem from charge exchange collisions, and at high energies, the characteristic collisionless bimodal peak is clearly established. In these cases, particularly, the scaling of the bimodal peak width can be observed. For both cases, Xe + ions establish a bimodal peak narrower than the bimodal peak formed by Ar + ions.

Revision and analysis of the energy balance model
For a fundamental understanding of the energy distribution within the system, the energy balance model resembled by eqs. (5) -(8) may be used, as shown in figure  4. We calculate all parameters and properties by means of a PIC/MCC simulation averaged over 3000 RF periods. The plot shows two bars for each of the chosen gas compositions. The grey bar on the left-hand side represents the total absorbed energy flux S abs,tot . The colored bars on the right-hand side resolve the different channels of energy dissemination in detail. The colors blue (electron energy lost at the electrode 9 ε e ), red (averaged energy consumption per e/Ar + -pair 9 ε c,e,Ar ), and green (averaged energy consumption per e/Xe + -pair 9 ε c,e,Xe ) represent the right-hand side of equation (6). The right-hand side of equation (7) is depicted in pink (Ar + ion energy loss at the electrode 9 ε i,Ar`) , cyan (energy loss by isotropic scattering 9 ε is,Ar`) , and purple (energy loss by backscattering 9 ε cx,Ar`) . The remaining colors olive (Xe + ion energy loss at the electrode 9 ε i,Xe`) , brown (energy loss by isotropic scattering 9 ε is,Xe`) , and orange (energy loss by backscattering 9 ε cx,Xe`) visualize the right-hand-side of equation (8).
At first, it is noticeable that figure 4 shows a roughly square-root-shaped increase of the absorbed energy flux density S abs as a function of the xenon fraction x Xe . This trend is a consequence of the boundary conditions in combination with the varied gas composition. The PIC/MCC simulations considered in this work use a singlefrequency voltage source as a boundary condition for calculating the electric field. The energy flux density is calculated self-consistently according to the plasma state. At low xenon fractions x Xe , xenon neutrals and Xe + ions successively provide additional loss mechanisms, and the energy consumption increases rapidly. Whereas at higher xenon fractions, xenon already dominates the discharge, and the energy consumption slowly saturates. Lieberman and Lichtenberg present the scaling law n s 9 S abs [1]. In section 3.1, we discussed that the trend of the plasma density n i,tot as a function of the xenon fraction x Xe (comp. fig. 2 a)) is approximated by a parabola. Combined with the square-root-shaped trend of the absorbed energy flux density S abs as a function of the xenon fraction x Xe , we see the resulting trend of n i,tot and S abs match the anticipated scaling.
The results calculated for pure argon (x Xe " 0.0) and pure xenon (x Xe " 1.0) discharges resemble the classical model given by equation (4). The results demonstrate that all individual loss terms sum up to the total energy flux and thus prove the models' exact energy conservation. Both the argon case and xenon case reveal that the energy loss due to colliding ions (argon: cyan and purple, xenon: brown and orange) has a significant contribution to the energy balance (argon: « 31.1 % of the total energy, xenon: « 35.6 %). These findings are similar to the study of Jiang et al. [77].
The remaining bars of figure 4 review the modified energy balance model presented in equation (5) to (8). It shows, for some exemplary gas mixtures, that the suggested balance for multiple ion species is complete and that each species' energy transfer can be traced individually. Furthermore, the results show that for a complete energy balance for plasmas with two ion species, colliding ions' energy transfers are at least as important as they are in mono ionic plasmas [77]. Especially, the energy losses due to charge exchange collisions (ε cx,Ar`-purple-or ε cx,Xe`-orange-, resp.) make up for a significant amount of the transferred energy.
Both the individual energy transfers of each particle species and the exact resolution of specific loss channels will in the following prove useful to understand and analyze the discharge. To make the results comparable, we decide to switch the representation of the energy flux density S abs from absolute units to relative units (comp. fig. 5). Thereby, Figure 5. The energy balance equations (5) -(8) applied for the background gas variation. All properties are calculated from a PIC simulation and referred to the total absorbed energy flux S abs,tot . All plots show the right-hand side of their corresponding equation in relative units. a) represents equation (5), b) equation (6), c) equation (7) and d) equation (8). (conditions: p gas " 3 Pa, l gap " 25 mm, V RF " 100 V, f RF " 13.56 MHz) we refer to the energy fluxes of each case individually with respect to the total energy flux S abs,tot . Figure 5 shows a rearrangement of the data of figure 4 in the relative representation explained before. Each of the subplots a) to d) respectively present the right-hand side of equations (5) to (8). The abscissa of all plots mark energy flux densities in relative units, and the ordinates are in units of the gas fractions (x Xe or x Ar , resp.). The color schemes for figures 5 b), c), and d) are similar to the ones used in figure 4. Figure 5 a) introduces a new color scheme for the total energy fluxes absorbed by electrons (bright blue), Ar + ions (fuchsia), and Xe + ions (lime green).
In section 3.1, we point out two observations. First, Xe + ions are for a wide range of mixtures the dominant ion species. Second, for constant gas pressure, collisional features of the IEDF depend on the gas composition, and even a collisional/collisionless hybrid regime can be reached. Both observations are confirmed and explained by the energy balance. Figure 5 a) shows that for a xenon fraction x Xe between 0.15 and 0.2, Ar + ions and Xe + ions absorb an equal amount of energy (30 % of S abs or « 3 W/m 2 ). Simultaneously, the production of Xe + ions is more effective than the production of Ar + ions. This increased effectiveness is due to the lower ionization energy of xenon (ε thr,Xe " 12.12 eV) compared to argon's ionization energy (ε thr,Ar " 15.8 eV).
The case for a xenon admixture of 20 percent (x Xe " 0.2) serves as the best example for this finding. There are, with a Xe + ion fraction x Xe`« 0.7 (comp. fig. 2), more Xe + ions than Ar + ions inside the discharge. Nevertheless, more energy per electron-ion pair is consumed to produce Ar + ions (red) than for the generation of Xe + ions (green) ( fig.  5 b)). This finding is explained by the lower excitation and ionization levels of xenon compared to argon. Simultaneously, these lower excitation and ionization levels open up new loss channels for the electrons inside the system. Raising the xenon fraction x Xe yields more and more electrons that are not energetic enough to participate in inelastic processes in an argon discharge. Thus, the averaged electron energy ε e drops when going from an argon discharge to a xenon discharge. The decreasing loss term ε e (blue) in figure 5 b) hints at the average electron energy of the system and gives evidence of this explanation. All in all, this shows that the production of Xe + ions fills an unoccupied energetic niche where numerous low energetic electrons can participate. Therefore, a significant production of Xe + ions is observed even for low xenon fractions x Xe and Xe + ions are the dominant ion species for the majority of the possible Ar/Xe mixtures.
The trends observed in the IEDFs ( fig. 3) and the conclusions drawn from this observation are confirmed by the energy balance as well ( fig. 5 c) and d)). Looking at the losses due to charge exchange collisions ε cx,i for both Ar + ions ( fig. 5 c), purple) and Xe + ions ( fig. 5 d), orange), it becomes apparent that the collisional features are switched between Ar + and Xe + ions when going towards more argon, or xenon respectively, dominated gas mixtures. The losses due to charge exchange for Ar + ions ε cx,Ar`m onotonically fall as a function of the xenon fraction x Xe ( fig. 5 c), purple) while the corresponding term for Xe + ions ε cx,Xe`m onotonically raises, when displayed as the same relation ( fig. 5 d), orange). The slight difference in the trends is explained by the dominance of the Xe + ions in the discharge. While the density of Xe + ions rapidly increases, when adding small amounts of xenon to an argon background (comp. fig.  2 a)), the density of Ar + ions vanishes as fast among the dominant Xe + ions. Hence, there are not enough Ar + ions present in discharges dominated by Xe + ions, so that the losses of Ar + ions in total cannot significantly contribute to the energy absorbed by the discharge (fig. 5 a), fuchsia).
In addition to this, the mean energy of Xe + ions at the electrode ε i,Xe`s hows a very different trend than all the collisional quantities ( fig. in 5 d), olive). Instead of monotonically rising with the xenon fraction x Xe as the corresponding Ar + term does as a function of the argon ratio x Ar (comp. fig. 5 c), pink), the Xe + curve shows a maximum at x Xe " 0.4. This maximum is closely connected to the dominance of Xe + ions. At 40 percent xenon admixture (x Xe " 0.4), Xe + ions already make up for about 90 percent of the ions in the discharge ( fig. 2 a)). At the same time, argon is the dominant background rendering Xe + ions more or less incapable of doing a relevant amount of charge exchange collisions. This lack of charge exchange collisions is seen in the IEDF of Xe + ions, that even for a xenon fraction x Xe " 0.5 shows a characteristic collisionless single bimodal peak ( fig. 3, right). For lower xenon fractions x Xe , the number density n Xe`a nd the flux density Γ Xe`a re lower and fewer Xe + ions reach the electrode. This decrease results in a lower energy loss. For higher xenon fractions x Xe , the charge transfer collision of Xe/Xe + becomes more and more probable. This trend manifests in the IEDFs ( fig. 3, right) and the trend of the loss term for charge exchange ε cx,Xe`( fig. 5 d), orange). Thus, the energy loss of Xe + ions to the surface finally drops because the energy gets dissipated more strongly to the neutral gas via charge exchange collisions.
The minimum of the total energy flux density absorbed by electrons S abs,e ( fig. 5 a), bright blue) has a similar explanation. For a xenon fraction x Xe " 0.5, electrons absorb the lowest amount of energy. Under these conditions, Xe + ions make up for almost all the ions in the discharge. Figure 2 b) shows that for a xenon fraction x Xe " 0.5 the Xe + ion fraction x Xe`i s approximately 0.9. At the same time, xenon atoms make up for just 50 percent of the background gas. The amount of collisions with argon or xenon particles respectively is, as argued before, significantly reduced compared to mixtures with a high amount of either of the gases. Thus, for xenon fractions x Xe ă 0.5, the production of Ar + ions causes electrons to absorb and invest more energy. For xenon fraction x Xe ą 0.5, collisions with xenon neutrals become successively more probable, and the production of Xe + ions consumes more energy (comp. fig. 5 b), green) without significantly changing the discharge conditions any more (comp. fig. 2). Additionally, figure 3 shows that x Xe " 0.5 is optimal for producing high energetic ions. Both ion species establish the characteristic collisionless bimodal peaks and impact the surface with high energies. Therefore, the relative amount of energy brought by ions to the surface is maximal. For lower xenon fraction (x Xe ă 0.5), the IEDF of Ar + ions is visibly affected by collisions and vice versa for higher xenon fraction (x Xe ą 0.5).

Influence of the driving voltage on the discharge
In terms of our simulation, a raised driving voltage equals, if all other parameters (gas composition, pressure, etc.) are kept constant, raising energy input to the system. Figure 6 a) shows a semi-logarithmic representation of the time and space averaged total plasma density n i,tot as a function of the gas fractions (x Xe or x Ar , resp.). The different colors differentiate the data for different RF amplitudes (black = 100 V, red = 250 V, blue = 500 V, green = 1 kV). The black curve shows the same data as figure 2 a). Due to the aforementioned higher input energy, the plasma density is raised in general while the several curves' general trend is kept. Independent of the driving voltage, argon discharges have a significantly lower plasma density than xenon discharges, and the transition while varying the gas composition shows the same non-linear trend. In sections 3.1 and 3.2, we discuss that in this context, non-linear means parabolic.
Apart from this, a varied driving voltage alters the dominance of Xe + ions. Figure  6 b) shows a similar plot to figure 2 b). The Xe + ion fraction is presented as a function of the gas fraction (x Xe or x Ar , resp.). The colors have the same meanings as in figure 6 a), and the black curve was also presented before (see fig. 2 b)). Figure 6 b) shows that, for a fixed xenon fraction x Xe , a raised voltage reduces the fraction of Xe + ions x Xep resent in the discharge. The case for x Xe " 0.2 is a good example of this observation. When increasing the driving voltage from 100 V to 1 kV, the ratio of Xe + ions x Xed rops from approximately 0.7 to roughly 0.6. Once again, the energy balance ( fig. 7) explains the discharge mechanisms governing how an increased driving voltage raises the plasma density. Similar to figure 5, terms on the right-hand side of the energy balance equations (5) -(8) are shown in relative units and as a function of the gas fractions (x Xe or x Ar , resp.). In contrast to figure 5, each panel of figure 7 represents just one term of the respective equation's right-hand side. The different curves represent data for different driving voltages V RF , ranging from V RF " 100 V to V RF " 1000 V. The color scheme is analog to figures 4 and 5. Figure 7 a) shows the total energy flux density absorbed by electrons (S abs,e ) in bright blue. Figure  7 b) depicts the total energy flux density absorbed by Ar + ions (S abs,Ar`) in fuchsia, and figure 7 c) presents the total energy flux density absorbed by Xe + ions (S abs,Xe`) in lime green. Together figures 7 a) -c) show the right-hand side of equation (5). Therefore, Figure 7. The energy balance equations (5) -(8) applied for both the variation of the background gas and the driving voltage V RF . All properties are calculated from a PIC simulation and referred to the total absorbed energy flux S abs,tot . All plots show one term on the right-hand side of their corresponding equation in relative units. a) shows the electron's part S abs,e of eq. (5). a 1 ) -a 3 ) show the three terms of equation (6) and add up to the respective curve of a). b) represents the Ar + ions' part S abs,Arò f equation (5). b 1 ) -b 3 ) present the three terms of equation (7) and sum up to the respective curve of b). c) shows the Xe + ions' part S abs,Xe`o f equation (5). c 1 ) -c 3 ) depict the three terms of equation (8), and their addition gives the respective curve of c). (conditions: p gas " 3 Pa, l gap " 25 mm, f RF " 13.56 MHz) the corresponding data points horizontally always add up to give 100 % (or the total energy flux density S abs,tot , resp.). Vertically the details of each particle species' power absorption are presented. Figures 7 a 1 ) -a 3 ) each show one term of the right-hand side of equation (6). The average energy loss of electrons at the electrodes ε e is shown in figure 7 a 1 ) in blue. The averaged amount of energy needed to create an electron/Ar + ion pair (ε c,e,Ar ) is found in panel a 2 ) in red, and the related term for electron/Xe + ion pairs (ε c,e,Xe ) is depicted in panel a 3 ) in green. The individual terms of the right-hand side of equation (7) are shown in figures 7 b 1 ) -b 3 ). They reveal the details of the Ar + ion dynamics by presenting the average energy loss by Ar + ions at the electrodes (ε i,Ar`, fig. 7 b 1 ), pink), the energy loss of Ar + ions caused by isotropic scattering (ε is,Ar`, fig.  7 b 2 ), cyan), and the energy loss of Ar + ions due to backscattering (ε cx,Ar`, fig. 7 b 3 ), purple). Similarly, figures 7 c 1 ) -c 3 ) show the right-hand side of equation (8). They unravel the details of the Xe + ion dynamics by showing the average impingement energy of Xe + ions at the electrodes (ε i,Xe`, fig. 7 c 1 ), olive), the energy lost by Xe + ions in isotropic scattering collisions (ε is,Xe`, fig. 7 c 2 ), brown), and the energy lost by Xe + ion in backscattering collisions (ε cx,Xe`, fig. 7 c 3 ), orange). Vertically, the sum of the data in the subscript labeled panels gives the curves of the non-subscript labeled one (e.g., panels a 1 ) -a 3 ) sum-up to panel a)).
In general, it is apparent that a raised driving voltage reduces the ratio of energy coupled to the electrons ( fig. 7 a)) and raises the fraction absorbed by both Ar + and Xe + ions ( fig. 7 b) or fig. 7 c), resp.). The increased energy consumption into the ion contribution mainly consists of two parts. First, a raised driving voltage V RF increases the voltage drop across the boundary sheaths, and ions gain higher impingement energies after crossing the sheath collisionlessly. This is shown in figure 7 b 1 ) for Ar + ions and in figure 7 c 1 ) for Xe + ions. Second, an increased energy gain for the ions inside the sheath goes along with an increased energy loss caused by charge exchange collisions. The corresponding terms ε cx,Ar`f or Ar + ions ( fig. 7 b 3 )) and ε cx,Xe`f or Xe + ions ( fig.  7 c 3 )) support this hypothesis. Furthermore, the cross sections for charge exchange dominate the ones for isotropic scattering at high energies (comp. fig. 1 c) and d)). Correspondingly, the already low energy losses by Ar + ions (ε is,Ar`, fig. 7 b 2 )) and Xe + ions (ε is,Xe`, fig. 7 c 2 )) caused by isotropic scattering decrease due to the increased driving voltage V RF . The maximum of figure 7 c 1 ) was discussed in section 3.2. The energy-efficient production of Xe + ions already creates a high amount of Xe + ions for small xenon fractions x Xe . Thus, there are optimal parameters for Xe + ions to bombard the surface with the least collisional loss (x Xe " 0.4 for V RF " 100 V, sec. 3.2). The aforementioned enhanced role of backscattering and decreased influence of isotropic scattering causes the optimal parameters for higher driving voltages V RF to shift to higher xenon fractions x Xe (e.g., X Xe " 0.5 for V RF " 1000 V, fig. 7 c 1 )).
In terms of ion production, the previous assessment shows that the higher the driving voltage is set, the smaller the fraction of the energy consumed for creating new electron/ion pairs becomes. Furthermore, the maximal amount of energy consumed for creating Xe + ions in a pure xenon background (x Xe " 1.0, fig. 7 a 3 )) is always lower than the corresponding maximum for Ar + ions in a pure argon background (x Ar " 1.0, fig. 7 a 2 )). As argued before, this finding correlates with the fact that the threshold energies of all inelastic processes involving xenon are significantly lower than those involving argon. This observation additionally reveals why Xe + ions dominate the discharge for most conditions. It becomes best visible by comparing the pure argon case (x Xe " 0.0, fig. 7 a 2 )) with the pure xenon case (x Xe " 1.0, fig. 7 a 3 )) for a driving voltage of 1 kV. Here, roughly eight percent of the total energy flux density is used to produce an electron/Ar + ion pair ( fig. 7 a 2 )). For the corresponding pure xenon case, only five percent of the total energy is used to produce an electron/Xe + ion pair ( fig. 7 a 3 )). Simultaneously, the xenon case's plasma density is more than one order of magnitude higher than in the argon case (comp. fig. 6 a)).
Since the production of Xe + ions remains more effective for all applied driving voltages, there has to be another reason why the dominance of Xe + ions is reduced. A close examination of figures 7 a 1 ) and 8 explains the observed. Both panels of figure  8 are similar in structure to figure 6 a), but show the individual ion densities (n Ar`i n fig. 8 a) or n Xe`i n fig. 8 b), resp.) as a function of the gas fraction (x Xe or x Ar , resp.). The different colors again mark different values of the driving voltage V RF , and the color scheme is the same as in figure 6 a). The trend of the Ar + ion density in figure 8 a) already reveals the underlying process responsible for the decreased dominance of Xe + ions. Even for the base case (V RF " 100 V), the maximum of the density of Ar + ions is found for a xenon fraction x Xe " 0.2 and not for a xenon fraction x Xe " 0.0 as it is vice versa the case for Xe + ions (comp. fig. 8). This maximum is shifted by a raised driving voltage to a xenon admixture of 30 percent (x Xe " 0.3, fig. 8 a)). Recalling figure 6 a), it was observed that adding xenon to an argon background is equivalent to monotonically raising the plasma density. Therefore, a small xenon admixture to an argon discharge means creating more electrons that will mostly collide with an argon atom. As a result, the probability of ionization of an argon atom is higher than it is in a case with no or lower xenon admixture, that is, without these additional electrons. Thus, the density of Ar + ions is higher than in a discharge without xenon admixture. This effect affects the Ar + ions for low voltages as long as most neutrals are argon atoms. A raised driving voltage shifts the maximum of the Ar + ion density and the benefits of this synergy effect to higher xenon fractions.
For Xe + ions, on the other hand, this synergy effect cannot be observed ( fig. 8 b)). This observation is due to the higher ionization energy of argon. Figure 7 a 1 ) helps to understand this observation by showing the energy lost by electrons at the electrodes ε e . The general trend of the curves for ε e is a reduction by a raised driving voltage ( fig.  7 a 1 )). An equivalent conclusion is that energy is dissipated more efficiently inside the volume of the discharge. In terms of our non-equilibrium low-pressure discharge, there are just two ways for electrons to lose energy. Either they interact inelastically with the background or transfer their energy to the surface by arriving at the electrodes. The first option was discussed before ( fig. 7 a 2 ) and a 3 )), and the second option is discussed here. Both processes similarly respond to the increased driving voltage, which means that a higher driving voltage increases the ion production efficiency. Simultaneously, the energy dissemination efficiency is increased the more xenon is added to the background gas. In section 3.1, we discuss that an increased amount of xenon atoms in the discharge provides lower energetic electrons with the opportunity to get involved in inelastic processes compared to a discharge with lower or no xenon addition (see fig. 5 d)). In figure 7 a 1 ), the same trend is observed for all depicted driving voltages. As a function of the xenon fraction x Xe , the energy lost by electrons at the electrode ε e is monotonically falling. Vice versa, argon has higher thresholds for inelastic processes, especially ionization, than xenon (comp. tab. 1). Thus, adding argon to a xenon background cannot produce a higher electron density that would cause more ionization of xenon. The synergy effect does not take place for Xe + ions that benefit from additional ionization of argon. Figures 9 and 10 show, similar to figure 3, IEDFs normalized to the respective particle flux densities at the electrode surface. The difference between figure 9 and 10 is in the gas composition ( fig. 9: x Xe " 0.1, fig. 10 x Xe " 0.9). Both figures contain IEDFs for Ar + ions in the left column panels and IEDFs for Xe + ions in the right one. The difference between each figure's four rows is the altered amplitude of the RF voltage V RF given on each panel's top. Panels of the same row share the same voltage.
The figures show that for the IEDF, a raised driving voltage, first of all, means that the averaged sheath voltage xφ s y increases. This increase manifests in the width of the characteristic collisionless single bimodal peak. Its width scales with V RF and a xφ s ptqy { xsptqy [28,78]. Kawamura et al. [28] give the averaged sheath width xsptqy in terms of the collisionless Child-Langmuir law: with e the elementary charge, m i the ion mass, ε 0 the vacuum permitivity, and xj i ptqy the averaged ion current inside the sheath. For argon, the increased bimodal peak's width is found in figure 10 (left) and for xenon in figure 9 (right). From top (V RF " 100 V) to bottom (V RF " 1000 ), the width of the highest energetic bimodal peak increases (Ar + ions: fig. 10 left, Xe + ions: fig. 9 right). At the same time, a higher driving voltage at a constant pressure produces higher energetic ions. Higher kinetic energy enlarges the mean free path of those ions since the mean free path is energy-dependent [2] and the collision cross sections fall at high energies (see fig. 1). Thus, the distance between the peaks in the low energetic part of the IEDFs that are connected to charge exchange collisions is increased with the driving voltage [26,79].
Furthermore, a raised driving voltage causes the emergence of multiple bimodal structures within the IEDFs. These structures have first been reported as double peaks by Wild et al. [79]. For our scenario, they become visible for V RF " 500 V and V RF " 1000 V and establish both for Ar + ( fig. 9) and Xe + ions ( fig. 10). Here, charge exchange collisions are responsible for the appearance of the low energetic peak. Section 3.1 discusses that low energetic peaks vanish for Ar + ions when the xenon fraction x Xe is raised and vice versa. For a second or third bimodal peak to establish, two requirements have to be met. First, ions have to be able to react to the sheath electric field. Second, there has to be some sort of hybrid regime that we discussed in section 3.1. Combining these requirements also means that only charge exchange collisions that happen clearly above the averaged sheath position can establish an additional bimodal structure. Under these conditions, the slow ions produced through charge exchange experience the sheath's modulation that eventually determines their impingement energy. A charge exchange collision inside the sheath during the collapsing phase causes the ions to gain slightly lower impingement energy than a charge exchange during the expanding sheath phase.
According to Lieberman and Lichtenberg [1], there is a weak dependency between the average position of the sheath edge and the voltage amplitude (s m 9 V 1{4 RF ). Thus, it is more likely for the collisional structures of IEDFs at higher voltages to show bimodal structures. The IEDFs of Xe + ions at a xenon fraction x Xe " 0.9 ( fig. 10) are the best example for this conclusion of our study. For V RF " 100 V, the results clearly show a single bimodal peak and several non-bimodal charge exchange peaks. For V RF " 250 V, the main bimodal peak is centered around « 110 eV, and at least one additional bimodal peak around 87 eV is visible. At 500 V, the IEDF has at least four bimodal peaks (centered around « 130 eV, « 170 eV, « 180 eV, and « 210 eV). The case for V RF " 1000 V shows at least four bimodal peaks as well (centered around « 190 eV, « 260 eV, « 325 eV, and « 410 eV). For that case, solely charge exchange collisions that take place deep inside the boundary sheath and close to the electrode do not show any sign of bimodal features.
The hybrid regime of the IEDFs itself is also influenced by a raised driving voltage V RF . For a voltage amplitude V RF " 250 V, a slightly higher voltage than that of the base case, the hybrid regime appears for lower admixtures of xenon ( fig. 9) or argon, respectively ( fig. 10). Here, the broadening and amplification effects of a raised driving voltage prevail. Thus, the hybrid regime establishes earlier than for lower voltages.
The IEDFs for even higher driving voltages (see V RF " 500 V and V RF " 1000 V in fig. 9 and 10) are again more collision-dominated and show a different trend. For 500 V, the bimodal part of the distribution function is less populated than the low energetic part. For 1000 V, the highest energetic peak for both Ar + ( fig. 9) and Xe + ions ( fig.  10) are damped compared to the lowest energetic peaks. This trend arises from the fact that the cross section for charge exchange collisions for high energies drops much slower than the cross sections for isotropic scattering (comp. fig. 1 c) and d)). Therefore, charge exchange is the preferred process at high energies. For driving voltages much higher than 100 V, the hybrid regime is shifted back to higher mixing ratios.

Conclusion
The objective of this work was to investigate the ion dynamics of plasmas containing two ion species. This investigation was conducted by simulating a low-pressure capacitively coupled plasma with a mixture of argon and xenon as the background gas. The overall result is that the gas composition serves as a means to control the collisionality of the ion species and thus the ion dynamics. Section 3.1 shows that the gas composition (more specifically the argon fraction x Ar or xenon fraction x Xe , respectively) significantly affects the discharge, especially the ion dynamics. The effect on the discharge resembles a parabolic function of the plasma density and the xenon fraction x Xe (comp. fig. 2). A complete energy balance that we self-consistently calculate based on a PIC/MCC simulation helps understand this effect. Inelastic processes in xenon (e.g., ionization with ε i,Xe " 12.12 eV) have significantly lower energetic thresholds. Thus, electrons distribute their energy more efficiently when the xenon fraction x Xe is raised. We show that especially the ionization process in xenon is energetically more favorable than in argon. This disparity leads to Xe + ions being the dominant ion species for a broad range of xenon fractions x Xe .
For the ion dynamics, we present that the gas composition controls the collisional characteristics of the IEDF. Between argon and xenon, only non-resonant charge transfer collisions are possible. Three-body collisions do not occur in relevant amounts in the low-pressure regime. Therefore, a varied xenon fraction x Xe shifts the multiple low energetic peaks (characteristic for charge exchange and a collision dominated regime) from argon (most pronounced at x Xe " 0, fig. 3 left) to xenon (most pronounced at x Xe " 1, fig. 3 right). Additionally, a collisional/collisionless hybrid regime is present for specific gas fractions. Some ions experience the discharge within this hybrid regime as collision dominated while others traverse the boundary sheath without collisions. The analysis of the energy balance helps to understand these effects as well. It reveals that charge exchange is, even at low-pressures, a relevant energy loss process for ions. A raised xenon fraction x Xe depletes (Ar + ions) or contributes to (Xe + ions) this process for the respective ions ( fig. 5). Thus, the addition of xenon increases (Ar + ions) or decreases (Xe + ions) the impingement energies of the respective ions. Furthermore, the energy balance reveals optimal parameters for the impingement energy of ions in this mixture. In this context, optimal refers to overall minimal collisional losses for the ions, thus desirable conditions for processes (e.g., ion-assisted etching). For x Xe " 0.4, the combined fraction of the total energy that ions lose at the surface is maximal. This example shows that the gas compositions allow tailoring the discharge to the requirements of specific applications.
A variation of the driving voltage V RF attenuates the dominance of the Xe + ions (sec. 3.3). The reason for this observation is a synergy effect. The argon's ionization process benefits from additional electrons created during the ionization of xenon. An extensive analysis of the energy balance is needed to understand this synergy effect and differentiate why it only occurs for Ar + but not Xe + . Furthermore, it is presented that the increased driving voltage V RF intensifies structures (e.g., broadens the width of bimodal peaks) and further complicates the IEDFs (e.g., by creating multiple bimodal peaks). The energy dependence of the cross section for charge exchange causes the hybrid regime to shift to different mixing ratios when solely varying the driving voltage. Both observations are supported by the analysis of the energy balance too. Overall, the energy balance has proven to be a practical and impactful diagnostic. The results of section 3.3 show that the gas composition controls the ion dynamics over a wide range of driving voltages. However, the effect of varied gas compositions is not entirely independent of the driving voltage.
Future work based on this study will develop in two directions. On the one hand, the model system Ar/Xe needs to be left behind. The presented basic principles have to be investigated in more complex and process relevant gas mixtures like Ar/CF 4 or CF 4 /H 2 . The energy balance model can be adapted to and should be tested for these gas mixtures. On the other hand, based on this work's findings, the influence of a combination of multi-frequency discharges and a varied gas composition on the ion dynamics should be investigated. For example, a multi-frequency approach could be used to optimize the ion production, which at V RF " 100 V was found to be optimal for a xenon fraction x Xe " 0.4, further.
Another open research question is: How does the addition of secondary electron emission and realistic surface coefficients alter the ion dynamics? The argon ionization's synergy effect, especially, could be significantly affected when secondary electrons cause an amplification of the ionization process. To our best knowledge, there are no published experimental results that analyze the influence of the gas mixture on the IEDFs. Nor are there studies that experimentally report about the hybrid regime or the synergy effect within the ionization of argon. All of these studies would be crucial to validate our findings and simulation.