Electron heating mode transitions in radio-frequency driven micro atmospheric pressure plasma jets in He/O₂: a fluid dynamics approach

Radio-frequency (RF) micro atmospheric pressure plasma jets (µ-APPJs) have become an attractive plasma source for surface treatment [1–6] and, in particular, for biomedical applications [7–11]. Such jets are usually operated in helium or argon with a small admixture of molecular gases, such as nitrogen, oxygen or combinations of both. The control of the production of reactive oxygen and nitrogen species (RONS) within the discharge volume is crucial for these applications. One method is the use of voltage waveform tailoring [12, 13], which has been shown the capability of tuning the dynamics of electron energy distributions (EEDs) [14–18]. Typically, the EED deviates from a Maxwellian distribution particularly for high energies even at atmospheric pressure since the energy relaxation length can be larger than length scale of the reduced electric field gradients. Furthermore, the energy relaxation frequency can be faster than the temporal change of the reduced electric field.

Initial investigations on electron power absorption dynamics in µ-APPJs were performed in nominally pure argon or helium. The electron heating modes were assumed to be similar to those of low pressure capacitive discharges, i.e. the α-mode and the γ-mode [19–22]. The α to γ transition was later discussed in the context of adding impurities [23, 24]. Hemke et al [25] firstly pointed out based on results of particle-in-cell/Monte Carlo collision (PIC/MCC) simulations for a pure helium discharge that the ionization dynamics is mainly produced by Ohmic heating in the so-called Ω-mode. It was demonstrated that even a small impurity in the noble gases can change the dominant ionization path to Penning ionization, leading to a decrease of the breakdown voltage [26, 27]. Experimentally, phase resolved optical emission spectroscopy (PROES) was used to study the dynamics of energetic electrons based on the wavelength integrated optical emission [28–31] and the use of selected emission lines (Ar: 750 nm [22], O: 844 nm [32, 33]). Bischoff et al [15] proposed that the 706.5 nm helium line can be used to probe the ionization dynamics in helium when using a small nitrogen admixture. Electron power absorption mode transitions were observed as well in that work by both experiments and PIC/MCC simulations for µ-APPJs operated in He/N₂ mixtures.

Low pressure electropositive radio frequency capacitively coupled plasmas (RF-CCPs) can operate in two modes, the α-mode and the γ-mode [34]. Similarly, RF driven µ-APPJs can also operate in two modes, the aforementioned Ω-mode and the Penning-mode. In the Ω-mode, electrons are accelerated by a bulk electric field while having a high neutral collision frequency leading to a decreased electron conductivity in spite of the electron conduction current being high. Although the spatio-temporal ionization dynamics of the Ω-mode is similar to those of the α-mode, the physical mechanism is different. In the α-mode in low pressure RF-CCPs, energetic electrons are mainly generated by sheath expansion. The electric field adjacent to the sheath can be either the ambipolar field or a drift field, or even a combination of both due to low electron density [35–37].

There is also a fundamental difference between the Penning-mode in µ-APPJs and the γ-mode in low pressure RF-CCPs. The largest electron impact ionization rate occurs inside the sheath in both modes at maximum sheath extension. In the γ-mode this ionization is produced by secondary electron emission from the electrodes. Whereas in the µ-APPJs, the maximum ionization in the Penning-mode is caused by energetic electrons producing highly excited states of the rare gas, followed by Penning ionization, for example, in He/O₂ mixtures, He* + O₂ $\rightarrow$ e + He + O$_{2}^{+}$.

Strongly electronegative capacitive RF discharges at low pressure, such as O₂ or CF₄, operate in a drift-ambipolar mode [38], dominated by a drift electric field in the bulk and an ambipolar electric field near the sheath edge. The drift field here results from a low electron density in the bulk due to depletion by attachment rather than the high collision frequency in the Ω-mode. The ambipolar field is a consequence of a local maximum in the electron density at the sheath edge.

Another important aspect of µ-APPJs for any kind of applications is the production of reactive species [39], which are highly relevant for chemical reactions. Complex plasma chemistries, such as He/O₂ mixtures, consist of hundreds of reactions. A sensitivity analysis [40], as well as a discussion of uncertainties and errors [41] induced by individual reactions were performed in the frame of a global model. Simulated species densities (atomic oxygen densities [33], ozone densities [42]) were benchmarked against experiments. The major reactions that lead to the generation and loss of these two species were discussed. Reactive species generations were also found to be affected by discharge parameters. The pressure dependence of O₂(a$^{1}\Delta_{g}$) densities [43], the effect of power pulse control on O₂(a$^{1}\Delta_{g}$) production [44], as well as the influence of the electrode configuration on RONS production [45] were investigated based on two-dimensional fluid simulations. Additional reactive neutral species are generated through interactions between the effluent and its surroundings. Parametric studies on RONS production in a APPJ flowing into humid air [46], or interacting with water [47] were performed based on a two-dimensional plasma hydrodynamics model.

In this work, we investigate a radio frequency capacitively coupled micro atmospheric pressure He/O₂ plasma jet based on a two-dimensional fluid dynamics approach. The purpose of this study is to address the effect of voltage amplitude and molecular reactive admixtures on the electron heating mode transition, the charged species dynamics and the neutral species densities. Based on the simulation results and comparisons to experimental results, it is shown that a fluid model can qualitatively describe the variation tendency induced by control parameters and capture the main physical mechanism in such micro atmospheric pressure discharges. The remaining content of this paper is structured as follows: a brief introduction of the simulation model is provided in section 2, followed by discussions of results in section 3. Finally, conclusions are drawn and prospects are made.

2.1. Description of the simulation model

The simulations in this investigation were performed using the computational modeling platform nonPDPSIM developed by Mark Kushner [48, 49]. It is a two-dimensional multi-species fluid dynamics code based on an unstructured grid, used for medium to high pressure weakly ionized plasmas. The model and prior applications have been previously discussed in detail [48, 49]. nonPDPSIM has been successfully applied to the study of µ-APPJs [21, 39, 43–47]. Here, only a brief description of nonPDPSIM is provided.

For each charged species, the particle conservation equation is solved. The particle flux is expressed in terms of the drift-diffusion approximation. In order to allow for the non-Maxwellian characteristics of the electron energy distribution function, the electron transport coefficients, as well as the rate constants in the source and loss terms, are obtained by solving a zero-dimensional Boltzmann equation based on the two-term approximation. These generated coefficients are firstly tabulated as a function of the reduced electric field, subsequently altered in dependence of the mean electron energy, or the effective mean electron temperature. The effective mean electron temperature is obtained by solving the electron energy conservation equation. The fluid equations of charged species are coupled to incompressible (or if needed compressible) Navier–Stokes equations, which are used to describe the neutral species transport. Poisson's equation is solved for the electric potential.

A schematic of our simulation domain is shown in figure 1. It is based on the COST reference micro-plasma jet introduced in [50]. In our simulations, a finite volume of the model is considered, whose edges are grounded. The powered electrode is at the bottom and the grounded electrode is at the top. Two dielectrics are placed between the respective electrode and the adjacent grounded wall, with a relative permittivity equal to 4. Gases flow in on the left side and are mixed in the gas mixing volume. It is 5 mm in height. Therefore, the discharge cannot be ignited in this region at voltages typically used to drive the jet. After the mixing region there is the discharge channel of 30 mm in length and 1 mm in height. A side chamber is located after the discharge domain for the effluent before gases flow out on the right. The unstructured mesh of design includes approximately 12 000 nodes for the plasma region.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The discharge is operated in He/O₂ mixtures at atmospheric pressure. The gas flow rate can significantly affect plasma discharges and species generations. A fast gas flow will carry most of the generated neutral species outside the jet, meanwhile, it will cause turbulence in the jet. Those together lead to unsteady discharges. While a slow gas flow will lead to reactive species fully built up in the jet so that the treatment distance is limited. As a result, the gas flow is influenced by the length of the jet, and it is required to generate high fluxes of reactive species at the nozzle in an energy efficient way. Considering the feature size of a COST jet, it usually works in a gas flow rate approximate 1 slm [15–17, 33, 39, 50, 51]. The gas flow is fixed at 1 slm in simulations. The oxygen concentration ratio is set to 0.05% , 0.25%, and 0.5% for different simulation cases. The voltage source is fixed at a frequency of 13.56 MHz and the voltage amplitude is varied from 400 to 600 V. The included species are ground state neutral species He, O₂, O, O₃, excited state neutral species O₂(v = 1–4) (first four vibrational levels of O₂), O₃(v), O₂(a$^{1}\Delta_{g}$), O₂(b$^{1}\Sigma_{g}^{+}$), O(¹D), He* (ensemble of He(2³S) and He(2¹S)), positive ions O$_{2}^{+}$, O⁺, He⁺, negative ions O⁻, O$_{2}^{-}$, O$_{3}^{-}$, and electrons. For electron impact collisions with helium atoms, cross sections of the elastic [52], the excitations [52, 53] and the ionization process [52] are considered. For electron impact collisions with oxygen molecules, the reactions and the corresponding cross sections proposed by Gudmundsson et al [54] are considered. These cross sections are used to generate transport coefficients by solving the Boltzmann equation. Reactions between electrons and other neutral species are neglected, since those neutral densities are at least one order of magnitude lower than the molecular oxygen density. Reactions between heavy species (ions and neutrals) are based on Turner [41]. The number of total reactions is reduced according to the sensitivity analyses [40]. The chemical reactions considered in this work are the same as those listed in [18]. Surface coefficients for neutral species are treated as a loss probability [33, 55]. The coefficients are identical to those listed in [18]. The ion induced secondary electron emission coefficients are chosen to be 0.2, 0.06, and 0.1 for He⁺, O$_{2}^{+}$, and O⁺ respectively based on the formula given in [56].

2.2. Description of the experimental set-up

The first row of figure 2 shows the computational spatio-temporally resolved He(3³S) excitation rates as a function of the driving voltage amplitude at 400, 500, and 600 V (different columns), with the gas flow fixed at 1 slm and the oxygen concentration kept constant at 0.05%. The results are taken at the center of the discharge channel (marked by the black dotted line in figure 1). The maximal excitation rates are different for each case. The results are normalized by the respective maximum. As the driving voltage amplitude increases, the spatio-temporal dynamics of the excitation rate changes, i.e. an electron heating mode transition is induced. At lower voltage, the majority of energetic electrons are generated in the plasma bulk when the sheaths are expanding; while at higher voltages, strong electron impact excitation rates appear inside the sheaths at the time of maximal sheath voltage. It should be noted that even though the electric field in the sheath is much stronger than in the bulk at 400 V, the electron density in the bulk is much higher compared to the sheath. As a consequence, most energetic electrons (above 22.7 eV) are generated by the acceleration due to the bulk drift field and the discharge is operated in the Ω-mode [15, 25].

When the voltage is increased to 600 V, the electric field in the sheath becomes strong enough to dominate the production of energetic electrons, causing the discharge to be operated in the Penning-mode [15, 25]. As can be seen in figure 2, the strong electron impact excitation rates occur when the sheaths expand to the maximum, i.e. when the electric field in the sheath is strongest. Those electrons mainly originate from Penning ionization inside the sheath. The same electron heating mode transitions are found via PROES measurements in a lower voltage amplitude range from 270 to 355 V, as shown in the second row of figure 2. Each case is normalized by the respective maximum. A larger grounded electrode compared to the powered electrode leads to the weak asymmetry of the patterns in the experimental results. The vertical stripes in experimental results are due to synchronization issues of the camera. As the camera is directly triggered by the applied waveform, at some moments the internal delay generator of the camera is not very stable, thus, leading to an inaccurate set of the delay time. The difference of the working voltage amplitude between simulations and experiments is due to the ignorance of the electron kinetic effects in fluid simulations, particularly for the electrons at high energy, such as the electrons generated from Penning ionization and surface emissions, which are then accelerated by the electric field. On the other hand, the total number of non-uniform unstructured cells is limited to a reasonable number to save the computational cost. Due to the far larger effective electrode length than the electrode gap, a relatively coarse meshing along the electrode gap is used, which can lead to a reduced precision of the spatial resolution between the electrodes. However, the simulation results still show a qualitative agreement with the experiments. It is very important to understand such a mode transition, since it results in a lot of plasma parameter variations. For example, as shown in the third row of figure 2, the time-averaged helium metastable density (simulated) is increased, with two peaks near the sheath edges at 600 V. Helium metastables are generated via electron impact excitations, while mainly destructed by Penning ionizations with O₂. Both reactions proceed fast. Since O₂ is uniformly distributed, the helium metastable density profile depends on the generation rate distribution, i.e. it is the highest at the center of the discharge gap in the Ω-mode, while two peaks near the sheath edges are formed in the Penning-mode.

The first row of figure 3 shows the computational spatio-temporally resolved He(3³S) excitation rates as a function of the O₂ concentration. The voltage amplitude is fixed at 500 V. The results are taken again at the center of the gas flow channel (marked by the black dotted line in figure 1). The maximal excitation rates are different for each case. The results are normalized by the respective maximum. For the 0.05% O₂ concentration case, the excitation rates are stronger in the sheaths when the sheaths expand to the maximum, indicating that the discharge is predominantly operated in the Penning-mode. Increasing the O₂ concentration induces the discharge to be predominantly operated in the Ω-mode, since the excitation rates are stronger in the plasma bulk when the sheaths are expanding. Such a transition is also shown by PROES measurements in the second row of figure 3 at 355 V. The slight asymmetry is caused by the larger area of the grounded electrode. The reason for the use of higher voltage amplitudes in the simulations has been discussed above. This transition is caused by the combination of two effects: electronegativity and collisions. More electronegative gas can lead to more negative ions and a lower electron density in the bulk (shown below), leading to a stronger bulk drift electric field, since the electron conductivity is inversely proportional to the electron density. This is similar to the drift pattern of the drift-ambipolar mode [38] in low pressure RF strongly electronegative capacitive discharges. Besides, the destruction rate of He* is enhanced by the increased O₂ density. However, the majority of electrons generated via the Penning ionization cannot be accelerated to high energy due to the more frequent inelastic collisions in the presence of more molecular gas, leading to a decreased population of He*. Correspondingly, the helium metastable density decreases significantly by adding more O₂ as shown in the third row of figure 3.

Zoom In Zoom Out Reset image size

As shown above, the helium metastable density is highly relevant to the electron heating dynamics, and so are the charged species densities. Figure 4 shows the simulated time-averaged charged species density profiles between the electrodes as a function of the voltage amplitude. The O₂ concentration is kept constant at 0.05%. All the results are taken at the center of the discharge channel (marked by the black dotted line in figure 1). In such cases, merely electrons, O$_{2}^{+}$ and O⁻ are dominant, while the other charged species are negligible. Their densities are increased by increasing the voltage amplitude, but the variation for O⁻ is weak. This is due to the competition between the enhanced major generation (electron impact dissociation attachment) and destruction (recombination with O$_{2}^{+}$ and reactions with oxygen neutrals). It can be seen that the time-averaged quasi-neutrality is broken, since the O$_{2}^{+}$ density is higher than the sum of electron and O⁻ densities, for example, at the center of the discharge gap. We believe that it is an averaging effect. The boundary loss for electrons is pronounced due to a small discharge gap (1 mm) and the strong oscillation driven by the RF electric field, which results in a narrow profile of the electron density at each moment. Ion densities are almost time independent in steady state, thus, the spatio-temporal profiles of ion densities are not shown. The instantaneous quasi-neutrality is fulfilled locally (shown in figure 5), while the mean electron density is decreased by averaging over one RF period.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 6 shows the simulated time-averaged charged species density profiles between the electrodes as a function of the O₂ concentration. The voltage amplitude is kept constant at 500 V. All the density profiles are shown at the center of the discharge channel (marked by the black dotted line in figure 1). At low O₂ concentration (0.05%), the electron density is higher than the negative ion densities. The major negative ion is O⁻, while O$_{2}^{-}$ and O$_{3}^{-}$ are negligible. As the O₂ concentration is increased, the electron density decreases, and the hump in the density profile becomes significant due to the enhancement of the negative ion population at the center of the electrode gap. At 0.5% O₂ concentration, the O⁻ and O$_{3}^{-}$ densities predominate over the electron density. The O⁻ density increases slightly compared to the significant enhancement of the O$_{2}^{-}$ and O$_{3}^{-}$ densities. To understand such a variation of the negative ion density as a function of the O₂ concentration, it is necessary to investigate the major chemical reactions leading to the construction and destruction of those negative ions. We estimate the production rate and loss rate induced by each relevant reaction based on the time-averaged density of each species and the corresponding rate constant. By comparing those estimated production and loss rates, we find that O⁻, O$_{2}^{-}$ and O$_{3}^{-}$ are mainly produced by the reactions

$$\begin{align} & \mathrm{e} + \mathrm{O}_{2} \rightarrow \mathrm{O} + \mathrm{O}^{-}, \end{align} \tag{ 1 }$$

$$\begin{align} & \mathrm{e} + \mathrm{O}_{2} + \mathrm{He} \rightarrow \mathrm{O}_{2}^{-} + \mathrm{He}, \end{align} \tag{ 2 }$$

$$\begin{align} & \mathrm{O}^{-} + \mathrm{O}_2 + \mathrm{He} \rightarrow \mathrm{O}_{3}^{-} + \mathrm{He}, \end{align} \tag{ 3 }$$

while the destruction rates, as mentioned above for negative ions, differ merely slightly between those different species. It can be seen that the generation of negative ions is proportional to the O₂ concentration. However, reaction (3) also corresponds to a destruction mechanism of O⁻, causing the O⁻ density to increase only slightly by increasing the O₂ concentration. Such variations of the negative ion and electron densities as a function of the voltage amplitude and the O₂ concentration can affect the electronegativity. Here the electronegativity is defined as the ratio between the sum of the time-averaged negative ion densities and the time-averaged electron density at the center of the discharge gap. As shown in figure 7, this ratio decreases as the voltage amplitude increases (red square). However, a significant increase of the electronegativity results from the increasing O₂ concentration (blue dot).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As a consequence of the variation of the electronegativity, the electric field perpendicular to the electrode changes with the O₂ concentration. In figure 8, the left figure shows the electric field component perpendicular to the electrodes at t₁ = 8.3 ns corresponding to the moment of the maximal excitation rate in the bulk (marked in the first row of figure 3). We restrict the region from 0.2 to 0.6 mm to magnify the differences among the bulk electric fields. The electric field in this region mainly causes the Ω-pattern of the electron heating dynamics. This drift electric field becomes larger as the O₂ concentration increases, since the electron density gets lower, inducing an attenuation of the electron conductivity. It is known that the strong bulk electric field in atmospheric pressure discharges is caused by the very high collisionality of electrons, mainly due to electron elastic collisions, since cross sections of electron elastic collisions are much larger than those of electron inelastic collisions. In our study cases, the predominant contribution is electron helium elastic collisions due to the high flow of the helium gas. The increasing O₂ concentration from 0.05% to 0.5% cannot lead to a significant increase of the total collision frequency because O₂ is still a very small admixture and the cross section of electron oxygen elastic collisions is comparable with the one of electron helium elastic collisions. Attentively, here we only point out that the enhancement of the bulk electric field is not mainly caused by the elastic collisions, but by the increased electronegativity induced by adding more O₂. The right plot of figure 8 shows the electric field perpendicular to the electrodes at t₂ = 19.2 ns corresponding to the moment of the maximal excitation rate in the top sheath (marked in the first row of figure 3). Again, only the region between 0.6 and 0.9 mm is shown to magnify the differences. The electric field in the top sheath becomes smaller as the O₂ concentration increases, which weakens the Penning-pattern of the electron heating dynamics there. Interestingly, the low electric field in the sheath leads to a self-amplification mechanism of this electron heating mode transition. Since the voltage is fixed across the discharge gap, the enhanced electric field in the bulk corresponds to a larger potential drop there. As a result, the voltage drop across the sheath is reduced, leading to the decreased electric field there.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

So far, we have discussed the electron heating mode transition dependence on the voltage amplitude and the O₂ concentration. We have also analyzed the charged species density variations induced by the control parameters and the effect of the electronegativity on the electron heating dynamics. Another important aspect for practical applications is the production of the desired reactive neutral species, such as atomic oxygen and ozone. We firstly focus on the computational distributions of neutral species densities along the direction of the gas flow at the fixed voltage amplitude of 500 V and the O₂ concentration of 0.05%, as shown in figure 9. The densities are taken at the center of the electrode gap (highlighted by the blue dotted line in figure 1). Regarding the coordinate system, the gas nozzle is located at 0, indicating that the coordinate is negative inside the jet (discharge domain), while the coordinate of the effluent is positive. The He*, O(¹D) and O₂(v = 1–4) densities are flat distributed inside the jet, and rapidly decrease to be negligible after the nozzle. The O, O₂(a$^{1}\Delta_{g}$) and O₂(b$^{1}\Sigma_{g}^{+}$) densities increase in the jet and reach their highest values in a short distance from the nozzle outside the jet, and decrease slightly in the effluent. The O₃ density is continuously increasing until the end of the effluent propagation. Similar results have been reported by Hemke et al [39]. Such neutral species distributions along the gas flow result from the major generation and destruction reaction rates, together with the effect of the gas flow. According to the sensitivity analyses based on chemical reactions in He/O₂ plasma jets in [33, 40], He*, O(¹D) and O₂(v = 1–4) are generated via electron impact excitation and quenched by He or O₂. Both processes are fast so that their density distributions are uniform in the jet. No source is provided outside the jet, so that they vanish rapidly. O, O₂(a$^{1}\Delta_{g}$) and O₂(b$^{1}\Sigma_{g}^{+}$) are also generated by electron impact molecular oxygen reactions at fast rates, but vanish slowly by reactions between neutrals. In such cases, the gas flow affects the transport of those species to show non-uniform distributions in the jet. Due to their relatively slow destruction rates, the species generated in the jet can be carried by the gas flow to the effluent. The decrease in the effluent is because of the lack of sources, but is not significant. O₃ is generated and destructed by neutrals slowly. The main generation reaction is

$$\begin{align} \mathrm{O} + \mathrm{O}_2 + \mathrm{He} \rightarrow \mathrm{O}_{3} + \mathrm{He}, \end{align} \tag{ 4 }$$

which provides a source in the effluent.

Figure 10 shows three reactive oxygen species (O₂(a$^{1}\Delta_{g}$), O and O₃) density distributions along the gas flow as a function of the voltage amplitude and the O₂ concentration. Densities are taken at the center of the electrode gap (highlighted by the blue dotted line in figure 1). Those species densities can be enhanced by either increasing the voltage amplitude or the O₂ concentration within the range discussed in this work. But increasing the O₂ admixtures is more effective.

Zoom In Zoom Out Reset image size

In this work, the spatio-temporal electron heating dynamics in a He/O₂ RF capacitive atmospheric pressure micro plasma jet is investigated by two-dimensional fluid dynamics simulations and PROES measurements. Two electron heating modes (the Ω-mode and the Penning-mode) are found and the electron heating dynamics can be changed by varying the voltage amplitude as well as the O₂ concentration. Based on the analyses of the construction and destruction mechanisms of the negative ions, it is pointed out that the generation of the negative ions is largely dependent on the O₂ concentration. The increased negative ion density induced by adding more O₂ enhances the drift bulk electric field, which contributes to the electron heating mode transition. The densities of different neutral species are found to show different distributions along the direction of the gas flow inside the jet and in the effluent due to the species relevant chemical reaction rates as well as the effect of the gas flow.

We notice that the electron heating mode transition can be induced within a lower range of the voltage amplitude in experiments (the second rows of figures 2 and 3). However, our simulated results clearly show that the fluid model is capable of investigating RF micro atmospheric pressure discharges qualitatively and providing guidance for practical applications. To achieve quantitative agreement of the results between the electrodes with experimental measurements, a kinetic treatment of electrons is needed. Complex chemistry between heavy species can hardly be treated by PIC/MCC simulations to the best of our knowledge. A real hybrid model would be a more effective solution for simulating sophisticated discharges with multiple species more accurately. This has been already shown by the authors elsewhere based on a coupling between a 1d3v (one dimension in displacement, three dimensions in velocity) PIC/MCC model and a 2D simplified fluid model [18]. However, in order to consider the lateral gas flow, a sophisticated hybrid model including 2d3v PIC/MCC algorithm and 2D fluid model would be ultimately needed. Moreover, a large number of cells are required to resolve the plasma dynamics precisely in both directions. The implementation and application of such model is a future goal. In the meantime, one has to rely either on a hybrid 1d3v PIC/MCC model that allows for a correct kinetic description of electrons but omits the actual gas transport, or on the fluid model which considers the transport of all species, but treats electrons not fully kinetically.

This work is supported by the German Research Foundation in the frame of SFB 1316 (Projects A4 and A5) and MU2332/11-1. We gratefully thank Prof Mark Kushner for his support and guidance of using nonPDPSIM. We also thank Dr Zoltán Donkó for fruitful discussions. Y. Liu and T. Mussenbrock acknowledge the hospitality of Brandenburg University of Technology Cottbus-Senftenburg where most of the simulation work was performed.

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