Simulation and modeling of radio-frequency atmospheric pressure plasmas in the non-neutral regime

Radio-frequency-driven atmospheric pressure plasma jets (RF APPJs) play an essential role in many technological applications. This work studies the characteristics of these discharges in the so-called non-neutral regime where the conventional structure of a quasi-neutral bulk and an electron depleted sheath does not develop, and the electrons are instead organized in a drift-soliton-like structure that never reaches quasi-neutrality. A hybrid particle-in-cell/Monte Carlo collisions (PIC/MCC) simulation is set up, which combines a fully kinetic electron model via the PIC/MCC algorithm with a drift-diffusion model for the ions. In addition, an analytical model for the electron dynamics is formulated. The formation of the soliton-like structure and the connection between the soliton and the electron dynamics are investigated. The location of the electron group follows a drift equation, while the spatial shape can be described by Poisson–Boltzmann equilibrium in a co-moving frame. A stability analysis is conducted using the Lyapunov method and a linear stability analysis. A comparison of the numerical simulation with the analytical models yields a good agreement.


Introduction
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 an insignificant amount of gas heating that proves essential in the treatment of biological surfaces in medical applications [10][11][12][13] and other heat-sensitive materials [1][2][3]6].
The above-mentioned reasoning is proper for both atmospheric pressure and low-pressure plasmas. However, lowpressure plasmas require a sophisticated vacuum setup while the most simple atmospheric pressure plasma can be realized with just two electrodes and a sufficient power supply. This difference 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 nonequilibrium plasmas have specific difficulties too. Most notable in this context is, due to the high collisionality, the prevalence to transit to a local thermodynamic equilibrium [3,4,6]. Numerous different plasma source concepts and even more implementations that take countermeasures against significant gas heating have been used to circumvent this difficulty [3,6,15].
Standardized plasma sources are required to make it easier to compare experiments. 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 deduction, 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 λ D compared to the discharge length L. Boundary effects (e.g. sheath formation) take place on the length scale of λ D ≈ 10 −1 -10 −2 mm (estimated for the expected parameters T e ≈ 1 eV and n e ≈ 10 15 -10 17 m −3 ) 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. Although there are some similarities to Townsend discharges, we chose to call this discharge regime where no quasi-neutral plasma bulk is formed a nonneutral discharge regime. The Townsend regime was reported as operation modes in direct-current (DC) and dielectric barrier discharges (DBDs). For the former device, the discharges are known as dark Townsend discharge [1,5], and for the latter, the so-called Townsend mode is the simplest operation mode [29][30][31]. However, the key characteristics of these discharges is a vacuum-like electric field that is not significantly disturbed by the presence of the charged particles. The nonneutral regime does not share this feature. Later sections (cf, sections 4 and 5) will elaborate on this comparison and the disparities.
Although we suspect to see a non-neutral discharge within a He/O 2 mixture in previous work [25], a non-neutral regime of the COST-Jet has not been analyzed in detail. Employing a one-dimensional hybrid particle-in-cell/Monte Carlo collisions (PIC/MCC) simulation, we explore a non-neutral discharge for the He/N 2 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. The prevalence of the structure distincts the non-neutral regime from other yet reported discharge modes like the γ, Ω or Penning mode [32][33][34][35][36].
We, according to the structure's tendency to keep and quickly regain its shape, characterize it as drift-soliton-like. The expression soliton describes wave phenomena with a high degree of dimensional stability that occur in many fields of physics [37,38]. In the context of plasma applications, actual drift-solitons have mainly been reported for magnetized plasmas [39,40]. We find the plasma properties such as the electron temperature T 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 understanding of the formed soliton. Therefore, an abstract and simplified yet analytical model for the electron dynamics is formulated. A detailed mathematical 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 [41][42][43] 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 accompanied 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 as a tool for diagnostics. The final section summarizes our findings, draws conclusions, and outlines perspectives for future work.

The COST reference jet
The COST-Jet is a capacitively coupled radio-frequencydriven 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 [15,[22][23][24][32][33][34][35]. 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 [22][23][24][32][33][34][35].
Although we firmly believe that a complete model of the Cost-Jet, a long tube (30 mm) with a tightly bounded 1 mm × 1 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 [22][23][24][32][33][34][35], we focus on describing the plasma and especially the electron dynamics between the electrodes along the z-axis and perpendicular to the gas flow. Moreover, these studies have shown that a one-dimensional simulations including a fully self-consistent and kinetic treatment of the electrons as presented in the following section reproduce all features observed in experiments [32][33][34][35]. Figure 1(b) shows a crosssection 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, 22-26, 33, 34], less commonly argon [18,52]) mixed with a molecular trace gas (usually nitrogen [24,33,34,52] and/or oxygen [15,18,19,22,23,25,26,52,53]). 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 and the applicability of the numerically demanding particle-in-cell (PIC) method. Both will be presented in following sections. The details of the chemistry set used for this study are mainly oriented on previous work [28] and are shown in table 1. The table reveals that except of the Penning ionization (No. 31) no follow-up reactions need to be considered for the helium/nitrogen mixture. This means that no other reaction requires its reactants to be formed previously in another reaction (e.g. the helium metastable He * needs to be formed by reaction No. 2 or 3 to than take part in reaction No. 31). A quick estimate (lifetime τ m ≈ 10 ns, flow velocity of the gas stream v f ≈ 100 m s −1 ) of the average distance a helium metastable drifts during its lifetime ends up in the range of 1 μm. Thus, the accumulation of helium metastables along the gas flow axes is a negligible effect. The one-dimensional representation of the jet's plasma is valid as long as the cooling effect of the gas flow is considered (e.g. by fixing the temperature of the simulated gas background).

Hybrid particle-in-cell/Monte Carlo collisions simulation
The PIC simulation is a numerically sophisticated tool that provides self-consistently calculated statistical representations of the distribution functions [54][55][56][57]. The foundation of the PIC algorithm traces back to the 1940s [54]. In the 1960s, a Monte Carlo collisions (MCC) scheme was added that lead to the commonly used PIC/MCC simulation [55]. In recent years, the PIC/MCC simulation has had great success in the research of non-equilibrium plasmas in general [1,[56][57][58][59], and the scheme has been applied in the atmospheric pressure regime [28,33,34]. 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,58,60].
Overall, hybrid concepts in plasma modeling are commonly used. The review paper of van Dijk et al [58] lists numerous techniques, variations, and applications. The hybrid PIC simulation has many applications in astrophysical plasma studies [61][62][63][64]. 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 PIC/MCC code Eehric. The following subsections present some details of the implementation.

Description of the electrons
As in the classical PIC/MCC scheme, the electrons are described as particles, represented by so-called superparticles. Their dynamics are traced by solving the equations of motion individually for each superparticle [54][55][56][57]. The code Eehric applies the Leap-Frog algorithm on a regular one-dimensional Cartesian discretization of the equations of motion. The width Δz of the resulting spatial grid is coupled to the Debye length by λ D 2 Δz to avoid numerical instabilities [55,57,59]. For the Cost-Jet geometry, the small electrode separation allows a low number of grid points to satisfy the requirement. We fixed the number of grid points for this study to 201.
The usual Monte Carlo technique is implemented to account for collisions [55,57]. Additionally, a null-collision technique [65] enhances the numerical efficiency of the collision routine. The necessary cross sections for the electron processes listed in table 1 are retrieved from the website of the LXCat project [44][45][46]. For helium, the cross section data originate from the database Biagi-v7.1 [47] based on version 7.1 of the simulation program Magboltz [66]. For the nitrogen cross sections, the Phelps database [48] provides data based on data by Phelps and Pitchford [49]. In terms of this study, as well as in previous work [28], the effects of gas heating are neglected.
The collision probability for electrons P e is found by where ν e is the electron collision frequency and Δt the discrete timestep [28,55,57,59]. The criterion for stability Δtω pe 0.2 that requires the time step Δt to sufficiently resolve the electron plasma frequency ω pe is overshadowed by the criterion ν e Δt 1 [55,57,59]. The latter criterion follows from a requirement that forbids particles to suffer consecutive collisions per time step. Under atmospheric pressure conditions, this criterion enforces a much smaller time step than the stability requirement. We use 1.2 million time steps for each RF cycle to limit the probability for a second (i.e. missed) collision per time step to approximately 1 percent. The results shown in later sections are averaged over 2000 RF periods and binned to 1000 diagnostic time steps within an RF period. This averaging process allows for satisfactory results with a relatively low number of about 10 000 superparticles. As Erden and Rafatov [67] present, the influence of the number of superparticles on the statistics always has to be considered.

Description of the ions
Eremin et al present that a drift-diffusion approximation is a valid representation of the ion dynamics at high pressures [27]. Thus, Eehric adapts this insight and simulates the He + , He 2 + , and N 2 + ions by evaluating the following form of the one-dimensional drift-diffusion approximation: The index s represents the individual ion species. n i is the ion number density, Γ i is the ion flux density, μ i is the mobility constant, E z denotes the electric field, and D i is the mobility constant.
As argued in previous work [28], the trace gas admixture is usually one percent or lower. Thus, the influence of nitrogen on the gas transport is negligible, and the binary diffusion coefficients and ion mobility in helium yield a sufficiently accurate description. The necessary mobility data have been measured by Frost [68] (He + /He), Beaty and Patterson [69] (He 2 + /He), and McFarland et al [70] (N 2 + /He). The numerical values for He 2 + /He and N 2 + /He are retrieved from the data collection of Ellis et al [71]. Values for the binary diffusion coefficients have been obtained by Deloche et al [72] (He + /He and He 2 + /He) and Walker and Westenberg [73] (N 2 + /He). For the numerical solution of equations (2) and (3), the scheme presented by Scharfetter and Gummel [74] is applied. Other authors [75] refer to this scheme as exponential scheme due to its exponential nature that provides an inherent switching between upwind and downwind differencing.
Ion-induced secondary electron emission is included by basically counting the number of secondary electrons N SE,r/l at each electrode separately and releasing electrons as long as N SE,r/l 1. To account for an in time average sufficient number of secondary electrons, the decimal part of N SE,r/l are treated by a random process. Whenever the Table 1. Plasma chemical reactions considered in the simulation. The last column uses the following abbreviations for the data sources of the corresponding cross sections or reaction rates: (1) data from the Biagi-v7.1 database obtained via the website of the LXCat project [44][45][46][47], (2) data based on [49] obtained from the Phelps database via the website of the LXCat project [44][45][46]48], (3) reaction rates retrieved from Brok et al [50] and Sakiyama and Graves [51]. The second column from the right gives threshold energies ε thr for electron reactions and reaction rates k r for ion reactions.
equally distributed pseudo-random number in the interval [0, 1) R 01 N SE,r/l , a secondary electron is generated at the respective electrode (left/right). The total number of secondary electrons is calculated by w e is the weight of the electron superparticles, γ i,s the ioninduced secondary electron emission coefficient, and Γ i,el,s the ion flux density towards either the right or left electrode surface. Whenever a secondary electron is launched from the surface, the respective counter N SE,e/l is reduced by one or for fractions reset to zero. The coefficients are estimated based on an empirical formula given by Raizer [5] that finds application in recent work [36,60]. The formula reads with ε thr,s the ionization threshold energy of species s (cf table 1), and ε φ the work function of the electrode surface material (cf Wilson [76]). For the He/N 2 mixture at hand, the relevant coefficients approximate to γ i,He + = 0.25, γ i,He 2 + = 0.25, and γ i,N 2 + = 0.1. These values may be rough estimates. However, the results show that the nonneutral discharge regime in insensitive towards a slight arbitrariness of the secondary electron emission coefficients (cf, section 4).
To account for Penning ionization according to process no. 31 (cf table 1), a procedure similar to previous work [28] is adapted. It is assumed that for both excitation reactions of helium (cf reaction no. 2 and 3, table 1), about 50 percent result in the formation of the metastable states He(2 1 S) and He(2 3 S), respectively. The time constant for diffusion is small compared to the lifetime of the metastables. Thus, diffusion is neglected, and a metastable atom is created by saving the position x m and generation time t m of the according excitation event. Each metastable is assigned with a lifetime k r is the reaction rate for the Penning ionization retrieved from [50,51], n N 2 denotes the neutral gas density of nitrogen, and R 01 is a pseudo-random number on the interval [0, 1). For t t m + τ m , the metastable decays via the Penning process by contributing to the ionization at position x m .

Simulation results and discussion
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 e (blue), the total ion density n i,tot (red) with spatial resolution. The panels (b) (n e ), (c) (n N 2 + )), and (d) (n i,tot ) provide spatially and temporally resolved representations of the respective number densities. First, the figure shows that the N 2 + ions are the clearly dominant ion species for this case (figures 2(c) and (d)). Panels (c) and (d) show no by naked eye visible difference. The density of He 2 + ions already is about three orders of magnitude lower than that of the N 2 + ions. For this case, the density of He + ions, which is even lower than the density of the He 2 + ions. Hence, both ion densities are not shown. The extremely low He + ion density stems from the efficient ion conversion He + → He 2 + (cf reaction 32 table 1). This result is in accordance to Martens et al [77]. Both the He 2 + and the N 2 + ion density profiles exhibit only a slight temporal modulation. The plasma is operated in the RF regime (ω p,i ω RF ω 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 e is approximately Gaussian-shaped (figure 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 (figure 2(b)). This behavior reminds us on two of the three key characteristics of solitons as given by Drazin and Johnson ([37], 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 discharge 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 [33,34] cannot be used in the context of a non-neutral regime.
The black lines in figures 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 δ z calculates from the second spatial moment. The black lines mark the interval Z ± 1.28δ 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 φ (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 (φ: figure 3(c), E z : figure 3(d)). Due to not having a quasineutral region, the potential of the discharge is, in first approximation, the electrical potential of a parallel plate capacitor bearing a positive space charge [78]. Therefore, the spatial profile of the electric potential resembles a compressed parabola ( figure 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. figure 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. figure 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 (figures 3(b) and (d)). This distinct field structure shares analogies to the dark Townsend discharges in DC discharges [1,5] and the Townsend mode of DBDs [29][30][31]. Both of these discharge regimes are characterized by an electric field barely disturbed by the charged particles. Yet, the plateau shape similar to classical discharge modes [36,57] stems from the interaction of charged particles with the field. Since we find characteristics of both the Townsend regime and 'classical' RF discharges in our non-neutral regime, we define the non-neutral discharge regime as an intermediate regime in between.
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 2 + ions: (c), N 2 + ions: (d) and the electron temperature T 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 e from the diagonal element p zz of the pressure tensor, which is defined as with m e the electron mass, n e the electron number density, v 2 z the electron's mean-square velocity, and u the mean velocity of the electrons. The main advantage of this method is that no assumptions on the electron energy distribution function are required. Details of these diagnostics are found in previous work [57]. 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 ( figure 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 ( figure 4(a)). Furthermore, the moments when the movement speed of the electron group drops to zero (t ≈ 35 ns and t ≈ 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 ( figure 4(a)) completely overshadows the energy delivered to the ions (figures 4(c) and (d)). The power absorption by ions (figures 4(c) and (d)) basically mimics the particular structure of the electrical field ( figure 3(b)). For the positive half-wave of the RF voltage (t ≈ 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 ≈ 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 2 + and N 2 + 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. figures 3(b) to 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 ( figure 3(b)) and the maxima of the ion energy absorption (figures 4(c) and (d)). At t ≈ 20 ns (grounded) and t ≈ 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.  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 discharge 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 (table 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 ( figure 5(a)) coincide with the patterns of electron power absorption ( figure 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 ≈ 35 ns, at grounded electrode: t ≈ 70 ns), electrons get lost, lose energy, and do not ionize. The light green structures adjacent to the respective electrodes (grounded: t ≈ 10-25 ns, powered: t ≈ 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 ( figure 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 (figure 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 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 discharge regime? And (ii) why do the electrons behave similarly?

Mathematical modeling
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 [58]. 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 i as a spatial and temporal constant. The electron model is given by the equation of continuity with the generation term neglected, The reason for this neglect is a time scale argument that can be deduced from figure 5. There is approximately one percent of the density created during one RF-cycle. Therefore, sources are negligible for the established plasma on the time scale of an RF period. The flux Γ e is calculated with the diffusion constant D e and the mobility μ e , The electric field E z (z, t) and its potential φ(z, t) are given by Poisson's equation: 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 [79]: The system of equations is completed by suitable asymptotic conditions for z → ±∞. The electron density is assumed to vanish and the electric field is asymptotically linear; the constants E ±∞ have the dimension of the electric field: 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 N (t), the center of mass (COM) Z(t), and the soliton-averaged field E(t): Based on the properties of the above quantities, we introduce a comoving frame of reference, which turns out to be currentless. Furthermore, we simplify the equations by introducing a dimensionless notation. After a few intermediate steps, the original system (equations (8)-(11)) collapses to a second-order differential equation known as the Poisson-Boltzmann equation [80,81,[83][84][85]: The skipped mathematical details of the derivation are found in appendix A. The Poisson-Boltzmann equation is of second order and therefore needs two conditions. The symmetry around the origin demands ∂φ/∂ξ = 0. The second degree of freedom can conveniently by fixed by setting φ(0) = ln(n 0 ), with n 0 ∈ [0, 1]. This parameter can be identified as the ratio between the electron density and the ion density at ξ = 0. Analytical solutions to (17) are well known but quite inconvenient to handle [84,85]. Figure 6 presents numerical solutions for the electron density n e , the electric potential φ, 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 section 4).
Furthermore, the Poisson-Boltzmann equation displays a crucial distinction between this non-neutral discharge regime, the Townsend regime, and 'classical' RF discharges. The right-hand side of equation (17) accounts for the influence of charged particles on the electric field. Although the exponential term (i.e. the influence of electrons) might be neglected for very non-neutral scenarios (e.g. the blue curve of figure 6), there is always a non-zero contribution of charged particles (i.e. a significant disparity from a Townsend scenario). On the other hand, recent work has shown that equation (17) alone does not suffice to describe the whole discharge domain of a 'classical' RF discharge (i.e. it is part of a more sophisticated model [80,82]). Therefore, the non-neutral regime sets itself apart from both of the other discharges.
Intuitively, the stationary solutions are stable. An effective method to show this formally is Lyapunov's direct method [41][42][43]. The starting point is a balance equation for the density of the Helmholtz free energy [83,[85][86][87][88]: (18) The Helmholtz free energy itself reads as follows, where the temporally constant last term in the integrand serves to render the integral finite: The time derivative of F is negative semi-definite, so that it is a suitable candidate for a Lyapunov functional [88]. Note that the equality sign is only assumed when the flux and with it the dissipation vanishes: A solution of equation (17) is stable when it represents a minimum of the free energy under the constraint of a fixed particle number N . A constraint-free variational problem can be formulated with the help of a Lagrangian multiplier Λ: + n e (log(n e ) − 1 + Λ)) dξ.
The first variation of L reads as follows, where the second identity is obtained by partial integration of the first term, taking Poisson's equation for δE z into account: (E z δE z + δn e log(n e ) + Λδn e ) dξ, The functional L has an extremum when δL vanishes for all δn e . This implies that the bracket in the second identity vanishes. The electrons must be in Boltzmann equilibrium. The fact that Λ = 0 is a consequence of the choice of units: The second variation (24) directly shows that all extrema of the functional are minima, the equilibria are therefore stable: To assess how quickly the equilibrium is restored, we analyze the evolution of a small perturbation of a stationary solution of equation (17). We assume that the perturbation is current-free and has an integrated density of zero. Expressed in the electrical field δE z , the evolution can be described by a parabolic differential equation, completed by suitable asymptotic conditions and an initial condition at t = 0: It is advantageous to interpret the dynamics in the language of functional analysis. We consider the space of all complex functions ψ on R which have a finite norm that is motivated by the second variation of the Helmholtz free energy, (26) This space can be turned into a Hilbert space by defining the scalar product The evolution operator is now defined as A short calculation verifies the identity from which we conclude that the evolution operator is selfadjoint and positive definite, so that its spectrum lies on the positive real axis: Figure 7 shows the first eigenvalues and eigenfunctions for the case n 0 = 0.75. Evidently, one eigenvalue is λ 1 = 1, the corresponding eigenfunction is The following identity (established by means of partial integration) implies that all other eigenvalues are larger than unity.
Written in this compact functional analytic notation, the evolution equation for the perturbation and its initial condition read: An application of the Laplace transform yields which can be solved in terms of the resolvent: The Laplace back transform gives an explicit solution: Utilizing that all eigenvalues of T are positive, the answer to the perturbation is a sum of decaying exponential modes: It can be stated, that the system returns to its equilibrium in an exponential fashion. The speed is governed by the dielectric relaxation time τ e . For the simulation of section 4, the dielectric relaxation time is τ e ≈ 38.7 ns, faster than the RF period (T RF ≈ 73.7 ns). The slowest time constant λ 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 (A.18). 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 original shape quite quickly.

Simulation and model: physical conclusions
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 i /n e and calculate the solution of equation (17). 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 e , electric potential φ, and electric field E z . The left column presents spatially and temporally resolved results of the analytical model. The electron density n e is shown in panel (a), the electric potential φ 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 e , (e): φ, (h): E z . The right column presents temporal snapshots of the comparison between the analytical model and the simulation results (c): densities, (f ): φ, (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 φ and field E z . The neglect of wall interactions causes the modeled soliton to get stretched at the walls (cf, figures 8(a)-(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 e and the electric field E z are connected by one integration, and the electric potential φ is calculated by integrating n 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 (cf, figure 8(c)), and the electric potential φ bears the best agreement (i.e. there is just a slight quantitative deviation-cf, figure 8(i)). Moreover, the qualitative agreement between the temporally and spatially resolved profiles of the electric potential φ and the field E z (compare figures 8(d) and (e) and (g) and (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 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 e is thus well justified.
The analytical model assumes that T 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?

Summary and outlook
The purpose of this study was to describe and analyze the nonneutral discharge regime of capacitively coupled atmospheric pressure plasma jets via the combined use of numerical simulation and analytical modeling. In the simulation section, a selfdeveloped 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,33,34]. 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. However, we suspect that the distinct structure of the electric field may be key to achieve this experimental evidence. The electric field induced second harmonic generation (E-FISH) technique that was just recently successfully applied to the COST-Jet [89] seems as a promising candidate to set up the required experiment. Due to the discussed similarities of the non-neutral regime to the Townsend mode of DBDs, we expect to find experimental proof in the low-power range.

Acknowledgment
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.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.
Appendix A. Discussion of the current-free reference frame Integration of the particle balance equation (8) over the spatial axis leads to the insight that the particle number 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: Integrating the continuity equation with the weight z, we derive a differential equation that describes the drift motion of the