On singlet metastable states, ion flux and ion energy in single and dual frequency capacitively coupled oxygen discharges

We apply particle-in-cell simulations with Monte Carlo collisions to study the influence of the singlet metastable states on the ion energy distribution in single and dual frequency capacitively coupled oxygen discharges. For this purpose, the one-dimensional object-oriented particle-in-cell Monte Carlo collision code oopd1 is used, in which the discharge model includes the following nine species: electrons, the neutrals O(3P) and O2(X3Σg−), the negative ions O−, the positive ions O+ and O2+, and the metastables O(1D), O2(a1Δg) and O2(b1Σg+). Earlier, we have explored the effects of adding the species O2(a1Δg) and O2(b1Σg+), and an energy-dependent secondary electron emission yield for oxygen ions and neutrals, to the discharge model. We found that including the two molecular singlet metastable states decreases the ohmic heating and the effective electron temperature in the bulk region (the electronegative core). Here we explore how these metastable states influence dual frequency discharges consisting of a fundamental frequency and the lowest even harmonics. Including or excluding the detachment reactions of the metastables O2(a1Δg) and O2(b1Σg+) can shift the peak electron temperature from the grounded to the powered electrode or vice versa, depending on the phase difference of the two applied frequencies. These metastable states can furthermore significantly influence the peak of the ion energy distribution for O2+-ions bombarding the powered electrode, and hence the average ion energy upon bombardment of the electrode, and lower the ion flux.


Introduction
Capacitively coupled plasma (CCP) oxygen discharges have various applications in plasma processing, such as ashing of photoresist, etching of polymer films, oxidation and deposition of thin film oxides. In these applications the interaction of energetic ions with solid surfaces is a key process. The positive ions, mainly created in the plasma bulk, are accelerated through the space-charge sheath adjacent to the electrodes. Etch profiles are defined by the ion energy distribution (IED) and the ion angular distribution (IAD) on the substrate. Therefore, it is important to understand how ions acquire energy as they travel through the sheath in order to be able to predict and to control the IED and IAD. The IED is an important property of the plasma. It contains information about the bombarding energy of ions striking the substrate or the discharge walls. In addition, ) and O 2 (b Σ + g 1 ), and an energy-dependent secondary electron emission yield for oxygen ions and neutrals, to the discharge model. We found that including the two molecular singlet metastable states decreases the ohmic heating and the effective electron temperature in the bulk region (the electronegative core). Here we explore how these metastable states influence dual frequency discharges consisting of a fundamental frequency and the lowest even harmonics. Including or excluding the detachment reactions of the metastables O ( ∆ a g 2 1 ) and O 2 (b Σ + g 1 ) can shift the peak electron temperature from the grounded to the powered electrode or vice versa, depending on the phase difference of the two applied frequencies. These metastable states can furthermore significantly influence the peak of the ion energy distribution for O + 2 -ions bombarding the powered electrode, and hence the average ion energy upon bombardment of the electrode, and lower the ion flux.
Keywords: capactively coupled discharge, oxygen, particle-in-cell/Monte Carlo collisions, detachment (Some figures may appear in colour only in the online journal) Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. one can determine the position where the ions are formed in the plasma and the collisionality of the discharge from the IED. The oxygen chemistry is rather involved, and the discharge includes both atomic and molecular ions, atoms and molecules and various metastable states. The two low lying singlet metastable states of the molecular oxygen, O ( ∆ a g 2 1 ) and O ( ) Σ + b g 2 1 , can be present in the discharge in significant amounts since the ∆ a g 1 state is relatively stable against deactivation by collisions with other molecules and chamber walls, while the Σ + b g 1 state is effectively produced by energy transfer from the metastable atom O( 1 D) [1]. The electron kinetics and the electron heating properties are essential for the charged particle reproduction that balance the particle losses in steady state discharge operation. In radio frequency (rf) driven discharges such as the CCP the oscillating electric field accelerates the electrons to energies that are high enough to ionize enough of the gas atoms or molecules to maintain the discharge. In these discharges there is power deposited by ohmic heating in the bulk plasma due to collisional momentum transfer between the oscillating electrons and the neutrals. For this mechanism to be dominant the electron neutral mean free path must be smaller than or comparable to the discharge dimensions. However, at low pres sure the main contribution to the electron heating is often due to the rapid movement of the electrode sheaths. The electron heating is via stochastic or collisionless electron heating by the expanding sheaths that leads to a generation of electrons beams that then enter the plasma bulk [2]. This sheath motion and thus the stochastic heating can be enhanced by self-excited non-linear plasma series resonance (PSR) oscillations in asymmetric discharges [3][4][5]. When a discharge is operated through this electron heating mechanism it is commonly referred to as the α-mode [6]. At high applied voltages and pressures secondary electron emission can contribute or even dominate the ionization. This operation regime is referred to as γ-mode [6]. In electronegative discharges large electron density gradients can develop within the rf period and lead to generation of ambipolar fields that can accelerate the electrons and is referred to as the drift-ambipolar (DA) mode [7]. A good summary of the fundamental mechanisms of the electron heating covering both simulations and experimental results is given by Donkó et al [8]. Earlier we have demonstrated how these singlet metastable molecular states influence the electron kinetics and the electron heating mechanism in the capacitively coupled oxygen discharge operated at a single frequency of 13.56 MHz [9][10][11]. We found that at higher pres sure (50-500 mTorr) the electron heating occurs mainly in the sheath region and at low pressure (10 mTorr), ohmic heating in the bulk plasma (the electronegative core) dominates. We found that the detachment by the singlet molecular metastable states is the process that has the most influence on the electron heating process in the higher pressure regime, while it has only a small influence at lower pressure [10,11].
The applied power, the discharge pressure, the driving frequency and the electrode separation all play a crucial role in determining the IED. When using a single frequency voltage source, sheath voltages are high even though the ion flux is relatively low, and an increase in ion flux will result in a further increase in sheath voltage. This leads to high ion bombarding energies, which in turn can cause damage to wafers placed on the driven electrode, lead to a loss of linewidth control or limit reproducibility in batch processing [12]. The ion bombardment energy and the ion flux on the substrate cannot be controlled independently in conventional single frequency CCP discharges. A separate control of ion properties, the ion flux and the mean ion energy as well as the shape of the IED at the electrodes, is of significant importance in plasma processing applications. Independent control of the ion flux and ion energy was suggested in dual frequency CCP discharges over two decades ago [13,14]. The application of two or even multiple frequencies provides an enhanced separate control of ion flux and mean energy as compared with single-frequency discharges and has been explored extensively by both simulations [8,[15][16][17][18][19] and experiments [8,13,14,[20][21][22][23] over the past two decades. In the case of two frequencies or dual-frequency (DF) excitation two approaches have been applied [8]: (i) The 'classical' DF discharges where the excitation is a combination of a high-frequency f H and a lowfrequency f L component. A separate control is effective when f f L H and f L is typically in the MHz range while f H is in the range of few tens to hundreds of MHz. Then the high-frequency voltage amplitude controls the plasma density, and thus, the ion flux to the electrodes and the low-frequency voltage controls the mean ion energy. This is due to the fact that the electron heating rate is proportional to the square of the driving frequency. However, interference between the high and low frequencies can cause loss of separate control of the ion energy and flux [16,24]. The other approach is the (ii) the electrical asymmetry effect (EAE) where the discharge is driven by a superposition of a fundamental and a second harmonic [18,23,25]. Heil et al [25] discovered that an asymmetry can be achieved in a CCP discharge by applying a fundamental frequency along with its even harmonics. They showed theoretically and by simulations that a dc-self bias develops in the discharge, which is dependent on the phase difference between the two applied frequencies. Thus, the dc-self bias, and hence also the energy of ions bombarding the electrodes, can be controlled by varying the phase angle. However, the ion flux is for the most part dependent on the applied voltage so a separate control of the ion energy and ion flux can be achieved using this method. It has been verified both experimentally and numerically, both for electropositive argon discharges and electronegative oxygen discharges [5,18,23,[26][27][28][29]. This method has proven more effective than applying a high frequency and low frequency component to have a separate control of the ion energy and the iox flux, as proposed by Goto et al [13]. More recently it has been explored to tailor the voltage waveforms by adding a number of harmonics and vary the shape of the waveform by individually vary the amplitude and phase of the harmonics to create arbitrary waveforms [30]. When these multi-harmonic waveforms are applied to the electrodes the plasma exhibits an asymmetric response and a dc-self bias can be generated even for perfectly geometrically symmetric systems. The electrically asymmetric plasma response arising from the application of such waveforms can be split up into [31,32]: the amplitude asymmetry effect (AAE) and the slope asymmetry effect (SAE). The AAE is due to the difference between the maximum and minimum values of the driving voltage waveform, while the SAE (sawtooth-like waveforms) arises from the different positive and negative slopes of the applied voltage waveform.
Recently there have been attempts to control the IED in oxygen CCP discharges. Derzsi et al [33] explored two approaches: 'Peaks' and 'valleys' waveforms, and electrical assymmetry effect (EAE) waveforms in order to control the IED in oxygen CCP both experimentally and by PIC/MCC simulations. However, the simulation results indicate that the discharge operates at high electronegativity. They also observe the key role the singlet metastable O 2 (a ∆ g 1 ) plays in the discharge. In a subsequent study they report a transition of the electron heating mode from the DA-mode to the α-mode as the number of harmonics that compose the voltage waveform are changed or the gas pressure [32]. Furthermore, they find that the number of harmonics has a strong influence on the discharge electronegativity and the generation of dc-bias. Schüngel et al [34]  ) was taken to be a constant fraction of the O 2 density, and the metastable O ( ) was not included in the discharge model. They found that the density profiles of the charged particles shifted to the right and left depending on the phase angle θ, but since the time-averaged power absorption of electrons was approximately constant, the peak densities and electronegativity varied only within ±10%. The dc self bias was found to be approximately linear with the phase angle θ for θ < < 0 90 and θ < < 90 180 . Furthermore, the PIC/MCC simulations agreed well with the experiments, apart from the fact that the experimental setup was inherently asymmetric so the dc-self bias profile was shifted compared to the symmetric case in the PIC/MCC simulations. Zhang et al [29] investigated dual frequency argon and oxygen discharges operated at 30 mTorr and 103 mTorr with PIC/MCC simulations. They found that the self-bias increased approximately linearly with the phase angle θ, for θ < < 0 90 for both discharge species, hence providing control of the ion energy. However, the ion flux significantly varied with θ for the oxygen discharge, and therefore limiting control of the plasma density, power absorption and electronegativity.
Here we use the 1d-3v particle-in-cell Monte Carlo collision code oopd1 to explore the ion energy distribution function (IED) for discharges including and excluding the detachment processes by the metastables O ( ∆ a g ) and . We both consider a discharge operated at a single frequency of f = 13.56 MHz, and a dual frequency discharge operated at a fundamental frequency of f = 13.56 MHz and its second harmonic, 2f = 27.12 MHz. The model setup, the plasma chemistry and the parameters used in the calculations are given in section 2. In section 3 we explore the variations in the IED while varying the electrode spacing and pressure, including and excluding the detachment process. We also discuss the influence of the singlet metastable states on the electron heating mechanism for the dual frequency discharge. Concluding remarks are given in section 4.

The simulation
We assume a capacitively coupled discharge operated between two electrodes, where one is grounded and the other is driven by an rf voltage in the single frequency case and in the dual frequency case. We study different values of the phase angle θ, from 0° to 90° in steps of 15°. The time step is taken to be ∆ = × − t 3.68 10 11 s and the simulation grid consists of 1000 uniformly distributed cells. The grid spacing and timestep are chosen such that the electron plasma frequency and the electron Debye length of the low-energy electrons are resolved where ω pe is the electron plasma frequency. The simulation was run for 2750 rf cycles or ).
The full oxygen reaction set and the cross sections used have been discussed in our earlier works and will not be repeated here [9,11,35]. The treatment of energy and scattering angles in the elementary collision processes such as electron-neutral, electron-ion, ion-neutral collisions, and neutral-neutral collisions is as it was implemented in the xpdp1 code and discussed by Vahedi and Surendra [36]. The only difference is that oopd1 and xpdp1 use different algorithms for the scattering of the incident and ejected electrons as oopd1 uses a relativistic treatment of the electrons and a revised differ ential cross sections for electron scattering as discussed elsewhere [35]. For this current work we added a feature that calculates the dc-self bias iteratively ensuring that the charged particle fluxes averaged over one period of the fundamental driving frequency are equal on the both the electrodes [37]. In oopd1 it is possible to implement different particle weights for each species, where particle weight is the ratio of the number of real particles to computational particles. The particle weights used in this study are stated in table 1. We furthermore used a subcycling factor of 16 for the heavy particles, and their initial density profiles were taken to be parabolic [38]. The neutral particles are only tracked kin etically if their energy reaches a preset threshold value, but the charged particles are tracked at all energies. The thresholds are chosen so that the number of simulated particles remains within a suitable range typically 10 4 -10 5 particles. Particles with energy less than the threshold energy are assumed to have a fixed density and a Maxwellian velocity distribution at the neutral gas temperature = T 26 n mV. The partial pres sures for the neutral species were determined with a volume averaged (global) model at 50 mTorr [39], and these values are used here for all pressures studied. These are the same partial pressures as used in our earlier works [9][10][11] and are kept the same regardless of pressure to ease comparison, however, it should be noted that these values are expected to vary with pressure. The preset threshold energies for each particle, wall quenching and recombination coefficients, and partial pressures as used in the model, are listed in table 1. The wall recombination coefficient of 0.5 for the neutral atoms O( 3 P) was measured by Booth and Sadeghi [40] for a pure oxygen discharge in a stainless steel reactor at 2 mTorr. Note that this is a rough assumption as it is known that the wall recombination coefficient drops significantly with increased pressure [41]. This leads to underestimation of the atomic oxygen density. However, the atomic oxygen density is low and is expected to decrease with increased pressure so this does not have significant influence on the results here. We use the same value for the recombination of O( 1 D). We use a quenching probability of 0.007 for the metastable ) as estimated by Sharpless and Slanger [42] for iron. This might be an overestimate, since their estimate for aluminum is <10 −3 , and Ryskin and Shub [43] measured a quenching coefficient of 0.0044 for iron and × − 5 10 5 for aluminum. Also Derzsi et al [32,33] find that a value of 0.006 results in the best overall agreement between experimental and simulation data for the ion fluxes at the electrodes in an oxygen discharge. They also point out that the value of this quenching coefficient is not well known and explore its effect on the excitation rates of O(3p 3 P) state and compare to experimental findings by phase resolved optical emission spectroscopy (PROES) to confirm their choice [32]. A higher quenching coefficient ) density, and thus less effective detachment by the O ( ∆ a g ) and higher electronegativity [44]. The quenching coefficient for the metastable molecule ) is estimated to be 0.1, assuming that this quenching coefficient should be around two orders of magnitude larger than for the O ( ∆ a g ) state [45]. However, this may be an overestimate as a value of 0.026 has been measured for nickel and a value of 0.01 for copper [46]. Adding secondary electron emission yield leads to increased electron density and decreases the sheath width [11]. Secondary electron emission is included in our discharge model, and we use the fit for secondary electron emission yield as a function of energy for oxygen bombarding oxidized metal surfaces developed elsewhere [11]. ) and O ( )

Results
. We see that detachment by the singlet metastable has more significant influence on the IED than detachment by the singlet metastable O ( ∆ a g ). This agrees with our earlier observations that the O ( ) state has more significant influence on the electron kinetics in the oxygen discharge [48]. At low pressure (10 mTorr), the mean free path of the O + 2 -ions is much longer than the sheath width, so most ions that get accelerated across the sheath region arrive at the electrodes with the mean sheath voltage drop or the time-averaged plasma potential of V 0 /2 = 111 V. This is due to the fact that under these conditions, most ions are formed in the bulk region through electron impact ioniz ation and then travel from the plasma bulk through the sheath to the electrode without experiencing collisions. We see that adding detachment by the O ( ) to the discharge model shifts the ion energy peak to slightly higher energy in a dual frequency discharge as seen in figure 1(d). As the pressure is increased the IED broadens and the mean kinetic energy shifts toward lower values. Furthermore, as the pressure is increased, secondary peaks appear in the IED due to the formation of lowenergy ions, created within the time-varying rf sheath, mainly due to charge-exchange collisions. The oscillating peaks of the IED are generated by the sheath dynamics due to ionization and charge exchange collisions as demonstrated by Babaeva et al [49]. When a slow ion is created at a certain location, the electron density front may be between this newly created ion and the electrode. Thus the ion is within the quasi-neutral region and experiences no electric field. As the electron density front moves away from the electrode (the sheath expands), the ion accelerates towards the electrode. The ions are bunched and each rf cycle creates a new bunch. These bunches appear as distinct peaks in the IED. Note that the locations of the peaks do not depend on the nature of the ion creation. A recent work of Schüngel et al [50] introduces a simple model that takes into account the ion creation in the sheath and demonstrate how it gives the IED, in a CCP driven by arbitrary voltage waveforms. They use the model to demonstrate how the primary ions lead to relatively narrow distribution around the maximum ion energy, while the secondary ions arrive at the electrode with energies between zero and the energy of the primary ions. Thus at intermediate pressure (50 mTorr), the ions participate in charge exchange collisions with the neutrals in the sheath region, as seen in figures 1(b) and (e) for single frequency and dual Table 1. The parameters of the simulation, the particle weight, the threshold above which dynamics of the neutral particles are followed, the wall recombination and quenching coefficients, and the partial pressures used.

Species
-frequency respectively. The IED profile is a superposition of the distribution formed due to the charge exchange processes and a distribution with a peak at approximately 111 V, which turns out to be a flat IED profile, as seen in figure 1(b). For both single and dual frequency discharges the peak ion energy is shifted to lower energy with the addition of the detachment processes. As the pressure is increased further, the IED becomes more heavily weighted toward lower energies. At the highest pressure value explored (200 mTorr), the mean free path of the ions is much shorter than the sheath width, so the ions will typically undergo charge exchange col lisions in the sheath before arriving at the electrode as seen in figures 1(c) and (f ) for single frequency and dual frequency discharges, respectively. Thus, the IED at 200 mTorr, shown in figure 1(c), resembles a roughly exponential decreasing ion flux with energy, with distinct peaks due to charge exchange collisions. The IEDs at 50 mTorr, as seen in figure 1(e), agree with the exper imentally determined IED for the valley case reported by Derzsi et al [33]. In particular the high energy peak is very sharp in the experimental results as is the case when the detachment process are included in the model (full reaction set). Similarly, the IEDs at 200 mTorr, as seen in figure 1(f), agree with the experimentally determined IED for the valley case reported by Derzsi et al [33]. The change in the IED with pressure also agrees qualitatively with the measurements of Wild and Koidl [51]. They measured IED in an asymmetric oxygen plasma for various pres sure values, and at a bias voltage of 500 V. They found that at low pressures, the IED is saddle shaped with a double peak centered around the bias voltage, due to the ions that accelerate in the sheath without engaging in collisions. At higher pressures they found the saddle shape to disappear, which was explained by the creation of thermal ions in the sheath caused by charge exchange collisions and rf modulation of the sheath potential. Kawamura et al [52] analyzed IED's in high and low frequency rf discharges theoretically. The two regimes are determined by the ratio of the rf period (τ rf ) and the time it takes an ion to traverse the sheath at the dc voltage V 0 (τ ion ), i.e. the crucial parameter is / τ τ rf ion . They find that where n 0 is the bulk plasma density. In the low frequency regime ( / τ τ 1 rf ion ) the ions respond to instantaneous sheath voltage, since they traverse the sheath in only a fraction of an rf cycle. The IED is broad and bimodal, and its width ∆E i approaches the maximum sheath voltage drop. However, in the high frequency regime ( / τ τ 1 rf ion ), the width of the bimodal IED curve is proportional to / τ τ rf ion , so it approaches zero with decreased discharge pressure. In this limit, the ions take much longer than one rf period to traverse the sheath so the phase of the rf cycle at which they enter the sheath is not important for the IED shape, and the ions only respond to average sheath voltages. Hence the IED curve has a narrow peak around the average sheath voltage, which is the same curve as we see in our simulations for the lowest pressure value of 10 mTorr.   ) and O ( ) is excluded. The potential drop at each sheath should be around V 0 /2 = 75 V, but as the phase angle is adjusted the maximum energy of the O + 2 -ions striking the powered electrode can be adjusted from around 40-80 eV. When the full reaction set is included in the discharge model, the peak of the IED curve is   ) and O ( ) is excluded from the model. Figure 3 shows the IAD at the powered electrode for the same case. We see that as the phase angle is increased from 0° to 90° the IAD is shifted to the right, causing a higher average impact angle at the powered electrode. Figures 6(a) and (b)  ) and O ( ) is included in the simulation, the average energy is lower compared to when it is excluded as seen in figure 6(a). Furthermore, the average impact angle is slightly larger when the detachment reactions are included in the discharge model as seen in figure 6(b). The x-component of the ion velocity is the only component that gets accelerated across the sheath region, while the other two components remain unchanged in a collisionless discharge and approximately equal to the thermal velocity in the bulk region. Thus, the ion velocity is anisotropic at the electrode, and the higher the energy of the ion bombarding the electrode, the narrower is the impact angle. In accordance with this, the average impact angle at the powered electrode gets significantly larger as the average energy of the bombarding ion gets smaller, when the phase angle θ is varied from 0° to 90°, as can be seen in table 2. However, for higher pressure discharges where the ions undergo charge exchange and elastic collisions with the thermal neutral background, the IED is shifted towards lower energies and hence the IAD broadens [52]. We see that the average flux decreases significantly when the detachment processes are included in the reaction set from × 11-12 10 17 m −2 s −1 to roughly × 8 10 17 m −2 s −1 , as seen in figure 6(c). When the detachment reactions are included in the discharge model, the average ion flux varies by around ±5%, but only by around ±3% when they are excluded. These values are much lower than the ones found by Zhang et al [29], who found the ion flux to vary with θ by ±12% in a 30 mTorr dual frequency discharge, and by ±15% in a 103 mTorr dual frequency discharge with a gap separation of 2.5 cm. They found dramatic changes in the discharge properties with θ and concluded that a separate control of ion flux and energy in oxygen discharges was more limited than in argon discharges, where the flux varied only by ±5% at 30 mTorr, and ±12% at 103 mTorr. Note that here the background partial pressures of the various neutral species is kept the same for all pressures so this effect is only due to detachment processes. In reality there are variations in the partial pressures with pressure that could also have similar effect. Also in this context it is important to note that the treatment of energy and scattering angles in the elementary collision processes such as ion-neutral col lisions, and neutral-neutral collisions of heavy particles is as discussed by Vahedi and Surendra [36].
In ) and O ( ) is excluded from the discharge model. We see that as the detachment reactions are added, the electron heating goes from being mostly ohmic bulk heating to col lisionless heating in the sheaths. Furthermore, the sheath width increases when the detachment reactions are added to the discharge model. This is in accordance with what we have demonstrated earlier for the single frequency oxygen discharge [9]. As the phase angle is increased from 0° to 90° the sheath width at the powered electrode decreases, while the sheath width at the grounded electrode increases. Figure 5 shows the IED when both the phase angle θ is varied from 0° to 90°, and the discharge gap L is varied from 2.5 cm to 6.5 cm.
In figure 5(a) we explore the case when the full reaction set is used in the discharge model, but in figure 5(b) we exclude the detachment reactions by the metastables O ( ∆ a g ) and O ( ) . We see that for a small discharge gap (L = 2.5 cm) the exclusion of the detachment reactions causes the IED to have a clearer peak, as we also observed in figure 2. However, for a larger discharge gap (L = 6.5 cm), the peak of the IED is much more apparent in the case where the full reaction set is used in the discharge model. The discharge gap does not affect the position of the IED peak for different values of θ; as previously observed the IED peak position varies from 40 eV to 75 eV depending on the phase angle. For θ = 90 , where the maximum peak of the IED is at its lowest, we see a secondary peak in the IED curve at around 65 eV. This effect is not observed for the other phase angles. Figure 7 shows a spatiotemporal plot for one rf period of the electron heating rate profile for the same case, for the values θ = 0 , θ = 45 and θ = 90 . The left column shows the case with the full reaction set but the right column shows the case when detachment by both O ( ∆ a g ) and O ( ) is excluded from the discharge model. In all cases the electron ) and O ( ) is excluded from the discharge model, the electron heating occurs both in the sheath and bulk region, but there is also strong electron cooling in the sheath region during the sheath collapse, which leads to a net very little net heating in the sheath region averaged over one period, compared to the bulk region. In the case where the full reaction set is included in the discharge model, the heating occurs mainly in the sheath region, which is wider than in the case when the detachment reactions are excluded. This transition coincides with a decrease in the electronegativity as discussed elsewhere [10,11]. So by adding the detachment by the metastables decreases the electronegativity to <1, and the discharge operating mode changes from the α-DA-mode to α-mode. Similarly, Derzsi et al [32] see a strong decrease in the electronegativity as harmonics are added to the voltage waveforms and an operation mode transition from α-DA-mode to α-mode is observed for 10 and 15 MHz driving frequency. We also see that when the sheath heating occurs mainly close to the powered electrode, for θ = 0 , the ohmic heating occurs primarily close to the grounded electrode. Conversely, when the heating occurs primarily close to the grounded electrode, for θ = 90 , the ohmic heating occurs closer to the powered electrode. In accordance with this, the effective electron temper ature is higher close to the powered electrode for θ = 0 and higher close to the grounded electrode for θ = 90 when the full reaction set is included, and conversely for the case when detachment by both O ( ∆ a g 2 1 ) and O ( ) is excluded from the discharge model. This can be seen in figure 8, which shows a spatiotemporal plot for one rf period of the effective electron temperature for this case at the phase angles θ = 0 , θ = 45 and θ = 90 . The left column corresponds to the case where the full reaction set is included in the discharge model, but the right column corresponds to the case where the detachment reactions are excluded. We note that the effective electron temperature is significantly higher when the detachment by the singlet metastable states is neglected. The time averaged effective electron temperature in the discharge center is = T 3.2 eff eV when the full reaction set is used in accordance with experimental findings in a single frequency 13.56 MHz and 75 mTorr discharge, operated at 150 V [53,54]. When the detachment processes are neglected the time averaged effective electron temperature is much higher or = T 4.4 eff eV.

Conclusion
The one-dimensional object-oriented particle-in-cell Monte Carlo collision code oopd1 was used to explore the influence of the metastable states on the IED and IAD in single and dual frequency oxygen discharges. We find that the detachment reactions by the metastable states O ( ∆ a g 2 1 ) and O ( ) have a significant influence on the discharge properties. In particular, the average energy of the O + 2 -ions bombarding the powered electrode, can be significantly affected depending on whether the detachment reactions are included in the discharge model. Furthermore, the electron heating changes from being mainly ohmic heating in the bulk to become sheath heating with the inclusion of these reactions. In dual frequency discharges this causes the peak electron temperature to shift from the powered electrode to the grounded electrode or vice versa, depending on the phase angle θ. The average ion energy is somewhat lower and the average ion flux is significantly lower when the full reaction set is included in the simulation.