The form of electron energy probability function (EEPF) in low pressure radio frequency (RF) bounded plasma has been a subject of considerable interest [1–7]. Godyak et al. [8] looked at EEPF of inductively coupled plasma (ICP) in Argon over a wide range of pressure and power combinations. Meige and Boswell [9] looked at EEPF in an ICP to investigate the origin of the break energy at the plasma space potential. Bhattacharjee et al. [10] looked at the effect of excitation frequency. In this article we investigate the effect of varying the length and pressure of the discharge. Our aim is to see the effect of these important parameters on the EEPF, which dictates the overall discharge properties such as ionization and other rate constants.

The present research is based on an earlier work by Bhattacharjee et al. [10], where the effect of the excitation frequency on EEPF was investigated. Here we investigate the effect of the axial length (L) and the pressure (p) on EEPF of the same system, and show that any change in one of these parameters while keeping the other fixed, produces quantitatively close changes in the discharge parameters—implying that these two parameters are “self-similar” for low energy electrons.

The motivation further arises from the fact that, normally in most wave assisted plasma devices, the heating frequency is maintained constant (e.g., 13.56 MHz for rf discharges and 2.45 GHz in microwave discharges [11–13]), and two important parameters that may be varied are the system length (size) and the pressure, which can determine the dynamics of the particles (electrons and ions) present in the discharge. In many actual experiments, it may be inconvenient and quite cumbersome to vary the length of the system. In particular, it is desirable to know whether instead of varying the length of the discharge, is it possible to achieve nearly equivalent effect (such as in the nature of EEPF) by varying the pressure of the discharge? This would be quite helpful in applications such as plasma processing or plasma sources for ion beam extraction which require optimization of discharge conditions [14, 15] to achieve a certain profile of the EEPFs.

On another viewpoint, Paschen's law shows that the discharge length and pressure are critical parameters (actually the product pL) that determine the breakdown threshold of a gas in a dc discharge. Lisovskiy et al. [16] derived a similar law for RF discharges. The interplay between the parameters L and p is important to determine the breakdown threshold for RF discharges. Therefore, the effect of L and p on the discharge, particularly on the dynamics of the electrons is the subject of the research presented in this article. The RF heating frequency is kept fixed at 13.56 MHz and the length and pressure has been varied over a wide range i.e., 0.1–0.5 m and 1–10 mTorr, respectively.

The article is organized as follows. In section II, the PIC simulation scheme is presented. The simulation results from the PIC code are presented in section III. The simulation results are compared with an analytical model in section IV. Finally, the results are summarized and conclusions from the research drawn in section V.

In the following, we discuss the PIC simulation scheme. The plasma is confined between two infinite parallel plates separated by a distance L along the spatial x direction, which is the discharge length. L is varied to investigate its effect on EEPF. Initially the neutral gas pressure is maintained at 1 mTorr, this is varied later to see its effects on EEPF. Argon has been used as the test gas. The simulation is based upon a well-known PIC scheme [17–20] with Monte Carlo Collisions [21]. The simulation is one dimension (1D) in space (x) and three dimensional (3D) in velocities (v_x, v_y, and v_z). Accurate electron-neutral and ion-neutral collision cross sections are employed to ensure realistic simulations; data can be found in References Tachibana [22] and Pack et al. [23] for electron neutral collisions (elastic, exciting, and ionizing) and in Reference Phelps [24] for ion-neutral collisions (elastic and charge exchange). Coulomb collisions are not included in the present case as their effect would be negligible in the present low density condition.

For the present investigation, the electrons are heated with a scheme intended to model “inductive” excitation similar to that described by Turner [25]. Here, the RF heating field is localized at the center of the discharge over a length 0.05 m for all situations of different L and p. The shape of the RF electric field along the x axis is a top hat and the amplitude is typically ~2.8 × 10⁴ V/m in each direction y and z, which corresponds to a heating current density of 5 A/m². The RF electric field is applied in the y and z directions perpendicular to the spatial x direction where the boundaries are located. This allows us to heat the electrons both in y and z directions, momentum and energy being transferred to the other directions via electron-neutral collisions. The amplitude of the RF electric field is taken as uniform over the length of the heating region and localized within it. In the present investigation, the heating scheme is very much similar to the previous work by Bhattacharjee et al. [10]. The schematic of the heating field has been shown in Figure 1.

The simulation is allowed to run for several thousands of RF cycles in order for the system to reach equilibrium. The initial temperatures of the ions and electrons are taken as 0.026 eV (room temperature) and 3 eV, respectively. It is found that the typical time required to reach equilibrium is ~40 μs after the initiation of the discharge. Each simulation particle is a “macroparticle” representing 2 × 10⁸ actual particles. The number of cells along the x axis is 250 and the time step is ~1 × 10⁻¹⁰ s, which satisfies the stability and accuracy criteria for PIC scheme Δt ≪ 0.2/ω_p and Δx ≪ λ_D where Δt and Δx are the temporal and spatial increments used in the simulation and ω _p and λ_D are the plasma frequency and the Debye length respectively. The typical parameter values are plasma density ~2.5 × 10¹⁴ m⁻³, electron temperature T_e ~ 5 eV, λ_D ~ 10⁻³ m, and ω_p ~ 8.94 × 10⁸ rad/s.

The EEPFs are presented for three different velocity components v_x, v_y, and v_z and obtained at the center of the discharge. The EEPFs are measured along the abscissa and ordinate is in natural log scale so that a Maxwellian distribution yields a straight line [26]. The electron energies have a resolution 2 × 10⁻³ eV and these have been acquired after the system is well deep into the equilibrium (typically 100 μs) and run for additional number of ~5 × 10⁵ time steps (~50 μs) and averaged. The system is finite in the x direction and infinite in y and z directions.

Keeping the heating frequency and the discharge pressure fixed at 13.56 MHz and 1 mTorr, respectively, the length of the one dimensional plasma system has been varied from 0.1 to 0.5 m. For varying pressure, the length of the system has been fixed at 0.1 m and the pressure varied in the range 1–10 mTorr. The plasma space potential ϕ_p has been first looked at for different L and p within the above range. Figures 2A,B represent ϕ_p under the above conditions.

The formation of sheath near the boundaries can be clearly seen from the profile of ϕ_p (Figures 2A,B). The values of ϕ_p at the center (x/L = 0.5) have been plotted against L and p in Figures 3A,B, respectively.

We note a decreasing trend of ϕ_p at the center of the discharge with increase in both L and p.

Next we investigate the variation of EEPF in the x (bounded) direction at different L and p at the center of the discharge which are presented in Figure 4.

It is observed from Figure 4 that the typical structure of the EEPFs in the x direction is two Maxwellian with a step and this feature is common for all L and p [10]. The break energy tracks the value of ϕ_p at the center of the plasma [9]. The EEPFs comprise of low energy trapped electrons whose energies are below the break energy (i.e., below ϕ_p) and high energy escaping electrons whose energies are above the break energy (i.e., above ϕ_p). The break energy is shown by arrows in both the figures. With increase in L, the population of low energy trapped electrons as well as high energy escaping electrons increase (cf. Figure 4A). With increase in p, the population of low energy trapped electrons also increases as was observed with increasing length, but the population of the high energy escaping electrons decrease. This is expected, because increase in pressure increases the collisionality of the plasma and the high energy electrons lose energy by transferring energy to the neutral atoms. From the slope of the curves, the “temperature” of the two types of electron populations have been evaluated. It has been found that the “temperature” of the low energy trapped electrons is higher than the “temperature” of the escaping high energy electrons as plotted against L and p in Figure 5. This typical structure of the EEPF has also been reported in experiments by Godyak et al. [8].

The variation of T_e of both trapped and escaping electrons in a magnetically expanding current free double layer (DL) plasma was also investigated by Takahashi et al. [27] The same structure of EEPF and the same behavior of T_e has been reported in that article too (Figure 3 of reference 27).

Next we will discuss the variation of EEPF in y and z (unbounded, where there is no sheath) directions for different values of L. All the measurements have been taken at the center of the discharge.

We present a comparative study of EEPF's in y direction only (Figure 6) as it has the same shape as that of the z direction [10]. From Figure 6, it is observed that for all L and p under consideration, the nature of the EEPF is more linear (in the energy range 0–30 eV) compared to that of the EEPF in the x direction where the EEPFs are a two Maxwellian with a step. However, a tail region populated by high energy (“hot”) electrons (with energies above 30 eV) is observed which from now on will be referred to as a “hot tail.” The break energy is not as demarcated as in x direction.

From Figure 6A, it can be observed that, with an increase in L, the population of both low energy bulk electrons and the high energy tail electrons increases indicating that the heating effect is greater at larger lengths. From Figure 6B, we observe that with increase in pressure p, the population of the low energy electrons increase in a similar manner with varying length, however, the hot tail diminishes with increase in pressure as was seen earlier for the case of EEPF in the x direction. The electron temperature in the two parts of the distribution (low energy electrons and the hot tail) has been found from the slope of the curve and is plotted in Figure 7.

We observe that with increase in L and p, the temperature of the high energy electrons in the tail region [labeled as (T_e)_high in Figure 7] increases quite rapidly. But the temperature of the low energy electrons [labeled as (T_e)_low in Figure 7] shows a slight decrease with increase in both L and p. (T_e)_high tends toward (T_e)_low at smaller lengths and pressure, indicating that the distribution becomes more and more Maxwellian with diminished hot tail at smaller L and p.

The discharge parameters behave in equivalent manner with changes in L and p. ϕ_p is less for both higher length and pressure (Figures 1, 2). This is evident because with increase in both L and p, the temperature of high and low energy electrons decreases. The observed effects on the EEPFs for both varying L and p have been verified by writing the particle balance equation for cylindrical plasmas [28] for our one dimensional (1-D) plasma system. As there is no radial loss, we neglect the radial terms.

For our 1-D system, we equate the total surface particle loss to total volume ionization as follows:

where K_iz, n_g, L, and u_B are ionization rate constant, gas density, plasma length and Bohm velocity respectively. h_L is the axial edge to center density ratio [28] with

where, λ_i = ion mean free path given by λi=1ngσi, σ_i is the ion-atom scattering cross section ~10⁻¹⁴ cm² (for low energy ions, 0.026 eV). The ideal gas equation is:

where p is the neutral gas pressure, n_g is the neutral gas density, T is the temperature of the gas taken to be 25°C or 298 K in the simulations.

Using this, λi=1330p cm, where p is in Torr.

From Equation (1) we obtain,

The Bohm velocity

T_e is the electron temperature, M is the ionic mass: 39.948 a.u. for Argon. Equation (4) can be rewritten as,

where deff=L2hL is the effective plasma size.

For ionization, K_iz(T_e) = < σ_iz (v) v >_v, where σ_iz is the ionization cross section. The average is done over the distribution found from the PIC simulations. Hence,

where f(v) is normalized as,

For ionization, we expand the Thomson cross section [σ_iz (v)] around the ionization threshold energy E_iz to obtain,

where σ₀ = (e/4πϵ₀E_iz)² and E=12mv2e. With e and m being the elctron charge and mass respectively.

Using the form of ionization cross section of Equation (9), K_iz has been evaluated from Equation (7) in terms of energy. For calculating u_B from Equation (5), we have used the weighted average of the electron temperature T_e of the two population (trapped and escaping) of electrons, found from the PIC EEPFs. After calculating K_iz and u_B, the quantity n_gd_eff has been evaluated from Equation (6).

For analytical calculations, the variation with pressure of the neutral gas density n_g has been found from the ideal gas equation (Equation 3). We can obtain d_eff from the equation, deff=L2hL where h_L can be found from Equation (2). For both varying L and p, we have calculated the product ngdeff=pL2KBThL. The analytical value of n_gd_eff is plotted with the values obtained from the simulations (Equation 6) (calculated from K_iz and u_B) against L and p in Figures 8A,B, respectively for varying L and p.

It is therefore seen that the simulated and the analytical results are nearly the same for both varying L and p.

Figure 8 shows an agreement of the simulated results with a global analytical model. Equation (4), obtained from the particle balance Equation (1) may be further explored for combining data obtained at different L and p. We rewrite Equation (4) as follows:

where h_L been taken from Equation (2) and n_g from the ideal gas equation (Equation 3) and λi=1330p cm =10−2330p m.

We look at the behavior of the left hand side of Equation (10) with the product pL for different p's and L's under consideration. Figure 9 presents the behavior. In Figure 9, KizuBkT in Equation (10) has been plotted against pLpLmax where pL_max is the maximum value of pL for the either work conditions. The product pL is thus normalized to match the range of pL values for varying p and L.

From Figure 9, we note that the behavior of KizuBkT in Equation (10) follows the same trend with varying L and p. The data sets for varying p and varying L while keeping the other parameter fixed, have been indicated by red circles and black squares respectively in Figure 9. Thus, even if we have used the same sets of values of the product pL while varying L and p keeping other fixed, the absolute value of the quantity KizuBkT is different for the two data sets because it has been obtained from the simulated EEPFs that differs for the two data sets. But both the variations can be fitted with a single universal curve (shown by blue curve in Figure 9) obtained from the right hand side of Equation 10. Thus, it behaves in a qualitatively similar manner with the product pL.

In the light of the above, we now look at the variation of the plasma parameters i.e., the space potential, and the electron temperature as a function of α (normalized pL).

Figure 10 shows such variation of ϕ_P.

We observe a qualitative agreement of the two variations as predicted by Equation 10 through Figure 9. Now we look at the variation of T_e for both low energy electrons found from EEPF in x and y direction.

From Figure 11 we note that the temperatures of low energy electrons (found both from EEPF in the x and y direction) behave in a qualitatively similar manner with varying L and p. However, we recall our observations from Figures 4, 6 that with increasing L and p, the high energy part of the EEPFs (population of high energy escaping electrons in EEPFx and the hot tail in EEPFy) do not behave in a similar manner. So by observation, we can say that the parameters p and L are approximately similarity parameters for the low energy electrons trapped in between the sheaths.

Finally, we have looked at the plasma electron density obtained by integrating the EEPFs over all velocities or energies. Figure 12 shows the results.

Both the results follow the same trend and is consistent with the behavior of electron temperature.

The similarity in the variation of T_e with both L and p essentially arises from the collisionality of the electrons with the neutral particles.

The changes in the electron temperature with L and p, determine the modification of the electron energy distributions with L and p as presented in the manuscript. It is seen directly from Figures 5, 11, and indirectly from Figures 9, 10 that the electron temperature increases as the system gets shortened or the pressure is decreased. This can be explained by the ratio of mean free path to the system length (λ /L: λ being the mean free path and L be the system length). λ is inversely proportional to n_g, the neutral gas density. From the ideal gas equation (Equation 3), we see that for a system with lower pressure, n_g is lower and this implies a larger λ. So with decrease in pressure, λ approaches the system length L and the ratio λ /L approaches one. This decreases the probability of collisions, and the electrons suffer less number of collisions and less energy loss happens. Therefore, the electron temperature is higher at a lower pressure. In a similar manner, if the pressure is kept fixed and the system is shortened, λ approaches L, and similarly, the electrons suffer fewer collisions leading to lesser loss of energy. A consequence of this is rise in electron temperature. This further confirms the role of similarity of the two parameters L and p. In summary in both cases, the ratio λ /L changes in a similar manner and thus changes the electron temperature which causes modification in the EEPFs.

The shape of the distributions in the two directions [bounded: x, and unbounded: y (or z)] may be understood as follows. The results indicate that for different L and p under consideration, the EEPFs in the x direction has a shape of a bi-Maxwellian with a step and that in the y (or z) direction it is a non-Maxwellian with a hot tail. The sheaths are formed only in the x direction. So, along the x direction, the EEPFs comprise primarily of two parts: the first part corresponds to electrons having energies less than the plasma potential which remains trapped in the bulk plasma. And the second part consists of electrons having energies higher than ϕ_p, most of which escapes. As the EEPFs have been obtained at the center of the discharge, the population of the electrons corresponding to the second part are quite small compared to the trapped ones of the first part. In the y and z directions, the system is unbounded. No sheaths are formed and therefore ϕ_p has no role in electron trapping and no break energy is observed. Therefore, the nature of the distribution changes.

We summarize the results obtained from the PIC simulations as follows.