Modelling the electric field in reactors yielding cold atmospheric–pressure plasma jets

The behavior of the electric field in Cold Atmospheric–Pressure Plasma jets (CAPP jets) is important in many applications related to fundamental science and engineering, since it provides crucial information related to the characteristics of plasma. To this end, this study is focused on the analytic computation of the electric field in a standard plasma reactor system (in the absence of any space charge), considering the two principal configurations of either one–electrode or two–electrodes around a dielectric tube. The latter is considered of minor contribution to the field calculation that embodies the working gas, being an assumption for the current research. Our analytical technique employs the cylindrical geometry, properly adjusted to the plasma jet system, whereas handy subdomains separate the area of electric activity. Henceforth, we adapt the classical Maxwell’s potential theory for the calculation of the electric field, wherein standard Laplace’s equations are solved, supplemented by the appropriate boundary conditions and the limiting conduct at the exit of the nozzle. The theoretical approach matches the expected physics and captures the corresponding essential features in a fully three–dimensional fashion via the derivation of closed–form expressions for the related electrostatic fields as infinite series expansions of cylindrical harmonic eigenfunctions. The feasibility of our method for both cases of the described experimental setup is eventually demonstrated by efficiently incorporating the necessary numerical implementation of the obtained formulae. The analytical model is benchmarked against reported numerical results, whereas discrepancies are commented and prospective work is discussed.


problem formulation and Mathematical Modeling
The experimental setup of the two plasma reactors that we consider is shown in Fig. 1. Both reactors have an electrode that is electrically biased, to create the necessary electric field for the initiation of the electrical discharges. In the plasma-jet applications the applied field is usually time-varying, where the waveform might be sinusoidal or pulsed. As a first step, in the present work we consider a steady state problem that corresponds to the beginning of the pulse in a pulsed voltage waveform. Moreover, at this stage and for modelling purposes, we neglect the effect of the dielectric material of the dielectric tube by considering only the electric permittivity of the working gas. This is expected to affect the results by predicting a higher peak value of the electrostatic potential close to the electrodes and a faster descent away from them. To assess the effect of this assumption on the results, we examine below the numerical simulation of the electric field in Jánský and Bourdon 15 , where it can be seen that the discrepancy in the axial distribution of the electric field for a dielectric tube with relative permittivity ε = 4 r is relatively small. It is noted that this is not the case in general, as the effect of the dielectric on the electric field is  26 showed that the axial component of the electric field undergoes a small change in the rate of decrement, across the dielectric tube. The radial component, on the other hand is significantly affected, especially in the presence of a charge density, during the streamer propagation. Thus, in the present work, we consider that the assumption of a negligible dielectric, which is necessary for the analytical solution, is a fair conjecture for the electrostatic field.
The two plasma reactors that we consider are based on the experimental configuration that is very common in the plasma-jet community [27][28][29] . The first reactor ( Fig. 1(a)) consists of an annular electrode placed around a dielectric tube. The solution area consists of two regions, as indicated in Fig. 1(c), which include the cylindrical area inside the electrode and extend from the half-length of the electrode on the left to infinite distance on the right. The second reactor ( Fig. 1(b)) consists of two annular electrodes placed around a dielectric tube, the first of which is grounded and the second is electrically biased. The solution domain in this case consists of four regions, as indicated in Fig. 1(d), which include the cylindrical area inside the electrodes and extend from the half-length of the electrode on the left to infinite distance on the right. Aiming to mathematically model the electrostatic activity within the plasma reactor system (Fig. 1), we discriminate it into properly adjusted subsectors as demonstrated in Fig. 2 (this scheme incorporates both configurations, with the first one to be a subcase of the second). In this concept, air is treated as a single species, so a binary gas-air fluid system is considered. To this end, we implement to our work the circular cylindrical coordinate system 23 , which in terms of the variables ρ ∈ +∞ [0, ) , ϕ π ∈ [0, 2 ) and ∈ −∞ +∞ z ( , ), it is defined via either the Cartesian or the cylindrical basis aŝˆ∑  www.nature.com/scientificreports www.nature.com/scientificreports/ m hd 3 4  in which, by definition of the gradient ∇ and the Laplacian Δ differential operatorŝˆρ for any ρ ∈ +∞ [0, ), ϕ π ∈ [0, 2 ) and ∈ −∞ +∞ z ( , ), Maxwell's equations 24 reduce to j j j j for ∈ Ω r j , where V j and E j for = j m 1, 2, 3, 4, are the corresponding to each region (3)-(7) electrostatic potentials and electric fields, respectively.
At this point, the system modeling has been hatched and an analytical solution will be presented. Though, it should be underlined that: The characteristics of common potential theory, incorporating harmonic functions as in (9), are well established and standard techniques have been developed when the domain boundaries involve a certain kind of conditions, for instance either Dirichlet or Neumann type. However, the situation becomes complicated when the solution of the given equation is required to satisfy different boundary conditions on disjoint parts of the boundary of the domain where the condition is stated 30 , as it is the case of the current research. Several mathematical methods have been proposed in manipulating such physical problems, indicating that there are no straightforward techniques towards the representation via certain solutions 30 . Hence, here we are obliged to solve several mixed-type boundary value problems, wherein the given solution of the Laplace's Eq. (9) is based upon the classical method of separation of variables and the produced cylindrical harmonic eigenfunctions. Additional terms in this main solution could be incorporated in future work as an optimization process.
On the other hand, due to the rotational symmetry and without loss of generality, we may proceed to our analysis retaining azimuthal independence and consequently to exclude the ϕ-dependence in any potential fieldcalculation. Under this aspect, the electrostatic potential functions admit c n n n n n n 1 1 , 1 c n n n n n n 2 2 , 1 c n n n n n n 3 3 , 1

Electrostatic Potential and Electric Field in the CAPP Jet Systems
Primarily, condition (15) is readily satisfied, since according to the definition of the Bessel's functions 25 . In the sequel, the limiting condition (30), applied on the relative field (14), forces us to calculate Next, symmetry condition (21), reinforced on the electrostatic potential field (10), reads n n n n n hd 1 which reveals that and, consequently, the implicated field (10) is rewritten as ( 1) with ≥ n 1 is a new constant coefficient. Our next task includes the manipulation of boundary conditions (16)- (20), which in view of (35), (11)- (13) and (32), they give c n n n hd n 1, 1 n n n hd n n n n 1 n n n hd n n n n 1  (36) and (38) is that while in order to avoid trivial zero potential fields and , we are obliged to evaluate the corresponding involved separation constants as n n hd n n hd n n hd e n n hd n n hd n n hd and , www.nature.com/scientificreports www.nature.com/scientificreports/ On the other hand, the manipulation of the rest of the transmission in the z-direction conditions (22)-(29), initially require the proper representation of the constant potential V e via a series expansion in terms of Bessel eigensolutions and to this end we adopt the following expansion 25 and www.nature.com/scientificreports www.nature.com/scientificreports/ n n m n n hd leading to the fact that potential fields (46) and (47) coincide to whereas the section of sets (6) and (7) is representing an important result that secures the matching of the fields due to the same Neumann condition for ρ ρ = hd . The transition conditions (49)-(54) involve different arguments of Bessel functions of zeroth order and they are treated in view of the orthogonality relation , C n (2) , D n , ′ C n , ′ D n (3) and D n (4) with ′ ≥ n n , 1 , which is handled with cut-off techniques in order to obtain quadratic systems, whenever the number of sets n′ ≥ 1 of relations (66)-(71) coincides with the number of terms n ≥ 1 of the infinite series, that is ′ ≅ ... , concerning (66)-(71). Notice that the specified value ∈ N  ⁎ is conveniently chosen so as to obtain the order of the desired accuracy of the numerical implementation for the evaluation of the unknowns. Doing so, we have to solve systems of type 6 × 6, 12 × 12, …, 6N × 6N, thus, defining the N × N blocks of matrixes, each one of type 6 × 6 (a total 6N × 6N matrix) of the coefficients of the unknowns as (72) as well as denoting the vector × N 6 1 of the unknown coefficients via and incorporating the vector × N 6 1 of the known constants through , C n (2) , D n (2) , C n , D n (3) and D n (4) for = ... n N 1, 2, , and in terms of (64) and (65) were solved using MATLAB. The number of terms, taken = N 4000 for all cases, was adequate to ensure a smooth transition of the solution between the domains, in the results area, which extends radially from the cylindrical symmetry axis to the internal surface of the dielectric tube. It is noted that closer to the upper boundary of the solution domain, i.e. to the electrode (see Fig. 1(b,d)), more terms would be required in order to obtain a smooth curve for the potential fields.
Once the unknown constant coefficients are evaluated, then the solution of the particular electrostatic potential boundary value problem is given by n n n hd n n 1 1

numerical implementation Results and Discussion
Single electrode reactor (two-region problem/ Fig. 1(a) and 1(c)). The first reactor that is considered is based on the work of Jánský and Bourdon 15 . The electric field is created by a single electrode wrapped around a dielectric tube with internal radius and wall thickness of 1.6 mm. The length of the electrode is 5 mm (in the present calculations the half-length is considered with symmetry boundary condition) and it is biased with a constant applied voltage of 12 kV. In the numerical simulation of Jánský and Bourdon 15 , grounded planes were considered at the edges of the computational domain, in order to ensure that the potential decreases down to zero far from the electrode. Specifically, on the z-axis, grounded planes were placed at 10 cm from left and right edges of the driven ring, i.e. at z = 10 cm with reference to Fig. 3 (see ref. 15 for details on the solution domain).
In the present analysis, the zero potential boundary condition is set at the boundary of the semi-infinite domain, so there is no grounded electrode. It is noted that in both cases, the tube permittivity is approximated close to 1, similar to that of the gas flowing through the tube (i.e. helium), according to the assumption discussed above. The differences in the axial component of the electric field along the symmetry axis are shown in Fig. 3(a), where the grey area corresponds to the region inside the annular electrode. The dashed blue curve was copied from figure No. 1 of Jánský and Bourdon 15 and the solid black line depicts the results obtained from our analytical approach. The present calculations predict a much higher peak value (almost 2.5 times higher) and a much steeper decrease, compared to the numerical results with the grounded electrodes. The electric field in the present analytical solution is zero at an axial distance of 5 mm from the electrode, while in the results with the grounded plates the field extends farther than 30 mm from the electrode, although it is very weak. It must be noted that the main differences are not due to the distance between the grounded planes and the electrode. It can be verified that placing the grounded planes at longer distances from the electrode has an effect that cannot alone account for the differences on the axial electric field distribution. Indeed, to assess the effect of increasing the distance of the grounded plates www.nature.com/scientificreports www.nature.com/scientificreports/ in numerical solutions and the discrepancy due to the dielectric constant of the tube, we solved numerically the Laplace equation, using grounded planes at different distances from the electrode. The solution of the Laplace equation was conducted with the OpenFOAM C++ library 31 , using multi-region support. Figure 3(b) shows the comparison of these numerical solutions with the results of Jánský and Bourdon 15 for two values of the dielectric constant and for grounded planes at different distances. The graphs support the above statement.
Consequently, the discrepancies between the numerical and the analytical results are synergistically attributed to: (i) The conventional choice to solve the mixed-type boundary value problems of our model (Laplace's Eq. (9)) by applying cylindrical harmonic eigenfunctions only. (ii) The fact that the analytical model developed here is not expanded to the radial domain considered by Jánský and Bourdon 15 . (iii) The different modelling of the boundary conditions. Specifically, in the present analytical solution with the zero-potential boundary condition at infinity, the decay of the electric potential in the axial direction is exponential. Otherwise, if the zero-potential boundary condition is placed at a finite distance, though far from the electrode, via grounded planes, the decay of the electric potential is not necessarily purely exponential. (iv) The non-inclusion of the dielectric layer in analytical model. These points determine the prospective work that should be realized as a continuity of the present report.
The second single electrode reactor is based on the experimental configuration reported recently by our group 13 (modified in the sense that the grounded electrode has been removed). It consists of an annular electrode 20 mm in length and internal radius of 3 mm (see Fig. 1(a)), fixed around a dielectric tube of 1 mm wall thickness. The applied voltage is maintained here at +5 kV. This value has been selected according to our previous experiments, i.e. when the reactor of reference 13 was driven by sinusoidal high voltage at 10 kHz, a peak-to-peak amplitude of around (2 × 7) kV was necessary to ignite the electrical discharge, whereas when the same reactor was driven by pulsed high voltage (2 kHz; rising/falling time ~100 ns; duty cycle 10%) 32 , a plateau of about +7 kV was the threshold for discharge ignition. Thus, the amplitude of +5 kV is here selected as a value being not far from a practical operating one but at the same time being low enough to prevent any space charge formation (electrostatic field case). The middle of the 20 mm electrode is positioned at axial coordinate z = 0 mm, where the symmetry conditions are imposed, so the electrode extends up to z = 10 mm. Figure 4 shows the contours of the electrostatic potential ( Fig. 4(a)), the radial electric field (Fig. 4(b)) and the axial electric field (Fig. 4(c)) in the area inside the dielectric tube, following the analytical solution. The radius ρ of the plotted area ranges from 0 to 2 mm (i.e. the inner radius of the dielectric tube) and the axial distance ranges from z = 0 mm, which corresponds to the center of the electrode, to z = 20 mm, where the electrostatic potential and the electric field are close to zero. The results show that the potential drops very steeply and its value becomes negligible after approximately 2 mm from the end of the electrode. Concerning the electric field, it is observed that the axial component is approximately double the radial one, and both components decrease down to zero at a few mm distance from the electrode. Figure 5 depicts the axial distribution of the electrostatic potential ( Fig. 5(a)), the radial electric field (Fig. 5(b)) and the axial electric field (Fig. 5(c)), at two radii, one at the center of the cylindrical geometry ρ = 0 mm (black solid line) and one at the surface of the dielectric ρ = 2 mm (blues dashed line). The grey area corresponds to the region which is inside the annular electrode. The electrostatic potential is almost constant in that region, due to fact that the electrode is relatively long, compared to the one considered in the previous case ( Fig. 3(a)). The electrostatic potential begins to decrease close to the edge of the electrode and, as expected, the drop is steeper close to the surface of the electrode (ρ = 2 mm) than farther from it, i.e. at the symmetry center-line. Beyond a distance of approximately 5 mm from the edge of the electrode, the potential becomes zero. The electric field plots ( Fig. 5(b,c)) reveal that the there is a significant decrease of the electric field close to the symmetry axis. The radial component vanishes on the symmetry axis due to the symmetry boundary condition. The peak value of www.nature.com/scientificreports www.nature.com/scientificreports/ the axial electric field on the symmetry line is approximately half of the corresponding value on the surface of the electrode. Hence, it is inferred that the initiation of the plasma will occur on the surface of the electrode, where the electric field is stronger.
Double electrode reactor (four-region problem/ Fig. 1(b) and 1(d)). The double electrode reactor is based on the experimental configuration of our previous report 13 and it is comprised of two annular electrodes with length 20 mm and internal radius 3 mm, which are positioned with a gap of 35 mm between them (see Fig. 1(b)). The thickness of the dielectric tube is 1 mm and the applied voltage is 5 kV. The middle of the grounded electrode is positioned at z = 0 mm, where the symmetry condition is imposed, and it extends up to z = 10 mm. The electrically biased electrode extends from z = 45 mm to z = 65 mm. In Fig. 6, the electrostatic potential ( Fig. 6(a)), the radial electric field ( Fig. 6(b)) and the axial electric field (Fig. 6(c)) inside the dielectric tube's area depict the analytic solution. The radius ρ of the plotted area ranges from 0 to 2 mm (i.e. the inner radius of the dielectric tube) and the axial distance ranges from z = 0 mm (center of the grounded electrode) to z = 100 mm.
The results show that the potential decreases very rapidly, almost symmetrically, in both sides of the driven electrode. As in the single electrode case, the axial electric field is approximately double the radial one, but both components fade away after a few mm from the edge of the driven electrode. The region inside the grounded electrode, i.e. from z = 0 mm to z = 10 mm has zero potential and electric field, thus the only electrically active area is that around the driven electrode. This fact leads to the conclusion that the existence of a grounded electrode does not affect the electric field significantly, compared to the single electrode case.
This observation can also be made from the plots of Fig. 7, which depict the axial distribution of the electrostatic potential ( Fig. 7(a)), the radial electric field (Fig. 7(b)) and the axial electric field (Fig. 7(c)), at two radii, one at the center of the cylindrical geometry ρ = 0 mm (black solid line) and one at the surface of the dielectric ρ = 2 mm (blues dashed line). The grey line-filled area corresponds to the region which is inside the grounded electrode and the grey area without the line-filling corresponds to the region inside the driven electrode. It can be seen from all the figures that the electrostatic potential and the electric field have non-zero values only close to the driven electrode. The electrostatic potential does not have significant differences between the symmetry line at ρ = 0 mm and the dielectric surface at ρ = 2 mm. The electric fields components, however, are higher on the surface of the dielectric, similarly to the single electrode case. Additionally, the peak values of the electric field components are similar to the ones found in the single electrode case, which supports the notion that in the specific reactor configuration with the specific dimensions, the existence of a grounded electrode does not affect the electric field. www.nature.com/scientificreports www.nature.com/scientificreports/ conclusions In this contribution we have presented an analytical model for calculating the electrostatic potential and the electric field in plasma-jet reactors (CAPP jets). The solution is based on a multiple-region domain within the working gas that enables the calculation of the electric field in different reactor configurations, comprised of one or two electrodes with different sizes and dimensions. The contribution of the dielectric tube to the field variation is considered neglected. This is necessary in order to obtain the analytical solution and, although it may affect the electric field (especially its radial component) it is considered as a fair approximation, based on physical argumentation.
The mathematical modelling of the particular problem is focused on the solution of elliptic-type boundary value problems with either Dirichlet or Neumann boundary conditions, accompanied by the proper behavior of the fields at infinity, which approximates the area at the exit of the jet. Thus, every field potential is written in terms of harmonic cylindrical eigenmodes and three-dimensional expressions are derived as closed-type series formulae. Once the electrostatic potential is obtained within the entire domain of interest, the electric field is  www.nature.com/scientificreports www.nature.com/scientificreports/ readily recovered and all cases are being numerically implemented and discussed, confirming the validity of the analytical approach.
Specifically, the results for the single electrode case have shown the important role of the location of the zeropotential boundary condition, in order to calculate the accurate electric field. The specific boundary condition is difficult to impose at infinity in numerical solutions, so the analytical solution presented in the present study is a useful approach for calculating the electric field in plasma simulations. Moreover, the results for the doubleelectrode reactor showed that for specific electrode configurations, the existence of a grounded electrode may not affect the electric field and consequently the conditions for the initiation of the plasma. Finally, the solution presented in the present paper could be a useful tool in the design of plasma-jet reactors for tailored applications. To this end, the present mathematical model must be enhanced by incorporating the effect of the dielectric tube, expanding the domain to radial distances approaching the zero-potential boundary condition, adding more terms in the present main analytical solution. Work under progress is set towards these directions.