Electron power absorption in capacitively coupled neon-oxygen plasmas: a comparison of experimental and computational results

Phase Resolved Optical Emission Spectroscopy (PROES) measurements combined with 1d3v Particle-in-Cell/Monte Carlo Collisions (PIC/MCC) simulations are used to study the electron power absorption and excitation/ionization dynamics in capacitively coupled plasmas (CCPs) in mixtures of neon and oxygen gases. The study is performed for a geometrically symmetric CCP reactor with a gap length of 2.5 cm at a driving frequency of 10~MHz and a peak-to-peak voltage of 350 V. The pressure of the gas mixture is varied between 15 Pa and 500 Pa, while the neon/oxygen concentration is tuned between 10% and 90%. For all discharge conditions, the spatio-temporal distribution of the electron-impact excitation rate from the Ne ground state into the Ne $\rm{2p^53p_0}$ state measured by PROES and obtained from PIC/MCC simulations show good qualitative agreement. Based on the emission/excitation patterns, multiple operation regimes are identified. Localized bright emission features at the bulk boundaries, caused by local maxima in the electronegativity are found at high pressures and high O$_2$ concentrations. The relative contributions of the ambipolar and the Ohmic electron power absorption are found to vary strongly with the discharge parameters: the Ohmic power absorption is enhanced by both the high collisionality at high pressures and the high electronegativity at low pressures. In the wide parameter regime covered in this study, the PROES measurements are found to accurately represent the ionization dynamics, i.e., the discharge operation mode. This work represents also a successful experimental validation of the discharge model developed for neon-oxygen CCPs.


Introduction
Capacitively Coupled Plasmas (CCPs) represent important reactor types widely used in plasma etching, sputtering, deposition, and cleaning processes providing the basis for semiconductor manufacturing and biomedical applications [1][2][3][4][5]. In such plasmas, efficient control of the particle properties at the electrodes (e.g., the fluxes of chemically reactive species and the ion flux-energy distribution), which determine the plasmasurface interaction, is a basic requirement. In order to achieve improved process control in such systems, a detailed understanding of the electron power absorption and ionization dynamics is necessary. The spatio-temporal distributions of the electron power absorption and the ionization within the radiofrequency (RF) period determine the operation mode of the discharge. In low-pressure CCPs, several different discharge operation modes can be identified: the α-mode and the γ-mode [6] are typical in electropositive gases, while the driftambipolar (DA) mode [7] and the striation (STR) mode [8][9][10][11][12] can be observed in electronegative gases. In the α-mode, the ionization, caused by energetic electrons accelerated by electric fields during the times of sheath expansion, is concentrated at the bulk side of the expanding sheath edge, while in the γ-mode, the ionization, dominated by secondary electrons (SEs) accelerated by the strong electric field inside the sheaths, is concentrated within the sheath region. In the DA-mode, the ionization can be observed across the whole bulk region and at the collapsing sheath edges, as it is generated by fast electrons accelerated by the strong drift electric field in the plasma bulk caused by the low conductivity of the plasma bulk, and by the ambipolar electric fields at the sheath edges caused by the strong gradients of the electron density. In the STR-mode, which develops when both positive and negative ions can react to the fast variation of the RF electric field, the ionization, concentrated within the bulk region, exhibits layered structures called "striations"; these are generated as a result of the formation of alternating space charges and the modulation of the electric field and the energy gain of electrons in the plasma bulk. The simultaneous presence of the ionization patterns characteristic of the different discharge operation modes can be observed in low-pressure CCPs under specific discharge conditions, as well as transitions between these modes by varying the external control parameters [13][14][15][16][17][18]. Experimental observation of the space and time-resolved plasma emission based on Phase Resolved Optical Emission Spectroscopy (PROES) [19][20][21] provides invaluable information on the electron power absorption and excitation/ionization dynamics in low-pressure CCPs [7,8]. For PROES measurements, Ne is often used as a tracer gas due to its favorable spectroscopic properties. By adding Ne in a small concentration (typically 5% -15%) to the background gas, and measuring the emission from a carefully selected atomic transition (e.g., Ne 2p 5 ( 2 P o 1/2 )3p 0 → 2p 5 ( 2 P o 1/2 )3s 1 with a wavelength of 585.25 nm) the spatio-temporal electron impact excitation rate from the ground state into the upper level can be derived. This way, by selecting an emission line resulting from an excited state with a high electron impact excitation threshold energy, information about the dynamics of high-energy electrons (which are typically responsible both for the excitation and the ionization processes in the discharge) can be obtained, making PROES an effective non-intrusive diagnostic technique. Although PROES provides information about the spatio-temporal distribution of the electron-impact excitation dynamics from the ground state into the selected excited atomic state in the discharge, it is generally considered to probe the discharge operation mode (which, in turn, is determined by the spatio-temporal distribution of the ionization dynamics) as well. Recently, the applicability of PROES to probe the discharge operation mode was investigated in low-pressure CCPs operated in pure neon [22]. In this work, a detailed comparison of computational and experimental results focusing on the spatio-temporal distributions of the excitation and the ionization rates in geometrically symmetric singlefrequency neon CCPs has been performed in a wide parameter regime. Particle-in-Cell/Monte Carlo Collisions (PIC/MCC) simulations [23][24][25][26][27][28][29][30][31][32][33][34] and PROES measurements were carried out at driving frequencies ranging from 3.39 MHz to 13.56 MHz and pressures from 60 Pa to 500 Pa, at a fixed peak-to-peak voltage of 330 V and an electrode gap of 2.5 cm. At fixed frequencies, transitions of the discharge operation mode from the α-mode to the γ-mode were observed by increasing the pressure. The electron impact excitation rates from the ground state into the Ne 2p 5 ( 2 P o 1/2 )3p 0 state (for which we will use the simplified notation Ne 3p 0 ) obtained from PROES measurements and PIC/MCC simulations were in good agreement in all cases. However, it was found that PROES results do not always probe the ionization and significant γ-mode ionization can take place in the discharge even in the cases when this is not seen in the spatio-temporal distribution of the Ne 3p 0 excitation. This was explained by the difference in the cross sections of the electron impact excitation into the observed level and the ionization as a function of the electron energy [22]. Although the threshold energy of the Ne 3p 0 excitation process is close to the one of the ionization, the shapes of the cross sections are different within the electron energy range of the discharges studied. As a consequence of this, it is more likely that the highly energetic SEs within the sheaths induce ionization rather than excitation. Due to their high relevance in material processing, oxygen CCPs have been studied extensively, both experimentally and by simulations [13,[35][36][37][38][39][40][41][42][43][44][45][46][47][48][49][50]. In pure oxygen CCPs, transitions between the α-mode and the DA-mode have been found by changing the gas pressure [13,14], the gap distance [14,17], the driving frequency [13,47], the driving voltage waveform [13,15,18,46,51,52], and the external magnetic field [49]. A transition from a hybrid DA-α mode to a pure α-mode was observed by increasing the driving frequency at a constant pressure [47], as well as by increasing the pressure or the electrode gap [14]. PIC/MCC simulations and PROES measurements showed a transition from the DA-mode to the α-mode in oxygen CCPs driven by tailored voltage waveforms by changing the number of consecutive harmonics included in the driving voltage waveform at high base frequencies [13]. In this work, we study the electron power absorption and excitation/ionization dynamics in capacitively coupled plasmas (CCPs) operated in mixtures of neon and oxygen. PROES measurements are combined with PIC/MCC simulations in a wide parameter regime. The study is performed in a geometrically symmetric CCP reactor with a gap length of 2.5 cm, at a driving frequency of 10 MHz and a peak-to-peak voltage of 350 V. The pressure of the gas mixture is varied between 15 Pa and 500 Pa, while the mixing ratio of neon-to-oxygen is between 10% and 90%. Transitions between different discharge operation modes are found to be induced by changing the gas pressure and varying the neon/oxygen concentration ratio in the discharge, as well as the development of localized bright emission/excitation features at the edges of the bulk region at high pressures and high O 2 concentrations. The good qualitative agreement between the PROES results and the PIC/MCC simulation results supports the validity of the discharge model developed for neon-oxygen CCPs. While experimental validation of simulations is of major importance for the reliability of the simulation results, this has only been done for some inert gases [22,53] and electronegative gases [48], but such efforts are barely found for reactive gas mixtures [15], which are required for process development. The paper is structured in the following way. In section 2, the experimental setup is presented. In section 3, the simulation method and the neon-oxygen discharge model is described, including details on the collision processes in section 3.1, the surface processes in section 3.2, the gas heating module in section 3.3, the calculation of the density of O 2 (a 1 ∆ g ) metastable molecules in section 3.4, the study of the electron power absorption based on the Boltzmann term method in section 3.5, and the simulation settings in section 3.6. The results are shown and discussed in section 4. The conclusions are drawn in section 5.

Experimental setup and discharge conditions
The space and time resolved optical emission of the discharge is measured experimentally, by PROES [19][20][21], in which the emitted light from the selected excited atomic/molecular state is spatio-temporally resolved. From that, the electron impact excitation rate from the ground state into the observed level, W 0,i (x, t), can be calculated, as introduced in [21]. In order to perform a successful PROES measurement on a CCP, some crucial conditions have to be fulfilled [21]: (i) Knowledge of the optical transition rates in the gas is needed. (ii) The contribution of cascades, excitation from metastable levels and quenching to the population of the measured excited state need to be negligibly low. (iii) The intensity of the emitted light at the measured line has to be high enough. (iv) No superposition with other optical lines is tolerated within the spectral resolution of the spectrometer or the interference filter. (v) The lifetime of the relevant excited state has to be short enough to temporally resolve the RF period (100 ns in the current study). All in all, the choice of the optical line to measure is critical. A typical line that satisfies these criteria well and is, therefore, often applied for PROES is the emission line corresponding to the Ne 2p 5 ( 2 P o 1/2 )3p 0 → 2p 5 ( 2 P o 1/2 )3s 1 transition, with a wavelength of 585.25 nm and lifetime of 16.26 ns [54]. Its threshold energy for electron impact excitation from the ground state is 18.966 eV [55]. For the current measurements, our geometrically symmetric "Budapest v.3" plasma reactor is used, which was first introduced in [22]. A sketch of the setup is shown in figure 1. The flat disk electrodes made of 304 grade stainless steel with identical diameters of 14.0 cm are situated within a quartz cylinder. Their separation is 2.5 cm. The upper electrode is driven by RF voltage, while the lower one is grounded. There is no active cooling applied on the electrodes. The chamber can be evacuated via a turbomolecular and a rotary pump, providing a base pressure of approximately 10 −5 Pa. The experiments are performed in a gas flow of ≈3 sccm set by two flow controllers. The driving voltage is provided by a waveform generator (Juntek JDS-2900), which is connected to a linear power amplifier (RM BLA-300) and an impedance matching box (MFJ-949E). The voltage drop across the discharge is measured by a high voltage probe (HP 10076A, 100:1) right at the power feeding point, and the pressure is monitored by capacitive gauges (Pfeiffer Vacuum CMR264 and MKS Baratron 631A). PROES measurements are performed with a fast-gateable ICCD camera (4 Picos, Stanford Computer Optics). The optical emission from the Ne 3p 0 state with a wavelength of 585.25 nm is measured. In order to limit the measured light to this specific line, an interference filter is applied with a central wavelength of 585 nm and a spectral full width at half maximum of ∼10 nm. The gate width of the camera is set to 2 ns. The camera is equipped with a Thorlabs MVTC23013 0.128x bi-telecentric lens, by which two-dimensional pictures can be taken. The spatial resolution is approximately 150 µm. Due to the lateral uniformity of the plasma, the data are averaged in the direction perpendicular to the discharge axis, which reduces the noise significantly.

Simulation method
The simulations are based on a one dimensional in space and three dimensional in velocity space (1d3v) Particle-in-Cell/Monte Carlo Collisions (PIC/MCC) simulation code [23][24][25][26][27][28][29][30][31][32][33][34]. This code (named PICit! ), recently developed in our group, is suitable to model geometrically symmetric capacitively coupled RF discharges in various mixtures of electropositive/electronegative gases. The particles traced in the simulations of neonoxygen gas mixture plasmas are electrons, Ne + ions, fast Ne atoms (Ne f ), O + 2 ions, O − ions and fast O 2 molecules (O f 2 ). The set of collision processes taken into account in the neon-oxygen discharge model is based on the sets of collision processes used previously for simulation studies of CCP discharges operated in pure neon [22] and oxygen [13,17,[50][51][52]56], complemented with "cross processes" between oxygen and neon species and collision processes for fast neutrals (see section 3.1). The heating of the background gas due to elastic collisions of fast neutrals and ions with thermal (background) atoms/molecules, as well as heating up of the electrodes due to inelastic collisions of plasma particles with the electrodes [57] is taken into account in the model (see section 3.3). The complexity of the model necessitates a relatively detailed description, which we provide below.

Collision processes
The elementary collision processes included in the discharge model are listed in table 1.

Electron collisions:
The collision processes between electrons and Ne atoms are elastic scattering (process 1), excitation (processes 2-10) and ionization (process 11). Nine electron impact Ne excitation processes are considered, including the Ne 3p 0 excitation from the ground state (process 10), whose population dynamics is captured experimentally by the PROES measurements. The cross sections for the e − + Ne collisions are adopted from the Biagi-v7.1 dataset [58] and are plotted in figure 2(a). Collisions between electrons and O 2 molecules include elastic scattering (process 12), rotational and vibrational excitation (processes 13-17), metastable excitation (processes 18 and 19), dissociative attachment (process 20), electronic excitation (process 21), dissociation (processes [22][23][24], ionization (process 25), and dissociative excitation (process 26). Electron impact detachment of O − ions (process 27) and dissociative recombination of O + 2 ions (process 28) are also considered. The e − + O 2 process list and cross sections adopted are largely based on the "xpdp1" set [60]. However, a few changes have been introduced [51]: we (i) replace the elastic collision (process 12) cross section with the elastic momentum transfer cross section of Biagi [59] and use, accordingly, isotropic electron scattering, (ii) replace the original "xpdp1" cross sections for the ionization (process 25), electron impact detachment (process 27) and dissociative recombination (process 28) with those recommended by Gudmundsson [61]. The cross sections for the electron impact collisions with oxygen species are plotted in figure 2(b) (processes 12-17) and figure 2(c) (processes 18-28).    and a backward scattering channel (process 30) is considered and the cross sections are taken from [62] in the functional forms σ 29 (ε com ) = 1.06 · 10 −19 ε −0.5 com m 2 and where ε com is the value of the kinetic energy in the ion-atom center-of-mass reference frame in units of eV. As Ne + + O 2 collisions, charge transfer (process 31) and isotropic elastic scattering (process 32) are included in the model. Process 31 is a non-resonant charge transfer reaction that is possible as the ionization potential of Ne is higher than that of O 2 . The rate coefficient of this reaction at thermal energies was reported in [63]. Schlumbohm [66] reported a measurement of this cross section in the energy range of 3-200 eV, and found it to be between 3.5 ×10 −20 m 2 and 4 ×10 −20 m 2 . We use a cross section of constant 3.7×10 −20 m 2 for this process. Process 32 is treated with its Langevin cross section: where µ is the reduced mass, g is the relative velocity and ε com is kinetic energy in the centre-of-mass frame in units of eV. The polarizibility values are: are also included in the model. The metastable O 2 (a 1 ∆ g ) singlet delta oxygen molecules are not traced in the simulation, however, their density is calculated self-consistently based on the balance of their creation, de-excitation, and diffusive transport as discussed in section 3.4 below. The cross sections of processes 36-38 are from [61], for process 39 the cross section is taken from [65]. As "cross-process" between O − ions and Ne species, isotropic elastic scattering with Ne atoms (process 40) and mutual neutralisation with Ne + ions (process 41) is considered. Process 40 is treated with its Langevin cross section. Process 41 is treated in a similar way as the mutual neutralisation O − + O + 2 −→ O + O 2 (process 38). Due to the lack of cross section data for O − + Ne + −→ O + Ne neutralization, we use the same cross section for process 41 as that used for process 38. The cross sections of the collision processes used for O − ions are plotted in figure 3 Fast neutral collisions: A background gas atom or molecule acquiring a kinetic energy greater than a threshold energy of E thr = 9 · 3 2 k B T after a collision is considered to be a fast neutral, and is traced in the simulation until its energy drops below E thr . Thermalization of fast neutrals due to collisions with atoms/molecules of the background gas contributes to the heating of the gas (see section 3.3). Elastic collisions of fast Ne atoms (Ne f ) and fast O 2 molecules (O f 2 ) with background O 2 molecules and Ne atoms (processes [42][43][44][45] are considered in the model in the following way: (i) the scattering cross sections are derived from σ V , the viscosity cross section and (ii) isotropic scattering in the COM frame is assumed. As reliable cross section data for Ne f and O f 2 for neonoxygen mixture are not available in the literature, the cross sections are calculated for processes 42-45 based on the pair-potential between the particles, for which the Lennard-Jones (LJ) type is assumed: where the energy LJ and the characteristic length σ LJ can be found in tables, e.g. in [68]. The calculation of the cross section starts with computing the scattering angle χ as a function of the impact parameter b. This is done as described in [69]. Having obtained the χ(b) function for a given LJ potential, the viscosity cross section is computed as [70]: The χ(b) function is computed for a set of energies (in the COM frame) and the above integral is evaluated for all these energy values, resulting in the energy-dependent elastic scattering cross section σ fast (ε) = 3 2 σ V (ε). To obtain the LJ parameters for the Ne-O 2 "cross collisions", the Lorentz-Berthelot combining rule is used [71]: where a and b denote the interacting species. The a ↔ b symmetry of these equations results in the same cross sections for Ne f +O 2 and O f 2 +Ne collisions. These computed cross sections (for processes [42][43][44][45] are plotted in figure 4. The calculation of the cross sections of fast neutrals was validated by comparing our results for argon gas with those given in [70].

Surface processes
The surface processes taken into account in the discharge model are electron reflection, secondary electron emission (SEE) induced by O + 2 ions and Ne + ions, and surface quenching of O 2 (a 1 ∆ g ) metastable molecules. Constant surface coefficients are specified for these processes. At the electrode surfaces, electrons are elastically reflected with a probability of η e = 0.7. This value is adopted based on the findings of Schulenberg et al. [53]. For O + 2 ions, a SEE coefficient of 0.015 is considered, while for Ne + ions the SEE coefficient is set to 0.1 [72]. The value of the surface quenching probability (surface destruction probability), which controls the loss rate of O 2 (a 1 ∆ g ) metastable molecules, is set to α = 8 · 10 −4 in this study (see discussion in section 3.4).

Gas heating
The collision frequency of elementary processes is strongly dependent on the local gas temperature (through the dependence on the background gas density). In order to cover the wide range of discharge parameters properly, gas heating needs to be included in the model. The coupling of the input RF electric power to the gas heating is realised by two processes: (i) in the gas phase, ions are accelerated by the electric field and energy is transferred to the background gas atoms/molecules through collisions, mediated by fast neutrals, and (ii) at the surfaces, the flux of particles interacting with the electrodes causes them to heat up.
Heat transfer: The thermal power density, P (x), accumulates the excess energy transferred to the background gas, originating from the particles accelerated by the electric field, per unit time and volume [73]. It acts as the source term in the steady-state heat equation, that is included in the model to calculate the equilibrium temperature distribution in the discharge gap: were κ is the thermal conductivity specific to the gas mixture. Analogously to the linear particle shape function in the PIC method, if a collision with the background gas happens within the i-th grid cell, it contributes to P (x) at the neighbouring grid points in the following way: where A is the surface of the electrodes, ∆x is the length of the grid cell normal to the electrodes, and ∆t H is the time of data collection (several RF cycles). ∆E is the energy gain of a background gas atom/molecule due to the collision. Collisions of ions, as well as fast neutrals with the background gas, contribute to P (x). When thermalized, the neutrals deposit all their energy to the P (x) thermal power.
Boundary condition: In order to solve eq. (6), the boundary values of the temperature distribution have to be defined. In regimes with Knudsen number (Kn) close but below unity the "slip flow" approximation can be applied [74]. The Knudsen number is defined as Kn = l/Λ, where l is the molecular mean free path and Λ is the characteristic size of the system, i.e. the length of the discharge gap in our case. The resulting boundary condition is of the third kind, introducing a temperature jump at the electrodes as described in [75,76]: where T wall is the temperature of the walls (i.e. the electrodes). λ is called the temperature jump distance, and is proportional to the mean free path perpendicular to the electrodes, and it depends on the gas mixture. For further details see [75]. Note that eq. (8) imposes conditions on the slope of the temperature distribution at the boundaries.
To satisfy these, the following iterative method is used: (i) Initially the temperature jump is neglected and wall temperatures (T wall ) are used as boundary values for the Thomas algorithm which solves eq. (6). (ii) Since the resulting temperature distribution does not satisfy eq. (8), the right hand side of eq. (8) is set as new boundary values to solve eq. (6) again. (iii) Step (ii) is repeated until condition (8) is fulfilled within pre-defined tolerance limits. The convergence of this method was demonstrated in [75]. Heat convection at the electrodes: The heating of the electrode surfaces is caused by inelastic collisions of plasma particles with the electrodes. In order to accurately describe the experimental system, this change of wall (electrode) temperature must be accounted for in the gas heating module, as it enters directly in eq. (8). The equilibrium wall temperature is determined by the balance of heat influx at the plasma side, thermal conduction of the electrode assembly, and heat dissipation to the environment at the outer surface of the system. In our calculations, we consider a simplified geometry of the electrode arrangement, as compared to the experimental realisation (consisting of the electrodes, spacers, and flanges), which is shown in figure 5. This simplification is permissible because of the use of a calibration procedure for the effective heat conductivity of the setting, to be revealed below. Letq 0 denote the heat flux originating from the discharge reaching the electrode and T ∞ the room temperature. k and L e are the thermal conductivity and the width of the electrode, respectively. h is the convection heat transfer coefficient of air. For the model system of figure 5, the steadystate homogeneous heat equation inside the electrode and its boundary conditions read: where x = 0 represents the plasma side and x = L e the external surface of the electrode assembly. Eqs. (10) and (11) are boundary conditions for boundaries at point 1. and 2. in figure 5. The solution of eq. (9) yields the wall temperature in terms of the heat flux from the plasma, assuming equalq 0 , k, and L e values for both electrodes: The value ofq 0 could, in principle, be calculated tracing individual particle-surface interactions, however, an even simpler and more precise estimation can be achieved by assuming that, in equilibrium, all electric power absorbed by the plasma (charged plasma particles only) from the RF field becomes finally absorbed by the electrodes and is dissipated in form of heat. This means that the heat flux reaching one electrode is the time average ( ... ) of the half of the total electrical input power per unit area: where the sum runs over all charged particles i with charge q i , weight factor W i , velocity v i and electric field at the particles' position E i . Having calculated this quantity, the wall temperature can readily be obtained from eq. (12), if (L e /k + 1/h) is known. While this factor (that consists of three components, each characterizing different aspects of the overall heat dissipation) could be approximated based on the exact dimensions and material properties of the electrode constructions, we have chosen an alternative way to determine it, using a calibration technique. Tunable diode laser absorption spectroscopy (TDLAS) measurements were performed in pure Ar discharges in the same experimental system. The details of this technique are described in [53]. The measured gas temperature was then compared with PIC/MCC results obtained for different values of the (L e /k + 1/h) coefficient, and its value giving the best agreement was determined and used throughout this study.

Oxygen metastable balance
Previous studies on oxygen CCPs have shown that at low pressures the discharge contains more O − ions than electrons, thus it is electronegative [46,48,50], which is quantified by the global parameter β = n O − / n e , i.e., the global electronegativity of the discharge. The O − ion density is determined by the balance of its creation (by means of electron impact dissociative attachment of O 2 molecules, process 20 in table 1) and the primary loss in the gas phase due to associative detachment in collisions with O 2 (a 1 ∆ g ) metastable molecules (process 39 in table 1). In order to accurately compute the O − ion density distribution, it is important to determine the density of the O 2 (a 1 ∆ g ) metastable species self-consistently. This is done by modeling the balance of creation, transport, and surface de-excitation of the O 2 (a 1 ∆ g ) molecules, which are incorporated as continuum species, described by the stationary diffusion equation where D is the diffusion coefficient, n m (x) and S m (x) are the density and time averaged (over several RF cycles) source rate distributions of the O 2 (a 1 ∆ g ) metastables. Like in the gas temperature calculation module, here we aim only for the stationary solution of the transport equation and assume that it does not depend on the actual time-evolution, which, on the other hand, can be significantly slower than the relaxation of other plasma parameters. The source term in eq. (14) is calculated analogously to the thermal power input in eq. (7) but with the microscopic contribution of process 18 in table 1 for every electron in each time step. The contribution of gas phase losses of metastables (through processes 39 in table 1) is not considered in the metastable density balance calculation. The validity of this simplification is supported by the collision statistical analysis, showing that the average probability of associative detachment is < 1% of the probability of metastable excitation for all conditions. To obtain proper pressure, temperature, and gas composition dependent diffusion coefficients we utilize the predictions of the Chapman-Enskog theory [77,78], resulting in the general form where D 0 is the gas composition dependent diffusion coefficient at a reference pressure The dominant loss process of O 2 (a 1 ∆ g ) metastables is the de-excitation at the electrodes, which is incorporated in the boundary conditions. Similarly to the gas heating calculation, due to the low pressure condition, where the collision mean free path is comparable to the characteristic size of the system, the boundary condition is of the third kind, relating the boundary value with its gradient as with κ = l(2 − α)/( √ 3α), where l is the collision mean free path of O 2 molecules, and α = 8 · 10 −4 is the surface destruction (de-excitation or quenching) probability. The mean free path is calculated from the diffusion coefficient using the formula derived from the kinetic theory of gases [79] With all terms evaluated eq. (14) is solved using the Thomas tridiagonal algorithm, the metastable density values at the boundaries n m (x = 0/L) are adjusted, and the solution is iterated until the boundary conditions (16) are met.

Electron power absorption
The characteristics of the electron power absorption have been investigated using the Boltzmann term method. This method, which provides spatio-temporally resolved information about the electric field and the power absorbed by the electrons, has been thoroughly described in earlier works [50,56,80]. Therefore, here only a short overview is given. The Boltzmann term analysis is based on the electron momentum balance equation, which can be rearranged for the electric field, which then can be divided into various terms, given by E tot = E in + E ∇p + E Ohm [56], where (n e u e ) + ∂ ∂x (n e u 2 e ) , Here m e and e are the electron mass and charge, respectively, n e is the electron density, u e is the mean velocity, Π c is the electron momentum loss and p denotes the diagonal element of the pressure tensor. Each of these electric field terms corresponds to a distinct physical mechanism: E in is the electric field term originating from inertial effects; E Ohm , the Ohmic electric field, is a consequence of electrons colliding with the particles of the background gas, and E ∇p is due to pressure effects. This electric field term can be split into two additional terms, according to E ∇p = E ∇n + E ∇T , where E ∇n = − T n e e ∂n e ∂x , where T = p /n e is the parallel electron temperature. E ∇n is, in regions where quasineutrality holds, identical to the "classical ambipolar" electric field [81] and E ∇T originates from the spatial gradient of the electron temperature. From the electric field terms the corresponding power absorption terms can be easily calculated by multiplying each of these terms with the spatio-temporally resolved electron conduction current density, j c . Table 2 lists the most important physical and numerical parameters of the simulations grouped according to the gas pressure. The physical parameters were chosen to cover a low to intermediate pressure regime and a wide mixing range. The numerical parameters were chosen to ensure the fulfillment of the usual stability and accuracy requirements imposed on the explicit electrostatic PIC/MCC scheme. For a recent detailed discussion on these criteria and the effect of different numerical parameters on the simulation results see e.g. [82]. Further parameters specific to surface processes, heat conductivity, and metastable molecule diffusion are introduced in the previous sections. Figure 6 shows the spatio-temporal distribution of the electron impact excitation rate from the ground state into the Ne 3p 0 state measured by PROES at six different pressures: 15 Pa (first row), 31 Pa (second row), 62 Pa (third row), 125 Pa (fourth row), 250 Pa (fifth row) and 500 Pa (sixth row). At a given pressure (panels in a given row), results obtained for different neon/oxygen concentration of the background gas mixture (with decreasing O 2 concentration in panels from left to right) are shown: 10% Ne-90% O 2 (first column), 30% Ne-70% O 2 (second column), 50% Ne-50% O 2 (third column), 70% Ne-30% O 2 (fourth column), 90% Ne-10% O 2 (fifth column). All panels of figure 6 cover one RF period on the horizontal axes, and the vertical axes show the distance from the powered electrode.

Experimental results
At the lowest pressure of 15 Pa (first row), excitation at both the expanding and the collapsing sheath edges can be observed, as well as in the bulk region. The strongest excitation is found at the bulk side of the collapsing sheath edge at both electrodes. At this pressure, measurements have been performed for O 2 concentrations between 50% and 90% in the background gas mixture. For all cases, the spatio-temporal distribution of the excitation rate exhibits similar patterns (panels (a1)-(c1)), independently of the Ne/O 2 mixing ratio, suggesting hybrid α-DA discharge operation mode (with dominant electron power absorption and excitation due to the DA-mode). The main excitation patterns observed at this pressure are labeled in panel (b1): I. indicates the excitation peak at the expanding sheath edge (α-peak), II. denotes the excitation in the central bulk region (drift feature), and III. signals the excitation at the collapsing sheath edge (ambipolar peak). At 31 Pa (second row), for O 2 concentrations decreasing from 90% to 30% (panels (a2)-(d2)), strong excitation at the expanding sheath edge (α-peak, I.) and significant excitation in the bulk region (drift feature, II.) with a weak excitation peak at the collapsing sheath edge (ambipolar peak, III.) are found, suggesting a hybrid α-DA discharge operation mode (with dominant α-mode). By decreasing the O 2 concentration at this pressure, the α-peak (I.) at the expanding sheath edge is enhanced. At the lowest O 2 concentration of 10%, the ambipolar peak (III.) at the collapsing sheath edge vanishes and the excitation plot suggests discharge operation in pure α-mode (see panel (e2) At the highest pressure of 500 Pa (sixth row), measurements have been carried out for up to 70% O 2 concentration. At the lowest O 2 concentration of 10%, a single excitation peak is seen at both electrodes at the expanding sheath edge (α-peak, I.), indicating discharge operation in pure α-mode (panel (e6)). Excitation in the bulk and development of two additional excitation peaks (IV. and V.) is visible starting from 30% O 2 concentration (see panels (b6)-(d6)). Compared to the 250 Pa case, at this high pressure, the development of patterns IV. and V. starts at a lower O 2 concentration (30% vs. 50% at 250 Pa). These two excitation patterns are intensified as the O 2 concentration is increased, clearly dominating the excitation at 70% and 50% O 2 concentrations (panels (b6) and (d6)), while only weak excitation in the close vicinity of the expanding sheath edge (α-peak, I.) can be observed in these cases. The observed mode transitions as a function of the pressure and O 2 admixture are related to changes of the electronegativity of the discharge, which will be discussed in detail later based on the simulation results.
In summary, at the lowest pressure, weak α-peak, strong ambipolar peak, and weak drift features are found in the excitation rate. With increasing pressure, the α-peak and the drift feature are enhanced, while the ambipolar peak is reduced. At intermediate pressures, the α-peak is the dominant excitation pattern in all mixtures. Further increase of the pressure leads to the formation of two distinct excitation peaks at the edges of the bulk region, which dominate the excitation at high O 2 concentrations. At the highest pressure, these excitation peaks are enhanced and the development of these features is found also in mixtures with lower O 2 concentration.

Simulation results
In order to reveal the physics behind the formation of the different excitation patterns observed experimentally by PROES in the wide pressure range and gas mixing range presented above, especially those found at high pressures and high O 2 concentrations in the bulk region, PIC/MCC simulations have been performed. The PIC/MCC simulations covered the whole parameter regime studied by PROES, i.e. pressures between 15 Pa and 500 Pa, and Ne/O 2 concentrations between 10% and 90%, at 10 MHz and 350 V peak-to-peak voltage.  the density of negative O − ions in the bulk. The electronegativity is very high in the discharge center, the ratio of the negative ion density and electron density is about 100 ( figure 8(a)). The global electronegativity (the ratio of the density of negative ions and electrons averaged over the electrode gap) is about 30 under these conditions (see figure 8(b)). The density of Ne + ions is about three orders of magnitude lower than the density of O + 2 ions in the discharge center and the density of fast neutral species is high over the whole discharge gap. In this gas mixture, local maxima in the time averaged electron density distribution at the edges of the bulk region and the high negative ion density in the bulk are found also at a higher pressure of 31 Pa (panel (b)). However, at 31 Pa the electron density is enhanced in the discharge center compared to the 15 Pa case, while the density of Ne + ions, as well as the density of fast neutrals, decreases in the bulk. The electronegativity is maximum in the discharge center ( figure 8(a)), the global electronegativity of the discharge is about 25 ( figure 8(b)). By increasing the pressure to 500 Pa for the 30% Ne-70% O 2 mixture (panel (c)), the time averaged electron density exhibits its maximum in the discharge center and local minima at the bulk edges. The density of Ne + ions as well as the density of fast neutrals drops in the bulk. The bulk region is wide at this high pressure. The O − density profile has a local minimum in the discharge center and local maxima near the edges of the bulk region. Consequently, the electronegativity has a local minimum in the center and local maxima at the bulk-sheath boundary ( figure 8(a)). In the discharge center, the density of negative ions is about one order of magnitude higher than the electron density. The global electronegativity is about 8 under these conditions ( figure 8(b)). By decreasing the O 2 concentration in the mixture to 10% at this high pressure ( figure 7(d)), the density of electrons increases, while the O − density remains high in the bulk region. The O − density is about four times higher than the electron density in the discharge center. The global electronegativity has a low value of about 2 under these conditions, i.e. the discharge remains electronegative also in case of the lowest O 2 concentration of 10% in the mixture. As it can be seen in figure 8(b), the global electronegativity decreases with increasing the pressure in all mixtures. For O 2 concentrations above 30% the electronegativity does not depend on the Ne/O 2 concentration ratio in the background gas mixture. The simulation results show that the self-consistently calculated density of singlet delta oxygen metastable molecules increases with the pressure. However, the O 2 (a 1 ∆ g ) density remains low in all cases, below 3% of the density of background O 2 molecules. Figure 9 shows the time averaged density distribution of O 2 (a 1 ∆ g ) molecules in the gap as a function of pressure for 30% Ne-70% O 2 gas mixture. At low pressures, the metastable density profiles are flat, while at high pressures, profiles with peak metastable density in the center of the discharge are obtained. This is due to the pressure dependence of the diffusion coefficient (see eq. 14, at high pressure the diffusion slows down). The change of the wall temperature as a function of pressure for various Ne/O 2 concentrations of the background gas mixture is shown in figure 10(a). The wall temperature, T wall , increases significantly compared to the initial wall temperature (which was set to 300 K in the simulations) with increasing the pressure. At a given pressure, T wall decreases with decreasing the O 2 content in the gas mixture. This is due to the different efficiency of the power absorption of Ne and O 2 at identical discharge conditions. In panel (b) of figure 10 the spatial distribution of the gas temperature with respect to the wall temperature (left vertical axis) and the power input field (right vertical axis) are shown at three different pressures for the 30% Ne-70% O 2 case. The temperature profiles show the maximum gas temperature in the discharge center at all pressures, which is only slightly higher than the temperature of the electrodes. Compared to the wall temperature, an increase by only about 1.5 K in the discharge center is found at the highest pressure. The simulations revealed that at the low electrical power levels of the present setup, being only a few Watts, the increased temperature of the gas is mainly due to the increased temperature of the electrodes, and the gas temperature does not significantly change within the gap. Figure 11 shows the spatio-temporal distribution of the electron impact excitation rate  their relative intensity). In the following, the simulation results are analysed in detail in four representative cases: (i) 15 Pa, 30% Ne-70% O 2 , (ii) 31 Pa, 30% Ne-70% O 2 , (iii) 500 Pa, 30% Ne-70% O 2 , and (iv) 500 Pa, 90% Ne-10% O 2 . For these cases, the spatio-temporal plots of the electron-impact excitation rate from the ground state into the Ne 3p 0 state are shown in panels (b1), (b2), (b6), and (e6) of figure 6 (PROES) and figure 11 (PIC/MCC simulation), respectively. Under these conditions, the main excitation patterns identified above (labeled I.-VI. in figure 6) are exhibited. The time averaged particle density distributions and the electronegativity obtained for these four cases were presented in figure 7 and figure 8(a), respectively, and discussed above.  In panel (c) the temporally averaged negative ion (O − ) density divided by a factor of 50 is also included (dashed black line) to illustrate that the density of O − ions is much higher than the electron density in the bulk region. As it can be seen in figure 7(a), the time averaged O − density is about two orders of magnitude higher than the time averaged electron density in the discharge center. As a consequence of this, the ratio of the O − ion and electron density, i.e. the electronegativity is high (about 100) in the discharge center ( figure 8(a)). Under these conditions the ratio of negative ions and electrons averaged over the electrode gap, i.e. the global electronegativity is also high (about 30, see figure 8(b)). The spatio-temporal plots of the Ohmic power absorption, P Ohm , and the ambipolar power absorption, P ∇n , are shown in panels (d) and (e), respectively. Due to the high electronegativity in the discharge center, the conductivity of the plasma is low in the bulk region. Consequently, the Ohmic power absorption, P Ohm , peaks in the discharge center at the times of maximum RF current within the RF period, as it can be seen in panel (d). The ambipolar power absorption, P ∇n , peaks at the edges of the bulk region at the time of sheath collapse at both electrodes (see panel (e)) due to the local maximum of the electron density and the corresponding large gradients in the electron density (panels (b) and (c)). Under these conditions, the ambipolar power absorption is the dominant power absorption mechanism. Panel (f) shows temporal snapshots of the Ohmic and ambipolar electric field, E Ohm and E ∇p , respectively, at a selected time instance indicated by black vertical dashed lines in panels (d) and (e). At this time, the Ohmic electric field has a negative sign in the bulk, while the ambipolar electric field changes the sign three times. The sum of these two terms results in a strong electric field at the grounded electrode side of the bulk (at the instantaneous collapsing sheath side) and a weak, almost zero electric field at the powered electrode side of the bulk (at the instantaneous expanding sheath side). Electrons accelerated by these electric fields induce strong excitation (III.) at the collapsing sheath edge and weaker excitation (II.) in the discharge center (see panel (a)). Weak excitation at the expanding sheath edge (I.) can be also observed. These are excitation patterns characteristic of a hybrid α-DA discharge operation mode with dominant DA-mode, specific to electronegative discharges (panel (a)). Discharge characteristics obtained from the simulations performed for the same gas mixture as in case of Figure 12, i.e. 30% Ne-70% O 2 , but at a higher pressure of 31 Pa, are shown in Figure 13. In this case, the spatio-temporal distribution of the excitation rate is dominated by the α-peaks (I.) in combination with drift features (II.) in the bulk and ambipolar peaks (III.) at the collapsing sheath edges. Similarly to the 15 Pa case, the electron density profiles plotted at different time instances show strong peaks at the edges of the bulk region (panel (c)). The electron density is low in the bulk (panels (b) and (c)), however it is enhanced in the discharge center compared to the 15 Pa case. Similarly to the 15 Pa case, the electronegativity is maximum in the discharge center ( figure 8(a)) and the global electronegativity has a high value of about 25 ( figure 8(b)). The Ohmic power absorption peaks in the discharge center (panel (d)) and its magnitude is similar to that obtained at 15 Pa. The ambipolar power absorption peaks at the edges of the bulk region (panel (e)) and its magnitude is significantly reduced compared to the 15 Pa case. The superposition of the Ohmic electric field and the ambipolar electric field at the selected time instance (panel (f)) results in a strong electric field at the powered electrode side of the bulk (at the instantaneous expanding sheath side) as well as at the grounded electrode side of the bulk (at the instantaneous collapsing sheath side). The electrons accelerated by these electric fields induce the strong excitation peak (I.) at the expanding sheath edge and the weaker excitation peak (III.) at the collapsing sheath edge. The nonzero electric field in the discharge center leads to excitation in the central bulk region (feature II.) These are excitation patterns characteristic of a hybrid α-DA discharge operation mode with dominant α-mode (panel (a)). The discharge characteristics obtained at a very high pressure of 500 Pa in 30% Ne-70% O 2 gas mixture are shown in figure 14. Under these conditions a weak α-peak (I.) and two additional excitation peaks (IV. and V.) in about the same time interval are found in both halves of the RF period, as well as weak excitation patterns (VI.) in the sheaths (panel (a)). The electron density distributions at different time instances show a local minimum close to the bulk edges, at about the position of maximum sheath expansion at both electrodes and maximum in the center of the discharge (panels (b) and (c)). The time averaged density of negative ions has a local minimum in the center and local maxima at the borders of the sheath region (note that the time averaged O − density is divided by 50 in panel (c)). The reason for this is the increasing locality of the charged particle transport as a function of pressure. During sheath expansion, energetic electrons are created at the expanding sheath edge, which induce a local formation of O − ions via high-energy threshold dissociative attachment (reaction 20 in table 1). At low pressures, these processes are not localized at the sheath edge. As a consequence of these characteristics, the electronegativity has a local maximum at about the positions of local minima in the electron density ( figure 8(a)), which leads to the development of drift electric fields in these regions. These are enhanced by the bulk electric field induced by the high collisionality of the plasma at this high pressure. Therefore, both the collisionality due to the high pressure and the electronegativity of the discharge contribute to the ohmic power absorption term. As a result of the superposition of these two effects, the Ohmic power absorption peaks at the edges of the bulk (panel (d)), in the same regions where the excitation patterns IV. and V. are found. Compared to the Ohmic power absorption, the ambipolar power absorption (panel (e)) is less significant under these conditions. This can be seen also in panel (f) which shows the Ohmic and ambipolar electric fields at the selected time instance: the ambipolar electric field is weak compared to the Ohmic field. The Ohmic field plotted at the time instance indicated in panel (d) shows local maxima at the bulk edges which are responsible for the excitation peaks IV. and V. at the bulk edges. The simulation results obtained at the highest pressure of 500 Pa by increasing the Ne concentration to 90%, i.e. for the 90% Ne-10% O 2 gas mixture, are shown in figure 15. Under these conditions strong excitation is found only at the expanding sheath edges (pattern I.) at both electrodes (panel (a)). The electron density profiles for the different time instances show only small local minima at the edges of the bulk region (panel (c)). The time averaged O − density has a local minimum in the center and local maxima at the borders of the sheath region. The discharge remains electronegative (the global electronegativity is about 2, see figure 8(b)). The spatio-temporal distribution of the Ohmic power absorption (panel (d)) shows peaks at the edges of the bulk, but these are less pronounced compared to the case of high O 2 concentration at the same pressure. The ambipolar power absorption (panel (e)) is mainly concentrated at the expanding sheath edges and has about the same magnitude as the Ohmic power absorption. At this low O 2 concentration of 10% in the mixture, there are no local maxima in the Ohmic electric field at the bulk edges at the selected time instance, and the ambipolar electric field exhibits only slight modulation in the gap (panel (f)). Under these conditions, the discharge operates in α-mode. In summary, based on the PIC/MCC simulation results, the mechanisms behind the development of the various excitation features revealed by the PROES measurements in Ne-O 2 CCPs in the wide pressure range for different gas mixing ratios have been successfully revealed. Finally, figure 16 shows the spatio-temporal distribution of the electron impact ionization rate of oxygen obtained from PIC/MCC simulations for different pressures and different Ne-O 2 mixtures (for the same conditions as those covered in case of the PROES results presented in figure 6 and the excitation rates in figure 11). The ionization patterns seen in panels of figure 16 reflect the excitation patterns of figure 11. This means that the PROES measurement results in mixtures of Ne and O 2 under the conditions studied probe both the excitation and the ionization dynamics in the discharge. Under these conditions, the PROES measurements provide correct information on the discharge operation mode.

Conclusions
We have performed Phase Resolved Optical Emission Spectroscopy (PROES) measurements combined with 1d3v Particle-in-Cell/Monte Carlo Collisions (PIC/MCC) simulations in low-pressure capacitively coupled neon-oxygen gas mixture plasmas. This study covered a wide pressure range and a wide mixing range of Ne and O 2 gases for a geometrically symmetric plasma reactor with a gap length of 2.5 cm, operated at a driving frequency of 10 MHz and a peak-to-peak voltage of 350 V. The pressure of the gas mixture was varied between 15 Pa and 500 Pa, while the neon/oxygen concentration was changed between 10% and 90%. For all discharge conditions, the spatio-temporal distribution of the electron-impact excitation rate from the Ne ground state into the Ne 3p 0 state was recorded by PROES. The measured electron-impact excitation rate from the Ne ground state into the Ne 3p 0 state was compared to the PIC/MCC simulation results on the Ne excitation rate, resulting in a good qualitative agreement in the whole parameter regime. At the lowest pressure, a weak α-peak at the expanding sheath edge, strong ambipolar peak, and weak drift feature in the bulk region were found in the excitation rates. With increasing pressure, the α-peak and the drift feature were found to be enhanced, while the ambipolar peak was reduced. At intermediate pressures, the α-peak was found to be the dominant excitation pattern in all mixtures. Further increase of the pressure resulted in the formation of two distinct excitation peaks at the edges of the bulk region, which dominated the excitation at high O 2 concentrations. Based on the emission/excitation patterns, multiple discharge operation regimes were identified. It was found that the localized bright emission features at the bulk boundaries at high pressures and high O 2 concentrations are caused by local maxima in the electronegativity. The relative contributions of the ambipolar and the Ohmic electron power absorption were found to vary strongly with the discharge parameters: the Ohmic power absorption was enhanced by both the high collisionality at high pressures and the high electronegativity at low pressures. In the wide parameter regime covered in this study, the PROES measurements were found to accurately probe the ionization dynamics in the discharge, i.e. the discharge operation mode. The simulation revealed that the temperature of the electrodes increases significantly compared to the initial wall temperature with increasing the gas pressure. It was found that the power deposition within the gas causes only a slight increase of the gas temperature above the temperature of the electrodes, which was, however, found to be significant due to the heating of the electrodes by the particles from the plasma. This finding points out the importance of the thermal balance of the electrode construction in determining the electrode and gas temperatures under operating conditions at moderate electrical power levels.