On the importance of excited state species in low pressure capacitively coupled plasma argon discharges

In the past three decades, first principles-based fully kinetic particle-in-cell Monte Carlo collision (PIC/MCC) simulations have been proven to be an important tool for the understanding of the physics of low pressure capacitive discharges. However, there is a long-standing issue that the plasma density determined by PIC/MCC simulations shows quantitative deviations from experimental measurements, even in argon discharges, indicating that certain physics may be missing in previous modeling of the low pressure radio frequency (rf) driven capacitive discharges. In this work, we report that the energetic electron-induced secondary electron emission (SEE) and excited state atoms play an important role in low pressure rf capacitive argon plasma discharges. The ion-induced secondary electrons are accelerated by the high sheath field to strike the opposite electrode and produce a considerable number of secondary electrons that lead to additional ionizing impacts and further increase of the plasma density. Importantly, the presence of excited state species even further enhances the plasma density via excited state neutral and resonant state photon-induced SEE on the electrode surface. The PIC/MCC simulation results show good agreement with the recent experimental measurements in the low pressure range (1–10 Pa) that is commonly used for etching in the semiconductor industry. At the highest pressure (20 Pa) and driving voltage amplitudes 250 and 350 V explored here, the plasma densities from PIC/MCC simulations considering excited state neutrals and resonant photon-induced SEE are quantitatively higher than observed in the experiments, requiring further investigation on high pressure discharges.


Introduction
Radio frequency (rf) capacitively coupled plasma (CCP) discharges have a wide range of applications, but in particular they are applied for material processing in the semiconductor industry. The CCPs are currently indispensable for etching and thin film deposition in integrated circuit fabrication. The discharge is created when a rf voltage is applied across a gasfilled region between two electrodes. The gas pressure can vary over a wide range and the discharge properties vary with the operating pressure. For example, at low operating pressure (a few Pa), the mean free path for both electrons and ions is larger than the electrode spacing, and the ion velocity is highly non-isotropic when ions arrive at the electrode surface after the acceleration across a high sheath field [1,2]. For the fabrication of solar panels and flat panel displays, which involves thin film deposition on a large area surfaces, the pressure is often higher, or in the intermediate pressure regime, on the order of 130 Pa [3][4][5][6][7].
Numerical simulations are of significant importance in understanding the fundamental mechanisms of capacitive discharges, especially, particle-in-cell simulations coupled with Monte-Carlo collisions (PIC/MCC), which self-consistently calculate electron/ion energy and velocity distribution functions, and have been proven to be a very powerful tool, applicable over a wide pressure range [8][9][10][11]. Discharges operated in noble gas, especially argon, are usually studied to gather understanding of the fundamental mechanisms in plasma discharges. In earlier studies, Roberto et al [12] and Lauro-Taroni et al [13] showed that the direct ionization from the ground state atom is the main ionization source for low pressure argon discharges, and the ionization from exited state species, like metastable pooling, and step-ionization processes, is negligible, by means of PIC/MCC simulations, but in the absence of secondary electron emission (SEE) from the surface. Therefore, in most of the particle-based modeling studies at the current time [14][15][16][17][18], the excited state atoms with a relatively low density compared to the feedstock gas are excluded [13,19]. Meanwhile, it is also known that with increasing pressure, the plasma gas composition may differ from the feedstock gas, and the reaction-generated excited state species, especially the metastable atoms, Ar m , which are long-lived, although having a lower density than the feedstock gas atoms, can alter the discharge properties through metastable pooling and stepionization. We know that while fluid modeling inevitably has certain assumptions on ion/electron energy distribution function and transport coefficients, particle-in-cell simulation tracing excited state neutrals as individual particles via a direct simulation of MCCs is computationally expensive.
In our recent works on capacitive argon discharges [20][21][22], the excited state neutrals are modeled as space-and time-evolving fluids incorporated with PIC/MCC simulations that treat charged particles as individual computer particles. This way, not only are the nonlocal dynamics of charged particles properly included as conventional particlebased modeling, but the effect of the excited state neutrals is also captured. It is found that the presence of metastable atoms enhances the plasma density by a factor of three at an intermediate pressure of 213 Pa, and the main ionization source alters from electron impact ionization of the ground state atom at low pressure (7 Pa) to metastable pooling and step-ionization at higher pressure (666-2000 Pa) [20]. The photon emission process of the Ar(4p) manifold, Ar(4p) → Ar m + hν, contributes a considerably large proportion of Ar m production at low pressure, while the electron impact excitation gradually dominates the Ar m production at high pressure (666-2000 Pa). The plasma density axial profile is also found to transition from parabolic shape at low pressure to a 'passive' flat bulk shape at high pressure due to the ionization occurring at the sheathplasma interface with/without excited state neutrals [7,22]. In absence of excited state neutrals, the ionization near the sheath-plasma interface is due to electron impact ionization; while, in presence of excited state species, the ionization near the sheath edge is mainly due to metastable pooling ionization because the metastable atoms Ar m are mainly located near the sheath edge [22].
In addition to the plasma processes in the bulk region, the surface processes of SEE have also attracted increasing attention in recent years [21,[23][24][25][26][27]. The simplest assumption for plasma surface processes is to set the ion-induced SEE coefficient to be a constant, often it is assumed to be 0.1, and the electron is simply assumed to be reflected with a probability of 0.2. Under such assumptions, for rf driven CCPs operated at low pressure, the secondary electrons are accelerated by the high sheath field to penetrate the bulk plasma region collisionlessly, and to become lost at the opposite electrode, even though the secondary electrons absorb power from the rf source. However, the effect of SEE on the plasma density is generally assumed to be less important at low pressure [1]. The γ-mode in which secondary electron-induced ionization breakdown within the sheath is dominant is known to appear at high pressure and high voltage [28].
In 1998, Gopinath et al [29] introduced a more realistic treatment for electron-electrode interaction in multipactor discharges, which are sustained within two parallel metal plates. In their approach the empirical Vaughan's formula describing electron-induced secondary electron yield (SEY) that depends on the impact energy and incident angle [30,31], is adapted. Different from the simple assumption of a constant electronreflection coefficient, the empirical Vaughan's model considers the incident electron energy and angle, allowing the electron-induced SEE yield to be above unity if the incident electron energy is above certain energy, i.e. tens of volts for most metal and dielectric materials. It implies that an ioninduced secondary electron with high energy, after the acceleration in high sheath field, is capable of producing more than one secondary electron at the opposite electrode if it can overcome the opposite sheath potential barrier and arrive at the electrode at a high energy. The newly ejected secondary electron moves toward the bulk plasma region under the action of the sheath field, that usually points toward the electrodes, and produce ionization. This regime was found by Horváth et al [24] to significantly enhance the plasma density in low pressure (0.5 Pa), high voltage (∼1000 V) rf driven argon CCPs.
The ion-induced SEE coefficient in reality relies on the ion impact energy and surface models. For a clean metal electrode surface, the SEY of argon ions is almost a constant 0.07, when the ion energy is below around 500 V, and slightly rises with increasing ion impact energy, while, for a dirty electrode surface, the argon ion-induced SEY is as low as 0.01 for ion energy below 10 V and increases up to 0.4 when the ion energy is 1000 V [32].
The discharge characteristics under realistic ion and fast atom-induced SEYs have been studied for single-frequency and multi-frequency low-pressure capacitive rf discharges driven by tailored voltage waveforms in argon [23,[25][26][27] and electronegative gas [33,34]. In our recent work, the excited state neutral-induced SEE was considered in intermediate pressure rf driven argon CCPs, and a γ-mode discharge is observed in presence of excited state neutral-induced SEE due to a large excited neutral-to-ion flux ratio of SEE [21]. In addition, the electrode material effects for SEE from SiO 2 , Cu, Si surfaces have also been examined by PIC/MCC simulations by several groups [16,[35][36][37].
Although fruitful knowledge on both bulk plasma and surface processes in low-pressure argon CCPs has been accumulated in the past years, there is still a long-standing issue on the quantitative validation between PIC/MCC simulations and experimental measurements. In the early studies, Vahedi et al [9] compared the modeling results from PIC/MCC simulations to laboratory measurements by Godyak et al [38] and showed that the electron energy probability functions from the PIC/MCC simulations are very similar to experimentally determined ones, over a gas pressure range from a few Pa to tens of Pa, but the plasma density from the simulations is lower than the experimentally determined ones by a factor of two. The deviations in the electron density were reevaluated and confirmed when revisiting the topic of electron power absorption in capacitive argon discharge by Lafleur et al [39]. In fact, lower plasma density in particle-based simulations is frequently observed when quantitatively comparing to experimental measurements, and they may differ from each other even by a few times [27,40,41], indicating that the physics considered in the PIC/MCC simulation models may be incomplete.
More recently, Schulenberg et al [42] conducted a multidiagnostic experimental validation of 1d3v PIC/MCC simulations for a low pressure rf driven capacitive argon discharge in typical discharge conditions, i.e. the operation gas pressure is in the range from 1 to 20 Pa, and the rf driving voltage amplitude is in the range from 150 to 350 V. It was found that a good agreement between simulation and experiments can be achieved only if the electron reflection coefficient is set to be an 'effective' constant up to 0.7. This 'effective' electron reflection coefficient was used by Derzsi et al [43], while studying a capacitive O 2 /Ne discharge.
In this work, we focus on the comparison of PIC/MCC simulation results against the recent experiments conducted by Schulenberg et al [42] and we will explore the physics that may be missing in modeling CCPs in the past by PIC/MCC simulations. The PIC/MCC code (oopd1) [20][21][22] that was strictly benchmarked and upgraded recently to be capable of considering the excited state neutrals-involving reaction processes such as metastable pooling ionization, step ionization, photon emission, production and loss of excited state neutrals, as well as SEE induced by electrons, ions, and neutrals (including excited atom states). We also explore the effect of resonant photon-induced SEE from the electrode surface in this work. The electron-induced SEE, modeled by the empirical Vaughan model is identified to be important, the excited state speciesrelated processes including the neutral-induced SEE on the electrode surface, and resonant photon-induced SEE, are also found to play a critical role in low pressure CCPs. The fluxes of different species like neutrals, ions, and photons flowing towards the electrodes are also analyzed in detail.
In section 2, we briefly introduce the PIC/MCC model. The modeling results and their comparison with recent experiments are shown in section 3, and finally, discussions and conclusions are drawn in section 4.
All the simulations conducted here are for capacitive argon discharges sustained between two planar stainless steel electrodes with the left hand electrode connected to a rf voltage source through a dc blocking capacitor and the right hand electrode grounded. The electrode spacing is kept at 4 cm, the driving frequency is 13.56 MHz, the feedstock gas pressure varies from 1 to 20 Pa, and the rf voltage amplitudes are 150, 250, and 350 V, respectively. Those conditions are the same as in the experiments of Schulenberg et al [42].
In the simulations, the charged particles (electrons and/or ions) are treated as individual particles, and a computer particle represents a cluster of 10 6 -10 12 real particles and moves in the direction normal to the electrode with its velocities updated in three directions using an explicit leap-frog integration scheme. When the system reaches a steady state, the particle number per discrete space cell width and per Debye length is much larger than unity to minimize the discrete particle noise, thus a large number of computer particles (∼10 5 ) are traced. The electric field is obtained by solving Poisson's equation where the net charge density is collected over all the charged particles by linearly weighting the charge within the cell into two neighboring discrete grid points. One thousand of discrete cells are used in the simulation with each cell width resolving the Debye length.
The collision dynamics of the charged particles is simulated by a null collision method [10]. For ions, the elastic and charge exchange collisions with neutrals are implemented by using the isotropic elastic scattering and backward scattering cross sections of Phelps [61]. The electron-neutral collisions in this work considers elastic scattering, electron-neutral excitation, and electron impact ionization. The argon atoms are excited to multiple levels of states, including the metastable-level atoms Ar m and the resonance radiative atoms Ar r of the 4s manifold, and the 4p-manifold (Ar(4p)). The corresponding energy thresholds and cross sections are given in detail in our previous works [20,21].
Since the electron elastic scattering influences the electron power absorption, and in turn affects the plasma density, here we introduce more details about the calculation of the elastic scattering angle. The elastic scattering between electrons and neutrals is implemented in the center-of-mass system, and the scattering angle χ follows the expression derived by Okhrimovskyy et al [62]: where ξ is the screening parameter and R is a random number (R ∈[0,1]). For a conventional screened Coulomb interaction, ξ = 4ϵ/(1 + 4ϵ) and ϵ = E inc /E 0 with E inc the incident electron energy in the center-of-mass system before collision, and E 0 = 27.21 eV, respectively. To make the elastic momentum transfer more consistent with experiments, the screening parameter ξ is also derived by Okhrimovskyy et al [62] from the normalized differential scattering cross section and expressed as a function of the ratio of the momentum transfer cross section and the total elastic scattering cross section recommended by Hayashi [63]. The comparison of different calculations of scattering angle, like the non-isotropic scattering angle by Vahedi et al [10], equation (1) and isotropic scattering based on the momentum transfer cross section, the screening parameter including Coulomb screening and Hayashi's screening parameter, and their effects on plasma density electron temperature, electron power deposition at low (6.7 Pa) and intermediate (213 Pa) pressure were discussed in more detail in an earlier work [22]. In this work, the screening parameter based on the experimental momentum transfer and total elastic scattering cross section recommended by Hayashi is adapted to calculate the scattering angle in equation (1). The algorithm for electron-neutral excitation and electron impact ionization in terms of energy partition and determination of the scattering angle of ejected electrons can be found elsewhere [20,53].
When an energetic electron arrives at the electrode surface, electron-induced SEE is considered and the emission coefficient is calculated by a modified empirical Vaughan's formula [24,[29][30][31]64]. The total SEE coefficient (γ e ) fits the experimental data provided by Baglin et al [65] and by Furman and Kirby et al [66,67] for normal incidence on stainless steel surface. The experimental data provided is for primary electron energy in the range 10-1000 eV [65], and 50-2000 eV [66,67], respectively. γ e based on Vaughan's model (Vau) used in this work and the experimental data (Exp a and b) [65][66][67] for the fitting are shown in figure 1. For low energy primary electrons incident on the surface, an elastic reflection process with the emission velocity specularly symmetric with respect to the surface normal vector is considered as done elsewhere [24,64,68]. For higher energy primary electrons, the true SEE is described by the conventional Vaughan's model [30,31], with 3% elastic reflection component, and 7% inelastic backscattered component [29]. The corresponding three components elastic reflection, inelastic backscattered and true secondary electron coefficients are plotted and labeled as η e , η i and η SEE , respectively, in figure 1. The parameters characterizing the conventional Vaughan's formula, elastic reflection and inelastic back-scattering that fit experimental measurements for stainless steel surface are summarized in table 1 [24].
The excited state atoms, the metastable-level Ar m (two levels treated as one), the radiative-level Ar r (two levels treated as one) and all the levels of the 4p-manifold Ar(4p) (taken as one species), are treated as time-and space-varying diffusing fluids. The feedstock gas argon atoms are assumed to be a fixed fluid that is spatially uniform within the two metal electrodes. The gas temperatures used in our simulations are the values experimentally measured by Schulenberg et al in [42] and increase from 300 to 330 K when the pressure varies in the range of 1-20 Pa. The effect of the gas temperature is comparable to the effect of the different plasma species induced SEE coefficients. For example, the plasma density for 150 V and The threshold for elastic reflection 0 eV ηe,max The maximum elastic reflection coefficient The energy corresponding to ηe,max 5.5 eV ∆e The control parameter for the decay of ηe 10 eV re Portion of elastically reflected high-energy electrons 3% [29] Portion of inelastically back-scattered high-energy electrons 7% [29] 20 Pa is around 10% higher at room temperature than that at 330 K, however, the excited-state neutral and photon induced SEE increases the plasma density by around 22% and 54%, respectively (see more detailed discussion later in figure 3(b)).
The excited state atoms are assumed to diffuse only against the background gas and not against each other as the excited state atom densities are much lower than the feedstock gas density. The governing diffusion equation for fluid species also includes the source and loss terms originated from the impact reactions between fluid and charged particles, as well as fluidfluid reactions. When an electron or ion is generated via fluidfluid reactions, the particle is added as an individual particle. On the electrode surface, the recombination coefficient is set to be 0.5, a value estimated by Stewart [69], implying that the excited state atoms are 50% quenched and 50% reflected on the electrode surface. The excited state atoms also lead to SEE on the surface with a coefficient of 0.21 for Ar m and Ar r , and 0.27 for Ar(4p) [70]. The resonance radiation of Ar r is partially imprisoned at low pressure, and the fraction of the radiation escaping depends on the specific gas pressure and electrode spacing. The Walsh model [71] is used to calculate the escape factor g. Typically, higher pressure leads to a lower value of g, and a smaller fraction of the resonant photons escape from the plasma to the surface. The escaped photons are traced and their impact on the electrode surface to produce SEE is calculated. For a photon energy of 11.62 eV, the resonant photon-induced secondary electron emission coefficient is set to 0.075, taken from the measurements of Feuerbacher and Fitton [72] for stainless steel.

Results
In this section, we present the PIC/MCC simulation results and the comparison to experimental measurements with emphasis on the plasma density and the effects of different surface processes. For comparison, the discharge conditions set in the PIC/MCC simulations are the same as in the experiments of Schulenberg et al [42] as follows: the electrode spacing is fixed at 4 cm, and the driving frequency is 13.56 MHz, the feedstock gas pressure is in the range of 1-20 Pa, and rf voltage amplitude is 150, 250, and 350 V. The surface processes that we explore focus on the electron reflection (constant coefficient), electron-induced real SEE (incident energy and angle-dependent SEY using the modified empirical Vaughan's formula), excited state species-induced SEE, and resonant photon-induced electron emission. These processes are considered incrementally in the PIC/MCC simulations (see the more detailed definitions in the following). The case of electron reflection with a constant coefficient 0.2 is the same as in the simulations of Schulenberg et al [42] except that more reactions involving excited state species like metastable pooling and step-wise ionization are considered. However, in this work we explore the more realistic surface process of SEE from energetic electrons, excited state species, and photons rather than using an 'effective' electron reflection coefficient.
In our simulations, the ion-induced SEE yield is a function of the ion impact energy and we assume dirty electrodes [21,32,73] for all the cases. Generally, a clean electrode means that the metal surface is flashed at a very high temperature for a long time and the properties may change after a duration of discharge operation. Therefore, the surface condition of a 'dirty' electrode is used to calculate the ion-induced SEE yield [32,73] in this work. For other parameters involving surface processes, even though these SEE coefficients may be imperfectly known, they are chosen as reasonable as possible based on the published literature.
The PIC/MCC simulations were conducted for four cases corresponding to four different surface models:   For a comparison of the cases explored an overview is given in table 2. Note that the energetic neutral state atoms are tracked as individual particles only when the energy is above 32 eV that is the threshold for fast neutral-induced SEE for a dirty metal surface. However, since the energetic neutral state atom-induced SEY is energy dependent and much lower than the ion-induced SEY in the energy range of interest (0-160 eV) where 160 eV is almost the highest ion energy at 350 V [42], the role of the fast ground state atoms induced SEY is indeed negligible, which is supported by simulations including and excluding energetic ground state neutrals.
For conciseness, we will refer to different electron emission processes and cases via the SEE coefficient, for example, if case I is discussed, then we will refer to 'γ e = 0.2; γ exc = 0; γ ph = 0', directly. It is worth noting that both the metastable state atoms Ar m and the resonant state atoms Ar r (also Ar(4p)) flowing toward the electrodes can produce SEE. Here we use γ exc to represent the excited state neutral-induced secondary electron coefficient. Note that the flux of metastable atoms is much higher than the flux of the other excited state species by two or three order of magnitudes in the pressure range investigated here. Figure 2 shows the electron density versus gas pressure for the three driving voltage amplitudes. The figures show the experimentally determined value and the PIC/MCC simulation results for cases I-IV. The black dotted lines with error bar in figures 2(a)-(c) show the experimental measurements of the plasma density at the discharge center by Schulenberg et al [42], versus pressure (1-20 Pa) at driving voltage amplitudes of 150, 250, and 350 V, respectively. We can see that the electron density n e monotonically increases with increased pressure for all driving voltage amplitudes except that a slight decrease is seen at 20 Pa for 250, and 350 V. Similar to the observations that were made by Schulenberg et al [42], the plasma density from the PIC/MCC simulations for the case 'γ e = 0.2, γ exc = 0, γ ph = 0' is significantly lower than the experimental measurements at the low pressure of 1 Pa. However, when the more realistic electron-induced SEE process described by modified Vaughan's model is considered, as shown by case II 'γ e = Vau, γ exc = 0, γ ph = 0' in figure 2(a), the plasma density is significantly enhanced and increases from 0.24 × 10 15 to 0.62 × 10 15 m −3 for 150 V. Incrementally, the presence of γ exc , as shown by case III: 'γ e = Vau, γ exc = 0.21, γ ph = 0', further increases the plasma density to 0.8 × 10 15 m −3 , and γ ph = 0.075 further increases the plasma density around 7.5% (case IV). This implies that the addition of excited state species, that produces excited state neutral and resonant photon impact on the surface and creates secondary electrons, enhances the plasma density by 35% in total.
At higher pressure, in the range 5-10 Pa, the plasma densities with or without the consideration of Vaughan's model are almost the same. However, the excited state neutral and resonant photon-induced secondary electrons still enhance the plasma density. To make the plasma density enhancement from different surface processes more visible quantitatively, the plasma density profiles at 1 and 20 Pa for 150 V are shown in figures 3(a) and (b), respectively. When the driving voltage amplitude is 250 or 350 V, while the pressure is kept at 1 Pa, the plasma density is also enhanced by a factor of two. Similar to the case of 150 V, as shown in figures 2(b) and (c), the change of γ e from 0.2 to Vaughan's model has minimal effects on the plasma density at higher pressures (5-20 Pa) for 250 and 350 V. The contribution of SEE from different surface processes depends on the electron dynamics present in different parameter ranges that will be further discussed in the following.
The external conditions controlling the discharge properties explored here in the PIC/MCC simulations are the SEE, mainly determined by the incident electron flux, energy, and angle, incident ion flux and energy, neutral (including excited states) flux and resonant photon flux. Firstly, we analyze the origin of the enhanced plasma density for Vaughan's model compared to a constant electron reflection coefficient. Comparing these two surface models, the only difference is the dependence of the SEY on the incident primary electron properties, i.e. the SEY in Vaughan's model relies on the impacting electron energy and incident angle. Although an oblique incidence produces a slightly higher SEY compared to a normal incidence at the same impact energy, the SEY strongly depends on impact energy, therefore, we plot the electron energy distribution, which mainly responds to the spatiotemporal electric field for a collisionless low pressure discharge (1 Pa).  figure 4(a)), the E x in figure 4(a) is limited within the range of −0.3 × 10 4 -0.3 × 10 4 V m −1 , and the temporally-evolving boundary of the gray area represents the sheath edge. We can see in figure 4(b) that the primary electrons with energy above the sheath potential arrive on the right electrode (x = 4 cm) and form a pattern during the sheath collapse. The color represents the fraction of electrons within the specific energy range in the colorbar, and all primary electrons incident on the electrode having energy below 25 eV will contribute a small number of emitted secondary electrons.
The ions respond to the time-averaged sheath field due to their larger mass than electrons, and flow towards the electrodes at an almost constant flux, thus, the ion-induced secondary electrons are always present over the whole rf period regardless of the sheath collapse or expansion. As a result, the high sheath field (see arrow A in the figure 4(a)) in the gray area accelerates these secondary electrons to a high energy and consequently they impact the opposite (right) electrode (x = 4 cm) in the form of an electron beam as shown in figure 4(c). The energy of the secondary electron beam corresponds to the rf sheath potential. As the potential drop across the bulk plasma is small, the sheath potential drop when the sheath approaches the maximum width is almost equal to the driving voltage amplitude, 150 V, and the maximum energy of the electron beam is around 150 V. The plasma potential with respect to the two electrodes is equal and has no net contribution to the secondary electron energy gain. Comparing the starting point of the secondary electron beam (t/T ≈ 0.5) in  figure 4(a), a phase delay can be found. The reason is that the accelerated secondary electrons emitted from the left electrode (x = 0 cm) at the initial phase duration of sheath expansion are reflected by the opposite sheath potential, and arrive at the opposite electrode (x = 4 cm). Note that the travel time of secondary electrons across the bulk region is much smaller than the rf period, and consequently it contributes very little to the time delay. The energetic secondary electrons striking the surface produce a large number of new secondary electrons, as the SEY is above unity for Vaughan's model (also called δ electrons in some literature [14,16,37]). The new secondary electrons are accelerated by the ambipolar electric field (see arrow B in figure 4(a)) to inject into the bulk plasma. Thus these new secondary electrons will have a similar properties as the primary electrons and are accelerated and decelerated back and forth by the sheath expansion and collapse, leading to additional ionization impacts in the bulk region. A considerable number of secondary electrons having the same properties as the primary electrons are accumulated in the simulation to form a pattern in figure 4(c) similar to that of primary electrons in figure 4(b). Ultimately, the plasma density is significantly enhanced by energetic secondary electrons that produce a large number of secondary electrons due to the surface process described by Vaughan's model.
At low pressure (1 Pa), the electron mean free path is much larger than the electrode spacing, and the secondary electrons can penetrate the bulk region without collisions. When the pressure is increased to 20 Pa, the electron mean path is significantly shortened, and the electrons lose their energy via inelastic excitation and ionization impacts, with the latter producing more electrons directly, enhancing the plasma density. As can be seen in figure 4(d), the behavior of the electron beam generation disappears for secondary electrons bombarding the right electrode. Therefore, the plasma density is almost the same for constant electron reflection coefficient and Vaughan's model at higher pressure (5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15)(16)(17)(18)(19)(20). In figures 4(e) and (f), we show more information on the ionization and secondary electron energy distribution at the midplane (x = 2 cm) for discharge conditions the same as of figures 4(b) and (c), respectively. As expected, the primary electron energy distribution oscillates continuously in response to the oscillating sheath.
The secondary electron energy distribution shows two beams due to the separate acceleration of the secondary electrons by the alternating high sheath fields adjacent to the left and right electrode. The two secondary electron beams at the mid-plane are almost symmetric for the first and second half rf period due to the symmetry of the discharge. The highest energy of the secondary electrons in figure 4(f) is above 150 V because of the presence of the plasma potential. The electron beams were also observed in previous PIC/MCC simulations where only ion-induced secondary electrons were considered and where the ionization solely produced by those ballistic electrons was analyzed [74].
In addition to electron and ion-induced SEE discussed above, the excited state neutrals and resonant photons also lead to SEE with the former and the latter controlled by the metastable atom Ar m and resonant photon flux flowing towards the electrodes, respectively. The resonant photon is from the escaped radiation of the resonant level atoms Ar r . As the Ar m and Ar r diffusion is not affected by the oscillating electric  field, their density profiles are less significantly varying over the whole rf period. This implies that the secondary electrons are also emitted continuously from both the left and right electrodes over time, similar to the dynamics of ion-induced secondary electrons.
Quantitatively comparing the time-averaged fluxes of ions, excited state neutrals, and photons can help understand the importance of different species, as shown in figures 5(a)-(c) for 150, 250, 350 V, respectively. The flux of excited state neutrals is defined as Γ n,exc = (1 − 0.5γ rec ) −1 n exc v n /4 with the mean thermal velocity v n = √ 8k B T/πm and recombination coefficient γ rec = 0.5 and n exc is the density of excited state neutrals on the electrode. We can see that the excited state neutral flux, Γ n,exc , increases at first and then decreases with increasing the feedstock gas pressure in the range of interest, and maximum Γ n,exc appears at 2 Pa. There are two competitive factors affecting Γ n,exc : (i) the excited state neutral densities; and (ii) the collisions with the feedstock gas. The higher excited state neutral densities at higher pressure tend to result in higher flux Γ n,exc . However, the more frequent collisions between excited state neutrals and feedstock gas reduce the flux toward the electrodes. Factor (i) is dominant over factor (ii) when the pressure is between 1 and 2 Pa, and factor (ii) dominates when the pressure is higher than 2 Pa. As γ exc is almost three times larger than γ ph , the back flowing secondary electron flux, γ exc Γ n,exc , from neutral atoms in excited state is larger than the secondary electron flux that is due to photons bombarding the electrodes γ ph Γ ph , leading to a higher density gain at 1 Pa from excited state neutrals than photons (see figure 3(a)).
The ion and resonant photon fluxes increase with increased gas pressure, and the density gain from photon-induced secondary electrons becomes greater than from excited state neutrals (see figure 3(b)). The resonant photon flux incident on the electrode is proportional to the product of the Ar r density and the escape factor g, i.e. fraction of the escaped photons after radiation trapping. We can see that the Ar r densities in figure 6(a) increase significantly at higher pressures due to the generation from the chemical reaction processes in the plasma bulk. The escape factor g shown in figure 6(b) decreases from 3 × 10 −3 to 0.8 × 10 −3 with pressure for a given electrode spacing. It is worth noting that, indeed, the Ar r loss in the plasma is affected by the radiation rate (linearly proportional to the escape factor g). Finally, the nonlinear combined effect results in an increased photon flux flowing toward the surface. This resonance photon-to-ion flux ratio is at the upper limit of what has been measured experimentally for inductively coupled plasmas [75], and so our results would be the upper limit of what influence γ ph might have.

Discussions and concluding remarks
Kinetic PIC/MCC simulations have been an important tool for the understanding of the physics of low pressure capacitive discharges in the past three decades. However, the plasma density from PIC/MCC simulations commonly shows quantitative deviations from the experimental measurements, even in argon discharges. Generally, the observed plasma densities from particle-based simulations are lower than what is observed experimentally, when experimental measurements are compared to simulation results [27,40,41], lower even by a few times [27], indicating that certain physics may be missing in previous modeling of low pressure rf capacitive discharges. More recently, the PIC/MCC simulations were validated against experimental diagnostics [42], and it was found that an 'effective' electron reflection coefficient of 0.7 is required in simulations to match the experimental measurements. However, the underlying physics supporting the coefficient of 0.7 is not fully understood.
In this work, we explored the potential surface processes that may result in deviations between PIC/MCC simulations and the recent experiments [42]. Using the same external conditions as the experiments, argon discharges are sustained within two stainless steel electrodes with a spacing of 4 cm. The feedstock gas pressure is varied in the range of 1-20 Pa, and the driving voltage amplitudes are 150, 250, and 350 V at a driving frequency of 13.56 MHz. It is found that both the energetic electron-induced SEE and excited state atoms play important roles in low pressure (1 Pa) rf capacitive argon plasmas. The ion-induced secondary electrons are accelerated by the high sheath field to bombard the opposite electrode, producing a considerable number of secondary electrons. These energetic electron-induced secondary electrons are accelerated by the ambipolar electric field within the collapsed sheath toward the bulk plasma and act similarly to the primary electrons, i.e. these electrons are accelerated and decelerated back and forth by the oscillating sheaths over time, leading to additional ionization impacts that finally increase the plasma density.
Importantly, the presence of excited state species further enhances the plasma density via SEE from the excited state neutrals and resonant state photons (originating from the escaped radiation of resonant level atoms Ar r ) bombarding the electrode surface. With the consideration of SEE from energetic electrons, excited state neutrals and resonant photons, the PIC/MCC simulation results agree with the recent experimental measurements at low pressure (1-10 Pa).
At higher pressures, above 5 Pa, the plasma density from the electron surface process described by Vaughan's model is almost the same as that considering a constant electron reflection coefficient of 0.2, because of the absence of energetic electron beams at higher pressure where the electron mean free path is shorter. The excited state neutral and resonant photoninduced secondary electrons still can enhance the plasma density through direct ionization impacts within the plasma bulk region. To further understand the relative importance of excited state neutrals and resonant photons, we examined their flux flowing toward the electrode surface. The excited state neutral-induced secondary electron flux dominates over the photons at 1 Pa. With increasing pressure in the range of interest (5-20 Pa), the resonant photons play a more important role in SEE than excited state neutrals.
We also note that the plasma density from PIC/MCC simulations with both excited state neutral and photon-induced SEE included in the model are quantitatively higher than determined in the experiments at the highest pressure of 20 Pa for 250 and 350 V. There are a few possible reasons for the deviation: (i) the plasma discharge is more local radially at high pressure and a density peak usually appears above the electrode edge at high pressure [40,76]; radial geometry effects are excluded in the current one-dimensional simulations; (ii) outgassing processes from the electrode and glass wall may exist under ion impacts at higher pressures due to contamination with air or water vapor in the experiments; these contaminants may affect the discharge characteristics; (iii) the photon-induced SEE yield may be different for different treatments of metal electrode; for example, a reduced SEY from experiments is reported for heat-treated metal surface (Au, Cu and Pt) compared to untreated metal surface as discussed by Phelps et al [32], and the coefficient 0.075 is an absolute upper bound that can be expected from nonoxidized, hydro-carbon free, defect free surfaces. In day-today usage, γ ph is likely smaller than 0.075. In addition, it was stated by Schulenberg et al [42] that an additional discharge around the Langmuir probe started to appear at pressure above 20 Pa, which may influence the experimental results at that pressure. Therefore, in the future, it is worth further investigating the possible missing physics in PIC/MCC modeling, and further developing improved measurement tools, for quantitative comparisons between simulations and experiments in high pressure discharges.