1. Introduction
Low pressure radio frequency (rf)-driven capacitively-coupled discharges have a range of material processing applications such as plasma etching and plasma enhanced chemical vapor deposition within the microelectronics industry. These discharges have been explored extensively over the past few decades. However, a few issues remain to be fully understood, including the electron heating mechanism, in particular when driven by multiple frequencies [
1], and the role of surfaces regarding recombination and quenching of various species and phenomena such as secondary electron emission and electron reflection [
2,
3]. The modern capacitively-coupled discharge consists of two parallel electrodes separated by a few cm and is driven by a radio-frequency power generator. The plasma forms when rf voltage is applied between the electrodes. The electrons that gain enough energy from the resulting electric field produce positive ions, negative ions, and electrons through electron impact ionization of neutral atoms and molecules and electron impact dissociative attachment of molecules, which forms the plasma. The plasma is separated from the electrodes by space charge sheaths. Multiple frequencies are commonly applied in order to achieve separate control of ion flux and ion energy, as the ion flux dictates the throughput of the process and the ion energy determines the etching and deposition parameters on the wafer surface.
The particle-in-cell (PIC) method, when combined with Monte Carlo (MC) treatment of collision processes, is a self-consistent kinetic approach that has become a predominant numerical approach to investigate the properties of the low pressure capacitively-coupled discharge. This approach is commonly referred to as particle-in-cell Monte Carlo collision (PIC/MCC) method. The basic idea of the PIC method is to allow typically a few hundred thousand computer-simulated particles (superparticles) to represent a significantly higher number of real particles (density in the range of
–
m
) [
4,
5,
6]. In a PIC simulation, the motion of each particle is simulated and the various macro-quantities are calculated from the position and velocity of these particles. The particle interaction is handled through a macro-force acting on the particles, which is calculated from the field equations at points on a computational grid. This method allows us to follow the spatio-temporal evolution of the various plasma parameters such as particle density, particle energy, particle fluxes, and particle heating rates.
The kinetics of the capacitively-coupled oxygen discharge have been studied for over two decades starting with the seminal work of Vahedi and Surendra [
7] using the 1D
xpdp1 PIC/MCC code. Since then, a number of PIC/MCC studies have been reported on oxygen and Ar/O
discharges using the
xpdx1 series of codes, in both symmetrical and asymmetrical geometry, performed over a range of pressures and compared to experimental findings [
8] and to analytical density profiles [
9], showing good agreement, to explore the formation of the ion energy distribution function in an O
/Ar mixture in an asymmetric capacitively-coupled discharge [
10], and the influence of the secondary electron emission on the density profiles and the electron energy distribution function (EEDF) [
11]. Other 1D PIC/MCC codes have been developed to explore the oxygen discharge. A 1D PIC/MCC model developed in Greifswald, that includes the metastable oxygen molecule O
as a fraction of the ground state molecule, was used to determine the ion energy distribution function (IEDF) in oxygen CCP [
12,
13]. Furthermore, they found by comparison with experiments that one sixth of the oxygen molecules are in the metastable singlet delta state. A 1D PIC/MCC code, developed in Dalian [
14,
15], was applied to explore the electrical asymmetry effect in a dual-frequency capacitively-coupled oxygen discharge. Similar to Bronold et al. [
12], this work assumed a constant density for the singlet metastable molecule O
. More recently, a 1D PIC/MCC code that was developed in Budapest was used to explore the heating mechanism in a capacitively-coupled oxygen discharge driven by tailored waveforms (composed of
N harmonics in addition to a fundamental frequency
) [
16,
17]. Furthermore, a PIC/MCC fluid hybrid model was applied to explore the electron power absorption and the influence of pressure on the energetics and particle densities [
18,
19]. In all of these works, only electrons, the positive ion O
, and the negative ion O
were treated kinetically, and the positive ion O
was neglected. Furthermore, none of the metastable states were treated kinetically. The one-dimensional object-oriented plasma device one (
oopd1) code allows having the simulated particles of different weights, which allows for tracking both charged and neutral particles in the simulation. Earlier, we benchmarked the basic reaction set for the oxygen discharge in
oopd1 to the
xpdp1 code [
20].
In recent years, the oxygen reaction set in the
oopd1 code was improved significantly [
20,
21,
22]. Using this improved discharge model, we showed that the singlet metastable molecular states have a significant influence on the electron heating mechanism in the capacitively-coupled oxygen discharge [
21,
22,
23,
24] as well as the ion energy distribution [
25]. We demonstrated that, when operating at low pressure (10 mTorr), the electron heating is mainly located within the plasma bulk (the electronegative core), while, when operating at higher pressures (50–500 mTorr), the electron heating appears almost solely within the sheath regions [
22,
23]. Furthermore, when operating at low pressure, the electron heating within the discharge is due to a hybrid drift-ambipolar-mode (DA-mode) and
-mode, and while operating at higher pressures, the discharge is operated in a pure
-mode [
26,
27]. We have also shown that 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 almost negligible influence at lower pressures [
22,
23,
24].
Secondary electron emission and electron reflection from the electrodes have often been neglected in PIC/MCC simulations. When it is included, it is common to assume the secondary electron emission to have a constant value (independent of the discharge conditions such as the energy of the bombarding ions), while only the ion-induced secondary electron emission is taken into account, and thus, the contributions of other species are neglected [
2,
3]. The effects of including a constant secondary electron emission yield are increased electron density, enhancement of the density profiles and the electron energy distribution functions (EEDFs), decreased sheath width, and the electron heating rate profiles changing significantly in both argon [
28] and oxygen [
11] discharges. Furthermore, it has been demonstrated that an asymmetry can be introduced by having electrodes with different secondary electron emission properties in a capacitively-coupled discharge [
29], which was later extended to also include the electrical asymmetry effect in a dual-frequency capacitively-coupled discharge driven by two consecutive harmonics with different electrode materials [
30]. In these studies, the secondary electron emission yield was set to be a constant. A few recent studies have emphasized using realistic secondary electron emission yields for both fast neutrals and ions bombarding the electrodes [
2,
3,
22,
28,
31,
32,
33].
In an earlier study, we explored the role of including an energy-dependent secondary electron emission yield for both O
and O
-ions and O and O
neutrals in an oxygen discharge [
22]. We noted that this had a significant influence on the discharge properties, including increased electron and ion densities and decreased sheath width. Here, we study systematically how the secondary electron emission and the electron reflection from the electrodes influence the charged particle profiles, the electron heating processes, the electron energy probability function (EEPF), and the effective electron temperature, in a single frequency voltage-driven capacitively-coupled oxygen discharge by means of numerical simulation, for a fixed discharge voltage, while the discharge pressure is varied from 10–50 mTorr. The simulation parameters and the cases explored are defined in
Section 2, and the simulation results found by including and excluding the ion-induced secondary electron emission and electron reflection are compared in
Section 3. We give a summary and concluding remarks in
Section 4.
2. The Simulation
The one-dimensional (1d-3v) object-oriented particle-in-cell Monte Carlo collision (PIC/MCC) code
oopd1 [
34,
35] is herein applied to a capacitively-coupled oxygen discharge. The
oopd1 code, like the well-known
xpdp1 code [
7], is a general plasma device simulation tool capable of simulating various types of plasmas, including breakdown, accelerators, beams, as well as processing discharges [
20].
The oxygen reaction set included in the
oopd1 code is rather extensive. Like
xpdp1, it includes the ground state oxygen molecule O
, the negative ion O
, the positive ion O
, and electrons [
7,
20]. In addition, oxygen atoms in the ground state O
P) and ions of the oxygen atom O
[
20], the singlet metastable molecule O
(a
), and the metastable oxygen atom O
D) [
21], and the singlet metastable molecule O
(b
) [
22] were added along with the relevant reactions and cross-sections. The full oxygen reaction set was discussed in our earlier works where the cross-sections used were also given [
20,
21,
22]. Furthermore,
oopd1 has energy-dependent secondary electron emission coefficients for oxygen ions and neutrals as they bombard both clean and dirty metal electrodes [
22]. Thus, for this current work, the discharge model contains nine species: electrons, the ground state neutrals O(
P) and O
, the negative ions O
, the positive ions O
and O
, and the metastables O(
D), O
, and O
(b
. We herein use the secondary electron emission yield for a dirty surface as given in our earlier work [
22].
We assume a geometrically-symmetric capacitively-coupled discharge where one of the electrodes is driven by an rf voltage at a single frequency:
while the other electrode is grounded. Here,
is the voltage amplitude,
f the driving frequency, and
t the time. The discharge operating parameters assumed are the voltage amplitude of
V, an electrode separation of 4.5 cm, and a capacitor of 1 F connected in series with the voltage source. The driving frequency is assumed to be 13.56 MHz. These are the parameters used in our earlier works using
oopd1 [
20,
21,
22,
24,
27] and in the work of Lichtenberg et al. [
9] using the
xpdp1 code. The discharge electrode separation is assumed to be small compared to the electrode diameter so that the discharge can be treated as one dimensional. We assume the electrode diameter to be 10.25 cm, which is needed in order to determine the absorbed power, and set the discharge volume for the global model calculations applied to determine the partial pressure of the neutral species. The time step
and the grid spacing
are set to resolve the electron plasma frequency and the electron Debye length of the low-energy electrons, respectively, according to
, where
is the electron plasma frequency, and the simulation grid is taken to be uniform and consists of 1000 cells. The electron time step is set to
s. The simulation was run for
time steps, which corresponds to 2750 rf cycles. It takes roughly 1700 rf cycles to reach equilibrium for all particles, and the time averaged plasma parameters shown, such as the densities, the electron heating rate, and the effective electron temperature, are averages over 1000 rf cycles. All particle interactions are treated by the Monte Carlo method with a null-collision scheme [
4]. For the heavy particles, we use sub-cycling, and the heavy particles are advanced every 16 electron time steps [
36]. Furthermore, we assume that the initial density profiles are parabolic [
36].
The kinetics of the charged particles (electrons, O
-ions, O
-ions, and O
-ions) was followed for all energies. Since the neutral gas density is much higher than the densities of charged species, the neutral species at thermal energies (below a certain cut-off energy) are treated as a background with fixed density and temperature and maintained uniformly in space. These neutral background species are assumed to have a Maxwellian velocity distribution at the gas temperature (here,
= 26 mV). The kinetics of the neutrals are followed when their energy exceeds a preset energy threshold value. The energy threshold values and the particle weights used here for the various neutral species included in the simulation are listed in
Table 1. The partial pressures of the background thermal neutral species were calculated using a global (volume averaged) model of the oxygen discharge, as discussed in Proto and Gudmundsson [
37]. The fractional densities for the neutrals O
, O(
P), O
(a
, and O
(b
, estimated using the global model calculations at 10, 25, and 50 mTorr, are listed in
Table 2. These values are used as input for the PIC/MCC simulation as the partial pressures of the neutral background gas. Note that not all the neutrals considered in the global model calculations are shown in
Table 2. Due to recombination of atomic oxygen and quenching of metastable atoms and molecules on the electrode surfaces, discussed below, there is a drop in the high energy (energy above the threshold value) atomic oxygen density and an increase in the high energy oxygen molecule densities next to the electrodes, as shown in our earlier work [
22]. Thus, assuming uniformity of the background gas is thus somewhat an unrealistic assumption.
The electrode surfaces have significant influence on the discharge properties. There are a few parameters regarding the surface interaction of the neutral species that have to be set in the discharge model. For a neutral species that hits the electrode, we assume it returns as a thermal particle with a given probability. Similarly atoms can recombine on the electrode surfaces to form a thermal molecule with a given probability. As the oxygen atom O(
P) hits the electrode, we assume that half of the atoms are reflected as O(
P) at room temperature, and the other half recombines to form the ground state oxygen molecule O
at room temperature. Thus, for a neutral oxygen atom in the ground state O(
P), we use a wall recombination coefficient of 0.5, as measured by Booth and Sadeghi [
38], for a pure oxygen discharge in a stainless steel reactor at 2 mTorr. Similarly, as the metastable oxygen atom O(
D) hits the electrode, we assume that half of the atoms are quenched to form O(
P) and that the other half recombines to form the ground state oxygen molecule O
at room temperature. For the surface quenching coefficients of the singlet metastables molecules on the electrode surfaces, we assume for the singlet metastable O
(a
) a value of
, and for the singlet metastable O
(b
), we assume a value of
, based on the suggestion by O’Brien and Myers [
39] that the surface quenching coefficient for the b
state is significantly larger than for the a
state. We explored the influence of the surface quenching coefficients of the singlet metastable molecule O
(a
) on the discharge properties in an earlier work [
37]. There, we demonstrated that the influence of
on the discharge properties and the electron heating mechanism can be significant indeed. The partial pressures listed in
Table 2 were calculated by a global model using these surface quenching and recombination parameters as discussed in our earlier study [
37].
In the simulations, we either neglect electron reflection from the electrode or assume that electrons are reflected from the electrodes with a probability of 0.2, which is the number of elastically-reflected electrons per incoming electron, independent of their energy and angle of incidence. This value is based on the summary of values presented by Kollath [
40] for various materials. This value has been used by others in PIC/MCC simulations of capacitively-coupled discharges [
2,
3]. However, in reality, the reflection of electrons is known to depend on the electrode material, incident electron energy and the angle of incidence [
40,
41]. Furthermore, for all the cases explored here, we neglect secondary electron emission due to electron impact of the electrodes. The four cases explored for each pressure are listed in
Table 3.
3. Results and Discussion
The choice of the surface quenching coefficient for the singlet metastable O
(a
) of
was based on our earlier study of the time averaged electron heating profile between the electrodes
[
37]. In this study, we found that at 10 mTorr, almost all the electron heating occurred within the plasma bulk (the electronegative core), and the electron heating profile was almost independent of the surface quenching coefficient for the singlet metastable molecule O
(a
, while the DA-heating mode dominated the time averaged electron heating over one rf cycle. At 25 mTorr, the time averaged electron heating occurred both in the bulk (the electronegative core) and in the sheath regions, and a hybrid DA- and
-mode heating was observed. When operating at 50 mTorr, electron heating in the sheath region dominated, and the discharge was operated in a pure
-mode. Thus, this choice of pressure values and
gave us three distinct operating regimes to analyze further.
The electron energy probability function (EEPF) in the discharge center is shown in
Figure 1, for the various combinations of pressures, including and excluding secondary electron emission and electron reflection from the electrodes, for a total of four cases for each pressure, as shown in
Table 3.
Figure 1a shows the electron probability function (EEPF) at 10 mTorr. At low electron energy, the EEPF curved outwards, and a high energy tail was apparent when secondary electron emission was excluded from the simulation. We see that adding secondary electron emission to the discharge model enhanced the EEPF. When including the ion-induced energy-dependent secondary electron emission yield, more electrons were created at the electrodes, which were subsequently accelerated to the plasma bulk across the sheath. Thus, more high energy electrons were created in the discharge, and the EEPF exhibits a high energy tail when secondary electrons were emitted from the electrodes. This high energy tail extended up to roughly 240 eV. At 10 mTorr, both cases (including and excluding electron reflection) including secondary electron emission overlapped, and both cases (including and excluding electron reflection) neglecting the secondary electron emission overlapped. Thus, including electron reflection from the electrodes had negligible effects on the EEPF.
Figure 1b shows the EEPF at 25 mTorr. We see that the shape changed for all four cases as the pressure increased. The electron reflection had a negligible effect on the EEPF. This can be seen from the overlap of the green dashed line on the black one, when secondary electron emission was excluded, and from the overlap of the red line on the blue one, where secondary electron emission was included. Furthermore, an overall reduction of the high energy part of the curve, compared to the 10 mTorr case, was observed when secondary electron emission was neglected. This means that, when the pressure was raised and the secondary electron emission was neglected, there were fewer hot electrons within the bulk.
Figure 1c shows the EEPF at 50 mTorr. Here, the transition, which already started at 25 mTorr, was fully accomplished and the shape of the EEPF now curved inwards or was bi-Maxwellian for all four cases. As before including secondary electron emission led to a high energy tail. Furthermore, now, the black dashed line with the green one and the blue dashed line with the red one overlapped almost perfectly, which indicates that the electron reflection from the electrodes had negligible effects.
Figure 2 shows the profile of the time averaged power absorption by the electrons over an rf cycle
. A predominance of the electron heating within plasma bulk was observed in all four cases at 10 mTorr. We see that, when including the ion-induced secondary electron emission, the difference between including and excluding the electron reflection at the electrodes was very small within the plasma bulk. The same occurred when the secondary electron emission was excluded. A maximum in the power absorption in the bulk and a minimum in the sheath region were observed when the secondary electron emission was excluded and the electron reflection was included in the simulation (black line in
Figure 2a). On the contrary, a maximum in the power absorption in the sheath edge and a minimum in the bulk were seen when the secondary electron emission was included and the electron reflection was excluded (red line in
Figure 2a). At the transition pressure of 25 mTorr, the situation was drastically changed. A combination of the electron heating in the plasma bulk and in the sheath region was observed in all four cases. Indeed
Figure 2b shows that there was a maximum in the power absorption in the bulk and a minimum in the sheath edge when both secondary electron emission and electron reflection were excluded (green line in
Figure 2b), while a minimum in the bulk and a maximum in the sheath edge were observed when the secondary electron emission and the electron reflection were included (blue line in
Figure 2b). At 50 mTorr, the transition was fully accomplished and the electron heating was almost solely in the sheath region or stochastic electron heating. This is clearly seen in
Figure 2c when averaged over the rf cycle. Indeed, the electron reflection did not play much of a role. The maximum in the power absorption was observed when secondary electron emission was included in the simulation and the sheath was slightly narrower.
In order to explore the observed transition further, we plot the center electron density as a function of pressure in
Figure 3a. The electron density increased with increased pressure. At 10 mTorr, all four cases exhibited a similar electron density, and the electron density was slightly enhanced when the electron reflection was included in the simulation. At 25 mTorr, we see that including both the secondary electron emission and the electron reflection gave the highest center electron density, while excluding both processes led to the lowest center electron density. The differences in electron density between including and excluding both secondary electron emission and electron reflection were bigger than at 10 mTorr. At 50 mTorr, we see that including the ion induced secondary electron emission increased the center electron density and that including electron reflection increased the electron density even further.
Further insights about the observed transition are shown in the plot of the center electronegativity as a function of pressure in
Figure 3b. At 10 mTorr, the discharge was the most strongly electronegative. The electronegativity decreased from ∼110 at 10 mTorr to ∼20 at 50 mTorr. The electronegativity was higher (lower) when electron reflection was excluded (included) in the simulation; however, all four cases were very close to each other. The maximum (minimum) value of the electronegativity was reached when both secondary electron emission and electron reflection were excluded (included). At 25 mTorr, we observe that the gap between including and excluding both secondary electron emission and electron reflection was bigger than at 10 mTorr. We observe that, when secondary electron emission was included, excluding the electron reflection enhanced the electronegativity. The same occurred when secondary electron emission was included. Indeed, in this case, excluding both secondary electron emission and electron reflection gave the highest electronegativity. At 50 mTorr, the electronegativity was drastically reduced. We observed that electronegativity was lowest when secondary electron emission was included in the simulation and that electron reflection did not play much of a role. On the other hand, the electronegativity was highest when secondary electron emission and electron reflection were excluded from the simulation.
Figure 4 shows the spatio-temporal behavior of the electron power absorption
, where
and
are the spatially and temporally-varying electron current density and electric field, respectively. The figures show the electron power absorption for the various combinations of pressures, including and excluding secondary electron emission, while excluding electron reflection from the electrodes. For each of the figures, the abscissa covers the whole inter-electrode gap, from the powered electrode on the left-hand size to the grounded electrode on the right-hand size. Similarly, the ordinate covers the full rf cycle. Note that each of the six figures may have different magnitude scales, represented by the color scales on the right-hand side of each figure. Therefore, there can be differences in the six figures, not only qualitative, but also quantitative.
Figure 4a,b shows the spatio-temporal behavior of the electron power absorption including and excluding
at 10 mTorr, respectively.
Figure 4c shows the difference between including and excluding the ion-induced secondary electron emission from the electrodes. In
Figure 4a,b, the most significant heating is observed in the sheath region, during the sheath expansion, and the most significant cooling is observed during the sheath collapse. Here, significant energy gain (red and yellow areas) and small energy loss (dark blue areas) were evident within the plasma bulk region. We observe electron heating during the sheath collapse on the bulk side of the edge of the collapsing sheath (next to the instantaneous anode), while there was cooling (electrons loose energy) on the electrode side (the lower left-hand corner and upper center on the right-hand side). This kind of electron heating structure was observed experimentally in a capacitively-coupled SF
/N
discharge [
42] and SiH
discharge [
43] using spatiotemporal optical emission spectroscopy. This heating mechanism was explored further using the relaxation continuum model [
44], where electron heating due to three processes was identified: sheath expansion (
-mode), high electric field within the bulk, and ionization due to formation of a double layer on the instantaneous anode side, which resulted in acceleration of electrons. Indeed, in electronegative discharges, this electron heating within the plasma bulk can be the dominating electron heating mechanism [
42,
44]. This heating mechanism, which is due to electrons that are accelerated by strong drift and ambipolar electric fields within the plasma bulk and at the sheath edges in strongly electronegative discharges, was later coined as drift ambipolar (DA) electron heating [
45]. In highly electronegative discharges, these electrons are often found to dominate the ionization processes. As seen in
Figure 3b, the electronegativity was high ∼110 at 10 mTorr, which is essential for the DA-heating to be effective. The electron heating occurred both within the bulk and in the sheath regions, and a hybrid DA- and
-mode heating was observed. By looking at the time averaged electron heating profile in
Figure 2a, we see that there was electron cooling in the sheath region and all the electron heating occurred in the bulk region averaged over one rf cycle.
Figure 4c shows that there were no significant differences in the power absorption between the two cases, and in fact, there was a slightly higher electron heating within the discharge when secondary electron emission was excluded.
Figure 4d,e shows the spatio-temporal behavior of the electron power absorption for
and
, respectively, at 25 mTorr.
Figure 4f shows the difference between including and excluding secondary electron emission from the electrodes. At 25 mTorr, a transition process was observed. Indeed,
Figure 4e shows that the heating and the cooling in the sheath regions were reduced, while
Figure 4d shows that a significant contribution to the electron heating in the bulk region was observed. Therefore, a hybrid DA- and
-mode heating was observed, where the DA-heating was more important when secondary electron emission was excluded (
Figure 4e) than when it was included (
Figure 4d). This is clearly seen in the difference plot shown in
Figure 4f when cooling was seen as the difference. We see in
Figure 2b that in this case, the time averaged power absorption was observed both in the plasma bulk, as well as in the sheath regions.
Figure 4g,h shows the spatio-temporal behavior of the electron power absorption for
and
at 50 mTorr, respectively.
Figure 4i shows the difference between including secondary electron emission (
Figure 4g) and excluding secondary electron emission (
Figure 4h). Here, the electron heating rate in the sheath regions reduced again, and there was almost no heating in the plasma bulk, as seen in
Figure 4g,h; a pure
-mode was observed for both plots. This is clearly seen in
Figure 2a when averaged over the rf cycle. As seen from the difference plot shown in
Figure 4i, the electron heating in the sheath region was quantitatively more important for
than for
. There was clearly higher electron heating at sheath expansion when secondary electron emission was included. However, including ion-induced secondary electron emission from the electrodes decreased the overall electron power absorption.
Figure 5 shows the spatio-temporal behavior of the effective electron temperature. It shows the effective electron temperature as a function of position between the electrodes within one rf cycle, for the various combinations of pressures including and excluding secondary electron emission from the electrodes. At 10 mTorr, the effective electron temperature was high within the plasma bulk, and no important difference was observed in the spatio-temporal behavior of the effective electron temperature between
and
, as seen in
Figure 5a,b, respectively. A peak in the effective electron temperature was observed within the bulk region on the instantaneous anode side and agrees with the region of peak electron heating seen in
Figure 4a,b. The difference in the effective electron temperature calculated excluding and including the secondary electron emission is seen in
Figure 5c. We see that the effective electron temperature within the plasma bulk was slightly higher when secondary electron emission was included. The peaks in the effective electron temperature were higher when secondary electron emission was included. We also observe that the effective electron temperature had a peak within the plasma bulk at the instantaneous anode side and at the sheath expansion at 25 mTorr for both
and
, as seen in
Figure 5d,e, respectively. This is clearly seen in the difference plot in
Figure 5f. The electron effective temperature was higher for
than for
in the bulk region. In particular, the peak in the bulk region on the instantaneous anode side increased when secondary electron emission was included. At 50 mTorr, we observe a peak in the effective electron temperature during the sheath expansion. We also see that there was an increase in the effective electron temperature within the bulk when secondary electron emission was included. For
, the effective electron temperature in the bulk and in the sheath region was lower than when secondary electron emission was included, as seen in
Figure 5g,h. This behavior is clearly manifest in the difference plot shown in
Figure 5i. At all pressures, we found that including secondary electron emission in the discharge model increased the electron energy.