Mode Transition of Filaments in Packed-Bed Dielectric Barrier Discharges

We investigated the mode transition from volume to surface discharge in a packed bed dielectric barrier discharge reactor by a two-dimensional particle-in-cell/Monte Carlo collision method. The calculations are performed at atmospheric pressure for various driving voltages and for gas mixtures with different N2 and O2 compositions. Our results reveal that both a change of the driving voltage and gas mixture can induce mode transition. Upon increasing voltage, a mode transition from hybrid (volume+surface) discharge to pure surface discharge occurs, because the charged species can escape much more easily to the beads and charge the bead surface due to the strong electric field at high driving voltage. This significant surface charging will further enhance the tangential component of the electric field along the dielectric bead surface, yielding surface ionization waves (SIWs). The SIWs will give rise to a high concentration of reactive species on the surface, and thus possibly enhance the surface activity of the beads, which might be of interest for plasma catalysis. Indeed, electron impact excitation and ionization mainly take place near the bead surface. In addition, the propagation speed of SIWs becomes faster with increasing N2 content in the gas mixture, and slower with increasing O2 content, due to the loss of electrons by attachment to O2 molecules. Indeed, the negative O− 2 ion density produced by electron impact attachment is much higher than the electron and positive O 2 ion density. The different ionization rates between N2 and O2 gases will create different amounts of electrons and ions on the dielectric bead surface, which might also have effects in plasma catalysis.


Introduction
Plasma catalysis is gaining increasing interest for various environmental applications, such as gaseous pollutant removal, the splitting of CO 2 , hydrogen generation and O 3 production [1][2][3][4][5][6][7][8][9].Plasma catalysis can be regarded as the combination of a plasma with a catalyst, and often results in improved performance, in terms of selectivity and energy efficiency of the process.Plasma is an ionized gas, consisting of various reactive species, like electrons, positive ions, negative ions and radicals.These reactive species are created by applying a potential difference to a gas.The gas itself remains at room temperature, which is beneficial in terms of energy saving compared to classical thermal catalysis.Plasma-based pollutant removal and gas conversion have recently gained increased attention, being possible alternatives for chemical reactors.The production of new molecules in the plasma can be selective when a catalyst is added to the plasma.Plasma catalysis can be realized by introducing dielectric packing beads (coated with catalyst material) in the discharge gap, forming a packed bed dielectric barrier discharge (PB-DBD) reactor.The DBD generally occurs in a filamentary mode by applying a high driving voltage [10][11][12][13][14], which induces a very fast ionization avalanche, propagating from powered electrode to grounded electrode, i.e., a so-called streamer [12,13,[15][16][17][18][19].Each streamer starts when the driving voltage passes a certain threshold, and will further polarize the dielectric surface [20].
The ionization avalanches can occur either within the bulk plasma (so-called volume discharge) or along the surface of a dielectric (so-called surface discharge), sustained by electron impact surface ionization waves (SIWs) [21].SIW discharges operating in N 2 at low pressure (5-100 Torr) and created by high voltage (10-15 kV) nanosecond pulses, traveling along a dielectric surface or a liquid surface, have been experimentally investigated by Intensified Charge Coupled Device (ICCD) camera images [22], with measured propagation speeds of ∼5 × 10 5 m/s.Furthermore, a SIW discharge in hydrogen was experimentally reported, by measuring the time-resolved electric field at special positions, within tens of nanosecond pulse scales, using a picosecond four-wave mixing technique [23].The results showed that the discharge developed itself as a SIW propagating along the dielectric surface at an average speed of ∼10 6 m/s with a maximum electric field of ∼2.3 × 10 6 V/m [23].
On the other hand, volume discharges can be sustained by filamentary microdischarges (MDs) in a PB-DBD reactor, by limiting the dimensions at atmospheric pressure, following the famous Paschen law [24][25][26].A filamentary MD may be tens of micrometers in space and can exist for a few nanoseconds [27], yielding a high concentration of reactive species in a narrow gap, which is beneficial for pollutant remediation.Furthermore, filamentary MDs can co-exist with SIWs along the surface of the dielectric beads in a PB-DBD reactor under proper conditions [28,29].This may yield high concentrations of chemically active species and radicals on the surface due to electric field enhancement along the dielectric surface, compared to an unpacked DBD.
The generation of reactive species at catalyst surfaces and in the gaps between the packing beads is very important for plasma catalysis.Therefore, it is crucial to understand the properties of surface and volume discharges and their formation mechanisms.However, only a few theoretical studies have been performed for PB-DBD reactors [11,[29][30][31][32][33][34][35][36][37].Russ et al. [30] adopted a two-dimensional (2D) hydrodynamic theory to study transient MDs in a PB-DBD reactor for a gas mixture 80% N 2 , 20% O 2 , and 500 ppm NO, with 23 species and 120 plasma reactions.However, this work was limited to a one-directional discharge, without considering discharge channels in the voids of the packed beads.Kang et al. [31] implemented a 2D fluid model to investigate the typical properties of MDs in a ferroelectric PB-DBD, focusing only on the surface discharge with strong electric field for two packing beads.Chang et al. [32] studied a ferroelectric packed bed N 2 plasma for a spherical void between two pellets based on the 1D Poisson equation and transport equations.This work focused on the electron density, electron temperature, and electric field as a function of the applied voltage, the discharge gap size, and the dielectric constants of the packing beads.
In our previous work, we employed a kinetic 2D particle-in-cell/Monte Carlo collision (PIC/MCC) model to investigate the formation and propagation of a filamentary MD in a PB-DBD, and the results were compared to an unpacked DBD [11].Van Laer and Bogaerts developed a comprehensive 2D fluid model for a PB-DBD in helium [33][34][35].They reported that a packing enhanced the electric field strength and electron energy near the contact points of the dielectric beads.Microplasma discharges in humid air in a PB-DBD reactor were experimentally characterized through ICCDimaging and numerically simulated by the 2D multi-fluid simulator non PDPSIM, with a negative driving voltage at atmospheric pressure [29].The behavior of three kinds of discharge modes, including positive streamers (restrikes), filamentary MDs and SIWs, were reported in their work.Wang et al. [36] studied the microplasma characteristics by means of fluid modelling and experimental observations using ICCD imaging, for packing materials with different dielectric constants, in dry air at atmospheric pressure.They also studied the same three types of discharges predicted in Ref. [29].In addition, they reported a transition in discharge mode from surface discharges to local filamentary discharges upon increasing the dielectric constant of the packing from 5 to 1000.Finally, Kang et al. [37] experimentally and numerically investigated surface streamer propagation on an alumina bead, and predicted three distinct phases, i.e., the generation and propagation of a primary streamer with a moderate electric field and velocity, rapid primary streamer acceleration with an enhanced electric field, and slow secondary streamer propagation.
Although several studies have been performed to better understand the microplasma properties in a PB-DBD reactor, as outlined above, the interaction mechanisms between the plasma and the catalyst are still poorly understood.Both physical and chemical effects can play a role in this interaction, and this will affect the discharge types and properties.Therefore, in the present work, we investigate the filament formation and mode transition between volume and surface discharges, as a function of driving voltage and O 2 /N 2 mixing ratio at atmospheric pressure, by using a 2D PIC/MCC model in a micro gap PB-DBD reactor.Specifically, we aim to obtain a better understanding of the production mechanisms and concentrations of reactive species sustained by the different modes.

Results and Discussion
The spatial and temporal evolutions of the species densities, electric field, and excitation and ionization rates are presented in this section, including the whole process of filament formation and transition to SIWs.

Effect of Driving Voltage
It has been reported that the SIW properties are dependent on both the direction and value of the driving voltage and the surface properties [21].In Figures 1-7, we present the effect of driving voltage on the density of the reactive species and the electric field for a dry air discharge (so only a mixture of 80% N 2 and 20% O 2 , no other air components).Figure 1 shows the electron density for four different driving voltages, i.e., −5 kV (a1-a4), −10 kV (b1-b4), −20 kV (c1-c4) and −30 kV (d1-d4), and four different time intervals.Note that the different time intervals at various driving voltages are based on the filament formation times under the different conditions.The driving voltages are all above the breakdown voltage of air (∼3 kV/mm) in glow mode.Since a filamentary discharge needs higher breakdown voltage than a glow discharge, no filament occurs at −5 kV, so only the seed electron density (∼10 18 m −3 ) is shown in Figure 1a1-a4, which corresponds with the values in Ref. [29].A limited volumetrically sustained filament appears at −10 kV, with an average speed of 3 × 10 5 m/s and a maximum electric field of 8.6 × 10 7 V/m.For −20 kV, the filament is formed in the gap at around 0.25 ns with an average speed of 2 × 10 6 m/s.Subsequently, the SIWs are developed and start dominating the discharge.The SIWs can be further classified as downward-and upward-propagating modes.The downward-directed SIWs are formed in the time interval of 0.25-0.6ns with an average velocity of 2.2 × 10 6 m/s, while the upward-directed SIWs are developed in the time range of 0.6-0.8ns with an average velocity of 2 × 10 6 m/s.It is worth to note that a multi-channel filamentary MD also develops with the upward-directed SIWs, i.e., the plasma extends away from the dielectric surface.However, the maximum density in the MD is nearly two orders of magnitude lower than in the SIWs, as shown in Figure 1c4, which indicates that the discharge is mainly governed by the SIWs.For the driving voltage of −30 kV, the filament is formed and sustained in the time range of 0-0.2 ns, as nearly pure SIWs propagating downward along the surface with a high speed of 5 × 10 6 m/s, while no upward-directed SIW develops.There is obvious asymmetry in the electron density for −30 kV, due to the inherent statistical character of the streamers.
Therefore, as the driving voltage increases, a mode transition will happen, i.e., the filament is mainly sustained by an MD at low driving voltage (−10 kV), it becomes a combination of a limited portion of MD and significant SIW at a moderate driving voltage (−20 kV), and finally it is almost completely dominated by SIWs.The physical mechanism is that the charged species will rapidly escape to the dielectric surface under the effect of high driving voltage (strong electric field) and charge the dielectric.The surface charging will induce a strong tangential component of electric field along the dielectric surface, and lead to the formation of SIWs.The SIWs will in turn yield high concentrations of reactive species, as well as high electron density on the surface, which will further charge the dielectric, and gradually form a filament along the dielectric surface.This filament may increase the surface activity of the beads, which is probably beneficial for plasma catalysis.
Our results are consistent with the theoretical and experimental predictions in Refs.[22,23,29,36].For −10 kV, the average speed of the volume streamer is 3 × 10 5 m/s, which is in the same order as in Refs.(∼5 × 10 5 and ∼10 5 -10 6 m/s) [29,36].For −20 kV, both volume and surface filaments co-exist.Therefore, the filament could propagate in the gaps and on the surface of the dielectric beads.This trend is in good agreement with recent experimental observations by fast-camera imaging for numerous packing beads with dielectric constant of 25.Once the SIWs are initiated in the PB-DBD, we calculate a propagation speed of the SIWs around 10 6 m/s, which is in the same order as in literature (∼5 × 10 5 and ∼10 6 m/s) [22,23].The reactive species are highly concentrated on the surface due to the presence of SIWs along the catalyst surface, which agrees with experimental observations in literature [22,36].
Figures 2-4 show the density of N + 2 , O + 2 and O − 2 ions, respectively, for the same conditions as in Figure 2.
As shown in Figure 2, the N + 2 ion density profile is a bit different from the electron density distribution at −20 kV, i.e., with a higher density in the voids between the beads.The maximum N + 2 ion density is a bit higher than the electron density, due to loss of electrons by attachment to O 2 molecules.However, the N + 2 ion density profile becomes very similar to the electron density distribution at −30 kV, because the discharge is now dominated by pure surface discharges.This indicates that the electron attachment is not significant for a pure surface discharge.Figure 3 clearly shows that the O + 2 ions are mainly formed along the dielectric surface by the SIWs, and the maximum density is about 20 times lower than the maximum electron density for all driving voltages, because it is more difficult for the discharge to become sustained in O 2 gas [38].This can be justified by the ionization rate shown later in Figure 5.   Therefore, in total, the electrons and N + 2 ions are present both along the dielectric surface and in the gaps between the beads, whereas the O + 2 and O − 2 ions only concentrate along the dielectric surface at −20 kV.The difference in spatial distributions of charged reactive species may affect the catalytic selectivity, which is important for plasma catalysis, as the interaction between the charged species and the surface might influence the morphology and work function of the catalyst, and accordingly affect the catalyst performance.However, the exact influence cannot be deduced from our model and can be quite complicated, which is beyond the scope of this work.
It has been reported [11,33,36,39] that the local electric field can be enhanced near the contact points and at the boundaries between beads and dielectric layers.This is because the dielectric beads and plates are strongly polarized by the applied electric field, which reduces the potential drop in the dielectric, but increases the potential drop in the gas gap.However, most of the previous works only focus on local electric field in the gas gap [29,40], but they do not study the variation and influence of the electric field in the dielectric materials.
Figure 6 shows the electric field amplitudes for the same conditions as in Figure 1, and Figure 7 shows the vertical and horizontal electric field components for −20 kV (at 0.8 ns) and −30 kV (at 0.2 ns).As clearly seen from Figures 6 and 7, the presence of dielectric beads and plates induces a sharp boundary with the gas phase, which strongly enhances the displacement electric field inside the dielectric material near this boundary.As a result, the electric field on the surface and in the voids of the packing beads is also enhanced, due to the significant charge accumulation in the narrow gaps (see Figures 1-4 above).Indeed, the local electric field inside the vertical gaps between two beads (or between bead and plate) is enhanced so much that the maximum value is much larger than the average external electric field, even for the lowest driving voltage of −5 kV without breakdown.The maximum values, presented in Figure 6, stay constant when there is no breakdown, but they increase with time when there is a discharge.
The electric field is more enhanced in the vertical gaps between the dielectric beads than inside the materials for the −10 kV case, due to the very limited filamentary MD, as displayed in Figure 1b.On the other hand, at −20 kV and −30 kV, the electric field on and near the surface of the dielectric materials is significantly enhanced, due to the high concentration of the charged species on the catalyst surfaces (see Figures 1-4 above).In addition, as shown in Figure 6c4, the electric field is also obviously enhanced in the head of the upward-directed SIWs (noted in Figure 1c4 above), where the electric field is characterized by an overlap of both vertical and horizontal electric field components.Behind the SIW head, the vertical component rapidly reduces, as shown in Figure 7a1.Moreover, at −30 kV, the electric field has its maximum values on the surface of the beads and plates, and the values gradually decrease away from the surface boundary, because of the pure surface discharge mode.The high electric field on the dielectric bead surface is important for plasma catalysis [8], because it could induce locally stronger electron impact reactions, thus higher reaction rates and a higher plasma density.The SIW discharges therefore play an important role in a PB-DBD reactor for plasma catalysis applications.The physical mechanism is related to the high concentration of reactive species on the surface of the beads, which yields tangential components of the electric field, and the latter in turn enhance the SIWs.
We calculated the tangential components of the electric field E x to be around 2 × 10 7 V/m, as shown in Figure 7b2, which is larger than the value of 6 × 10 6 V/m from literature [29].The difference is mainly due to the much smaller discharge gap of 1 mm, and bead size of 508 µm (bead diameter) used in this work, compared to the large discharge gap of 1 cm and bead size of 1.8 mm (bead diameter) in Ref. [29], as the same driving voltage will induce a much stronger electric field in a smaller gap.In addition, we consider ZrO 2 dielectric beads with a dielectric constant of 22, which yields stronger polarization than for quartz material with a dielectric constant of 4, used in Ref. [29].This stronger polarization also induces a higher electric field in the gap.
In order to elucidate the mechanisms giving rise to the SIWs on the dielectric surfaces, we need to distinguish where the reactive species are generated, either on the dielectric surfaces or in the gap.To clarify this, the electron impact excitation rate (a1-b1) and ionization rate (a2-b2), for N 2 (a1-a2) and O 2 (b1-b2) components in the mixture, are plotted in Figure 8, for 0.8 ns and a driving voltage of −20 kV.
As seen in Figure 5, both excitation and ionization mainly take place on the surface of the beads.The maximum excitation rate is about one order of magnitude higher than the maximum ionization rate, due to the lower threshold energy needed for excitation reactions (see Table 1 below).The excitation rate for the N 2 component is almost uniformly distributed on all surfaces, and fills the gaps between beads 3, 4 and 5, indicating that the excited N * 2 molecules will be fairly uniformly distributed on these surfaces and in the gaps.The ionization rate for the N 2 component is locally enhanced on the top surface of dielectric bead 4 and the bottom dielectric plate, due to the high local electric field there, which corresponds to the higher N + 2 ion density at these positions, shown in Figure 2c4.Although the excitation and ionization rates for the N 2 component spread out from the surface and cause restrikes [29], which may further give rise to multi-channel filamentary MD (see Figure 1c4 above), the maximum values in the restrikes are about one order of magnitude smaller than on the surface, as shown in Figure 5a1-a2.This again demonstrates that the SIWs are the main discharge mode at the driving voltage of −20 kV, while the MDs sustained by column-like filaments are more or less negligible.
Table 1.Reaction set of a N 2 /O 2 gas mixture used in the model, as well as the threshold energy for electron impact ionization and excitation collisions.The cross sections for all reactions are adopted from Refs.[41][42][43][44], and are downloaded from the LXCat database [45].

Reactions Threshold Energy (eV)
Electron On the other hand, both the excitation and ionization rates for the O 2 component show discontinuous distributions on the bead surfaces, except for the top surface of central bead 4, as shown in Figure 5b1-b2.This indicates that the production of the excited O * 2 molecules and the O + 2 ions is not a continuous process, but rather a discrete process, as predicted in Ref. [29].This discontinuous ionization rate leads to a rather discretely distributed O + 2 ion density (see Figure 3 above).The difference between the ionization rate for N 2 and O 2 gas components may be attributed to the different mean free path for ionization collisions.The mean free path for ionization collisions with N 2 is about 250 µm (hence similar to the radius of the beads, i.e., 254 µm), while it is ∼2 mm (similar with the perimeter of the beads) for O 2 .A shorter ionization mean free path results in more intensive ionization collisions, whereas an ionization mean free path longer than the bead radius and the gap distance between beads results in weaker ionization and may yield a discrete distribution, as shown in Figure 5b2.Although the densities of the excited species (N * 2 and O * 2 ) are not displayed in this work, the electron impact excitation rate reveals that these excited species are highly concentrated on the surface.The electron impact ionization rate profile can be considered as evidence for SIW formation on the surface.

Effect of the Gas Composition
We vary the O 2 /N 2 mixing ratio in this subsection, to obtain better insight in the different character between electropositive and electronegative gases.The effect of gas composition on the production of MDs inside porous ceramics was investigated experimentally by Hensel et al. [38], using photographic visualization and optical emission spectroscopy.It was predicted that a higher O 2 content resulted in the redistribution of the MD channels inside the porousceramics [38], while the total number of MD channels was reduced, since the breakdown voltage was increased.Their results can be partially correlated to the excitation and ionization rates, which are stronger in N 2 than in O 2 , as shown in Figure 5 above.Therefore, it is important to understand the effect of gas composition on the mechanisms of filament formation and mode transition between volume and surface filaments in a PB-DBD reactor, which has not been reported before.
In this section, we will demonstrate that the gas composition is a critical parameter for mode transition.Figure 8 shows the electron density.The filament is initiated from t = 0 with a few seed electrons.The filaments travel for a similar distance within 0.2 ns, and arrive at the central bead 4 at 0.25 ns, for the different gas compositions, except in the case of 10% O 2 and 100% O 2 .However, after the filaments reach the central bead, their propagation speed becomes faster in N 2 gas and gradually becomes slower with increasing O 2 contents, as shown in Figure 8 at 0.4 and 0.7 ns.Indeed, the average speed of the filaments is 2 × 10 6 m/s during the first 0.25 ns, while during 0.25-0.4ns, the average speed is 3 × 10 6 m/s for pure N 2 gas, indicating the fastest SIW, 2.3 × 10 6 m/s for 10% O 2 , 2 × 10 6 m/s for 50% O 2 , and 1.7 × 10 6 m/s for pure O 2 .
The physical mechanism is that the discharge is more difficult to be sustained, and thus it needs a higher breakdown voltage and leads to a slower discharge evolution, for electronegative gas, due to the loss of electrons by attachment to O 2 molecules.This feature was experimentally confirmed by Hensel et al. [38], through measurements of the breakdown voltage in a N 2 /O 2 gas mixture.Therefore, for the same driving voltage, the discharge can more easily be created in N 2 than in O 2 , resulting in gradually slower propagation speeds of the filaments upon higher O 2 fractions in the mixture.
For pure N 2 gas, our calculations predict a dominated surface discharge with a negligible volume discharge.Indeed, the maximum electron density (1.51 × 10 24 m −3 ) on the surface is approximately five orders of magnitude higher than in the gap (3.51 × 10 19 m −3 ).The electron density is almost uniformly distributed on the bead surfaces after the filament reaches the central bead 4.
The maximum electron density on the bead surfaces is represented by the red color in Figure 8.It is smooth for the first two cases, indicating that the electron density is indeed uniformly distributed, while it becomes discontinuous with increasing O 2 contents.However, the maximum electron density still occurs on the surface and not in the gap, giving rise to a high concentration of electrons on a specific part of the surface in pure O 2 .Again it shows that the SIWs can produce a high density of charged species on the surface of the beads, which might be beneficial for gas treatment by plasma catalysis [46].Finally, while adding more O 2 to the mixture makes the SIWs to become discontinuous, the gas mixture almost does not affect the volume discharge, as seen in Figure 8.

Smooth Smooth discontinuous
Figure 9. N + 2 ion density (m −3 ) for the same conditions as in Figure 8.

Figure 10. O +
2 ion density (m −3 ) for the same conditions as in Figure 8.
The N + 2 ion density profile is similar to the electron density profile, but the maximum N + 2 ion density is two times larger than the electron density, again due to the loss of electrons by attachment to the O 2 molecules.The number of lost electrons can be determined by the O − 2 ion density (see Figure 11  below).The maximum N + 2 ion density (2.91 × 10 24 m −3 ) is always found on the surfaces, and the speed of the N + 2 ion filament is almost the same as for the electron filament discussed above.

Figure 11. O −
2 ion density (m −3 ) for the same conditions as in Figure 8.
The O + 2 ion density exhibits different profiles with increasing O 2 content, as shown in Figure 10.The O + 2 ion density is at least two times smaller than the electron density, for the same reason as above, i.e., it becomes harder for the discharge to be created for high O 2 content at a fixed driving voltage, resulting in a lower O + 2 ion density.As a result, the O + 2 ion density is at least four times smaller than the N + 2 ion density at equal gas fractions (50%/50% mixture).In pure O 2 , the propagation speed of the SIW is about two times slower than in pure N 2 , and no MD channel is formed above the surface of beads 3 and 5, due to the absence of an upward-directed SIW.
The O − 2 ions are again mainly formed on the surface by the SIWs, with a maximum value of 4.01 × 10 24 m −3 in 100% O 2 , which is 2.6 times higher than the maximum electron density and 5.8 times higher than the maximum O + 2 ion density, indicating that electron attachment is quite significant.Indeed, attachment is easier than ionization, because it requires no threshold energy, while the ionization threshold is 12.06 eV for O 2 and 15.58 eV for N 2 .Therefore, the O − 2 ion density sharply increases with rising O 2 content, as expected.
To summarize, the O + 2 and O − 2 ion densities are enhanced, while the electron and N + 2 ion densities drop on the bead surface, for higher O 2 content, as expected.This might affect the surface reactions in plasma catalysis [47,48].

Model Assumptions
The dimensions of the whole reactor are 1.65 mm in the y direction and 10 mm in the x direction.The discharge is sustained between two parallel plate electrodes covered by two dielectric layers of 0.3 mm thickness, separated by a gap of 1 mm in the y direction.Dielectric beads are inserted in the gap, forming a packed bed reactor.The dielectric constant of the layers and the beads is r = 22, characteristic for ZrO 2 .A schematic illustration of the reactor is presented in Figure 12, showing in the x direction only the central part including the dielectric beads, but the entire simulation region is taken from 0 to 10 mm in the x dimension.The dielectric layers and beads are colored in blue.We consider five dielectric spheres with diameter d = 508 µm.In order to leave some gas space for the filament formation between the packing beads, these five beads are packed in a non-strict spherical packing manner [11], with distance between adjacent bead centers of (0.1 + √ 3/2)d in the y dimension and 1.1d in the x dimension.Note that some gaps are present between the dielectric beads in Figure 12.This was intentional to show filament formation.Indeed, streamers cannot propagate in a 2D system without gaps between the beads, while they can propagate either along the bead surfaces or the gaps between the beads in a 3D case.We thus assume a certain gap between the beads, to allow the streamers propagate in our 2D model.Furthermore, even in experimental (3D) packed-bed reactors, there are some gaps between the dielectric beads, i.e., the cross sections between the beads are similar with the 2D model, when a lot of beads are packed together, as the beads cannot touch each other perfectly.This can be observed from Ref. [49].On the other hand, it needs extremely long calculation times to model a PB-DBD with 5 beads in an entire 3D geometry, owing to severe mesh requirements and very large number of macro particles.Therefore, we consider a 2D model.The upper electrode (y = 1.65 mm) is powered by a pulsed voltage with a rise time of 0.1 ns and then kept constant for the duration of the simulation, considered as the cathode.The lower electrode (anode; y = 0 mm) is grounded.The dielectric surfaces are considered as absorption boundaries for the reactive plasma species, i.e., all reactive species will be removed from the simulation if they hit the surface of the dielectrics and they cannot participate in the discharge anymore.The absorbing boundary condition is a reasonable and widely used boundary condition in PIC model [50].When charged species are absorbed on the dielectrics, they can emit a secondary electron and they will deposit their charge on the dielectric surface.The deposited charge is determined by the charge of the ion striking the dielectric surface and the secondary electron emission charge and coefficient, accounting for the formalism Here, Q D is the deposited charge with initial value of zero, Q ion is the charge of the ions striking the dielectric surfaces of the beads or layers, and Q se is the charge of the secondary electron, respectively.The effect of secondary electron emission is self-consistently coupled in the PIC/MCC model, assuming a constant ion impact secondary electron emission coefficient of 0.15 [29].
Photoionization is neglected in this study.Indeed, we found in our previous work [39] that the results are nearly the same with and without photoionization in a gap of tens of µm, which is the typical gap size between the packing beads and between the dielectric layers and the beads in a PB-DBD.Furthermore, the photoionization rate was nearly two or three orders of magnitude lower than the electron impact ionization rate, even in a large gap (hundreds of µm ) [29,36].Indeed, as will be shown below, the filaments in the present model are mainly sustained by SIWs, and the volume discharges will be negligible in a narrow gap of ∼10 µm.
The simulations are applied to a mixture of N 2 and O 2 at atmospheric pressure with a temperature of 300 K.Only charged species, i.e., electrons (e) and ions (N + 2 , O + 2 and O − 2 ), are simulated in this PIC/MCC model, as the total ionization degree is typically less than 10%.Filaments are initiated from initial seeding electrons emitted from the surface of the top dielectric plate in the region x ∈ [4.9, 5.1] mm.A few seed electrons can be generated by cosmic radiation or external emission.The emission current density is set as 10 5 Am −2 .Subsequently, the filament sustains itself through anode-directed avalanches due to the electron impact ionization and secondary electron emission.Dissociation process, and the behavior of atomic ions, are not included in the model, because the cross sections of these are relatively small, thus omitting these reactions will almost not affect the kinetics of the streamer.

Simulation Method
We developed a 2D PIC/MCC model based on the VSIM simulation software [51], which has been widely used and validated [11,39,51].This PIC/MCC simulation is based on an explicit and electrostatic method, which was first introduced and described in detail in Ref. [52].The PIC/MCC model takes advantage of accounting for the detailed kinetic behavior of charged particles, compared to a fluid model.There are four steps in a PIC/MCC model: (1) pushing the particle velocity based on the previous electric field; (2) weighting from the positions of all charged particles to obtain the charge density; (3) summing and extrapolating the charge density to achieve a new electric field based on the Poisson equation; and (4) using a standard MCC method to describe electron impact elastic collisions, excitations, ionizations and electron attachment.
The particle pushing procedure is described by the Newton equations where Here n represents the nth time step (δt).The new electric field in time t n+1 is solved by the Poisson where = r 0 is the permittivity, and ρ n+1 is the total charge density, given by ρ n+1 = ∑ p Z p eS(r − r n+1 p ).
with the light speed c, and the space steps in x and y directions δx and δy.As mentioned above, the simulation region is 1.65 mm in the y direction and 10 mm in the x direction, and we use a mesh of 1000 × 500 uniform grid points.Dirichlet boundary conditions are adopted in the y direction, and Neumann conditions are used in the x direction.
A standard MCC method [53] is used to account for the electron impact elastic, excitation, ionization and attachment collisions with N 2 and O 2 gas molecules, as listed in Table 1.The cross sections and threshold energies used for these reactions are adopted from the Refs.[41][42][43][44] and downloaded from LXCat database [45].
We consider only three types of ions, since they are the most important ones, with the lowest ionization threshold and largest density compared to other ions.We consider a simulation time up to 0.8 ns, which is enough to develop both the volume and surface discharges in a PB-DBD reactor.Thus, the effect of electron-ion and negative ion-positive ion recombination is negligible.Indeed the recombination reactions have a larger relaxation time up to the microsecond time scale, as predicted by Kruszelnicki et al. [29], as it requires a sequence of two-body reactions or three body reactions.

Conclusions
We have applied a 2D PIC/MCC model to study the formation and mode transition of filamentary discharges for a PB-DBD reactor with various driving voltages and N 2 /O 2 gas mixtures at atmospheric pressure.The discharge is sustained between two parallel plate electrodes covered by two dielectric plates, separated by a gap distance of 1mm.Dielectric beads are inserted in the gap to form a PB-DBD reactor.
As the driving voltage increases, a pure surface discharge gradually dominates, because the charged species can more easily escape to the beads and charge the bead surface due to the strong electric field at high driving voltage.This significant surface charging will enhance the tangential component of the electric field along the bead surface, yielding SIWs.The SIWs, in turn, yield a high concentration of reactive species on the bead surface, which will affect the chemical reactions (and energy efficiency) in plasma catalysis.Indeed, electron impact excitation and ionization mainly take place on the bead surfaces.
The propagation speed of SIWs becomes faster in N 2 gas and slower with increasing O 2 content, because it is more difficult for the discharge to be created, and it yields a slower discharge evolution in electronegative gas, due to the loss of electrons by attachment to O 2 molecules.This trend can also be understood from the significant difference in ionization rates between N 2 and O 2 gases.These different ionization rates will create different amounts of electrons and ions on the dielectric bead surface, which might have consequences for plasma catalysis.

Figure 2 .
Figure 2. N + 2 ion density (m −3 ) for the same conditions as in Figure 1.

Figure 3 .
Figure 3. O + 2 ion density (m −3 ) for the same conditions as in Figure 1.

Figure 4
Figure 3 clearly shows that the O +2 ions are mainly formed along the dielectric surface by the SIWs, and the maximum density is about 20 times lower than the maximum electron density for all driving voltages, because it is more difficult for the discharge to become sustained in O 2 gas[38].This can be justified by the ionization rate shown later in Figure5.Figure4presents the O − 2 ion density distribution.The O − 2 ions are generated by electron impact attachment.Thus, the space and time distributions of the O − 2 ion density are nearly the same as for the electron density, but the maximum density is about two times smaller than the maximum electron density due to the smaller fraction of O 2 in air.

Figure 4 .
Figure 4. O − 2 ion density (m −3 ) for the same conditions as in Figure 1.

Figure 5 .
Figure 5. (a1,b1) Electron impact excitation rate (m −3 s −1 ) and (a2,b2) electron impact ionization rate (m −3 s −1 ), at 0.8 ns for a driving voltage of −20 kV, in an air discharge, for the nitrogen (a1,a2) and oxygen (b1,b2) components.The same color scale is used in all panels to allow comparison.The maximum values are noted in the figure.

Figure 6 .
Figure 6.Electric field amplitude |E| (V/m) for the same conditions as in Figure 1.The maximum values are noted in each panel.They occur at the contact points between the beads, and are larger than the color scale, but the color scale is chosen in such a way to better illustrate the general behavior.

Figure 7 .
Figure 7. Electric field components in the y and x direction, (a1,b1) E y (V/m) and (a2,b2) E x (V/m), (a1,a2) for −20 kV (0.8 ns), and (b1,b2) for −30 kV (0.2 ns), in an air discharge.The maximum values are noted under the figure.They occur at the contact points between the beads, and are larger than the color scale, but the color scale is chosen in such a way to better illustrate the general behavior.

12 .
Geometry used in the packed-bed dielectric barrier discharge (PB-DBD) reactor.The entire simulation domain is 10 mm in the x dimension and 1.65 mm in the y dimension, but only a smaller part in the x direction is shown here, for clarity.The discharge gap is 1 mm, and the top and bottom electrodes are covered by 0.3 mm thick dielectric plates.The top electrode is 0.05 mm thick, and it acts as powered electrode.The bottom electrode (x = 0) is grounded.Dielectric packing beads with diameter of 508 µm are placed in the gap.The numbering is used later in the paper.

( 5 )
Here, p represents the electrons (e) and ions (N + 2 , O + 2 and O − 2 ), and r p , v p , Z p , m p and S are the particle position, velocity, charge, mass and shape function in space.The shape function is chosen to be a first order b-spline (b l ) function, S(r − r p ) = b l ( r−r p r ), with r = δxδy.The time step δt fulfills the Courant condition cδt < 1 -impact ionization e + O 2 → 2e + O + O 2 → e + O 2 e + N 2 → e + N 2