Effect of Electron Energy Distribution on the Hysteresis of Plasma Discharge: Theory, Experiment, and Modeling

Hysteresis, which is the history dependence of physical systems, is one of the most important topics in physics. Interestingly, bi-stability of plasma with a huge hysteresis loop has been observed in inductive plasma discharges. Despite long plasma research, how this plasma hysteresis occurs remains an unresolved question in plasma physics. Here, we report theory, experiment, and modeling of the hysteresis. It was found experimentally and theoretically that evolution of the electron energy distribution (EED) makes a strong plasma hysteresis. In Ramsauer and non-Ramsauer gas experiments, it was revealed that the plasma hysteresis is observed only at high pressure Ramsauer gas where the EED deviates considerably from a Maxwellian shape. This hysteresis was presented in the plasma balance model where the EED is considered. Because electrons in plasmas are usually not in a thermal equilibrium, this EED-effect can be regarded as a universal phenomenon in plasma physics.

in plasma physics. Therefore, a full quantitative experiment and improved model are required to provide a universal solution to plasma hysteresis.
Here, we present theory, experiment, and modeling about origin of the hysteresis phenomena. It is found that, as a new theory, evolution of the electron energy distribution function (EEDF) creates hysteresis in the plasmas. Significantly, this EEDF-effect on the hysteresis can be generalized to a universal phenomenon in gas discharge plasmas. This is because electrons are not in a thermal equilibrium in most plasma discharges, and EEDFs deviate considerably from Maxwellian shape due to wave-electron resonant interaction like Landau damping [25][26][27][28][29][30][31] and energy-dependent electron heating by the Ramsauer-Townsend scattering minimum as quantum-mechanical effect [32][33][34][35][36][37] .

Experimental results in Ramsauer and non-Ramsauer gas discharges
The experiment is done in an inductively coupled plasma with a planar type antenna coil. The radio-frequency (13.56 MHz) power is applied to the two-turn antenna coil through automatically controlled impedance matching network. To measure the plasma power, which indicates the transferred power to plasma excluding the matching network and antenna coil, the current monitoring of the antenna coil and measurement of the system resistance are performed. To measure plasma density and EEDF, a single Langmuir probe is used. Detailed descriptions for the experimental setup and the plasma power measurement can be seen in the Method section. To see the plasma mode transition and the hysteresis, input power increases or decreases on 1 W basis step by step, and following results are obtained.
When the input power or the plasma power increases and decreases at Ar gas of 40 mTorr, the plasma density shows an identical trajectory and thus, the plasma has no hysteresis (Fig. 1a). However, a huge hysteresis loop of the plasma density is present at Ar gas of 250 mTorr (Fig. 1b). When the input power increases from 26 W to 27 W at the high gas pressure of Ar, an abrupt increase in the plasma density is observed, indicating the E-H mode transition ( Supplementary Fig. 1 When the input power or the plasma power decreases, the plasma density still remains at the H mode until plasma power of 12.56 W at input power of 23 W. A slight decrease in the input power from 23 W to 22 W results in an abrupt change of the discharge characteristic to an E mode with low plasma density. Through this transition of the discharge mode, there is a huge hysteresis with the trajectory [I→ II→ III→ IV] in case of the Ar gas of 250 mTorr (Fig. 1b). In He gas discharges, it is noted that the hysteresis is not observed at either a high gas pressure of 300 mTorr (Fig. 2c) or a low gas pressure of 40 mTorr (not presented here). This direct experimental evidence implies that the plasma has bi-stability characteristic with strong nonlinearity only at the high pressure Ramsauer gas discharge.
Why then does this hysteresis occur at the high pressure Ramsauer gas? A crucially different plasma property on the discharges (Fig. 1) is the shape of the electron energy probability function (EEPF), which is caused by the quantum-mechanical effect of the Ar gas, called the Ramsauer-Townsend scattering minimum. The Ar (Ramsauer gas) has a minimum cross section for electron-neutral collisions at low electron energy range below 1 eV, while the He (non-Ramsauer gas) has no minimum cross section in this electron energy range as indicated in Fig. 2e 38,39 . Due to the collision cross section of Ar, the electron heating rate and EEPF strongly depend on the electron energy at the high pressure experiment, as follows.
In the low plasma density mode (E mode) with an Ar gas pressure of 250 mTorr where the electron-neutral collision frequency v m is much higher than the driving angular frequency, the Ohmic power P ohm transferred to electrons per unit volume is inversely proportional to v m at given E-field as P ohm ≈ e 2 E 2 /(v m 2m). Thus, the low energy electrons are more efficiently heated than the high energy electrons at the Ramsauer gas, and the EEPF shows Druyvesteyn distribution rather than Maxwellian distribution at the high gas pressure 32 (Fig. 2a). On the other hand, the EEPF transits into Maxwellian distribution as shown in Fig. 2b due to the heating mode transition and the enhancement of the electron-electron collisions in the high plasma density (H) mode. This occurs because the electrons are accelerated by a circular induced electric field in the skin depth, and the confinement of electrons is enhanced in the discharge mode.
This evolution of the EEPF with the heating mode transition can exhibit hysteresis in plasma physics, as follows. The electron heating and electron-electron collisions replenish high energy electrons, which were originally depopulated by inelastic collisions at the high gas pressure. This enables the replenished electrons to participate in ionizations, causing collisional energy loss ε c to create an ion-electron pair to be strongly reduced. Therefore, ionization efficiency can be dramatically increased with the evolution of the EEPF during the E-H mode transition. With decreasing plasma powers, the EEPF still remains a Maxwellian distribution at a certain threshold regime such as region (III) in Fig. 1b. In the regime with Maxwellian EEPF, the ionization efficiency is sufficient to produce the H mode plasma, and thus H-E transition occurs at lower plasma power than that in E-H transition. Therefore, the hysteresis becomes possible due to the change of the EEPF. Note that the EEPFs on both E and H modes have a Maxwellian distribution at the high He pressure of 300 mTorr (non-Ramsauer gas) and at the low Ar pressure of 40 mTorr (Supplementary Fig. 2,3). Hysteresis was not observed in these conditions (Fig. 1a,c).

Discussion
To evaluate the effect of the EEPF on the hysteresis, theoretical calculations of the ε c and the stable plasma density were investigated by using stepwise global model considering EEPF for the first time. In the stepwise model, we consider argon atom energy levels of resonance states (4s r ), metastable states (4s m ), and 4p excited states Fig. 5. To obtain the ε c , power and particle balance equations (eqs. (6) and (7)) are used. In the theoretical calculation, EEPF should be considered because the reaction rate constant is strongly dependent on the shape of the EEPF. The rate constants are obtained by using collision    Figure 3 shows the ε c curves obtained from self-consistent calculation using equations (2), (6) and (7) in the Method section. When only Maxwellian EEDF is considered, ε c,Maxw decreases slightly from 71.2 V to 60.1 V, and such a minor variation in ε c,Maxw is not enough to explain the density jump and the hysteresis of Fig. 1b. However, when the EEDF-modification is included in the stepwise balance model, ε c changes significantly and thus, the hysteresis can be clearly explained. At the region (I) in the E mode (Fig. 1b), the ε c,Druy is 165.2 V, while ε c,Maxw showed a remarkable reduction to 61 V at the region (II) with the E-H mode transition. With decreasing plasma power, the high plasma density is still maintained (ε c,Maxw : 61.6 V). With further decreasing the plasma power, ε c increased significantly to 148.5 V on the H-E mode transition with the modification of the EEDF. From these variations of the EEDF, ε c has a huge hysteresis loop in Fig. 3, which shows a consistent result with the plasma density-hysteresis of Fig. 1b.
To precisely show the EEDF-effect on the plasma-hysteresis, stable plasma density was obtained theoretically. The plasma density is determined by the balance between transferred power P trans and power dissipation P diss . Using equations (2-7) in the Method section, we can obtain the P diss curves for each EEDF (Fig. 4). The P trans curves were calculated from equation (8). For the stable working point of the plasma density, the intersection point between P trans and P diss should satisfy ∂P trans /∂n e < ∂P diss /∂n e . Figure 4 shows the stable points (■: evolution of the EEPF considered, □ : only Maxwellian EEPF considered).
If the EEDF effect is not considered, and a conventional analysis with the assumption of Maxwellian EEDF is only used in the stepwise model, then there is little hysteresis during the E-H and H-E transition: trajectory B [a→ b→ c→ II→ III→ d→ e]. When the evolution of the EEPF is considered, however, the hysteresis is clearly found. At the E mode, a stable working point is built between P trans and P diss,Druy , while stable working point is made between P trans and P diss,Max at the H mode. Until a certain point (region III), stable working point is built between P trans and P diss,Max . After reduction of the plasma power below a threshold value, the discharge is transited into E mode, and stable point is presented between P trans and P diss,Druy (region IV). From these discharge dynamics with evolution of the EEPF, trajectory A (I→ II→ III→ IV) is made and thus, the hysteresis loop is definitely solved in Fig. 4. Therefore, our theoretical approach for the plasma balance equation considering the EEDF gives a result consistent with the hysteresis experiment in Fig. 1b.
Our experiment demonstrates that the plasma hysteresis occurs at high pressure Ramsauer gas discharge where the evolution of the EED occurs with the heating mode transition. The plasma balance modeling clearly explains that the EED-modification creates a strong plasma nonlinearity and bi-stability of plasma. Based on our experiment, theory, and modeling, we conclude that the hysteresis is caused by the evolution of the EEPF. Because electrons in plasmas are usually not in a thermal equilibrium, this EED-effect can be regarded as a universal phenomenon in plasma physics. This research is of prime interest to the plasma physics and application community as it contributes to study of plasmas for fundamental physics, industrial nano-processing plasma, fusion reactor, and bio-medical plasma. The suggested theory in this study will aid the resolution of incompletely solved plasma phenomenon, as well as the plasma hysteresis. Method Experimental device description. The experiment was performed in conventional planar type inductively coupled plasma (ICP) reactor with a cylindrical shape of an inner diameter of 26 cm and a discharge gap of 12 cm. RF power of 13.56 MHz was applied to two-turn antenna coil through automatically controlled impedance matching network, and the reflection of the input power P in was under 1%. The plasma power P p , which implies the transferred power to plasma, was obtained using current I measurement of the antenna coil. When the chamber is in base pressure under 3 × 10 −6 Torr, plasma is not ignited. In this case, the system resistance R s including antenna coil and matching network can be found as R s = P in /I 2 where I is the coil current. When Ar or He gas is inserted and gas pressure is adjusted to the experimental condition, the plasma is turned on with applying RF power. With the plasma sustainment, the P in is consumed to two components (system power loss P s = I 2 R s and P p ), and the P p can be obtained using P p = P in − I 2 R s . EEPF and plasma density measurements. The absolute values of the EEPF are obtained by using the AC superposition method 40,41 . Since second derivatives (d 2 I p /dV 2 ) of the I-V probe current is proportional to the EEPF (g p (ε) = ε −1/2 g e (ε)) in isotropic plasmas 32 , the second harmonic current I 2ω measurement in the AC superposition method gives the EEPF g p (ε) and the EEDF g e (ε) as Theory details for plasma stability and hysteresis. The collisional energy loss ε c is obtained from the stepwise plasma balance equation. In this work, argon atom energy levels of resonance states (4s r ), metastable states (4s m ), and 4p excited states 42,43 are considered, as illustrated in following Table. 1 and Fig. 5. Then, the ε c is given from where K iz,i , K ex,j , and K el,k , are the rate constants for the ionizations, excitations, and elastic collisions, and the last term of equation (2) is reduced as K el n g n e (3 m/M)T e 42,43 . In equation (2), the EEDF should be considered because rate constant K is strongly dependent on the EEDF shape. The rate constant for the collision reaction between a and b is defined as where σ is the collision cross section for each reaction process in Table 1, and the collision cross sections data for the Ar gas are referred to in refs 38,44-49. The general EEDF can be written as follows 50 :  (3) and (4)  where A eff , ε i , ε e , P diss , and P trans are the loss area as A eff = 2πR[0.86R(3 + h/2λ i ) −1/2 + 0.8h(4 + R/2λ i ) −1/2 ] 18 , the ion energy loss, the electron kinetic energy loss, the dissipated power, and the total transferred power, respectively. From equations (2-7), P diss curves for each EEDF are obtained (P diss,Druy curves for the Druyvesteyn distribution and P diss,Max curves for the Maxwellian distribution). In the inductive discharges, plasma power is transferred by two components: E coupling (P cap ) and H coupling (P ind ). From the Maxwell equation and integrating the Poynting vector over the surface area of the interface of the plasma and dielectric wall, P trans in the plasma is presented as 51  where α = ((ω pe /c) 2 /(− 1 + iv m /ω)) 1/2 , J 0 and J 1 are the zeros and first order Bessel functions, V c /l is the voltage difference, λ = κ p 1/2 ω/c where κ p = 1− ((ω pe /c) 2 /(1− iv m /ω)), and λ′ = κ d 1/2 ω/c. Finally, the plasma density is determined by the balance equation between the P trans and the P diss (eqs. (7) and (8)).