Hysteresis between gas breakdown and plasma discharge

Hysteresis occurs when a system undergoes different processes depending on the direction of the mode transition. In low-temperature plasmas, hysteresis has been widely observed in direct-current (DC) discharges,^1–4 radio frequency (RF) inductive plasmas,^5–8 and microwave plasmas.⁹ Understanding the stability and mode transition of plasma discharge plays an important role in controlling the plasma behavior in engineering applications, such as microelectronics processing and space propulsion.

In DC discharges, hysteresis is mainly observed between the Townsend (breakdown) and glow discharge^10,11 and during the glow-to-arc transition.¹¹ Recent study also shows the transition between stationary and oscillatory modes in DC glow discharges.¹² Gas breakdown is a phenomenon in which electrically insulating gas becomes conductive by generating a self-sustaining plasma. DC breakdown is typically characterized by the breakdown voltage, V_b, as a function of pd, where p is the pressure and d is the distance of the anode–cathode gap, and the ion-induced secondary electron emission (cf. Paschen theory¹³). For the Townsend-to-glow transition, self-sustaining discharge is ignited at the breakdown voltage (cf. Townsend discharge) and a normal glow discharge can be sustained at a lower voltage than the breakdown voltage. The current–voltage characteristics shows that the plasma resistance decreases during the Townsend-to-glow discharge, resulting in a negative differential resistance,^14–17 i.e., $d V d / d I d < 0$⁠, where V_d is the discharge voltage between the anode and cathode gap and I_d is the discharge current. Some studies suggest the existence of multiple solutions assuming steady-state conditions, often called bifurcation, in DC plasma systems.^3,18–20

In actual experiments, typically a ballast resistor is used to limit the current drawn from the discharge by varying the discharge voltage. Hence, the discharge voltage can be given as $V d = V app − R I d$⁠, where V_app is the applied voltage and R is the resistance of the circuit. The applied voltage is gradually increased, i.e., during the forward voltage sweep, until breakdown is observed. Several fluid simulations are used to obtain the self-sustaining conditions using the forward sweep.^21–23 Once the plasma is formed, increasing the applied voltage leads to an increase in the discharge current I_d and thus a decrease in the discharge voltage V_d (cf. glow discharge). In Ref. 24, a backward sweep of the voltage was performed from the glow discharge regime. By decreasing the applied voltage, the current decreases until the plasma is extinguished. For a DC discharge with a large gap, the current–voltage characteristics between the forward and backward sweep are different, thereby generating a hysteresis. The voltage at which the plasma extinguishes, here we call the plasma sustention voltage, is lower than the breakdown voltage. Additionally, Ref. 24 reports the hysteresis between breakdown and plasma discharge is mitigated in a microgap where field emission from the cathode is dominant. While computational models have been used to study such hysteresis using different electrical resistances,²⁵ the detailed mechanisms of the hysteresis and the fundamental plasma characteristics are still not well understood.

Paschen curve for gas breakdown is typically obtained in the limit of a small initial space charge by searching for the voltage at which self-sustaining plasma is observed. When an initial seed plasma is placed and if the voltage is below the breakdown voltage, i.e., V_d < V_b, the plasma decays. On the other hand, when the voltage is above the breakdown voltage, i.e., V_d > V_b, the seed plasma starts to grow in time. Thus, there is a specific voltage at which the plasma density stays time-invariant. This is the condition for a self-sustaining plasma and the voltage required is called the breakdown voltage, i.e., V_d = V_b. A variety of computational models have been developed to predict the breakdown characteristics in both DC and RF discharges.^26–30 Recently, a fluid moment approach was benchmarked against a particle-in-cell Monte Carlo collision (PIC-MCC) model and drift-diffusion (DD) model, demonstrating the importance of electron transport within the DC discharge.³¹ 

In this paper, a one-dimensional PIC-MCC model is used to study the mechanisms behind the hysteresis between the Townsend and glow regimes. We performed two kinds of simulations using the PIC-MCC model. First, a ballast resistor is considered for the PIC-MCC model to demonstrate forward and backward voltage sweeps. The forward voltage sweep is performed from the breakdown voltage in the absence of field emission. The simulation results suggest that transition from the Townsend regime to a steady-state glow regime is observed, resulting in a negative differential resistance, which is consistent with the experimental observations.^14,24 Then, backward voltage sweep is performed from the steady-state glow discharge. The discharge voltage obtained from the backward voltage sweep can be smaller than the breakdown voltage, exhibiting a hysteresis similar to experiments.^14,24 Second, a finite initial plasma density (⁠ $n p 0$⁠) is assumed within the anode–cathode gap to study the plasma characteristics during the backward voltage sweep (without a ballast resistor). In the limit of low plasma density, V_s = V_b is obtained. The simulation results suggest that V_s decreases with an increased initial plasma density and saturates above a certain initial plasma density. The minimum sustaining conditions coincide with the discharge voltage obtained from the backward sweep, demonstrating that the minimum sustaining results are unique solutions.

The PIC-MCC model is described in Sec. II. Section III discusses the hysteresis obtained in the current–voltage characteristic during a forward and backward applied voltage sweep for the small pd case. Section IV shows the Paschen curves for the gas breakdown as well as the minimum sustaining conditions in the presence of an initial plasma density. In Sec. V, the steady-state plasma profiles for low and high pd regimes are analyzed to investigate the mechanisms of the hysteresis between breakdown and plasma discharges. Finally, Sec. VI shows that the field emission can mitigate such hysteresis.

In this study, an explicit, electrostatic, particle-in-cell (PIC) simulation coupled with Monte Carlo collisions (MCC)^32,33 is developed for ions and electrons to investigate the hysteresis.

Both ions and electrons are resolved in one dimension of physical space and in three dimensions of velocity space (1D3V). The leapfrog algorithm is employed to update the particle position and velocity. Argon is considered as the working gas. The electron-neutral collisions are accounted for, whereas ion-neutral collisions are not considered. Null collision algorithm is employed for the simulation of electron-neutral collisions.^34,35 The Biagi cross-sectional database is used for electron-argon collision.³⁶ The database is the collision between electron and the ground-state neutral including singly charged ionization, first three electronic excitation, and elastic collision. The neutral atom density is assumed to be constant, as the plasma discharge is weakly ionized. It is to be noted that fast neutrals, metastables, and photoemission effects are found to be critical to validate theory and experiments in Ref. 37. Reference 38 also reports that the fast heavy species play an important role in the ionization and excitation process in the abnormal regime. However, these effects are not considered for simplicity in the present study. Also, stepwise ionization from metastable states is not considered in the PIC-MCC model, which might be important for the low-temperature case (i.e., high pd).^39,40 Inclusion of these effects is reserved for future work.

The ions and electrons are absorbed when colliding with the anode or cathode, contributing to the current. Ion-induced secondary electron emission (SEE) is accounted for, i.e., when an ion hits an electrode, an electron is emitted with a certain probability, namely, the ion SEE coefficient, γ_i. In this paper, constant $γ i = 0.1$ is assumed for simplicity, while it has been reported that γ_i has ion energy dependency.³⁷ 

The PIC-MCC simulation code is parallelized using message passing interface (MPI). The PIC-MCC simulation code is benchmarked against the Turner test case based on 1D capacitively coupled plasma (CCP)⁴⁶ and 1D cross field discharge test case.³³ Also, the Paschen curve for the gas breakdown condition under $d = 1 cm$ is also benchmarked with a full-fluid moment model and drift-diffusion (DD) model.³¹ 

Based on the experimental setup in Ref. 24, we assume that the neutral gas is argon at atmosphere pressure, i.e., p = 1 atm, and room temperature of 300 K, resulting in $n g ≈ 2.5 × 10 25 m − 3$⁠, where n_g is the neutral number density. The distance between the anode and cathode d is varied $d ∈ [ 1 μ m , 100 μ m ]$ to obtain the self-sustaining voltage as a function of pd. The cathode and anode are located at $x / d = 0$ and $x / d = 1$⁠, respectively. It is important to note that the self-sustaining condition (cf. Paschen curve) is a function of pd, and it is confirmed that by changing the pressure to a smaller value, hence changing d to a larger value, generate a self-similar solution under the influence of same physical processes. In addition, streamers are not present in, or not the focus of, this work. The streamer mechanism becomes important when pd is the order of $10 3 Pa · m$ (e.g., Meek condition),^10,47 while this paper focus on $p d = 0.1 – 10 Pa · m$⁠. Here, note that 1 Pa · m is approximately 0.75 Torr · cm. A streamer might transition to the arc discharge due to the heating of neutral gas and electrodes,⁴⁸ and the glow-to-arc transition⁴⁹ is not observed in this work.

Dirichlet boundary conditions are employed for the anode and cathode. $ϕ = V d$ is fixed at the anode, i.e., $x / d = 1$⁠, while $ϕ = 0$ V is fixed at the cathode, i.e., $x / d = 0$⁠. Potential and electric field are self-consistently obtained by solving the Poisson equation using a tridiagonal matrix solver.

The grid size $Δ x$ is chosen by taking the minimum of two condition: $Δ x = min ( 0.1 λ D 0 , d / 50 )$⁠, where $λ D 0$ is the Debye length calculated based on the initial plasma condition. We confirmed that all simulation results satisfy $Δ x / λ D s < 1$⁠, where λ_Ds is the Debye length calculated based on the steady-state plasma properties. The simulation time step $Δ t$ is determined as $Δ t = min ( Δ t CFL , Δ t MCC )$⁠, where $Δ t CFL = Δ x / v e , max$ (here $v e , max$ is the maximum electron velocity inside the domain) is the time step based on the Courant–Friedrichs–Lewy (CFL) condition and $Δ t MCC$ is the time step based on the collision frequency. The time step restriction for particle collisions, i.e., $Δ t MCC$⁠, is necessary to avoid having a particle undergo more than one collisional event over a few time steps. In the present simulations, $ν e n , max Δ t MCC = 0.02$ is employed, where $ν e n , max$ is the maximum electron-neutral collision frequency.

Because of the small current and large voltage in the Townsend discharge, while the current is large and the voltage is small in the glow discharge, a ballast resistor is needed to obtain the Townsend-to-glow transition. In this section, the current–voltage characteristics is obtained from increasing and decreasing the applied voltage, i.e., forward and backward voltage sweep, using a ballast resistor. In these results, the field emission is not taken into account. Thus, the results are applicable for low pressure (p < 1 atm) with a longer anode–cathode gap.

First, the breakdown voltage is obtained when the initial plasma seed density, $n p 0$⁠, is small. It is confirmed that when the initial plasma density seed is small enough, the self-sustaining voltage remains the same (see Sec. IV). The voltage sweep starts from the breakdown voltage (⁠ $V app ≈ V d ≈ V b$⁠), as I_d is negligibly small in the breakdown regime. Then, the applied voltage is slightly increased (e.g., +1 V) to perform the forward voltage sweep.

Figure 1 shows current–voltage characteristics obtained for the forward and backward voltage sweep obtained from the PIC-MCC model. It is to be noted that the reduced current density is given by $j ∗ = j / p 2$⁠, where j is the current density (cf. scaling law^10,38,50). From Ref. 50, the scaling-law is established from the nondimensionalization of Poisson equation considering the drift velocity. The validity of scaling law is confirmed using the PIC-MCC model when $p = 0.01$ and $0.1 atm$ are used. In this work, we focus on pd = 0.24 Pa · m, which corresponds to the left branch as shown later.

Another observation from Fig. 1 is that there exist multiple voltages for the same current density, e.g., during the forward and backward sweep. From a mathematical point of view, this could be considered bifurcation because of the existence of multiple steady-state solutions. Bifurcation has been used to study self-organization of cathode spots in glow and arc discharges.^3,18–20 Figure 1 shows that multiple solutions exist because of the different paths that the plasma takes to reach a steady-state solution. For instance, changing the circuit resistance may change the steady-state current and voltage in the glow regime via the forward sweep from the Townsend regime. During the backward sweep, the existence of space charge modifies the plasma characteristics, thus leading to the steady-state plasma profile being different from the forward sweep. How such circuit effects contribute to self-organization in 2D and 3D configurations is reserved for future investigations.

The plasma properties at the breakdown condition, i.e., V_b = 450 V, are obtained from the PIC-MCC simulation model, as shown later, and are used as the initial condition for a simulation, where V_app is set to $V b + 1$ V $= 451$ V. During the transient, when I_d is still small, V_d can be larger momentarily than V_b, after which the plasma density and thus the discharge current rapidly grow. Then, the voltage drop between the resistor increases and in turn V_d decreases. This is consistent with the Townsend-to-glow transition,^14,24 as I_d is larger but V_d is smaller in the glow discharge regime than in the Townsend regime. In Fig. 1, the red line corresponds to a transition from Townsend discharge, i.e., $V app ≈ V d ≈ V b ≈ 450$ V, $j d ∗ ≈ 17 μ A / ( Pa · m ) 2$⁠, to glow-like discharge, i.e., V_app = 451 V, V_d = 332 V, and $j d ∗ ≈ 2.3 × 10 3 μ A / ( Pa · m ) 2$⁠. In addition, the reason why this transition exhibits a negative differential resistance is because with a small change of V_app, $d V d / d I d = − R$ based on Eq. (2) at steady state. Thus, the slope is dependent on the circuit resistance.

In an actual discharge, it is well known that the current density remains constant in the glow discharge regime. The current increases while the voltage remains the same because of the increase in the plasma column area.^1,50,52–55 If the plasma area reaches the electrode area, then the glow discharge transitions to abnormal glow.⁵⁶ In a 1D simulation, such plasma column expansion cannot be modeled, and hence, constant current density cannot be obtained in a glow-like discharge. In fact, Ref. 24 shows that a glow discharge with constant voltage cannot be obtained in an experiment with a pin-to-plate geometry. Under such a case, the plasma moves into a regime that is similar to an abnormal glow, where both current and voltage increase.

Once the steady-state condition is obtained in the glow-like regime using the forward sweep, the backward sweep is performed, i.e., the applied voltage is decreased. The most important observation from Fig. 1 is that the plasma can be sustained at $V app < V b$ during the backward sweep, which clearly shows a hysteresis similar to the experimental observations.^14,24

Figure 1 shows that decreasing V_app from V_app = 451 V results in decrease in both the current density and discharge voltage, which is indeed a characteristics of the abnormal glow.^14,50 This trend continues down to $V app ≈ 393$ V. The discharge voltage is minimum at V_app = 393 V: V_d = 300 V and $j d ∗ ≈ 1.8 × 10 3 μ A / ( Pa · m ) 2$⁠. However, further decreasing V_app (i.e., V_app < 393 V) results in a decrease in I_d but an increase in V_d. The plasma resistance can be written as $R p = V d / I d$⁠. Alternatively, from the plasma perspective, the electron resistivity is given by $η = ( m e ν m ) / ( e 2 n e )$⁠, where m_e is the electron mass, ν_m is the electron momentum-transfer collision frequency, e is the elementary charge, and n_e is the electron density. The change in the slope during the backward sweep, i.e., 385 V $≤ V app ≤ 393$ V and 393 V $≤ V app ≤ 451$ V, indicates that the plasma resistance R_p changes. The discussion about the plasma resistance will be presented in Sec. V A 3.

If further decreasing V_app below 385 V, the plasma cannot be sustained. Note that at V_app = 385 V, V_d = 355 V and $j d ∗ = 612 μ A / ( Pa · m ) 2$⁠. As can be seen, V_app approaches V_d. With a small amount of change in V_app, I_d keeps decreasing to a point where $V d = V app$ and I_d = 0. An important observation here is that V_app = 385 V is already well below the breakdown condition, i.e., V_b = 450 V. Hence, while this is not surprising, the results show the current–voltage characteristics will not return to the breakdown voltage after the backward voltage sweep. This observation is also consistent with the experiments in Ref. 24, in which the breakdown voltage (forward sweep) and the sustaining voltage (backward sweep) were found to be different.

In this section, instead of varying the applied voltage and modeling the system with a resistance to determine the anode–cathode voltage for self-sustaining discharge, we vary the initial plasma density to obtain the self-sustaining voltage conditions for different initial plasma densities. This approach is motivated by the investigation of a Paschen curve to study the breakdown conditions when the seed initial density is low. In Secs. IV and V, we study the self-sustaining voltage without the field emission, corresponding to the large gap case in Ref. 24. In Sec. VI, we investigate the effects of field emission on the hysteresis.

The self-sustaining voltage is obtained from the PIC-MCC model in the presence of initial plasma density, $n p 0$⁠. Here, the resistor and applied voltage are not considered. A self-sustaining condition refers to the condition when the plasma is at steady state, i.e., plasma is not growing nor decaying.^26,27 Thus, mathematically, the self-sustaining condition is achieved when $∂ ⟨ n i ⟩ / ∂ t = ∂ ⟨ n e ⟩ / ∂ t ≈ 0$⁠, where $⟨ n i ⟩ and ⟨ n e ⟩$ are the spatially averaged ion and electron densities, respectively. When $n p 0$ is small, i.e., $n p 0 ≈ 0$⁠, the voltage at which the discharge is self-sustaining corresponds to the breakdown voltage, denoted as V_b. When $n p 0$ is increased, i.e., $n p 0 ≠ 0$⁠, we repeat the process of obtaining the self-sustaining voltage, denoted as V_s. Here, V_s can be considered the minimum sustaining condition because the plasma decays when V_d < V_s. It might be possible to achieve a steady-state discharge at the same current level with V_d > V_s and depending on the resistor, as discussed in Sec. III.

The self-sustaining voltage, V_d = V_s, is obtained for different initial plasma densities, $n p 0 = 10 18 , 10 19 , 10 20 , 4 × 10 20$⁠, and $10 21 m − 3$⁠. These plasma densities correspond to an ionization fraction of $n p 0 / n g = 10 − 7 – 10 − 4$⁠. Quasineutral and uniform plasma profile is assumed as the initial condition. The initial electron and ion VDFs are assumed to be Maxwellian with electron temperature of T_e = 3 eV and ion temperature of 300 K, respectively. Additionally, the number of initial electron and ion macroparticles is set to 200 000 each, based on a particle convergence study.²⁶ The particle weight of ions and electrons is adjusted based on the initial plasma density. The plasma density in the experiment of Ref. 24 can be estimated from $I d ≈ e n p U e S t$⁠, where $I d ≈ 10 − 100 μ A$⁠, S_t is the surface area of the cathode pin, and U_e is the electron bulk velocity at the cathode. Here, we approximated $U e ≈ k B T e / 2 π m e$ with T_e = 3 eV and $S t ≈ π r t 2$ where k_B is the Boltzmann constant and r_t is the radius of pin (⁠ $r t = 1 μ m$⁠). Using these values, the plasma density is expected to be $10 20 − 10 21 m − 3$⁠, which is consistent with the range of initial plasma density considered for the PIC-MCC simulations in this section.

The self-sustaining voltage V_s is obtained for a given pd using an automatic determination method, which is similar to the Newton's root-finding method.³¹ First, the simulation starts with an arbitrary, small applied voltage $V d k = 0$⁠, where k is the number of iterations, resulting in a plasma profile, $n i k = 0$⁠. If $∂ ⟨ n i k ⟩ / ∂ t < 0$⁠, the applied voltage for the next simulation is set to be $V d k + 1 = 2 V d k$⁠. If $∂ ⟨ n i k + 1 ⟩ / ∂ t > 0 , V d k + 2$ is chosen to be the average of $V min$ and $V max$⁠, where $V min = V k$ and $V max = V k + 1$⁠. If $∂ ⟨ n i k + 2 ⟩ / ∂ t > 0$⁠, then $V max = V d k + 2$⁠. Otherwise, $V min = V d k + 2$⁠. Then, the next iteration value is chosen as $V d k + 3 = ( V max + V min ) / 2$⁠. The procedure is continued until the residual reaches 1% unless stated otherwise to obtain $∂ ⟨ n i ⟩ / ∂ t ≈ 0$⁠. $V d ≈ V b$ can be achieved after sufficient number of iterations (e.g., 10–15). The breakdown voltage obtained from the automatic method is verified with the analytical solution of DC Paschen curve of microgaps where the field emission is important.^57,58

Figure 2 shows the self-sustaining voltage curves obtained from the PIC-MCC model for different initial plasma densities, $n p 0$⁠, in the absence of field emission. The breakdown voltage V_d = V_b is obtained in the case of $n p 0 ≤ 10 18 m − 3$⁠, which corresponds to $n p 0 / n g ≤ 10 − 7$⁠. The minimum pd is observed at $p d ≈ 0.24 Pa · m$⁠, at which the breakdown voltage is V_b = 450 V. At $p d < 0.24 Pa · m$⁠, breakdown does not occur. To further complete the Paschen curve, the solution at higher voltage is obtained by searching for the self-sustaining condition of pd for given constant applied voltages. This multi-valued Paschen curve in DC has been observed in a wide range of experimental observations⁵⁹ and is observed as a consequence of neglecting the fast neutrals, radiative state, and metastable state.³⁷ It is also observed for the $d = 1 cm$ case including PIC-MCC, full fluid moment (FFM), and drift-diffusion (DD) models.³¹ The ionization cross section, σ, generally has the maximum value at a certain electron energy (e.g., ε = 92 eV for argon) and then decreases as the electron energy increases.³⁶ At higher voltage, due to the decrease in $σ ( ϵ )$⁠, where ε is the electron energy, higher pd is required to obtain sufficient ionization, resulting in a multi-valued Paschen curve at the left branch as shown in Fig. 2(b).

Most importantly, Fig. 2 shows that the self-sustaining voltage becomes smaller than the breakdown voltage for finite initial plasma densities, i.e., $n p 0 ≥ 10 19 m − 3$⁠. Four observations can be made. First, in the left branch (low pd), V_s deviates from V_b for $n p 0 > 10 19 m − 3$⁠, while in the right branch (high pd), V_s deviates from V_b for $n p 0 > 10 18 m − 3$⁠. Second, the difference between gas breakdown voltage and the self-sustaining voltage is less sensitive to the initial plasma density for the pd value close to the minimum voltage condition, i.e., $p d ≈ 1 Pa · m$⁠. Third, as shown in Fig. 2(b), the presence of initial plasma density shifts the minimum pd to a smaller value, which means that it is easier to self-sustain the discharge due to the presence of initial charge density. Finally, the decrease in the self-sustaining voltage saturates at $n p 0 ≥ 10 21 m − 3$ in the left branch and at $n p 0 ≥ 4 × 10 20 m − 3$ in the right branch.

Figure 3 shows the time history of space-averaged ion density for (a) $p d = 0.24 Pa · m$ and (b) $p d = 10 Pa · m$⁠, respectively. For the left branch, in Fig. 3(a), the simulation typically reaches steady state around 5 ns. To obtain converged results, the initial number of macroparticles is increased by ten times for $n p 0 = 10 22 m − 3$⁠. For $n p 0 = 10 22 m − 3$⁠, a significant portion of the initial plasma seed is lost to the electrodes during the initial transient. However, it can be seen that the steady-state results for $n p 0 = 10 21$ and 10²² m^–3 are indeed identical, as shown in Fig. 2.

For the right branch (⁠ $p d = 10 Pa · m$⁠), in Fig. 3(b), the simulation time is approximately 10 times larger because of the smaller Knudsen number, i.e., the flow is more collisional. Similar to the simulation results of the left branch, for large $n p 0$⁠, e.g., $n p 0 ≥ 4 × 10 21 m − 3$⁠, the initial number of macroparticles is increased by four times for $n p 0 = 4 × 10 21 m − 3$ and ten times for $n p 0 = 10 22 m − 3$ in order to obtain converged results. Also, it is observed that the growth or decay of electron density is very sensitive to the anode–cathode voltage for $n p 0 ≥ 10 19 m − 3$ at the right branch (⁠ $p d = 10 Pa · m$⁠). Here, V_d is obtained using the automatic method until the difference in V_d reaches 0.05 V, which corresponds to approximately 0.03% residual. In fact, as shown in Fig. 3(b), there is some evidence of temporal oscillations for the space-averaged ion density, which could be correlated with the spatial fluctuation of plasma profile that is shown in Sec. V B. This could be a feature of striations observed in the plasma column.

As shown in Fig. 3, the self-sustaining condition is equivalent to the steady-state condition. Additionally, as shown in Fig. 2, the self-sustaining conditions for different $n p 0$ result in different steady-state solutions, yielding different discharge currents.

Figure 4 shows the current–voltage characteristics for the self-sustaining conditions with assumed initial plasma densities together with the results of Fig. 1. The steady-state results during the backward voltage sweep coincide with the steady-state results of self-sustaining conditions for different $n p 0$⁠. For instance, V_app = 393 and 385 V cases during backward voltage sweep are similar to the self-sustaining results for $n p 0 = 10 21 m − 3$ and $4 × 10 20 m − 3$⁠, respectively. This is an indication that the self-sustaining conditions shown in Fig. 2 are “unique” solutions, despite the use of a uniform initial plasma density. It can therefore be concluded that $V d ≈ V s$ condition for $n p 0 = 10 20 m − 3$ cannot be obtained during the backward sweep for this selected resistance.

In fact, as shown in Figs. 1 and 4, the slope of the current–voltage curve during the forward sweep is based on the circuit resistance. In other words, if a larger circuit resistance is selected, the forward sweep will result in a larger slope in the current–voltage plot. In such a case, it might be possible for the forward sweep to intersect the minimum sustaining condition (i.e., the black line in Fig. 4). Hence, with a larger circuit resistance, the steady-state voltage after the forward sweep may result in a voltage closer to V_b. This is consistent with some experimental observations, where the circuit resistance provides different current–voltage responses.⁵¹ The simulation results suggest that the combination of R and current–voltage characteristics determines whether $V d ≈ V s$ is obtained or not during the backward sweep.

In Sec. IV C, it is observed that the minimum sustaining conditions, V_d = V_s, are unique conditions to the plasma discharge. To understand the change in the plasma resistance, the plasma profiles for the minimum sustaining conditions at the left branch (low pd; pd = 0.24 Pa · m) and at the right branch (high pd; pd = 10 Pa · m) are discussed in this section.

Figure 5 shows the steady-state (self-sustaining) profiles of (a) potential $ϕ$⁠, (b) electric field E, (c) electron axial velocity u_e, and (d) electron temperature T_e for $n p 0 = 10 19 , 10 20 , 4 × 10 20 , 10 21$⁠, and $10 22 m − 3$ at pd = 0.24 Pa · m. The self-sustaining voltages are $V d = 450 , 412 , 355 , 299 ,$ and 299 V, respectively.

In the case of $n p 0 ≤ 10 19 m − 3$⁠, the potential and electric profiles are close to the Laplace solution (⁠ $E 0 = − V d / d$⁠), which is a feature of the Townsend discharge, as shown in Figs. 5(a) and 5(b). Figure 5(c) shows that the electron bulk velocity monotonically increases toward the anode, and Fig. 5(d) shows that the electron temperature T_e reaches approximately 50 eV near the anode, as the electron mean free path is on the order of the anode–cathode gap.

Figure 6 shows the steady-state profiles of ion density n_i and electron density n_e for $n p 0 = 10 19 , 10 20 , 4 × 10 20 , 10 21$⁠, and $10 22 m − 3$⁠. The ion density is approximately uniform except for near the anode and increases as $n p 0$ increases. For Figs. 6(a)–6(c), the ion density is at least an order of magnitude larger than the electron density. Assuming constant u_e, the steady-state continuity equation provides $n e ( x ) = n e , c exp ( α x )$⁠, where $n e , c$ is the electron density at the cathode, $α = ν ion / u e$ is the first Townsend coefficient, and $ν ion$ is the ionization frequency.

As $n i − n e ≈ n i$ for these cases due to non-neutrality, it can be seen from the Gauss's law that $ϵ 0 d E / d x ≈ e n i$⁠, where ϵ₀ is the vacuum permittivity and e is the elementary charge. When n_i is uniform, i.e., $n i ≈ n 0 , E ( x ) = E c + ∫ 0 x e n i / ϵ 0 d x ≈ E c + e n 0 x / ϵ 0$⁠. This linear dependence of electric field can be seen in Fig. 5(b). Indeed, as shown in Figs. 5(a) and 5(b), the space charge effects become important, affecting the self-sustaining condition. It can be seen that the electric field becomes positive near the anode (i.e., reversal of the electric field) for $n p 0 = 10 21 m − 3$⁠, which is associated with the formation of a quasineutral region at $x / d ∈ [ 0.8 , 1.0 ]$⁠. In particular, Figs. 5(c) and 5(d) show a sharp decrease in u_e and T_e for $n p 0 = 10 21 m − 3$ near the anode, i.e., at $x / d ∈ [ 0.9 , 1.0 ]$⁠. Here, since higher electron density is observed as shown in Fig. 6(d), the ionization rate increases rapidly toward the anode. The electrons generated by an ionization event in this region will not have sufficient acceleration as they are already close to the anode, resulting in the decrease in u_e. In particular, in the quasineutral region (⁠ $x / d ∈ [ 0.9 , 1.0 ]$⁠), as the electric field is close to 0, resulting in the rapid decrease in u_e.

As $n p 0$ increases, the location of the maximum electron temperature is shifted toward the cathode and the electron temperature decreases near the anode due to Joule heating, i.e., $| u e E |$⁠. In particular, at $x / d ∈ [ 0.9 , 1.0 ]$⁠, the smaller Joule heating and the larger energy losses due to ionization can contribute to the sharp decrease in electron temperature. It can also be seen that the results of $n p 0 = 10 22 m − 3$ are nearly identical to the $n p 0 = 10 21 m − 3$ results, indicating the saturation of the self-sustaining condition.

As shown in Fig. 6(d), for $n p 0 = 10 21 m − 3$⁠, the electron density increases more rapidly than the exponential increase near the anode (⁠ $x / d ∈ [ 0.8 , 1.0 ]$⁠). This results in the plasma to approach quasineutrality, significantly affecting the potential profile near the anode, as shown in Fig. 5(a). The ion density is still 1.5–2 times larger than the electron density at $x / d = 0.9 – 1.0$ but is much closer to a quasineutral plasma compared to the other lower $n p 0$ cases. For $n p 0 = 10 22 m − 3$⁠, the minimum Debye length at steady state can be estimated as $λ D , min ≈ 0.2$ μm by considering $n e ≈ 5 × 10 21 m − 3$ and $T e ≈ 1 eV$ at $x / d = 0.9$⁠. Here, d = 1 μm. In addition, at $x / d ∈ [ 0 , 0.8 ]$⁠, the Debye length is larger than the gap size. Therefore, most of the region in the anode–cathode gap is non-neutral (i.e., $n i ≫ n e$⁠), which is a feature of the left branch of the Paschen curve, i.e., rarefied regime.

Table I shows the energy balance at the left branch (rarefied regime) for each initial plasma density. Here, we show the volume-averaged collisional loss due to the electron-neutral collisions, $ε c$⁠, the averaged ion energy loss to the cathode, $ε i , c$⁠, and the averaged electron energy loss to the anode, $ε e , a$⁠. It can be seen from Table I that the total energy is conserved, i.e., $V d ≈ ε c + ε i , c + ε e , a$⁠. It is to be noted that the ion energy loss to the anode and the electron energy loss to the cathode are negligibly small.

Several observations can be made from Table I. First, as the Knudsen number is large in the left branch (i.e., low pd regime), electrons need to be sufficiently energized to contribute to generating an electron–ion pair. Thus, $ε c$ stays more or less constant. Second, the electron energy loss to the anode $ε e , a$ decreases as $n p 0$ increases because of the formation of a quasineutral plasma near the anode. The electron velocity and temperature at the anode decrease as $n p 0$ increases, as shown in Figs. 5(c) and 5(d). Third, the ion energy loss to the cathode $ε i , c$ decreases slightly, as stationary ions are already present in the domain as an initial condition. Therefore, it can be concluded that the decrease in self-sustaining voltage for a larger $n p 0$ is associated with the decrease in $ε i , c$ to some extent but, more importantly, the decrease in $ε e , a$⁠. The self-sustaining voltage saturates due to the formation of a quasineutral plasma near the anode minimizing the electron energy loss to the anode.

From the plasma profiles as shown in Figs. 5 and 6, the plasma resistance can be discussed to understand the current–voltage characteristics as shown in Fig. 4. From drift-diffusion formulation, the plasma resistivity¹¹ can be approximated as $η ≈ m e ν e n / ( e 2 n e )$⁠. Using the plasma profiles, let us discuss the plasma resistance in 385 V $≤ V app ≤ 393$ V and 393 V $≤ V app ≤ 451$ V of the backward sweep in Figs. 1 and 4. First, to discuss the transition in 385 ∼ V $≤ V app ≤ 393$ ∼ V, let us consider the increasing $n p 0$ cases of the self-sustaining discharge. For argon, the electron momentum-transfer collision frequency slightly decreases for electron temperatures above 10 eV but remains overall constant, assuming a Maxwellian electron energy distribution function. Thus, starting from the breakdown condition, V_b = 450 V, where n_e is small at $n p 0 = 10 18$ m^–3, the plasma resistivity is large. As $n p 0$ increases the electron density increases, resulting in a smaller plasma resistivity, i.e., $η ∝ n e − 1$⁠. This can explain the change in the slope of the current–voltage characteristics of the self-sustaining conditions. Second, in the abnormal regime (393 V $< V app < 451$ V), the overall structure of the electron temperature and electron density is similar to the result of V_d = V_s at $n p 0 ≥ 10 21 m − 3$ [as shown in Fig. 6(d)], but the near-anode region becomes more quasineutral and spatially larger. At the anode, the electron temperature is almost unchanged. The peak electron density at the anode increases by approximately 20% from V_app = 393 to 451 V. Consequently, the plasma resistance indeed decreases.

However, an important difference between the two regions is the differential resistance, i.e., how much the voltage changes over the change in applied current. The key difference between the two regions that can be clearly observed from the PIC-MCC simulations is the change in the plasma characteristics. The plasma is non-neutral, similar to Townsend discharge, in the 385 V $≤ V app ≤ 393$ V range, as can be seen from Fig. 6. It can be considered that increasing the electron density does not require more voltage until it reaches a quasineutral plasma. On the other hand, the plasma becomes more quasineutral in the range of 393 V $≤ V app ≤ 451$ V. Once a quasineutral plasma is established, the voltage needs to be increased in order to get higher plasma density.

Figure 7 shows the profiles of (a) potential and (b) electron temperature for $n p 0 = 10 18 , 10 19 , 10 20 , 4 × 10 20 , 10 21 , 4 × 10 21$⁠, and $10 22 m − 3$ in the right branch (⁠ $p d = 10 Pa · m$⁠). The self-sustaining voltage is $V d = 250 , 217.25 , 193.15 , 188.95 , 188.85 , 188.80$⁠, and 188.80 V, respectively. Additionally, $Δ x = d / 200$ is used to obtain better resolved plasma profiles for $n p 0 = 10 18 m − 3$ for visualization purpose, i.e., it is confirmed that similar results can be obtained with a coarser grid size, e.g., $Δ x = d / 50$⁠.

The potential profile for $n p 0 ≤ 10 18 m − 3$ is similar to the Laplace solution as the space charge effect is negligible, as shown in Fig. 7(a). The ion density is not equal to the electron density in the domain, i.e., the plasma is non-neutral and the electron density exponentially increases from the cathode to the anode for the breakdown condition, as shown in Fig. 8(a). These two features are characteristics of the Townsend discharge regime and are consistent with the smallest $n p 0$ cases in the left branch as shown in Figs. 5 and 6. Additionally, it was observed in Ref. 31 that the electron velocity distribution function near the cathode becomes a combination of beam electrons accelerated from the cathode and thermalized electrons due to electron-neutral collisions. The oscillating T_e near the cathode can be attributed to the relaxation between such two electron populations.^60,61

Figures 7 and 8 show a few interesting features for the right branch. First, while the self-sustaining voltage decreases between $10 18 m − 3 ≤ n p 0 < 10 20 m − 3$⁠, the self-sustaining voltage saturates at $n p 0 ≥ 10 20 m − 3$⁠. Second, in spite of the saturation of self-sustaining voltage, V_d, at higher $n p 0$⁠, the plasma structure is affected. There are a few distinct features that can be seen for the right branch, particularly for the cases where $n p 0 ≥ 10 20 m − 3$⁠. As $n p 0$ increases, (1) the cathode fall, i.e., the potential drop near the cathode, increases, (2) a plateau (E = 0) region forms and expands, and (3) the anode sheath is consistently electron attracting. Figures 7(b) and 7(c) show that T_e increases near the cathode (⁠ $x / d ∈ [ 0 , 0.05 ]$⁠) proportional to the potential drop in the cathode sheath, the electron temperature becomes small in the plateau region (about 3 eV), but increases near the anode region up to 5 eV. Finally, it can be seen from Figs. 8(b)–8(e) that the plasma becomes quasineutral for $n p 0 ≥ 10 19 m − 3$ in the right branch. The electron temperature oscillates by $1 − 2$ eV in the plasma column region (⁠ $E ≠ 0$⁠) along with an oscillating electric field.⁶² This is likely related to the fact that the spatially averaged plasma density fluctuates in time, for the high pd results, as shown in Fig. 3. The plasma density is maximum near the cathode sheath edge and the diffusion flux (⁠ $d n e / d x < 0$⁠) carries the electron current toward the anode where the drift flux is small (E = 0). For $n p 0 ≥ 4 × 10 21 m − 3$⁠, the quasineutral region is expanded to the anode, and then, the steady-state plasma profiles are saturated. The formation of the quasineutral plasma is indicative of a transition from Townsend discharge to a normal glow discharge.

Table II shows the energy losses in the right branch (collisional regime). In the collisional regime, the electron energy is mainly used for the inelastic collision (excitation, ionization) and the electron energy loss to the anode is small. The energy loss mechanisms in the right branch (collisional regime) are different from the left branch (rarefied regime), which is shown in Table I.

Two interesting observations can be made from Table II. First, the initial decrease in the self-sustaining voltage (up to $n p 0 = 10 20$ m^–3) is primarily driven by the decrease in the ion energy loss toward the cathode, $ε i , c$⁠, while a gradual decrease in collisional loss, $ε c$⁠, can also be observed. It is likely that the existing seed plasma excites sufficient electron emission from the cathode to self-sustain a plasma, lowering the requirement for the voltage. The other notable observation is that the self-sustaining voltage saturates and remains nearly constant for $n p 0 ≥ 4 × 10 20$ m^–3. It can be seen from Table II that the decrease in $ε c$ cancels out with the increase in $ε i , c$⁠. This is consistent with the observation in Fig. 7: the cathode fall increases and the plateau (⁠ $E ≈ 0$⁠) region expands when the initial plasma density is increased from $n p 0 = 10 20$ m^–3 to $n p 0 = 4 × 10 21$ m^–3.

Figure 9 shows the profile of collisional energy loss due to the elastic collisions, $S coll , el = n e ν e l Δ ε el$⁠, ionization, $S coll , iz = n e ν i z Δ ε iz$⁠, and excitation, $S coll , ex = n e Σ j ν e x Δ ε ex$⁠. For $n p 0 ≤ 10 18 m − 3$⁠, both $S coll , ex$ and $S coll , iz$ increase toward the anode according to the electron density profile, as shown in Fig. 8(a). In the case of $n p 0 ≥ 10 20 m − 3$⁠, both ionization and excitation processes are concentrated around the cathode (⁠ $x / d ∈ [ 0 , 0.05 ]$⁠) due to higher electron temperature. These high-energy electrons can ionize gas with less electronic excitation compared to the ionization (i.e., small ionization cost⁶³) from the rate coefficient point of view.³¹ In the quasineutral region, because $T e ≈ 3$ eV, the energy loss due to the elastic collision is comparable to both the excitation and ionization energy losses. This is particularly interesting because this indicates that the minimum sustaining condition is obtained as a condition for minimum energy required to sustain a discharge. In the non-zero electric field region, as $T e ≈ 5$ eV as shown in Fig. 7(b), the electronic excitation increases significantly, while the ionization is still small. Note that the excitation profile obtained from the PIC-MCC model is consistent with the luminescence pattern obtained from the experiment.^1,2 For $n p 0 ≥ 4 × 10 21 m − 3$⁠, the increase in the excitation energy loss is only observed inside the anode sheath.

The expansion of the quasineutral region results in a decrease in the total energy loss due to collisions (ionization and excitation), $ε c$⁠. Simultaneously, ions generated by the ionization travel to the cathode and get accelerated by the cathode sheath potential. As shown in Table II, the ion energy loss, $ε i , c$⁠, increases due to the increase in the cathode sheath potential, as shown in Fig. 7(a) for larger $n p 0$⁠. As a result, the decrease in V_d = V_s can be attributed mainly to the larger decrease in $ε c$ for $10 18 m − 3 < n p 0 < 10 20 m − 3$⁠. Additionally, V_d = V_s saturates, since the increase in the ion energy loss to the cathode cancels out the decrease in the total collisional energy loss for $n p 0 ≥ 10 20 m − 3$⁠. For $n p 0 = 4 × 10 21 m − 3$⁠, as the quasineutral region reaches near the anode, the collisional loss due to excitation is minimized, while the cathode sheath potential $ϕ s$ is close to V_d (i.e., $ϕ s ≈ V d$⁠). This results in the saturation of steady-state plasma profile at $n p 0 ≥ 4 × 10 21 m − 3$⁠. From the perspective of the global energy balance,¹¹ the maximum plasma density near the anode sheath is also saturated at $n p 0 ≥ 10 20 m − 3$ shown in Figs. 8(c)–8(e) when the total energy loss in the system $ε loss ≈ ε c + ϕ s$ is approximately constant.

For smaller gaps (e.g., $d ≈ 1 μ m$⁠), if the discharge voltage is on the order of 100 V, the electric field is on the other of 10⁸ V/m, in which field emission becomes important. Hence, field emission is introduced according to the Fowler–Nordheim equation, shown in Eq. (1), assuming β = 20. From the current–voltage characteristics, as shown in Fig. 4, the self-sustaining voltage can be smaller than the breakdown voltage during the backward voltage sweep. It is further shown that the self-sustaining condition assuming an initial plasma density is consistent with the steady-state solution during the backward voltage sweep. Hence, instead of conducting a voltage sweep to study the hysteresis, we investigate the self-sustaining voltages for different initial plasma densities in the presence of field emission. If the breakdown voltage and the self-sustaining voltage are the same, this indicates that the voltage remains constant, while the current increases. Thus, it can be expected that the self-sustaining voltage also remains constant during the backward sweep, resulting in mitigation of the hysteresis.

Figure 10 shows the effects of field emission on self-sustaining voltages. In particular, the breakdown conditions without initial space charge effects can be seen for the results when $n p 0 ≤ 10 19 m − 3$⁠, which approximately corresponds to $n p 0 / n g ≤ 10 − 6$⁠. In the presence of field emission effects, the required voltage for a self-sustaining condition decreases monotonically as the gap size, d, decreases.^30,57,58,64 In addition, microgap breakdown (with field emission) allows for the self-sustaining plasma at small pd values, below the minimum pd values for cases in the absence of field emission as shown in Fig. 2 and even to pd = 0.

As shown in Fig. 10, the self-sustaining voltage curves can be divided into three regimes (e.g., for the $n p 0 = 10 21 m − 3$ case): pd > 0.28 Pa · m, 0.12 Pa · m $< p d ≤ 0.28$ Pa · m, and $p d ≤ 0.12$ Pa · m. (I) At $p d ≥ 0.28$ Pa · m, the self-sustaining voltage is identical to the cases without field emission for all cases. As discussed in Sec. IV B, one can expect a hysteresis due to the initial plasma density effects when forward and backward sweeping the voltage. (II) Field emission is dominant at $p d ≤ 0.12$ Pa · m, as the self-sustaining voltage becomes identical over all initial plasma densities. In this case, it can be considered that hysteresis is mitigated as the breakdown voltage (in the absence of initial plasma density) and the self-sustaining voltage (in the presence of initial plasma density) are identical, which is consistent with the experimental observations in Ref. 24.

(III) At the intermediate regime, i.e., 0.12 Pa · m $< p d ≤ 0.28$ Pa · m, the self-sustaining voltages in the presence of field emission (solid lines) deviate from those when the field emission is negligible (dashed lines). In this regime, the self-sustaining voltages still depend on the initial plasma density, $n p 0$⁠. Additionally, the reason why the self-sustaining voltage for $n p 0 > 10 19 m − 3$ is lower than the breakdown voltage is due to the space charge effects. In the presence of an initial plasma, a sheath is formed near the cathodes, which enhances the local electric field and hence the field emission current density at the cathode. A higher initial plasma density provides larger electric field near the cathode, resulting in larger field emission. Consequently, the self-sustaining voltage decreases as the initial plasma density becomes high as shown in Fig. 10. This enhancement is called the ion-enhanced field emission in Refs. 65 and 66.

In this paper, the hysteresis of a DC plasma discharge is investigated and analyzed using a 1D PIC-MCC simulation. First, the forward and backward voltage sweeps are demonstrated accounting for a constant ballast resistor in the PIC-MCC model without field emission. The simulation results show an evidence of hysteresis in the current–voltage characteristics, in which the steady-state voltage during the backward voltage sweep can be smaller than the breakdown voltage. In the forward voltage sweep, with a small increase in the applied voltage, transition from the Townsend regime to the glow regime is observed, resulting in a decreased discharge voltage and an increased current. When the applied voltage is decreased from the steady-state result in the glow regime obtained from the forward voltage sweep, the decrease in the discharge voltage linearly decreases with as the current density decreases until it gets to a minimum voltage. Further decreasing the applied voltage results in an increase in discharge voltage, while the discharge current decreases, indicating a change in the plasma resistance during the backward voltage sweep. The PIC-MCC simulation results indicate that the change in the differential resistance of the plasma is due to the plasma being non-neutral (feature of Townsend discharge) or quasineutral (feature of glow discharge). Finally, when the applied voltage is decreased further, the plasma cannot be sustained and the current goes to zero. Consequently, the applied voltage when the current becomes zero via the backward voltage sweep is smaller than the breakdown voltage, which is consistent with the experimental observation.²⁴ 

Second, to understand the discharge conditions of the backward sweep, the self-sustaining voltages are obtained by considering several initial plasma densities, $n p 0$⁠, without any ballast resistor. The breakdown voltage and the self-sustaining voltage are obtained when the initial space charge density is small and when a finite initial plasma density is considered, receptively. It is observed that $n p 0 / n g ≈ 10 − 6 – 10 − 7$⁠, where n_g is the background gas density, is the threshold between the breakdown and self-sustaining discharge conditions. The 1D PIC-MCC simulations in the absence of field emission (representative of large gap experiments) suggest that the self-sustaining voltage in the presence of initial plasma density is smaller than the breakdown voltage. The saturation of the voltage decrease at larger $n p 0$⁠, together with the observation that the plasma within the anode–cathode gap becoming quasineutral, indicates a transition from a Townsend discharge to a glow discharge. From the current–voltage analysis, it is shown that self-sustaining conditions for different initial plasma densities result in unique steady-state solutions, yielding different discharge current. In addition, the discharge current and voltage obtained from the self-sustaining cases with different $n p 0$ coincide with the steady-state solutions during the backward sweep, indicating that the discharge distinguishes by the backward sweep differently from the generation of glow discharge by the forward sweep.

To investigate the plasma behavior associated with the transition between the breakdown and self-sustaining discharges for the minimum sustaining voltage condition V_d = V_s, steady-state plasma profiles obtained from the PIC-MCC model are analyzed. In both rarefied (left branch; small pd) and collisional (right branch; large pd) regimes, the effects of the space charge and formation of a quasineutral plasma region are key to explain the transition. In the rarefied regime, the self-induced electric field near the anode causes a decrease in the electron velocity, Joule heating, and thus electron temperature. Consequently, as $n p 0$ increases, the electron energy loss to the anode decreases, which leads to a smaller V_d required to sustain the plasma. When the electron energy loss to the anode is minimized, the minimum sustaining voltage and steady-state plasma profiles are saturated at larger $n p 0$⁠. The transition from Townsend to glow discharges is consistent with the change in the plasma resistance during the backward voltage sweep. On the other hand, in the collisional regime, the space charge due to $n p 0$ influences the cathode sheath potential (cf. cathode fall), which increases the electron temperature and the ionization rate near the cathode. As $n p 0$ increases, a quasineutral plasma and the formation of a plateau (E = 0) region are observed. In this region, the Joule heating is small so that the electron temperature decreases. However, the electric field near the anode is still finite, and the anode sheath is electron attracting. Due to the expansion of the quasineutral plasma region, the total energy loss of electrons due to collisions decreases. However, simultaneously, the ion energy loss to the cathode increases due to the increase in the cathode sheath potential. The decreased electron energy loss and the increased ion energy loss result in the saturation of the self-sustaining voltage in the right branch for larger $n p 0$⁠.

Finally, self-sustaining voltage curves are obtained in the presence of field emission, which is observed for microgaps (e.g., on the order of 1 μm). Field emission allows for the formation of a self-sustaining discharge at a small pd value, even down to pd = 0. In the pd regime where field emission is dominant, the breakdown voltage (with forward voltage sweeping) and the self-sustaining voltage (with backward voltage sweeping) converge to the same value for a wide range of initial plasma density cases. This indicates that the hysteresis will not occur in this regime and is consistent with the experimental observation.²⁴ In the absence of the hysteresis, the same current–voltage characteristics can be obtained for the forward and backward sweep. Additionally, the present 1D PIC-MCC simulations show the effects of initial plasma density in an intermediate pd regime, which can contribute to a hysteresis. The space-charged enhanced field emission is similarly observed in experiments.^64,67

This material is based on work supported by the U.S. Department of Energy, Office of Science, Office of Fusion Energy Sciences, under Award No. DE-SC0020623, the Office of Naval Research under Award No. N00014-21-1-2698, and Lam Research Corporation. Y.Y. acknowledges the Japan Society for the Promotion of Science (JSPS) for the postdoctoral fellowship. Some of the computing for this project was performed on the Sherlock cluster. We would like to thank Stanford University and the Stanford Research Computing Center for providing computational resources and support that contributed to these research results.

The authors have no conflicts to disclose.

Yusuke Yamashita: Conceptualization (equal); Formal analysis (lead); Methodology (lead); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). Kentaro Hara: Conceptualization (equal); Formal analysis (supporting); Funding acquisition (lead); Methodology (supporting); Project administration (equal); Software (supporting); Validation (supporting); Visualization (supporting); Writing – original draft (supporting); Writing – review & editing (equal). Saravanapriyan Sriraman: Funding acquisition (supporting); Project administration (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.