A practical guide to modeling low-current quasi-stationary gas discharges: Eigenvalue, stationary, and time-dependent solvers

The work is concerned with the modeling of low-current quasi-stationary discharges, including the Townsend and corona discharges. The aim is to develop an integrated approach suitable for the computation of the whole range of existence of a quasi-stationary discharge from its inception to a non-stationary transition to another discharge form, such as a transition from the Townsend discharge to a normal glow discharge or the corona-to-streamer transition. This task includes three steps: (i) modeling of the ignition of a self-sustaining discharge, (ii) modeling of the quasi-stationary evolution of the discharge with increasing current, and (iii) the determination of the current range where the quasi-stationary discharge becomes unstable and the non-stationary transition to another discharge form begins. Each of these three steps is considered in some detail with a number of examples, referring mostly to discharges in high-pressure air


I. INTRODUCTION
The physics of many gas discharge systems has been understood reasonably well by now. High-quality data for evaluation of transport and kinetic coefficients and tools performing such evaluation are publicly available, e.g., LXCat 1 and LoKI. 2 Sophisticated numerical models have been developed for the simulation of gas discharge systems. Such models are virtually universally based on time-dependent solvers, which give detailed information on spatiotemporal distributions of plasma parameters and are indispensable for studies of discharges with fast temporal variations, such as high-frequency discharges, pulsed discharges, streamer and spark discharges, etc, e.g., Refs. 3 and 4.
Time-dependent solvers can also be used for the computation of steady-state (or, more precisely, quasi-stationary) gas discharges: an initial state of a discharge is specified and its relaxation over time is followed until a steady state has been attained. An alternative is to use stationary solvers, which solve steady-state equations describing a steady-state discharge by means of an iterative process unrelated to time relaxation. Stationary solvers offer important advantages in simulations of steady-state discharges. In particular, they are not subject to the Courant-Friedrichs-Lewy criterion or analogous limitations on time stepping. This allows one to speed up simulations, with improvement by orders of magnitude in many cases, and is particularly important for modeling discharges with strongly varying length scales, e.g., corona discharges, where a variation of the mesh element size by orders of magnitude is indispensable. Moreover, stationary solvers allow decoupling of physical and numerical stability. Another useful feature of stationary solvers is their ability to compute patterns of complex behavior that can manifest itself in the modeling of gas discharges even in apparently simple quasi-stationary situations. Time-dependent solvers may not provide important information in such cases. Several such examples referring to the modeling of glow discharges and thermionic arc discharges can be found in Ref. 5.
Although most of the popular ready-to-use toolkits for gas discharge simulation employ time-dependent solvers, e.g., nonPDPSIM 6 and Plasma module of commercial software COMSOL Multiphysics®, stationary solvers for gas discharge modeling are provided by Plasimo; 7 COMSOL Multiphysics® provides stationary solvers for general partial differential equations; and although the Plasma module of COMSOL Multiphysics® is intended to work with time-dependent solvers, it can still be used with stationary solvers. 5 This work is concerned with modeling of low-current discharges, including the Townsend and corona discharges, the aim being to develop an integrated approach suitable for the computation of the whole range of existence of a quasi-stationary discharge from its inception to a non-stationary transition to another discharge form, such as the transition from the Townsend discharge to a normal glow discharge or the corona-to-streamer transition. It is convenient to divide the task into three steps: (i) modeling of the ignition of a self-sustaining discharge, (ii) modeling of the quasi-stationary evolution of the discharge with increasing current, and (iii) the determination of the current range where a quasistationary discharge ceases to exist and the above-mentioned nonstationary transition begins. Each of these three steps is considered in some detail with a number of examples, referring mostly to discharges in atmospheric-pressure air.
From the mathematical perspective, the problem of ignition of self-sustaining discharges is an eigenvalue problem. 8 In this work, several methods for its numerical solution are discussed and compared. The method of choice, the so-called resonance method, is physics-based and requires solving a boundary-value problem for steady-state linear partial differential equations, which may be routinely done by means of ready-to-use solvers, including commercial ones. Note that, in addition to being of theoretical interest, understanding the ignition of self-sustaining discharges is important for applications: it is a useful reference point in the investigation of breakdown in high-voltage electrical equipment in low-frequency, e.g., 50Hz, electric fields, where the time of variation of the applied voltage is much longer than the ion drift time. (In this work, the term "breakdown" means a transition from low-to high-current discharge, usually spark or arc discharge, accompanied by a sharp decrease of the discharge voltage, examples being streamer or leader breakdown; the meaning that is usual in electrical engineering.) Moreover, in certain conditions, the breakdown voltage coincides with the voltage of ignition of a self-sustaining discharge, e.g., Ref. 9.
The solution describing the ignition of a self-sustaining quasi-stationary discharge, obtained at the first step, may be conveniently extended to higher currents by means of stationary solvers. The most time-consuming step when using stationary solvers is usually finding a suitable initial approximation, which requires intelligent guesswork. Fortunately, in simulations of low-current self-sustaining discharges, this step can be performed in a routine way using the resonance method. In this work, such integrated approach is discussed in some detail and examples of its application to corona discharges of different configurations and both polarities are shown.
As the current of a quasi-stationary discharge increases, the discharge will lose stability and a non-stationary transition into another discharge form occurs. The loss of stability against small perturbations may be studied by means of solving the eigenvalue problem resulting from linear stability theory. In Ref. 10, this approach was used to study the stability of the Townsend and glow discharges. An alternative approach to investigation of stability is to apply a perturbation to a steady-state solution and to follow the development of the perturbation by means of a time-dependent solver. This approach allows studying stability against both small and finite perturbations. In this work, this approach is illustrated by the example of the stability of a positive point-to-plane corona against perturbations of various amplitudes.
The outline of the paper is as follows. A model of low-current discharges in high-pressure gases is briefly described in Sec. II. Computation of initiation of self-sustaining gas discharges is considered in Sec. III. The eigenvalue problem, which governs the discharge initiation, is formulated and three different methods for its solution are given and compared: direct solution of the eigenvalue problem for the self-sustainment voltage, investigation of stability of the no-discharge solution, and the resonance method. An integrated approach for the calculation of low-current quasi-stationary discharges, which is based on a combination of the resonance method and stationary solvers, is discussed in Sec. IV. Investigation of stability of low-current quasi-stationary discharges by means of eigenvalue and time-dependent solvers is considered in Sec. V. A brief summary is given in Sec. VI. In order not to overload the paper, some material has been combined into three appendixes: Appendix A, where the boundary conditions for drift-diffusion equations are briefly discussed; Appendix B, concerned with plasmachemical processes and transport coefficients of low-current discharges in high-pressure air; and Appendix C, where the effective reduced temperature of a pair of ion species in high electric fields is briefly discussed.

II. A MODEL OF LOW-CURRENT DISCHARGES IN HIGH-PRESSURE GASES
Mathematical ideas discussed in this work apply to both hydrodynamic and kinetic models of gas discharges. For brevity, here the consideration is restricted to the conventional system of hydrodynamic equations describing low-current discharges in highpressure gases, which comprises equations of conservation and transport of charged particles, excited states, and radicals produced in the discharge, and the Poisson equation, Here, subscript α identifies different species produced in the discharge (positive and negative ions, the electrons, excited states, and radicals); n α , J α , D α , μ α , S α , and Z α are, respectively, number density, density of transport flux, diffusion coefficient, mobility, net volume rate of production, and charge number of species α; f is the electrostatic potential; e is the elementary charge; and ε 0 is the permittivity of free space. The source terms S α in the equations of conservation for electrons and positive ions include, in addition to terms describing production of these species in collisional processes, the photoionization term S ph . The transport equation (2) for the charged particles is written under the so-called driftdiffusion approximation. This system of equations also includes other relevant equations, such as equations governing the photoionization and the electron and neutral-gas energy equations. If motion of the neutral gas plays a role, then the convective terms have to be added to the lhs of Eq. (1) and the system of equations includes also the Navier-Stokes equations.
The system of equations is supplemented by usual boundary conditions. In particular, boundary conditions for densities of charged particles on the surfaces of electrodes and dielectric surfaces (in case they border the active zone of the discharge) may be written as described in Appendix A.
Most examples given in this work refer to low-current discharges in high-pressure dry air and have been computed with the use of transport coefficients and a kinetic scheme of plasmachemical processes described in Appendix B. The local-field approximation is employed, so all the transport and kinetic coefficients, including those for the electrons, are assumed to depend on the local reduced electric field E=N and the neutral-gas temperature T. (Here, E ¼ ∇f j j is the electric field strength and N is the number density of the neutral gas.) The kinetic scheme does not consider excited states or radicals and takes into account one species of positive ions, which are designated A þ , the electrons, and three species . The photoionization is evaluated by means of the three-exponential Helmholtz model, 11 with each of the terms satisfying the Helmholtz partial differential equation, Here, A j and λ j are constants (parameters of the three-exponential fit function) given in Ref. 11, p O2 is the partial pressure of molecular oxygen, and I(r) is the product of the probability of ionization of a molecule at photon absorption and the local photon production rate. The latter is assumed to be proportional to S i (r), the rate of ionization of neutral molecules by electron impact, and I r ð Þ is written as 12 where p is the neutral-gas pressure and p q = p þ p q À Á is a quenching factor that accounts for the non-radiative de-excitation of radiating states of nitrogen molecules due to collisions with other molecules. The quenching pressure p q is set equal to 30 Torr. 13,14 The examples reported in this work are limited to low discharge currents, where the discharge-induced heating and motion of the neutral gas are negligible, and refer to T ¼ 300 K.
The boundary conditions for electron and ion densities are those specified at the end of Appendix A. The electron emission flux is related to the flux of incident positive ions by the effective secondary emission coefficient γ, which is assumed to characterize all mechanisms of secondary electron emission (due to ion, photon, and excited species bombardment). 3 The rate of photoionization is set to zero at all solid surfaces, S (j) ph ¼ 0 (j ¼ 1, 2, 3), similarly to Refs. 15 and 16.
The numerical modeling reported in this work was performed with the use of commercial software COMSOL Multiphysics®. The following interfaces were used: transport of diluted species or TDS [equations of conservation and transport of charged species, Eqs. (1) and (2)], electrostatics [the Poisson equation, Eq. (3)], and coefficient form partial differential equations [the Helmholtz equations, Eq. (5)]. The streamline (Galerkin-Petrov) and crosswind diffusion stabilizations, which are default options of the TDS interface, were kept activated in all cases except where otherwise indicated. It should be stressed that both stabilizations are consistent, i.e., the corresponding terms vanish when the iterations have converged.

III. COMPUTING INITIATION OF SELF-SUSTAINING GAS DISCHARGES
A. The eigenvalue problem The condition of initiation of a self-sustaining gas discharge, where the discharge voltage is just sufficient for the electron impact ionization to compensate losses of the charged particles, is well known for wide parallel-plate electrodes, where the applied electric field is uniform and diffusion of the charged particles is of minor importance. The condition reads where α is the Townsend ionization coefficient, d is the discharge gap width, and γ is the effective coefficient of secondary electron emission from the cathode. An approximate relation for other discharge configurations is obtained by replacing αd by the so-called ionization integral, which is the line integral of α evaluated along the electric field line that ensures the biggest value of this integral. Such approach is theoretically incomplete and cannot be generalized, without invoking arbitrary assumptions, to account for potentially important effects such as the photoionization or the presence of dielectric surfaces. Nevertheless, it remains the main tool used in industrial applications, e.g., recent work. 17 However, the problem of ignition of self-sustaining discharges may be solved by accurate mathematical means and this solution is relatively simple. 8 At the ignition, the charged particle densities are very low and the applied electric field is not perturbed by plasma space charge, so the rhs of Eq. (3) may be dropped and this equation assumes the form of the Laplace equation, Let us first assume that no dielectric surfaces border the active zone of the discharge, then perturbation of the applied electric field by surface charges deposited on the dielectric need not be considered as well. (The case of a discharge along a dielectric surface will be considered in Sec. IV C.) Therefore, the distribution of the electric field in the gap, for a given gap geometry, is governed solely by the applied voltage U and may be found by means of standard electrostatic simulations disregarding the presence of charged particles in the gap. Of course, the electric field distribution is linear with respect to U and it is sufficient to perform the electrostatic simulation only once and then to scale the computed electric field distribution to each required value of the applied voltage U. For definiteness, it is assumed that the applied voltage is defined as the potential of the anode with respect to the cathode, thus U . 0. At the ignition, there is no need to consider nonlinear processes, such as reactive collisions involving two or more particles produced in the discharge (ions, electrons, excited states, or radicals), which includes the stepwise ionization and the recombination of charged particles. Thus, there is no need to consider excited states and radicals and it is sufficient to consider the equations governing the charged-particle distributions in the gap, Eqs. (1) and (2) with α referring to the charged species.
According to the conventional definition, the ignition voltage is the value of the applied voltage U such that the applied electric field produces an ionization intensity just sufficient to compensate for the losses of charged particles, so a steady-state balance of the charged particles can be reached. Therefore, the system of equations should be considered in the stationary form, i.e., without derivatives with respect to time in Eq. (1). The resulting equations governing distributions of the ion and electron densities read where subscript α refers to the charged species (the ions and the electrons).
Since there is no need to consider nonlinear processes at the ignition, the source terms S α are linear with respect to the chargedparticle densities. Note that this implies that the dependence of the photoionization rate S ph on the electron density n e is linear. This is indeed the case as exemplified by Eqs. (4)-(6): I r ð Þ the rate of ionization of neutral molecules by electron impact is proportional to the local electron density n e r ð Þ; therefore, the dependence of S ph r ð Þ on n e r ð Þ, while being non-local, is linear. Note that this is the case also in models where the photoionization rate is computed by evaluation of an integral (e.g., Refs. 6,14,18,and 19) rather than by solving partial differential equations.
Equations (9) and (10) are supplemented by Eqs. (4)-(6) or similar and boundary conditions Eq. (A6) with ξ a ¼ ξ e ¼ 1=2 or similar. The obtained boundary-value problem, which is considered for a given distribution of the electric field in the gap, governs distributions of the ion and electron densities at inception and is linear (with respect to these densities).
We consider conditions where there is no external ionizing radiation and the ionization mechanisms are direct ionization by impact of electrons, accelerated by the applied electric field, and photoionization by photons produced in the discharge; there is no thermionic, thermo-field, and field electron emission from the cathode, as well as no electron photoemission caused by external radiation. Then, the above-described linear boundary-value problem governing distributions of the ion and electron densities in the gap is homogeneous. Since the problem is considered for a given distribution of the electric field in the gap, the applied voltage U is a control parameter of the problem. The aim is to find a value U ¼ U 0 such that the problem describes the inception of a self-sustaining discharge, i.e., a low-current steady-state discharge.
From the mathematical perspective, this is a boundary-value problem for a system of partial differential equations (or integrodifferential equations, if the photoionization rate is computed by evaluation of an integral) and this problem is linear (with respect to the charged-particle densities) and homogeneous (no external ionization terms). For all values of the applied voltage U, the problem admits a trivial solution: the ion and electron densities are zero at all points in the gap, which correspond to a situation where no discharge has been ignited. The task is to find a value of U such that the problem admits, in addition to the trivial solution, also a nontrivial one. In mathematical terms, this is an eigenvalue problem and U is the eigenparameter.
Thus, the ignition of a self-sustaining discharge is described by an eigenvalue boundary-value problem for a system of partial differential or integrodifferential equations governing distributions of the ion and electron densities. The physical sense of the problem is that the applied voltage should be such that direct ionization by the impact of electrons accelerated by the applied electric field and photoionization by photons generated in the discharge are just sufficient to compensate for the losses of charged particles. In agreement with the above, in the special case of wide parallel-plate electrodes, the formulated eigenvalue problem may be reduced to one dimension and leads to the well-known self-sustainment condition (7), provided that similar simplifications are applied: no diffusion, one ion species, and no photoionization.
There is one more special case where the applied electric field can be considered uniform and the formulated problem may be reduced to one dimension: the problem of self-sustainment field in the cross section of a plane or cylindrical positive column of a long low-current discharge. The formulated problem takes a form similar to the well-known eigenvalue problem for an ordinary differential equation describing the ambipolar diffusion of charged particles in the cross section of a positive column under the assumption of quasineutrality, e.g., Ref. 4. The solution to the latter problem is well known: the cosine for a plane column and the zero-order Bessel function for a cylindrical column, respectively, Eqs. ( 4. At the discharge ignition, the transversal electric field is zero, the diffusion is free and not ambipolar, and there is no quasi-neutrality. Therefore, the above-cited solution 4 is valid for the distribution of the election density n e at the discharge ignition provided that the ambipolar diffusion coefficient is replaced by the electron diffusion coefficients D e . The ion density exceeds n e by the factor D e =D i , where D i is the diffusion coefficient of the ions.
In the general case, the formulated eigenvalue problem for the ignition of a self-sustaining discharge is multidimensional and its numerical solution is not quite trivial. Three numerical methods for solving this problem are discussed in Secs. III B-III E.
B. Direct solution of the eigenvalue problem for the self-sustainment voltage The most direct approach to solving the eigenvalue problem, formulated in Sec. III A and governing the ignition of a self-sustaining discharge, is to apply an appropriate standard eigenvalue solver directly to the eigenvalue problem as it was formulated, considering U as the eigenparameter. The solver will return a set of eigenvalues (spectrum) and eigenfunctions associated with each of the eigenvalues.
Each of the eigenfunctions is a vector with each of the components describing a distribution of the number density of one of the ion species or the electrons. Some of the eigenfunctions have components all of which are non-negative (or nonpositive, which is equivalent since eigenfunctions are defined to the accuracy of constant factors) at all points of the computation domain, i.e., the discharge gap. Eigenvalues associated with such eigenfunctions may be physically meaningful. Other eigenfunctions have one or more components that change sign in the discharge gap. Since the components of the eigenfunctions describe distributions of chargedparticle species, the eigenvalues, with which such eigenfunctions are associated, have no physical meaning. There can exist also eigenvalues associated with eigenfunctions that have some of their components positive and others negative. Such eigenvalues have no physical meaning as well.
In addition, most eigenvalues are complex (and appear in pairs of complex conjugate values), and others are real. Since the eigenparameter is the applied voltage U, complex eigenvalues have no physical meaning, while real ones may be physically meaningful.
One of the eigenvalues must simultaneously be real and associated with an eigenfunction with all the components being nonnegative at all points of the computation domain. This eigenvalue will give the voltage of ignition of a self-sustaining discharge (for brevity, the term "self-sustainment voltage" will be also used).
This approach requires the use of an eigenvalue solver for boundary-value problems for systems of partial differential or integrodifferential equations. For example, an eigenvalue solver for boundary-value problems for systems of partial differential equations provided by commercial software COMSOL Multiphysics® is used in this work. A more laborious option is to first manually discretize the partial differential or integrodifferential equations, thereby converting the problem into an eigenvalue problem for a system of linear algebraic equations, for which solvers are more easily accessible.

C. Investigation of stability of no-discharge solution
The second approach is indirect. Let us return to the full problem considered in Sec. II, which includes, in particular, the Poisson equation. This problem admits two steady-state solutions. First, there is a solution with the ion and electron densities equal to zero at all points in the gap and with the electrostatic potential governed by the Laplace equation. This solution describes a situation where the voltage is applied to the electrodes; however, no discharge has been ignited. Of course, the discharge current I is zero for this solution. It should be stressed that although the ion and electron densities equal to zero do not represent a solution to the eigenvalue problem formulated in Sec. III A (which, by definition, must be non-trivial), the zero densities do represent, jointly with the Laplacian potential distribution, a formal solution of the full problem for any value of the applied voltage U. Second, there is a solution with non-zero charged-particle densities and the electrostatic potential governed by the Poisson equation. This solution describes the self-sustaining discharge and exists for all values of the discharge current I, starting from infinitesimal values at the inception.
A graphic illustration is shown in Fig. 1. The current-voltage characteristic (CVC), described by the no-discharge solution, i.e., the one with the zero ion and electron densities, coincides with the positive part of the ordinate axis The solution with non-zero densities, describing the self-sustaining discharge, branches off from the former solution at the ignition point In mathematical terms, this is a bifurcation point. Points of bifurcation of steady-state solutions are neutrally stable. Therefore, the spectrum of stability of the no-discharge solution at U ¼ U 0 must include the zero eigenvalue. It is known from experiment that the no-discharge solution is stable for U , U 0 and unstable for U . U 0 . Note that the no-discharge solution, of course, cannot be realized experimentally and cannot be calculated with a time-dependent solver in the range U . U 0 , where it is unstable. However, it should be emphasized once again that the no-discharge solution satisfies the full problem for any applied voltage U and can be trivially computed by means of a stationary solver for any U including U . U 0 .
The stability of the no-discharge solution is governed by a problem coinciding with the eigenvalue problem formulated in Sec. III A except that the time derivative terms @n α =@t in the charged-particle conservation equations (9), which were omitted in the formulation of the eigenvalue problem, are taken into account. The time dependence of all charged-particle densities n α is assumed to be exponential and described by the same factor e λt , where λ is a parameter to be determined (the perturbation growth FIG. 1. Schematic of current-voltage characteristics at inception. The solution with non-zero ion and electron densities, describing the self-sustaining discharge, branches off from the solution with the zero densities, describing the no-discharge situation, at the bifurcation point The obtained time-independent boundary-value problem for a system of linear homogeneous partial differential (or integrodifferential) equations for n α1 r ð Þ represents an eigenvalue problem with λ being the eigenparameter. Its spectrum, i.e., the set of λ values, computed for a given U value governs the stability of the solution with the zero ion and electron densities, describing the no-discharge situation for this U. It should be stressed that this eigenvalue problem does not coincide with the eigenvalue problem formulated in Secs. III A; in particular, the eigenparameter is different: λ instead of U.
Most eigenvalues are complex. Such eigenvalues describe perturbation modes that grow or decay in time in an oscillatory manner. Other eigenvalues are real and describe perturbation modes that grow or decay monotonically in time.
The spectrum is computed, by means of an appropriate eigenvalue solver, for a number of values of the discharge voltage U, starting from a value that is surely lower than the self-sustainment voltage U 0 . Then, U is increased in small increments.
At low applied voltages U, all the eigenvalues must be positioned in the left half-plane; in other words, real parts of all the eigenvalues must be negative, which means that all perturbation modes decay in time. As U increases, one of the eigenvalues must cross the origin while the others still remain in the left half-plane. All the components of the eigenfunction associated with the zero eigenvalue must be non-negative. This crossing represents the point of neutral stability or, in other words, the bifurcation point being sought, and the corresponding U value is the self-sustainment voltage U 0 .
The spectrum of stability can be computed using the same two options that were described in the last paragraph of Sec. III B. Note that in Ref. 20, the second option was employed in order to determine the self-sustainment voltage of the parallel-plate discharge.
Note that this approach is in essence an alternative way of solving the eigenvalue boundary-value problem formulated in Sec. III A; in contrast to the direct approach, which is described in Sec. III B and in principle gives the whole spectrum (set of U values) of the eigenvalue boundary-value problem formulated in Sec. III A, this approach is designed to find only one U value, namely, the one that corresponds to the self-sustainment voltage.

D. The resonance method
The third approach is physics-based and may be described as follows. The stationary boundary-value problem governing distributions of the ion and electron densities at inception, formulated in Sec. III A, is considered and the equations of conservation of the electrons and one of the positive ion species are supplemented by a term describing an external ionization source. In other words, Eq. (9) for the electrons (α ¼ e) and one of the positive ion species (α ¼ i) are replaced by the following equations: where W is the external ionization term. This term does not depend on the particle densities, in particular, on the electron density, nor on the applied voltage and is specified more or less arbitrarily. For example, it may be set constant in space in one-dimensional (1D) problems and a Gaussian function with a maximum on the discharge axis is a natural choice in axially symmetric problems. The obtained stationary boundary-value problem describes the non-self-sustained discharge, generated in the same plasmaproducing gas and the same electrode configuration. The problem is solved for different values of the applied voltage U. On physical grounds, one can expect that a kind of resonance will appear when U becomes equal to the self-sustainment voltage.
Note that the problem is linear and the external ionization source W is the only inhomogenous term. Therefore, the absolute values of W are irrelevant as far as the task is restricted to the calculation of the self-sustainment voltage: the scaling of W affects the scaling of the number densities of the charged particles and does not affect the self-sustainment voltage. On the other hand, some care in the choice of the scaling of W is useful if this method is used as a part of the integrated approach for modeling quasistationary low-current discharges in the whole range of its existence starting from the inception; and a comment on this point is provided in Sec. IV A.
The stationary boundary-value problem describing the non-self-sustained discharge, being linear (and inhomogenous), may be routinely solved by means of ready-to-use solvers; hence no need for manual discretization. For example, a linear solver for boundary-value problems for systems of partial differential equations provided by COMSOL Multiphysics® is used in the modeling reported in this work. Note that COMSOL automatically activates the nonlinear solver option in the case where the default streamline and/or crosswind diffusion stabilization is kept activated in the TDS interface since the stabilization terms appearing in the species conservation equations are nonlinear. However, the number of iterations required is very small, typically no more than two, so the procedure remains the same whether stabilization is activated or not.
Similarly to the previous one, this method is in essence another way to solve the eigenvalue boundary-value problem for the system of partial differential or integrodifferential equations, formulated in Sec. III A, which finds not the whole spectrum but rather only the self-sustainment voltage. The method was introduced and termed the resonance method in the preceding work. 8 Of course, the very idea of finding a real eigenvalue by varying the eigenparameter and searching for this or that kind of resonance is quite obvious. For example, a change of sign of the determinant of the system of algebraic equations obtained by a finite-difference discretization of a boundary-value problem was used in Refs. 21-23 as an indication of bifurcation of regimes of current transfer to, respectively, cold cathodes of glow discharges and hot thermionically emitting cathodes of arc discharges.

E. Comparing different methods
In order to illustrate the above three methods, it is natural to apply them to the most basic 1D model of a parallel-plate discharge. Let us consider, as an example, a microdischarge in xenon at the pressure of 30 Torr. The interelectrode gap is d ¼ 0:5 mm.
Calculations have been performed in the framework of a simple model accounting for a single positive ion species and the electrons.
Transport and kinetic coefficients are the same as in Ref. 24. The secondary electron emission coefficient γ was set equal to 0:03. The simulations have been performed both with and without the streamline and crosswind diffusion stabilization on two different uniform numerical meshes, one with 170 and the other with 3000 mesh elements.
Let us first consider results obtained by means of the resonance method. The CVC of the non-self-sustained discharge, j U ð Þ (here, j is the current density), was computed numerically in the framework of this method for an external ionization source term that was uniform in space and is depicted in Fig. 2 by solid lines. There is a singularity at U ¼ 175:05 V, after which the solution loses its physical significance: the ion and electron densities turn negative, as well as the discharge current. This is the expected resonance-like phenomenon, and the U value at which it happens is the self-sustainment voltage.
An approximate analytical description of the CVC of the parallel-plate non-self-sustained discharge may be obtained by neglecting the diffusion of the charged particles and is given by the following formula: The CVC described by this equation duly conforms to the CVC given by the numerical calculations, shown in Fig. 2. In particular, there is a singularity, after which j turns negative. The value of U at which the singularity occurs is obtained by setting the denominator of Eq. (13) equal to zero, which leads to the conventional condition of self-sustainment Eq. (7). The value of U determined by this condition is approximately 177:20 V, which duly conforms to the above-mentioned value of the self-sustainment voltage of 175:05 V, computed numerically. Note that the 1:2% difference is due to the neglect of diffusion in Eq. (13). Thus, the expected resonance-like behavior of characteristics of the non-self-sustained discharge, which should occur at an applied voltage coinciding with the self-sustainment voltage, is indeed present and may be readily identified as a singularity in the current-voltage characteristic I U ð Þ, accompanied by a change of sign of current. Note that 170 mesh elements were sufficient to obtain the self-sustainment voltage with five significant digits, and hence the method is not too sensitive with respect to the mesh. Moreover, the computed self-sustainment voltage was the same with five significant digits regardless of whether stabilization was activated or not.
The real spectrum computed in the framework of the direct solution of the eigenvalue problem for the selfsustainment voltage (the method described in Sec. III B) comprised two eigenvalues, U ¼ 15:633 V and U ¼ 175:05 V. The ion and electron densities described by the eigenfunctions associated with both of the eigenvalues were non-negative at all points of the gap and both eigenvalues up to five significant digits were not affected by either the change in the number of mesh elements or the activation of stabilization. (The calculations also returned at least one more real eigenvalue; however, that eigenvalue was constantly increasing with refinements of the numerical grid.) The bigger one of the two eigenvalues conforms very well to the above value of the self-sustainment voltage computed by means of the resonance method. Note that the product αd, evaluated using the bigger eigenvalue, equals 3:47, while ln 1=γ þ 1 ð Þ¼3: 54 for γ ¼ 0:03. Thus, the bigger eigenvalue approximately satisfies the conventional self-sustainment condition Eq. (7) as it should. (The small discrepancy is due to the neglect of diffusion in this condition.) On the other hand, αd ¼ 1:36 Â 10 À3 when evaluated using the smaller eigenvalue. Thus, the eigenvalue U ¼ 15:633 V is way too low to fulfill the self-sustainment condition. Therefore, this eigenvalue can only be an artifact. Surprisingly, no warning signs appeared in the course of the numerical calculations; it should be stressed, in particular, that this eigenvalue did not depend on the mesh. Additional comments on this point are provided later in this section.
The spectrum of stability of the solution with the zero ion and electron densities (the no-discharge solution), computed in the framework of the method described in Sec. III C, includes a real eigenvalue associated with an eigenfunction with both components non-negative at all points of the gap. There is also a large number of complex conjugate eigenvalues, all of which were associated with eigenfunctions whose components change sign in the computation domain. The real eigenvalue is plotted in Fig. 2 by the dotted line. As expected, it is negative at low voltages and then changes sign. The latter occurs at U ¼ 175:05 V and the corresponding point is marked in Fig. 2 by the circle; this is the desired point of neutral stability or, in other words, the bifurcation point, and the corresponding U value is the self-sustainment voltage U 0 . This value was not affected by the change in the number of mesh elements; on the other hand, the method appeared not to work with the stabilization activated.
Results obtained in this work for the spectrum of stability of the no-discharge solution conform to results of Ref. 20. Only one real eigenvalue was found in Ref. 20, which is plotted in Fig. 5 of Ref. 20. The components of the eigenfunctions associated with the real eigenvalue, which are plotted in Fig. 6  Thus, all three methods described in Secs. III B-III D gave the same value for the self-sustainment voltage, U ¼ 175:05 V, as they should, with the exception of the extraneous root given by the direct solution of the eigenvalue problem.
As another example, consider the results of numerical simulation of the ignition of a positive wire-to-plane corona discharge under conditions of experiment: 25 air, 1 atm, wire diameter 89 μm, and wire-to-plane distance 10 mm. The self-sustainment voltage was computed by means of the direct numerical solution of the eigenvalue problem and by means of the resonance method, on two different numerical meshes, one with 1 65 240 and the other with 5 34 184 degrees of freedom, with and without stabilization. The direct numerical solution of the eigenvalue problem gave the self-sustainment voltage of 4278:5 and 4278:3 V on the two meshes with stabilization and 4284:4 and 4285:3 V without stabilization. No extraneous roots have been detected, in contrast to the above example of parallel-plate discharge. In the framework of the resonance method, a singularity in the current-voltage characteristic of the non-self-sustained discharge with a change in the sign of the current, similar to that shown by the solid lines in Fig. 2, was found. The self-sustainment voltage of 4278:0 and 4278:2 V with stabilization and 4282:2 and 4280:7 V without stabilization was obtained. All the above values are close to each other and to the experimental ignition voltage, which is approximately 4:3 kV as seen in Fig. 7 of Ref. 25. Let us now compare the three methods from the methodological point of view. The methods are mathematically equivalent as far as computation of the self-sustainment voltage is concerned, so the question is which method is easier to implement and use and which method has wider applicability.
The method based on the direct numerical solution to the eigenvalue problem requires computing and analyzing the spectrum only once, while the stability method requires doing this several times (for different U values). Thus, the stability method is more laborious. On the other hand, the method of direct solution of the eigenvalue problem must be used with caution in view of the potential appearance of extraneous roots, as in the example of parallel-plate discharge. This is an unfortunate drawback of a theoretically simple method. Perhaps, this difficulty could be overcome by further fine-tuning the eigenvalue solver used or by using a different solver. On the other hand, this difficulty is a manifestation of a general trend: finding the spectrum by means of an eigenvalue solver, required in the framework of the first and second methods, is a nonlinear task and, as such, requires care in certain cases; in particular, distinguishing between physical and artificial eigenvalues is not always easy.
In contrast, the stationary boundary-value problem, describing the non-self-sustained discharge in the framework of the third (resonance) method, is linear and its solution is straightforward. The self-sustainment voltage U 0 in the resonance method has to be identified from computations of the non-self-sustained discharge for different U values rather than directly as in the first method. It should be stressed, however, that the pattern of the CVCs of the non-self-sustained discharge, which is shown by the solid lines in Fig. 2 and comprises a singularity at U ¼ U 0 with the subsequent change of sign of current, has been observed in every case where the resonance method has been applied. This includes also simulations not shown in this paper. In all the cases, the resonance and, thus, the U 0 value were readily identifiable and the procedure was straightforward, fast, and reliable. Typical calculations for 2D geometries take about 10 min on a desktop computer.
An important advantage of the resonance method is its physical transparency. In particular, it allows the use of this method not only for simulations of ignition of a self-sustaining discharge but also as a module of a more general code for the modeling of quasi-stationary (self-sustained) discharges. This is described in Sec. IV.
In summary, the resonance method is the method of choice for the calculation of the ignition of self-sustaining discharges at least for now. The method has given useful results in a wide range of conditions and complex geometries, e.g., it was used for investigation of discharge ignition in air between concentric-sphere and concentric-cylinder electrodes with microprotrusions of different shapes on the surface of the inner electrode in a wide range of air pressures. 26 Investigation of discharge ignition along a dielectric surface is considered in Sec. IV C.

IV. INTEGRATED APPROACH FOR THE CALCULATION OF LOW-CURRENT QUASI-STATIONARY DISCHARGES
As discussed in Sec. I, stationary solvers represent an appropriate tool for the modeling of quasi-stationary discharges and offer important advantages over time-dependent solvers. The most time-consuming step when using stationary solvers is usually finding a suitable initial approximation, which requires intelligent guesswork. Fortunately, in simulations of low-current selfsustaining discharges, this step can be performed in a routine way using the resonance method. The resonance method provides a very accurate first approximation for a steady-state solution at a very low current, describing the onset of the self-sustaining discharge. Once this solution has been computed, it can be used as a starting point for the simulations with the discharge current being gradually increased. The solution for each current serves as an initial approximation for simulations with the next current value.
In this way, one obtains an integrated approach for modeling quasi-stationary low-current discharges in the whole range of its existence starting from the inception. This integrated approach is described in this section, and examples of its application are given.

A. Combining the resonance method and stationary solvers
The procedure is illustrated by the flow chart shown in Fig. 3 and may be described as follows. At the first step, the non-self-sustained discharge is computed as described in Sec. III D.
Remind that, as discussed in Sec. III A, the plasma space charge is neglected at this step, i.e., the rhs of the Poisson equation (3) is dropped or, in other words, the Poisson equation (3) is replaced by the Laplace equation (8). Nonlinear processes, such as reactive collisions involving two or more particles produced in the discharge, are neglected as well. The discharge voltage U is increased in small increments ΔU until the discharge current becomes negative. Let us designate by U 1 the highest voltage for which the current is still positive and by I 1 the corresponding current.
The second step consists in recalculation of the solution obtained at the first step for U ¼ U 1 but with the discharge current being the control parameter instead of the voltage. This amounts to solving the same problem as at the first step but considering the discharge voltage as an unknown that has to be found from the condition that the discharge current (an integral characteristic of the solution) takes the given value I 1 . This problem is nonlinear and, in principle, requires iterations. However, since the initial approximation being used represents the exact solution for I ¼ I 1 , only one iteration is needed without damping.
At the third step, the solution for I ¼ I 1 is recalculated with the external ionization term W in Eqs. (11) and (12) switched off. The convergence is very fast provided that ΔU is sufficiently small (typically, no more than three iterations are needed without damping if ΔU does not exceed 1% of U ). The discharge voltage will slightly increase (by an amount smaller than ΔU) and the obtained value will represent a little more accurate estimate of the self-sustainment voltage.
At the fourth step, the solution for I ¼ I 1 is recalculated with account of the space charge, i.e., with the Laplace equation replaced by the Poisson equation, and with account of all relevant nonlinear processes and with the stabilization activated as needed (if it has not been activated right from the first step). The change in U will be very small and only one iteration is needed provided that I 1 is small enough, say, in the order of nanoamps or smaller, which can be ensured by rescaling the external ionization term W, used in the resonance method. As the outcome of this step, one obtains a complete solution describing the onset of the selfsustained discharge.
At the fifth and final step, the quasi-stationary low-current discharge is computed in the whole range of its existence. The simulations are performed with the discharge current being gradually increased. The solution describing the onset of the self-sustained discharge, obtained at the previous step, is used as a starting point, and the solution for each current is used as an initial approximation for simulations with for next current value.
We conclude this section with a few practical hints. In gas discharge modeling, the species conservation and transport equations are frequently solved in the logarithmic formulation (e.g., this option is available in the Plasma module of COMSOL Multiphysics®), where the dependent variables are logarithms of the species number densities. The logarithmic formulation ensures that the number density of any of the species is never negative. However, such formulation introduces additional nonlinearity and was found to be less efficient for steady-state modeling than the original formulation (the one with the dependent variables being the species number densities). The modeling reported in this work has been performed in the original formulation and no sizable negative values of species densities have appeared provided that the numerical mesh is not too coarse.
It is essential that the code allows specifying not only the discharge voltage as a control parameter but also the discharge current with the possibility of an easy and seamless switching between the two. A standard way to specify the discharge current in gas discharge modeling is an implicit one, which involves introducing an external circuit comprising a voltage source and a ballast resistance. In COMSOL Multiphysics®, an alternative is available: one can use the "Global Equation node" option to specify discharge current directly without introducing an external circuit.
The modeling of the quasi-stationary low-current discharge in the whole range of its existence (the fifth above-described step) starts with the discharge current I being the control parameter. As I gradually increases, it may be helpful to switch the control parameter to U, in order to accelerate convergence. When simulating corona discharges, this can usually be done when the discharge voltage has increased by about 200 V from the ignition voltage.
It often happens in modeling that iterations converge for a value of the control parameter, but fail to converge for the next value, no matter how small the increment of the control parameter. Since a solution can turn back or join another solution but cannot just disappear, such a break-off represents a failure of the method. The most frequent reason is that there is a region of fast variation, which is not adequately resolved by the numerical mesh being used (e.g., the corona attachment to the electrode has expanded and is no longer adequately resolved). A refinement of the mesh solves the problem. Adaptive mesh refinement is a powerful tool.
The second most frequent reason is that an extreme point of the CVC or a turning point has been encountered: a code cannot pass through these points if operated with, respectively, U or I as a control parameter. An obvious fix is to switch the control parameter. If the CVC has a complex form, the control parameter has to be switched several times in order to compute the whole range of existence of the steady-state solution. On the other hand, not all sections of a complex-form CVC are stable and stability analysis is indispensable. An example of such a situation is encountered in Sec. V.

B. Coaxial and point-to-plane corona discharges
As examples, CVCs of positive and negative corona discharges in atmospheric-pressure air, computed by means of the integrated approach described in Sec. IV A, are shown in Figs. 4 and 5 for the concentric-cylinder and point-to-plane corona configurations, respectively. Also shown are time-averaged CVCs obtained experimentally. Agreement between the modeling and the experiment is quite good.
Computed distributions of the excitation rate of radiating nitrogen states N 2 C ð Þ, which are the main source of visible radiation emitted by air coronas, are shown in Fig. 6 for the point-to-plane corona and various applied voltages of both polarities along with the corona images recorded in the experiment. It can be seen that the modeling correctly reproduces the features of the radiation intensity distributions for both polarities.
The CVCs shown in Figs. 4 and 5 have been computed by means of stationary solvers in the framework of the approach described in Sec. IV A. In principle, the CVC of the positive coronas, as shown in Figs. 4 and 5(a), could be computed also by means of time-dependent solvers, which are the usual tool in gas discharge modeling. However, the switching from a stationary solver to a time-dependent one causes a very significant increase of the computation time, typically by several orders of magnitude, as documented in Ref. 28. This is unsurprising given that the mesh element size in the modeling was in the order of fractions of micrometer near the surface of the corona electrode, so the Courant-Friedrichs-Lewy criterion or analogous limitations require extremely fine time stepping.
On the other hand, an attempt to compute the time-averaged CVCs of the negative coronas shown in Fig. 5(b) by means of a time-dependent solver, reported in Ref. 28, was unsuccessful. This is an indication that the glow negative corona is unstable in these conditions. This conclusion is consistent with the wellknown fact that the negative DC corona in air usually operates in a pulsed regime: the current waveform reveals the so-called Trichel pulses. Therefore, the agreement between the steady-state modeling and experimental data seen in Figs. 4, 5(b), and 6 shows that a steady-state model provides a reasonably accurate description of time-averaged characteristics of pulsed negative coronas. Another example of an agreement of such type may be found in Ref. 29: the electric field distribution in a negative corona in an axially symmetric concentric-cylinder gap in air, computed by means of a steady-state model, was shown to be in good agreement with the time-averaged measurements, 30 performed in a regime with repetitive Trichel pulses. 31 Thus, in the case of a pulsed corona, time-dependent solvers, while being an adequate tool for studies of spatiotemporal evolution of individual pulses (e.g., Refs. 31-35), cannot be used for the direct calculation of time-averaged characteristics of the discharge. The latter can be conveniently computed by means of stationary solvers in the framework of the integrated approach being used in this section.

C. Discharge along a dielectric surface
Consider now a case where the discharge active zone is in contact with a dielectric surface. The formulation of the eigenvalue problem describing the ignition of a self-sustaining discharge, given in Sec. III A, needs to be slightly modified in such cases, as well as the procedure of solving it by means of the resonance method.
Continuing to consider low-current quasi-stationary discharges (which implies, in particular, that the time of variation of the applied voltage is much longer than the ion drift time), one can treat surface charges on a dielectric surface as quasi-stationary. Therefore, the surface is under the floating potential and the boundary condition for the electric field follows from the condition of the normal component of the local current density being equal to zero at each point of the surface. This boundary condition is linear and homogenous with respect to the ion and electron densities. Therefore, while distributions of the ion and electron densities in low-current quasi-stationary self-sustained discharges vary proportionally to the discharge current I, the electric field distribution is independent from I. This is similar to what happens in the case where no dielectric surfaces border the active zone of the discharge, considered in Sec. III A. It is this independence of the electric field distribution from I that makes the concept of self-sustaining (ignition) voltage applicable in the presence of dielectric surfaces. The difference from the no-dielectric case is that the electric field in the presence of a dielectric is coupled to distributions of the charged particles and cannot be computed independently.
Thus, one needs to consider a stationary boundary-value problem involving partial differential and possibly also integrodifferential equations, governing distributions of ions and electrons in the gap, and the Laplace equation governing the electrostatic potential with the aim to determine the value U ¼ U 0 of the applied voltage such that the problem admits a non-trivial solution for the particle densities. Note that this problem becomes ill-posed in the no-discharge situation, where the ion and electron densities are zero and the boundary condition of zero normal component of the current density at the dielectric becomes trivially satisfied and brings no information concerning the electric field. Therefore, it is unclear whether this problem can be solved by means of standard eigenvalue solvers. However, it can be efficiently solved by means of the resonance method. The procedure is as follows. The first and second steps of the procedure described in Sec. IV A are performed with the boundary condition of zero normal component of the electric field at the dielectric surfaces (instead of zero normal component of the current density). As in the case where no dielectric is present, the electric field distribution is decoupled from the ion and electron distributions and is linear with respect to U, so it is sufficient to perform the electrostatic simulation only once.
At the next step, the obtained solution for I ¼ I 1 is recalculated with the boundary condition of zero normal electric field at the dielectric surfaces being replaced by the condition of zero normal component of the current density. This involves solving the equations governing the distributions of the ions and the electrons in the gap and the Laplace equation for the electrostatic potential, and these equations are now coupled through the boundary condition of zero normal current density. The problem is nonlinear, and the distributions computed at the previous steps are used as an initial approximation.
Then, the third to fifth steps of the procedure described in Sec. IV A are performed, and, thus, the discharge is computed in the whole range of its existence. These steps do not need to be modified due to the presence of the dielectric and are performed in the same way as described in Sec. IV A.
As an example, Fig. 7 shows the self-sustainment voltage computed by means of the above procedure (excluding the fourth and fifth steps of Sec. IV A) for an axially symmetric configuration shown in the inset: two disk electrodes separated by a cylindrical dielectric surrounded by atmospheric-pressure air. Note that the dielectric is depicted as a massive cylinder, however this is irrelevant, as is its dielectric constant, since the equations, being stationary, are solved only in the air-filled region. The coefficient of secondary electron emission from the cathode and the dielectric was 0:03. The simulations have been performed for the dielectric radius R in the range from 3 to 8:2 mm. The upper edge of the positive electrode and the lower edge of the negative electrode were rounded with a 0:2 mm radius as indicated in the inset. For R ! 7:7 mm (i.e., when the dielectric sticks out by 0:2 mm or more, which is the case shown in the inset), the edges of the dielectric were rounded with a 0:2 mm radius as well. For R 7:3 mm (i.e., in the case where the dielectric is recessed by 0:2 mm or more), there were roundings at the lower edge of the positive electrode and the upper edge of the negative electrode but not at the edges of the dielectric. Finally, for 7:3 mm , R , 7:7 mm, there were roundings at the dielectric edges as well as at the lower edge of the positive electrode and the upper edge of the negative electrode. The radius of the roundings was 0:2 mm, and the centers of curvature were positioned in such a way that the surfaces of the electrode and dielectric formed an angle of 90 at the triple junction. Thus, the contact angle between the metal and dielectric surfaces at the triple junction was 90 in all the geometries with the aim to limit the local electric field.
It is seen from Fig. 7 that there is little difference between the values of the self-sustainment voltage computed with the boundary condition at the dielectric surface being zero normal current density or zero normal electric field. For R & 5 mm, when the dielectric is away from the edges of the electrode disks, the selfsustainment voltage is virtually constant and equals approximately 10 kV. This value is in good agreement with the experimental breakdown voltage for R ¼ 3:5 mm reported in Ref. 36, which is represented in Fig. 7 by a triangle. Note that this value is lower than the corresponding Paschen value (approximately 11 kV), which is consistent with the discharge being not plane-parallel (the breakdown path computed for R & 5 mm represents an arc attached to the electrode edges). As R increases, i.e., as the dielectric approaches the active zone of the discharge, the self-sustainment voltage decreases and attains a minimum when R becomes equal to the electrode radius. For higher R, the self-sustainment voltages starts increasing again.
There seem to be no quantitative experimental data published on the effect of dielectric surfaces on self-sustainment voltages. It is known from experiment that the presence of a dielectric facilitates the breakdown (e.g., review 37 ), and it is indeed seen in Fig. 7 that the self-sustainment voltage decreases as the dielectric approaches the active zone of the discharge. It would be very interesting to check experimentally the rapid increase of the self-sustainment and/or breakdown voltage when the dielectric sticks out, predicted by the modeling and seen in Fig. 7 in the range R ! 7:5 mm.

V. STABILITY OF LOW-CURRENT QUASI-STATIONARY DISCHARGES
The integrated approach to the modeling of quasi-stationary low-current discharges, described in Sec. IV A, offers the possibility FIG. 7. Discharge between two disk electrodes separated by a cylindrical dielectric. Air, 1 atm. Circles, crosses: self-sustainment voltage, modeling with the boundary condition at the dielectric surface being zero normal current density (circles) or zero normal electric field (crosses). Triangle: breakdown voltage, experiment. 36 of calculating all existing stationary solutions. In other words, the integrated approach makes it possible to compute all the steady states of the discharge that are theoretically possible, regardless of their stability and whether they can be realized in a particular experiment. This is a consequence of the decoupling of physical and numerical stability, which is characteristic of stationary solvers. An example of the usefulness of this feature is the possibility of computation of time-averaged characteristics of pulsed negative coronas by means of a stationary solver, discussed at the end of Sec. IV B.
The stability of the computed steady states may be studied separately. Stability against small perturbations may be studied by means of solving the eigenvalue problem resulting from linear stability theory. Note that the eigenvalue problem governing stability of 2D (axially symmetric or planar) states against 3D perturbations may be formulated and solved in the 2D geometry, which greatly facilitates the task.
An alternative approach to the investigation of stability is to employ the same code that was used for the calculation of the steady-state solution with the stationary solver replaced by a timedependent one. A perturbation is applied to the steady-state solution, and the development of the perturbation is followed by means of the time-dependent solver. This approach does not involve solving an eigenvalue problem and allows studying stability against both small and finite perturbations. However, one should keep in mind that numerical stability of time-dependent solvers is not necessarily equivalent to physical stability: in Ref. 5, a number of examples were shown where a time-dependent solver of COMSOL Multiphysics was able to compute some steady states that are unstable in the framework of linear stability theory and/or was unable to compute some stable states. Moreover, time-dependent solvers cannot give any information concerning stability against perturbations of symmetries different from the one of the calculation domain of the solver. For example, 2D time-dependent simulations say nothing about stability of 2D steady states against 3D perturbations, which are frequently the most dangerous ones. Consequently, care should be employed in drawing conclusions about the stability of steady states on the basis of time-dependent simulations.
In Ref. 10, the main features of stability of axially symmetric steady states of the Townsend and glow discharges were investigated by means of solving the eigenvalue problem in the framework of linear stability theory, considering as an example conditions of microdischarges in xenon. Both real and complex eigenvalues have been detected, meaning that perturbations can vary with time both monotonically and with oscillations. The fundamental mode of the axially symmetric glow discharge is stable when the discharge is operated in the abnormal regime. There is a rather wide window of stability of the axially symmetric normal discharge, which is characterized by a normal spot at the center of the cathode. The subnormal discharge is unstable, the Townsend discharge is stable at low currents.
As another example, let us apply the approach based on the use of time-dependent solvers to study the stability of steady states of a positive point-to-plane corona discharge against perturbations of various amplitudes. It is expected that in this specific geometry, the most dangerous perturbations are axially symmetric ones, hence stability may, in principle, be studied by means of 2D timedependent modeling.
In Fig. 8, the CVC of a steady-state positive corona under conditions of experiment, 28 computed in a wide current range by means of a stationary solver, is shown. Also shown are the experimental data from Ref. 28, which refer to a more narrow current range and are the same as the data shown by triangles in Fig. 5(a). The CVC manifests a maximum, shown by the circle B, followed by a turning point C. Leaving analysis of reasons of such behavior beyond the scope of this work, let us consider stability of different steady states.
Let U 0 be the voltage of the steady state the stability of which we want to investigate. One can consider different forms of the perturbation imposed over the steady-state solution. In this work, the initial condition for the time-dependent modeling was obtained by running the steady-state code for U ¼ U 0 þ ΔU with a certain value of ΔU. (In other words, the perturbation was set equal to the difference between the steady-state solutions for U ¼ U 0 þ ΔU and U ¼ U 0 .) The time-dependent code was invoked with this initial condition and was run with the fixed corona voltage U ¼ U 0 . If in the course of temporal evolution the perturbation decays and the discharge relaxes to the steady state U ¼ U 0 , then this state is stable against the considered perturbations and for the voltagecontrolled discharge. Otherwise, the state is unstable.
Results of the calculations for ΔU ¼ +20 V are schematically shown in Fig. 9. Here, f designates a state positioned on the ascending branch of the CVC (branch AB in Fig. 8). On being perturbed with +20 V (states h and i, respectively), the discharge relaxes to the state f as shown by the arrows, which means stability. The same outcome is obtained in calculations for other values of ΔU, provided that ΔU j j is not too big. The same outcome is obtained also for other steady states on the ascending branch. The states after the maximum of the CVC (on the branch BCD in Fig. 8) are unstable. In some cases, the discharge relaxes to the steady state with the same U on the ascending branch on the CVC. An example is shown in Fig. 9: on being perturbed with +20 V (states j and k, respectively), the discharge relaxes not to the state g but rather to f . In other cases, no steady state was attained and streamer-like structures appeared to start forming.
The above fits very well into the conventional pattern of stability of voltage-controlled discharges against small perturbations: the states on the ascending branch of the CVC are stable, and the states beyond the maximum of the CVC are unstable.
On the other hand, the states on the ascending branch of the CVC are unstable against perturbations of a bigger amplitude: simulations for higher values of ΔU j j have shown that stability is lost (and streamer-like structures appear to start forming) as ΔU j j reaches a threshold value, which is between 150 and 200 V for lower U 0 and decreases as U 0 increases (e.g., overvoltages of 80 V or more in the state f triggered the formation of streamer-like structures). In other words, the states on the ascending branch of the CVC are stable if the perturbations of the corona voltage are within approximately 5% for lower voltages and about 1% for higher voltages. This may explain why a stable positive corona discharge is more difficult to realize for high corona voltages. In any case, it would be very interesting to investigate experimentally the effect of voltage stabilization on positive corona.

VI. SUMMARY
The work aims at developing an integrated approach for the computation of low-current quasi-stationary discharges from the inception to a non-stationary transition to another discharge form, such as a transition from the Townsend discharge to a normal glow discharge or the corona-to-streamer transition. This task involves three steps: (i) modeling of the ignition of a self-sustaining discharge, (ii) modeling of the quasi-stationary evolution of the discharge with increasing current, and (iii) the determination of the current range where the quasi-stationary discharge ceases to exist and the abovementioned non-stationary transition begins. Each of these three steps is considered in some detail with a number of examples, referring mostly to discharges in atmospheric-pressure air.
The ignition of self-sustained gas discharges is governed by a multidimensional boundary-value eigenvalue problem for a system of stationary partial differential equations (and possibly also integrodifferential equations), formulated in Sec. III A. The physical sense of the problem is that the applied voltage should be such that direct ionization by the impact of electrons accelerated by the applied electric field and photoionization by photons generated in the discharge are just sufficient to compensate for the losses of charged particles. There are two special cases where the applied electric field may be considered as uniform and the formulated problem is reduced to one dimension. The first one is the case of discharge ignition between wide parallel-plate electrodes. The formulated eigenvalue problem leads to the well-known self-sustainment condition (7) in this case, provided that similar simplifications are applied. The second special case is the one of self-sustainment field in the cross section of a plane or cylindrical positive column of a long low-current discharge. In this case, the formulated problem takes a form similar to the well-known eigenvalue problem for an ordinary differential equation describing the ambipolar diffusion of charged particles in the cross section of a positive column under the assumption of quasi-neutrality. The difference is that at the discharge ignition, the transversal electric field is zero, the diffusion is free and not ambipolar, and there is no quasi-neutrality.
In the general case, the formulated eigenvalue problem for the ignition of a self-sustaining discharge is multidimensional and its numerical solution is not trivial. Three different methods for its solution are described in Secs. III B-III D: direct solution of the eigenvalue problem for the self-sustainment voltage, investigation of stability of the no-discharge solution, and the resonance method. The methods are mathematically equivalent as far as computation of the self-sustainment voltage is concerned, so the question is which method is easier to implement and use and which method has wider applicability. In the authors' experience, the resonance method is the preferred one. The method is based on solving linear partial differential equations and can be implemented with the use of standard solvers, including commercial ones. It is robust and fast and has given useful results in a wide range of conditions and complex geometries.
An important advantage of the resonance method is its physical transparency. In particular, it allows one to use this method not only for simulations of ignition of self-sustaining discharges but also as a module of a more general code for the modeling of quasi-stationary self-sustained discharges. Such an integrated approach, based on a combination of the resonance method and stationary solvers, is described in Sec. IV A and allows modeling of quasi-stationary lowcurrent discharges in the whole range of their existence starting from the inception. The use of stationary solvers instead of time-dependent ones dramatically reduces the computation time, and this reduction is especially large in discharges with strongly varying length scales, such as corona discharges, where a variation of the mesh element size by orders of magnitude is indispensable. The integrated approach to the modeling of quasi-stationary low-current discharges offers the possibility of calculating all existing stationary solutions, or, in other words, all steady states of the discharge that are theoretically possible, regardless of their stability and whether or not they can be observed in a particular experiment. An example of the usefulness of this feature is the possibility of computation of time-averaged characteristics of pulsed negative coronas by means of a stationary solver, discussed at the end of Sec. IV B. Note that time-dependent solvers, while being an adequate tool for studies of spatiotemporal evolution of individual pulses, cannot be used for the direct calculation of time-averaged characteristics of a pulsed corona.
The stability of the computed steady states may be studied separately. Stability against small perturbations may be studied by means of solving the eigenvalue problem resulting from the linear stability theory. An alternative approach to the investigation of stability is to employ the same code that was used for the calculation of the steady-state solution with the stationary solver replaced by a time-dependent one; and a perturbation is applied to the steady-state solution and the development of the perturbation is followed by means of the time-dependent solver. An example of application of a time-dependent solver for the investigation of stability of a positive point-to-plane corona discharge against perturbations of various amplitudes is given in Sec. V.
A challenging task is to develop, on the basis of the resonance method, an accurate and fast tool for evaluation of the selfsustainment voltage in conditions of industrial interest. Such tool would provide a reference point for the optimization of hold-off capabilities of high-voltage switchgear operating in low-frequency fields. It would be a useful supplement to the technique based on evaluation of the Townsend ionization integral, which is the main tool used in engineering practice. Another challenging task is to establish a connection between the self-sustainment and breakdown voltages. Note that in certain conditions the two voltages coincide, e.g., Ref. 9. It would be very interesting to check experimentally the rapid increase of the self-sustainment and/or breakdown voltage in cases where the dielectric protrudes into the plasma, predicted by the modeling and seen in Fig. 7 in the range R ! 7:5 mm. Another interesting unresolved issue is to investigate, through an experiment supported by a modeling similar to that described in Sec. V, the effect of voltage stabilization on the positive corona. The authors have no conflicts to disclose.

APPENDIX A: BOUNDARY CONDITIONS FOR DRIFT-DIFFUSION EQUATIONS
The hydrodynamic (drift-diffusion) equations (2) are justified provided that the local macroscopic length scale L is much larger than the relevant microscopic scale λ α . When α refers to an ion species, λ α is represented by the ion mean free path. When α ¼ e, λ e is represented by the electron energy relaxation or maxwellization length. On distances of the order of λ α or smaller from solid surfaces, the hydrodynamic equations lose their validity. Therefore, a kinetic analysis is indispensable in order to formulate boundary conditions on solid surfaces for hydrodynamic equations in a rigorous way: a kinetic equation for the velocity distribution function on distances of the order of λ α must be solved in the first approximation in the small parameter λ α =L and the obtained solution must be asymptotically matched with a solution of the hydrodynamic equations on the scale L; and the procedure of the asymptotic matching will provide the boundary condition being sought.
Let us denote by the index n the projections of the corresponding vectors onto the local normal directed from the solid into the plasma. Let us first consider the case where (i) E n , the local normal electric field at the surface, is weak enough so that its work over the charged particles of species α over the distance λ α is much smaller than the characteristic translational energy of the particles, e E n j jλ α ( kT α ; (ii) all particles of species α coming to the surface from the plasma are absorbed and none are reflected; and (iii) the surface does not emit the particles of species α. (Here, T α is the species effective temperature and k is Boltzmann's constant.) In this case, the distribution function on distances of the order of λ α is of the order of λ α =L, i.e., asymptotically small, compared to values of the distribution function on distances of the order L. Therefore, the above-described asymptotic matching procedure requires that n α the particle number density governed by the hydrodynamic equations vanishes at the surface, thus giving the trivial boundary condition n α ¼ 0 for the hydrodynamic equations. Note that this boundary condition has exactly the same physical meaning as the standard no-slip boundary condition on solid surfaces for the Navier-Stokes equations. If the local electric field is moderate or strong, e E n j jλ α * kT α , then the asymptotic structure of the solution on the distances of the order λ α is rather complex for both the electrons 38-40 and the ions; 41 no simple exact solution is possible, and, therefore, there is no unique way to formulate simple boundary conditions. This explains why the number of different existing variants of boundary conditions is so large, e.g., Refs. 42, 43, and references therein. On the other hand, the effect of the boundary conditions over the distribution of particles attracted to the surface (the positive ions in the case E n , 0 and the negative ions and the electrons in the case E n . 0) is localized on the length scale kT α =e E n j j in the case of moderate or strong electric fields. This scale is comparable to, or smaller than, λ α , and the hydrodynamic (drift-diffusion) equations are anyway inapplicable on this scale. The effect of the boundary conditions on the distribution of the repelled particles will be qualitatively correct if the flux of particles emitted by the surface is described correctly. Hence, the exact form of the boundary conditions is not very important in the case of moderate or strong electric fields, and this is consistent with what is found in the modeling practice. Therefore, the use of simple approximations is advisable.
The simplest variants may be summarized as follows. The hydrodynamic transport flux of particles of species α in the direction from a solid surface into the plasma may be represented as where J αÀ is the density of the flux of the particles coming from the plasma to the surface, r α is the reflection coefficient, and J α ð Þ em is the density of the flux of the particles emitted by the surface. Considering r α and J α ð Þ em as known, one needs to express the kinetic quantity J αÀ in terms of hydrodynamic parameters in order to obtain an explicit boundary condition for hydrodynamic equations.
The simplest approximation is to assume that J αÀ is proportional to n α the local number density and the mean speed of chaotic motion, where ξ α is a numerical coefficient and m α is the species particle mass. Equation (A1) assumes the form It is well known that if the velocity distribution of particles of species α is isotropic and Maxwellian, then the density of particle flux in any direction, in particular, in the direction to the surface, equals n α C α =4, so ξ α in Eq. (A2) may be set equal to 1=4. The condition obtained in this way is equivalent to the so-called Thomson-Loeb formula, which has been widely used in fluid modeling of gas discharges, e.g., Ref. 44 and references therein.
If the number of particles moving from the surface is low, then the velocity distribution near the surface is strongly anisotropic. Approximating this distribution by a half-Maxwellian function, one obtains ξ α ¼ 1=2. Note that the obtained boundary condition in the case of the electrons and r e ¼ 0 coincides with the corresponding boundary condition employed in Plasma module of COMSOL Multiphysics®.
Another variant of boundary condition may be obtained by representing the velocity distribution near the surface as a superposition of two distributions, describing the particles moving to and from the surface. Assuming that each of these distributions is half-Maxwellian, one can write where n αþ and n αÀ are the number densities of the particles moving into the plasma and to the surface, respectively. Solving these equations for n αþ and n αÀ and substituting into the expression for the net particle number density, n α ¼ n αþ þ n αÀ , one finds Solving this equation for J αÀ and substituting into Eq. (A1), one obtains the boundary condition, This condition has been known for a long time, e.g., in the particular case r α ¼ 0, it coincides with Eq. (38) of Ref. 45; see also recent work. 43 The second term on the rhs of Eq. (A5) is consistent with the second term on the rhs of Eq. (A2) if one sets ξ α ¼ 1=2 1 þ r α ð Þ. However, the first term on the rhs of Eq. (A5) involves the coefficient 2= 1 þ r a ð Þ, which is absent in Eq. (A2) and which exceeds unity except for r a ¼ 1.
From now on, let us restrict the consideration to the usual situation where the reflection of the charged particles is insignificant, i.e., r α ¼ 0 for all charged-particle species. Let us apply each of the above boundary conditions to three limiting cases. The first one is the case of local equilibrium, where most of the emitted particles return to the surface: J αn j j ( J a ð Þ em . Equation (A2) with ξ α ¼ 1=4 and Eq. (A5) give the correct result n α ¼ 4J The second limiting case is the one where most of the emitted particles are swept away by the electric field into the bulk of the plasma and only a small number return to the surface. The drift velocity of the particles in the vicinity of the surface is directed into the plasma and is much greater than C α , so the second term on the rhs of Eqs. (A2) and (A5) is small. Equation (A2) assumes the form J αn ¼ J The former expression would be adequate for an isotropic velocity distribution function. However, the distribution near an absorbing non-emitting surface is strongly anisotropic, and the latter expression appears to be more appropriate in such cases.
Thus, each of the boundary conditions correctly reproduces two of the limiting cases but not the third one. Since the boundary condition (A5) gives an unrealistic result in the important limiting case where most of the emitted particles are swept away by the electric field into the plasma, this condition can be hardly recommended. Equation (A2) with ξ a equal to 1=4 or 1=2 appears to be more physical.
In many cases, the ion emission is neglected, i.e., J α ð Þ em for all ion species may be set equal to zero. Then, Eq. (A2) assumes the form for the ions and the electrons, respectively. The electron emission flux J e ð Þ em must be evaluated with account of all relevant mechanisms, including the electron emission under the effect of charged and excited particles and photons produced in the discharge (secondary electron emission), photoemission caused by external radiation, thermionic emission, thermo-field, and field electron emission from the cathode surface. It should be stressed that in the case of weak electric field, where e E n j jλ α ( kT α or, equivalently, the drift speed is much smaller than C α , the first boundary condition in Eq. (A6) is equivalent to the trivial boundary condition n α ¼ 0 as it should.
The natural choice is to set ξ a ¼ 1=2 in the first equation (A6). ξ e in the second equation may also be set equal to 1=2 except at hot surfaces in high-pressure (arc) plasmas, where the dominating electron emission mechanism is thermionic emission and a significant part of the emitted electrons return to the surface, hence the value ξ e ¼ 1=4 may be more appropriate. The boundary conditions (A6) with ξ a ¼ ξ e ¼ 1=2 were used in most of the simulations described in this work.
For computational reasons, it may be convenient to replace the boundary condition (A6) for the positive ions at the cathode and negative ions and the electrons at the anode by the condition of zero normal derivative of n α , which amounts to neglecting the diffusion flux of the attracted particles to the electrode in comparison with the drift flux. This eliminates difficulties with the computation of the distribution of the attracted particles on the length scale kT α =e E n j j, where the solution given by the hydrodynamic equations is non-physical anyway. This simplification was employed in most of the simulations described in this work.
As mentioned above, the effect of the precise form of boundary conditions in the modeling is not strong. As an example, one can refer to the parallel-plate discharge, considered in Sec. III E: the discharge initiation voltage, computed numerically using the boundary conditions Eq. (A6) with ξ i ¼ ξ e ¼ 1=2, differs from the value that follows from Eq. (13), which was obtained without account of diffusion with the boundary conditions J en ¼ ÀγJ in at the cathode and n i ¼ 0 at the anode, by mere 1:2%. (The index i here refers to the positive ions.) As another example, one can mention that the change of the boundary conditions for the negative ions at the cathode and for the positive ions at the anode from the condition (A6) with ξ α ¼ 1=2 to n α ¼ 0 produced little effect in the modeling of corona discharges reported in Ref. 28.
3:3 Â 10 À4 and 4:9 Â 10 À4 m 2 V À1 s À1 over the reduced field range up to 120 Td. In principle, variations of reduced mobilities with the reduced field can be readily introduced into numerical models. Since, however, a constant value is used for the mobility of O À 2 (and complex/cluster ions) and given that the above variations are not huge, a constant value of 5:2 cm 2 V À1 s À1 is chosen for the reduced mobility of O À in order to be consistent. Note that this value corresponds to the value given in Table IIf (p. 81) of Ref. 48 for the reduced field of 100 Td.
The O À 3 ions are present mostly in the drift zone at low reduced fields. According to Ref. 48, the reduced mobility of O À 3 varies between approximately 2:7 Â 10 À4 and 3:4 Â 10 À4 m 2 V À1 s À1 over the reduced field range up to 100 Td and between approximately 3:5 Â 10 À4 and 3:3 Â 10 À4 m 2 V À1 s À1 over the reduced field range between 100 and 200 Td. According to Ref. 50, the reduced mobility of O À 3 varies between 2:5 Â 10 À4 and 3:1 Â 10 À4 m 2 V À1 s À1 over the reduced field range up to 100 Td and between 3:1 Â 10 À4 and 2:9 Â 10 À4 m 2 V À1 s À1 over the reduced field range between 100 and 240 Td. In order to be consistent, we have chosen for the reduced mobility of O À 3 ions a constant value of 2:7 Â 10 À4 m 2 V À1 s À1 , which corresponds to the value shown in Fig. 2 of Ref. 50 for the reduced field of 60 Td.
The diffusion coefficients of all ion species are related to the mobilities through Einstein's relation with the corresponding effective ion temperatures evaluated by means of the Wannier formula, Eq. (C1) of Appendix C. We remind that both Einstein's relation and the Wannier formula are accurate in the case of ions with a constant mobility. The mobility of the electrons was taken from Ref. 51, and the longitudinal and transversal diffusion coefficients of the electrons were evaluated with the use of the online version of the Bolsig+ solver 52 and the cross sections. 53 The kinetic scheme and relevant kinetic data used in this model are summarized in Table I. Reactions (1)-(4) and (6)- (8) are the same that were considered in Ref. 12 where M is any of the molecules N 2 and O 2 . However, the contribution of the process O À 2 þ N 2 is small; therefore, the collisional detachment from O À 2 is written in this model with account of collisions only with O 2 . The approximations of rate constants of reactions 5-8 were taken from Table 2 of Ref. 54. (The rate constant of the collisional detachment from O À 2 was multiplied by the factor of 5, in agreement with the above.) It should be stressed that these approximations are valid for variable gas temperature T, in contrast to the approximations given in the Appendix of Ref. 54 and used in Ref. 12, which are valid only for T ¼ 300 K.
Note that, although numerical results reported in this paper refer to the constant gas temperature T ¼ 300 K, the applicability of the modified model for variable T is an improvement that will be exploited in subsequent work. It is expected that the modified model is applicable under conditions where the degree of dissociation of oxygen molecules is sufficiently low and oxygen atoms do not significantly affect the balance of charged particles (in particular, the rate of destruction of negative ions). At air pressures of the order of 1 atm, this corresponds to gas temperatures of the order of 1000 K and lower.
Another modification to the kinetic scheme 12 is an account of recombination. Again, this modification is not very relevant to the results reported in this paper; however, it will be useful for subsequent work. It should be stressed that due to the lack of sufficient experimental information, the account of recombination cannot be introduced in an accurate way and should be considered rather as Note that the third factor in the first term on the rhs of Eq. (C2) is the so-called reduced temperature of the species α and β, so the physical meaning of this equation is clear. Note also that Eq. (C2) is consistent with the well-known fact that the mean kinetic energy of relative motion of ions and neutrals is characterized by the effective ion temperature: setting in (C2) v dβ ¼ 0, m β ¼ M, T β ¼ T, one obtains T αβ ¼ T α as it should be.
Strictly speaking, the use of the effective reduced temperature T αβ for the evaluation of rate constants is justified in the case of binary ion-ion reactions. However, in the absence of better options it is natural to use these temperatures also in the case of three-body reactions, where the third body is a neutral particle, in the same way as the effective ion temperature is used for the evaluation of rate constants of three-body ion-molecular reactions with the third body being a neutral particle. Note that in the particular case of ion-ion recombination reactions, the factor v dα À v dβ À Á 2 in (C2) may be replaced by v dα þ v dβ À Á 2 .

DATA AVAILABILITY
The data that support the findings of this study are available within the article.