Numerical approaches in simulating Trichel pulse characteristics in point-plane configuration

In this work, a detailed comparison is made of a few different approaches to numerical modeling of non-equilibrium gas discharge plasmas in dry ambient air at atmospheric conditions, leading to Trichel pulse discharge. Simulation models are based on a two-dimensional axisymmetric finite element discretization of point-plane geometry. The negative corona discharge and the hydrodynamic approximation for generic ionic species (electrons, positive and negative ions) are used. The models account for the drift, diffusion, and reactions of the species. They comprise continuity equations coupled to Poisson’s equation for the electric field. Three different formulations were used to specify the ionic reaction rate coefficients. In the first one, the reaction coefficients are approximated by the analytical expressions as a function of the electric field intensity. Two others extract the reaction coefficients from the solution of the Boltzmann equation as a function of the reduced electric field or the electron energy. The effect of gas flow and heating on the pulse characteristics is also investigated. The accuracy of the models has been validated by comparing them with the experimental data.


Introduction
Corona discharge is a locally self-sustained low energy electric discharge that exhibits different modalities at different polarities and voltage levels. Negative corona discharge in a certain voltage range and operating conditions generates a regular train of current pulses, called the Trichel pulses [1]. * Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
With at least a two-electrode setting and given the right conditions, Trichel pulses appear in the discharge current waveform, initiating a transient regime that lasts up to a few microseconds, gradually falling to the quasi-stationary train of regular pulses. The first pulse in the transient regime is always the most energetic, rising up to a few milliamps, and is followed by less intensive transient pulses until they stabilize at about 100 µA. The pulse parameters (magnitude, frequency and others) are affected by the electrode configuration, external electrical circuit [2], voltage level, and gas operating conditions, such as humidity, gas temperature, pressure and chemical compositions [3][4][5][6][7][8].
For the applied voltages in the range of −5 to −30 kV and air gaps in the order of cm, the Trichel pulses in quasi-steady-state regime have a frequency range of a few kilohertz to a few hundred kilohertz, an average current in the range of 5-100 µA, the pulse rise time in the range of nanosecond and a short duration pulse width, usually in the range of 5-200 ns. It has been recently claimed [9] that the pulse front is directly associated with the electron avalanche close to the corona electrode. A different explanation was proposed by Cernák et al [10], who postulated that the positive streamers towards the discharge electrode are propagated first. During the pulse plateau, the discharge is weakly sustained until the next pulse forms. The pulse rise and fall times are highly dependent on the ionization coefficient, and the velocities of the negative and positive ions, as discussed in [7]. The nature, physics and regimes of electric discharges were also discussed in [11][12][13].
An interesting dynamics and multi-physics nature of the pulsed electric discharges opens additional lines of research.
To name a few, Zhang et al [14] studied Trichel pulses' electromagnetic wave signatures and discharge characteristics. Instead of supplying the discharge system with a negative dc voltage, Mizeraczyk et al [15]. used a series of precisely controlled ramp voltage pulses to trigger the system manually and reproduce Trichel pulses. The controlled Trichel pulse generation can find potential applications in electrospray thrusters, or where there is a need for precise manipulation of the production of the charge species, the electric field, and heat and momentum transfer to the objects. The issue of emission of ultrasonic waves was addressed in [16] considering the linearized form of Navier-Stokes' (NS) equations. Finally, the effect of ionic wind and transverse flows were discussed in [17][18][19], respectively.
Self-pulsation can occur in both the negative and positive discharges. It has also been observed in a variety of electrode configurations and geometries, such as in the hollow cathode discharges (HCDs), and the parallel-plate configurations. A brief comparison between HCD, parallel-plate discharge pulses, and regular Trichel pulses generated in curved geometries was reported in [20].
Numerical simulation is inevitable in understanding the underlying nonlinear dynamics, and microphysical processes of pulsed discharges. Existing approaches include fluid (or hydrodynamic) models, the kinetic models based on Boltzmann equation (BE) or Monte Carlo simulations, and the particle-in-cell methods. The kinetic models are the most accurate, but rarely used, due to computational costs. Instead, the hydrodynamic models are adopted most often. Mathematical modeling of pulse discharges date as early as 1966 [21]. However, this field gained momentum only after the papers were published by Morrow [22,23] on the numerical simulation of pulsed discharges using the finite difference method (FDM) with the flux corrected transport (FCT). FCT is a conservative shock-capturing scheme which is proposed for solving certain class of equations with discontinuities that occur mostly in gas dynamics and magnetohydrodynamics. The accuracy of this scheme is maintained by rejecting spurious diffusion and oscillations. More recently, FCT method in combination with the finite element method (FEM) was also used to study the Trichel pulses in air [6].
Higher order linear approximations such as central difference, high order upwind, hybrid difference, and power law give false predictions by introducing spurious oscillations. Low order approximations on the other hand, introduce numerical diffusion. FCT and the total variation diminishing (TVD) schemes work based on combining the low and high order solution approximations [24,25]. FCT always calculates pre-solutions with both low and high order approximations, before merging them into the final solution. TVD does it in a single-step procedure using special functions called flux (or slope) limiters, by a switching mechanism between the low and high order solutions and limiting the solution gradient near the sharp changes or discontinuities. This mechanism maintains the high order approximation in the area of a low solution gradient, but blends both low and high order approximations in areas of sharp changes. Scharfetter-Gummel method [26] is another class of nonlinear, accurate and robust discretization schemes for stiff drift-diffusion transport systems with shocks, and was originally proposed for semiconductor devices. Hagelaar and Pitchford [27] demonstrate an effective utilization of this scheme for solving the electron energy distribution function formulated from the BE. Very recently, a high-order Scharfetter-Gummel scheme has been proposed for gas discharge equations [28]. The coupling between charged particle transport and space charge density cannot be treated explicitly because of the timestep restrictions. However, even if implicit or semi-implicit methods are used, the time step would become restricted by the explicit evaluation of the electron energy source term. In [29], authors presented a seminal work for an implicit treatment of the electron energy source term, which overcomes time-stepping restrictions to speed up the computations by two orders of magnitude.
Another strategy to accelerate the computations adopts a segregated two-step procedure, based on a quasi-steady state approximation for electron concentration [30]; first, a stationary electron continuity equation is solved and then timedependent continuity equations for ions are treated. The simplified model speedup has been three orders of magnitude over a fully coupled model. In [31], authors circumvented the computational burden by introducing an inertial term into the electron transport equation and suppressing the process of pulse formations. In [32], authors used stationary solvers instead of approximate or time-dependent methods, and this made it possible to do a wide range of numerical experiments on the timeaveraged characteristics of corona discharges. Equations for electric discharges not only need longer computational times; they are also stiff and highly nonlinear. This can be mediated to some extent, by logarithmic transformation of the electron energy and density equations [33].
Most numerical studies on the pulsed electric discharges were restricted to two-dimensional Cartesian and twodimensional axisymmetric models. An example of a onedimensional model was reported by Durán-Olivencia et al [34], using a multi-species hydrodynamic approximation. Another one-dimensional approximation of Trichel pulses in air, however using fluid particle-in-cell method, was attempted by Soria-Hoyo et al [35], where they accurately reproduced a sequence of a hundred pulses for a sphere-to-plane geometry and discussed the effects of the reaction coefficients. Study done by Georghiou et al [13] is one of very few, on the three-dimensional hydrodynamic approximation of corona discharge. They present their findings with a detailed review of the numerical approximation of gas discharges in air at standard conditions. The rest of the papers cited here are 2D or 2D axisymmetric, unless stated otherwise.
In the hydrodynamic descriptions of electric discharges, it is widely accepted that considering just a few generic charge species (i.e. electrons, negative and positive ions) could sufficiently capture a general characteristics of the discharge and the process dynamics. In this case, air is usually assumed as a composition of oxygen and nitrogen molecules only. The density of the charge species is predicted by solving a set of charge continuity equations, which account for the drift and diffusion under Poissonian electric field, and some number of most important chemical reactions. The Poisson equation is used to calculate the electric field distribution. The discharge is sustained by including photo-ionization and/or the secondary electron emission from the cathode surface. The last step before solving the governing equations is to extract the swarm parameters: the reaction rates, electron mobility, and diffusion coefficients. These coefficients can be expressed by analytical expressions using one of a few versions available in literature or using the electron energy distribution function obtained by solving the BE. The freeware BOLSIG+ solver [36] developed by Hagelaar and Pitchford [27] is used most often as a tool for determining the swarm parameters. The swarm parameters are calculated as functions of local reduced electric field E/N (local field approximation-LFA), or average electron energy (local energy approximation-LEA). Therefore, in the case of LEA, the electron energy conservation needs to be solved for the electron temperature, in addition to the charge conservation equations and the Poisson equation.
Morrow [22] presented one of the earliest attempts of developing three-species hydrodynamic models of Trichel pulses, but he successfully captured only the first pulse. Napartovich et al [37] developed a three-species quasi-one-dimensional model to describe Trichel pulse formation in air that can reasonably capture the pulse shape and average characteristics of the discharge. In [38], authors studied the impact of reaction coefficients on the pulse characteristics, such as the frequency, average current, and the pulse rise time. The reaction and transport coefficients were defined using analytical expressions taken from [6,13], where the data in [13] were initially calculated by Kang et al [39].
Considering a higher number of ionic species or more complex chemical gas composition, gives more details about physics of the process, but is much more computationally demanding. Air is a composition of O 2 and N 2 , but other gases are also present. Neutral molecules, free radicals and excited states may be considered, too. The number density of the heavy species is sometimes calculated based on the multi-component diffusion transport equation, a simplified form of the Maxwell-Stefan equation, which uses Maxwell-Stefan diffusivities. Often a larger number of species is considered to investigate the role of particular chemical reactions, to study different modes of discharge and transition regimes, or to explore in more detail the impact of some species on the spatio-temporal characteristics of the discharge. For example, Yanallah et al [40] investigated ozone generation by glow discharge in oxygen, for a wire-to-cylinder configuration, supplied with an applied voltage bellow the initiation of pulses. The model included 38 reactions between electrons, ions and neutral species, in ground and excited states, involving 15 species. The reactions were selected based on the work of Soria et al [41], and the reaction rates were extracted using the LFA method. In addition to the transport equations, in [42,43], the electron energy conservation was also considered with the swarm parameters obtained by LEA method. Chen et al [33] examined different stages of the pulses in a coaxial cylindrical gap geometry by adopting a hydrodynamic model with 15 species and 83 reactions: 61 involving electrons and 22 ionion reactions. The electron mobility and diffusion coefficients were obtained using LFA, and the reaction rates were calculated based on LEA. Kokovin et al [44] studied two modes of discharge and their transition from Trichel pulses to stationary glow discharge in air. Their model featured 9 species and 23 reactions, where the reaction rates and transport coefficients obtained using LEA method. Salah et al [45], used two different 2D axi-symmetric models to study Trichel pulses in air for a pin-plate configuration: a three-species model and a more detailed one. Reaction rate and transport coefficients were extracted using analytical expressions in their general three-species model, and using LEA method, in the multispecies model.
The existence of seed electrons is vital in the development of self-sustained discharges. These electrons are generated either by photons reacting with the air particles, leading to the photoionization effect (also known as impact ionization), or as a results of positive ions, or high-energy photons, hitting the discharge electrode. In the discharge models, interaction of photons with the cathode surface that leads to the photoemission is usually neglected and the production of the secondary electrons due to positive ions hitting the surface is implemented by setting proper boundary conditions on the cathode surface. However, modeling the photoionization is computationally much more demanding.
The physical process and mathematical formulations of the photoionization are described in [46][47][48][49][50][51][52]. An efficient implementation of the photoionization, incorporated by Georghiou et al [13], calculates the source term on a coarser secondary grid, and interpolates the values onto the fine primary grid used for the rest of the discharge computations. Even though the photoionization occurs in the discharges of both polarities, its effect can be ignored in the negative discharges [30,53,54]. Therefore, there are only few studies in which the photoionization is considered in the modeling of negative current pulses [7,16,32,33,55]. Pulses can also occur for the positive discharge. Studies on the positive discharges include the photoionization effect, as their main electron source both for the pulsed positive supply voltages [56][57][58], or the positive dc voltages [43].
LEA is less commonly used; most often when a detailed chemical model is needed, or higher energy discharges and discharge transition regimes need to be investigated [33, 42-45, 59, 60]. Swarm parameters are extracted from the solution of BE as discussed earlier. In obtaining these swarm parameters, most studies assume that heavy particles have the same temperature as the ambient air and use a fixed gas temperature in conjunction with the electron temperature.
LFA is a compromise between accuracy and computational expense [13,17,18,31,40,61]. In LFA the reaction parameters are calculated from BE and the swarm parameters are tabulated as a function of the reduced electric field [31]. These parameters can also be approximated by analytical expressions. Georghiou et al [13] summarized a few resources, in which the parameters have been obtained from experimental results or by solving the kinetic model. They used Morrow and Lowke [62] in their own study of hydrodynamic modeling of non-thermal plasmas, and obtained the parameters by solving the stationary BE. Morrow and Lowke [62] modeled positive streamers in air by solving 1D continuity equations using FDM with FCT, and a 2D Poisson's equation. Transport data used in their study were fitted to the analytical expressions in terms of reduced electric field and presented. Yanallah et al [40] used the data reported in [63] for the reaction coefficients in nitrogen-oxygen mixtures.
The speed of the bulk gas flow is normally two orders of magnitude smaller than that of the charge species due to the electric migration. Therefore, the convective flux is usually ignored in the species transport equations. An increased computational cost is another reason for this. If the NS equations are solved, one-way coupling between the flow and the discharge is usually assumed. Chen et al [17,18,61] solved the NS equations with incompressible, weakly compressible, and fully compressible assumptions. The flow models were coupled with the energy equation to account for the gas heating due to the energy deposition from discharge current to neutral gas molecules. In [59], reactive compressible NS equations were solved for the gas heating effect, as was first discussed in [64]. In the study of glow discharges carried out by Yanallah et al [40], the energy conservation and inviscid-incompressible form of the NS were considered. In some scenarios, bulk motion is not negligible and must be considered. Instead of solving the NS equations, Guo et al [19] considered the convective flux with a constant air velocity, to investigate the effect of the transverse airflow. Assuming the ionization and attachment coefficients are influenced by the air humidity and pressure, as it was discussed earlier in [27], the heating, humidity and pressure effects were discussed in [54,55]. Similarly, the effect of air pressure was studied in [65].
The most common platform for the numerical study of the electric discharges is COMSOL Multiphysics software package, a commercial simulation tool based on FEM [66]. As discussed earlier, numerical instabilities are created by discontinuities, or sharp spatio-temporal gradients of rapidly developing discharge fronts. Adding an artificial diffusion to the solution at the expense of a lower local spatial accuracy, can stabilize the solution. This method is used in the streamline upwind Petrov-Galerkin (SUPG) method, which was proposed in [67]. SUPG is similar to the flux limiters used for the finite volume or FDMs [25,[68][69][70][71]. Other stabilization techniques that can be considered, which are independent of the numerical schemes or platforms used, include stabilizing the ionization process [72], changing the shape of the electrodes, for example cathode from conical to hyperboloid, and anode from plane to cylinder, introducing a convective external gas flow, or changing the material. Some authors use either selfdeveloped codes, or platforms other than COMSOL [13,40,41,73].
The present work incorporates three approaches to model the Trichel pulses dynamics and investigate their modeling capabilities. They use different formulations to specify the ion reaction coefficients: expressed by analytical equations (AEs) as a function of electric field, or by solving the BE. In this case the reaction coefficients can be expressed as a function of the local reduced electric field or the local electron energy. The numerical results for the average electric power and pulse frequency rate are compared with the experimental data found in the literature. The numerical models include continuity equations of charge species, Poisson equation for the electric field, and electron energy conservation for the electron temperature. As for the fluid flow, a full set of NS equations (conservation of mass, momentum, and energy) is solved for neutral gas incompressible bulk motion. The effect of gas heating is also examined.

Model description
The negative corona discharge is studied in the configuration shown in figure 1. This is a typical point-plane type discharge consisting of two electrodes: a needle electrode with a hyperboloid profile supplied with a high negative voltage V ap and the plate electrode, which is grounded. The gap between both electrodes is filled with dry air of relative permittivity ϵ r . Air is composed out of N 2 and O 2 molecules with 72:28 ratio. Initially, the atmospheric pressure and temperature 300 K are assumed, but in some models the gas pressure and temperature are affected by the discharge. The discharge, energy and the flow problems are treated as 2D-axisymmetric.

Mathematical model
The discharge model considers three generic charge species: electrons, positive ions and negative ions. A set of three drift-diffusion equations account for the drift, diffusion, and destruction of ionic species. The electric field is produced by the electrodes and the space charge, and is governed by the  Charge species existing or generated in the computational domain either contribute to the ionic reactions and get neutralized or destroyed, or leave the domain. Species deposited on the metal surfaces are removed instantly.
A set of charge continuity equations is solved to predict the concentration of charge species [74]. The source terms of the equations describe chemical reactions, and the flux terms take into account drift and diffusion under electric field. In this study, five chemical reactions are considered: ion-ion and ion-electron recombination, electron attachment, and electron impact ionization. Having both voltage polarities require considering a wider range of reactions and inclusion of cluster ions, heavier and metastable molecules, excited states, and stepwise ionization. For a negative DC discharge, a much smaller set of reactions needs to be considered. The reactions are given in table 1.
The source reaction terms S are defined as: where ⃗ E is the electric vector field, µ is the mobility, α is the ionization coefficient, β is the recombination coefficient, η is the attachment coefficient, and subscripts e, p and n refer to the electrons, and positive and negative ions, respectively.
Three-body attachments with neutral third body: Seed electrons are essential for self-sustained positive and negative discharge processes. When the positive ions collide on the metal surface of the discharge electrode, seed electrons are produced from the surface. This process, called the secondary electron emission, is responsible to sustain the negative discharges, and treated as a boundary condition on the discharge electrode according to equation Γ e = γc p µ p | ⃗ E|, where Γ e is the electron flux and γ is the secondary electron emission coefficient [75]. Photoionization which is most important phenomenon for sustaining the positive discharges, is not considered in the discussed model.
Swarm parameters are treated differently in the models discussed in this paper, except for the mobility of the positive and negative ions, which are equal to 2.43 × 10 −4 m 2 (V · s) −1 and 2.70 × 10 −4 m 2 (V · s) −1 , respectively [13], and the diffusivity D of the positive and negative ions, which are equal to 2.8 × 10 −6 m 2 s −1 and 4.3 × 10 −6 m 2 s −1 [76], respectively.
LFA method: The mobility, ionization and attachment coefficients of the electrons were extracted in terms of the reduced electric field E/N, by solving the two-term BE using Bolsig+ [36], where the electron collision cross sections were obtained from the SIGLO and PHELPS databases [79,80]. Once the electron mobility is obtained, the electron diffusion coefficient is calculated as D e /µ e = 1 V, known as the Einstein relationship [81].
LEA method: The mobility, ionization, attachment, and diffusion coefficients of the electrons were extracted in terms of the electron mean energyε, by solving the two-term BE using Bolsig+ [36], where the electron collision cross sections were obtained from the SIGLO and PHELPS databases [79,80].
A continuity equation governs the electron energy density distribution: (2) where t is the time, c ϵ = c eε = 3 2 k B c e T e is the electron energy density, k B is the Boltzmann constant, T e is the electron temperature, S jl is the electron Joule heating, S el and S inel are heating source terms from electron-neutral elastic and inelastic collisions, respectively. The source terms are defined by: where c i is the number density of neutral heavy species i, η i is the rate coefficient for energy loss from electron-neutral inelastic collisions, and ∆ ϵi is the threshold energy loss of electron-neutral inelastic collisions. m e is the mass of electron, m i is the mass of neutral species i, T e and T are the electron and gas temperatures, respectively, and v e,i is the electron momentum transfer collision frequency to the neutral species i, extracted from Bolsig+ [36].
The Laplacian electric field ⃗ E L also needs to be determined to calculate the discharge current pulse train in air using Sato's formula [82]: where V ap is the applied voltage, Γ s is the flux of charge species s, and A is the discharge area.

Initial and boundary conditions.
Initially, a neutral plasma is assumed with zero concentration for the negative ions, and a Gaussian distribution c 0 for the electrons and positive ions [83]: where N max = 10 6 m −3 , s 0 = 25 µm the radius of initial distribution, r 0 and z 0 representing the coordinates of the discharge tip.
On the remote boundaries, zero normal derivative for electric displacement is assumed. Depending on the direction of the convective flux of species, zero species density distribution, or zero normal derivative is set.
Zero normal density of positive ion and zero negative ion density is assumed on the discharge electrode. On the ground electrode, the voltage and the density of the positive ions are zero. Zero normal derivative density is set for the electrons and negative ions.
As for the electron energy density equation, an energy density flux is imposed on the discharge electrode, zero normal derivative of energy density is set on the ground electrode, and zero energy density is assumed on the remote boundary. The electron energy density flux is: Where v e,th is the electron thermal velocity defined by: The mass, momentum and energy conservation equations of neutral gas molecules are governed by the incompressible NS and energy equations: where ρ a , ν and p are the density, kinematic viscosity and pressure of the air, respectively, ⃗ u is the air flow velocity vector, ⊗ is the outer product operator, ⃗ g is the gravitational vector field, and ⃗ F ehd is the electrohydrodynamic body force. In the energy conservation equation, T is the air temperature, α T = kT ρaCp is the air thermal diffusivity, C p is the air heat capacity, k T is the thermal conductivity of air, and the heating source term S heat is due to the energy transfer from discharge to air.
It is assumed that the kinetic energy is instantly transferred into the thermal energy, therefore, instantaneous gas heating occurs upon collisions between ions and neutrals. The gas heating by energy transfer from electrons to the neutrals is divided into the fast and slow timescale processes. On the shorter timescale, fast energy transfer occurs though elastic collisions and rotational excitation, and though the excitation of N 2 molecules. The energy transfer on the longer timescale happens through the process of vibration-translation V-T collisions, after the relaxation time of vibrational energy of N 2 molecules. The power S heat deposited by the discharge is: (9) where the current densities of the electrons, positive ions and negative ions are denoted by ⃗ j e , ⃗ j p and ⃗ j n respectively, the fractional energy transferred in elastic collisions and rotational excitation of air molecules, electronic excitation and vibrational excitation of N 2 molecules, are η (el−r) , η ex and η v respectively. The three fractions are calculated using Bolsig+ [36]. It was hypothesized in [84], that the efficiency of the energy transfer through the electronic excitation is ζ ex = 30%. The V-T energy transfer is obtained from a phenomenological relationship [85]: where the energy density S vt stored in vibrational states, is relaxed in τ vt = 20 µs [86]. The equation of states for perfect gases, p = (ρ g /m g )k B T, closes equations (8)-(10).

Initial and boundary conditions.
The air flow is initially at rest, therefore the velocity vector is set to zero. The temperature is set to T 0 = 300 K. No-slip Dirichlet boundary conditions are assumed on the solid surfaces of both electrodes, (⃗ u = 0). The remote boundary is an open boundary, imposing zero normal stress, or constant pressure on the boundaries. This condition ensures that the air can enter or leave the domain freely. On the axis of symmetry, zero flux condition is applied for both temperature and flow velocity. For the energy equation, zero normal derivative of temperature is applied on the remote boundary, and constant temperature on both electrodes.

Model parameters
The parameters of the discharge model studied in this paper are listed in table 3. The hyperboloid profile of the discharge electrode is defined by a parametric hyperbolic curve of a form r = 2Rt, and z = Rt 2 + l 1 , where R is the radius of the discharge electrode tip, r and z the spatial coordinates, and t the parameter of the hyperboloid. A range of voltages from −5 to −10 kV was applied to the discharge electrode. The gap length between the electrodes is l 1 and the computational domain is a rectangular box of dimension l 2 × l 3 , truncated far away from the discharge region for a reasonable balance between accuracy and computation time.
The commercial FEM platform COMSOL 5.6 was used for numerical modeling. Based on the three approaches introduced earlier for treating the swarm parameters, three different cases for defining the reaction rate coefficients were considered: AE, LFA and LEA. For all three cases, the Poisson and Laplace equations were solved in the electrostatics module. The Laplacian filed is needed only for the discharge current calculations using (4). Furthermore, the charge transport equations were solved in the transport of diluted species module. A second-order accurate interpolation scheme was used for the solution of the Poisson and Laplace equations. However, a first-order scheme was used for highly advective transport equations, dominated by the first derivatives of the species concentrations. The streamline diffusion method was activated to stabilize the transport equations. This method adds a small numerical diffusion term, negligible compared to the advection term, yet highly effective in stabilizing the solution and facilitating the convergence. For the LEA case, equations (2) and (3) were solved in the transport of diluted species module with a first-order discretization scheme. Additionally, to study the effect of ambient air, equations (8)-(10) were also solved. The laminar flow and heat transfer in fluids modules were used for the NS and energy equations, and the domain ordinary differential equations and differential algebraic equations module (domain ODEs and DEAs) was used for (10). Linear discretization scheme was used for the NS and energy equations, and a thirdorder scheme for the ODEs. After the FEM-based spatial discretization, the implicit backward differentiation formulation was adopted with adaptive time-stepping method, and all the equations were solved in the time domain.
Spatial discretization of the computational domain was done using unstructured triangular grid with about 150 000 elements, leading to about 0.7 million degrees of freedom. Very fine elements were used near the discharge tip, starting from approximately 500 nm, and gradually growing coarser farther from the discharge region. Iteratively refining the mesh elements, different mesh densities were investigated until sufficient accuracy was achieved.

V-I and frequency characteristics
The voltages ranging from −6 to −10 kV were applied for AE and LFA models. Comparison between numerical and experimental results for the pulse frequency and average current, is shown in figure 2. The experimental data has been taken from [87]. LFA switches form the Trichel pulses regime to the glow discharge regime at voltages above −9.5 kV. Both models capture a linear relationship for the DC current and frequency, versus applied voltage, whereas in the experiment, a parabolic relationship holds. Overall, results from AE model are closer to those from the experiment. The agreement for the DC current is much closer compared to that of the frequency, for both numerical models, confirming the earlier findings [87].

Pulse characteristics
A quantitative comparison between the pulse shapes for the AE and LFA at different voltage levels is shown in figures 3 and 4, respectively. Four different voltages have been considered for each model. At each voltage, two pulses were sampled and aligned at t = 0. Numbers on both axes locate the peaks of the pulses. An increased voltage increases both the minimum and the maximum value of the discharge current for   both models, but this has a higher impact on the pulse peaks. The minima and maxima in the AE model increase approximately linearly when the voltage magnitude increases. However, this trend is different for the LFA model, where the minima increase faster and the maxima increase slower, as the voltage magnitude increases.
The discharge stops between two pulses. A plateau between each two bursts of the pulse is flatter at lower voltages for both models. The pulses appear sharper and the plateau appears wider and flatter for LFA than that for AE. Depending on the input voltage, the plateau can be roughly divided into up to three distinct levels, on which the currents stay for some time, before they rise again onto the next pulse peaks. In the LFA model (figure 4), three levels can be identified at lower voltages, and two levels can be identified at higher voltages. In the AE model (figure 3), two levels can be identified at lower voltages, but only one at higher voltages.
A comparison between the single pulse shapes obtained from the AE, LFA, and LEA models at −7 kV is presented in figure 5. The inset plot shows the pulse currents in the logarithmic scale. Numbers on the y axis indicate the peaks of the current pulses.
For the LEA model, the discharge peak is an order of magnitude higher than for two other models, rising to −1119 µA. The current pulse for this model is also significantly sharper, and its rise time is shortest. For the LFA model, the peak of the discharge current is equal to −141.84 µA. The rise time is longer than that of the pulses given by the LEA and shorter than that by the AE models. The pulse width is widest for the AE model and the discharge peak is equal to −65.91 µA.

Species distributions
The spatial distribution of the charge species concentrations, c e , c p , and c n in the discharge region at the instant of the discharge current peak are shown in figure 6. Figures 6(a)-(c) show these concentrations obtained from the AE model. Very close to the discharge tip, the concentration of electrons is negligible. It is exponentially increasing until it reaches the maximum value at about 80 µm from the tip. From this point to the ground electrode, it exponentially decreases with a lower gradient than it was increasing before, and expands in the radial direction. In contrast, the concentration of positive ions is considerable only close to the discharge tip. The length of the ionization region is about 100 µm. The concentration of the negative ions is distributed very similarly to that of the electrons. The concentration of the electrons is one order of magnitude lower than that of the positive and negative ions. The distributions given by the LFA model, shown in figures 5(d)-(f), are qualitatively similar to that given by the AE model. Although they are an order of magnitude higher and distributed more locally. The distributions given by the LEA model are rather different. c e and c n expand over a large area and their maxima locate at about 30 µm from the tip and 80 µm from the axis of symmetry. Their values decay exponentially in both directions. However, c p is localized very near the discharge tip, the ionization region is extended within 50 µm from the discharge tip.

Reaction coefficients
The differences between the results obtained from different discharge models are surprisingly high. The only differences between these models are the reaction coefficients. In order to investigate this effect, the electron mobility, the ionization coefficient and the attachment coefficient for −7 kV and points on the axis of symmetry, at the instant of the discharge current peak are shown in figures 7-9, respectively. The electron mobilities are rather close for the AE and LFA models, starting from 0.022 and saturating to about 0.088. For the LEA model, the mobility rises exponentially and saturates to about 1.8. An order of magnitude higher electron mobility is a reason for very sharp discharge pulses given by the LEA model. Another reason for difference in the pulse characteristics is the difference in the ionization coefficients along the axial direction from the discharge tip. Faster changes correspond to sharper current pulses, and vice versa.
3.6. The effect of ambient gas using LEA All results presented above were obtained assuming air at atmospheric pressure and 300 K. Moreover, the ion convection due to air motion were neglected. While these simplifying   assumptions are commonly accepted, the actual effect of gas heating and motion were checked for the discussed model. This was done for LEA model and just one voltage level equal to −7 kV.
To investigate the effect of air motion, the convection has been added to the species flux: However, the gas flow also generates non-uniform pressure, which affects the reaction rates. While some of these effects are probably negligible, the variable pressure can affect: • The positive and negative ion mobility: where µ 0 is the mobility at normal conditions, µ p is the mobility at a given pressure, p 0 is the normal air pressure (1 atm) and p is the actual pressure. • Concentration of the neutral species where N p is the gas concentration for a given air pressure, k B is the Boltzmann constant, and T is air temperature. This concentration modifies the electron mobility, and ionization and attachment coefficients. The Boltzmann coefficient determines the reduced values of these coefficients and they should be multiplied by the actual concentration of the neutral molecules.
In addition to its effect on the concentration of neutral molecules, the gas temperature also influences the ion diffusion coefficients where D and T are actual diffusion coefficient and temperature at any point, respectively, and D 0 and T 0 are the reference diffusion coefficients and temperature, respectively. The gas temperature also impacts the electron mobility, and ionization and attachment coefficients. This effect is not strong. Moreover, it would be impractical to solve the BE at each point of space. As an approximate solution, the reaction rate coefficients were calculated for two gas temperatures (300 K and 900 K). Then, for a given temperature the reaction coefficients were interpolated between these two extreme values.
The gas temperature also influences the density ρ g , and dynamic viscosity µ g of air. The density is updated by the equation of states for perfect gases. The dynamic viscosity is approximated by the piece-wise quartic polynomial function:  (15) which is constructed by fitting empirical data over the interval of 200 K-1600 K. The Sutherland equation may be used alternatively.
The procedure adopted for LEA including gas motion and heating (LEA * ) algorithm is shown in figure 10. The flowchart demonstrates one and two-way couplings between different parts of the problem. It shows how the discharge model interacts with the flow model and how different parameters are updated. Figures 11(a) and (b) show the air flow pattern for two voltage levels (−7 kV and −9 kV) supplied to the discharge electrode and LEA * approximation. The flow is generated due to the electrohydrodynamic body force close to the discharge tip. Near the axis of symmetry air moves from the discharge tip towards the plate. The flow direction is reversed farther from the axis. At both voltage levels, the air jet forms vortices which can be identified by the streamlines. The flow patterns in both cases are qualitatively very similar. The maximum velocity magnitude for −7 kV and −9 kV are about 4.5 ms −1 and 6.0 ms −1 , respectively.
Gas heating also strongly depends on the electrode voltage (figures 11(c) and (d)). The maximum temperatures for −7 kV and −9 kV, are 550 K and 750 K, respectively. The heat is generated locally very close to the discharge tip, and is convected away, spreading over a larger area in the air gap towards the plate. Insets of the figures cover an area, extending 300 µm in axial and 800 µm in radial directions from the discharge tip. The temperature distributions in both cases are also qualitatively very similar.
All models discussed so far, AE, LFA, LEA and LEA * , yielded surprisingly different Trichel pulse characteristics. The basic pulse parameters for all these models are summarized in table 4. Figure 11. The air flow patterns ((a), (b)) and the temperature distributions ((c), (d)) for two voltage levels −7 kV (left column) and −9 kV (right column) using LEA * model. Spatial dimensions are in millimeters, the flow velocities in ms −1 , and the temperatures in K.

Conclusions
This work presents a detailed comparison between four different numerical approaches, analytical expressions, LFA, LEA, and LEA including gas motion and heating, to model nonequilibrium gas discharge plasmas in dry ambient air at atmospheric conditions. All simulation models were performed for a point-plane geometry with 2d-axisymmetric finite element discretization. The negative corona discharge formulation is considered, leading to Trichel pulse discharge. The discharge models are based on the hydrodynamic approximation for three generic reactive species (electrons, positive and negative ions). The discharge models comprise continuity equations for the reactive species coupled to Poisson's equation for the electric field. The continuity equations account for the drift, diffusion, and reaction terms. AE, LFA, LEA and LEA * formulations were used to specify the ionic reaction rate coefficients. In the AE model, the reaction coefficients are approximated by the analytical expressions as a function of the electric field intensity. In the other models, the reaction coefficients are approximated from the solution of the BE. In the LFA model, the coefficients were extracted in terms of the reduced electric field E/N and in the LEA models, the coefficients were extracted in terms of the electron mean energyε. It was concluded that the DC current, frequency, and the pulse parameters given by the AE model is closer to the experimental data published in [87]. The comparison made between the AE and LFA models for the voltage range from −6 kV to −10 kV suggests that the models capture DC current more closely to the experimental data than the frequency. The frequency is overestimated by a factor of 4-8. The pulse width is smallest for the LEA model, and largest for the AE model. The pulse peak for the LEA model is about an order of magnitude higher than that for the other models. Capturing high current pulses with very short pulse widths given by the LEA model requires a very fine adaptive time-stepping method, and leads to very long computations.
The effect of gas flow and heating on the discharge morphology and the pulse characteristics were also investigated for voltage levels (−7 kV and −9 kV) using the LEA formulation. The convective term was considered in the continuity equations in addition to the drift, diffusion, and reaction terms. The resulting gas flow generates a non-uniform pressure field, which affects the reaction rates. The effect of the gas temperature on the diffusion coefficient was also considered, which in turn affects the reaction rates, as well as the air density and dynamic viscosity. The discharge dynamics and the gas flow and heating patterns produced by the altered reaction coefficients and the added convection term were studied using this two-way coupling between the discharge and the flow models. Maximum temperature and velocity magnitude of (550 K and 4.5 ms −1 ), and (750 K and 6.0 ms −1 ) have been obtained for voltage levels of −7 kV and −9 kV, respectively.

Data availability statement
The data cannot be made publicly available upon publication because they are not available in a format that is sufficiently accessible or reusable by other researchers. The data that support the findings of this study are available upon reasonable request from the authors.