Self-pulsing of Dielectric Barrier Discharges at Low Driving Frequencies

This paper investigates the self-pulsing of Dielectric Barrier Discharges (DBDs) at low driving frequencies. In particular, (a) the dependence of current on the product pd of gas pressure p and the gas gap length d, (b) the effects of lossy dielectrics (in resistive discharges) and large dielectric permittivity (in ferroelectrics) on current dynamics, (c) the transition from Townsend to a dynamic Capacitively Coupled Plasma (CCP) discharge with changing pd values, and (d) the transition from Townsend to a high-frequency CCP regime with increasing the driving frequency. A one-dimensional fluid model of Argon plasma is coupled to an equivalent RC circuit for lossy dielectrics. Our results show multiple current pulses per AC period in Townsend and CCP discharge modes which are explained by uncoupled electron-ion transport in the absence of quasineutrality and surface charge deposition at dielectric interfaces. The number of current pulses decreases with an increasing applied frequency when the Townsend discharge transforms into the CCP discharge. The resistive barrier discharge with lossy dielectrics exhibits Townsend and glow modes for the same pd value (7.6 Torr cm) for higher and lower resistances, respectively. Finally,we show that ferroelectric materials can amplify discharge current in DBDs. Similarities between current pulsing in DBD, Trichel pulses in corona discharges, and subnormal oscillations in DC discharges are discussed. 1


I. INTRODUCTION
The Dielectric Barrier Discharge (DBD) is an alternating current (AC) discharge between electrodes covered by dielectrics 1 .They have certain properties of Capacitively Coupled Plasma (CCP) sources at low gas pressures and glow discharges at atmospheric gas pressures.DBDs find numerous applications for ozone generation 2 , decontaminating of biological samples 3 , thin-film deposition 4,5 , excimer emission 6,7 , cold plasma processing 8,9 , and many more 10,11 .
Depending upon the operating conditions, gas type, and pressure, a DBD can be spatially homogeneous or filamented.Also, they can operate in a glow or Townsend mode 12 .The glow mode, also known as an atmospheric pressure glow discharge (APGD), is distinguished by a higher discharge current reaching hundreds of milliamperes, non-uniform distribution of the electric field within the discharge gap, and plasma quasi-neutrality in the discharge bulk 13 .The glow discharge has been studied in air, nitrogen, helium, and other noble gases under atmospheric pressure and is usually characterized by one current pulse per half cycle [14][15][16] .APGDs were used for nanomaterial synthesis 17 , thin film deposition 18 , inactivation of bacteria 19 , and microplasma thrusters 20 .Conversely, the Townsend mode of discharge at atmospheric pressure, also known as an atmospheric pressure Townsend discharge (APTD), has a much lower discharge current of a few milliamperes, an almost uniform electric field along the discharge axis, and a lower electron density compared to the ion density 13 .In this mode, the electron and ion densities within the gap could differ substantially, and the Debye length is larger than the gas gap 13 .The Townsend discharge mode has applications for detecting volatile organic compounds 21 , ultraviolet light emission 22 , and thin film deposition 23 to name a few.
The trains of current pulses (or self-pulsation) in Townsend discharges have been observed in experiments and numerical simulations (see Ref. 24 for further references).The theory presented by Nikandrov and Tsendin, 25 explained the formation of current pulses in Townsend mode by deposition of surface charges on dielectric surfaces.The Townsend mode occurs when the ion drift time τ i through the discharge gap is much less than the inverse frequency of the applied voltage, i.e., where E br is the breakdown electric field, d is gas gap length, and µ i is the ion mobility.The opposite extreme, ωτ i > 1, corresponds to typical conditions of CCP operation with a plasma-sheath discharge structure 26 .
The self-pulsing of current in filament-free discharges was initially found by Bartnikas et   al. 27,28 while studying corona discharges in the 1960s.Akishev et al. [29][30][31] studied the self-pulsing regime of DBD numerically and corona discharges (both numerically and experimentally) and suggested that negative differential resistance is responsible for self-pulsation of current.Levko et al. 32 simulated the self-pulsing of subnormal DC discharges using a 2D fluid model and attributed the observed oscillation to ion transit time instability.Raizer et al. 33 studied the self-pulsing current in a Townsend DBD with semiconductor layers.They emphasized the dependence of the secondary emission coefficient on the electric field strength as a source of instability.Zhang et al. 34,35 studied the evolution mechanism of self-pulsing and the influence of impurities.
This paper aims to clarify the physics of self-pulsing of DBDs at low driving frequencies and its possible relation to subnormal oscillations in DC discharges 36 and Trichel pulses in corona discharges 37,38 .We illustrate that the self-pulsing of the filament-free DBD occurs because of the decoupling of electron and ion transport in the absence of quasineutrality in Townsend discharge.
We show that self-pulsing can also occur in CCP operating in a dynamic regime due to the absence of quasineutrality in the sheath.For the purpose of this paper, a one-dimensional fluid model of DBD is adequate.We solve the appropriate set of fluid equations using COMSOL 6.0 and attach an equivalent RC circuit to naked electrodes to model lossy dielectrics.We identify different regimes of discharge operation including quasi-DC (Townsend) discharges and CCP discharges operating in dynamic and high-frequency modes.
In this paper, we also discuss the dynamic regime of discharge operation.The dynamic regime is characterized by an almost constant (in time) value of plasma density in the gap but a substantial oscillation in electron temperature (mean energy) over an AC period.The dynamic regime was first introduced in Ref. 39 for the collisionless sheath operation in CCP.Later, it was identified for Inductively Coupled Plasma (ICP) at low driving frequencies 40 and more recently, for a positive column in AC discharged 41 and CCP 42 .The dynamic regime occurs when the applied frequency is between the inverse of the energy relaxation time and the ambipolar diffusion time (see Section III below for details).
The DBD characteristics depend also on the type of barrier materials used.The use of Resistive barrier discharge (RBD) and ferroelectric barrier discharge (FBD) have been discussed (see Ref. 11,43,44 ).The RBD uses high-resistive layers (aka a lossy dielectric), to prevent arcing.The FBD has ferroelectric layers with very high dielectric permittivity for applications in thin film deposition 45 , plasma actuators 46 , seed treatment, and plasma medicine 47 .
This paper is organized as follows.Section II describes the numerical model, which includes model geometry, fluid-Poisson model, boundary conditions, and reaction process.Section III discusses the characteristic scales and introduces a classification of AC discharges.Section IV contains our study's results, which are followed by the conclusions in Section V.

II. NUMERICAL MODEL A. Model Geometry
A schematic diagram of the DBD used in our simulations is shown in Figure 1.Two dielectric layers of thickness L = 0.1 mm and dielectric constants ε r are connected to a powered and the ground metal electrodes.An Ar gas at 1 atmospheric pressure (760 Torr) is placed between the two dielectrics of thickness L.An external sinusoidal voltage V (t) of amplitude V 0 is connected to two metal electrodes.There are a total of 268 grid points spanning between the electrodes, with 200 grid points within the gas gap.We solve the Fluid-Poisson model for 1D geometry.

B. Fluid-Poisson Model
For weakly ionized plasma, the time evolution of every plasma particle species p is obtained by solving the continuity equation 48 : where n p is the number density, Γ p is the density flux, S p is the source term, and index p denotes electron, ion, or excited atoms.
The density flux is given by where µ p is the mobility, q p is the particle charge, E is the electric field, and D p is the diffusion coefficient.For neutral particles, the first term on the right-hand side of the Eq. 3 is zero as q p = 0.
The source term S p is written in terms of a rate coefficient given by: where index r refers to a reaction type; c p,r is the net number of particles of species p created in the reaction of type r, and R r is the reaction rate for reaction r which is proportional to the densities of the reacting particles, i.e., R r = k f n n n p for two-body reactions 26 , n n is neutral gas density.
For electron-induced reactions, the rate coefficients have been computed from cross-section data by the following integral: where m is the electron mass, ε is electron kinetic energy, σ k is the collision cross-section, and f 0 is an electron energy distribution function (EEDF), which is assumed Maxwellian in the present paper.
The continuity equation for electron energy density (n ε = n e ε), is calculated from the energy balance equation 48 ∂ where n ε is the electron energy density, S ε is the energy loss or gain due to inelastic collisions, and Γ ε is electron energy flux calculated as 49 where µ ε and D ε are the electron energy mobility and the energy diffusion coefficients.The electron energy loss is obtained by summing the collisional energy loss overall reactions: where the first term on RHS represents heating by the electric field and the second term represents energy loss in collisions.Here, n r is the density of the target particles and the ∆ε r is the threshold energy.
The electrostatic potential V is calculated from the Poisson equation: where ε m is the permittivity of the medium.

C. Boundary Conditions
The boundary and surface conditions used in our simulations are described below.
The net surface charge σ s at dielectric surfaces is obtained from the particle fluxes as where j e and j i are total electron and ion current densities.
The boundary condition for the electric potential on a dielectric surface is where the subscripts 1 and 2 represent the gas gap and solid dielectric material, respectively.
The electron flux normal to the electrodes or walls, assuming no reflection, is given by where the electron thermal velocity is defined as v th,e = ( 8K B T e πm ), and γ is the secondary electron emission coefficient.
The normal component of the electron energy flux is given by where ε p is the mean energy of the secondary electron emitted by p th species.
For ions, and excited neutral species, which are lost to the wall due to surface reactions 48 , the flux is: The boundary conditions for the electric potential at the two electrodes are and where φ powered is the potential at the powered electrode and φ ground is the potential at the ground electrode.
The initial conditions used for the simulations are a gas temperature of 400 K, an electron density of 10 6 m −3 , a mean electron energy of 5 eV, and a mean energy of secondary electrons is 2.5 eV.We define the electron temperature as where ε denotes the electron mean energy 50 and 1 eV = 11,600 K.Note that the results of the simulation do not change with changes in the initial conditions such as changing the initial electron density to 10 6 m −3 from 10 5 m −3 .

D. Transport Coefficients and Reaction Processes
The set of volume and surface reactions considered in our simulations are shown in Table I and Table II.
The electron mobility is calculated by using the expression, µ e = e mν ea , and the ion mobility is given by µ i = 0.12 m 2 V s 1Torr p , taken from book 41,52 .

III. CHARACTERISTICS SCALES AND CLASSIFICATION OF AC DISCHARGES
The ion drift time across the gap is characterized by the time (τ i ) defined in Equation (1).Ion diffusion can usually be neglected because the ion thermal velocity is small compared to its drift velocity when the electric field is not close to zero.The time scale corresponding to the electron is much shorter because of the large electron mobility.
The time scale for electron diffusion across the gap is given by 52 where D e is the electron diffusion coefficient.
The characteristics time scales coming from the electron energy equation include the electron thermal diffusion time scale (τ T e ) which is of the order (τ D e ) and the energy relaxation time (τ ε ) which describes the electron cooling rate in collisions with neutrals.The latter is calculated by using the expression 53 where M is the mass of the atom.AC discharges can operate in different regimes depending on the value of the driving frequency with respect to the time scales introduced above.
In quasi-neutral plasma, the electric field forces the electrons and ions to diffuse at equal rates.
The charged particle dynamics are characterized by the slow ambipolar diffusion time scale τ a , defined as: where D a is the ambipolar diffusion coefficient.As a result, the plasma density does not change substantially during the AC period, but the electron temperature oscillates substantially over time.
The electron temperature is nearly uniform over space because the electron thermal diffusion occurs fast (at the free electron diffusion time scale) in both the sheath and plasma.
The dynamic regime of a collisionless sheath in CCP was first considered in detail by Nikandrov and Tsendin 39 for the frequency interval ωτ i < 1 < ωτ B , where . In our case of the collisional sheath, the motion of the electrons and ions is controlled by drift.Hence, the dynamic regime for the collisional sheath occurs in the frequency range ωτ b e < 1 < ωτ i .In the plasma, the dynamic regime corresponds to the frequency interval ωτ ε < 1 < ωτ a 40 .In CCP discharge with a plasma-sheath structure, both the sheath and the plasma can be in the dynamic regime at low driving frequencies.
By increasing ω, the high-frequency discharge regime corresponds to (ωτ ε > 1).In this regime, the electrons and ions respond to the average field in plasma, but electrons respond to the instantaneous field in the sheath.Electron temperature and the ionization rate are controlled by the time-average field in the plasma.
We performed simulations for different regions of Figure 2 by changing gas gap length (d), pressure (p), applied voltage (V), and driving frequencies (f ).The results of these simulations are presented in section sections IV.

IV. RESULTS AND DISCUSSIONS
We investigate the characteristics of an Argon DBD at atmospheric pressure and the transition from Townsend mode to CCP by changing pd.In addition, we also investigate the quasi-DC discharge, lossy DBD discharge, and ideal DBD discharge modes using the circuit equivalent of a resistive barrier discharge.Moreover, we study the transition from the Townsend mode of discharge to a dynamic CCP mode with a change in frequency in an ideal DBD.We use 1D simulations using the fluid-Poisson model in COMSOL and illustrate the physics of current multipulses and discharge characteristics.The details of each investigation are described below.

A. Oscillations of Current Density in the Townsend Mode
In this section, we discuss the oscillating behavior, i.e., the self-pulsing nature of current density (measured near the powered electrode) in the 1D DBD simulation of Ar gas.The discharge parameters are a frequency of 5 kHz, applied voltage amplitude 1000 V, gas gap distance 0.1 mm, secondary emission coefficient 0.05, and relative permittivity 5, i.e., ε r = 5.Note that the amplitude of the current density pulses decreases from the first to the sixth peak.The higher amplitude of the current density during pulses is because of the electron current and the lower value current density between the pulses comes from the ion current.Also, there is a phase shift for current density and applied voltage because of the capacitive coupling.Now we investigate the mode of discharge in our simulation (Townsend vs. glow mode).For this, we theoretically estimate the ion transit time using equation ( 1) and ion mobility (µ i ) as a function of E N in Td (1 Td = 1 Townsend = 10 −21 V m 2 ) following Basurto et.al 54 .The breakdown electric field required for equation ( 1) is given by E br = U br d , where d is gas gap distance and U br is the breakdown voltage, which is calculated from the simulation.The breakdown voltage is calculated as the difference between voltages measured at the dielectric-plasma interface node of the computational domain, which is 200 V for this case.The theoretically calculated value of the ion transit time is 24 µs.The ion transit time satisfies the condition, ωτ i < 1, and the discharge is the Townsend mode 25 .This is a quasi-static regime (quasi-stationary state) when the conduction current flows as a train of current pulses associated with the deposition of surface charges on dielectrics.
An interesting phenomenon that we observe in the Townsend discharge mode is the lowfrequency oscillation of the discharge current density (self-pulsing), as shown in Figure 3.These oscillations appear because of the time lag between ion formation near the instantaneous cathode (ground electrode for positive half cycle) and the current due to ion-electron emission.These pulses are a result of electron avalanches initiated by the background electrons and secondary electrons emitted from the instantaneous cathode due to ion bombardment.A gas breakdown occurs when the gap voltage exceeds the breakdown voltage.Rapid multiplication of seed electrons is linked to neutral gas species being ionized by electron impact and the growth of an electron avalanche.The newly generated electron-ion pairs move to the electrodes, driven by an electric field and diffusion.In the absence of plasma, the transport of electrons and ions occurs independently.Electrons drift and diffuse much faster than ions, and electron drift dominates over diffusion.As the breakdown voltage is much greater than the electron temperature, electrons drift to the instantaneous anode, and a net negative surface charge is deposited at the anode surface.
The gap voltage decreases below the breakdown voltage, and the electron multiplication process stops.The above process corresponds to the first current pulse in the discharge curve of Figure 3.This is how a current pulse is formed.At low driving frequencies (ωτ i < 1) ions have enough time to reach the instantaneous cathode during the half-period of the applied voltage increasing the electric field in the gap.The second current pulse is formed when the electric field again exceeds the breakdown value.Thus, the distance between the peaks is about the ion transit time in the Townsend discharge mode.Self-pulsing of DBD is similar to the self-pulsing of DC and corona discharge.It happens in the Townsend regime when the electron and ion transport is decoupled because the electric field is weakly perturbed by space charges 55 .We discuss further details of the discharge characteristics of this mode in Subsection IV B below.
In Figure 4, we present the validation of our model for helium discharge.As experimental data for current pulsing in filament-free DBD in the Townsend mode was unavailable, we conducted a simulation based on the experimental results of Shin et al., 56 and compared our results with the theoretical results presented by Nikandrove and Tsendin 25 .Though our model is oversimplified in terms of chemical and ionization kinetics and neglects photo-electron emission, we observe that our simulation results have reasonable agreement with both experiment and theory.We obtained four peaks at 500 Hz and ten peaks at 100 Hz.Additionally, we determined that the separation between peaks for 500 Hz applied frequency is approximately 2.8 µs, which is the same order as ion transition time, i.e., τ i = 3.3 µs.This validates our statement that the distance between the peaks is about the ion transition time in a Townsend discharge.Note that the ion transit time in Argon is larger than the ion transit time in Helium because of the mass difference.Therefore, achieving a Townsend regime in Argon is more difficult than in Helium.towards the anode.After reaching the peak value, the electron density starts to decrease over time.The increase in electron density is seen during all six current density pulses.Similarly, before the discharge, the ion density is uniform along the discharge axis and has a very small value as in the electron density.The ion density then swiftly increases during the discharge, near the instantaneous cathode, reaches a peak value, and then starts to decrease again.The increase in ion densities is seen during all six current density pulses.some oscillations along the discharge axis that are seen as spikes in the ion density profile near the instantaneous cathode.The ion density is around two orders of magnitude larger than the electron density during the pulse.This is a typical characteristic of the Townsend mode discharge where quasi-neutrality (n i = n e ) is not observed.This behavior strongly supports the conclusion that the discharge is in the Townsend mode.

Temporal Variation of Surface Charge Density
The surface charge accumulated at each plasma dielectric interface is calculated from Equation (11). Figure 8 illustrates the surface charge accumulation on the left (powered dielectric) and right (grounded dielectric) interfaces.The surface charge plays an important role in the gas discharge.This is because the presence of surface charge modifies the value of the gas gap voltage (V g ) by altering the overall voltage drop across the dielectric (V d ).The relationship between the gas gap voltage and surface charge density is given by 57 where V a is the applied voltage and V d is given by the relation where C d stands for the capacitance of each dielectric.The σ 1 and σ 2 are the surface charge densities on the left and right dielectrics, respectively.
It is evident from Figure 8 that before the beginning of the current pulses, the sum of the absolute values of the surface charge density on both dielectrics is a higher constant value, indicating that V d is larger and consequently V g is lower than the breakdown voltage.This situation prevents any discharge of the gas.As the surface charge density on both dielectrics begins to decrease over time, V d also decreases, leading to an increase in V g to a value larger than the breakdown voltage.Consequently, a breakdown in the gas gap is initiated.As discharge occurs, the surface charge density decreases first and then remains constant for a relatively longer time, causing V g to rise again, which in turn triggers the onset of the next discharge process as shown in Figure 8.

C. Effect of Discharge Parameters on DBD
Consider now the effect of different discharge parameters, such as pd values, driving frequency, and properties of dielectric materials on discharge dynamics.We also classify the discharge based on different characteristic scales, especially in the case of pd values and driving frequency.

Dependence of Discharge Dynamics on the pd Values
The dependence of the gas breakdown voltage on the product of pressure (p) and gap length (d), i.e., U br,pd = f (pd), is known as a Paschen curve 52 .The experimental Paschen curves for a static breakdown exhibit a minimum at a certain pd, which corresponds to optimal breakdown conditions.The Paschen curve for a static breakdown can be approximated as 26 , where γ is the secondary electron emission coefficient, and A and B are constants that describe the electron multiplication process and depend on the gas type.The experimental Paschen curve 52 for Argon is shown in Figure 9, illustrating that it has a minimum breakdown voltage of approx- velocity distribution is nonlocal in space and time, and electron "runaway" may occur 58 .Hence the fluid model does not work well near the Paschen minimum.Therefore, in this paper, we perform simulations for the pd values on the right branch of the Paschen curve.At large values of pd, ambipolar drift could become important, and a thermocurrent instability may occur 59 .
Figure 10 shows the simulation results for a gas gap distance of 3 mm (pd value of 228 Torr cm).Panel (a) illustrates the multi-pulse current density throughout the considered third cycle.
This case exhibits six current density pulses per half cycle.The spatial variation of the electron density (solid lines) and ion density (dashed lines) during five different phases of the half cycle of the considered cycle is shown in panel (b).As we can see, the electron density is equal to the ion density in most of the gas gap, except near the sheath region.The difference between the electron and ion density in the sheath region is more pronounced in the instantaneous cathode.
In the sheath region, the ion density is greater than the electron density, and a sharp increase in the electric field near the instantaneous cathode is observed, as shown in panel (c).This implies the existence of the cathode fall region, where the electric field increases almost linearly towards the cathode 60 .The quasineutrality of the plasma in the gas gap and the non-uniform nature of the electric field along the gas gap implies that the discharge is in the CCP glow mode.The spatiotemporal evolution of electron temperature over the considered cycle is shown in panel (d).
The electron temperature is high in the sheath because of the Joule heating 61 , i.e., the electric field in the sheath accelerates electrons, leading to increased kinetic energy and, consequently, higher temperatures.The temporal variation of the electron temperature at the center of the gas gap (r = 1.5 mm) for the considered first half cycle is shown in panel (e).Here, the electron temperature exhibits a significant oscillation over the half-period of the driving AC source 41 .Now we calculate the characteristic time scale to classify the discharge dynamics for the pd value of 228 Torr cm.We obtain that ωτ ε = 0.20 and ωτ a = 1.04 × 10 4 , see table III.These values satisfy the condition, ωτ ε < 1 < ωτ a , implying that the discharge is in the dynamic regime of CCP discharge.
Figure 11 shows the simulation results for a gas gap distance of 10 mm (pd value of 760 Torr cm).Panel (a) displays the current density for multiple pulses over the considered third cycle.In this case, six pulses are observed with relatively higher current density than the pd value of 228 Torr cm.The spatial variation of the electron density (solid lines) and ion density (dashed lines) during five different phases of the considered half cycle is shown in panel (b).As we can see, the electron density is equal to the ion density in most of the gas gap, except in the sheath region.As before, the difference between the electron and ion density in the sheath region is more pronounced in the instantaneous cathode.In the sheath, ion density is greater than electron density, and a sharp increase in the electric field near the instantaneous cathode is observed, as shown in panel (c).The sheath region, in this case, is even narrower than the previous case of pd value 228 Torr cm.The quasineutrality of the plasma in the gas gap and the non-uniform nature of the electric field imply that the discharge is in CCP glow mode.The current density multi-pulses are observed in this case as well.The spatiotemporal evolution of electron temperature over the cycle is shown on panel (d).As before, the electron temperature is high near the dielectrics instead of the center of the gap.The temporal variation of the electron temperature at the center of the gas gap (r = 5 mm) for the first half cycle of panel (a) is shown on panel (e).Again, the electron temperature exhibits a significant oscillation over the half-period of the driving AC source.In addition, the overall electron temperature is slightly lower than in the previous case.As before, the characteristic time scales, ωτ ε = 0.047 and ωτ a = 4.4 × 10 5 (see Table III) satisfy the condition 'ωτ ε < 1 < ωτ a ', implying that the discharge is in the dynamic regime of CCP discharge mode.The other discharge parameters are kept the same as the base-case simulation.
The physical reasoning for the dynamic regime of CCP glow mode for higher pd values (228 and 760 Torr cm) can be explained as follows.Due to the slow ion time scale, ions do not have a Note that all the plasma parameters: electron temperature, electron density, ion density, and Debye lengths are calculated at the gas gap center at the normalized time of τ = 0.25 of the considered cycle enough time to reach the cathode during the half-period of the applied voltage.They accumulate in the gap, and plasma forms during several AC periods.The ions and electrons are "glued" together by the electric field in the plasma.The electric field becomes inhomogeneous, and the typical plasma-sheath structure of a capacitively coupled plasma (CCP) is formed, with the potential drop in the plasma being about several electron volts.A slow, ambipolar time scale now characterizes the motion of electrons and ions in the plasma.However, electron heat transport does not depend on quasineutrality and occurs in the plasma and sheath at the same time scale as the free electron diffusion.Therefore, the plasma operates in a dynamic regime with small modulations of the plasma density and substantial modulations of the electron temperature and electron-induced ionization and chemical reactions.The current pulses are still observed in the dynamic regime of CCP operation as the processes in the anode sheath remain like those in Townsend discharges, i.e.
electrons and ions are decoupled in the sheath, and current pulsing may still occur because of this decoupling.

Dependence of Discharge Dynamics on Driving Frequency
To study the effect of applied frequency on discharge dynamics, we perform the simulations of a 1D DBD discharge in Ar for frequencies of 10 kHz and 25 MHz, keeping the other parameters the same as for the base case.These results are classified based on the characteristic scales of the different regions of Figure 2.  Figure 12 shows the results for a driving frequency of 10 kHz.Panel (a) shows the number of current pulses for the considered cycle and panel (b) shows the spatial variation of electron and ion density during five different phases of the half cycle of the AC period.The number of current density pulses in this case is four which is less than for the base case of 5 kHz (see Figure 3).This is expected because the number of current density multi pulses is proportional to ω −1/2 for the Townsend mode of discharge 25 .Also, we observe that the electron density during the first pulse is much less than the ion density and the plasma scale length (0.01 cm) is less than the Debye length (0.0454 cm).This supports the conclusion that the discharge is in the Townsend mode.
Figure 13 shows the results for a driving frequency of 25 MHz.Panel (a) shows the terminal This value is smaller than the previous 10 kHz case.In addition, the spatial variation of electron and ion density during the pulse shows that the electron and ion densities are equal during the pulse except at the sheath.The sheath region is most pronounced on the instantaneous cathode (right side).In addition, the plasma scale length (0.01 cm) is greater than the Debye length (0.0001 cm) and ωτ a > 1 implies that the discharge characteristic corresponding to the 25 MHz case is in a high-frequency regime.The characteristic scale sizes for both 10 kHz and 25 MHz cases are shown in Table IV.

Dependence of Pulses on Lossy Dielectrics in a Resistive Discharge
To study a lossy dielectric in resistive discharges, we used an equivalent circuit that emulates the resistive barrier discharge, as illustrated in where C is the capacitance of each capacitor in Farad, ε r = 5 is the relative permittivity of the dielectric, and ε 0 is the relative permittivity of the vacuum.The A and L are the area and width of each dielectric.We use the same dielectric on both sides, which gives C 1 = C 2 = 3.447 nF.When the resistance is very low (5 Ω), no current passes through the capacitors, and all the electron current passes through the resistor as depicted in Figure 15 (a).The lossy dielectric circuit in this case behaves as a naked electrode DBD.As we increase the resistance from 5 Ω to 5 kΩ and 50 kΩ, the electron current flowing through the resistor decreases while allowing the displacement current to pass through the capacitor (panels (b) and (c)).The circuit in this case behaves as a lossy dielectric.With a further increase in the resistance (to 500 kΩ), the electron current through the resistance ultimately reaches zero (panel (d)), and all the displacement current passes through the capacitor.The lossy dielectric circuit in this case behaves as a DBD with a perfect dielectric.The temporal variation of the current density for all four resistance values described above is shown in Figure 16.For the lower value of resistance (R = 5 Ω, panel (a)), the current density and applied voltage are in phase, and there are no current pulses in this case.This behavior is similar to a typical quasi-DC discharge.A similar glow mode of discharge at a low pd value of 1.27 Torr cm is obtained by Denpoh et al. 42 for a CCP discharge.For higher resistance (R = 5 kΩ, 50 kΩ, and 500 kΩ) in panels (b)-(d), the current leads the applied voltage by 90 • , because of CCP coupling and multiple pulses are observed in all three cases.The magnitude of the current density is slightly increased with the increase in resistance.In short, we can conclude that a dielectric material with a low conductivity (higher resistance) is required to produce discharge in DBD more easily.
To investigate the mode of discharge in lossy dielectrics in a resistive discharge, we plot the spatiotemporal evolution of the electric field (Figure 17) and spatial evolution of the electron and ion densities during the first pulse of the positive half cycle (Figure 18).For the lowest resistance (5 Ω) case, the electric field is not uniform along the gas gap and the electron (solid curve) and ion density (dotted curve) during the pulse are nearly equal (n e = n i ) in a gas gap, except within the sheath region.Moreover, the electron density is almost zero within the sheath region and the ion density is significantly higher than the electron density.This indicates that the discharge in this case is similar to a quasi-DC glow mode.For the higher value of resistance (5 kΩ, 50 kΩ, and 500 kΩ cases), the electric field is uniform in the gas gap.In addition, the electron density is less than the ion density in all three cases.This indicates that the discharge in all three cases of higher resistance is in the Townsend mode.

Current Pulses for High Permittivity Ferroelectrics
In this section, we study the dynamics of the atmospheric pressure barrier argon plasma using ferroelectric materials as a dielectric layer ferroelectric barrier discharge (FBD) 44 .Outstanding properties of ferroelectrics are high permittivity and a large value of the secondary electron emis- sion coefficient.For ferroelectric materials with relative permittivities from 25 to 1000, we observed a significant enhancement in the current density, even for the same value of γ as shown in Figure 19 ((a)-(e)).The higher discharge current density for a large value of permittivity is due to their stronger polarization charge 64 that facilitates a much higher charge accumulation onto the electrode surface as seen on panel (f) causing an enhancement in the discharge process and hence current density.Note that the ferroelectric with relative permittivity of 85.2 in Figure 19 (c) is lithium niobate (LiNbO 3 ) and this is close to the permittivity value of water (80.10).Hence the result obtained in this case could be similar to using water as liquid dielectrics.The temporal evolution of surface charge density on the powered electrode (left) for all three choices of permittivity in the ferroelectrics.Other simulation parameters are kept the same as in the base case.Note that the surface charge density is scaled by a factor of 0.2 for ε r = 1000 to highlight its greater magnitude within the Y-axis range.minimum).For the higher pd values of 228 and 760 Torr cm, the dynamic regime of the CCP has been observed with significant oscillations in electron temperature over an AC period.Our study Based on our analysis, self-pulsing of DBD appears similar to the subnormal oscillations in DC discharges and Trichel pulses in corona discharge.They all happen in the regime, where quasineutrality is absent and the electron and ion transport is decoupled because the electric field is weakly perturbed by space charges 55 .However, additional studies are desirable because the simplest fluid model with a Maxwellian EEDF used in our studies has limited applicability and accuracy.We expect that using fluid models with non-Maxwellian EEDF will not change our results quantitatively.The analysis of non-local electron kinetic effects on self-pulsing phenomena could be the subject of future studies.

FIG. 1 .
FIG. 1. Schematic diagram of a DBD cell.Two dielectric layers of permittivity ε r and thickness L = 0.1 mm (for both) are placed in between two metal electrodes (gray rectangle).Ar gas with a separation length d is placed between the two dielectrics at atmospheric pressure.A sinusoidal voltage source is connected to the electrodes.

FIG. 2 .
FIG. 2. A schematic diagram representing different characteristic scales used for the classification of AC discharges where τ b e the electron drift time scale, τ D e is the electron diffusion time scale, τ i ion time scale, and τ a is the ambipolar diffusion time scale.The abbreviation QSS is for Quasi-stationary state and HF is for High frequency.The self-pulsing is a common feature that we observed in dynamic and QSS regimes in our paper.

Figure 3 (
FIG. 3. (a) The time evolution of applied voltage (solid blue curve), gas gap voltage (dashed light blue curve), and terminal current density (red curve) over six cycles for Argon plasma at frequency 5 kHz, applied voltage amplitude 1000 V, and gas gap distance of 0.1 mm.(b) A single cycle over the time 0.45 − 0.65 ms showing current density multi-pulses in positive and negative half cycles with a breakdown voltage of approximately 500 V, and (c) time-resolved pulses of the first half cycle

Figure 5 (
Figure 5 (a) illustrates the spatiotemporal evolution of the electric field during the six cycles shown in Figure 3.The electric field appears to be constant throughout the gas gap over all six cycles.There is a change in the magnitude of the electric field from a positive to a negative value over every half cycle.Panel (b) represents the spatial distribution of the electric field during five different phases of the half cycle of the AC period from time 0.45 − 0.55 ms.The electric field remains constant throughout the gas gap for all five different phases.The uniform electric field over the gas gap distance indicates that the discharge is in Townsend mode.

Figure 6
Figure 6 displays the spatiotemporal evolution of the electron density (left panels) and ion density (right panels).The top panels show the electron and ion densities over the six cycles and the bottom panels show the same quantities over the half-cycle from 0.45 − 0.55 ms.Before the discharge, the electron density is uniform and has a very small value along the discharge axis (see panel (c)).During gas discharge, seed electrons become energized, which initiates an avalanche.The electron multiplication in the gap causes a gradual increase in the electron density

Figure 7
Figure 7 displays the spatial evolution of the electron and ion densities during the five different phases of the half cycle of the AC period from 0.45 − 0.55 ms.During each phase of half AC period, the electron density is maximum near the instantaneous anode (left side) and gradually decreases towards the instantaneous cathode (right side) whereas the ion density gradually increases towards the instantaneous cathode.In addition, the spatial ion density profile shows

FIG. 6 .
FIG. 6.The spatiotemporal evolution of the (a) electron and (b) ion density over the six cycles.The bottom panels show the spatiotemporal evolution of the (c) electron and (d) ion density over the half-cycle from 0.45 − 0.55 ms.

FIG. 7 .
FIG. 7. The spatial distribution of electron (solid curves) and ion densities (dashed curves) during the five different phases of the half cycle of the AC period from 0.45 − 0.65 ms.The left ordinate corresponds to the electron density and the right ordinate to the ion density.
FIG. 9.The Paschen curve for argon with A = 4.53, B = 78.33,and γ = 0.05 (solid blue curve).The Paschen curve with the same values of A, and B but with different electron emission coefficients (γ = 0.5) is also shown to illustrate how the Paschen curve can be modified based on γ.The experimental data (solid black curve) is extracted from the book 52 .The dashed vertical gray lines with increasing thickness represent the pd values for different gas gap distances at which simulations are performed.

FIG. 10 .
FIG. 10.Variation of the terminal current density showing multi-peaks (panel (a)), spatial variation of electron and ion density during five different phases of the half cycle of the considered third cycle (panel (b)), spatial variation of the electric field during five different phases of the half cycle of the AC period of the third cycle (panel (c)), the spatiotemporal evolution of electron temperature over one cycle (panel (d)), and spatial variation of the electron temperature over a half cycle (panel (e)) at 3 mm (pd = 228 Torr cm).

FIG. 11 .
FIG. 11.Variation of the terminal current density showing multi-peaks (panel (a)), spatial variation of electron and ion density during five different phases of the half cycle of the AC period from 0.45 − 0.65 ms (panel (b)), spatial variation of the electric field during five different phases of the half cycle of the AC period from 0.45 − 0.65 ms (panel (c)), the spatiotemporal evolution of electron temperature over one cycle (panel (d)) and spatial variation of the electron temperature (panel (e)) at 10 mm (pd = 760 Torr cm).The other discharge parameters are kept the same as the base-case simulation.

bFIG. 12 .
FIG. 12. Time evolution of current density multi-peaks (panel (a)), and spatial evolution of the electron (solid curve) and ion density (dashed curve) during five different phases of the half cycle of the AC period of the third cycle (panel (b)) at the applied frequency of 10 kHz.The other discharge parameters are kept the same as for the base case simulation.

FIG. 13 .
FIG. 13.Time evolution of current density multi-peaks (panel (a)), and spatial evolution of the electron (solid curve) and ion density (dashed curve) during five different phases of the half cycle of the AC period of the third cycle (panel (b)) at frequency 25 MHz.The other discharge parameters are kept the same as for the base case simulation.

Figure 14 .
FIG. 14.A schematic diagram of the circuit equivalent to the lossy dielectrics in the DBD simulations.Two RC parallel circuits are connected on each side of the gas gap to mimic the lossy dielectric 62 .

FIG. 16 .
FIG. 16.Time evolution of the applied voltage and terminal current density (measured near the plasma domain) for different values of resistance in the RC parallel circuit for lossy dielectrics: (a) 5 Ω, (b) 5 kΩ, (c) 50 kΩ, and (d) 500 kΩ.

FIG. 17 .
FIG. 17.The spatiotemporal evolution of the electric field at five different phases of the half cycle of the AC period of the third cycle for different values of resistance in lossy dielectrics: (a) 5 Ω, (b) 5 kΩ, (c) 50 kΩ, and (d) 500 kΩ.

FIG. 18 .FIG. 19 .
FIG. 18.The spatial distribution of the electron density (solid curve) and ion density (dotted curve) at five different phases of the half cycle of the AC period of the third cycle for different values of resistance in lossy dielectrics: (a) 5 Ω, (b) 5 kΩ, (c) 50 kΩ, and (d) 500 kΩ.

TABLE I .
Table of Collisions and Reactions Modeled (electron impact cross section are obtained from

TABLE III .
Table of plasma parameters for Gas Discharge at different pd values a .

TABLE IV .
Table of plasma parameters for Gas Discharge at different frequencies b