Controlling plasma properties under differing degrees of electronegativity using odd harmonic dual frequency excitation

The charged particle dynamics in low-pressure oxygen plasmas excited by odd harmonic dual frequency waveforms (low frequency of 13.56 MHz and high frequency of 40.68 MHz) are investigated using a one-dimensional numerical simulation in regimes of both low and high electronegativity. In the low electronegativity regime, the time and space averaged electron and negative ion densities are approximately equal and plasma sustainment is dominated by ionisation at the sheath expansion for all combinations of low and high frequency and the phase shift between them. In the high electronegativity regime, the negative ion density is a factor of 15–20 greater than the low electronegativity cases. In these cases, plasma sustainment is dominated by ionisation inside the bulk plasma and at the collapsing sheath edge when the contribution of the high frequency to the overall voltage waveform is low. As the high frequency component contribution to the waveform increases, sheath expansion ionisation begins to dominate. It is found that the control of the average voltage drop across the plasma sheath and the average ion flux to the powered electrode are similar in both regimes of electronegativity, despite the differing electron dynamics using the considered dual frequency approach. This offers potential for similar control of ion dynamics under a range of process conditions, independent of the electronegativity. This is in contrast to ion control offered by electrically asymmetric waveforms where the relationship between the ion flux and ion bombardment energy is dependent upon the electronegativity.


Introduction
Non-equilibrium capacitively coupled plasmas (CCP) at low pressure are commonly used for etching and deposition of nanoscale structures in the semiconductor industry. However, the rapid advance of the industry towards smaller feature sizes places a strong emphasis on plasma control in order to achieve the atom-scale precision required for the realisation of novel device designs. In this context, numerous plasma parameters must be optimised and controlled in order to provide favourable process conditions and produce highquality components. Two of the most important are the ion flux and ion bombardment energy at the substrate electrode. The ion flux at the substrate is one of the principal parameters in determining the rate of a process and as such it is one of the key parameters that plasma control strategies aim to optimise. The ion bombardment energy required is largely process dependent with certain processes requiring low ion energies [1,2] and other processes requiring ion energies above certain thresholds in order to break bonds on the substrate.
In single frequency CCPs, there exists a rather defined relationship between the ion flux and bombardment energy [3], as such the range of independent control over the two quantities in such plasma sources is limited. In this context, much effort has been invested to develop plasma control strategies allowing for independent control of the ion flux and bombardment energy. Many investigations into this control have focussed on driving CCPs with more than one frequency. The original concept was introduced by Löwe et al [4] and subsequently investigated by numerous authors [5][6][7][8][9][10][11]. The basic premise of the conventional dual frequency approach is to achieve independent control over the ion flux and bombardment energy by using a very high frequency to control the former and a low frequency to control the latter. However, in order for truly independent control to be achieved functional separation of the two frequencies is required in order to ensure that the effects of the two frequencies are decoupled [6]. This generally necessitates the use of a high frequency, on the order of 100MHz or more. The disadvantage of using these very high frequencies is the introduction of standing wave and skin effects which affect the radial uniformity of the plasma and, consequently, the ion flux [12][13][14][15][16][17] which can be detrimental for process outcomes.
As a result of the difficulties caused by radial non-uniformities at very high driving frequencies, numerous investigations have been carried out on multiple frequency plasma sources using lower frequencies [9][10][11][18][19][20][21][22][23][24]. In these systems, functional separation of the frequencies is generally not achieved and as a result separate control of the ion flux and ion bombardment energy is limited. Furthermore, it has been shown that in many such systems the coupling between the multiple frequencies is highly nonlinear and is mediated by complex spatio-temporal electron heating [9-11, 18, 19, 25].
In recent years, specific excitation schemes have been developed to harness the nonlinear coupling of multiple frequencies in order to gain independent control of the ion flux and ion bombardment energy, along with other plasma properties. These schemes fall under the general heading of voltage waveform tailoring. The most common implementation of this technique relies on the exploitation of the electrical asymmetry effect. For a detailed discussion of progress in this area the reader is referred to the recent review of Lafleur [26]. Briefly, it was determined in 2008 by Heil et al [27] that the use of a voltage waveform composed of a fundamental frequency and its second harmonic (for example, 13.56 and 27.12 MHz) allows for the generation of an asymmetric plasma (i.e. one with a dc self-bias) in a geometrically symmetric system. This effect is possible due to the fact that the positive and negative amplitudes of the voltage waveform are not equal, i.e. the waveform has an amplitude asymmetry. Furthermore, the dc self-bias, and therefore the ion bombardment energy at the electrode, in such a system can be controlled by varying the phase shift between the two harmonics, effectively changing the relative extent of the maximum and minimum of the voltage waveform. The ion flux is comparatively insensitive to the change in the phase shift meaning that the two quantities can be controlled independently of one another. Since the original study many investigations have demonstrated this effect through both simulation and experiment and extended it to the use of more than two frequencies to increase the amplitude asymmetry of the waveform [28][29][30][31][32][33][34][35][36][37][38][39][40][41]. It has also been demonstrated that an electrical asymmetry can be generated using sawtooth-like waveforms where the rise and fall times of the voltage waveform differ [41][42][43][44][45][46][47][48]. In this case, the dc self-bias is generated as a result of a so-called slope asymmetry as opposed to the previously discussed amplitude asymmetry.
In addition, the concept of voltage waveform tailoring also incorporates the use of voltage waveforms which do not induce an electrical asymmetry and hence maintain the symmetry of the plasma. Such waveforms have been predicted to offer enhanced control over the electron dynamics in atmospheric pressure plasma jets [49,50] and the sheath dynamics in low-pressure oxygen plasmas [51]. The fact that these symmetric tailored waveforms allow for control of the plasma properties without the need to change the symmetry of the plasma may be advantageous under certain conditions. For example, it has been shown that the dc self-bias formed in electrically asymmetric plasmas can change over a wide range dependent on the gas in which the plasma is formed [44] and can thus be difficult to predict. The electronegativity of the plasma, which in part determines the electron heating mode [52], has been found to be particularly important in this regard [31,37,44,48].
Overall, voltage waveform tailoring has proven to be a very effective method to achieve control of various parameters in plasmas of industrial interest. A topic of particular importance in such plasma sources is the role played by surface interaction probabilities in determining the overall discharge dynamics and the electron heating. Voltage waveform tailoring has previously been predicted to be a more effective method of controlling the ion flux and ion bombardment energy than conventional dual frequency approaches in cases where the secondary electron emission probability at the electrodes is high and electron heating can be dominated by secondary electrons emitted from the electrodes [53]. Furthermore, there have been a number of investigations detailing the importance of secondary electron emission coefficients in such plasma sources and the problems and opportunities presented for plasma processing applications by their variation from material to material [54][55][56][57].
It is also known that surface interaction probabilities for reactive neutral species can strongly influence the properties and charged particle dynamics of low pressure plasma sources [38,[57][58][59][60][61][62]. A prominent example is the role played by the surface quenching probability of the molecular oxygen metastable molecule O 2 (a g 1 D ) in determining the electronegativity (i.e. the ratio of the negative ion density to the electron density) of oxygen plasmas [38,48,57,59,63]. This dependence occurs as a result of the high rate coefficient for electron detachment from Oions in collisions with O 2 (a g 1 D ) molecules in comparison with the equivalent reaction with ground state molecular oxygen [64]. Thus, the overall electron detachment rate, and consequently the negative ion density, is strongly related to the O 2 (a g 1 D ) density, which is determined by its surface quenching probability. The electronegativity subsequently plays an important role in determining the electron heating dynamics in the plasma and as such it is possible for these dynamics to change significantly with variations in the surface quenching probability of O 2 (a g 1 D ). Such effects are important to consider as they can lead to process drifts in industrial applications and thus limit reproducibility from process to process. As a result, techniques to control plasma properties in plasmas which are susceptible to changing electronegativity are highly desired.
Presented here is a numerical study of the dynamics of a geometrically symmetric capacitively coupled oxygen plasma where the powered electrode is driven by tailored voltage waveforms described by the following relation: where f lf and f hf are the low and high frequencies, 13.56MHz and 40.68MHz respectively, V lf and V hf refer to the amplitudes of the voltage waveform corresponding to each of the individual frequencies and θ is the phase shift between them. Examples of typical voltage waveforms taking the described form are shown in figure 1. It should be noted that as the high frequency used in this case is an odd multiple of the low frequency these waveforms are symmetric, both in amplitude and slope, and hence do not induce a dc self-bias.
The aim of this work is to determine if the use of symmetric voltage waveforms can allow for similar levels of plasma control under both low and high electronegativity conditions, overcoming some of the difficulties associated with electrically asymmetric systems applied to electronegative plasmas. To do this, plasmas driven by the waveform described in equation (1) are characterised under the variation of (i) high and low frequency components (V hf and V lf ) and phase shift (θ) and (ii) the contribution of O 2 (a g 1 D ) to the background gas. Two values of the O 2 (a g 1 D ) density are chosen which represent the two extremes predicted in the literature through the variation of the O 2 (a g 1 D ) surface quenching probability [38,57]. The plasma electronegativity resulting for each of these O 2 (a g 1 D ) densities thus represent approximate extremes in the electronegativity of oxygen plasmas. Subsequently, the full range of control over the ion flux and the average sheath voltage at the powered electrode by variation of the waveform parameters, V hf , V lf and θ is investigated for both low and high electronegativity conditions.

Numerical simulation
In this work, a self-consistent one-dimensional (1D) fluid model, incorporating a semi-kinetic treatment of electrons and energy dependent ion mobilities is used [65]. The domain is a 1D slice of a typical geometrically symmetric CCP reactor with a discharge gap of 25mm, similar to various commercial and research systems [22,34,66]. The simulation approach and governing equations are described in detail in [65]. Briefly, the simulation solves the mass and momentum (via the drift-diffusion approximation) conservation equations for electrons and ions. In addition, the energy conservation equation is solved for electrons while their transport properties and electron impact rate coefficients are calculated in advance using the two-term approximation Boltzmann solver BOLSIG+ [67] and incorporated into the fluid simulation as lookup tables. Increases in ion temperature as a result of strong electric fields, such as in the plasma sheath, are approximated using a semi-empirical formula describing ion heating in an electric field. Furthermore, the dependence of the ion mobility on the reduced electric field is accounted for by interpolation of experimentally measured ion mobilities in the low field limit and extrapolation of these values to the high field limit assuming a rigid sphere model [65].
In comparison to particle-in-cell (PIC) simulations, which are commonly used for such investigations, the semikinetic fluid simulations used in this work have the advantage of shorter solution times, at the expense of additional assumptions. In this context, it is important to consider if these additional assumptions influence the accuracy of the simulation results. Insight into such effects can be gained from the works of Becker et al [68], Turner et al [69] and Surendra [70] who have undertaken detailed comparisons between PIC simulations and fluid simulations similar to those used in this work. In the work of Becker et al, the authors found that the investigated parameters, for example, electron density and ion flux at the electrodes, agreed to between 10% and 50% between the different fluid approaches studied and PIC simulations for He and Ar CCPs. These results were for a pressure range of 10-80Pa and are therefore consistent with the pressure of 40Pa used in this work. The agreement was found to be poorest for Ar, an issue that was attributed in part to the Ramsauer minimum in the momentum transfer cross section and its influence on electron  [71,72]. The reaction mechanism used is kept as simple as possible, incorporating eight gas phase reactions for species production and loss, as shown in table 1. The densities of neutral species are not solved for self-consistently, rather the background gas is assumed to be of fixed density, given by the ideal gas law for a pressure of 40Pa and a temperature of 300K. The assumption of a constant gas temperature of 300K, both across the simulation domain and for the different plasma conditions, is justified for the low power conditions studied in this work. In particular, the gas temperatures and electron densities measured by Wegneret al [73] for an E-mode inductively coupled plasma in O 2 , which exhibits similar plasma dynamics to CCPs, imply that the gas temperature remains at approximately 300K over the range of electron densities studied in this work i.e. 10 14 -10 15 m −3 . The neutral density in this work is composed of ground state molecular oxygen, O 2 with a variable fraction of the metastable, O 2 (a g 1 D ). This allows for the simulation of the oxygen discharge in either low electronegativity , as discussed previously. It should be noted that the aim here is not to accurately describe any particular oxygen plasma, but rather to demonstrate the different dynamics possible under two possible extremes in the electronegativity and how these may be controlled.

Charged species density profiles
In order to address the control of oxygen plasma dynamics at different electronegativities we first characterise the basic plasma phenomena under both low and high electronegativity conditions. Figure Examining the low frequency only cases, shown in figures 2(a) and (c), a pronounced change in charged particle density profiles is observed when the O 2 (a g 1 D ) content is changed. In the case where O 2 (a g 1 D ) is 16% of the ground state O 2 density the charged particle density profiles are typical of a weakly electronegative plasma where both the peak and averaged electron and negative ion densities are similar. In this case, the average electronegativity is approximately 1. As the O 2 (a g 1 D ) content is decreased to 0.5% of the ground state O 2 density, shown in figure 2(c), the electron density in the plasma is much lower than the negative ion density and hence the electronegativity rises significantly, in this case to approximately 19. Here, the negative ion density is almost entirely responsible for the quasi-neutrality of the bulk plasma.
Under dual frequency excitation when the O 2 (a g 1 D ) content is 16% of the ground state O 2 density, shown in figure 2(b), the charged particle profiles are similar to the low frequency only case, shown in figure 2(a). The main differences are that the charged species densities are higher, due to the increased efficiency of electron heating with the addition of the high frequency waveform [51] and the central plateau in each profile is wider, due to the reduction in sheath width at higher plasma densities. The positive and negative ion density profiles show a Units: rate coefficients in m 3 s −1 ; gas temperature T g in K; relative gas temperature T T 300 f(ò) indicates that the rate coefficients are obtained from the EEDF calculated with the two-term Boltzmann equation solver BOLSIG+ [67] using electron impact cross section data. Additionally, electron-impact excitation of O 2 into rotational, vibrational and electronically excited states according to the cross sections given in [74,75] is accounted for in the EEDF calculation to properly simulate electron energy losses. c It is worth noting that several values exist in the literature for this rate coefficient, however the recent measurements and discussion published by Midey et al [82] show that most of these measurements agree within the combined experimental errors.
flat central region, with slight peaks on either side. This is particularly visible in the case of the negative ions. These slight peaks in the negative ion density are a result of competing production and destruction mechanisms. Negative ions are produced in the simulation by electron attachment to ground state O 2 (reaction R4 in table 1), the rate for this reaction is maximum at the sheath edge, where electron energies are highest as the electron attachment process in oxygen has a threshold of several eV [77]. The dominant destruction of negative ions occurs by electron detachment collisions with O 2 (a g 1 D ), (reaction R8 in table 1), which is equally distributed throughout the simulation domain. The result of these two competing processes is that the net production of negative ions is slightly greater just inside the sheath edge than in the centre of the plasma. These effects cause the peaks at the edges of the density profile, and its central minimum. Overall, this has a small effect on the average electronegativity in this case which is also approximately 1.
When the O 2 (a g 1 D ) content is decreased to 0.5% of the ground state O 2 density, shown in figure 2(d), the electron density is once again significantly lower than the negative ion density meaning the average electronegativity has increased, in this case to approximately 17, similar to the low frequency only case shown in figure 2(c). Here, the peaks in negative ion density on either side of a central minimum do not occur. This is because, in the case of high electronegativity, electrons are heated throughout the plasma bulk. The more uniform spatial profile of the electron heating means that any peak in electron energy near the sheath edge is not as pronounced as in the low electronegativity case and thus approximately equal production and destruction of negative ion occurs throughout the plasma bulk. Time averaged charged species density profiles for an O 2 (a g 1 D ) content of 16% driven by (a) low frequency only excitation (V lf =200 V, V hf =0 V) and (b) dual frequency excitation (V lf =100 V, V hf =100 V) and an O 2 (a g 1 D ) content of 0.5% driven by (c) low frequency only excitation (V lf =200 V, V hf =0 V) and (d) dual frequency excitation (V lf =100 V, V hf =100 V). For all plots the phase shift θ=0°.

Time and space resolved electron dynamics
The time and space resolved electron density profiles for the same four cases as figure 2 are shown in figure 3. The powered electrode is at the bottom of each plot and the instantaneous position of the plasma sheath edge, defined according to the condition given in [83], is represented by the dashed black lines. In both cases where the O 2 (a g 1 D ) content is 16% and the electronegativity is around 1, shown in figures 3(a) and (b), the time and space resolved electron density profiles are similar to those predicted for purely electropositive plasmas [52]. In this regime of low electronegativity the electrons are concentrated in the centre of the discharge gap throughout the entire radio-frequency cycle in order to preserve quasi-neutrality while a small number of electrons approach the electrodes during periods of sheath collapse in order that the fluxes of positive and negative charges to the electrodes are equal throughout the rf cycle. In general, the same is true for both single frequency and dual frequency cases, with the main differences being a more complex temporal evolution of the sheath structure and overall higher electron density in the dual frequency case.
In the high electronegativity cases, shown in figures 3(c) and (d), the time and space resolved electron density profiles are significantly different to the low electronegativity cases. In the low frequency only case, shown in figure 3(c), the electron density peaks in the regions where the sheath is collapsed, in contrast to the weakly electronegative case under the same conditions, where the electron density peaks in the centre of the discharge gap. This occurs because, in the strongly electronegative case, the negative ion density largely compensates the positive ion density in the centre of the discharge gap, as can be observed in figure 2(c). However, because of their high inertia the negative ions cannot respond instantaneously to the rapidly varying electric fields and are largely confined to the region between the maxima in sheath edge position. These are marked with horizontal dashed white lines in figures 3(c) and (d). The electrons, which can respond to the rapidly varying fields then accumulate between the position of maximum sheath extension and the electrode in order to maintain quasi-neutrality in these regions. Similar electron density profiles have been predicted by PIC simulations for CCPs produced in highly electronegative CF 4 under comparable operating conditions [24,52].
In the dual frequency case, shown in figure 3(d), the build-up of electrons in the regions where the sheath is collapsed is less significant than in the low frequency only case. This results from the fact that the sheath maxima are closer to the electrodes in the dual frequency case, due to the higher plasma density and the sheath being fully collapsed for shorter time periods, due to the dual frequency modulation of the sheath. These factors mean that the electron density in these near-electrode regions deviates only slightly from the central electron density.
The differences in the time and space resolved electron density profiles between the low and high electronegativity cases are reflected in the time and space resolved electron impact ionisation rates for ground state O 2 , shown for the same conditions as figures 2 and 3 in figure 4. Here, the instantaneous position of the plasma sheath edge is marked by dashed white lines. The electron impact ionisation dynamics in the low electronegativity cases, shown in figures 4(a) and (b), confirm that the plasma is operated in α mode where the dominant electron heating occurs at the sheath expansion, in both the low frequency only and dual frequency cases. The structures corresponding to this sheath expansion electron heating are marked with I.
In the low frequency only case, shown in figure 4(a), one of these structures occurs per radio-frequency cycle at each electrode which is characteristic of a symmetric discharge. In the dual frequency case, shown in figure 4(b), there are two visible sheath expansion ionisation structures at each electrode as a result of the more complex temporal evolution of the plasma sheath edge. That the intensity of the first sheath expansion structure, occurring at the powered electrode at around 8ns, is greater than that of the second, occurring at the powered electrode around 33ns, is a result of the differing ion densities in the two regions of space and time. For both of these periods of sheath expansion the change in voltage per unit time is similar, however, the sheath expansion occurring around 8ns occurs close to the electrode in a region of lower ion density than the sheath expansion at 33ns. The result of this is that the velocity of the expanding sheath, which determines the electron heating as a result of the sheath expansion, is greater at 8ns than at 33ns leading to a higher electron heating rate. This gives rise to a higher rate of electron impact ionisation during the sheath expansion occurring at 8ns than that occurring at 33ns. In addition, the maximum and time averaged rates of electron impact ionisation are significantly higher in the dual frequency cases as the maximum sheath velocities are greater due to the influence of the high frequency component of the waveform. This in turn leads to higher current densities and power depositions which result in the higher densities of charged species in the dual frequency cases for both low and high electronegativies, as shown in figures 2 and 3.
In the high electronegativity cases, shown in figures 4(c) and (d), the time and space resolved electron impact ionisation dynamics are significantly different from the low electronegativity cases, shown in figures 4(a) and (b). Under low frequency only excitation, shown in figure 4(c), there is evidence of ionisation in three distinct spatial regions; at sheath expansion (I), near the collapsing sheath edge (II) and in the plasma bulk (III). Ionisation at sheath expansion is less significant compared to the low electronegativity case, and ionisation near the collapsing sheath edge and in the plasma bulk are dominant. This reflects the fact that the discharge is operating in the drift-ambipolar heating mode [52] as opposed to the α mode discussed previously. In this mode strong electric fields are present in the plasma bulk in order to accelerate the relatively small number of electrons across the plasma bulk so that current continuity is fulfilled. These 'drift' electric fields result in the ionisation observed in the bulk plasma (III). Near the collapsing sheath edge the strong ionisation features are caused by ambipolar [52] or electron pressure induced [57] electric fields which arise from the sharp electron density gradient that results from the build-up of electrons between the maximum of the sheath edge position and the electrode when the sheath has collapsed.
Under dual frequency operation at high electronegativity, shown in figure 4(d), the dominant electron impact ionisation feature at the powered electrode occurs during the initial sheath expansion from the fully collapsed sheath around 8ns, analogous to the dual frequency low electronegativity case ( figure 4(b)). However, in contrast to the low electronegativity case, there are also ionisation features visible in the bulk plasma and at the collapsing sheath edges as a result of the aforementioned drift-ambipolar electron heating. The transition, from the most intense ionisation occurring at the collapsing sheath edge to the expanding sheath edge, with the introduction of the high frequency component to the waveform occurs because of the increasing expanding sheath edge velocity. This results in increased collision-less heating of electrons at the expanding sheath edge in comparison to the low frequency only case. Furthermore, the sheath expansion ionisation feature is extended further into the plasma bulk than is observed in the low electronegativity case as a result of drift electric fields formed to preserve current continuity in the plasma bulk, as with the high electronegativity, low frequency only case, shown in figure 4(c). The change in position of the dominant ionisation rate feature with the introduction of the high frequency waveform is notable as the electron dynamics in this case resemble those of a weakly electronegative plasma, while the electronegativity is still high, and comparable to that of the low frequency only case.

Control of plasma properties
It has been demonstrated that the charged particle and ionisation dynamics can change significantly as a function of both the high and low frequency components of the driving voltage waveform and the O 2 (a g 1 D ) content of the background gas. As a result, it is important to compare the level of control over the plasma properties of industrial interest that can be achieved under both low and high O 2 (a g 1 D ) content using the presented dual frequency approach. Figure 5 shows the relative changes in the average sheath voltage (used as an estimate of the change in average ion bombardment energy) and the ion flux at the powered electrode for an O 2 (a g 1 D ) content of (a) 16% and (b) 0.5% for different driving voltage waveforms. In each plot the data points have been normalised to the respective values for the low frequency only waveform i.e. 0% high frequency component. The black squares represent a variation of the high frequency component percentage with the phase shift between the low and high frequency components of the waveform, θ=0°, going from 0% high frequency at the top left, to 100% high frequency on the top right of each plot. The red circles represent the equivalent data points for θ=180°, in this case going from 10% high frequency to 90% high frequency component. The area between both sets of points represents approximately the full range of control that can be achieved utilising both frequency and phase variations. This is because varying the phase between 180°and 360°has the opposite effect as varying the phase between 0°and 180°. This is shown in detail in [51].
Considering first the black squares, where θ=0°, the simulation predicts that the ion flux at the powered electrode changes nonlinearly as the percentage high frequency component contributing to the voltage waveform is increased. This is the case for both low (figure 5(a)) and high ( figure 5(b)) electronegativity. This is as a result of the increased efficiency of electron heating, and consequently electron impact ionisation, as the sheath velocity increases at higher percentages high frequency component. The increased ionisation rate leads to higher positive ion densities and therefore higher ion fluxes to the electrode. Additionally, in both cases, the average sheath voltage exhibits a minimum, relative to the low or high frequency only waveforms, around 50% high frequency component. This minimum is related to the shape of the dual frequency voltage waveform and the corresponding modulation of the plasma sheath motion and is determined only by the percentage of each frequency contributing to the waveform [8]. As a result, the trend in average sheath voltage with increasing high frequency component is similar at both low and high electronegativity, which is not necessarily the case when electrically asymmetric waveforms are used as the sheath voltage is also affected by the electrical asymmetry and consequently the nature of the background gas [44].
The red circles in figures 5(a) and (b), which represent dual frequency waveforms where θ=180°, exhibit lower average sheath voltages and slightly higher ion fluxes relative to the equivalent waveform where θ=0°. In these cases the lower average sheath voltages result from the lower amplitudes of the overall voltage waveforms where θ=180°, this can be seen by comparing the blue and red waveforms in figure 1. These lower voltage amplitudes result from destructive interference between the two frequencies which decreases the amplitude of the overall voltage waveform steadily as the phase shift is varied between 0°and 180°. The slight increase in the ion flux when θ=180°compared to when θ=0°results from a small increase in electron heating as a result of a higher sheath expansion velocity in the case where θ=180°, as discussed in [51].
Overall, the control of the average sheath voltage and ion flux to the powered electrode, relative to the low frequency only case, is predicted to be similar at both low and high electronegativity. In the low electronegativity case, shown in figure 5(a), the minimum sheath voltage achieved using the dual frequency waveform with θ=180°is approximately 25% lower than that of the low or high frequency only waveforms for the same voltage amplitude. This minimum has an average ion flux approximately 1.7 times greater than that for a low frequency only driving voltage. In the high electronegativity case, shown in figure 5(b), the minimum in average sheath voltage is also around 25% lower than for low or high frequency only driving voltages. The average ion flux at this point is approximately 1.2 times that for the low frequency only driving voltage. This parameter space, where the average sheath voltage is lower and the ion flux is higher than the low frequency only case may represent a favourable process window for applications requiring high ion fluxes but low ion bombardment energies at the substrate, but where radial non-uniformities may begin to become significant with the use of high frequency only waveforms [84]. The full range of control of the average ion flux is slightly reduced in the high electronegativity case compared to the low electronegativity case, with the maximum relative increases in the ion flux compared to the low frequency only case being 4.2 and 5.8 respectively.
The above explanations are valid for both high and low electronegativity cases as a result of the fact that the relationship between the average ion flux and average sheath voltage, in symmetric multiple frequency plasmas, is largely determined by the dynamics of the plasma sheath. The dominant effects of the high electronegativity are by contrast largely confined to the bulk plasma. The similar control of the two quantities under both conditions utilising this technique is partly as a result of the fact that the symmetry of the plasma is maintained using this dual frequency approach. This means that shifting the ionisation dynamics from sheath collapse dominated to sheath expansion dominated, as the percentage high frequency contribution to the driving frequency waveform is increased in the high electronegativity case, has little effect on the trend in the ion flux and sheath voltages. As a result, these trends are similar to the case where the electronegativity is low. As such, this technique may be useful for processes which are susceptible to process drifts as a result of changes in the plasma electronegativity. Furthermore, the results presented here may be useful in optimising multiple frequency plasma sources currently used in industry that do not exhibit a significant electrical asymmetry.

Conclusions
The charged species dynamics in low pressure oxygen plasmas in differing regimes of electronegativity driven by dual frequency waveforms composed of odd harmonics have been investigated using 1D numerical simulations. The fundamental phenomena underpinning the plasma structure and sustainment have been characterised through examination of the time averaged charged particle profiles and the time and space resolved electron dynamics. It has been demonstrated that oxygen discharges can operate in distinctly different modes dependent upon the O 2 (a g 1 D ) content of the background gas and the form of the driving voltage waveform. However, despite these differing modes of operation the relationship between the average sheath voltage and the average ion flux at the powered electrode is similar as a function of the high frequency component contributing to the voltage waveform and the phase shift between the two frequency components.
In general, plasmas produced in other electronegative gases, for example CF 4 , can operate at even higher electronegativity than the O 2 plasmas considered in this work. Further work remains to be done to understand if similar control can be attained in these plasmas using the odd harmonic dual frequency approach presented here. In any case, the discussed technique may be of value for industrial plasma processing applications operating in regimes where the electronegativity is susceptible to drifts as a result of changing surface conditions or other phenomena and as such provides a potentially valuable addition to the toolbox for industrial plasma processing already offered by the concept of voltage waveform tailoring.