Electron power absorption dynamics in a low pressure radio frequency driven capacitively coupled discharge in oxygen
aa r X i v : . [ phy s i c s . p l a s m - ph ] F e b Electron power absorption dynamics in a low pressure radio frequency drivencapacitively coupled discharge in oxygen
A. Proto and J. T. Gudmundsson
1, 2, ∗ Science Institute, University of Iceland, Dunhaga 3, IS-107 Reykjavik, Iceland Department of Space and Plasma Physics, School of Electrical Engineering and Computer Science,KTH Royal Institute of Technology, SE-100 44, Stockholm, Sweden (Dated: February 8, 2021)We use the one-dimensional object-oriented particle-in-cell Monte Carlo collision code oopd1 toexplore the properties and the origins of both the electric field and electron power absorption withinthe plasma bulk for a capacitively coupled oxygen discharge operated at 10 and 100 mTorr for45 mm of gap distance. The properties of the electric field at three different time slices as well astime averaged have been explored considering the moments of the Boltzmann equation. The electronpower absorption is distinctly different at these operating pressures. The most relevant contributionsto the electric field at different time steps come from the pressure terms, the ambipolar and theelectron temperature gradient terms, along with the ohmic term. The same applies for the electronpower absorption. At both 10 and 100 mTorr the relative ohmic contribution to the electron powerabsorption remains roughly the same, while the ambipolar term contributes to power absorptionand the temperature gradient term to electron cooling at 100 mTorr, and the opposite applies at10 mTorr. At 100 mTorr the discharge is weekly electronegative, and electron power absorption ismainly due to sheath expansion while at 10 mTorr it is strongly electronegative, and the electronpower absorption occurs mainly within the electronegative core and drift-ambipolar mode dominates.The agreement between the calculated values and the simulations is good for both the electric fieldand the electron power absorption within the plasma bulk and in the collapsed sheath region for allthe cases considered.
I. INTRODUCTION
The low pressure radio frequency (rf) driven capaci-tively coupled discharge has been applied in integratedcircuit manufacturing for a few decades. The capacitivelycoupled discharge consists of two parallel electrodes, typ-ically with a radius of few tens of cm, separated by afew cm and driven by a power generator. These dis-charges have been explored extensively over the past fewdecades. The power transfer mechanism, which is com-monly referred to as ’electron heating’ or ’electron powerabsorption’ in the literature [1], is still a topic ratherpoorly understood. Although the electron power absorp-tion mechanism is a topic widely studied and discussedover the past decades, a fully consistent and generalmathematical-physical explanation of the several phys-ical mechanisms involved in the power transfer mecha-nism are still lacking. This is in particular true for theelectronegative capacitively coupled discharge.It is widely accepted that the electron heating can bedivided into two components: the ohmic heating (col-lisional) and the stochastic heating (collisionless) whileseveral operating modes have been identified in the ca-pacitively coupled discharge including the stochastic elec-tron heating due to the sheath motion ( α -mode) [2], sec-ondary electron emission due to ion and neutral bom-bardment of the electrodes [3], the drift ambipolar DA-mode [4], non linear electron resonance heating (NERH) ∗ [email protected] [5–9], electron bounce resonance effect [10, 11] and thegeneration of series resonance oscillations [5, 7]. In astrongly electronegative discharge the electrical conduc-tivity tends to be low and, due to large ion inertia highelectric field is induced within the plasma bulk (elec-tronegative core). Furthermore, ambipolar fields appearnear the sheath edges.The particle-in-cell (PIC) method, when combinedwith Monte Carlo (MC) treatment of collision processes,is the most frequently used numerical approach to inves-tigate the properties and the operating modes of low pres-sure capacitively-coupled-discharges. The combination ofthe particle-in-cell (PIC) method and Monte Carlo colli-sion (MCC) treatment of collision processes is commonlyreferred to as the PIC/MCC method. The PIC/MCCmethod is a self-consistent kinetic approach that has be-come a predominant numerical approach to investigatethe properties of the low pressure capacitively-coupleddischarge.The one-dimensional-object-oriented plasma deviceone ( oopd1 ) code allows having the simulated particlesof different weights so that in principle both chargedand neutral particles can be tracked during the simu-lation. Earlier we benchmarked the basic reaction set forthe oxygen discharge in the oopd1 code to the xpdpl code [12]. In recent years the oxygen reaction set inthe oopd1 code has been improved significantly [12–14].The oopd1 code has been applied to explore the electronpower absorption in the capacitively coupled oxygen dis-charge while varying the various external parameters andoperating conditions such as discharge pressure [14–16],driving voltage amplitude [17], driving frequency [18], thesecondary electron emission [14, 19], the surface quench-ing of the metastable states [20] and the electrode gapdistance [21].During the past decades several attempts to describecorrectly the behaviour of the electron heating using theBoltzmann equation have been made. Surendra andDalvie [22] were the first to set up a mathematical modelto describe the electron power absorbed using the Boltz-mann equation for both electrons and ions using the PICresults as input. In the years that followed several au-thors used the formulation set by Surendra and Dalvie[22] to develop similar models inspired by the their re-sults [23–29]. Among these Brinkmann [27] derived aunified description of electron power absorption in ca-pacitively coupled discharges using a mathematical for-mulation where the electron density profile has been ap-proximated by a smooth step function, finding that thetotal time averaged electron power absorption is the sumof four terms, each one corresponding to one of the heat-ing mechanism knows from separate previous theories,i.e. NERH, stochastic heating (hard wall model), ambipo-lar/pressure heating and ohmic heating. Brinkmann alsodemonstrated that a time dependent temperature is nec-essary to obtain a non zero time averaged electron powerabsorption. More recently, Schulze et al used a simplifiedmoment analysis of the Boltzmann equation (the Boltz-mann term analysis) where the electron temperature gra-dient was both neglected and considered [30, 31] in orderto describe both the electric field and the electron powerabsorbed in an electropositive low pressure capacitivelycoupled argon discharge. They found that the time av-eraged ambipolar electron power absorption completelyvanishes for a temporally independent electron tempera-ture. This approach has recently been applied to explorethe electron power absorption mechanisms in a capaci-tively coupled oxygen discharge by Vass et al. [32]. Usingthe Boltzmann term analysis, they found that the ohmiccontribution to the electron power absorption is small atdifferent time steps at low pressure while it becomes im-portant at higher pressures. Finally, they observed thatat low pressure the space-time averaged electron powerabsorbed was entirely given by the ohmic term and thatthe pressure term contribution increases as the pressureis increased. Here we use the Boltzmann term analysis toinvestigate the origins of both the electric field and thepower absorbed by the electrons at different time stepswithin both the plasma bulk and the sheath collapsed re-gion and the related time averaged quantities within theplasma bulk in a capacitively coupled oxygen dischargeat 10 and 100 mTorr for 45 mm of gap distance. The elec-tron power absorption in the capacitively coupled oxygendischarge is distinctly different when the discharge is op-erated at 100 mTorr than when it is operated at 10 mTorr[14, 15, 21]. When operating at 100 mTorr with a gapsize of 45 mm the discharge is operated in α -mode andstochastic electron heating dominates while at 10 mTorrthe discharge is more electronegative and the DA-modedominates. The main task of the current work is to per- form a Boltzmann term analysis of a capacitively cou-pled oxygen discharge, in order to shed light on the theunderlying physical mechanism behind the electric fieldand the electron power absorbed, to gain understandingof the electron power absorption in an capacitively cou-pled oxygen discharge operated at different pressures (10and 100 mTorr), that represent hybrid DA- α mode andpure α mode, respectively. We will follow the frameworkof the Boltzmann term analysis given in the recent workof Schulze et al. [31] with some modifications, since thephysical conditions and gas considered are different. Thecurrent work is structured as follows. In Section II wegive a brief overview of the simulation set up. In SectionIII we show the spatio-temporal profiles of both the totalcharge density and the quasineutrality deviation alongwith both the total charge density and the density pro-files for all the species involved at different time steps andtime averaged. Section IV discusses the model. In sub-section IV A a simple fluid model based on the Schulzeet al. [31] work is employed to explore the behaviour ofboth the electric field and the electron power absorbed atboth 100 and 10 mTorr. The results from both the simu-lations and the calculations at both 10 and 100 mTorr arediscussed and compared in Section V. Finally, Section VIsummarizes our findings. II. THE SIMULATION
The one-dimensional (1d-3v) object-oriented particle-in-cell Monte Carlo collision (PIC/MCC) code oopd1 [12]is applied to a capacitively coupled oxygen discharge. In1d-3v PIC codes the model system has one spatial dimen-sion and three velocity components. The oopd1 code, likethe well known xpdp1 code, is a general plasma devicesimulation tool capable of simulating various types of de-vices, where the plasma is the main actor, such as particlebeams, electrical breakdown, particle accelerators as wellas processing discharges [12]. The oxygen reaction setincluded in the oopd1 code is rather extensive and ninedifferent species are considered: electrons, the groundstate neutrals O( P) and O (X Σ − g ), the negative ionsO − , the positive ions O + and O +2 , and the metastablesO( D), O (a ∆ g ) and O (b Σ +g ). The basic reaction setincluded O (X Σ − g ), O +2 and O − . In our earlier work weadded oxygen atoms in the ground state O( P) and ionsof the oxygen atom O + to the oopd1 code [12]. In a laterwork the singlet metastable molecule O (a ∆ g ) and themetastable oxygen atom O( D) were added [13], as wellas the singlet metastable molecule O (b Σ +g ) [14]. Thefull oxygen reaction set together with the cross sectionsused have been discussed in our earlier works and willnot be repeated here [12–14, 21]. We assume a geomet-rically symmetric capacitively coupled discharge whereone of the electrodes is driven by an rf voltage at a singlefrequency V ( t ) = V sin (2 πf t ) (1)while the other electrode is grounded. Here, V is thevoltage amplitude, f the driving frequency and t thetime. For this current study we assume the discharge tobe operated at the pressure of 10 mTorr and 100 mTorrwith voltage amplitude V =400 V with an electrode sep-aration of 4.5 cm. A capacitor of 0.1 µ F is connected inseries with the voltage source. The electrode diameterand the driving frequency are assumed to be 10.25 cmand 13.56 MHz respectively. These are the same param-eters as assumed in our previous work [21]. The timestep ∆ t and the grid spacing ∆ x are set to resolve theelectron plasma frequency and the electron Debye lengthof the low energy electrons respectively, according to ω pe ∆ t < .
2, where ω pe is the electron plasma frequencyand the simulation grid consists of 1000 equal cells. Theelectron time step is 3 . × − s. The simulation wasrun for 5 . × time steps, which corresponds to 2750 rfcycles as it takes roughly 1700 rf cycles to reach equilib-rium for all particles. Time averaged plasma parametersshown, such as the densities, the electron power absorp-tion, and the effective electron temperature, are averagedover 1000 rf cycles. All particle interactions are treatedby the Monte Carlo method with a null-collision scheme[33]. For the heavy particles, we apply sub-cycling, wherethe heavy particles are advanced every 16 electron timesteps [34] and an initial parabolic density profile has beenassumed [34].The kinetics of the charged particles (electrons, O +2 ions, O + ions and O − ions) was followed for all energies.Since the neutral gas density is much higher than thedensities of charged species, the neutral species at ther-mal energies (below a certain cut-off energy) are treatedas a background with fixed density and temperature andmaintained uniformly in space. The main challenge whenPIC/MCC simulations are applied to simulate moleculargases has to due with the timescale difference betweenthe processes of dissociation and the processes involv-ing charged particles. Therefore, a global model [35] isused beforehand to determine the partial pressure of thevarious neutrals created in the discharge as discussed inProto and Gudmundsson [20], i.e. the ground state neu-tral atoms O( P) and the metastables O( D), O (a ∆ g )and O (b Σ + g ) under certain control parameters includ-ing the discharge pressure, the absorbed power and thegap separation between the two electrodes, etc. The ab-sorbed power determined by the PIC/MCC simulationis used as an input parameter in the global model cal-culations, iteratively. The partial pressure of the atomsand metastable species obtained from the global modelcalculation is then used as the partial pressure of thesespecies in the neutral background gas in the simulation.Note that, a global model is mainly developed to model ahigh density low pressure discharges such as inductivelycoupled discharges, rather than capacitively coupled dis- charges and the proportion of the power absorbed by theelectrons in the former is much larger than in the latter.Therefore, the global model may overestimate the atomand metastable density within the discharge, especiallywhen operating at low pressure. The fractional densitiesfor the neutrals O (X Σ − g ), O (a ∆ g ), O (b Σ g ), O( P),O( D), estimated using the global model calculations, arelisted in Table I. These values have been used as inputfor the PIC/MCC simulation as the partial pressures ofthe neutral background gas. These neutral backgroundspecies are assumed to have a Maxwellian velocity dis-tribution at the gas temperature (here T n =26 meV).The kinetics of the neutrals are followed when their en-ergy exceeds a preset energy threshold value. The energythreshold values used here for the various neutral speciesare listed in Table II. The thresholds were chosen in or-der to keep the number of simulated particles within asuitable range, typically 10 − particles. Particleswith energy below this threshold energy are assumed tobelong to the neutral background.Note that the background neutrals are assumed to beuniform within the discharge. However, we are awarethat the electrode surfaces have a significant influenceon the neutral density profiles. The density profiles forfast neutrals indicate that the oxygen atom density de-creases and the molecular metastable density increasesin the electrode vicinity [14]. As an oxygen atom O( P)hits the electrode, it is assumed that half of the atomsare reflected as O( P) at room temperature and theother half recombines to from the ground state oxygenmolecule O (X Σ − g ) at room temperature. Similarly, asa metastable oxygen atom O( D) hits the electrode, halfof the atoms are quenched to form O( P) and the otherhalf, is assumed to recombine to form the ground stateoxygen molecule O (X Σ − g ) at room temperature. Thesurface quenching coefficients for the singlet metastablemolecules on the electrode surfaces are assumed to havea value of γ wqa = 0 .
007 and γ wqb = 0 . (a ∆ g )and O (b Σ g ), respectively. The influence of the surfacequenching coefficients of the singlet metastable moleculeon the electron heating mechanism has been explored indetail in an earlier work [20], where it has been demon-strated that the influence of γ wqa on the overall dischargeproperties can be rather significant. The surface quench-ing and recombination coefficients used in this currentwork are listed in Table II. Note that the oxygen reac-tion set used in this current study is signficantly moreextensive than the one used by Vass et al. [32] where theonly metastable state included is the singlet metastablemolecule O ( ∆ g ). To estimate the density of the sin-glet metastable molecule O ( ∆ g ) they assume a homo-geneous spatial density where a balance between elec-tron impact excitation and quenching at the electrodes[36, 37].Secondary electron emission and electron reflection hasbeen incorporated into the discharge model [14]. Theenergy dependent secondary electron emission yield fora dirty surface has been employed [14] and the electronsare assumed to be elastically reflected from the electrodesindependently of their energy and their angle of incidencewith a probability of 0.2 [19]. The secondary electronemission due to the electron impact on the electrodeshas been neglected here as in all our previous works onthe oxygen discharge. TABLE I: The relative partial pressures of the thermalneutrals at 10 and 100 mTorr calculated by a global (volumeaveraged) model for 1.8 W power.O (X Σ − g ) O (a ∆ g ) O (b Σ g ) O( P)10 mTorr 0 . . . . . . . . III. RESULTS FROM THE SIMULATION
Through recent PIC/MCC simulations of a capaci-tively coupled oxygen discharge it has been demonstratedthat the singlet metastable molecular states have a sig-nificant influence on the electron power absorption mech-anism [13–16] as well as the ion energy distribution [38].At low (high) pressure and high (low) electronegativity,i.e. 10 mTorr (50 – 500 mTorr), the electron power ab-sorption is mainly located within the plasma bulk (thesheath regions) [14, 15]. Furthermore, when operatingat low pressure, the time averaged electron power ab-sorption within the discharge is due to a hybrid drift-ambipolar mode, (DA-mode) and α -mode, and while op-erating at higher pressures, the electron power absorptionis due to stochastic heating and the discharge is oper-ated in a pure α -mode [17, 18]. It has also been demon-strated that detachment by singlet molecular metastablestates is the process that has the most influence on theelectron power absorption process in the higher pressureregime, while it has almost negligible influence at lowerpressures [14–16]. All the quantities returned by the sim-ulations and involved in the calculations for both theelectric field and the electron power absorption in thefollowing sections are arrays extended along the x − axis,i.e. the discharge gap length. In particular, every singlecomponent of the electron temperature has been calcu-lated as T e ,ii = e E e ,ii − m e e u ,i where E e ,ii and u ,ii , with i = x, y, z , are the mean electron energy density and theelectron mean velocity respectively. Since E e ,i = m e h v ,i i by definition, the expression for the electron temperaturegiven above is the same as the one shown by Wilczek et al.[1], when the particle mean velocity is not negligible.Figure 1 shows the density profiles for O +2 ions, O + ions, O − ions and electrons at 100 and 10 mTorr. At100 mTorr the center electronegativity is 3.55 and at 10mTorr it is 93.64. The electronegative discharge consistsof an electronegative core connected to electropositiveedge plasma regions [39, 40]. At 100 mTorr (Figure 1 (a))both O − ion and O +2 ion density profiles have a similar shape within the bulk region while the O − ion densitydecreases more steeply than the O +2 ion density profileapproaching the sheath edges. We also see that the elec-tron density profile is somewhat lower than both the O − ion and the O +2 ion density profiles within the bulk re-gion while it decreases sharply within both the sheathregions. At 10 mTorr (Figure 1 (b)) the situation is dif-ferent. First of all we see that O − and the O +2 ion densityprofiles perfectly overlap within the bulk region and thatthe O − ion density profile decreases more steeply thanthe O +2 density profile beyond the sheath edges. Thetime averaged value of both the O − and the O +2 ion den-sity profiles within the bulk region is slightly lower thanin the 100 mTorr case. The electron density in the bulkplasma is roughly two orders of magnitude lower thanthe O +2 and O − densities. Moreover, we observe that theO + density is higher than the electron density within thebulk region, contrary to the 100 mTorr case. We also ob-serve that the O + density profile is more flattened withinthe bulk region and that it decreases more steeply thanfor the 100 mTorr case. The electron density profile isflat and constant within the bulk region and it has equalabsolute maxima on both the sheath edges. Figure 2 (a) FIG. 1: The density profiles for charged particles at (a) 100mTorr and at (b) 10 mTorr in a parallel plate capacitivelycoupled discharge with a gap separation of 45 mm driven bya 400 V voltage source at 13.56 MHz.
TABLE II: The parameters of the simulation, the particle weight, the energy threshold above which kinetics of the neutralparticles are followed and the wall recombination and quenching coefficients for the neutral species on the electrode surfaces.Species particle weight energy threshold coefficient(range) [meV] recomb./quenchingO (X Σ − g ) 5 × (a ∆ g ) 5 ×
100 0.007O (b Σ g ) 5 ×
100 0.1O( P) 5 ×
500 0.5O( D) 5 ×
50 0.5 recomb/0.5 quenchingO +2 − -O + − -O − × − -e 10 − -FIG. 2: The spatio-temporal behaviour of the (a) totalcharge density and of the (b) quasineutrality deviation,defined by Eqs. (2) and (3), respectively, over the full gaplength for a parallel plate capacitively coupled oxygendischarge at 100 mTorr for 45 mm of gap separation drivenby a 400 V voltage source at driving frequency of 13.56 MHz. shows the spatio-temporal behaviour of the total chargedensity at 100 mTorr over the full gap length for a full period defined as followsTotal Charge Density = e (cid:16) n O +2 + n O + − n O − − n e (cid:17) (2)Firstly, we observe a net zero charge density within thebulk region and the fully collapsed sheath regions. Sec-ondly, the right (left) sheath region is positively chargedand reaches its maximum extension at t/τ rf = 0 . t/τ rf = 0 . (cid:16) n O +2 + n O + − n O − − n e (cid:17) n O +2 + n O + . (3)We observe that the quasineutrality deviation uniquelyidentifies the sheath region. Indeed, the quasineutralitydeviation value is 1 within the expanded sheaths whileit is 0 within the plasma bulk and the sheath collapseregions and it has an intermediate value on the bulk-sheath time varying interface. Figure 3 (a) shows thespatio-temporal behaviour of the total charge density at10 mTorr over the full gap length for a full period de-fined by Eq. (2). Firstly, we observe a net positive (neg-ative) charged stripe that appears on the sheath side (onthe bulk side) of both the sheath edges over the full rf-cycle (on both the collapsing sheath edges). The positivecharged stripe density strongly increases on the sheathside of the expanded sheath edge. This positive chargedstripe was absent in the 100 mTorr case as shown in Fig-ure 2 (a). Such a difference with respect to the 100 mTorrcase is due to the fact that the electron mean free pathis longer at low pressures, so that they leave the positiveions behind while crossing the sheath edge during thesheath collapse. For example, at t/τ rf = 0 .
25, on very
FIG. 3: The spatio-temporal behaviour of the (a) totalcharge density and the (b) quasineutrality deviation, definedby Eqs. (2) and (3), respectively, over the full gap length fora parallel plate capacitively coupled oxygen discharge at 10mTorr for 45 mm of gap separation driven by a 400 Vvoltage source at driving frequency of 13.56 MHz. short time scales, a net negative ambipolar field buildsup (Figure 3 (a)), which induces a recall force on the bulkelectrons and a pushing force on the bulk positive ionstoward the collapsing sheath edge. For this reason, anexcess of positive charges on the immediate sheath sideof both the collapsed sheath edges is observed. It’s worthnoting that, once crossed the sheath side of the collapsedsheath edge, the electrons are free to accelerate towardsthe left electrode due to the flux compensation effect. Onthe other hand, a positive peak in the ambipolar field isobserved on the bulk side of the collapsed sheath edgeat 10 mTorr (Figure 3 (a)). Such a peak is absent at100 mTorr (Figure 2 (a)). This can be explained con-sidering that, at longer time scales, when the electronsare repelled from the sheath side of the collapsed sheathedge, the ambipolar field changes the sign and becomespositive on the bulk side of the collapsed sheath edge.Such a behaviour is responsible for the negative chargeexcess observed on the bulk side of the collapsed sheathedge. Moreover, under the action of the ambipolar force, the electrons are confined within the bulk and the col-lapsed sheath region. The same reasoning can be alsoapplied at both t/τ rf = 0 .
50 and t/τ rf = 0 .
75. We alsoobserve that the sheath regions reach a larger extensionat both t/τ rf = 0 .
25 and t/τ rf = 0 .
75 with respect tothe 100 mTorr case, while the net positive charge presenton the sheath side of both the sheath edges has a lowervalue when the sheath approaches its maximum exten-sion. The net charge in the fully expanded sheath re-gions at 10 mTorr is globally lower than in the 100 mTorrcase. Figure 3 (b) shows the spatio-temporal behaviourof the quasineutrality deviation at 10 mTorr over the fullgap length defined as in Eq. (3). As in the 100 mTorrcase shown in Figure 2 (b), the quasineutrality devia-tion uniquely identifies the sheath region. We observe anon-quasineutral stripe on both the sheath edges whichis related to the positive charged stripe observed in Fig-ure 3 (a). At 100 mTorr the sheath has a smooth contourover the full rf cycle (Figure 2 (b)) while at 10 mTorr itends abruptly on the non-quasineutral stripe on both thesheath edges. Figure 4 (a) shows the total charge den-sity profile at 100 mTorr defined by Eq. (2) at differenttime steps and time averaged over all the full gap length.We observe that the total charge density profile is flatand very close to zero over the full bulk width for allthe cases considered. At t/τ rf = 0 .
25 ( t/τ rf = 0 .
75) thetotal charge density steeply increases while approachingthe sheath edge of the right (left) bulk-sheath interfaceand sharply decreases toward the right (left) electrode.On the other hand, the total charge density is approx-imatively constant at t/τ rf = 0 .
25 ( t/τ rf = 0 .
75) andzero within the fully collapsed sheath region and slowlyincreases toward the right (left) electrode. The chargedensity profile at t/τ rf = 0 .
25 is a mirror image of thecharged density profile at t/τ rf = 0 .
75. At t/τ rf = 0 . t/τ rf = 0 .
25 and the t/τ rf = 0 .
75 cases. Inthe time averaged case the total charge density steeplyincreases once passed both the sheath edges reaching alower maximum with respect to all the other cases con-sidered. We observe the time averaged case to be almostzero within the plasma bulk and to have an almost con-stant profile within the sheath regions while approachingboth the electrodes.Figure 4 (b) shows the total charge density profile at 10mTorr defined by Eq. (3) at different time steps and timeaveraged over the full gap length. The total charge den-sity is zero within the discharge center ( x = 0) for all thefour cases considered even thought there is always somecharge density within the plasma bulk (either positive ornegative), and the time average is zero. Also, the totalcharge density profile at t/τ rf = 0 .
25 ( t/τ rf = 0 .
75) showsan additional local maximum on the left (right) sheathedge with respect to the 100 mTorr case (Fig. 4 (a)). Thesame applies to the t/τ rf = 0 .
50 and the time averaged
FIG. 4: The total charge density at (a) 100 mTorr and at(b) 10 mTorr over the full gap length at t/τ rf = 0 .
25 (greenline), t/τ rf = 0 .
50 (red dashed line), t/τ rf = 0 .
75 (bluedotted dashed line) and time averaged (black dashed line)for a parallel plate capacitively coupled oxygen discharge for45 mm of gap separation driven by a 400 V voltage source atdriving frequency of 13.56 MHz. case, with the presence of two almost equal local max-ima on both the sheath edges, which are absent in the100 mTorr case. At t/τ rf = 0 .
25 the total charge densityprofile sharply increases (decreases) while approachingthe bulk side of the right (left) sheath edge and it has anabsolute maximum (minimum) on the sheath side (bulkside) of the right (left) sheath edge. Then the total chargedensity profile sharply decreases (steeply increases) oncepassed the right (left) sheath edge reaching a positivevalue (a local maximum) toward the right electrode (onthe sheath side of the left sheath edge). We also observethat the total charge density profile is approximativelyconstant within the left sheath region.
IV. MODEL DESCRIPTION AND RESULTS
During the past years several attempts have been madeto describe correctly the behaviour of the electric field using the Boltzmann equation. Surendra and Dalvie[22] used the first momentum Boltzmann equation to de-compose the electric field into a sum of different terms,each one corresponding to different physical mechanisms.They were able to isolate all the single terms contribut-ing to both the electric field and the electron power ab-sorbed. Moreover, they found the electron pressure termto be important for the collisionless heating and that, fora constant electron temperature, the collisionless electronheating vanishes upon time average. Since then severalattempts have been made to explain the behaviour of theelectric field at different time steps and upon time aver-age. In recent years Surendra’s framework has been em-ployed to explain the behaviour of the electric field withinthe bulk and in the sheath regions. In particular Schulzeand coworkers have used the zeroth momentum Boltz-mann equation (stationary continuity equation), with astationary density profile [30] and with a temporally de-pendent density profile together with a non zero ioniza-tion rate [31], combined with the first momentum Boltz-mann equation, with and without the change in the mo-mentum term [30, 31] to derive a space- and time-resolvedexpression for the different electric field terms involved.The Surendra-Dalvie framework has improved our knowl-edge of the physical mechanisms behind the origins of theelectric field within the bulk region and within the col-lapsed sheath region but has not given a general consen-sus on the origin of the electric field within the expandedsheath region.The DA-mode is associated with the creation of electricfield within the plasma bulk. The electric field within aplasma discharge is built up by several different phenom-ena, depending on the gas considered. The electroneg-ative discharges present a bigger number of phenomenathan the electropositive discharges, and the situation ismuch more complicated. For both electropositive andelectronegative discharges sheaths form near the elec-trodes, a positive net charge within the sheath regionbuilds up, leading to a potential profile that is positivewithin the bulk region and falls to zero near both elec-trodes [41]. However, a strong electric field within thebulk region has been observed, both experimentally andby simulations in electronegative discharges. The highvalue of the electric field has been related to the low dcconductivity within the bulk as discussed by Schulze et al.[4]. Furthermore, strong peaks in the electric field atthe sheath edges have been observed [42]. The observedpeaks have been related to the corresponding local max-ima of the electron density at the sheath edges which arecaused by the ambipolar field built up by a net chargeseparation between the positive charges as they are accel-erated towards the electrode, and the electrons, togetherwith the negative ions, confined within the bulk region.This is completely different from the situation observed inthe electropositive discharges, where the ambipolar fieldaccelerates the electrons toward the discharge center [4].
A. Simple Fluid Model
When operated at 100 mTorr (10 mTorr) the elec-tronegativity is low (high) and the discharge operates inpure α -mode (hybrid DA-mode and α -mode). Irrespec-tively of the different discharge modes and conditions, asimple fluid model for an electronegative discharge is suf-ficient to describe the physics of a such system. In thissubsection the simple fluid model applied to a dischargeoperated at both 100 mTorr 10 mTorr is discussed. Themodel describes the behaviour of the electric field andof the electron power absorption within the bulk region.This model is based on the approach used by Schulzeet al. [31], with the only difference that here both theionization rate and the change of momentum terms areassumed to be negligible and are set equal to zero. Themodel is valid within the bulk region and in the col-lapsed sheath regions only, since the electron density inthe expanded sheath region is very small. Moreover, thequasineutrality condition has not been imposed and theideal gas law has been employed in the first momentumBoltzmann equation. At 100 mTorr the pressure ten-sor (the temperature) is taken as (not) isotropic. Set-ting p e ≡ p e , xx = p e , yy = p e , zz , one finds Tr( p e ,ij ) =3 p e , xx = 3 p e = 3 en e T e = en e ( T e , xx + T e , yy + T e , zz ), sothat T e ≡ ( T e , xx + T e , yy + T e , zz ) /
3, i.e. the electron tem-perature is direction averaged. Accordingly to the cur-rent set up and in order to make the physical systemconsistent, the ideal gas law has to be seen as an ap-proximation. On the other hand at 10 mTorr neither thepressure tensor, nor the temperature are isotropic. Since p e , xx = p e , yy = p e , zz and T e , xx = T e , yy = T e , zz , we are leftwith p e ≡ p e , xx = en e , xx T e , xx , p e , yy = en e , yy T e , yy , and p e , zz = en e , zz T e , zz . The zeroth and the first momentumBoltzmann equation for electrons in a plasma dischargein the absence of magnetic field are the continuity equa-tion ∂n e ∂t + ∂∂x ( u e n e ) = G − L (4)where G and L are the reaction rates involving the cre-ation and the destruction of electrons, respectively, andthe momentum balance equation ∂∂t [ m e n e u e ] + ∂∂x (cid:2) m e n e u e (cid:3) + ∂∂x [ en e T e ]+ en e E + Π c = 0 , (5)respectively. According to this set up, there is no need forkeeping the continuity equation (Eq. (4)). Now, accord-ing to Lieberman and Lichtenberg [41], the momentumchange term Π c can be approximated by a Krook colli-sional operator as followsΠ c = X β m e n e ν e β ( u e − u β ) − m e ( u e − u G ) G + m e ( u e − u L ) L (6) where the summation is over all species, u e and u β are themean velocities of the electrons and the species β , respec-tively, and ν e β is the momentum transfer frequency forcollisions between electrons and species β . Now, neglect-ing the reactions involving the creation and destructionsof particles (e.g., ionization, recombination) and consid-ering only the O neutral species, with a negligible ve-locity compared to the electrons, the momentum changeterm becomes Π c = m e ν e n e u e (7)along with the continuity equation ∂n e ∂t + ∂∂x ( u e n e ) = 0 (8)Solving Eq. (8) with respect to the velocity gradient onefinds ∂u e ∂x = − u e n e ∂n e ∂x − n e ∂n e ∂t (9)Combining Eqs. (9), (7) and (5) together with the idealgas law one finds an expression for the electric field E = − m e e ∂u e ∂t | {z } I + m e e u n e ∂n e ∂x | {z } II + m e e u e n e ∂n e ∂t | {z } III − T e n e ∂n e ∂x | {z } IV − ∂T e ∂x |{z} V − m e u e ν c e | {z } VI (10)Each electric field term in Eq. (10) has its own origin.The first and the third term (I and III) are electron in-ertia terms due to the temporal variation in the electronvelocity and density, respectively. The second term (II)corresponds to an electric field due to the normalizedelectron density gradient. The fourth (IV) term corre-sponds to diffusion (ambipolar field) [4, 30]. The fifthterm (V) corresponds to the electron temperature gradi-ent. Therefore, terms IV and V represent electron heat-ing due to pressure effects which is a collisionless mecha-nism [23]. The sixth term (VI) is due to electron collisionswith atoms and molecules (drift field). Equation (10)has been applied to a given set of input parameters. Theinput parameters are the electron density and the elec-tron temperature from the simulation. The collision term(Term VI) was taken from the reaction rate given by thesimulation for an electron neutral elastic collision. Theelectron collision frequency values at 100 and 10 mTorrare ν c = 8 . × s − and ν c = 5 . × s − withinthe discharge center, respectively. Multiplying the elec-tric field coming from Eq. (10) times the electron currentdensity J e = − en e u e it is possible to find the electronabsorbed power as follows J e · E = m e u e n e ∂u e ∂t | {z } I − m e u ∂n e ∂x | {z } II − m e u ∂n e ∂t | {z } III + eu e T e ∂n e ∂x | {z } IV + en e u e ∂T e ∂x | {z } V + m e n e ν c u | {z } VI (11)Each electron power absorption term that constitutesEq. (11) has its own origin which is strictly related to theelectric field given by Eq. (10). The first and the thirdterm (I and III) are electron inertia power absorptionterms. The second term (II) corresponds to the powerabsorption term related to the electron density gradient.The fourth (IV) term is related to the ambipolar fieldorignating from the electron density gradient [4, 30]. Thefifth term (V) is related to the electron temperature gra-dient term for the electric field (Eq. (10)). The fourth andfifth terms are usually know in the literature as pressureheating terms, respectively [22, 23, 31]. The sixth term(VI) is related to the collisions and represents ohmic heat-ing. It’s worth noting that the electron power absorptionformula shown in Eq. (11) can be split as follows [22]( J e · E ) = ( J e · E ) Nonohmic + ( J e · E ) ohmic (12)where( J e · E ) Nonohmic = Term I + Term II + Term III+ Term IV + Term V (13)( J e · E ) ohmic = Term VI (14)In turn the Non ohmic contribution can be split up asfollows [25]( J e · E ) Nonohmic = ( J e · E ) Inertia + ( J e · E ) Pressure (15)where( J e · E ) Inertia = Term I + Term II + Term III (16)( J e · E ) Pressure = Term IV + Term V (17)This will be useful later when we identify the differentcontributions to the electron power absorption. We un-derline that the same split applied to the electron powerabsorption can be applied to the electric field formulashown in Eq. (10).
V. RESULTS AND DISCUSSION
The sheath location is determined by assuming thatthe density of negatively charged species has fallen tohalf the density of the positive charged species. Inmore detail the sheath edge position, which is defined as x sh ( t ), is taken to be the position x , where the condition n e ( x, t ) /n i = 1 / t/τ rf = 0 .
25 from the left electrode to the rightsheath edge, at t/τ rf = 0 .
50 from the left to the rightsheath edge, at t/τ rf = 0 .
75 from the left sheath edge tothe right electrode. For all the three time slices an almostperfect match between the overall terms summation andthe result from the simulations is observed.At t/τ rf = 0 .
25 (Figure 5 (a)) we see that the con-tribution from Terms II and III is negligible. On theother hand, the main contribution comes from Terms I,IV, V and VI. Term I is almost zero within the bulk re-gion up to the right sheath edge while it decreases tonegative values as it approaches the left sheath edge upto the left electrode. The electric field inertia term dueto temporally varying electron velocity becomes negativeapproaching the left electrode, indicating that the elec-tron velocity gradient is positive near the left electrode(Eq. (10)). Term IV is flat and zero within the bulk re-gion and sharply increases (sharply decreases) while ap-proaching the right (left) sheath edge. An absolute min-imum in the Term IV profile is observed on the sheathside of the left sheath edge. Moreover, once passed thisminimum, the profile of Term IV is approximatively con-stant over the left sheath region and it decreases while ap-proaching the left electrode. Term V is flat and zero overthe full bulk gap length. A small local minimum (maxi-mum) in the Term V profile is observed on the sheath sideof the right (left) sheath edge. Then it steeply increases(sharply decreases) while approaching the sheath side ofthe right sheath edge (the left electrode). Term VI is flatand zero within the bulk region up to the right sheathedge while it slightly increases as it approaches the leftsheath keeping an almost constant value while approach-ing the left electrode. Finally, Figure 5 (a) shows thatthe only important contributions to the electric field at t/τ rf = 0 .
25 comes from the inertia term related to thetemporal gradient of the electron velocity (Term I), fromthe pressure gradient related terms (Term IV and V) andfrom the ohmic heating term (Term VI). At t/τ rf = 0 . t/τ rf = 0 .
50 comes from the pressure gradient relatedterms (Terms IV and V) and from the ohmic contribu-0
FIG. 5: The electric field profile of the terms that constituteEq. (10) and their summation compared with the result ofthe simulations at (a) t/τ rf = 0 .
25 from the left electrode tothe right sheath edge, at (b) t/τ rf = 0 .
50 from the left to theright sheath edge, at (c) t/τ rf = 0 .
75 from the left sheathedge to the right electrode, for a parallel plate capacitivelycoupled oxygen discharge at 100 mTorr for 45 mm of gapseparation driven by a 400 V voltage source at drivingfrequency of 13.56 MHz. tion (Term VI). At t/τ rf = 0 .
75 (Figure 5 (c)) we seealmost a mirror image of Figure 5 (a). The contribu- tions from Terms II and III are negligible and Term VIis small. The main contribution comes from Terms I, IV,and V. The only significant contributions to the electricfield at t/τ rf = 0 .
75 comes from the inertia term relatedto the temporal gradient of the electron velocity and fromthe pressure gradient related terms. Therefore, the mostsignificant contribution to the electric field within thebulk plasma at 100 mTorr is due to the pressure gradientterms. An almost perfect match between the calculatedelectric field profile and the result from simulation canbe observed for all the time steps considered as shown inFigures 5 (a), (b) and (c).Figure 6 shows the electric field profile at 10 mTorrof the terms that constitute Eq. (10) and their summa-tion compared with the result of the simulations at (a) t/τ rf = 0 .
25 from the left electrode to the right sheathedge, at (b) t/τ rf = 0 .
50 from the left to the right sheathedge, at (c) t/τ rf = 0 .
75 from the left sheath edge tothe right electrode. At t/τ rf = 0 .
25 (Fig. 6 (a)) we seethat the contribution from terms II and VI is negligiblewhile Term III is small. We observe that the main con-tribution to the electric field comes from terms IV and V.Term IV is zero within the discharge center and approx-imatively flat and zero within the inner bulk region andsharply increases while approaching the bulk side of boththe sheath edges. A lower (higher) maximum on the bulkside of the left (right) sheath edge is observed. On theother hand, Term V has a similar behaviour except thatit increases less steeply while approaching the bulk side ofthe right sheath edge. The local maximum in Term V onthe left sheath edge overlaps almost perfectly with the lo-cal maximum in Term IV, total terms summation and theresult from the simulations placed in the same locationas well as the respective profiles within the inner bulk re-gion. Finally, Figure 6 (a) shows that the only importantcontributions to the electric field at t/τ rf = 0 .
25 comesfrom the pressure gradient related terms (Terms IV andV). At t/τ rf = 0 .
50 (Fig. 6 (b)) we see that the contribu-tion from Term I is negligible and Terms II, III and VIare small. We also observe that the main contribution tothe electric field comes from terms IV and V. Term IVis flat and zero within the inner core of the plasma bulkand decreases (increases) while approaching the bulk sideof the right (left) sheath edge. Then it increases againonce passed the right sheath edge building an absoluteminimum. Finally, Figure 6 (b) shows that the only im-portant contributions to the electric field at t/τ rf = 0 . t/τ rf = 0 . t/τ rf = 0 .
75 (Fig. 6 (c))we see a mirror image of the t/τ rf = 0 .
25 case. TermsII, III and VI are negligible and Term I is small. Themain contribution to the electric field comes from termsIV and V. The local minimum in Term V on the rightsheath edge overlaps almost perfectly with the local min-1
FIG. 6: The electric field profile of the terms that constituteEq. (10) and their summation compared with the result ofthe simulations at (a) t/τ rf = 0 .
25 from the left electrode tothe right sheath edge, at (b) t/τ rf = 0 .
50 from the left to theright sheath edge, at (c) t/τ rf = 0 .
75 from the left sheathedge to the right electrode (c), for a parallel platecapacitively coupled oxygen discharge at 10 mTorr for 45mm of gap separation driven by a 400 V voltage source atdriving frequency of 13.56 MHz. imum in Term IV, total terms summation and the resultfrom the simulations placed in the same location as wellas the respective profiles within the inner bulk region. To summarize the only important contributions to the elec-tric field at t/τ rf = 0 .
75 comes from the pressure gradientrelated terms (Terms IV and V).An almost perfect match between the calculated elec-tric field profile and the result from simulation can be ob-served over the full gap length up to the bulk side of theexpanding sheath edge for all the three time slices con-sidered. Moreover, at t/τ rf = 0 .
25 and t/τ rf = 0 .
75, thecalculated electric field sharply understimates the elec-tric field coming from the simulations while approachingthe bulk side of the fully collapsed sheath edge, while at t/τ rf = 0 .
50 the difference is very small. We observe thatthe inertia term (Term I), for all the three time slicesconsidered, is negligible compared to the 100 mTorr case(Figure 5). Such a term is absent due to the presence ofthe negative charged stripes placed on bulk side of thecollapsing sheath edge, which prevents the electrons fromcrossing the collapsing sheath edge (see discussion in Sec-tion III). Moreover, due to the presence of the positivecharged stripes placed on the sheath side over the full rfperiod, the electrons are prevented from increasing theirown velocity. Since at 10 mTorr there is higher number ofelectrons within the collapsing sheath region than at 100mTorr (Figure 2 (a) and Figure 3 (a) respectively), thedisplacement current is lower and the temporal changein the electron velocity due to the time varying electricfield is also lower.Figure 7 (a) shows the time averaged electric field pro-file at 100 mTorr for the three main contributing termsto the calculated electric field using Eq. (10) from theleft to the right (time averaged) sheath edge: Term IV(red line), V (blue dashed line) and VI (green dotteddashed line). We see that all the three terms consideredare flat and zero within the bulk region. Term IV steeplyincreases (decreases) while approaching the right (left)sheath edge. On the other hand, Term V sharply in-creases (decreases) while approaching the sheath side ofthe left (right) edge. Finally, Figure 7 (a) shows that theonly important contribution to the time averaged electricfield comes from the pressure terms (Terms IV and V).Figure 7 (b) shows the time averaged electric field pro-file in a discharge operated at 10 mTorr for the threemain contributing terms to the calculated electric fieldusing Eq. (10) from the left to the right (time averaged)sheath edge: Term IV (red line), V (blue dashed line) andVI (green dotted dashed line). We see that all the threeterms considered are zero within the discharge center.Term IV sharply increases (decreases) while approachingthe bulk side of the right (left) sheath edge. Term Vsteeply decreases (increases) while approaching the bulkside of the right (left) sheath edge and it steeply increases(decreases) while crossing the right (left) sheath edge.Finally, Figure 7 (b) shows that the only important con-tribution to the time averaged electric field comes fromthe pressure terms (Terms IV and V), just like at 100mTorr.Figure 8 (a) shows the time averaged electric field at100 mTorr calculated using Eq. (10) (blue dashed line)2
FIG. 7: The time averaged electric field profile at (a) 100mTorr and at (b) 10 mTorr of Term IV (red line), Term V(blue dashed line), Term VI (green dotted dashed line) fromEq. (10) from the left to the right sheath edge for a parallelplate capacitively coupled oxygen discharge for 45 mm ofgap separation driven by a 400 V voltage source at drivingfrequency of 13.56 MHz. and the result from simulations (red line) from the leftto the right sheath edge. An almost perfect match is ob-served over the gap length considered. Figure 8 (b) showsthe time averaged electric field at 10 mTorr calculatedusing Eq. (10) (blue dashed line) and the result fromsimulations (red line) from the left to the right sheathedge. We observe an almost perfect match within theinner bulk gap length. However, the calculated electricfield overstimates (understimates) the electric field fromthe simulation on the bulk side of both the sheath edges,resembling the observed slight difference between the cal-culated and simulated electric field at different time stepson the bulk side of the expanding sheath edge (Figure 6).Figure 9 shows the spatio-temporal behavior of theelectron power absorption J e · E at 100 mTorr. The fig-ures show the electron power absorption calculated us-ing Eq. (11) (Figure 9 (a)) and from the simulation overthe full gap length (Figure 9 (b)). The ordinate cov- FIG. 8: The time averaged electric field profile calculatedusing Eq. (10) (blue dashed line) and the result fromsimulations (red line) from the left to the right sheath edgeat (a) 100 mTorr and at (b) 10 mTorr for a parallel platecapacitively coupled oxygen discharge for 45 mm of gapseparation driven by a 400 V voltage source at drivingfrequency of 13.56 MHz. ers the full rf cycle. We see that almost all the elec-tron power absorption occurs during the sheath expan-sion and that the electron power absorption mode is apure α -mode. The calculated electron power absorptionalmost perfectly matches the result from the simulationon both the sheath edges. We recall that the theoreti-cal model doesn’t take the sheath dynamics into account,since it has been built for the bulk and collapsed sheathregion only. Within the gap length we see an almost per-fect match between the electron power absorption comingfrom the theoretical model and the simulation.Figure 10 shows the spatio-temporal behavior of theelectron power absorption J e · E at 10 mTorr, where J e and E are the spatially and temporally varying electroncurrent density and electric field, respectively. The fig-ures show the electron power absorption for the theoret-ical model (Eq. (11)) (Figure 10 (a)) and from the sim-ulation results over the full gap length (Figure 10 (b)).The ordinate covers the full rf cycle. Firstly, we observe3 FIG. 9: The spatio-temporal behaviour of the electronpower absorption over the full gap length (a) calculatedusing Eq. (11) and (b) the result from simulation for aparallel plate capacitively coupled oxygen discharge at 100mTorr for 45 mm of gap separation driven by a 400 Vvoltage source at driving frequency of 13.56 MHz. that at 10 mTorr a significant power absorption (red andyellow areas) and some electron cooling (dark blue areas)are evident within the plasma bulk region. The electronpower absorption appears during the sheath expansion(collapse) on the sheath side (on the bulk side) of thesheath edge, while there is electron cooling during thesheath expansion (collapse) on the bulk side of the sheathedge (on the electrode side). Indeed, at 10 mTorr theelectron heating mechanism is a combination of a driftambipolar (DA) heating in the bulk plasma and stochas-tic heating due to the sheath oscillation ( α -mode), as ithas been shown in our previous works [19–21]. Secondly,we observe that the calculated electron power absorptionresembles well the result from the simulation which isslightly overstimated within the sheath region only.Figure 11 shows the electron power absorption profileat 100 mTorr calculated using Eq. (11) (blue dashed line)and the result from simulations (red line) at (a) t/τ rf =0 .
25 from the left electrode to the right sheath edge, at(b) t/τ rf = 0 .
50 from the left to the right sheath edge,
FIG. 10: The spatio-temporal behaviour of the electronpower absorption over the full gap length calculated (a)using Eq. (11) and (b) the result from simulation for aparallel plate capacitively coupled oxygen discharge at 10mTorr for 45 mm of gap separation driven by a 400 Vvoltage source at driving frequency of 13.56 MHz. at (c) t/τ rf = 0 .
75 from the left sheath edge to the rightelectrode. An almost perfect match is observed for allthe time steps considered as shown in Figures 11 (a), (b)and (c).Figure 12 shows the electron power absorption profileat 10 mTorr calculated using Eq. (11) (blue dashed line)and the result from simulations (red line) at t/τ rf = 0 . t/τ rf = 0 .
50 from the left to the right sheath edge (b),at t/τ rf = 0 .
75 from the left sheath edge to the rightelectrode (c). An almost perfect match is observed forall the time steps considered as shown in Figures 12 (a),(b) and (c).Figure 13 (a) shows the time averaged electron powerabsorption profile at 100 mTorr of Term IV (red line), V(blue dashed line) and VI (green dotted dashed line) cal-culated using Eq. (11) from the left to the right sheathedge. All the three terms considered are flat and zerowithin the bulk gap length. We observe that Term IV(Term V) steeply increases (decreases) while approaching4
FIG. 11: The electron power absorption profile fromEq. (11) (blue dashed line) at (a) t/τ rf = 0 .
25 and the resultfrom simulations (red line) from the left electrode to theright sheath edge, at (b) t/τ rf = 0 .
50 from the left to theright sheath edge, at (c) t/τ rf = 0 .
75 from the left sheathedge to the right electrode, for a parallel plate capacitivelycoupled oxygen discharge at 100 mTorr for 45 mm of gapseparation driven by a 400 V voltage source at drivingfrequency of 13.56 MHz. the sheath side of both the sheath edges. On the otherhand, we observe that Term VI slightly increases while
FIG. 12: The electron power absorption profile fromEq. (11) (blue dashed line) at (a) t/τ rf = 0 .
25 and the resultfrom simulations (red line) from the left electrode to theright sheath edge, at (b) t/τ rf = 0 .
50 from the left to theright sheath edge, at (c) t/τ rf = 0 .
75 from the left sheathedge to the right electrode, for a parallel plate capacitivelycoupled oxygen discharge at 10 mTorr for 45 mm of gapseparation driven by a 400 V voltage source at drivingfrequency of 13.56 MHz. crossing both the sheath edges reaching small positivevalues on the sheath side of both the sheath edges. Fi-5
FIG. 13: The time averaged electron power absorptionprofile of Term IV (red line), Term V (blue dashed line),Term VI (green dotted dashed line) from Eq. (11) from theleft to the right (time averaged) sheath edge at (a) 100mTorr and at (b) 10 mTorr for a parallel plate capacitivelycoupled oxygen discharge for 45 mm of gap separationdriven by a 400 V voltage source at driving frequency of13.56 MHz. nally, Figure 13 shows that the main contributions to thetime averaged electron power absorption at 100 mTorrcomes from the pressure gradient related terms (TermIV and V) and from ohmic term (Term VI).Figure 13 (b) shows the time averaged electron powerabsorption profile at 10 mTorr of Term IV (red line), V(blue dashed line) and VI (green dotted dashed line) cal-culated using Eq. (11) from the left to the right sheathedge. We observe that all the three terms considered havea parabolic behaviour over the bulk gap length. More-over, Term V is higher (lower) on the bulk side of both thesheath edges (within the inner electronegative core) thanTerm VI, while Term IV is sharply lower over the full gaplength. In more detail, we see that Term IV is almost flatand zero within the discharge center up to the bulk sideof both the sheath edges, and that it slightly increaseswhile approaching both the sheath edges building equallocal maxima. Then it steeply decreases while crossing both the sheath edges, building two almost equal localminima. On the other hand, Term V slightly increaseswhile approaching the bulk side of both the sheath edges,building two almost equal local maxima, and it sharplydecreases while approaching both the sheath edges. Asregards to Term VI we see that it builds an absolutemaximum at the discharge center and that it sharply de-creases while approaching both the sheath edges.
FIG. 14: The time averaged electron power absorptionprofile calculated using Eq. (11) (a) (blue dashed line) andthe result from simulations (red line) over the full gap lengthat (a) 100 mTorr and at (b) 10 mTorr for a parallel platecapacitively coupled oxygen discharge for 45 mm of gapseparation driven by a 400 V voltage source at drivingfrequency of 13.56 MHz.
Figure 14 (a) shows the comparison between thetime averaged electron power absorption calculated usingEq. (11) (blue dashed line) and the result from simula-tion (red line) over the full gap length at 100 mTorr. InFigure 14 (a) we see that the calculated time averagedelectron power absorption overlaps almost perfectly withthe result from the simulation over the bulk gap lengthup to the sheath side of both the sheath edges. Moreover,it slightly understimates the result from the simulationwithin the inner core of both the sheath regions up tothe respective electrodes.6Figure 14 (b) shows a comparison between the timeaveraged electron power absorption calculated usingEq. (11) (blue dashed line) and the result from simu-lation (red line) over the full gap length at 10 mTorr. InFigure 14 (b) we see that the calculated time averagedelectron power absorption calculated using Eq. (11) over-laps almost perfectly with the result from the simulationover the bulk gap length up to the sheath side of boththe sheath edges. Moreover, the closer to the electrodes,the bigger is the difference between the calculation andthe simulation.
FIG. 15: The space-time averaged electron power absorptionprofile terms calculated using Eq. (11) at 10 mTorr (bluebars) and 100 mTorr (red bars) for a parallel platecapacitively coupled oxygen discharge for 45 mm of gapseparation driven by a 400 V voltage source at drivingfrequency of 13.56 MHz.
Figure 15 shows a comparison between the space-timeaveraged electron power absorption within the plasmabulk at 10 mTorr calculated using Eq. (11) (blue bars)and at 100 mTorr calculated using Eq. (11) (red bars).The space average has been taken within the bulk regiononly using the time averaged sheath locations as alreadydiscussed at the beginning of Section V. The histogramhas been constructed considering the calculated space-time averaged electron absorbed power and then buildingthe following quantity( J e · E ) X Term Percentage = 100 × ( J e · E ) X Term ( J e · E ) (18)where ( J e · E ) X Term labels the space-time averaged elec-tron power absorption related to the X term, where X refers to term I, II ... VI in Eq. (11). In Fig. 15 we ob-serve that in the 100 mTorr case the space-time averagedelectron power absorption comes from the pressure terms(Term IV and V) and the ohmic Term (Term VI). More-over, we see that Term IV is positive while Term V is neg-ative and they share almost the same magnitude in theabsolute value, while Term VI is the smallest. Therefore,the ambipolar term is a power absorption term while the electron temperature gradient presents power loss (elec-tron cooling). The main electron power absorption at100 mTorr is due to the pressure gradient terms. At 10mTorr the situation has changed drastically. First of allTerm IV and Term V flip their sign with respect to the100 mTorr case and are sharply smaller in the absolutevalue compared to the 100 mTorr case. In this contextthe ohmic term’s magnitude (Term VI) has not signifi-cantly changed compared to the 100 mTorr case but nowshares the same order of magnitude with respect to bothTerm IV and Term V. Moreover, in Figure 15 we recog-nize a general pattern where the ohmic term (pressureterm) contribution in the space-time averaged total elec-tron power absorption increases (decreases) when the thetotal pressure decreases (increases). Such a behaviour isexpected and has also been observed by Vass et al. [32].The only difference is that Vass et al. [32] find the pres-sure term contribution in the space-time averaged elec-tron power absorption to be negligible at low pressure.Like Vass et al. [32] we find the ohmic power absorp-tion to be important even at low pressure. Finally, at10 mTorr we observe the presence of small contributionscoming from Term I, II and III respectively.The inertiaterms I and III provide electron cooling while the electrondensity gradient term II contributes to electron power ab-sorption. VI. CONCLUSION
The one-dimensional object-oriented particle-in-cellMonte Carlo collision code oopd1 was applied to explorethe properties of the electric field and the electron powerabsorption at different time steps and time averaged overa full rf cycle within the plasma bulk in a capacitivelyoxygen coupled discharge at both 100 and 10 mTorr for45 mm gap distance. At both 100 and 10 mTorr the fluidmodel presented by Schulze et al. [31] was applied.At 100 mTorr at both t/τ rf = 0 .
25 and t/τ rf = 0 .
75 themain contributions to both the electric field and the elec-tron power absorption are due to the electron inertia termrelated to the temporal gradient of the electron velocity(Term I), the gradient pressure related terms (Term IVand V) and the ohmic heating term. At t/τ rf = 0 . t/τ rf = 0 .
25 and t/τ rf = 0 .
75 the maincontributions to both the electric field and the electronpower absorption come from the pressure gradient relatedterms (Term IV and V) only, while at t/τ rf = 0 .
50 asmall additional but not negligible contribution from thedrift field (Term VI) has been observed. Moreover, in7the time averaged case, the main contributions to theelectron power absorption come from both the drift field(Term VI) and the pressure gradient related terms (TermIV and V). We have also shown that the pressure gradientrelated terms and the ohmic term contribute to the timeaveraged electron power absorption.
Acknowledgments
This work was partially supported by the Icelandic Re-search Fund Grant No. 163086 and the University of Ice- land Research Fund.
Data Availability
The data that support the findings of this study areavailable from the corresponding author upon reasonablerequest. [1] S. Wilczek, J. Schulze, R. P. Brinkmann, Z. Donk´o,J. Trieschmann, and T. Mussenbrock, Journal of AppliedPhysics , 181101 (2020).[2] M. Lieberman and V. Godyak, IEEE Transactions onPlasma Science , 955 (1998).[3] V. A. Godyak and A. S. Khanneh, IEEE Transactions onPlasma Science , 112 (1986).[4] J. Schulze, A. Derzsi, K. Dittmann, T. Hemke, J. Meich-sner, and Z. Donk´o, Physical Review Letters , 275001(2011).[5] T. Mussenbrock and R. P. Brinkmann, Applied PhysicsLetters , 151503 (2006).[6] U. Czarnetzki, T. Mussenbrock, and R. Brinkmann,Physics of Plasmas , 123503 (2006).[7] T. Mussenbrock, R. P. Brinkmann, M. A. Lieberman,A. J. Lichtenberg, and E. Kawamura, Phys. Rev. Lett. , 085004 (2008).[8] M. A. Lieberman, A. J. Lichtenberg, E. Kawamura,T. Mussenbrock, and R. P. Brinkmann, Physics of Plas-mas , 063505 (2008).[9] Z. Donk´o, J. Schulze, U. Czarnetzki, andD. Luggenh¨olscher, Applied Physics Letters ,131501 (2009).[10] G. Y. Park, S. J. You, F. Iza, and J. K. Lee, PhysicalReview Letters , 085003 (2007).[11] Y.-X. Liu, Q.-Z. Zhang, W. Jiang, L.-J. Hou, X.-Z. Jiang,W.-Q. Lu, and Y.-N. Wang, Physical Review Letters , 055002 (2011).[12] J. T. Gudmundsson, E. Kawamura, and M. A. Lieber-man, Plasma Sources Science and Technology , 035011(2013).[13] J. T. Gudmundsson and M. A. Lieberman, PlasmaSources Science and Technology , 035016 (2015).[14] H. Hannesdottir and J. T. Gudmundsson, PlasmaSources Science and Technology , 055002 (2016).[15] J. T. Gudmundsson and B. Vent´ejou, Journal of AppliedPhysics , 153302 (2015).[16] J. T. Gudmundsson and H. Hannesdottir, AIP Confer-ence Proceedings , 120001 (2017).[17] J. T. Gudmundsson and D. I. Snorrason, Journal of Ap-plied Physics , 193302 (2017).[18] J. T. Gudmundsson, D. I. Snorrason, and H. Hannesdot-tir, Plasma Sources Science and Technology , 025009(2018).[19] A. Proto and J. T. Gudmundsson, Atoms , 65 (2018). [20] A. Proto and J. T. Gudmundsson, Plasma Sources Sci-ence and Technology , 074002 (2018).[21] J. T. Gudmundsson and A. Proto, Plasma Sources Sci-ence and Technology , 045012 (2019).[22] M. Surendra and M. Dalvie, Physical Review E , 3914(1993).[23] M. M. Turner, Physical Review Letters , 1312 (1995).[24] G. Gozadinos, M. M. Turner, and D. Vender, PhysicalReview Letters , 135004 (2001).[25] T. Lafleur, P. Chabert, and J. P. Booth, Plasma SourcesScience and Technology , 035010 (2014).[26] Y. Liu, J.-P. Booth, and P. Chabert, Plasma SourcesScience and Technology , 025006 (2018).[27] R. P. Brinkmann, Plasma Sources Science and Technol-ogy , 014001 (2016).[28] M. J. Grapperhaus and M. J. Kushner, Journal of Ap-plied Physics , 569 (1997).[29] I. Kaganovich, V. Kolobov, and L. Tsendin, AppliedPhysics Letters , 3818 (1996).[30] J. Schulze, Z. Donk´o´o, D. L. B G Heil1, T. Mussenbrock,R. P. Brinkmann, and U. Czarnetzki, Journal of PhysicsD: Applied Physics , 105214 (2008).[31] J. Schulze, Z. Donk´o, T. Lafleur, S. Wilczek, and R. P.Brinkmann, Plasma Sources Science and Technology ,055010 (2018).[32] M. Vass, S. Wilczek, T. Lafleur, R. P. Brinkmann,Z. Donk´o, and J. Schulze, Plasma Sources Science andTechnology , 025019 (2020).[33] C. Birdsall, IEEE Transactions on Plasma Science , 65(1991).[34] E. Kawamura, C. K. Birdsall, and V. Vahedi, PlasmaSources Science and Technology , 413 (2000).[35] E. G. Thorsteinsson and J. T. Gudmundsson, PlasmaSources Science and Technology , 055008 (2010).[36] A. Derzsi, T. Lafleur, J.-P. Booth, I. Korolov, andZ. Donk´o, Plasma Sources Science and Technology ,015004 (2016).[37] L. Wang, D.-Q. Wen, Q.-Z. Zhang, Y.-H. Song, Y.-R.Zhang, and Y.-N. Wang, Plasma Sources Science andTechnology , 055007 (2019).[38] H. Hannesdottir and J. T. Gudmundsson, Journal ofPhysics D: Applied Physics , 175201 (2017).[39] A. J. Lichtenberg, V. Vahedi, M. A. Lieberman, andT. Rognlien, Journal of Applied Physics , 2339 (1994).[40] A. J. Lichtenberg, I. G. Kouznetsov, Y. T. Lee, M. A. Lieberman, I. D. Kaganovich, and L. D. Tsendin, PlasmaSources Science and Technology , 437 (1997).[41] M. A. Lieberman and A. J. Lichtenberg, Principles ofPlasma Discharges and Materials Processing (John Wi- ley & Sons, New York, 2005), 2nd ed.[42] V. Georgieva, A. Bogaerts, and R. Gijbels, Journal ofApplied Physics93