Self-bias voltage formation and charged particle dynamics in multi-frequency capacitively coupled plasmas

In this work, we analyze the creation of the discharge asymmetry and the concomitant formation of the DC self-bias voltage in capacitively coupled radio frequency plasmas driven by multi-frequency waveforms, as a function of the electrode surface characteristics. For this latter, we consider and vary the coefﬁcients that characterize the elastic reﬂection of the electrons from the surfaces and the ion-induced secondary electron yield. Our investigations are based on Particle-in-Cell/Monte Carlo Collision simulations of the plasma and on a model that aids the understanding of the computational results. Electron reﬂection from the electrodes is found to affect slightly the discharge asymmetry in the presence of multi-frequency excitation, whereas secondary electrons cause distinct changes to the asymmetry of the plasma as a function of the phase angle between the harmonics of the driving voltage waveform and as a function the number of these harmonics.


I. INTRODUCTION
Capacitively coupled plasma sources (CCPs) driven by radio-frequency (RF) waveforms have been aiding plasma processing industry for decades.As RF current can flow through dielectric substances as well, the electrode materials are not restricted to conducting ones, which tremendously widens the range of applications.Utilising the bombardment by high energy ions, e.g., material can be removed from surfaces (plasma etching in microelectronics).At low bombarding energies, deposition from the plasma prevails (e.g. in fabrication of photovolatic devices).2][3][4] .
3][14] At high pressures, the ions flying through the sheaths collide several times with the atoms/molecules of the background gas and have, consequently, a low energy upon arrival at the electrode surfaces.In contrast, at low pressures, the ions have a long free path and can gain high energies while flying through the sheaths.When the ions cross the sheaths in a fraction of the RF period, their energy is determined by the instantaneous sheath voltage.When the ion transit time is much longer than the RF period the ion energy is largely determined by the time-averaged sheath voltage.][17] During the past decades, a number of approaches have been developed to provide additional "degrees of freedom" to control ion properties at the electrodes.9][20][21][22][23][24] The DF excitation utilises the "functional separation" of these two excitation signals: the highfrequency signal is responsible for the creation of the plasma whereas the low-frequency voltage is responsible for the acceleration of the ions.This way the amplitude of the highfrequency component controls the plasma density and, consequently, the ion flux at the surfaces while the amplitude of the low frequency signal controls the energy of the ions.The functional separation is most efficient when the two excitation frequencies are significantly different, but even in this case frequency coupling effects hinder the efficient control of the ion properties for this "classical" dual-frequency excitation. 25,26nother major step has been the discovery of the Electrical Asymmetry Effect (EAE) that allows to make geometrically symmetric plasma sources electrically asymmetric. 5This is achieved by applying a base RF and its second harmonic for the excitation of the plasma, which leads to a development of a DC self-bias voltage.It was explained by theory 27 and subsequently confirmed by both simulations 28 and experiments 29 that the self-bias voltage can be controlled by the phase angle between the driving voltage harmonics.As the self-bias voltage influences the voltage drops over the sheaths, the ion energy can readily be controlled whereas the ion flux as shown by subsequent studies can be kept at a reasonably constant level.The EAE also develops when geometrically asymmetric discharges are driven with specific waveforms and allows controlling the discharge properties within a wide range. 30,31tudies of the EAE were also extended to a higher number of harmonics (N > 2), various special waveforms like peaksand valleys-waveforms, 32 as well as sawtooth-waveforms 33 have been introduced and investigated both experimentally and computationally.These waveforms, known as "Tailored Voltage Waveforms" (TVW) 34 , have been shown to provide large flexibility for controlling charged particle dynamics, the spatio-temporal distribution of the rates of elementary processes (e.g.ionization and excitation), the electron energy distribution function, as well as the ion properties.As peaksand valleys-waveforms have markedly different positive and negative peak amplitudes these cause an amplitude asymmetry effect in the discharge.Sawtooth-type waveforms, on the other hand, have equal positive and negative peak amplitudes, however, notably different rising and falling slopes.These result in different sheath expansion velocities and, consequently, different rates of excitation and ionization at the two sides of the discharge, which generates an asymmetry and a DC selfbias voltage, termed as the slope asymmetry effect. 35n electrically asymmetric discharges, the DC self-bias voltage (η) has a direct effect on the ion flux-energy distribution (IFED) at the electrodes.7][38] Moreover, in the presence of η = 0, the IFED-s at the two electrodes will be different, which may be advantageous in applications when a high ion energy is required at one electrode, while this is to be prevented at the other electrode.The dependence of the self-bias voltage on the amplitudes and the phases of multi-frequency waveforms has thoroughly been investigated. 27,39,40A discharge asymmetry was also found to be induced by differing materials of the two electrodes, represented by, e.g., different electron reflection probabilities 41 or different secondary electron yields 42,43 or the combination of these. 44][47][48] The investigation of the nonlinear coupling of these various asymmetry effects clearly warrants further studies.Note, that most plasma reactors used in industrial applications are geometrically asymmetric.When such a discharge is driven by a multi-frequency waveform and has different electrode materials, three types of asymmetry effects are present simultaneously.If a magnetic field is applied as well, then four non-linearly coupled types of asymmetry mechanisms will be present.Moreover, many of these applications use an electronegative gas or gas mixture, in which the asymmetry effects may differ significantly [49][50][51][52][53][54] from those in thoroughly investigated electropositive discharges.
The formation of the DC self-bias voltage and its dependence on the properties of the applied waveform have been studied experimentally and via simulations in a number of studies.6][57][58] This particle based approach is fully capable to capture kinetic effects, 59 therefore is well suited for the description of plasma sources operated at low pressures where non-local particle transport [60][61][62][63] appears.An analytical model 27 based on the voltage balance of the discharge has aided the understanding of the observations.While a lot of knowledge has accumulated in the previous studies (part of which has been reviewed above) some details of the discharge dynamics and the self-bias formation in CCPs require further studies.Our aim here is to understand the effects of the number of harmonics used to construct the excitation wavefrom and to reveal how these vary as a function of the parameters of the surface processes: (i) the reflection coefficient of the electrons at the electrodes and (ii) the ioninduced secondary electron yield.Here, the same values of these parameters will be assumed for both electrodes.
Following the introduction of the physical setting considered, in Section II the analytic model and the basics of the computational method will be described in Sections II A and II B, respectively.Section III is devoted to the presentation of the results, where we provide a very detailed analysis that goes beyond the details covered in previous studies.In particular, we examine the effects of the floating sheath potentials and the finite voltage drop over the plasma bulk on the discharge asymmetry and the self-bias voltage.Subsequently, the effects of the surface processes are discussed.A brief summary is given in Section IV.

II. PHYSICAL SYSTEM AND METHODS
In this work, we consider a capacitively coupled plasma source that has two parallel, planar electrodes.The diameter of the electrodes is assumed to be much larger than the gap between them, allowing the use of a one-dimensional model.The discharge is excited by a voltage waveform 39 where ω 1 = 2π f 1 , with f 1 being the "base" radio frequency, φ k and θ k , respectively, the amplitude and the phase of the k-th harmonic.The amplitudes of the individual harmonics are set according to φ * is usually called the peak-to-peak voltage of the waveform.Note, however, that this is true only if Eq. (1) generates peaksor valleys-types waveforms.The first of these cases is realised by setting θ k = 0, ∀k, while the second case is realised by setting θ k = 0 for the odd values of k and θ k = 180 • for even values of k.By keeping the phase angles of all odd harmonics 0 • and varying the common value, denoted by θ , for all even harmonics, various waveforms (which include both the peaks and valleys cases) can be realised as shown in Figure 1 for N = 2 and N = 4.One can note that for arbitrary values of θ , the peak-to-peak amplitude of the waveform specified by (1) indeed varies, at θ = 90 • , e.g., φ pp ≈ 1.15φ * for N = 2 and φ pp ≈ 1.23φ * for N = 4. Peaks-type waveforms (θ = 0 • ) have sharp positive peaks and nearly flat negative parts between these peaks.Correspondingly, the sheath at the powered electrode is expanded for a relatively long part of the fundamental RF period, whereas the sheath at the grounded electrode is expanded for a short time.For valleys-type waveforms (θ = 180 • ) the scenario is reversed.

A. Model for the DC self-bias voltage formation
Figure 2 shows the equivalent electrical circuit consisting of the RF generator (G), the blocking capacitor (C), as well as the plasma that is represented by three circuit elements corresponding to the three regions of the discharge: the sheath at the powered side of the plasma (or "powered sheath"), the plasma bulk, and the sheath at the grounded side of the plasmas (or "grounded sheath"). 27Two of these elements, the sheaths, exhibit capacitive impedance, while the impedance of the bulk region consists of a resistive and an inductive part, originating from, respectively, the finite conductivity due to electron-atom collisions and from the inertia (finite mass) of the electrons. 64he balance equation for the voltage components marked in Figure 2 is Here, φ C is the DC voltage drop over the blocking capacitor, we assume the AC voltage drop over this element is negligible due to its high capacitance.In this case, the DC voltage drop φ C is the opposite of the DC self-bias voltage η that develops over the plasma due to the EAE. 27As a consequence, Eq. ( 3) can be rewritten as The model of the EAE, which assumes that (i) the sheath are fully collapsed at one side of the plasma at times of the extrema of the applied voltage waveforms and that (ii) there is no voltage drop over the bulk region of the plasma, predicts the dc self-bias voltage, based on the voltage balance of the circuit, to be where φ max and φ min are, respectively, the maximum and the minimum of the applied voltage waveform, φ (t).The more general expression, which considers the nonzero sheath voltages upon sheath collapse (i.e. the floating potentials) and the finite voltage drop over the plasma bulk, 53 is Here, we already introduced notations for the contributions of different origin: due to the waveform, η w , to the floating potentials, η f , and to the bulk voltage drop, η b .
In the above expressions, ε is the symmetry parameter, which is the magnitude of the ratio of the peak values of the sheath voltages at both sides of the plasma 27 : Calculations, which are not repeated here, express the extrema of the sheath voltages as 5 Here, e is the elementary charge, ε 0 the permittivity of free space, Q mp/mg the maximum charges within the sheaths, A p/g the surfaces, I sp/sg the sheath integrals, and n sp/sg the mean charged particle densities within the sheaths at the powered (p) and grounded sides (g) of the system.Note, that upon the original derivation of these expressions, 5 the electron front was assumed to exhibit a step profile at the sheath edge and the ion density profile was taken to be static within the sheath regions.Here, however, the data obtained from the PIC/MCC simulations include the slight penetration of a finite electron density into the sheaths, i.e.Q and n represent, respectively, the net charge and charged particle density.
The ratio of the sheath integrals appearing in the above expressions is customarily approximated 27 by a value of 1.0.In the symmetric system considered here, A p = A g = A. With these two simplifications, Eq. ( 7) becomes: Upon the presentation of the results, the most precise of ε will be computed from this equation, however, the validity of simplified approaches will also be tested, as follows: • The simplest approximation for ε is represented by a value of ε = 1, which neglects the differences between the magnitudes of the peak sheath voltages at the two electrodes.This approximation is termed as 'Model 1'.
• As a refinement, one may calculate ε with neglecting the difference between Q mp and Q mg in Eq. ( 10), i.e. taking which is the form that was used in the first model of the EAE. 5 This is our 'Model 2'.
The origin of any discharge asymmetry can also be approached from the most important elementary process in the plasmas: the ionization.This process is the primary source of the charged particles and under the conditions studied here, is driven by high-energy electrons that represent a minor fraction of the electron population. 59The two basic ways of gaining enough energy for ionization (relevant at our conditions) are: acceleration of the electrons (i) near the edges of the expanding sheaths ("α-heating") 66 and (ii) within the sheaths in the strong electric field whenever (secondary) electrons are emitted from the electrodes, due to, e.g., ion bombardment ("γ-heating"). 67he causes of asymmetry effects discussed above, like different positive vs. negative values or different rising vs. falling slopes of the driving voltage waveform, as well as different secondary electron yields at the two electrodes can also be viewed to act via establishing an imbalance of the ionization at the two sides of the plasma.A faster sheath expansion (controlled by the driving waveform), e.g., gives rise to a higher energy gain of the electrons and, generally, creates a higher charge density at the corresponding side of the plasma.In the presence of secondary electrons, the magnitude of the sheath voltages (accelerating these electrons) and the duration of the expanded phase of the sheath are the important factors.A higher sheath voltage and/or a longer expanded phase of the sheath gives rise to higher ionization rate.All these effects can couple in a complicated nonlinear way in a CCP.

B. Computational method
Our numerical results are obtained from one-dimensional (1D3V) bounded electrostatic PIC/MCC simulations.As this is a well-established method, the description of its details is omitted here, only some details specific to the current study are outlined below.More information about the approach can be found in the literature. 68ur code considers electrons and Ar + ions and follows their motion in an electric field that is defined by the potentials of the electrodes and the presence of the charged particles in the electrode gap.The powered electrode (situated at x = 0) is at a potential φ (t) + η, while the other electrode (situated at x = L) is grounded (φ (t) is defined by Eq. ( 1)).
The equation of motion of the charged particles is integrated using the leapfrog scheme, with a time step of ∆t = T /3000.The computational grid for the potential, the electric field, and the charged particle densities (that has a spatial resolution of ∆x) comprises 500 points.These parameters fulfil the relevant stability criteria of the PIC/MCC method. 69,70t the electrode surfaces, as already mentioned in Section I, two processes are considered.(i) Ar + ions arriving at the surface induce the emission of a secondary electron with a probability that is expressed by the secondary electron yield, γ. (ii) Electrons arriving at the electrode surfaces undergo an elastic reflection event with a probability R (of which the dependence on energy and angle of incidence is not taken into account).
The DC self-bias voltage of the discharges driven by N > 1 harmonics is determined in an iterative manner. 28At the initialization of the simulation, η = 0 V is set.After executing the simulation for a given number (typically 50) of RF cycles, the currents of the electrons and argon ions reaching each electrode are compared.Depending on the balance of these currents, the self-bias voltage is changed by a small quantity.This procedure is continued until η reaches a converged value and the time-averaged charged particle currents to each of the two electrodes balance over (within the noise level).
For our studies the identification of the position of the RF sheath edge as a function of time, s(t), is a crucial task.It is carried out based on computing the spatially and temporally resolved distributions of the electron and ion densities. 71The position of the sheath sheath edge, e.g., at the electrode at x = 0 (i.e.s = s p ) can be found from Here, h is a position where quasineutrality holds, we set this value as h = L/2.Solving the above equation for each time step within the RF cycle, the s(t) function can be determined at both sides of the discharge.When s p (t) and s g (t) are known, the voltage drops over the sheaths (φ sp (t) and φ sg (t)), the net space charges (Q p (t) and Q g (t)) and mean net charged particle densities (n sp (t) and n sg (t)) within the sheaths can readily be determined.

III. RESULTS
The simulations are carried out for Ar gas, at fixed values of the pressure, p = 10 Pa, the electrode gap, L = 2.5 cm, and the base frequency, f 1 = 13.56MHz.Driving voltage waveforms specified by Eq. ( 1) will be used with φ * = 300 V and with up to N = 4 harmonics, with phase angles over the whole domain of interest (0 • ≤ θ ≤ 360 • ).For the surface coefficients we adopt the following values.(i) For the secondary electron yield we take γ = 0, 0.2, and 0.4.For metal surfaces, γ is expected to be rather small, 0.1, while for dielectric surfaces the actual values may be well approximated by the higher γ-s adopted.(ii) For the elastic reflection coefficient of the electrons we take R = 0 and 0.2.
We start the presentation of the results by discussing the temporal behavior of the various quantities that appear in the voltage balance equation ( 4), for the cases of single-(N = 1) and dual-frequency (N = 2) excitation.The N = 1 case is displayed in Figure 3(a), for the base conditions of φ 1 = 150 V, f 1 = 13.56MHz, p = 10 Pa, L = 2.5 cm, at zero values of the surface coefficients.At this electrically symmetric excitation waveform the plasma is symmetric and no self-bias voltage develops.The magnitudes of the sheath voltages (φ sp and φ sg ) vary with opposite phase.Note, that φ sp ≤ 0 and φ sg ≥ 0 (cf.Eqs. ( 8) and ( 9)).The minima of |φ sp | and |φ sg | at the extrema of the applied voltage amount a few Volts.These are the socalled floating potentials, φ f sp and φ f sg , respectively, at the powered and at the grounded electrodes, which limit the losses of the electrons to the electrode where the sheath momentarily collapses.The sum of the sheath voltages approximates quite well the discharge voltage (eq.( 4)) as the voltage drop over the bulk of the plasmas, φ b , amounts ± few Volts only, due to the high conductivity of the plasma.
The results for the N = 2 case (with keeping all other parameters the same) are shown in Figure 3 1), with amplitudes given by Eq. ( 2), φ * = 300 V. T is the period of the fundamental RF frequency, f 1 .Note, that the bulk voltage drop is multiplied by a factor of 10.
A small additional contribution is provided by the nonzero voltage drop over the bulk region.As compared to the N = 1 case, now the behavior of the two sheaths is quite different.Due due to the specific applied voltage waveform, the sheath at the powered electrode collapses once within the "principal" RF cycle (T = 1/ f 1 ) at t/T = 0, while at the grounded side the sheath collapses twice during this period.Moreover, the sheath stays collapsed in the latter case for a longer time.As a consequence, the floating potential has a higher value of φ f sg ≈ 4.6 V as compared to the φ f sp ≈ −1.1 V found at the powered electrode.(Note, that these values are not resolved in the figure).These potentials ensure the compensation of electron and ion currents over an RF period by regulating the electron fluxes that reach the electrodes.
The time dependence of the length of the sheaths for this case is presented in Figure 4.The sheath at the powered side has a minimum length of about 0.06 cm, whereas at the grounded side the minimum of s g is ≈ 0.1 cm.This figure also shows the net charge contained within the sheaths, for a unit electrode area of 1 cm 2 .The temporal change of Q p and Q g FIG. 4. Time dependence of the sheath lengths (chain lines, left scale) and net charge in the individual sheaths (Q p and Q g ) and their sum (Q tot = Q p + Q g ), per unit area of A = 1 cm 2 (thick solid lines, right scale), for the N = 2, θ = 0 • case.The other conditions are the same as in Figure 3.
closely follows the variation of the length of the corresponding sheath.The sum of the charges in the two sheaths, Q tot is almost invariant of time, as Figure 4 reveals.The slow drop (small negative slope) of Q tot (t) is due to the continuous ion flux to the electrodes, while the temporary increase upon times of sheath collapses are due to the losses of the electrons to the electrodes.
After illustrating the time-dependent behavior of the relevant physical quantities for a few selected cases, next we address the behavior of the various voltages as a function of the phase angle θ .Figure 5 shows the maximum and minimum values of the applied voltage (φ max and φ min ), the peak sheath voltages ( φ sp and φ sg ), the floating potentials (φ f sp and φ f sg ), the bulk voltage drops at the times of the maximum and minimum of the applied voltage (φ b max and φ b min ), as well as the DC selfbias voltage, η, obtained directly from the simulation.(Recall that Eq. ( 6) formulates a connection between these quantities based on theory, but the comparison of these η values with the simulation results is presented later.) The difference between φ max and φ min equals φ * (being fixed at a value of 300 V) only at θ = 0 • and 180 • , at other values φ max − φ min > φ * .This disparity between φ max and |φ min | is higher in the case of N = 4, as compared to the N = 2 case, according to the increasing asymmetry of the applied waveform (see Figure 1).The maxima of the magnitudes, | φ sp | and φ sg as a function of θ , are on the other hand, very similar in the two cases with different number of harmonics (N).The floating sheath potentials show also very similar patterns in the N = 2 and N = 4 cases, while the voltage drop over the bulk plasma (although this is a small value), grows to almost a factor of two higher when the number of harmonics is doubled from N = 2 to N = 4.As to the DC self-bias voltage η, the highest values are obtained near, but not exactly at θ = 0 • and 180 • . 65 on the balance between the electron and ion currents to the electrodes.
Next, we address the question how well the model of the EAE, outlined in Section II A, reproduces these results for the self-bias voltage.For this, (i) the importance of the different terms in the expression ( 6) are examined, and (ii) various approximations ("modeling levels") for the calculation of the symmetry parameter ε (see the end of Sec.II A) are tested.This analysis is aided by Figure 6.The findings for N = 2 (panel (a)) and for N = 4 (panel (b)) are very similar, only the magnitude of η is higher in the N = 4 case.When Model 1 is used, i.e. ε is taken to be 1.0, a triangular shape for η(θ ) is obtained, which approximates reasonably the simulation results, although somewhat smaller |η| values are found at the extrema of the self-bias voltage.A similar η(θ ) dependence is found when Model 2 is adopted, however, the peak amplitudes of η obtained this way are higher than those predicted by the PIC/MCC simulations.Finally, Model 3 provides a very good agreement with the simulation results.
For this latter, the different contributions to the self-bias voltage, as specified in Eq. ( 6) are also displayed in Figure 6.The contributions of the floating potentials and the bulk voltage drop prove to be small and acting against each other over the whole range of the phase angle θ .The dominant term, η w , is thus hardy distinguishable from the sum of the three terms that yields η.The above observations make us conclude that consideration of the first term only in Eq. ( 6) is sufficient, however, for the calculation of ε the more precise form of Eq. ( 10) is required.
At this point it is useful to analyze the ε values obtained from the different models, together with its values taken directly from the simulation (via Eq. ( 7)).The corresponding data are shown in Figure 7 as a function of the phase angle θ (panel (a) for N = 2 and (b) for N = 4).Additionally, the terms involved in the calculation of ε via Eq.( 10) are also displayed in Figure 7.As in the case of the self-bias voltage (see Figure 6), the most accurate model (i.e.Model 3) reproduces very well the symmetry parameter as a function of θ resulting from the simulation.The ε = 1 approximation of Model 1 is clearly a bad choice, as ε varies with θ about ±10% in the N = 2 case and more than ±20% in the case of N = 4. Model 2 (which considers only the difference of the FIG. 7. Symmetry parameter ε obtained from the PIC/MCC simulation and from the different models, as well as the terms involved in the calculation of ε from Eq. (10).Note that "ε Model 2" is equivalent to n sp /n sg , and that "ε Model 3" is computed as (Q mg /Q mp ) 2 (n sp /n sg ) (see Eq. ( 10)).Discharge conditions: Ar at p = 10 Pa, L = 2.5 cm, f 1 = 13.56MHz, V pp = 300 V. mean charge densities within the sheaths), on the other hand, largely (by a factor of ≈ 2) overestimates the variation of ε with θ .These findings confirm that the charge dynamics, represented by the term (Q mg /Q mp ) 2 plays an important role. 65he inclusion of this term in the calculation of ε ("ε Model 3" in Figure 7) yields a very good agreement with the simulation data ("ε PIC").
Next, we analyze the behavior of the discharge in the vicinity of phase angles where the self-bias voltage is (i) zero and where (ii) it is maximised.The first domain corresponds to phase angles near θ = 90 • , while for the second domain θ is around 180 • .These two domains, respectively, are shown in Figures 8(a) and (b), for various values of the number of harmonics.Here we also include data obtained with different values of the secondary electron yield and the electron reflection coefficient.2][43] These studies have, however, focused on establishing an asymmetry by using unequal coefficients at the two electrodes, while here we study the effects of these parameters with equal values at both electrodes, i.e., they are not the causes of the asymmetry, but may modify it by influencing the plasma behavior.Based on Figure 8 the following observations can be made: 1.The electron reflection has a marginal effect on η, at all values of θ , N, and γ considered.
3. At N = 2, the angle where η becomes zero changes with γ, while for N = 4 no such dependence is observed.
4. The increase of γ decreases the maximum self-bias voltage at both N values.
In the following, we provide explanations for these observations.The first of these can be explained by the fact that although reflected electrons generally increase the plasma density, at the relatively low value of R considered here this increase is small. 41We observe only a slight effect of electron reflection at γ = 0, comparison of the pairs of data sets for Regarding the second observation, it is important to realise that at θ = 90 • the maximum and minimum values of the driving voltage waveform have the same magnitudes (φ max = −φ min = φ m ), i.e. no amplitude asymmetry is present at this specific θ .According to Eq. ( 5), any deviation of η from zero can only be attributed to ε = 1, as η/φ m = −(1 − ε)/(1 + ε).(Using Eq. ( 5) instead of the more precise Eq. ( 6) is justified by our earlier conclusion that the floating potential and the bulk voltage drop act against each other in the latter equation.)Having ruled out the effect of different (positive vs. negative) voltage amplitudes, the observed asymmetry of the discharge can only originate from the specific shapes of the driving voltage waveform.Indeed, it turns out that it is the slope of the waveform, i.e., dφ (t)/dt that is responsible for the observed effect.At N = 2, as shown (by the red solid line) in Figure 9(a), the driving waveform has a long falling slope and a short rising slope, i.e. it resembles a sawtooth-down waveform. 35his shape is expected to result in higher excitation/ionization rate at the grounded side of the discharge and this is confirmed by the computed ionization rate function shown in Figure 9(b) for N = 2 and γ = 0.The theory of sawtooth waveforms 34 predicts a negative self-bias voltage, in accordance with our observation in Figure 8(a).For N = 4, we find the fastest sheath expansion around times of t/T ≈ 0.4 at the powered side, and the ionization source (see Figure 9(b)) exhibits a peak at this side of the discharge.This explains the observation of a small positive self-bias voltage (that differs from the N = 2 case).TABLE I. Self-bias voltage (η) at θ = 90 • , at N = 2 and 4, as a function of γ, the values of the symmetry parameter ε and the terms involved in ε, as well as the peak sheath voltages.The explanation of the third observation is aided by tabulated values of the DC self-bias voltage, the symmetry parameter and its relevant terms, as well as the peak sheath voltages for various values of γ, at θ = 90 • , included in Table I.The γ = 0 case was discussed above and the origin of η < 0 at N = 2 was clarified.It is important to realise that at γ > 0 the ionization balance is also influenced by secondary electrons (emitted from the electrode surfaces).The energy gain and the multiplication of these electrons, and the consequent ionization are sensitive functions of the accelerating voltage, i.e. the peak sheath voltage drop.For N = 2, we find that φ sg < | φ sp | at γ = 0. Consequently, when secondary emission sets on at γ > 0, ionization at the powered side of the discharge increases by a higher amount as compared to that at the grounded side.This is the reason why we find an increasing value of n sp /n sg (seen in Table I), which, in turn results in an increase of ε with γ.The other factor involved in ε, (Q mg /Q mp ) 2 , also increases with increasing γ, further contributing to the change from ε < 1 to ε > 1 while γ reaches 0.4.This change of ε results in a switch of the sign of η.
To understand the behavior of (Q mg /Q mp ) 2 , the time dependence of the total charge Q tot (t) in the plasma (per unit area) is plotted in Figure 10. for various values of θ , N, and γ.All the curves seen in this figure exhibit slowly decaying parts that correspond to the continuous losses of ions at the electrodes and short, rapidly increasing segments where Q tot increases because of the losses of the electrons.Losses of electrons occur during the times of sheath collapse at either side of the plasma.Taking the case of N = 2 as an example, according to Figure 9(a) the grounded sheath collapses at t g /T ≈ 0.58, whereas the powered sheath collapses at t p /T ≈ 0.9.At the time of collapse of the grounded sheath Q tot resides in the powered sheath, thus the peak of Q tot at t g in Figure 10(a) can be associated with Q mp (illustrated for the N = 2, γ = 0.4 curve).Similarly, Q tot at t p peaks at a value of Q mg .
Having understood this, we can look now at the changes of Q mp and Q mg resulting from an increase of γ.When γ is increased from zero, more ionization occurs near the powered electrode (because of the higher peak sheath voltage there, see above) and, consequently, the ion flux to the powered electrode increases more than the ion flux to the grounded electrode.As the flux of the ions and the electrons to either of the electrodes must compensate each other on time average, an increased ion flux at one electrode also requires an increased electron flux at the same electrode.Recall that the electron flux can flow only during the collapse of the sheath.Due to the higher ionization at the powered side of the plasma the electron flux at the powered electrode will get enhanced more than at the grounded electrode.Consequently, the total charge in the plasma, Q tot , will be increased more at the time of sheath collapse at the powered side, as compared to that at the grounded side.According to the explanations about the behavior of Q tot (t) provided above, Q mg will get more enhanced than Q mp when γ is increased, leading to the increase of (Q mg /Q mp ) 2 , as observed.
The scenario at the higher number of harmonics, N = 4, is somewhat different.In this case, a significantly faster sheath expansion is observed in Figure 9(a) after the sheath collapse at the powered side of the discharge as compared to the N = 2 case (compare the falling slopes of the red solid vs. chain lines, around t/T = 0).This fast sheath expansion, in contrast with the N = 2 case, gives rise to a nearly symmetrical ionization function, which is a result of a compensation between the faster sheath expansion at the powered side and a higher sheath voltage at the grounded side (see Table I).Therefore, in the N = 4 case the arguments presented above do not apply for the changes of Q mg and Q mp with γ.Moreover, the explanation presented above clearly assumed a certain degree of locality of the electron kinetics, i.e., fluxes at a given electrode were assumed to be largely determined by the ionization rate near the same electrode.However, at the conditions considered, energy gain of the electrons at one side of the discharge can contribute to the ionization and particle fluxes at the other side of the discharge as well.This property of the discharge is obviously confirmed by the ("flat top") shape of the ionization source functions plotted in Figure 9(b), especially at γ > 0. As a consequence of these effects the small counter-acting changes of the n sp /n sg and (Q mg /Q mp ) 2 on ε, with increasing γ at N = 4 are difficult to explain.This interplay of these two terms results in a nearly constant ε and DC self-bias in these cases (see Table I).Finally, we address our fourth observation, i.e., the question: why does the self-bias voltage at its extremum (near θ = 180 • ) decrease with increasing γ (see Figure 8(b))?Corresponding data for the DC self-bias voltage, the symmetry parameter and its relevant terms for various values of γ are given in Table II.At γ = 0, ionization is maintained by electron energy gain at the phase of sheath expansion.At θ = 180 • (see the thick black lines in Figure 1), the expansion of the sheath is much faster at the powered electrode and therefore the ionization rate is also higher at that side of the discharge.This results in the high n sp /n sg values seen in Table II for both N = 2 and N = 4. (Note that the value is higher for the higher N due to the faster sheath expansion induced by the steeper voltage waveform.)As a consequence of this a large η is created.
When γ is increased from zero, ionization by the secondary electrons starts to play a role.For this contribution the duration of the expanded phase of the sheaths is a key parameter.For a longer period of expansion, that is actually found at the grounded side of the discharge due to the specific ap-plied waveform, a higher enhancement of the ionization rate is expected.In accordance with this, n sp /n sg decreases with increasing γ.
Using again the argument that a higher ionization rate at the grounded side of the discharge leads to an increase of the total uncompensated charge contained within powered sheath upon the collapse of the grounded sheath, the higher increase of Q mp with respect to that of Q mg (confirmed in Figure 10(b)), and the consequent decrease of (Q mg /Q mp ) 2 can be understood.
As both terms contributing to the symmetry parameter decrease, a strong suppression of the self-bias voltage appears.At the pressure and electrode gap values considered here, energetic electrons created at one side of the plasma also contribute to ionization at the other side of the plasma, as confirmed by the shape of the ionization source functions shown in Figure 9(b).Therefore, the discharge has a tendency to become symmetrical at high γ values.

IV. SUMMARY
In this work, we have examined the establishment of a discharge asymmetry and the concomitant formation of a DC self-bias voltage in capacitively coupled RF discharges driven by multi-frequency voltage waveforms.Computations have been carried out with various values of the coefficients that characterize the electrode surfaces, i.e., (i) the coefficient of elastic electron reflection and (ii) the ion-induced secondary electron yield.The latter ranged between zero and a high value of γ = 0.4, which can characterize high-electron-yield dielectric surfaces.
The understanding of the computational results has been aided by an analytical model that is based on a voltage balance of the RF discharge.We have shown that this model, in its more complete form when it also includes the charge dynamics (by accounting for the ratios of the net charges in the two sheaths), is able to successfully reproduce and explain the behavior of the DC self-bias voltage as a function of the phase angle between the harmonics of the driving voltage waveform.
The investigations of the surface coefficients indicated that the elastic reflection of the electrons, as long as equal values are used at both electrodes, has a minor influence on the discharge asymmetry and the self-bias voltage.The secondary electron emission coefficient (for which also the same values were adopted for both electrodes) was found to influence the discharge asymmetry and the self-bias voltage in a complicated manner, depending on the phase angle and/or the number of harmonics (N).These effects were understood based on the differences of the maximum sheath voltages and the durations of the expanded phases of the sheaths at the two sides of the discharge, as well as on the charge dynamics that is influenced by the ion fluxes to the electrodes.
At our choice of the gas pressure and the electrode gap, the ionization source function was found to be non-local and a high secondary electron yield induced a tendency to restore the symmetry of the discharge at the conditions of the highest amplitude asymmetry (i.e. in the case of peaks-and valleys-type excitation waveforms).Further studies could examine such effects at conditions of lower and higher pressures and/or electrode gaps when the nonlocality of the ionization could be enhanced or suppressed.

FIG. 2 .
FIG. 2. Equivalent electrical circuit of the system investigated.The shaded area marks the plasma region, the external circuit consists of the generator G and the coupling capacitor C.
(b), for the choice of θ = 0 • .The simulation reveals that for these conditions a self-bias voltage of η ≈ −54.4 V forms (indicated by the horizontal dashed line in Figure3)(b)).The substantial contributions to the discharge voltage, which is the sum of η and the generator voltage φ (t), are the sheath voltage drops, as above.

FIG. 3 .
FIG. 3. Time dependence of the quantities involved in the voltage balance of the discharge, for single-(a) and dual-frequency (b) excitation.For N = 2, the phase angle is θ = 0 • .Other discharge conditions: f 1 = 13.56MHz, p = 10 Pa, L = 2.5 cm, R = 0, γ = 0.The driving voltage waveform is defined by eq.(1), with amplitudes given by Eq. (2), φ * = 300 V. T is the period of the fundamental RF frequency, f 1 .Note, that the bulk voltage drop is multiplied by a factor of 10.
FIG. 5. Maximum and minimum values of the applied voltage (φ max and φ min ), peak sheath voltages ( φ sp and φ sg ), floating potentials (φ f sp and φ f sg ), bulk voltage drops at the time of the maximum and minimum of the applied voltage (φ b max and φ b min ), as well as the self-bias voltage η computed from the simulation.(a) N = 2, (b) N = 4.Note that some quantities are multiplied by 10.Discharge conditions: Ar at p = 10 Pa, L = 2.5 cm, f 1 = 13.56MHz, φ * = 300 V, R = 0, γ = 0.

FIG. 8 .
FIG. 8. Self-bias voltage obtained from the PIC/MCC simulations in the vicinity of θ = 90 • (a) and θ = 180 • (b), for different number of harmonics (N), secondary electron yields (γ), and electron reflection coefficients (R).The legend in (a) also holds for panel (b).In (b), the data corresponding to N = 2 are given in the left axis, while data for N = 4 are shown in the right axis.Discharge conditions: Ar at p = 10 Pa, L = 2.5 cm, f 1 = 13.56MHz, V pp = 300 V.

FIG. 10 .
FIG. 10.Total uncompensated charge per unit area (of 1 cm 2 ) as a function of time within a fundamental RF period for (a) θ = 90 • and (b) θ = 180 • , for various values of N and γ.Parts of the curves with a slow decay correspond to the continuous losses of positive ions, while parts with a steep rise correspond to sheath collapses when electrons are lost to the electrodes.The maximum charges in the two sheaths are illustrated for one of the cases in both panels.For more explanation see text.Discharge conditions: Ar at p = 10 Pa, L = 2.5 cm, f 1 = 13.56MHz, φ * = 300 V.

TABLE II .
Self-bias voltage (η) at θ = 180 • , at N = 2 and 4, as a function of γ, as well as the values of the symmetry parameter ε and the terms involved in ε.