A new 2D fluid-MC hybrid approach for simulating nonequilibrium atmospheric pressure plasmas: density distribution of atomic oxygen in radio-frequency plasma jets in He/O2 mixtures

A spatially two dimensional fluid-MC hybrid (fluid-kinetic) simulation method is developed and applied to the COST reference microplasma jet operated in helium with an oxygen admixture of 0.5%, excited by a single frequency voltage waveform with f = 13.56 MHz and ϕrms=275 V. The simulation approach is based on a fluid model augmented by a Monte Carlo module that generates electron impact rates for the continuity equations solved by the fluid module. This method is capable of providing the same level of accuracy as PIC/MCC simulations with an agreement within 5%–10% at atmospheric pressure, while being significantly faster (with a speedup factor of 30 for serial to 50 for parallel implementation). The simulation results are compared to previous measurements of atomic oxygen densities (Steuer et al 2021 J. Phys. D: Appl. Phys. 54 355204), and show a very good agreement. It is found that the buildup and saturation of the atomic oxygen density distribution along the jet are due to the interplay of chemical and electron impact reactions as well as of the gas flow. Comparing the simulation results to that of Liu et al 2021 J. Phys. D: Appl. Phys. 54 275204, it is inferred that fluid models where a 2-term BE solver is used, fail to describe the COST jet in an accurate manner due to the underestimation of the electron impact rates.

There has been a wide range of simulation methods applied to µAPPJs, such as 0d models, which are able to provide volume-averaged densities and are useful if the number of reactions is high [23][24][25][26], fluid models in one [27][28][29] and two spatial dimensions [30][31][32][33], where fluid equations are solved with predefined transport coefficients and reaction rates.These are usually calculated from e.g., 0d Boltzmann solvers for electrons, either as a function of the local reduced electric field, E/N, called the mean field approximation, or as a function of the mean energy, ε, called the mean energy approximation [33,34].Fluid models provide reasonably accurate spatiotemporally resolved data for the species densities and fluxes.A widely-used simulation technique for µAPPJs is provided by hybrid models, which make use of the timescale separation present among the physical processes taking place in such systems [35,36]: here, different processes are described by different submodels in the simulation: this can be done (i) within a fluid model, where e.g., a time-slicing technique is used for the neutrals, i.e., they are advanced in time with a significantly larger time step as compared to charged species [37][38][39][40], or (ii) using a fluid model for heavy charged species and neutrals, but a full PIC/MCC model for the electrons (these models are also called fluid-PIC/MCC models) [41][42][43].In this case, electron transport coefficients are not needed, since the full spatio-temporal dynamics can be readily calculated for electrons.The latter approach is capable of providing accurate data, close to the level of PIC/MCC, but due to the high number of collisions and the correspondingly large number of time steps, it is computationally very expensive in two dimensions.A complete (kinetic) description is provided by 1d3v Particle-In-Cell/Monte Carlo (MC) Collisions simulations [44][45][46][47][48][49][50][51], which are, however, computationally very expensive.Accounting for kinetic effects in case of electrons is, however, important for an accurate description of the plasma, in particular, if energetic electron groups are present, as is the case in RF microplasma jets operated in helium mixed with reactive gases, such as N 2 and/or O 2 [40,45,47].In these cases Penning-ionization can create electrons within the sheath regions, which are then accelerated by the high local electric field, thereby altering the high-energy tail of the energy distribution function (EEDF).Since this change in the EEDF cannot be captured accurately by fluid models using a 2-term 0d BE solver, as shown in [39,40], the experimental observations can only be reproduced using significantly higher voltages as compared to those applied in the experiments.In other words, these models significantly underestimate the electron density due to their inability to capture the physical processes caused by electrons with energies near the high-energy tail of the EEDF [39].
In this work, we propose a hybrid approach that is able to describe electrons kinetically while at the same time offers a significant speedup compared to PIC/MCC simulations.In this approach the fluid model, where the continuity equation is solved for each species (electrons, charged heavy species as well as neutral species) based on the drift-diffusion approximation, along with Poisson's equation, is augmented by a MC module where the spatio-temporally resolved electron impact rates are calculated based on the (spatiotemporally resolved) electric field obtained from the fluid module.The calculated rates are then transferred to the fluid model where they act as sources for the respective continuity equations.The gas flow field is assumed to be a parabolic Poiseuille flow.The two modules (fluid and MC) are iterated until convergence is achieved.Additionally, the electron transport coefficients (diffusion coefficient and mobility) are calculated by another MC module as a function of the local reduced electric field, E/N from the same set of cross sections that is used in the calculation of the electron impact rates.From these values a lookup-table is generated, which is then used in the simulations.This way, a speedup factor of 30 (in case of a serial implementation with smoothing of the rates) to 50 (in case of a parallel implementation of the MC module on 2-3 threads) can be achieved in 1d compared to (serial) PIC/MCC simulations.The performance of two dimensional simulations cannot be compared due to the lack of PIC/MCC data, but based on the large increase in computational requirements of the PIC/MCC method with increasing geometric complexity, an even higher acceleration factor can be safely assumed.Note, that with appropriate changes (e.g., solving the momentum balance/energy equations instead of assuming the drift-diffusion approximation), this method is applicable to simulating low pressure discharges as well.
As an example, we investigate the atmospheric pressure COST reference microplasma jet [52][53][54] operated in a He/O 2 mixture excited by a sinusoidal voltage waveform.The formation of atomic oxygen in this plasma source has drawn much attention [55][56][57][58][59][60].Therefore, in this work, we compare Two Photon Absorption Laser Induced Fluorescence (TALIF) measurements of the two dimensional spatial distribution of atomic oxygen from [59] to results obtained from our hybrid approach.We show, that there is an excellent agreement between simulation and experiments, which suggests that the simulation method is able to capture the relevant physical processes leading to the generation of atomic oxygen.Furthermore, we analyze how the atomic oxygen density builds up along the jet, which has not been done in two dimensions on a kinetic level before.We explain why the atomic oxygen density is saturated after an initial buildup near the end of the nozzle.This will shed light onto the role of the proper description of electron kinetics in atmospheric pressure RF plasmas.
The paper is structured as follows: in section 2 the computational approach is described in detail.Then, in section 3 the comparison between experimental and simulation results is presented and discussed.Finally, in section 4, conclusions are drawn.

Description of the COST-jet and the computational method
The schematic of the jet and the corresponding spatially two dimensional simulation domain is shown in figure 1.The jet consists of two plane parallel electrodes with a gap of L x = 1 mm, where the powered electrode (at x = 0) is driven by radio-frequency voltage waveforms.The length of the jet is L z = 30 mm, at its inlet (z = −30 mm), a given gas mixture is blown through the jet with a predefined mass flow rate.The width of the jet (in the y-direction) is 1 mm, this third dimension is however, not resolved in our simulations, i.e., the spatial extent of the plasma source in the y direction is considered to be infinite.
From a computational perspective, in order to be able to simulate the buildup of various species along the jet, a spatially two dimensional simulation method is needed, which is able to provide a kinetic description of electrons [39,40].A simulation method such as PIC/MCC would be desirable, however, due to the accuracy criteria of this method (high number of timesteps within each RF-period due to the atmospheric pressure conditions and a high number of superparticles due to the high number of computational grid cells in two dimensions) makes it impractical due to the enormous computational cost [61,62].
Therefore, the computational method used in this work is based on a fluid-MC hybrid (i.e., fluid-kinetic) approach, which makes use of the fact, that the underlying physical processes take place on separate timescales: the fastest being the electron dynamics (ps-ns), followed by chemical reactions and ion dynamics (ns-µs), while the neutral dynamics constitutes the slowest processes present in the plasma (on the order of a few ms).In order to mimic this, the computational method consists of separate modules that communicate with each other through various physical quantities, as indicated in figure 2.
The main part of the simulation is a fluid module, where for the density of each species (charged, including electrons, as well as neutral), n, the continuity equation is solved, based on the drift-diffusion approximation: where R k is the source function for species k as a result of chemical reactions (for which the rate constants are input parameters), while S k denotes the source function as a result of electron impact processes.Note, that losses for species k as a result of chemical reactions are included into R k .u is the flow field of the gas, which is given as input to the fluid module.The flux, Γ, is, according to the drift diffusion approximation, given by [63]  Here q k denotes the charge of the species (which, for neutrals is 0), µ k and D k are the mobility and diffusion coefficient of species k.
The fluid module calculates the electric field as well, by solving Poisson's equation for the electric potential, ϕ: where N c denotes the number of charged species considered.The fluxes are calculated using the Scharfetter-Gummel scheme [64].The boundary conditions for the charged particles at the electrodes are n| x=0,x=Lx = 0 for negatively charged species and ∇n| x=0,x=Lx = 0 for positively charged species, i.e., the particle flux at the electrodes consists of only drift (for positively charged particles) or both drift and diffusion (for negatively charged particles).For neutrals, surface loss probabilities are used in order to determine the given boundary condition.At the inlet (z = −L z ) the particle densities are set to zero, while at the outlet ∇n| z=0 = 0 is used for all species, since the diffusion flux is negligible compared to the flux from the gas flow in the z-direction [41].Similarly, for the potential, ∇ϕ| z=0,z=−Lz = 0 at both the inlet and the outlet.The continuity equations for the charged species as well as Poisson's equation are discretized based on a Crank-Nicolson scheme and solved together using an implicit SLOR algorithm on a spatio-temporal grid [65,66].Due to the fact, that the electron and ion densities are coupled by the Poisson's equation, the fluid module has to resolve the RF-cycle, however, a few hundred timesteps within an RF cycle proved to be sufficient to resolve all relevant physical processes in the single frequency case.For the neutral species a time-slicing method is used [41]: after a few hundred RF-cycles in the fluid module, the equations for the neutrals are solved until convergence is achieved, while the charged particle densities are held constant.This convergence usually occurs after a simulated time of a few ms.
The gas flow velocity profile, u, is given to the fluid module as an input.For usual inflow rates (on the order of 1 slm) the flow is laminar and compressibility effects are negligible (since the Mach number is on the order of 0.01) [40].The heating of the gas is, to a first approximation, negligible in the gas mixture studied in this work (He/O 2 [67]).We can assume that the plasma does not disturb the flow field itself.Therefore, it is sufficient to calculate the flow field once, before the simulation is started.Due to the simplicity of the system, the flow field can be very well described by a (parabolic-profile) Poiseuille-flow [68], which was verified using OpenFOAM (v10), an opensource computational fluid dynamics package [69].
In order to account for the kinetic nature of electrons [40], a separate MC module is applied for the generation of electron impact rate functions (S k ) based on cross sections.In this module, the equations of motion are solved for the electrons based on an electric field, E, which is calculated by the fluid module.Collisions are handled based on a MC procedure [70,71], using (total/momentum transfer) cross sections.The electron impact rate, S k for the kth species is calculated as where n g,i is the density of the background species an electron collides with during process i, f e is the electron velocity distribution function, v the velocity of the electron, σ i the cross section of the ith process and ξ i,k is an integer factor denoting how many particles of species k are created during the process i.In order to get the full S k , one needs to sum over all such processes.Then, the resulting spatio-temporally resolved electron impact rates (not rate coefficients), S k (r, t) for each species k are given to the fluid module as shown in figure 2. Note, that other physical quantities, that could in principle be calculated by the MC module for the electrons are not used during the iteration process.
In order to properly describe the underlying physics in the MC module, surface processes originating from heavy charged particles that result in the generation of electrons, such as secondary electron emission, are calculated based on the ion fluxes (Γ i , cf. figure 2) taken from the fluid module.Additionally, electron reflection is also considered.In order to account for the generation/destruction of electrons within one RF-cycle due to chemical reactions, those rate functions, which result in electron creation/destruction are given to the MC module (as seen in figure 2).In this module, the number of timesteps within an RF cycle is typically a few million to keep the probability of collision low enough [47].In this work, the collision probability threshold was set to 5%.The two modules (fluid and MC) are iterated until convergence is achieved.Typically, after a few hundred RF-cycles in the fluid module, the MC module is run for 2-5 RF-cycles depending on the acceptance level of statistical fluctuations for the spatiotemporal distribution of S k .The number of superparticles in the MC module is changed according to the required statistics: typically, 10-20 particles per grid cell are sufficient.
This approach requires more input data than fully kinetic models, such as PIC/MCC simulations.Apart from the cross sections for the electron impact processes, chemical rate coefficients and surface coefficients, the transport coefficients of the various species (diffusion coefficients and, in case of charged species, also mobilities) are also needed.In order to make the approach as close to a kinetic simulation as possible, a separate MC code was used to calculate the mobility as well as the (longitudinal) flux diffusion coefficients for electrons based on a method suggested in [72].We assume that the local field approximation is valid and the transport coefficients are a function of the local reduced electric field, E/N [34].The transport coefficients (diffusion coefficient and mobility) are only calculated once for each gas mixture, prior any actual simulation of the discharge.Based on the values a lookup table is generated (as a function of E/N), which is then used in the simulations.Calculating the transport coefficients instead of taking them from external sources makes the simulation more consistent since the transport coefficients are calculated from the same cross section set that is used in the MC module.Whenever reliable cross sections are available for other charged particles, they are used to calculate the transport coefficients.For ions, in case there are neither cross sections nor transport coefficients available (such as in case of, e.g., O + 2 +He, see later), a Langevin cross section is assumed and thus a constant mobility is used [73].The polarizabilities used were 0.208 Å 3 for He and 1.562 Å 3 for O 2 [74].The diffusion coefficient is then calculated by the Einstein relation, based on a constant background gas temperature, T g .In calculating the transport coefficients, due to the very low admixture of reactive gases, pure helium is used as a background gas.
With this method, a speedup factor of 30 in case of a serial implementation with smoothing of the rates to 50 (in case of a parallel implementation of the MC module on 2-3 threads) can be achieved compared to PIC/MCC simulations in 1 spatial dimension, based on our benchmark with the (serial) code presented in [47] in a He/N 2 mixture.This speedup originates from various factors: (i) the fluid module is much faster than the corresponding kinetic simulations due to the significantly lower number of arithmetic operations required (the fluid module usually takes up only 20% of the total runtime), (ii) the MC module does not have to fulfill the accuracy constraint of PIC/MCC simulations, viz.that the number of superparticles be sufficiently high, since that is needed so that numerical fluctuations in the electric field, after solving Poisson's equation, is minimized [61].Since in this model the electric field is supplied by the fluid module, the only reason for increasing the number of superparticles is that the signal-to-noise ratio of the electron impact source functions, S, is increased.Finally, (iii) the MC module is easily parallelizable.Thus, this method is more easily scalable to two dimensions than conventional PIC/MCC simulation codes.The results obtained from both methods agree within 5%-10%, where the discrepancy mainly originates from the transport coefficients and the assumption of the local field approximation, which breaks down in the sheath region.We need to mention, however, that this method is only applicable to steady-state discharges, since, due to the time-slicing, i.e., that within the iterative process, neutrals are advanced several ms in time (even though it would only take ≈3 ms for the gas to go from inlet to outlet) it is not possible to exactly characterize the elapsed physical time (above one RF-cycle) in the simulation.
This work focuses on He/O 2 mixtures.The charged species considered in the simulation are e − , He Here He * refers to an aggregate state of the 2 1 S and 2 3 S metastable states of helium [47].The electron impact cross sections (4 for He and 16 for O 2 ), the chemical reaction set (consisting of 91 reactions), the neutral diffusion coefficients as well as the surface coefficients are those presented in [41].The surface processes include secondary electron emission coefficients for positive ions and electron reflection, with values of r = 0.5, γ He + = 0.2, γ O + 2 = 0.05 and γ O + = 0.1 [41].The mobility of He + is from [75].For the oxygen ions a Langevin cross section is assumed for the calculation of the transport coefficients.The background temperature is set to T g = 350 K [60].For the simulations, a rectangular, uniform spatial grid of 60 × 60 proved to be sufficient.The number of timesteps in the fluid module was set to 100 per RF-period for neutral and 300 per RF-period for charged species, while in the MC module to 2.8 × 10 6 .The number of superparticles for electrons was set to 20 per cell.The grid spacing along with the timesteps used in the fluid module satisfy the CFL-condition, and all relevant physical processes were found to be properly resolved.The simulation converged after ≈2000 (fluid) RF-cycles.The requirement for convergence was that the densities do not fluctuate more than 1% over ≈100 RF-cycles.Using 12 threads for the MC module, the simulation took ≈6 days on an Intel ® Xeon ® Gold 6132 CPU.

Results
In this section, simulation results for the density of oxygen atoms are presented and compared to measurements from [60] in a mixture of He with 0.5% O 2 , excited by a single frequency waveform with f = 13.56MHz and an RMS voltage of ϕ rms = 275 V (corresponding to 1 W in the experiment) for various helium mass flow rates.The experimental results were obtained using TALIF.
Figure 3 shows experimental (top row) as well as simulation results (bottom row) of the two dimensional spatial distribution of the atomic oxygen density in the case of three different helium mass flow rates: 0.2 slm (a, d), 0.6 slm (b, e) and 1.0 slm (c, f).
As seen in the figure, the general trends and the buildup of the atomic oxygen density are properly captured by the simulation.As the mass flow rate (and with this, the maximum value of the velocity, u) increases, the residence time of the gas in the active plasma volume of the jet decreases and thus there is less time for the respective chemical reactions to take place.The absolute values obtained from the simulations agree with that of the experiment within ≈30%, however, one has to note that the estimated experimental error is around ≈50% [60].Nevertheless, this already shows the advantage a hybrid simulation can provide over e.g., fluid codes: due to the proper description of the electron dynamics, the computed values are within the range of the experimental results, which for fluid codes is not the case [39].Although the buildup of the oxygen atom density along the z-direction is properly captured, one needs to note, that the density profiles in the x-direction are slightly different in the simulation and in the experiment: this is most likely due to the choices for the surface loss probability and diffusion coefficient of atomic oxygen, which are both input parameters [41].
In order to investigate the buildup of the O atom density along the jet (i.e., along the z-direction) in more detail, figure 4 shows the normalized atomic oxygen density distribution averaged over the x-direction, for the three different mass flow rate values considered, as a function of the z-coordinate, obtained from the simulation along with experimental values.This representation eliminates the effect of systematic experimental uncertainty and allows side-by-side comparison of the data.
The experimental and simulation data show excellent agreement.In the case of the lowest mass flow rate value, 0.2 slm, the atomic oxygen density first increases reaching a local maximum at around z=−22 mm, beyond which a slight, but steady decrease can be observed, which is verified by the experimental results.
The same behavior cannot be observed for higher mass flow rates, since due to the increased velocity, the residence time is decreased such that the density is still in the buildup phase along the entire length of the jet.Although the effect of the gas flow is easily understood, the other two processes (electron dynamics and chemical reactions) are needed to explain the shape of the density distribution.In the following, we will investigate the case with 0.6 slm mass flow rate.
In order to understand what role the charged particle dynamics plays in the formation of atomic oxygen, figure 5 shows the spatial distribution of the time-averaged total positive ion density (a), the total negative ion density (b) and the electron density (c) along the direction of the jet.Panel (d) shows the spatial distribution of the electron and total negative ion density between the electrodes at specific z-values, indicated by the vertical dashed lines in panels (b, c).The dotted lines correspond to the negative ions density, while the full lines indicate the electron density.
As seen in panels (a, b) of figure 5, the positive as well as the negative ion density decreases along the jet.Furthermore, from the fact that both density profiles 'shrink' along the x-direction (i.e., between the electrodes), it follows that the spatio-temporal distribution of the electric field also changes along the flow.The electron density distribution (c, d) has a more interesting behavior: in the first ≈5 mm of the jet, the  electron density along the x-direction shows characteristics that are typical of electronegative discharges: since the negative ion density is appreciable near the inlet as per panel (b), the electron density profile will develop 'electropositive edges', as seen in panel (d).As we move along the discharge, and the negative ion density decreases, the electron density will increase in the center between the electrodes, and the electropositive edge will gradually disappear.This suggests that the electronegativity of the discharge decreases along the flow.
In order to understand what effects lead to the change in the charged particle densities observed above, figure 6  As seen in panels (a, b), the spatio-temporal distribution of the electric field changes as one moves along the jet: near the inlet, the 'bulk' electric field (i.e., the field between ≈0.2 mm and ≈0.8 mm) is higher in panel (a) compared to panel (b), which suggests that the respective electron impact reaction rates will also be higher near the inlet.The vertical lines in panels (a, b) indicate the time instances where the bulk electric field is maximum, which will correspondingly lead to maxima in the electron excitation/ionization rates (as part of the socalled Ω-peak, as shown later).
Panel (c) shows the temporal average of the x-component of the electric field, which exhibits similar patterns as figure 5 previously: the region where the temporal average of the electric field is small, shrinks, as one moves along the jet, which leads to the shrinking of the density profiles in the x-direction in figure 5. Figure 6(d) shows the spatial distribution of E x at the time instances indicated by the dashed/full vertical lines in panels (a, b), respectively.In the region where the 'bulk' electric field is high (in the range of ≈ (0.4 mm, 0.8 mm) for the blue curves and ≈ (0.2 mm, 0.6 mm) for the green curves), the absolute value of the electric field in panel (a), i.e., at a position near the inlet, is higher, which leads to an increased electron power absorption dynamics and consequently, higher electron impact rates.The reason for the higher electric field is the depleted electron density, as shown in figure 5(c): a smaller electron density leads to a smaller plasma conductivity and thus a higher Ohmic electric field [45].
In order to investigate how the electronegativity changes along the jet, figure 7 shows the two dimensional spatial distribution of the electronegativity.Here the electronegativity is defined as β(x, z) = n ion,-(x, z)/n e (x, z), where the densities are averaged over time.As seen in the figure, the electronegativity decreases from ≈ 3 near the inlet to below 1 in the first 10 mm of the jet.Furthermore, since, based on figure 6(c) the bulk width decreases along the jet, the region of high electronegativity gradually 'shrinks' in this direction as well.
Note, however, that the electronegativity is not as high as in, e.g., low pressure discharges (which can be as high as 100 [76,77]) and thus its effect on e.g., the electron power absorption dynamics will not be easily visible.
Since, ultimately, electrons are responsible for the generation of various species through electron impact collisions, it is important to understand how electrons are created and how that leads to their spatial density distribution seen in figure 5(c).Figure 8 shows the source functions for electrons originating from chemical reactions (first row) and from electron impact processes (i.e., ionization of helium and oxygen, second row), as a function of x and t, i.e., between the electrodes over one RF-cycle (a, c) and as a function of x and z (b, d).
The patterns seen in (a), the chemical reaction source function of the electrons, are due to Penning-ionization of oxygen, which is the dominant chemical reaction for the creation of electrons [47]: Helium metastables are primarily created by electron impact excitation.Since the threshold energies for the metastable states are relatively high (≈20 eV), their creation will take place where electron power absorption is high.The same is true for panel (c), i.e., the electron impact source function for electrons, which mainly originates from the ionization of oxygen (with a threshold energy of ≈12 eV).The numbers in panel (c) denote the most important electron power absorption patterns in atmospheric pressure RF discharges [47]: the marks 1 and 2 denote the Ω-peaks, which originate from the high Ohmic electric field in the bulk, due to the low conductivity as a result of the high collisionality of the plasma (cf.figures 6(a) and (b) and the corresponding vertical lines).The peak with the mark 3 corresponds to the Penning-peak, where electrons (mainly originating from Penning-ionization) ionize the  background gas due to the high 'sheath' electric field, which is due to the space charge present in this spatio-temporal region (as shown in figures 6(a) and (b)).These patterns can be correlated with the local maxima seen in panel (a), since the same mechanisms are responsible for the generation of helium metastables, which, in turn, through Penning-ionization, participate in the generation of additional electrons.Note, that there is a slight shift in time between the maxima in panel (c) and those in panel (a): the reason for this is, that the helium metastable density has to build up first in order for Penning ionization to be appreciable.The spatial evolution of the chemical electron source function is shown in panel (b).In the first 5 mm, where the electronegativity is high, so is the absolute value of the source function.This is in light with the behavior of the electric field (cf.figure 6), where the Ohmic electric field is higher near the inlet due to the depleted electron density, which leads to an increased production of He * , and thus, indirectly, of electrons.
Beyond the first 5 mm, the source function reaches a relatively steady value, where the two local maxima are closer to each other due to the change in the electric field distribution.The same trend can be observed in panel (d), showing the spatial distribution of the electron impact ionization rate: since it is primarily determined by the electric field, it also verifies that the reason for its decrease is the decrease of the electric field itself, which is intimately connected to the decrease of the electronegativity ( and thus the increase of the plasma conductivity).Note, that both the chemical and electron impact rates are important to take into account, since they are of the same order of magnitude (although the chemical rate is 3 times higher than the electron impact rate).
In order to understand the reason for the change in the electronegativity, figure 9 shows the density profiles (averaged over time and the x-direction) of the charged particles and of He Since only O − is created by electron impact, through dissociative attachment (with a mean rate of ≈2.5 × 10 23 m −3 s −1 ), the densities of all other negative ions decay rapidly, due to the generation of atomic oxygen and ozone, which destroy negative ions.This is why after a certain time a decrease can be observed even for O − .
Based on this, the charged particle dynamics is greatly affected by some of the neutral species, in particular, ozone and atomic oxygen.The converse is not true, i.e., in most cases, charged particles (except for electrons) do not affect the neutral density buildup appreciably due to their low density compared to that of the neutrals.As we shall see, electrons are important in properly describing the neutral density profiles, however, due to the well-separated timescales, over the timescale of the flow, the electron impact source functions can be averaged over one RF-cycle, and thus they can be taken into account as a temporally averaged background (in the sense that temporal modulations within an RF-cycle do not have an immediate effect on the neutral dynamics).We note, that in order to test whether additional charged species affect the plasma characteristics, we included in our simulations O + 4 ions with 7 additional chemical reactions taken from [26].Although the final density of O + 4 was ≈10% of that of O + 2 , apart from decreasing the densities of the other positive ions so that quasineutrality is fulfilled, no other significant changes were observed compared to the data presented in figure 9. Therefore, in order not to complicate the model further, we decided not to include O + 4 as an additional species.In order to understand how atomic oxygen builds up along the jet, figure 10 shows the averaged (over both time and the x-direction) density profiles for all neutral oxygen species (a) and the spatial distribution of the chemical rates for the most dominant reactions involving atomic oxygen (b).
As shown in panel (a), the most important species are atomic oxygen, ozone and the two metastable oxygen species, O 2 (a 1 ∆ g ) and O 2 (b 1 Σ g ).As mentioned above, the electron impact rates are almost constant along the jet, and act as a background.Furthermore, for O 3 and O 3 (v), there are no electron impact processes considered, and thus their density profiles are entirely determined by chemical reactions.Based on the absolute values of the mean chemical rates involving the creation/destruction of oxygen atoms, the following ones proved to be dominant: IV.He whose rates are shown in figure 10(b) spatially resolved along the jet.Reactions I-II produce, while reactions III-VI destroy atomic oxygen.The electron impact processes act as a positive background (of the order of ≈10 25 m −3 s −1 ).Reaction I constitutes an almost constant positive source throughout the jet, since O( 1 D) is primarily produced by electron impact processes.Reaction II has the highest absolute value among the reactions considered.In the first ≈10 mm of the jet, atomic oxygen has a sharp increase due to electron impact processes and reaction I-II.However, as the atomic oxygen density is increased, reactions III-VI become important, especially reaction III, which leads to the production of ozone.This slows down the increase of the atomic oxygen density, and together with diffusion losses to the walls (which is not shown here) eventually leads to the saturation in the density profile as shown in figure 4. We note, that reaction II is also the dominant loss mechanism for O 2 (b 1 Σ g ), which leads to the local maximum of its density near z ≈ −25 mm in panel (a).
Figure 11 shows the two dimensional spatial distribution of the chemical reaction rate, R O (a) and electron impact rate, S O (b) of atomic oxygen creation (both quantities are timeaverages).In light of the discussion above, it is clearly visible in panel (a), that the total chemical rate reaches zero beyond z ≈ −25 mm near the electrodes, and beyond z ≈ −20 mm in the center (i.e., at x = 0.5 mm).This is due to the onset of reactions III-VI.The reason for the nonuniform distribution of R O as a function of x is reaction II: since O 2 (b 1 Σ g ) is primarily produced by electron impact processes, it will attain its highest value in the center between the electrodes.But as per reaction II, a higher density for this metastable molecule will contribute positively to atomic oxygen, and therefore, a negative R O is reached at a later position in the center between the electrodes.
As shown in figure 11(b), beyond z ≈ −25 mm, the electron impact rate for atomic oxygen creation attains a uniform value (as a function of z), and thus acts as a background source.Note, however, that in this case the contribution of the electrons is higher compared to that of the chemical reactions (which is the opposite of what was observed in case of electrons, cf.figure 8).This also shows, that without a proper description of the electron dynamics, the absolute values of certain species can be severely underestimated, and this is why a hybrid approach is desirable for simulating the COST jet [39].We note, that the results of the charged and neutral particle dynamics presented for the 0.6 slm case can easily be transferred to the cases with flow velocities of 0.2 slm and 1.0 slm, respectively.Since the change in the flow velocity means a change in the residence time, the results for these two cases look 'selfsimilar' in the sense that the same patterns would be visible only at different z-coordinates, depending on the flow velocity (larger z for larger velocities).

Conclusions
In this work, a new hybrid (fluid-kinetic) simulation approach was introduced for the description of the COST microplasma jet.The approach is based on a fluid model, which solves the continuity equation for all species assuming the validity of the drift-diffusion approximation, and Poisson's equation for the electric potential.The fluid model is augmented by a MC code, which generates electron impact rates based on the electric field calculated by the fluid model, and the electron cross section set.The two modules are iterated until convergence is achieved.This approach can achieve a significant speedup compared to kinetic simulation methods (a factor of 30 (serial implementation) to 50 (parallel implementation)) with respect to a serial implementation of PIC/MCC), while providing almost the same accuracy (with an agreement within 5%-10% with respect to PIC/MCC) at atmospheric pressure.As an example, a discharge in He/O 2 mixture was simulated and the results obtained for the spatial density distribution of oxygen atoms were compared to the results of TALIF measurements [60].The experimental and simulation results show good quantitative agreement (within 30%), which verifies the correctness of the simulation technique.Note, that based on the work of Liu et al [39], a fluid model, which uses a two-term BE solver is not enough to describe electron impact processes accurately.
It was shown that the profile of the atomic oxygen density along the jet, where, after an initial buildup phase a saturation sets in, is due to the competition of two effects: (i) the electron impact processes, which, due to the separation of timescales between the electron dynamics and that of other chemical reactions, act as a a temporally averaged (over an RF-cycle) background over timescales of the neutrals, and (ii) chemical reactions, in particular, the ones that destroy atomic oxygen by converting it into ozone.This competition and the diffusion losses together result in a saturation of the atomic oxygen density.The results compared to that of [39] show, that although the electron dynamics, over the timescale of the neutral dynamics can be taken into account as a temporally averaged background source (over one RF-cycle), their proper description is crucial for a quantitative description of the plasma, which is absent from e.g., fluid models, where a two-term 0d BE solver is used.

Figure 1 .
Figure 1.Schematic of the simulation domain for the COST jet.

Figure 2 .
Figure 2. Schematic of the computational method.
shows the spatio-temporal distribution of the x-component of the electric field, E x , at two z-positions: z = −29 mm (a) and z = −5 mm (b), (these values correspond to the left-and rightmost vertical lines in figure 5), along with the two dimensional spatial distribution of the x-component of the temporally averaged electric field (c), and the spatial distribution of the x-component of the electric field at specific time instances, denoted by the vertical dashed/full lines in panel (a, b) (d).We only show the x-component of the electric field, since the y-component of the electric field is much smaller compared to that of E x .

Figure 5 .
Figure 5. Spatial distribution of the time-averaged total positive ion density, n ion,+ (a), total negative ion density, n ion,-(b) and electron density, ne (c) along the plasma channel in units of m −3 , as well as the spatial distribution of the electron density (solid lines) and the total negative ion density (dotted lines) between the electrodes (d) at positions of the vertical lines in panels (b, c) for a helium flow of 0.6 slm.Discharge conditions: ϕrms = 275 V, f = 13.56MHz, p = 10 5 Pa, He + 0.5% O 2 .

Figure 6 .
Figure 6.Spatio-temporal distribution of the x-component of the electric field, Ex at z = −29 mm (a) and z = −5 mm (b), two dimensional spatial distribution of the x-component of the temporally-averaged electric field (c) in units of kV m −1 , and spatial distribution of the x-component of the electric field at times indicated by the vertical dashed/full lines in panels (a, b) (d) in case of a helium flow of 0.6 slm.Discharge conditions: ϕrms = 275 V, f = 13.56MHz, p = 10 5 Pa, He + 0.5% O 2 .

Figure 7 .
Figure 7. Spatial distribution of the electronegativity, β, defined as the ratio between the total negative ion density and the electron density in case of a helium flow of 0.6 slm.Discharge conditions: ϕrms = 275 V, f = 13.56MHz, p = 10 5 Pa, He + 0.5% O 2 .

Figure 8 .
Figure 8. Electron source function as a result of chemical reactions, R e − (x, z, t), averaged over the z-direction (a) and averaged over time (b), electron source function due to ionization processes, S e − (x, z, t), averaged over the z-direction (c) and averaged over time (d) in units of m −3 s −1 in case of a helium flow of 0.6 slm.Discharge conditions: ϕrms = 275 V, f = 13.56MHz, p = 10 5 Pa, He + 0.5% O 2 .

Figure 9 .
Figure 9. Charged species and helium metastable density distributions (averaged over time and the x-direction) along the jet in case of a helium flow of 0.6 slm.Discharge conditions: ϕrms = 275 V, f = 13.56MHz, p = 10 5 Pa, He + 0.5% O 2 .

Figure 10 .
Figure 10.Neutral oxygen species density distributions along the jet (a) and spatial distribution of the dominant chemical reactions involving atomic oxygen along the jet (b) (both quantities are averaged over time and the x-direction) in case of a helium flow of 0.6 slm.In panel (b), reactions creating/destroying atomic oxygen are represented by full/dashed lines.Discharge conditions: ϕrms = 275 V, f = 13.56MHz, p = 10 5 Pa, He+0.5% O 2 .

Figure 11 .
Figure 11.Spatial distribution of the time-averaged total chemical rate (a) and electron impact rate (b) of atomic oxygen creation in case of a helium flow of 0.6 slm.Discharge conditions: ϕrms = 275 V, f = 13.56MHz, p = 10 5 Pa, He + 0.5% O 2 .