A non-neutral regime of radio-frequency atmospheric pressure plasma jets:
Simulation and modeling

Non-equilibrium plasmas have become a versatile and commonly used tool in modern society [1, 2, 3, 4, 5]. The origins of gas discharges trance back to discharges ignited under atmospheric pressure conditions. These first human-created plasmas root back to the 19th century [5, 6, 7]. Recently, those discharges were re-established as the focus of many works and considered to be the solution to fundamental challenges of the 21st century. The proposed applications range from environmental and energy-related topics such as CO₂ conversion [8, 9] to medical applications such as cancer treatment [10, 11] or wound healing [12, 13]. (A far more complete list of applications is given by Adamovich et al. [14]). The foundation of all applications of atmospheric pressure plasmas is their characteristic to enable complex chemistry without thermal reactions [1, 4, 6, 7, 14]. For the biological surfaces in medical applications [10, 11, 12, 13] and other heat-sensitive materials [1, 2, 3, 6], this characteristic of non-equilibrium plasmas allows for treatment at all.

While the above-mentioned reasoning is proper for both atmospheric pressure and low-pressure plasmas, the most significant difference between both types of discharges is the vacuum chamber. The opportunity to be installed in ambient air renders an atmospheric pressure plasma application, in most cases, cheaper and the technical requirements less sophisticated [6]. Additionally, applications that involve biological probes, bodily tissues, or even a whole patient [6, 10, 11, 12, 13, 14] cannot be realized under conditions anyhow nearby to vacuum. Nevertheless, atmospheric pressure non-equilibrium plasmas have specific difficulties too. Most notable in this context is the prevalence to transit, if provided with sufficient power, to a local thermodynamic equilibrium [3, 4, 6]. Numerous different concepts for plasma sources and even more realizations have been utilized to circumvent this difficulty [3, 6, 15].

Standardized plasma sources are called for to limit the vast parameter space. Therefore, the European COST (Cooperation in Science and Technology) Action MP1011 on ’Biomedical Applications of Atmospheric Pressure Plasma Technology’ provided Golda et al. [15] with the opportunity to build the COST Reference Microplasma Jet (COST-Jet) as a reference source [16]. In the following, the COST-Jet has proven to be highly reproducible [17], suitable for both experimental and theoretical studies [18], and surpassed the status of ”just a reference” source [19]. Hence, we choose the geometry of the COST-Jet for our models and simulations.

Despite many advances and developments in recent years [6, 14], some basic principles, reaction pathways, in short, a complete and fundamental understanding of atmospheric pressure non-equilibrium plasmas remains desirable. The intrinsic complexity of the non-linear system of non-equilibrium plasmas is coupled to transient physics and inherently complex chemistry at ambient pressure. As a result of this contraction, theory and simulations struggle to keep up with the challenges. To completely describe the tightly bounded capacitively coupled microplasmas present at atmospheric pressures, a multi-time-scale, kinetic, and three-dimensional model that considers all chemical details is needed. Nevertheless, all these at once are currently simply not feasible. A complete description of the chemistry involves dozens of species and hundreds of reactions [3, 4, 20, 21, 22, 23]. Models trying to resolve as many chemical details as possible are either global models with no spatial resolution [20, 21] or at maximum one-dimensional fluid models [22, 23]. Models prioritizing the spatial resolution result in two-dimensional fluid simulations [24, 25, 26] or one-dimensional kinetic models [27, 28] and restrict their chemistry to the essential minimum.

As the topical review by Bruggeman et al. [6] stresses, it is common knowledge that atmospheric pressure discharges do not necessarily form a quasi-neutral plasma bulk. This effect stems from the scaling of the Debye length $\lambda_{\mathrm{D}}$ compared to the discharge length $L$. Boundary effects (e.g., sheath formation) take place on the length scale of $\lambda_{\mathrm{D}}\approx 10^{-1}-10^{-2}\,$mm and disturb the quasi-neutrality the bulk plasma strives for. For low-pressure plasmas, $L$ is in the range of several centimeters and boundary effects govern just a small part of the plasma. However, atmospheric pressure plasmas use small discharge lengths $L$ to avoid thermalization. Boundary effects may consequently affect the whole discharge and prevent the formation of a quasi-neutral bulk region. We chose to call a discharge regime where no quasi-neutral plasma bulk is formed a non-neutral discharge regime.

Although we suspect to see it within a He/O₂ mixture in previous work [25], a non-neutral discharge regime of the COST-Jet has not been analyzed in details. Employing a one-dimensional hybrid particle-in-cell/Monte Carlo collisions (PIC/MCC) simulation, we explore a non-neutral discharge for the He/N₂ mixture. Based on the simulation data, this work aims to describe a non-neutral discharge regime for the COST-Jet. An emphasis is made on the electron dynamics of this regime. We find the electrons organized within a moving Gaussian-shaped spatial profile, which features an immense form-stability.

We, according to its tendency to keep and quickly regain its shape, characterize this structure as drift-soliton-like. The expression soliton describes wave phenomena with a high degree of dimensional stability that occur in many fields of physics [29, 30]. In the context of plasma applications, actual drift-solitons have mainly been reported for magnetized plasmas [31, 32]. We find the plasma properties such as the electron temperature $T_{\mathrm{e}}$ dominated by the dynamics of the soliton-like structure. However, the hybrid PIC/MCC simulation functions in this regard more like a numerical experiment and is too sophisticated to develop a fundamental of the formed soliton. Therefore, an abstract and simplified analytical model for the electron dynamics is formulated. A detailed analysis of this model helps to explain the formation and spatial mold of the electron group and provides insights to its dominant influence on the electron dynamics. Additionally, deliberations on the stability of the analytical solutions that include the ideas of Lyapunov [33, 34, 35] will be used to argue for the soliton-analogy.

The rest of the manuscript is organized as follows. Section 2 gives a brief overview of the COST-Jet and describes the chemistry under investigation. The following section introduces our simulation framework and essential details of its realization. The simulation-data-based analysis of the electron dynamics follows in section 4. Section 5 will provide the aforementioned analytical model that provides generalization and further insight. The analysis of the model is acompanied by deliberations on the stability of the soliton solution. A comparison between the simulation and our analytical model is made in section 6 and aims to visualize the reliability of the analytical model. The final section summarizes our findings, draws conclusions, and outlines perspectives for future work.

The COST-Jet is a capacitively coupled radio-frequency-driven plasma jet specifically designed to function as a reference plasma source. Nevertheless, recent work [19] shows that the COST-Jet is a more than decent plasma source for many applications. While Golda et al. [15] provide a complete characterization and definition of the COST-Jet, this section gives a brief introduction to the device’s essential features. By design, the COST-Jet is optimized for diagnostic purposes, especially optical access and is suitable for theoretical description. This choice is reflected in the two quartz plates glued to two stainless steel electrodes to form the rectangular discharge channel. The quartz plates allow for optical diagnostics along the horizontal gas-flow axis $x$ and vertical axis $z$. Figure 1 a) shows a sketch of the cross-section of the COST-Jet along the gas-flow-axis $x$. In addition, the capacitively coupled power delivery at 13.56 MHz combined with the rectangular geometry allows one-dimensional simulations to give meaningful results [24, 22, 23, 36, 37].

Although we firmly believe that a complete model of the Cost-Jet, a long tube ($30\,\mathrm{mm}$) with a tightly bounded $1\,\mathrm{mm}\times 1\,\mathrm{mm}$ lumen, has to regard all three spatial dimensions somehow, technical limitations currently force us to focus our study on a one-dimensional model. Similar to previous studies [24, 22, 23, 36, 37], we focus on describing the plasma between the electrodes along the $z$-axis and perpendicular to the gas flow. Figure 1 b) shows a cross-section of the setup, and the green line marks the position of the $z$-axis. For now, the gas flow is neglected and the plasma is assumed to be invariant in both the $x$- and $y$- direction.

In terms of chemistry, the COST-Jet is operated using a chemically inert carrier gas (usually helium [15, 24, 22, 23, 25, 26, 36, 37], less commonly argon [18, 44]) mixed with a molecular trace gas (usually nitrogen [24, 36, 37, 44] and/or oxygen [15, 22, 23, 25, 26, 18, 19, 44, 45]). For this study, we decided to examine the helium/nitrogen mixture. The resulting electropositive plasma yields a moderate amount of species and chemical complexity. Both are beneficial for the analytical model that will be presented in a following section. The details of the chemistry set used for this study are mainly oriented on previous work [28] and are shown in table 1.

The particle-in-cell (PIC) simulation is a numerically sophisticated tool that provides self-consistently calculated statistical representations of the distribution functions [46, 47, 48, 49]. The foundation of the PIC algorithm traces back to the 1940s [46]. In the 1960s, a Monte Carlo collisions (MCC) scheme was added that lead to the commonly used particle-in-cell/Monte Carlo collisions (PIC/MCC) simulation [47]. In recent years, the PIC/MCC simulation has had great success in the research of non-equilibrium plasmas in general [1, 50, 48, 49, 51], and the scheme has been applied in the atmospheric pressure regime [36, 37, 28]. Nevertheless, the above-discussed demands of atmospheric pressure discharges severely drive PIC/MCC simulations to their limits. Therefore, hybrid PIC/MCC simulation schemes for RF-driven plasma jets are utilized to lower the computational burden while keeping all essential pieces of physics [27, 52, 50].

Overall, hybrid concepts in plasma modeling are commonly used. The review paper of van Dijk et al. [50] lists numerous techniques, variations, and applications. The hybrid PIC simulation has many applications in astrophysical plasma studies [53, 54, 55, 56]. These studies reduce the numerical load of their calculations by describing the electrons as a fluid and solely resolving the ion dynamics kinetically. As Eremin et al. [27] point out, the situation in RF plasma jets is precisely inverted. Electrons bear dynamics that require a fully kinetic description, while a fluid description of the ion dynamics leads to a significantly more efficient model.

For this study, we adapt the hybrid PIC/MCC simulation concept and develop the one-dimensional explicit electrostatic Hybrid particle-in-cell/Monte Carlo collisions code eehric. The following subsections present some details of the implementation.

The most notable characteristics of the non-neutral regime of atmospheric pressure capacitively coupled plasma jets are seen in the density profiles. Figure 2 a) shows a snapshot of the electron density $n_{\mathrm{e}}$ (blue), the total ion density $n_{\mathrm{i,tot}}$ (red) with spatial resolution. The panels b) ($n_{\mathrm{e}}$), c) ($n_{\mathrm{N\textsubscript{2}\textsuperscript{+}}}$)), and d) ($n_{\mathrm{i,tot}}$) provide spatially and temporally resolved representations of the respective number densities. First, the figure shows that the N₂⁺ ions are the clearly dominant ion species for this case (fig. 2 c), and d)). Panels c) and d) show no by naked eye visible difference. The density of He₂⁺ ions already is about three orders of magnitude lower than that of the N₂⁺ ions. For this case, the density of He⁺ ions, which is even lower than the density of the He₂⁺ ions. Hence, both ion densities are not shown. The extremely low He⁺ ion density stems from the efficient ion conversion He⁺ $\to$ He₂⁺ (c.f. reaction 32 tab. 1). This result is in accordance to Martens et al. [70]. Both the He₂⁺ and the N₂⁺ ion density profiles exhibit only a slight temporal modulation. The plasma is operated in the RF regime ($\omega_{\mathrm{p,i}}\ll\omega_{\mathrm{RF}}\ll\omega_{\mathrm{p,e}}$).

For electrons, on the other hand, figures 2 a) and b) show a highly time-dependent structure. The spatial profile of the electron density $n_{\mathrm{e}}$ is approximately Gaussian-shaped (fig. 2 a)). The temporally and spatially resolved profile shows a narrow band of electrons pushed back and forth between the electrodes while approximately maintaining their mold (fig. 2 b)). This behavior reminds us on two of the three key characteristics of solitons as given by Drazin and Johnson ([29] p. 15): (i) a soliton has a permanent shape, (ii) a soliton is localized within a specific region, and (iii) solitons can interact without changing their shape. The dynamics of the electron group of figure 2 clearly meets the first two criteria. The third criterion cannot be met by a group of electrons. If, for some reason, two of these electron groups existed and met, they would interact and eventually form one bigger group of electrons. Hence, the group of electrons is not a real drift-soliton. However, we decide to refer to this structure as a drift-soliton in a figurative way.

The elementary characteristic of the non-neutral operation regime is visible by simultaneously looking at the electron density and ion densities. The peak electron density locally never resembles the total ion density. Quasi-neutrality is violated along the whole discharge region. This non-neutral behavior of the discharge has two implications. First, the electron dynamics feature specific characteristics that will be discussed in the following. Second, the classical segregation of capacitively coupled plasmas in a sheath and a bulk region is rendered meaningless for this operation regime. The concept of the boundary sheath edge, especially, and the analysis of its dynamics that in other studies provide orientation and valuable insight [36, 37] cannot be used in the context of a non-neutral regime.

The black lines in figure 2 c) and d) are introduced to compensate this loss of reference. They refer to the rough position of the soliton-like structure in space and time. The shape of the electron group is assumed to be Gaussian to calculate the black lines. Based on this assumption, the position of the symmetry axis $Z$ is calculated as the first spatial moment of the electron density, and the width of the structure $\delta_{z}$ calculates from the second spatial moment. The black lines mark the interval $Z\pm 1.28\,\delta_{z}$, in which 90 percent of the electrons are located.

The non-neutral behavior of the discharge reflects in the electric field and potential. Figure 3 shows the electric potential $\phi$ (a)) and the electric field in $z$-direction $E_{z}$ (b)) with spatial and temporal resolution. The remaining panels show temporal snapshots of the individual profiles with spatial resolution ($\phi$: fig. 3 c), $E_{z}$: fig. 3 d)). Due to not having a quasi-neutral region, the potential of the discharge is, in first approximation, the electrical potential of a parallel plate capacitor bearing a positive space charge [71]. Therefore, the spatial profile of the electric potential resembles a compressed parabola (fig. 3 c)), and the maxima and minima of the potential are primarily located at the electrodes. Accordingly, the electric field $E_{z}$ does not show the bulk-specific plateau around zero. When the electron group is around in the middle of the discharge (e. g. fig. 3 d), blue and red lines), the spatial profile of the electric field becomes almost linear (with a small plateau at the position of the maximal electron density). For the brief moments when the soliton-like structure is closest to either electrode, a tiny area develops that is closest to being quasi-neutral. A corresponding plateau close to the respective electrode becomes visible (e.g., fig. 3 d), orange and purple lines). Additionally, and in contrast to a quasi-neutral regime, the electric field often has a single direction and does not cross zero (fig. 3 b) and d)).

The, in the context of a capacitively coupled plasma, unique electric field structure and the strict localization of the electrons within a soliton-like group reflect in the whole electron dynamics. Figure 4 shows the power density $P_{s}$ for different particle species (electrons: a), He₂⁺ ions: c), N₂⁺ ions: d)) and the electron temperature $T_{\mathrm{e}}$ (b)). The species-specific power density $P_{s}$ is calculated as the product of the individual current density $j_{s}$ and the electrical field $E_{z}$. We retrieve $T_{\mathrm{e}}$ from the diagonal element $p_{\mathrm{zz}}$ of the pressure tensor, which is defined as

$\displaystyle p_{\mathrm{zz}}=m_{\mathrm{e}}\,n_{\mathrm{e}}\,\left(\langle v_{z}^{2}\rangle-u^{2}\right)$

(7)

with $m_{\mathrm{e}}$ the electron mass, $n_{\mathrm{e}}$ the electron number density, $\langle v_{z}^{2}\rangle$ the electron’s mean-square velocity, and $u$ the mean velocity of the electrons. Details of these diagnostics are found in previous work [49].

The dynamics of the soliton-like structure are reflected in the electron power absorption patterns. As discussed before, 90 percent of the electrons are located inside the “soliton area” enclosed by the black curves (fig. 4 a)). Accordingly, energy loss and gain almost exclusively happen inside this region. Whenever the structure moves, electrons are accelerated along its path and gain energy from interacting with the electrical field. The most substantial energy gain is observed when the electron group approaches the middle of the discharge gap, and the movement speed of the group directly correlates to the energy gain. While changing its direction at the electrodes, the soliton-like structure becomes decelerated and accelerated in the opposing direction afterward. The described behavior is, for example, seen between $25$ and $45\,$ns (fig. 4 a)). Furthermore, the moments when the movement speed of the electron group drops to zero ($t\approx 35\,$ns and $t\approx 70\,$ns) are the sole opportunities for a significant amount of electrons to leave the discharge. At these times, the soliton-like structure is closest to the electrodes, and the loss of electrons coincides with the only visible net loss of energy.

In summary, the electron dynamics of the non-neutral regime are coupled to the soliton-like structure rather than the oscillating boundary sheath (that in this case -at least in its classical sense- does not even exist). As a result, energy gain patterns along the movement of the electron group outweigh energy loss patterns at a standstill by a lot. Additionally, the energy absorbed by electrons (fig. 4 a)) completely overshadows the energy delivered to the ions (fig. 4 c) and d)). The power absorption by ions (fig. 4 c) and d)) basically mimics the particular structure of the electrical field (fig. 3 b)). For the positive half-wave of the RF voltage ($t\approx 0-37\,$ns), the electrical field is approximatively positive over the whole discharge gap, and ions are pushed towards the grounded electrode. For the negative half-wave of the RF voltage ($t\approx 37-74\,$ns), the situation is reversed. Ions are pulled towards the powered electrode. The energy absorption linked to these dynamics is simple. Both He₂⁺ and N₂⁺ ions show a similar absorption pattern. The stronger the electric field, the more energy is absorbed. Accordingly, the maxima of the absorbed power are in front of the electrodes and coincide with the maxima of the electrical field (comp. fig. 3 b) to fig. 4 c) and d)).

Another impact of the field structure shows inside the temporal and spatial behavior of the electron temperature displayed in figure 4 b). The figure, first of all shows, that the majority of the electrons (the ones inside the “soliton area”) share the same temperature. Outside of the “soliton area”, structures are visible that correlate to the secondary electron emission. The maxima of the electron temperature are located in front of the electrodes and switch sides each half-wave. These maxima coincide with the maxima of the electrical field (fig. 3 b)) and the maxima of the ion energy absorption (fig. 4 c) and d)). At $t\approx 20\,$ns (grounded) and $t\approx 55\,$ns (powered), the amount of ions reaching the respective electrode is highest, and so is the production of secondary electrons. Simultaneously, the absolute value of the field strength is maximal, and the secondary electrons gain an equivalent amount of energy.

Nevertheless, figure 4 b) is misleading regarding the role of secondary electrons. First, their number is low. Recalling the meaning of the black curves, it becomes evident that a tiny fraction of a few percent of the electron population consists of secondary electrons. Second, secondary electrons absorb an insignificant amount of the energy. Figure 4 a) shows structures (in light red color) in front of the respective electrodes (grounded: positive half-wave, powered: negative half-wave). However, the simultaneous energy gain by electrons inside the soliton-like structure outweighs these patterns. Third, another kind of electrons, the ones created by Penning ionization, is more crucial to the non-neutral regime. An investigation of the ionization dynamics will prove this point.

First, the term “Penning electron” must be defined. For this study, we chose to define an electron as a “Penning electron” when created by the Penning ionization process (tab. 1: No. 31). Penning electrons keep their status until their first ionizing collision. Please note that this definition is somehow arbitrary. According to the above definition, an electron avalanche or ionization cascade started by a “Penning electron” will be missed. Therefore, the importance of Penning ionization may be underestimated by our means. As the following section will point out, this source of arbitrariness cannot affect our results, and any representation of different ionization channels is sufficient for our arguments.

Panels a) to c) in figure 5 present temporally and spatially resolved profiles of the ionization rate of nitrogen differentiated by the diagnostics as mentioned above. The ionization patterns of the electrons inside the soliton-like structure (fig. 5 a)) coincide with the patterns of electron power absorption (fig. 4 a)). As discussed before, electrons mainly gain energy while the electron group as a whole is moving. Accordingly, there is ionization during the same time. When the electron group stops (at powered electrode: $t\approx 35\,$ns, at grounded electrode: $t\approx 70\,$ns), electrons get lost, lose energy, and do not ionize. The light green structures adjacent to the respective electrodes (grounded: $t\approx 10-25\,$ns, powered: $t\approx 50-65\,$ns) stem from ionization associated with the secondary electrons. Their contribution is insignificant compared to the ionization inside the soliton-like structure and the ionization pattern of Penning electrons.

The latter group of electrons has a diffuse temporal structure (fig. 5 b)). Due to its nature and the high pressure, Penning ionization introduces a delay between the initial electron impact producing the metastable helium state and the actual ionization when the metastable decays via collision with a nitrogen molecule. As a result, the ionization pattern of Penning electrons is vaguely connected to the energy absorption (the more energy, the more likely an inelastic collision). However, the distinct temporal structure becomes smeared out. The remainder is a spatial profile with an ionization maximum between $z=0.4-0.6\,$mm. Nonetheless, Penning ionization happens always and anywhere among the discharge gap (except a small region in front of the electrodes where no high energetic electrons are present at all).

In the total ionization pattern (fig. 5 c)), the structure caused by the movement of the electron group is dominant. Nevertheless, the ionization band in the region $z=0.4$ to $0.6\,$mm caused by Penning ionization is still visible. This observation and the comparable scales of figures 5 a) and b) lead to the conclusion that Penning ionization has a significant share on the total ionization. Anyhow, figure 5 d) shows that the density made up by Penning electrons is insignificant, and there is no sign of different behavior of the Penning electrons. Apparently, electrons are rapidly incorporated into the soliton-like structure regardless of their origin. The discussion of figure 4 b) regarding the electron temperature $T_{\mathrm{e}}$ painted a similar picture. All electrons more or less share the same fate, and neither diagnostic nor analysis of the simulation data could provide an answer. For now, two key questions remain open: (i) How does the soliton-like structure form inside the non-neutral regime?, and (ii) Why do the electrons behave similarly?

To develop a better physical understanding of the discharge dynamics and to answer the questions raised above, we now analyze an elementary drift-diffusion model. Such models trade, in a way, quantitative accuracy for mathematical transparency [50]. We are interested in the behavior of a soliton (the term used in the sense of section 4) that is not disturbed by the boundaries, and extend the solution domain to the real axis. Moreover, we focus on the dynamics of the electrons and the electric field and treat the ion density $n_{\mathrm{i}}$ as a spatial and temporal constant. The electron model is given by the equation of continuity with the generation term neglected,

$\displaystyle\frac{\partial n_{\mathrm{e}}}{\partial t}+\frac{\partial\Gamma_{\mathrm{e}}}{\partial z}=0.$

(8)

The flux $\Gamma_{\mathrm{e}}$ is calculated with the diffusion constant $D_{\mathrm{e}}$ and the mobility $\mu_{\mathrm{e}}$,

$\displaystyle\Gamma_{\mathrm{e}}$

$\displaystyle=-D_{\mathrm{e}}\,\frac{\partial n_{\mathrm{e}}}{\partial z}-\mu_{\mathrm{e}}\,E_{z}\,n_{\mathrm{e}}.$

(9)

The electric field $E_{z}(z,t)$ and its potential $\phi(z,t)$ are given by Poisson’s equation:

$\displaystyle\varepsilon_{0}\,\frac{\partial E_{z}}{\partial z}=-\varepsilon_{0}\,\frac{\partial^{2}\phi}{\partial z^{2}}=e\,(n_{\mathrm{i}}-n_{\mathrm{e}}).$

(10)

The total discharge current density $j_{z}(t)$ is assumed to be periodic in time, average-free and spatially homogeneous. It is related to the surface charge density, $Q(t)$, as [72]:

$\displaystyle-e\,\Gamma_{\mathrm{e}}+\varepsilon_{0}\,\frac{\partial E_{z}}{\partial t}=j_{z}(t)=-\frac{dQ}{dt}.$

(11)

The system of equations is completed by suitable asymptotic conditions for $z\to\pm\infty$. The electron density is assumed to vanish and the electric field is asymptotically linear; the constants $E_{\pm\infty}$ have the dimension of the electric field:

	$\displaystyle n_{\mathrm{e}}\big{|}_{z\to\pm\infty}=0,$		(12)
	$\displaystyle E_{z}\,-\frac{1}{\varepsilon_{0}}\,e\,n_{\mathrm{i}}\,z\Big{|}_{z\to\pm\infty}=-\frac{1}{\varepsilon_{0}}\,Q(t)+E_{\pm\infty}.$		(13)

The global dynamics of the electron soliton can be described by a closed model. We define three spatial moments of the electron density, namely the total number $\mathcal{N}(t)$, the center of mass (COM) $\mathcal{Z}(t)$, and the soliton-averaged field $\mathcal{E}(t)$:

	$\displaystyle\mathcal{N}(t)$	$\displaystyle=\int_{-\infty}^{\infty}n_{\mathrm{e}}(z,t)\,\mathrm{d}z,$		(14)
	$\displaystyle\mathcal{Z}(t)$	$\displaystyle=\frac{1}{\mathcal{N}(t)}\int_{-\infty}^{\infty}z\,n_{\mathrm{e}}(z,t)\,\mathrm{d}z,$		(15)
	$\displaystyle\mathcal{E}(t)$	$\displaystyle=\frac{1}{\mathcal{N}(t)}\int_{-\infty}^{\infty}E_{z}(z,t)\,n_{\mathrm{e}}(z,t)\,\mathrm{d}z.$		(16)

Integration of the particle balance equation (8) over the spatial axis leads to the insight that the particle number $\mathcal{N}$ is constant. Applying the same integral to Poisson’s equation (10) results in an identity which states that the offsets in the asymptotically linear electric field forms are related to the electron number:

$\displaystyle\mathcal{N}=-\frac{\varepsilon_{0}}{e}\,\left(E_{+\infty}-E_{-\infty}\right).$

(17)

Integrating the continuity equation with the weight $z$, we derive a differential equation that describes the drift motion of the COM with regard to the average field:

$\displaystyle\frac{\mathrm{d}{\mathcal{Z}}}{\mathrm{d}{t}}=-\mu_{\mathrm{e}}\,\mathcal{E}.$

(18)

Further, by manipulating the model equations and utilizing the asymptotic conditions, we obtain two additional identities. The first identity states

$\displaystyle\mathcal{E}=\frac{e\,n_{\mathrm{i}}}{\varepsilon_{0}}\,\mathcal{Z}-\frac{1}{\varepsilon_{0}}\,Q(t)+\frac{E_{+\infty}+E_{-\infty}}{2}.$

(19)

The second identity directly relates the current and the density averaged field:

$\displaystyle\varepsilon_{0}\frac{\mathrm{d}\mathcal{E}}{\mathrm{d}t}+\mu_{\mathrm{e}}en_{\mathrm{i}}\mathcal{E}=j_{z}$

(20)

Equations (18) and (19) combine to a closed differential equation for the COM:

$\displaystyle\frac{\mathrm{d}\mathcal{Z}}{\mathrm{d}t}=\frac{\mu_{\mathrm{e}}}{\varepsilon_{0}}\,Q(t)-\frac{e\,n_{\mathrm{i}}\,\mu_{\mathrm{e}}}{\varepsilon_{0}}\,\mathcal{Z}-\mu_{\mathrm{e}}\,\frac{E_{+\infty}+E_{-\infty}}{2}.$

(21)

When viewed as an initial value problem, the solution $\mathcal{Z}(t)$ of this differential equation depends on the initial value $\mathcal{Z}_{0}$ at $t_{0}$. Asymptotically, however, a periodic solution $\mathcal{Z}_{\mathrm{as}}(t)$ is reached which solely depends on the modulation $Q(t)$ and the asymptotic conditions. The first term is constant and the second term is periodic and average-free:

$\displaystyle\mathcal{Z}_{\mathrm{as}}(t)=-\frac{\varepsilon_{0}}{2\,e\,n_{\mathrm{i}}}\,\left(E_{+\infty}+E_{-\infty}\right)+\frac{\mu_{\mathrm{e}}}{\varepsilon_{0}}\int_{-\infty}^{\mathrm{t}}\exp\left(-\frac{\mathrm{e}n_{\mathrm{i}}\mu_{\mathrm{e}}}{\varepsilon_{0}}(t-\tau)\right)Q(\tau)\,\mathrm{d}\tau.$

(22)

The corresponding density averaged field is

$\displaystyle\mathcal{E}_{\mathrm{as}}(t)=\frac{e\,n_{\mathrm{i}}}{\varepsilon_{0}}\,\mathcal{Z}_{\mathrm{as}}-\frac{1}{\varepsilon_{0}}\,Q(t)+\frac{E_{\infty}+E_{-\infty}}{2}.$

(23)

This point $\mathcal{Z}_{\mathrm{as}}(t)$ is now taken as the origin of new coordinates which we, for a reason that will become clear shortly, call the current-free system:

$\displaystyle z=\mathcal{Z}_{\mathrm{as}}(t)+\xi.$

(24)

The electron density, the electron flux, and the electric field are transformed as follows, where $n_{\mathrm{e}}^{\prime}(\xi,t)$, $\Gamma_{\mathrm{e}}^{\prime}(\xi,t)$, and $E_{z}^{\prime}(\xi,t)$ are measured in the new coordinates:

	$\displaystyle n_{\mathrm{e}}(z,t)$	$\displaystyle=n_{\mathrm{e}}^{\prime}(\xi,t),$		(25)
	$\displaystyle\Gamma_{\mathrm{e}}(z,t)$	$\displaystyle=n_{\mathrm{e}}^{\prime}(\xi,t)\,\frac{\mathrm{d}\mathcal{Z}_{\mathrm{as}}}{\mathrm{d}t}+\Gamma_{\mathrm{e}}^{\prime}(\xi,t),$		(26)
	$\displaystyle E_{z}(z,t)$	$\displaystyle=\mathcal{E}_{\mathrm{as}}(t)+E_{z}^{\prime}(\xi,t).$		(27)

Immediately dropping the prime, the original equations (8) to (10) remain unchanged. For further streamlining, we adopt a dimensionless notation. We introduce the Einstein temperature $T_{\mathrm{e}}=eD_{\mathrm{e}}/\mu_{\mathrm{e}}$, the Debye length $\lambda_{\mathrm{D}}=\sqrt{\varepsilon_{0}T_{\mathrm{e}}/e^{2}n_{\mathrm{i}}}$, and the electric relaxation time $\tau_{\mathrm{e}}=\varepsilon_{0}/e\mu_{\mathrm{e}}n_{\mathrm{i}}$. Measuring length in $\lambda_{\mathrm{D}}$, time in $\tau_{\mathrm{e}}$, density in $n_{\mathrm{i}}$, flux in $n_{\mathrm{i}}\lambda_{\mathrm{D}}/\tau_{\mathrm{e}}$, electric field in $T_{\mathrm{e}}/e\lambda_{\mathrm{D}}$, and electric potential in $T_{\mathrm{e}}/\lambda_{\mathrm{D}}$, we obtain:

	$\displaystyle\frac{\partial n_{\mathrm{e}}}{\partial{t}}+\frac{\partial\Gamma_{\mathrm{e}}}{\partial{\xi}}=0,$		(28)
	$\displaystyle\Gamma_{\mathrm{e}}=-E_{z}\,n_{\mathrm{e}}-\frac{\partial{n_{\mathrm{e}}}}{\partial{\xi}},$		(29)
	$\displaystyle\frac{\partial{E_{z}}}{\partial{\xi}}=-\frac{\partial^{2}{\phi}}{\partial{\xi}^{2}}=1-n_{\mathrm{e}},$		(30)
	$\displaystyle 0=-\Gamma_{\mathrm{e}}+\frac{\partial{E_{z}}}{\partial{t}}.$		(31)

The asymptotic conditions translate to

	$\displaystyle n_{\mathrm{e}}\big{|}_{\xi\to\pm\infty}=0,$		(32)
	$\displaystyle E_{z}\,-\xi\big{|}_{\xi\to\pm\infty}=\pm(E_{\infty}-E_{-\infty})/2.$		(33)

In the new coordinates, the equation of motion for the COM is a simple relaxation. The COM converges exponentially to the reference point 0. The same statement holds for the particle-averaged electric field:

	$\displaystyle\frac{\mathrm{d}{\mathcal{Z}}}{\mathrm{d}{t}}$	$\displaystyle=-\mathcal{Z},$		(34)
	$\displaystyle\frac{\mathrm{d}{\mathcal{E}}}{\mathrm{d}{t}}$	$\displaystyle=-\mathcal{E}.$		(35)

Thus, it is reasonable to assume that the structure as a whole relaxes to an equilibrium. The stationary solution of eqs. (28) - (31) has zero flux, obeys Boltzmann equilibrium, and is described by the a second-order differential equation for the potential $\phi$ known as the Poisson-Boltzmann equation [74, 75, 76, 77, 73]:

$\displaystyle-\frac{\partial^{2}{\phi}}{\partial{\xi^{2}}}=1-\exp{\left(\phi\right)}.$

(36)

The Poisson-Boltzmann equation is of second order and therefore needs two conditions. The symmetry around the origin demands $\partial{\phi}/\partial{\xi}=0$. The second degree of freedom can conveniently by fixed by setting $\phi(0)=\ln(n_{0})$, with $n_{0}\in\left[0,1\right[$. This parameter can be identified as the ratio between the electron density and the ion density at $\xi=0$. Analytical solutions to (36) are well known but quite inconvenient to handle [76, 77]. Figure 6 presents numerical solutions for the electron density $n_{\mathrm{e}}$, the electric potential $\phi$, and the electric field $E_{z}$. The character of the solutions depends on the parameter $n_{0}$. For $n_{0}$ close to $1$, the bulk-sheath structure of low-pressure discharges is recovered. For values of $n_{0}$ that are considerably smaller, strongly non-neutral structures appear. Compared to figures 2 b) and 3, the orange and green curves already approximate the previously discussed simulation results (cf. sec. 4).

Intuitively, the stationary solutions are stable. An effective method to show this formally is Lyapunov’s direct method [33, 34, 35]. The starting point is a balance equation for the density of the Helmholtz free energy [78, 79, 80, 75, 77]:

$\displaystyle\frac{\partial}{\partial{t}}\left(\frac{1}{2}\,{E_{z}}^{2}+n_{\mathrm{e}}\,\log(n_{\mathrm{e}})-n_{\mathrm{e}}\right)+\frac{\partial}{\partial{\xi}}\left(\Gamma_{\mathrm{e}}\,\log(n_{\mathrm{e}})\right)=-\frac{{\Gamma_{\mathrm{e}}}^{2}}{n_{\mathrm{e}}}.$

(37)

The Helmholtz free energy itself reads as follows, where the temporally constant last term in the integrand serves to render the integral finite:

$\displaystyle\mathcal{F}=\int_{-\infty}^{\infty}\Biggl{(}\frac{1}{2}\,{E_{z}}^{2}+n_{\mathrm{e}}\,\log(n_{\mathrm{e}})-n_{\mathrm{e}}-\frac{1}{2}\left(\left|\xi\right|+\frac{\mathcal{N}}{2}\right)^{2}\Biggr{)}\,\mathrm{d}\xi$

(38)

The time derivative of $\mathcal{F}$ is negative semi-definite, so that it is a suitable candidate for a Lyapunov functional [80]. Note that the equality sign is only assumed when the flux and with it the dissipation vanishes:

$\displaystyle\frac{\mathrm{d}\mathcal{F}}{\mathrm{d}t}=-\int_{-\infty}^{\infty}\frac{{\Gamma_{\mathrm{e}}}^{2}}{n_{\mathrm{e}}}\,\mathrm{d}\xi\leq 0.$

(39)

A solution of equation (36) is stable when it represents a minimum of the free energy under the constraint of a fixed particle number $\mathcal{N}$. A constraint-free variational problem can be formulated with the help of a Lagrangian multiplier $\Lambda$:

$\displaystyle\mathcal{L}=\mathcal{F}+\Lambda\,\mathcal{N}=\int_{-\infty}^{+\infty}\left(\frac{1}{2}\,\left({E_{z}}^{2}-\left(\left|\xi\right|+\frac{\mathcal{N}}{2}\right)^{2}\right)+n_{\mathrm{e}}\bigl{(}\log(n_{\mathrm{e}})-1+\Lambda\bigr{)}\right)\,\mathrm{d}\xi.$

(40)

The first variation of $\mathcal{L}$ reads as follows, where the second identity is obtained by partial integration of the first term, taking Poisson’s equation for $\delta E_{z}$ into account:

	$\displaystyle\delta{\mathcal{L}}=$	$\displaystyle\int_{-\infty}^{+\infty}\left(E_{z}\,\delta{E_{z}}+\delta{n_{\mathrm{e}}}\,\log(n_{\mathrm{e}})+\Lambda\,\delta{n_{\mathrm{e}}}\right)\,\mathrm{d}\xi,$		(41)
	$\displaystyle=$	$\displaystyle\int_{-\infty}^{+\infty}\delta{n_{\mathrm{e}}}\bigl{(}-\phi+\log(n_{\mathrm{e}})+\Lambda\bigr{)}\,\mathrm{d}\xi.$

The functional $\mathcal{L}$ has an extremum when $\delta{\mathcal{L}}$ vanishes for all $\delta n_{\mathrm{e}}$. This implies that the bracket in the second identity vanishes. The electrons must be in Boltzmann equilibrium. The fact that $\Lambda=0$ is a consequence of the choice of units:

$\displaystyle n_{\mathrm{e}}=\exp(\phi-\Lambda)\equiv\exp(\phi)$

(42)

The second variation (43) directly shows that all extrema of the functional are minima, the equilibria are therefore stable:

$\displaystyle\delta^{2}{\mathcal{L}}$

$\displaystyle=\int_{-\infty}^{+\infty}\Biggl{(}{\delta{E_{z}}}^{2}+\frac{{\delta{n_{\mathrm{e}}}}^{2}}{n_{\mathrm{e}}}\Biggr{)}\,\mathrm{d}\xi\geq 0.$

(43)

To assess how quickly the equilibrium is restored, we analyze the evolution of a small perturbation of a stationary solution of equation (36). We assume that the perturbation is current-free and has an integrated density of zero. Expressed in the electrical field $\delta E_{z}$, the evolution can be described by a parabolic differential equation, completed by suitable asymptotic conditions and an initial condition at $t=0$:

	$\displaystyle\frac{\partial\delta E_{z}}{\partial t}=\frac{\partial^{2}{\delta{E_{z}}}}{\partial{\xi}^{2}}+\bar{E}_{z}\frac{\partial\delta E_{z}}{\partial\xi}-\exp\left(\phi\right)\delta{E_{z}},$		(44)
	$\displaystyle\delta E_{z}\big{|}_{\xi\to\pm\infty}=0,$
	$\displaystyle\delta E_{z}\big{|}_{t=0}=\delta E_{0}(\xi).$

It is advantageous to interpret the dynamics in the language of functional analysis. We consider the space of all complex functions $\psi$ on $\mathbb{R}$ which have a finite norm that is motivated by the second variation of the Helmholtz free energy,

$\displaystyle\mathbf{H}=\left\{\psi\;\Big{|}\;||\psi||^{2}=\int_{-\infty}^{\infty}\Biggl{(}|\psi|^{2}+\exp\left(-\phi\right)\left|\frac{\partial\psi}{\partial\xi}\right|^{2}\Biggr{)}\,\mathrm{d}\xi<\infty\right\}.$

(45)

This space can be turned into a Hilbert space by defining the scalar product

$\displaystyle\langle\psi,\psi^{\prime}\rangle=\int_{-\infty}^{\infty}\Biggl{(}\psi^{*}\psi^{\prime}+\exp\left(-\phi\right)\frac{\partial\psi^{*}}{\partial\xi}\,\frac{\partial\psi^{\prime}}{\partial\xi}\Biggr{)}\,\mathrm{d}\xi.$

(46)

The evolution operator is now defined as

$\displaystyle T\psi=-\frac{\partial^{2}{\psi}}{\partial{\xi}^{2}}-\bar{E}_{z}\frac{\partial\psi}{\partial\xi}+\exp\left(\phi\right)\psi.$

(47)

A short calculation verifies the identity

$\displaystyle\langle T\psi,\psi^{\prime}\rangle=\int_{-\infty}^{\infty}\Biggl{(}\exp(-\phi)\frac{\partial^{2}\psi^{*}}{\partial\xi^{2}}\frac{\partial^{2}\psi^{\prime}}{\partial\xi^{2}}+\left(3-\exp\left(-\phi\right)\right)\frac{\partial\psi^{*}}{\partial\xi}\frac{\partial\psi^{\prime}}{\partial\xi}+\psi^{*}\psi^{\prime}\Biggr{)}\,\mathrm{d}\xi,$

(48)

from which we conclude that the evolution operator is self-adjoint and positive definite, so that its spectrum lies on the positive real axis:

$\displaystyle\langle T\psi,\psi^{\prime}\rangle=\langle\psi,T\psi^{\prime}\rangle.$

(49)

Fig. 7 shows the first eigenvalues and eigenfunctions for the case $n_{0}=0.75$. Evidently, one eigenvalue is $\lambda_{1}=1$, the corresponding eigenfunction is

$\displaystyle\psi_{1}(\xi)=\exp\left(\phi(\xi)\right)\sim\frac{\partial\bar{E}_{z}}{\partial\xi}-1.$

(50)

The following identity (established by means of partial integration) implies that all other eigenvalues are larger than unity.

	$\displaystyle\langle T\psi,\psi\rangle-\langle\psi,\psi\rangle=$	$\displaystyle\int_{-\infty}^{\infty}\exp(\phi)\left(\frac{\partial^{2}}{\partial\xi^{2}}\left(\exp\left(-\phi\right)\psi\right)\right)^{2}\,\mathrm{d}\xi$		(51)
		$\displaystyle+\int_{-\infty}^{\infty}\exp\left(\phi\right)\left(2-\exp\left(\phi\right)\right)\left(\frac{\partial}{\partial\xi}\left(\exp\left(-\phi\right)\psi\right)\right)^{2}\,\mathrm{d}\xi\geq 0.$

Written in this compact functional analytic notation, the evolution equation for the perturbation and its initial condition read:

	$\displaystyle\frac{\partial\delta E_{z}}{\partial t}+T\,\delta E_{z}=0,$		(52)
	$\displaystyle\delta E_{z}\big{|}_{t=0}=\delta E_{0}.$

An application of the Laplace transform yields

$\displaystyle p\underline{\delta E_{z}}+T\underline{\delta E_{z}}=\delta E_{0},$

(53)

which can be solved in terms of the resolvent:

$\displaystyle\underline{\delta E_{z}}=\left(p+T\right)^{-1}\delta E_{0},$

(54)

The Laplace back transform gives an explicit solution:

$\displaystyle\delta E_{z}(t)=\int_{-i\infty}^{i\infty}\left(p+T\right)^{-1}\exp(pt)\,\mathrm{d}p\,\delta E_{0}.$

(55)

Utilizing that all eigenvalues of $T$ are positive, the answer to the perturbation is a sum of decaying exponential modes:

$\displaystyle\delta E_{z}(t)=\sum_{n=1}^{\infty}\langle\Psi_{n},\delta E_{0}\rangle\,\exp\left(-\lambda_{n}t\right)\,\Psi_{n}.$

(56)

It can be stated, that the system returns to its equilibrium in an exponential fashion.The speed is governed by the dielectric relaxation time $\tau_{\mathrm{e}}$. For the simulation of section 4, the dielectric relaxation time is $\tau_{\mathrm{e}}\approx 38.7\,$ns, faster than the RF period ($T_{\mathrm{RF}}\approx 73.7\,$ns). The slowest time constant $\lambda_{1}=1$ corresponds to a motion of the soliton as a whole: If displaced, the soliton returns to its original position within a few relaxation times. This is already captured by equation (34). The higher modes govern the restoration of the shape and act even faster. When distorted, for example by contact with the walls, the soliton assumes its orginal shape quite quickly.

The previous sections have described the non-neutral discharge regime in terms of a hybrid PIC/MCC simulation (section 4) and a simplified analytical model (section 5). How do the descriptions compare? For reference, we extract from our hybrid PIC/MCC simulation at each time $t$ the ratio $n_{0}=n_{\mathrm{i}}/n_{\mathrm{e}}$ and calculate the solution of equation (36). All quantities are transformed back into physical units and the laboratory frame.

Figure 8 compares the two models in terms of the temporally and spatially resolved electron density $n_{\mathrm{e}}$, electric potential $\phi$, and electric field $E_{z}$. The left column presents spatially and temporally resolved results of the analytical model. The electron density $n_{\mathrm{e}}$ is shown in panel a), the electric potential $\phi$ is shown in panel d), and the electric field $E_{z}$ is shown in panel g). The mid column displays temporally and spatially resolved simulation results for the same quantities (b): $n_{\mathrm{e}}$, (e): $\phi$, (h): $E_{z}$. The right column presents temporal snapshots of the comparison between the analytical model and the simulation results (c): densities, (f): $\phi$, (i): $E_{z}$. Overall, an excellent qualitative and a satisfactory quantitative agreement is observed. This agreement and details of the following discussion are best seen when looking at the animated version of figure 8 given in the attached movie 4. Deviations are most pronounced when the soliton-like structure approaches either of the surfaces. All deviations are plausibly explained by the fact that the analytical model assumes a constant ion density and neglects the electron absorption at the walls. The assumption of a constant ion density leads to lesser accurate numerical values for the electric potential $\phi$ and field $E_{z}$. The neglect of wall interactions causes the modeled soliton to get stretched at the walls (c.f., fig. 8 a), b), and c)).

A detailed analysis of figure 8 reveals that the difference of the model behaves as one would expect based on Poisson’s equation (10). The electron density $n_{\mathrm{e}}$ and the electric field $E_{z}$ are connected by one integration, and the electric potential $\phi$ is calculated by integrating $n_{\mathrm{e}}$ twice. Each integration by its mathematical nature results in a smoother quantity and differences between the integrands become increasingly irrelevant. Accordingly, the electron density shows both qualitatively and quantitatively the strongest deviations from the simulation results (c.f., fig. 8 c)), and the electric potential $\phi$ bears the best agreement (i.e., there is just a slight quantitative deviation - c.f., fig. 8 i)). Moreover, the qualitative agreement between the temporally and spatially resolved profiles of the electric potential $\phi$ and the field $E_{z}$ (compare fig. 8 d) to e) and fig. 8 g) to h)) makes the respective contours nearly indistinguishable.

In summary and supported by both simulation and theoretical analysis, the following important physical statements can be made:

• •

  There is indeed a non-neutral regime of RF-driven plasma jets. In the whole system, the phase-averaged electron density is considerably smaller than the ion density. A quasi-neutral bulk does not develop. Instead, the electrons organize themselves in the form of a soliton-like structure with a maximum that lies below the ion density. This reflects that the electron soliton is not much larger than the Debye length, but scales with it.

• •

  The observed soliton structure is stable. In the simulation, this is evident from the fast recovery of the form after distortions by the boundary. In the analytical model, this was formally proved by means of stability analysis. Physically, the soliton is (in the co-moving frame) an example of a Poisson-Boltzmann structure.

• •

  Due to being related to the dynamics of an equilibrium point (i.e., the COM $\mathcal{Z}(t)$), the soliton dynamics govern the electron dynamics. The previously discussed fast recovery from distortions means that electrons rapidly close up to the soliton.

• •

  One remarkable prediction of the simulation is the weak spatial gradient of the electron temperature. The assumption of a spatially constant $T_{\mathrm{e}}$ is thus well justified. The analytical model assumes that $T_{\mathrm{e}}$ it is also temporally constant; this is clearly an oversimplification that deserves attention.

In conclusion, the physical statements extracted from the model answer the previously raised key questions: (i) How does the soliton-like structure form?, and (ii) Why do the electrons behave in a uniform fashion?

The purpose of this study was to describe and analyze the non-neutral regime of capacitively coupled atmospheric pressure plasma jets via the combined use of numerical simulation and analytical modeling. In the simulation section, a self-developed hybrid PIC/MCC code was employed. It was found that the discharge does not develop a stationary quasi-neutral bulk. Instead, the electrons form a mobile Gaussian-shaped structure that we termed a soliton. A characteristic feature of this soliton is its ability to recover its form when perturbed by an interaction with the walls. Another peculiar observation was that the electron temperature within the soliton is temporally modulated but spatially constant.

For a better insight into the underlying dynamics, we analyzed a simplified and thus transparent model of the electron dynamics. In a current-free (i.e., co-moving) frame, the equilibrium solutions (solitons) were found to follow a Poisson-Boltzmann equation. Employing Lyapunov’s second argument, we could show that the solitons are stable. Linear stability analysis allowed to investigate the dynamic behavior of perturbations of the equilibrium solution. We found that all distortions decay on timescales comparable to the dielectric relaxation time, way below the RF period. The Boltzmann equilibrium character of the solitons explains also the observed uniform heating.

The work closed with a comparison between the results of the hybrid PIC/MCC simulation and a parameter fitted solution of the analytical description. Despite its simplified nature, the analytical model performed quite well. Deviations between the models could be directly traced back to the simplifying assumptions of the analytical model. In the near future, we plan to expand our work on the non-neutral discharge regime. We plan to establish precise criteria for the non-neutral regime, to distinguish it from the classical discharge modes [24, 26, 36, 37]. During this mode transition, close attention will have to be paid to classical concepts such as the boundary sheath and quasi-neutrality. Both have no meaning within the non-neutral regime, and preliminary evidence suggests that their definition will be tricky in the transition regime when $n_{0}$ approaches unity. Also the other limit poses open questions: The analytical model yields solutions for arbitrarily small parameters $n_{0}$. Is there a lower limit to the value of $n_{0}$ below which stable discharge solutions cannot exist?

As a final remark, we stress that we would be interested in experimental validation of our theoretical speculation. To the best of our knowledge, there are currently no published experimental studies that confirm the existence of a non-neutral discharge regime of RF driven atmospheric pressure plasma jets.

Funding by the German Research Foundation DFG, project 327886311 (CRC 1316), and the National Research, Development and Innovation Office (Hungary), project K134462, is gratefully acknowledged.