In direct-current (DC) discharge, it is well known that hysteresis is observed between the Townsend (gas breakdown) and glow regimes. Forward and backward voltage sweep is performed using a one-dimensional particle-in-cell Monte Carlo collision (PIC-MCC) model considering a ballast resistor. When increasing the applied voltage after reaching the breakdown voltage (*V _{b}*), transition from Townsend to glow discharges is observed. When decreasing the applied voltage from the glow regime, the discharge voltage (

*V*) between the anode–cathode gap can be smaller than the breakdown voltage, resulting in a hysteresis, which is consistent with experimental observations. Next, the PIC-MCC model is used to investigate the self-sustaining voltage (

_{d}*V*) in the presence of finite initial plasma densities between the anode and cathode gap. It is observed that the self-sustaining voltage coincides with the discharge voltage obtained from the backward voltage sweep. In addition, the self-sustaining voltage decreases with increased initial plasma density and saturates above a certain initial plasma density, which indicates a change in plasma resistivity. The decrease in self-sustaining voltage is associated with the electron heat loss at the anode for the low

_{s}*pd*(rarefied) regime. In the high

*pd*(collisional) regime, the ion energy loss toward the cathode due to the cathode fall and the inelastic collision loss of electrons in the bulk discharge balance out. Finally, it is demonstrated that the self-sustaining voltage collapses to a singular value, despite the presence of a initial plasma, for microgaps when field emission is dominant, which is also consistent with experimental observations.

## I. INTRODUCTION

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.^{9} 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.^{11} Recent study also shows the transition between stationary and oscillatory modes in DC glow discharges.^{12} 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

^{13}). 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*is the discharge voltage between the anode and cathode gap and

_{d}*I*is the discharge current. Some studies suggest the existence of multiple solutions assuming steady-state conditions, often called bifurcation, in DC plasma systems.

_{d}^{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 \u2212 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*and thus a decrease in the discharge voltage

_{d}*V*(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

_{d}*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,

^{25}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*, the plasma decays. On the other hand, when the voltage is above the breakdown voltage, i.e.,

_{b}*V*>

_{d}*V*, 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.,

_{b}*V*=

_{d}*V*. A variety of computational models have been developed to predict the breakdown characteristics in both DC and RF discharges.

_{b}^{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.

^{31}

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*is obtained. The simulation results suggest that

_{b}*V*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

_{s}*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.

## II. SIMULATION MODEL

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.

### A. PIC-MCC simulation model

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.^{36} 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 $ \gamma i = 0.1$ is assumed for simplicity, while it has been reported that

*γ*has ion energy dependency.

_{i}^{37}

^{41}

^{,}

*J*is the current density due to field emission,

_{FN}*A*and

*B*are Fowler–Nordheim constants,

^{42}

^{,}

*τ*and

*v*are the modified coefficients, $ \varphi W$ is the work function,

*β*is the field enhancement parameter, and

*E*is the electric field at the cathode, which is self-consistently obtained by solving the Poisson equation. Here,

_{c}*τ*and

*v*can be approximated by $ \tau 2 \u2248 1.1$ and $ v \u2248 0.95 \u2212 y 2$, where

*y*is the non-dimensional parameter defined as $ y = \Delta \varphi W / \varphi W$ and $ \Delta \varphi W = e \beta | E c | / ( 4 \pi \epsilon 0 )$ is Schottky lowering of the work function barrier.

^{43}The field enhancement parameter is determined by the surface roughness of the electrode, which is typically $ \beta = 1 \u2013 50$.

^{24,44,45}In this paper,

*β*= 20 is assumed and $ \varphi W = 4.5$ eV is used for tungsten based on the experimental work.

^{24}

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)^{46} and 1D cross field discharge test case.^{33} Also, the Paschen curve for the gas breakdown condition under $ d = 1 \u2009 \u2009 cm$ is also benchmarked with a full-fluid moment model and drift-diffusion (DD) model.^{31}

### B. Numerical conditions

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 \u2248 2.5 \xd7 10 25 \u2009 m \u2212 3$, where *n _{g}* is the neutral number density. The distance between the anode and cathode

*d*is varied $ d \u2208 [ 1 \u2009 \mu m , 100 \u2009 \mu 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 \u2009 Pa \xb7 m$ (e.g., Meek condition),

^{10,47}while this paper focus on $ p d = 0.1 \u2013 10 \u2009 Pa \xb7 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,

^{48}and the glow-to-arc transition

^{49}is not observed in this work.

Dirichlet boundary conditions are employed for the anode and cathode. $ \varphi = V d$ is fixed at the anode, i.e., $ x / d = 1$, while $ \varphi = 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 $ \Delta x$ is chosen by taking the minimum of two condition: $ \Delta x = min ( 0.1 \lambda D 0 , d / 50 )$, where $ \lambda D 0$ is the Debye length calculated based on the initial plasma condition. We confirmed that all simulation results satisfy $ \Delta x / \lambda D s < 1$, where *λ _{Ds}* is the Debye length calculated based on the steady-state plasma properties. The simulation time step $ \Delta t$ is determined as $ \Delta t = min ( \Delta t CFL , \u2009 \Delta t MCC )$, where $ \Delta t CFL = \Delta 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 $ \Delta t MCC$ is the time step based on the collision frequency. The time step restriction for particle collisions, i.e., $ \Delta t MCC$, is necessary to avoid having a particle undergo more than one collisional event over a few time steps. In the present simulations, $ \nu e n , max \Delta t MCC = 0.02$ is employed, where $ \nu e n , max$ is the maximum electron-neutral collision frequency.

## III. HYSTERESIS BETWEEN TOWNSEND AND GLOW DISCHARGES

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.

### A. Circuit-plasma model

*R*is the resistance and

*C*is the capacitance. The capacitance is used to mitigate the oscillations and achieve a steady-state discharge solution. In the present work,

*I*is obtained by the electron current and ion current at the anode. To reduce the statistical noise,

_{d}*I*is obtained by time-averaging over 10

_{d}^{4}time steps and then,

*V*is updated using Eq. (2) with a time step of $ \Delta t c = 10 4 \Delta t$, where $ \Delta t$ is the plasma time step for the PIC-MCC model. Consider that

_{d}*S*is the surface area of the plasma discharge. In the present 1D work, $ j = I d / S d$ is the current density, and $ R S d = 5 \u2009 \mu \Omega \xb7 m 2$ and $ C / S d = 11 \u2009 \mu F / m 2$ are considered. It is to be noted that if considering $ S d \u223c 10 \u2212 12 \u2009 m \u2212 2$ based on the experimental condition in Ref. 24,

_{d}*R*results in an order of $ M \Omega $, which is comparable to the resistor that is used in the experiment.

### B. Hysteresis in current–voltage characteristics

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 \u2248 V d \u2248 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 \u2217 = 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 \u2009 \u2009 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.

#### 1. Forward sweep

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*is set to $ V b + 1$ V $ = 451$ V. During the transient, when

_{app}*I*is still small,

_{d}*V*can be larger momentarily than

_{d}*V*, after which the plasma density and thus the discharge current rapidly grow. Then, the voltage drop between the resistor increases and in turn

_{b}*V*decreases. This is consistent with the Townsend-to-glow transition,

_{d}^{14,24}as

*I*is larger but

_{d}*V*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 \u2248 V d \u2248 V b \u2248 450$ V, $ j d \u2217 \u2248 17 \mu A / ( Pa \xb7 m ) 2$, to glow-like discharge, i.e.,

_{d}*V*= 451 V,

_{app}*V*= 332 V, and $ j d \u2217 \u2248 2.3 \xd7 10 3 \u2009 \mu A / ( Pa \xb7 m ) 2$. In addition, the reason why this transition exhibits a negative differential resistance is because with a small change of

_{d}*V*, $ d V d / d I d = \u2212 R$ based on Eq. (2) at steady state. Thus, the slope is dependent on the circuit resistance.

_{app}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.^{56} 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.

#### 2. Backward sweep

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*= 451 V results in decrease in both the current density and discharge voltage, which is indeed a characteristics of the abnormal glow.

_{app}^{14,50}This trend continues down to $ V app \u2248 393$ V. The discharge voltage is minimum at

*V*= 393 V:

_{app}*V*= 300 V and $ j d \u2217 \u2248 1.8 \xd7 10 3 \mu A / ( Pa \xb7 m ) 2$. However, further decreasing

_{d}*V*(i.e.,

_{app}*V*< 393 V) results in a decrease in

_{app}*I*but an increase in

_{d}*V*. The plasma resistance can be written as $ R p = V d / I d$. Alternatively, from the plasma perspective, the electron resistivity is given by $ \eta = ( m e \nu m ) / ( e 2 n e )$, where

_{d}*m*is the electron mass,

_{e}*ν*is the electron momentum-transfer collision frequency,

_{m}*e*is the elementary charge, and

*n*is the electron density. The change in the slope during the backward sweep, i.e., 385 V $ \u2264 V app \u2264 393$ V and 393 V $ \u2264 V app \u2264 451$ V, indicates that the plasma resistance

_{e}*R*changes. The discussion about the plasma resistance will be presented in Sec. V A 3.

_{p}If further decreasing *V _{app}* below 385 V, the plasma cannot be sustained. Note that at

*V*= 385 V,

_{app}*V*= 355 V and $ j d \u2217 = 612 \u2009 \mu A / ( Pa \xb7 m ) 2$. As can be seen,

_{d}*V*approaches

_{app}*V*. With a small amount of change in

_{d}*V*,

_{app}*I*keeps decreasing to a point where $ V d = V app$ and

_{d}*I*= 0. An important observation here is that

_{d}*V*= 385 V is already well below the breakdown condition, i.e.,

_{app}*V*= 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.

_{b}## IV. EFFECTS OF INITIAL PLASMA DENSITY ON SELF-SUSTAINING DISCHARGES

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.

### A. Numerical setup

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 $ \u2202 \u27e8 n i \u27e9 / \u2202 t = \u2202 \u27e8 n e \u27e9 / \u2202 t \u2248 0$, where $ \u27e8 n i \u27e9 \u2009 \u2009 and \u2009 \u2009 \u27e8 n e \u27e9$ are the spatially averaged ion and electron densities, respectively. When $ n p 0$ is small, i.e., $ n p 0 \u2248 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 \u2260 0$, we repeat the process of obtaining the self-sustaining voltage, denoted as

*V*. Here,

_{s}*V*can be considered the minimum sustaining condition because the plasma decays when

_{s}*V*<

_{d}*V*. It might be possible to achieve a steady-state discharge at the same current level with

_{s}*V*>

_{d}*V*and depending on the resistor, as discussed in Sec. III.

_{s}The self-sustaining voltage, *V _{d}* =

*V*, is obtained for different initial plasma densities, $ n p 0 = 10 18 , 10 19 , 10 20 , 4 \xd7 10 20$, and $ 10 21 \u2009 m \u2212 3$. These plasma densities correspond to an ionization fraction of $ n p 0 / n g = 10 \u2212 7 \u2013 10 \u2212 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

_{s}*T*= 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.

_{e}^{26}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 \u2248 e n p U e S t$, where $ I d \u2248 10 \u2212 100 \u2009 \mu A$,

*S*is the surface area of the cathode pin, and

_{t}*U*is the electron bulk velocity at the cathode. Here, we approximated $ U e \u2248 k B T e / 2 \pi m e$ with

_{e}*T*= 3 eV and $ S t \u2248 \pi r t 2$ where

_{e}*k*is the Boltzmann constant and

_{B}*r*is the radius of pin ( $ r t = 1 \u2009 \mu m$). Using these values, the plasma density is expected to be $ 10 20 \u2212 10 21 m \u2212 3$, which is consistent with the range of initial plasma density considered for the PIC-MCC simulations in this section.

_{t}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.

^{31}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 $ \u2202 \u27e8 n i k \u27e9 / \u2202 t < 0$, the applied voltage for the next simulation is set to be $ V d k + 1 = 2 V d k$. If $ \u2202 \u27e8 n i k + 1 \u27e9 / \u2202 t > 0 , \u2009 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 $ \u2202 \u27e8 n i k + 2 \u27e9 / \u2202 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 $ \u2202 \u27e8 n i \u27e9 / \u2202 t \u2248 0$. $ V d \u2248 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}

### B. Breakdown voltage *V*_{b} (for $ n p 0 \u2248 0$) and self-sustaining voltage *V*_{s} (for $ n p 0 \u2260 0$) as a function of *pd*

_{b}

_{s}

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*is obtained in the case of $ n p 0 \u2264 10 18 \u2009 m \u2212 3$, which corresponds to $ n p 0 / n g \u2264 10 \u2212 7$. The minimum

_{b}*pd*is observed at $ p d \u2248 0.24 \u2009 Pa \xb7 m$, at which the breakdown voltage is

*V*= 450 V. At $ p d < 0.24 \u2009 \u2009 Pa \xb7 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

_{b}*pd*for given constant applied voltages. This multi-valued Paschen curve in DC has been observed in a wide range of experimental observations

^{59}and is observed as a consequence of neglecting the fast neutrals, radiative state, and metastable state.

^{37}It is also observed for the $ d = 1 \u2009 cm$ case including PIC-MCC, full fluid moment (FFM), and drift-diffusion (DD) models.

^{31}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.

^{36}At higher voltage, due to the decrease in $ \sigma ( \u03f5 )$, 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 \u2265 10 19 \u2009 m \u2212 3$. Four observations can be made. First, in the left branch (low *pd*), *V _{s}* deviates from

*V*for $ n p 0 > 10 19 \u2009 m \u2212 3$, while in the right branch (high

_{b}*pd*),

*V*deviates from

_{s}*V*for $ n p 0 > 10 18 \u2009 m \u2212 3$. Second, the difference between gas breakdown voltage and the self-sustaining voltage is less sensitive to the initial plasma density for the

_{b}*pd*value close to the minimum voltage condition, i.e., $ p d \u2248 1 \u2009 \u2009 Pa \xb7 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 \u2265 10 21 \u2009 m \u2212 3$ in the left branch and at $ n p 0 \u2265 4 \xd7 10 20 \u2009 m \u2212 3$ in the right branch.

Figure 3 shows the time history of space-averaged ion density for (a) $ p d = 0.24 \u2009 \u2009 Pa \xb7 m$ and (b) $ p d = 10 Pa \xb7 \u2009 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 \u2009 m \u2212 3$. For $ n p 0 = 10 22 m \u2212 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^{22} m^{–3} are indeed identical, as shown in Fig. 2.

For the right branch ( $ p d = 10 \u2009 \u2009 Pa \xb7 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 \u2265 4 \xd7 10 21 m \u2212 3$, the initial number of macroparticles is increased by four times for $ n p 0 = 4 \xd7 10 21 \u2009 m \u2212 3$ and ten times for $ n p 0 = 10 22 m \u2212 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 \u2265 10 19 m \u2212 3$ at the right branch ( $ p d = 10 \u2009 \u2009 Pa \xb7 m$). Here, *V _{d}* is obtained using the automatic method until the difference in

*V*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.

_{d}### C. Relation to the hysteresis

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 \u2212 3$ and $ 4 \xd7 10 20 m \u2212 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 \u2248 V s$ condition for $ n p 0 = 10 20 \u2009 m \u2212 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.

^{51}The simulation results suggest that the combination of

*R*and current–voltage characteristics determines whether $ V d \u2248 V s$ is obtained or not during the backward sweep.

## V. PLASMA CHARACTERISTICS OF SELF-SUSTAINING DISCHARGES

In Sec. IV C, it is observed that the minimum sustaining conditions, *V _{d}* =

*V*, 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

_{s}*pd*;

*pd*= 0.24 Pa · m) and at the right branch (high

*pd*;

*pd*= 10 Pa · m) are discussed in this section.

### A. Left branch: Rarefied regime

#### 1. Plasma profiles

Figure 5 shows the steady-state (self-sustaining) profiles of (a) potential $\varphi $, (b) electric field *E*, (c) electron axial velocity *u _{e}*, and (d) electron temperature

*T*for $ n p 0 = 10 19 , 10 20 , 4 \xd7 10 20 , 10 21$, and $ 10 22 \u2009 m \u2212 3$ at

_{e}*pd*= 0.24 Pa · m. The self-sustaining voltages are $ V d = 450 , \u2009 412 , 355 , \u2009 299 ,$ and 299 V, respectively.

In the case of $ n p 0 \u2264 10 19 \u2009 m \u2212 3$, the potential and electric profiles are close to the Laplace solution ( $ E 0 = \u2212 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*for $ n p 0 = 10 19 , 10 20 , 4 \xd7 10 20 , 10 21$, and $ 10 22 \u2009 m \u2212 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

_{e}*u*, the steady-state continuity equation provides $ n e ( x ) = n e , c \u2009 exp ( \alpha x )$, where $ n e , c$ is the electron density at the cathode, $ \alpha = \nu ion / u e$ is the first Townsend coefficient, and $ \nu ion$ is the ionization frequency.

_{e}As $ n i \u2212 n e \u2248 n i$ for these cases due to non-neutrality, it can be seen from the Gauss's law that $ \u03f5 0 d E / d x \u2248 e n i$, where *ϵ*_{0} is the vacuum permittivity and *e* is the elementary charge. When *n _{i}* is uniform, i.e., $ n i \u2248 n 0 , \u2009 E ( x ) = E c + \u222b 0 x e n i / \u03f5 0 d x \u2248 E c + e n 0 x / \u03f5 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 \u2009 m \u2212 3$, which is associated with the formation of a quasineutral region at $ x / d \u2208 [ 0.8 , \u2009 1.0 ]$. In particular, Figs. 5(c) and 5(d) show a sharp decrease in

*u*and

_{e}*T*for $ n p 0 = 10 21 \u2009 m \u2212 3$ near the anode, i.e., at $ x / d \u2208 [ 0.9 , \u2009 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

_{e}*u*. In particular, in the quasineutral region ( $ x / d \u2208 [ 0.9 , \u2009 1.0 ]$), as the electric field is close to 0, resulting in the rapid decrease in

_{e}*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 \u2208 [ 0.9 , \u2009 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 \u2212 3$ are nearly identical to the $ n p 0 = 10 21 m \u2212 3$ results, indicating the saturation of the self-sustaining condition.

As shown in Fig. 6(d), for $ n p 0 = 10 21 \u2009 m \u2212 3$, the electron density increases more rapidly than the exponential increase near the anode ( $ x / d \u2208 [ 0.8 , \u2009 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 \u2013 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 \u2009 m \u2212 3$, the minimum Debye length at steady state can be estimated as $ \lambda D , min \u2248 0.2$ *μ*m by considering $ n e \u2248 5 \xd7 10 21 \u2009 m \u2212 3$ and $ T e \u2248 1 \u2009 eV$ at $ x / d = 0.9$. Here, *d* = 1 *μ*m. In addition, at $ x / d \u2208 [ 0 , \u2009 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 \u226b n e$), which is a feature of the left branch of the Paschen curve, i.e., rarefied regime.

*u*decreases near the anode, yielding a gradient of

_{e}*u*, as can be seen in Fig. 5(c). If $ d u e / d x \u2260 0$, the steady-state continuity equation can be written as

_{e}*u*. If $ \u2202 ( ln \u2009 u e ) / \u2202 x$ is constant, $ \alpha eff$ is constant, which results in $ n e ( x ) = n e , c \u2009 exp ( \alpha eff x )$. In fact, as shown in Fig. 5(c), $ u e > 0$ and $ d u e / d x < 0$ as $ n p 0$ increases. Thus, it can be seen from Eq. (3) that $ \alpha eff > \alpha $, which can reduce the space charge, i.e., $ n i \u2212 n e$, near the anode, promoting a quasineutral plasma near the anode, as can be seen in Fig. 6(d).

_{e}#### 2. Energy balance

*k*th type, $ k = el , \u2009 ex , \u2009 iz$ represents elastic collision, excitation, and ionization, respectively, and $ N \u0307 k = \u222b 0 d n \u0307 k d x$ is the volume-integrated collisional rate. The quantity $ \epsilon c$ indicates how much energy is required to generate one electron-ion pair. As the energy loss due to elastic collisions is small in the regime of interest, $ \epsilon c$ is mainly determined by the excitation and ionization rates, which is a function of the electron mean energy.

^{31}

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, $ \epsilon c$, the averaged ion energy loss to the cathode, $ \epsilon i , c$, and the averaged electron energy loss to the anode, $ \epsilon e , a$. It can be seen from Table I that the total energy is conserved, i.e., $ V d \u2248 \epsilon c + \epsilon i , c + \epsilon 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.

$ n p 0 ( m \u2212 3 )$ . | $ \epsilon c$ (eV) . | $ \epsilon e , a$ (eV) . | $ \epsilon i , c$ (eV) . | Sum (eV) . |
---|---|---|---|---|

$ \u2264 10 19$ | 21 | 127 | 303 | 451 |

10^{20} | 21 | 113 | 279 | 413 |

$ 4 \xd7 10 20$ | 22 | 71 | 263 | 356 |

$ \u2265 10 21$ | 23 | 4 | 271 | 298 |

$ n p 0 ( m \u2212 3 )$ . | $ \epsilon c$ (eV) . | $ \epsilon e , a$ (eV) . | $ \epsilon i , c$ (eV) . | Sum (eV) . |
---|---|---|---|---|

$ \u2264 10 19$ | 21 | 127 | 303 | 451 |

10^{20} | 21 | 113 | 279 | 413 |

$ 4 \xd7 10 20$ | 22 | 71 | 263 | 356 |

$ \u2265 10 21$ | 23 | 4 | 271 | 298 |

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, $ \epsilon c$ stays more or less constant. Second, the electron energy loss to the anode $ \epsilon 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 $ \epsilon 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 $ \epsilon i , c$ to some extent but, more importantly, the decrease in $ \epsilon 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.

#### 3. Plasma resistance

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^{11} can be approximated as $ \eta \u2248 m e \nu e n / ( e 2 n e )$. Using the plasma profiles, let us discuss the plasma resistance in 385 V $ \u2264 V app \u2264 393$ V and 393 V $ \u2264 V app \u2264 451$ V of the backward sweep in Figs. 1 and 4. First, to discuss the transition in 385 ∼ V $ \u2264 V app \u2264 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*is small at $ n p 0 = 10 18$ m

_{e}^{–3}, the plasma resistivity is large. As $ n p 0$ increases the electron density increases, resulting in a smaller plasma resistivity, i.e., $ \eta \u221d n e \u2212 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*at $ n p 0 \u2265 10 21 m \u2212 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

_{s}*V*= 393 to 451 V. Consequently, the plasma resistance indeed decreases.

_{app}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 $ \u2264 V app \u2264 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 $ \u2264 V app \u2264 451$ V. Once a quasineutral plasma is established, the voltage needs to be increased in order to get higher plasma density.

### B. Right branch: Collisional regime

#### 1. Plasma profiles

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

The potential profile for $ n p 0 \u2264 10 18 \u2009 m \u2212 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 \u2009 m \u2212 3 \u2264 n p 0 < 10 20 \u2009 m \u2212 3$, the self-sustaining voltage saturates at $ n p 0 \u2265 10 20 \u2009 m \u2212 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 \u2265 10 20 \u2009 m \u2212 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*increases near the cathode ( $ x / d \u2208 [ 0 , \u2009 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 \u2265 10 19 \u2009 m \u2212 3$ in the right branch. The electron temperature oscillates by $ 1 \u2212 2$ eV in the plasma column region ( $ E \u2260 0$) along with an oscillating electric field.

_{e}^{62}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 \u2265 4 \xd7 10 21 \u2009 m \u2212 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.

#### 2. Energy balance

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.

$ n p 0 \u2009 ( m \u2212 3 )$ . | $ \epsilon c$ (eV) . | $ \epsilon e , a$ (eV) . | $ \epsilon i , c$ (eV) . | Sum (eV) . |
---|---|---|---|---|

$ \u2264 10 18$ | 101 | 10 | 141 | 252 |

10^{19} | 95 | 11 | 112 | 218 |

10^{20} | 85 | 8 | 101 | 194 |

$ 4 \xd7 10 20$ | 70 | 6 | 115 | 191 |

10^{21} | 48 | 7 | 136 | 191 |

$ \u2265 4 \xd7 10 21$ | 33 | 8 | 149 | 190 |

$ n p 0 \u2009 ( m \u2212 3 )$ . | $ \epsilon c$ (eV) . | $ \epsilon e , a$ (eV) . | $ \epsilon i , c$ (eV) . | Sum (eV) . |
---|---|---|---|---|

$ \u2264 10 18$ | 101 | 10 | 141 | 252 |

10^{19} | 95 | 11 | 112 | 218 |

10^{20} | 85 | 8 | 101 | 194 |

$ 4 \xd7 10 20$ | 70 | 6 | 115 | 191 |

10^{21} | 48 | 7 | 136 | 191 |

$ \u2265 4 \xd7 10 21$ | 33 | 8 | 149 | 190 |

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, $ \epsilon i , c$, while a gradual decrease in collisional loss, $ \epsilon 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 \u2265 4 \xd7 10 20$ m^{–3}. It can be seen from Table II that the decrease in $ \epsilon c$ cancels out with the increase in $ \epsilon i , c$. This is consistent with the observation in Fig. 7: the cathode fall increases and the plateau ( $ E \u2248 0$) region expands when the initial plasma density is increased from $ n p 0 = 10 20$ m^{–3} to $ n p 0 = 4 \xd7 10 21$ m^{–3}.

Figure 9 shows the profile of collisional energy loss due to the elastic collisions, $ S coll , el = n e \nu e l \Delta \epsilon el$, ionization, $ S coll , iz = n e \nu i z \Delta \epsilon iz$, and excitation, $ S coll , ex = n e \Sigma j \nu e x \Delta \epsilon ex$. For $ n p 0 \u2264 10 18 \u2009 m \u2212 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 \u2265 10 20 \u2009 m \u2212 3$, both ionization and excitation processes are concentrated around the cathode ( $ x / d \u2208 [ 0 , \u2009 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^{63}) from the rate coefficient point of view.^{31} In the quasineutral region, because $ T e \u2248 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 \u2248 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 \u2265 4 \xd7 10 21 \u2009 m \u2212 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), $ \epsilon 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, $ \epsilon 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*can be attributed mainly to the larger decrease in $ \epsilon c$ for $ 10 18 \u2009 m \u2212 3 < n p 0 < 10 20 \u2009 m \u2212 3$. Additionally,

_{s}*V*=

_{d}*V*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 \u2265 10 20 \u2009 m \u2212 3$. For $ n p 0 = 4 \xd7 10 21 \u2009 m \u2212 3$, as the quasineutral region reaches near the anode, the collisional loss due to excitation is minimized, while the cathode sheath potential $ \varphi s$ is close to

_{s}*V*(i.e., $ \varphi s \u2248 V d$). This results in the saturation of steady-state plasma profile at $ n p 0 \u2265 4 \xd7 10 21 \u2009 m \u2212 3$. From the perspective of the global energy balance,

_{d}^{11}the maximum plasma density near the anode sheath is also saturated at $ n p 0 \u2265 10 20 \u2009 m \u2212 3$ shown in Figs. 8(c)–8(e) when the total energy loss in the system $ \epsilon loss \u2248 \epsilon c + \varphi s$ is approximately constant.

## VI. MITIGATION OF HYSTERESIS DUE TO FIELD EMISSION IN MICROGAPS

For smaller gaps (e.g., $ d \u2248 1 \u2009 \mu m$), if the discharge voltage is on the order of 100 V, the electric field is on the other of 10^{8} 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 \u2264 10 19 \u2009 m \u2212 3$, which approximately corresponds to $ n p 0 / n g \u2264 10 \u2212 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 \u2009 m \u2212 3$ case): *pd* > 0.28 Pa · m, 0.12 Pa · m $ < p d \u2264 0.28$ Pa · m, and $ p d \u2264 0.12$ Pa · m. (I) At $ p d \u2265 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 \u2264 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 \u2264 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 \u2009 m \u2212 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.

## VII. CONCLUSION

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.^{24}

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 \u2248 10 \u2212 6 \u2013 10 \u2212 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*, steady-state plasma profiles obtained from the PIC-MCC model are analyzed. In both rarefied (left branch; small

_{s}*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*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 (

_{d}*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.^{24} 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}

## ACKNOWLEDGMENTS

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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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).

## DATA AVAILABILITY

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

## REFERENCES

*Gas Discharge Physics*

*Principles of Plasma Discharges and Materials Processing*