Simple computation of ignition voltage of self-sustaining gas discharges

A robust, fast, and accurate numerical method is proposed for finding the voltage of the ignition of DC self-sustaining gas discharges in a wide range of conditions. The method is based on physical grounds and builds up from the idea that the ignition of a self-sustaining gas discharge should be associated with a resonance that would occur in a non-self-sustained discharge in the same electrode configuration. Examples of the application of the method are shown for various configurations: parallel-plate discharge, coaxial and wire-to-plane corona discharges, and a discharge along a dielectric surface. The results conform to the conventional Townsend breakdown condition for the parallel-plate configuration and are in good agreement with existing experimental data for the other configurations. The method has the potential of providing a reference point for optimization of the hold-off capability of high-power switchgear operating in low-frequency fields.


Introduction
The condition of initiation of a self-sustaining DC gas discharge, where the discharge voltage is just sufficient for the electron impact ionization to compensate losses of the charged particles, is well known in gas discharge physics, but goes under different names. For example, the terms 'self-sustaining condition' and 'breakdown condition' are encountered on p 544 of the book [1] and the terms 'condition for initiating a self-sustaining discharge' and 'criterion of ignition (selfsustainment)' are found on pp 131 and 175 of another classical book [2]. Note that in electrical engineering the term 'breakdown' has a different meaning: 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. In this work, the term * Author to whom any correspondence should be addressed.
'ignition of self-sustaining discharge' is used following [2] and 'breakdown' refers to transitions from low-to high-current discharge.
In addition to being of theoretical interest, understanding of ignition of self-sustaining discharges is important for applications: it is a useful reference point in investigation of breakdown in high-voltage electrical equipment in low-frequency (e.g., 50 Hz) electric fields, where the time of variation of the applied voltage is much longer than the ion drift time. Moreover, in certain conditions the breakdown voltage coincides with the voltage of ignition of self-sustaining discharge; e.g., [3].
An accurate derivation of the ignition condition exists for the parallel-plate geometry, where the applied electric field is uniform and diffusion of the charged particles is of minor importance. The condition reads αL = ln 1/γ + 1 , where α is the Townsend ionization coefficient, L 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 αL by the so-called ionization integral, which is the line integral of α evaluated along an electric field line; the field line that ensures the biggest value of the integral is chosen. However, such approach is theoretically incomplete and cannot be generalized to account for other potentially important effects without invoking more or less arbitrary assumptions. One of potentially important effects disregarded by the above approach is photoionization; for example, in positive coronas in electronegative gases the role of the secondary electron emission from the cathode is minor and seed electrons, which are subsequently multiplied by the impact ionization in the region of high electric field, are generated by photoionization. Another effect disregarded in the above approach is diffusion losses of the charged particles to walls, which may play a role, in particular, in the ignition of discharges along dielectric surfaces.
The aim of this work is to show that the problem of ignition of self-sustaining discharges may be solved by accurate mathematical means, and this solution is relatively simple. At the ignition, the charged particle densities are very low and the applied electric field is unperturbed by plasma space charge. There is no need to consider excited states and multi-step processes. Hence, the equations describing charged particle distribution in the gap are linear. One needs to find an applied voltage value where the volume ionization balances losses of the charged particles. From the mathematical perspective, this is a boundary-value eigenvalue problem for a system of stationary linear partial differential equations (and possibly also integro-differential equations) describing the distribution of the charged particles in the gap, the ignition voltage being the eigenparameter.
Thus, from the mathematical point of view, the problem of ignition of self-sustaining discharges is an eigenvalue problem. Acknowledging and exploiting this fact is appropriate and should be beneficial. However, this fact does not seem to have been recognized in the gas discharge theory, although the eigenvalue problem has the same physical meaning as the eigenvalue problem for a system of ordinary differential equations governing the electric field in the plane or cylindrical positive columns of low-current discharges, which is well known and described in textbooks; e.g., [1,4]. Related is the work [5], which proposed to determine the breakdown voltage by computing stability spectra for various values of the applied voltage.
In this work, an efficient physics-based method is proposed for numerical solution of the eigenvalue problem governing the ignition voltage of self-sustaining discharges. The method consists in modelling of a non-self-sustained lowcurrent discharge, generated by an external ionization source in the same plasma-producing gas in the same discharge configuration. The applied voltage is gradually increased until a resonance-like behaviour appears, and the voltage value at which this happens represents the sought-for ignition voltage of the self-sustaining discharge.
One can hope that the accurate and fast tool for evaluation of the ignition voltage, developed in this way, will be a useful supplement to the technique based on evaluation of the Townsend ionization integral, which is the main tool in studies of the hold-off voltage in engineering practice. Moreover, this approach can be used to produce an initial approximation for the modelling of the (already ignited) self-sustaining discharge by means of a stationary solver.
The outline of the paper is as follows. In section 2, the eigenvalue problem is stated and a method of its numerical solution is proposed. Calculation results for a parallel-plate discharge, coaxial and wire-to-plane corona discharges, and a discharge along a dielectric surface are given in section 3. The results conform to the conventional Townsend breakdown condition for the parallel-plate configuration and are in good agreement with existing experimental data for the other configurations. Conclusions are summarized in section 4.

The eigenvalue problem
The aim is to theoretically determine the ignition (inception) voltage of gas discharges of arbitrary configuration and at arbitrary pressure. The space charge in the interelectrode gap is negligible at the discharge inception. 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 the section 3.3 below.) 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 simulations only once and then rescale 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.
The distribution of ions and electrons in the gap is described by a system of transport and conservation equations, written for each of the ion species and the electrons, coupled with other relevant equations, such as species energy equations or equations describing the photoionization. Nonlinear effects, such as recombination of charged particles, or reactions involving ions or electrons and excited particles or radicals produced in the discharge, or effects caused by the Joule heating of neutral gas, are negligible at the discharge inception. There is no external (entering into the discharge from outside) 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. According to the conventional definition, the ignition voltage is a value of the applied voltage such that the applied electric field produces an ionization intensity just sufficient to compensate losses of the charged particles, so a steady-state balance of the charged particles can be reached. Therefore, the above system of governing equations should be considered in the stationary form, i.e., without derivatives with respect to time.
The system of equations is supplemented by usual boundary conditions at electrode surfaces, such as those describing absorption and/or reflection of ions and electrons and secondary electron emission from the cathode or photoemission under the effect of radiation produced in the discharge (we consider conditions where thermionic, thermo-field, and field electron emission are negligible, as well as electron photoemission caused by external radiation).
The resulting boundary-value problem governing distributions of the ion and electron densities is considered for a given distribution of the electric field, so 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, i.e., a low-current steady-state discharge.
In mathematical terms, the above problem, which governs distributions of the ion and electron densities, represents a boundary-value problem for a system of partial differential or integro-differential equations and this problem is linear (with respect to the 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 corresponds 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.

Method of numerical solution
Several methods of solving the above eigenvalue problem have been tested, including those based on the use of a standard eigenvalue solver. The method of choice is physics-based and may be described as follows. The stationary boundary-value problem governing distributions of the ion and electron densities, described in section 2.1, 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. 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 1D problems and a Gaussian function with a maximum on the discharge axis is a natural choice in 2D problems. The obtained stationary boundaryvalue problem, which describes the non-self-sustained discharge, 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 ignition voltage of the self-sustaining discharge.
The stationary boundary-value problem describing the nonself-sustained discharge is linear and inhomogenous and may be routinely solved by means of ready-to-use solvers; a solver provided by commercial software COMSOL Multiphysics is used in this work.
Note that the idea of finding a real eigenvalue by means of varying the eigenparameter and looking for some or other kind of resonance is pretty obvious; e.g., 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 [6,7] and [8] 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.

Parallel-plate configuration
Let us consider a parallel-plate microdischarge in xenon at the pressure of 30 Torr. The interelectrode gap is L = 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 equations for the ions and the electrons are written in the drift-diffusion approximation, transport and kinetic coefficients are the same as in [9]. The secondary electron emission coefficient γ was set equal to 0.03.
The current-voltage characteristic (CVC) of the non-selfsustained discharge, j (U) (here j is the current density), was computed numerically for an external ionization source term that was uniform in space, and is depicted in figure 1 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 occurs is the voltage of ignition of the self-sustaining discharge.
An approximate analytical description of the CVC in this case may be obtained by neglecting the diffusion of the charged particles and is given by the formula where W is the external ionization source and α is the Townsend ionization coefficient. The CVC described by this equation conforms, as is should, to the CVC given by the numerical calculations, shown in figure 1. 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 equation (1) equal to zero, which is the well-known condition of ignition of the self-sustaining parallel-plate discharge mentioned in the introduction. This value is approximately 177.20 V, which duly conforms to the above-mentioned value of the ignition voltage of 175.05 V. Thus, the expected resonance-like behaviour of characteristics of the non-self-sustained discharge, which should occur at an applied voltage coinciding with the ignition voltage of the self-sustaining discharge, 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 a uniform mesh with 170 elements was sufficient to obtain the ignition voltage with five significant digits, hence the method is not too sensitive with respect to the mesh. Moreover, the computed ignition voltage was the same with five significant digits regardless of whether the default stabilization was activated or not in the Transport of Diluted Species module of COMSOL Multiphysics, which was used for solving the species conservation equations.
The pattern of the CVCs I(U) of the non-self-sustained discharge, seen in figure 1, has been observed in every case where the method has been applied, including examples not shown in this paper, so the resonance and thus the ignition voltage were readily identifiable and the method was straightforward, fast, and reliable in all cases.

Coaxial and wire-to-plane corona discharges
Appropriate tools for the modelling of quasi-stationary gas discharges are stationary solvers, which, in particular, are not subject to the Courant-Friedrichs-Lewy criterion or analogous limitations on the mesh element size. Although most of the popular ready-to-use codes for gas discharge simulation employ time-dependent solvers, e.g., nonPDPSIM [10] and Plasma module of COMSOL Multiphysics, stationary solvers for gas discharge modelling are provided by Plasimo [11]; COMSOL Multiphysics provides stationary solvers for general partial differential equations; although the Plasma module of COMSOL is intended to work with time-dependent solvers, it can still be used with stationary solvers [12]. Several examples of gas discharge modelling where stationary solvers offer important advantages compared to time-dependent solvers can be found in [12].
The most time-consuming step in modelling by means of stationary solvers in many cases is finding a suitable initial approximation. Fortunately, for the modelling of low-current self-sustaining discharges this step can be greatly accelerated by using the resonance method: the resonance method provides an accurate first approximation for a steady-state solution at a very low current, which describes the onset of the self-sustaining discharge, and this solution can be used as a starting point for the simulations with the discharge current being gradually increased, the solution for each current serving  as an initial approximation for simulations for the next current value.
As examples, computed CVCs of (self-sustaining) positive and negative coaxial corona discharges and a positive wireto-plane corona discharge in atmospheric-pressure air in the experimental conditions [13,14] are shown in figures 2 and 3, respectively. The kinetic scheme comprises an effective positive ion species, electrons, O − , O − 2 , and O − 3 . The kinetic and transport coefficients are the same as in [15] with a few minor changes, the photoionization is evaluated by means of the three-exponential Helmholtz model [16]. The coefficient of secondary electron emission from the negative electrode was set equal to 10 −4 .
Also shown in figures 2 and 3 are experimental data, which have been taken from figure 3 of [13] and figure 7 of [14], respectively, with values of the discharge current in the latter figure converted to the linear current density j L . Agreement between the modelling and the experiment is excellent for the coaxial corona and within approximately 7% for the wire-toplane corona.
It should be stressed that the resonance method was used for computation of only one point for each of the positive and negative coronas, shown in figure 2, and of one point for the positive corona, shown in figure 3: in each case it is the state with the lowest current density, which represents the corona ignition state. All the other states have been computed in a more or less standard way by means of a stationary solver, with account of all relevant nonlinear processes, including the distortion of the applied electric field by the plasma space charge, with gradual increase of the linear current density j L . Note that the possibility of using the current density j L (instead of the discharge voltage U) as a control parameter was achieved by using the 'weak form' formulation available in COMSOL, without expressly introducing an external circuit comprising a voltage source and a ballast resistance.
Note that values of the ignition voltage of the wire-toplane corona, computed on two different numerical meshes, one with 165 240 and the other with 534 184 degrees of freedom, with and without stabilization, were quite close: 4278.0 and 4278.2 V, respectively, with stabilization and 4282.2 and 4280.7 V without stabilization.

Discharge along a dielectric surface
Deposition of charged particles on a dielectric surface occurs on particle drift times and is considered as a fast process in the context of this work, so surface charges on a dielectric surface are considered as quasi-stationary. Therefore, the surface is under the floating potential and the boundary condition for the electric field is that its normal component is such that the normal component of the local current density is zero at each point of the surface.
The latter boundary condition is linear and homogenous with respect to the ion and electron densities. Therefore, while distributions of the ion and electron densities in a very lowcurrent self-sustaining discharge in the presence of dielectric vary proportionally to the discharge current I, the electric field distribution is independent of I. This is similar to what happens in the case where no dielectric surfaces border the active zone of the discharge, considered in preceding sections, and the independence of the electric field distribution on I is a necessary condition for the concept of ignition (inception) voltage to be meaningful. The difference 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 boundaryvalue problem involving partial differential and possibly also integro-differential equations, governing distributions of ions and electrons in the gap, and the Laplace equation governing the electrostatic potential, and the aim is to determine the value U = U 0 of applied voltage such that the problem admits a nontrivial solution for the particle densities. Note that this problem becomes ill-posed in the no-discharge situation, where the Discharge between two disc electrodes separated by cylindrical dielectric. Circles, crosses: ignition voltage of self-sustaining discharge, modelling with the boundary condition at the dielectric surface being zero normal current density or zero normal electric field, respectively. Triangle: breakdown voltage, experiment [17].
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 ignition voltage is computed firstly with the boundary condition at the dielectric surface being zero normal electric field, and at the second step the latter condition is replaced with the condition of zero normal current density. As an example, let us consider an axially symmetric discharge between two disc electrodes separated by a cylindrical dielectric surrounded by atmospheric-pressure air, which is shown in the insert in figure 4 and is relevant to mediumvoltage circuit breakers. (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 simulations have been performed for a range of values of the dielectric radius R, starting from 7.3 mm, which corresponds to the dielectric bordering the rounded edges of the metal discs, down to 3 mm. The kinetic model of atmospheric-pressure air was the same as in section 3.2, the coefficient of secondary electron emission from the negative electrode was set equal to 0.03.
The computed dependence of the ignition voltage on the dielectric radius is shown in figure 4. Circles and crosses represent results obtained with the boundary condition at the dielectric surface being zero normal current density or zero normal electric field, respectively; the difference is hardly visible. For R 5 mm, when the dielectric is away from the edges of the electrode discs, the ignition voltage is virtually constant and equals approximately 9.7 kV. 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 is attached to the edges, as seen in figure 4. As R increases, i.e., as the dielectric approaches the active zone of the discharge, the ignition voltage decreases. Leaving a detailed discussion outside the scope of this work, we only note that for R = 7.3 mm the discharge zone is adjacent to the dielectric; an understandable electrostatic effect.
The account of secondary electron emission from the dielectric surface with the secondary electron coefficient equal to 0.03 has not affected the ignition voltage to the graphic accuracy.
There seem to be no published measurements of the effect of dielectric surfaces on ignition voltages. However, one can perform a qualitative comparison with measurements of the breakdown voltage: although a precise relationship between the ignition voltage, which ensures inception conditions of a (self-sustaining) pre-breakdown discharge, and the breakdown voltage is unknown, one can expect that the ignition voltage is lower than the breakdown voltage and follows similar trends.
The triangle in figure 4 represents the experimental breakdown voltage for a dielectric with 3.5 mm radius reported in [17]. As expected, it is somewhat higher than the calculated ignition voltage. It is known from experiment that the presence of a dielectric facilitates the breakdown (e.g., review [18]), and it is indeed seen in figure 4 that the ignition voltage decreases as the dielectric approaches the active zone of the discharge.

Conclusions
From the mathematical perspective, the problem of ignition of self-sustaining gas discharges is an eigenvalue problem. Acknowledging and exploiting this fact is appropriate and beneficial. In this work, this eigenvalue problem is considered and an efficient physics-based numerical method for its solution is proposed. The method consists in the modelling of a non-self-sustained low-current discharge, generated by an external ionization source in the same plasma-producing gas in the same discharge configuration; the applied voltage is gradually increased until a resonance-like behaviour appears. The stationary boundary-value problem, describing the nonself-sustained discharge, is linear and may be routinely solved by means of ready-to-use solvers, including commercial ones.
Examples of solution results are shown for several typical discharge configurations: parallel-plate discharge, coaxial and wire-to-plane corona discharges, and a discharge along a dielectric surface. The results conform to the conventional Townsend breakdown condition for the parallel-plate configuration and are in good agreement with existing experimental data for the other configurations.
The method has the potential to provide a reference point for optimization of hold-off capabilities of high-voltage switchgear operating in low-frequency (e.g., 50 Hz) fields. Moreover, the method is useful for the initiation of the modelling of (already ignited) self-sustaining discharges by means of stationary solvers, as shown by the examples of coaxial and wire-to-plane corona discharges.