Bifurcations of current transfer through a collisional sheath with ionization and self-organization on glow cathodes

A bifurcation analysis is performed of a dc glow discharge between parallel electrodes and of a dc nearcathode space-charge sheath bordering a uniform plasma column. A model of plasma is considered with a single ion species and motion of the charged particles dominated by drift. Bifurcation points are found at which steady-state modes with spots on the cathode branch off from the abnormal mode or from the mode corresponding to the falling section of the current-density–voltage characteristic. In both discharge configurations, bifurcations in the abnormal mode have been detected; an unexpected result given that loss of stability and pattern appearance in dc gas discharges are usually associated with a negative differential resistance of the discharge. The conclusion is drawn that the two most important mechanisms governing appearance of patterns on glow cathodes, which are electrostatic mechanism and diffusion, produce competing effects: the former favors appearance of modes with multiple spots, while the latter favors appearance of a mode with one spot. This may explain the appearance in experiments of a normal spot or, alternatively, of patterns with multiple spots.


I. INTRODUCTION
The existence of multiple modes of current transfer to electrodes of gas discharges is rather the rule than the exception.As far as cathodes of glow discharges are concerned, it has been known for many decades that steady-state current transfer can occur in an abnormal mode or in a mode with a normal spot ͑e.g., Ref. ͓1͔; some further references can be found in Refs.͓2,3͔͒.Recently, also steady-state modes with patterns of more than one spot have been observed ͓4-9͔.Patterns of a number of spots have been observed also on anodes of dc glow discharges ͑e.g., Ref. ͓10͔ and references therein; some further references can be found in Ref. ͓11͔͒.A diffuse mode and a mode with one spot can occur on cathodes of high-pressure arc discharges; e.g., Ref. ͓12͔ and references therein.Patterns of several spots have been observed on cathodes of vacuum arcs; e.g., Refs.͓13-15͔.Diffuse, constricted, and multiple-spot modes can occur on anodes of high-pressure arc discharges; e.g., Refs.͓16-18͔.A wealth of different patterns has been observed in dc-driven planar discharge systems in which one metal electrode is separated from the gas by a high Ohmic barrier or is made of a high-Ohmic or semiconductor material, and in dielectric-barrier discharges; e.g., Ref. ͓19͔ and references therein and Refs.͓20-22͔ and references therein, respectively.Patterns of cathode spots have been observed in non-self-maintained discharges; e.g., Refs.͓23,24͔ and references therein.Structures have been observed also on cathodes of a pulsed glow discharge; e.g., Ref. ͓25͔.
Multiple modes of current transfer to electrodes of gas discharges are of considerable scientific interest, both intrinsically and as an important example of a self-organization phenomenon.The understanding of these modes is of critical importance for design of a number of technical devices.
Basically, there are two approaches to theoretical description of these modes.The first approach, which may be called phenomenological, is based on similarities between the pat-terns observed on electrodes of gas discharges and patterns appearing in nonlinear dissipative systems in other fields, such as biology, chemistry, optics, and semiconductor science.In the framework of this approach, the distribution of parameters along the electrode surface is assumed to be governed by a reaction-diffusion equation ͑which is a diffusion equation with a nonlinear source term͒ or by a system of coupled reaction-diffusion equations; e.g., Refs.͓10,26-29͔ and references therein.For special situations, attempts have been made to derive reaction-diffusion-type equations by means of a two-scale asymptotic technique from basic equations that govern a particular discharge, such as equations of conservation of the ions and the electrons and the Poisson equation; Refs.͓30-32͔.However, in most situations, reaction-diffusion equations for the distribution of parameters along the electrode surface are just postulated on the basis of the above-mentioned similarities.
The other approach to theoretical description of multiple modes of current transfer to electrodes of gas discharges is based on a direct numerical solution of the basic equations describing a particular discharge.There are works concerned with a two-dimensional numerical simulation of a normal spot on glow cathodes and of the simplest axially symmetric patterns ͑a solitary spot; a central spot surrounded by a ring spot͒ on glow anodes; e.g., Refs.͓3,33͔ and references therein and Ref. ͓11͔ and references therein, respectively.The transition from a Townsend discharge to a normal glow discharge was investigated in Ref. ͓34͔.Temporal and spatiotemporal patterns in a thin glow-discharge layer sandwiched with a semiconductor layer between planar electrodes have been studied in Refs.͓35,36͔ by means of nonstationary one-and two-dimensional numerical simulations, respectively.A three-dimensional self-organized pattern in a dielectric-barrier discharge was calculated by means of a direct numerical simulation of the discharge in Ref. ͓37͔.Steady-state modes with different spot patterns on cathodes of high-pressure arc discharges have been computed and their stability studied; e.g., Refs.͓38,39͔ and references therein.
Let us restrict further discussion to the case where all existing modes of current transfer to electrodes are stationary, setting aside non-dc discharges and dc discharges with nonstationary patterns or patterns that appear stationary but are accompanied by oscillations of electrical parameters of the discharge ͑e.g., spots in an obstructed planar glow discharge; see Refs.͓40,41͔͒.The ultimate goal of a theory in this case is the capability of predicting what modes are possible under particular discharge conditions and which of these modes are stable.In mathematical terms, the task is to find multiple steady-state solutions to a strongly nonlinear multidimensional boundary-value problem which describes a particular discharge, and to investigate the stability of each of these solutions.Obviously, this task cannot be solved without advanced numerical modeling.On the other hand, it can hardly be solved by means of numerical modeling alone: qualitative information on the character and the range of existence of each of the solutions is needed in order to ensure the convergence of iterations to a desired solution.͑In particular, one should know in advance whether the desired solution does exist at the specified value of discharge current; a question rather difficult to answer by purely numerical means.Indeed, if iterations have diverged or converged to another solution, one would hardly know whether this is a numerical problem or whether the solution being sought simply does not exist under the specified conditions.͒Such information can be obtained by resorting to ideas and methods of the theory of self-organization in nonlinear dissipative systems, as is done in the framework of the first approach described above to theoretical description of multiple modes of current transfer to electrodes of gas discharges.
At present, the above goal has been achieved only in the theory of cathodes of high-pressure dc arc discharges.The general pattern of different steady-state modes of current transfer to high-pressure arc cathodes turned out to be rather complex: a large number of modes with different spot patterns may exist at the same discharge current and more than one of these modes may be stable; Refs.͓38,39͔.The theory is much less advanced for multiple steady-state modes on cathodes of dc glow discharges.The only aspect that is developed relatively well is numerical simulation of a normal spot on glow cathodes; e.g., Refs.͓3,33͔ and references therein.As far as modes with multiple steady-state spots are concerned, opinions diverge even on such basic question as to what mechanisms are responsible for formation of multiple spots.The authors of Refs.͓4,6,9͔ believe that basic mechanisms of glow discharge are insufficient to describe modes with multiple spots and additional mechanisms are needed, such as an increasing dependence of the effective secondary emission coefficient on the reduced electric field, or heating of the gas.This point of view was disputed in Ref.
͓42͔ with a reference to two-decade-old works ͓43-48͔.On the other hand, no time-independent modes structured along the electrode surface have been obtained in the recent modeling of a glow discharge sandwiched with a semiconductor layer, which was reported in Ref. ͓36͔ and was performed with account of only basic glow-discharge mechanisms: either a stationary state homogeneous along the electrode sur-face was found, or an oscillating state, homogeneous or structured.
Bifurcation theory is a powerful tool of analysis of selforganization in nonlinear dissipative systems.It provides important qualitative information on patterns that are possible in the system considered, and this information, in addition to being useful by itself, facilitates finding the patterns numerically.Application of the bifurcation theory has triggered the development of a theory of multiple modes of current transfer to cathodes of high-pressure dc arc discharges; results of bifurcation analysis have guided numerical modeling ͓38,49,50͔ and have been of critical importance for understanding stability of different steady-state modes ͓39,51͔.
In the present work, a bifurcation analysis is performed of a fluid system of equations describing steady-state current transfer to cathodes of high-pressure glow discharges, in the spirit of the approach proposed in Ref. ͓47͔.Bifurcation points are sought at which stationary multidimensional solutions branch off from the one-dimensional ͑1D͒ stationary solution describing states homogeneous along the cathode surface.The results show that time-independent solutions structured along the cathode surface do exist, and allow one to explain, at least in principle, the appearance of different steady-state spot modes on dc glow cathodes, namely, of modes with multiple spots as reported in Refs.͓4-9͔, or of the normal mode which is observed in most cases.Questions left beyond the scope of the work include the nonlinear behavior of bifurcating stationary solutions ͑which can be found by running a 3D code in the vicinity of each bifurcation point with an initial approximation constructed using the corresponding eigenfunction found in the course of bifurcation analysis, similarly to how multidimensional solutions describing stationary spots on high-pressure arc cathodes have been found in Ref. ͓38͔͒; stability of stationary solutions; bifurcations of nonstationary solutions ͑Hopf bifurca-tions͒, similar to those found in Refs.͓35,36͔ for a glow discharge sandwiched with a semiconductor layer.
The outline of the paper is as follows.The model is described in Sec.II.In Sec.III, a procedure is developed for finding bifurcation points at which steady-state modes with spots branch off from the abnormal mode or the ͑unstable͒ mode corresponding to the falling section of the currentdensity-voltage characteristic ͑CDVC͒.Results of numerical calculations for high-pressure xenon plasma are reported in Sec.IV and placed into a more general context in Sec.V. Conclusions of the work are summarized in Sec.VI.

II. THE MODEL
The bifurcation approach used in this work is applicable to virtually any self-consistent model of a glow discharge, including models taking into account multiple ion and neutral species with a complex chemistry and a nonlocal electron energy distribution.However, as a first step it is advisable to consider the simplest self-consistent model, which is the one with a single ion species and motion of the ions and the electrons dominated by drift.Note that this model has already been employed for investigation of patterns in Refs.͓35,36,47͔.

A. The system of equations
Let us consider a high-pressure dc glow discharge in a discharge vessel in the form of a right cylinder, not necessarily circular; see Fig. 1.The origin of Cartesian coordinates x , y , z is on the cathode surface and the z axis is perpendicular to the cathode surface, i.e., parallel to the lateral surface of the discharge vessel.
At current densities of interest for this work, the bulk plasma is quasineutral; the space charge is localized in thin layers ͑sheaths͒ adjacent to the electrodes and to the lateral surface of the discharge vessel.Two geometries will be treated.In the first geometry, the calculation domain includes the whole of the discharge between parallel electrodes.In the other geometry, the calculation domain includes only the near-cathode space-charge sheath bordering a uniform plasma column.
The system of equations governing distributions in the plasma of the number densities of the ions and the electrons, n i and n e , and of the electrostatic potential includes equations of conservation of the ions and the electrons and the Poisson equation: where i and e are the mobilities of the ions and the electrons ͓known functions of the modulus of the local electric field: i = i ͑E͒, e = e ͑E͒, E = ٌ͉͉͔ and S is the net rate of production of the charged particles.
The equations of conservation of the ions and the electrons, Eqs.͑1͒ and ͑2͒, are written in the fluid approximation; the left-hand sides of these equations account for drift of the ions and the electrons but not for diffusion.In part, this is justified by the high values of the electric field: the voltage between the electrodes or the voltage drop across the sheath, being of the order of several hundred volts, is much higher than the ion and electron temperatures, which means that diffusion of the ions and the electrons in the axial direction is a minor effect.On the other hand, diffusion of the ions and the electrons in the lateral direction will be taken into account approximately, by means of a ͑nondifferential͒ term that will be included in the production rate S; cf.Eq. ͑8.20͒ of Ref. ͓1͔.
The net production rate represents the difference between the ionization rate and the rate of losses of charged particles due to diffusion to the lateral surface of the discharge vessel and dissociative recombination: Here i is the ionization frequency, dif is the effective frequency of diffusion losses ͑Ref.͓1͔, pp.194-195͒ and ␤ dis is the coefficient of dissociative recombination.The ionization frequency is related to Townsend's ionization coefficient ␣, i = ␣ e E, and will be considered as a known function of E.
In discharges in narrow tubes at low pressures and moderate currents, the charged particle density may be low enough for diffusion to be the dominating loss mechanism, dif ␤ dis n i .The third term on the right-hand side ͑RHS͒ of Eq. ͑4͒ can be omitted in such cases; this is the model of diffusion-controlled discharge.

B. Boundary conditions
Let us proceed to the formulation of boundary conditions for Eqs.͑1͒-͑3͒.Fluxes of the ions and electrons to the cathode surface are related through ␥, the effective secondary emission coefficient.This results in the boundary condition where ␤ = i / e .The potential of the cathode surface is known and may be set equal to zero, = 0. ͑6͒ There is no electric current to the ͑insulating͒ lateral surface of the discharge vessel, which results in the boundary condition Here n is a direction locally orthogonal to the lateral surface.
If the calculation domain includes the whole of the discharge between parallel electrodes, then boundary conditions must also be specified at the surface of the anode.Taking into account that there is no ion flux from the anode into the plasma and assuming that the potential of the anode surface is known, one obtains Here U is the discharge voltage, which will be considered as a given ͑control͒ parameter.Note that Eqs.͑1͒ and ͑2͒ may be rewritten as e n e ͑n e − n i ͒, ͑10͒ and may be viewed as differential equations of the first order for n i and n e , while Eq.͑3͒ may be viewed as a second-order differential equation of the elliptic type for .Each firstorder equation requires one boundary condition, and this condition cannot be specified on a lateral surface of the discharge vessel, which is parallel to ٌ.An elliptic equation requires one boundary condition at every boundary of the calculation domain.Thus, the above-specified set of boundary conditions ͑one condition at the lateral surface and two conditions at each of the electrodes͒ conforms to the type of system of differential equations being considered.
If the calculation domain includes only the near-cathode space-charge sheath, the two lacking boundary conditions must be specified on the plasma side of the sheath.Boundary conditions that allow one to describe a smooth ͑asymptotic͒ transition from the sheath to the plasma column are derived in Appendix A. These conditions apply at large distances from the cathode surface ͑much larger than the scale of thickness of the sheath͒ and read ͑11͒ Here E ϱ is the electric field in the plasma column and U is the sheath voltage.U will be considered as a given parameter and E ϱ should be found as a part of the solution.

C. Transforming the equations
It is possible to eliminate the unknowns n i and n e from the system of governing equations ͑1͒-͑3͒, thus transforming the system to a single differential equation for .To this end, let us first subtract Eq. ͑1͒ from Eq. ͑2͒, arriving thus at the equation of continuity of electric current: Eliminating from this equation n i with the use of Eq. ͑3͒, one can rewrite it in the form Equation ͑2͒ may be rewritten as Eliminating n i from the LHS of Eq. ͑1͒, one can rewrite this equation as

͑15͒
Multiplying Eq. ͑14͒ by ␤ and adding Eq. ͑15͒, one obtains after transformations where K = ␣E − dif / e .This equation does not involve ٌn e and may be considered as an algebraic ͑quadratic͒ equation for n e .Solving this equation, one will obtain an explicit expression for n e in terms of .Substituting this expression into Eq.͑13͒, one will arrive at the desired equation containing only .However, an explicit form of this equation is unnecessary for the purposes of the present work, and Eqs.͑13͒ and ͑17͒ suffice.

A. General
If the current in a dc glow discharge is high enough, then it is distributed over the cathode surface more or less uniformly; the abnormal mode.As the current decreases, only part of the cathode surface remains active; a spot mode.͑A spot mode here and below means a normal mode or any of the modes with more than one spot observed in Refs.͓4-9͔.͒Both the abnormal mode and spot mode are stationary.
The nonlinear boundary-value problem formulated in the preceding section admits a 1D solution, F = F͑z͒.This solution is well known and describes steady-state modes that are homogeneous along the cathode surface: the abnormal mode; the mode corresponding to the falling section of the CDVC U͑j͒, which is unstable and is not realized; and the Townsend discharge, which represents the limiting case of very low current densities.͑In the following, these three modes will be jointly referred to as the distributed mode.͒The question is whether the boundary-value problem admits also solutions describing steady-state modes with spots.If such solutions exist, they are necessarily multidimensional, F = F͑x , y , z͒.While the 1D solution exists at all current values, although at low currents it is unstable and is not realized, multidimensional solutions should exist at low currents but may cease to exist at currents high enough.
It happens frequently in problems with multiple solutions of different symmetries that different solutions coincide at certain values of control parameters.In the problem considered, this means that the stationary multidimensional solutions may exactly coincide with the stationary 1D solution at certain values of the discharge current ͑or the voltage U͒.One can say that multidimensional solutions branch off from ͑or join͒ the 1D solution, and values of the discharge current ͑or of the voltage U͒ at which this occurs may be called bifurcation points.In terms of stability theory, these bifurcation points represent states at which the 1D solution becomes neutrally stable against nonoscillatory multidimensional perturbations.
The present section is concerned with finding these bifurcation points or, in other words, with finding values of the discharge current or of the voltage U at which modes with stationary spots branch off from the distributed mode.The analysis is performed in the framework of the boundaryvalue problem formulated in Sec.II.
The simplest bifurcations of stationary solutions governed by a single control parameter are of the following types ͑e.g., Refs.͓52-54͔͒: saddle-node ͑fold͒ bifurcations, transcritical bifurcations, and pitchfork bifurcations.A saddle-node bifurcation occurs when a solution reaches a limit of its existence region and turns back.Thus saddle-node bifurcations are unrelated to branching of essentially different solutions, and the latter can occur only through transcritical or pitchfork bifurcations.If a bifurcation is transcritical, one could think of a quasistationary transition from one stationary solution to another.The same applies to a pitchfork bifurcation if it is supercritical.There can be no quasistationary transition between stationary solutions if a pitchfork bifurcation is subcritical.
A discussion of these bifurcations given in Refs.͓49,51͔ in the context of steady-state current transfer to cathodes of dc arc discharges is expected to largely apply to the glow cathodes being treated in the present work.The saddle-node bifurcations occur at extreme points of the CDVC of the distributed mode, at which the 1D ͑distributed͒ mode reaches a limit of its existence region and turns back.Branching of multidimensional ͑spot͒ modes occurs at bifurcation points positioned outside extreme points of the CDVC of the 1D mode, and the corresponding bifurcations are transcritical or pitchfork.The analysis to be developed in the present work will allow one to find all bifurcation points at which steadystate multidimensional modes branch off from the steadystate 1D mode, whether the bifurcation is transcritical, or supercritical pitchfork, or subcritical pitchfork.Also, one will know what kind of bifurcation occurs at each point, transcritical or pitchfork.On the other hand, the present analysis does not include investigation of nonlinear behavior of bifurcating multidimensional solutions; therefore in the case of a pitchfork bifurcation it does not provide means for distinguishing between supercritical and subcritical bifurcations.

B. Formulating the eigenvalue problem
Let us designate by U 0 the value of the voltage which corresponds to one of the bifurcation points and by n i0 ͑z͒, n e0 ͑z͒, and 0 ͑z͒ the corresponding stationary solution ͑which is 1D͒.Stationary solutions in the vicinity of this point, i.e., at where n i1 , n e1 , and 1 are small perturbations of the solution at the bifurcation point, Substituting ͑18͒ into the governing equations, expanding in n i1 , n e1 , and 1 , and retaining terms of the zeroth order, one arrives at a problem governing the 1D solution.This problem can be easily solved, for example, as described in Appendix B. In the following, the functions n i0 ͑z͒, n e0 ͑z͒, and 0 ͑z͒ will be considered known.
Equations for perturbations are obtained by retaining terms of the first order.It is convenient to derive these equations by expanding Eqs.͑13͒ and ͑17͒: where the prime denotes d / dz, E 0 = 0 Ј, j 0 is the value of current density corresponding to U = U 0 , and ͑22͒ The ion and electron mobilities are not as strong functions of the electric field as is the ionization frequency; in this section i and e for brevity are treated as constant.In the derivation of Eq. ͑20͒, a two-term expansion of E in 1 was used, which reads as follows: Applying a similar procedure to the boundary conditions ͑5͒ ͑with n i eliminated in terms of n e and ͒, ͑6͒, and ͑7͒, one arrives at the following boundary conditions for the perturbations at the cathode surface: and at the lateral surface: If the calculation domain includes the whole of the discharge between parallel electrodes, then boundary conditions for the perturbations apply at the surface of the anode and are obtained from Eq. ͑8͒ with n i eliminated in terms of n e and : If the calculation domain includes only the near-cathode space-charge sheath, then boundary conditions for the perturbations apply at large distances from the cathode surface and are obtained from Eq. ͑11͒: Here E ϱ1 is a perturbation of the electric field in the plasma column, which should be found as a part of the solution.
Since the electric field in the plasma column is spatially uniform, E ϱ1 is constant.Now a problem governing perturbations has been formulated.Note that it involves only 1 and n e1 , i.e., it does not involve n i1 .We recall that the functions with the index 0, being governed by the 1D problem, may be considered known while the problem for perturbations is dealt with.As should have been expected, the problem for perturbations is linear and inhomogeneous, the inhomogeneity being represented by the term with U 1 .
Since the point U = U 0 under consideration is a bifurcation point, more than one steady-state solution exists in the vicinity of this point, i.e., at U 1 0. Hence, the problem governing perturbations has a nonunique solution for U 1 0. The latter means that the corresponding homogeneous problem, i.e., the one with U 1 = 0, has a nontrivial solution.The form of the homogeneous problem allows one to apply the method of separation of variables, similarly to how it is done in, e.g., many eigenvalue problems for the Schrödinger equation ͑e.g., ͓55͔͒, and to seek this nontrivial solution as n e1 ͑x,y,z͒ = g͑z͒⌽͑x,y͒, 1 ͑x,y,z͒ = f͑z͒⌽͑x,y͒.͑28͒ Substituting these expressions into Eqs.͑19͒ and ͑20͒ and dividing by ⌽, one arrives at

͑31͒
The independent variables in Eqs.͑29͒ and ͑30͒ are separated and the variables x and y may be involved only through the ratio ٌ Ќ 2 ⌽ / ⌽.Hence, this ratio is constant.Designating this constant by −k 2 , one can write an equation for the function ⌽͑x , y͒ in the form This equation must be solved in the cross section of the discharge vessel.Note that k has the meaning of a wave number characterizing variation of the bifurcating multidimensional solution in the plane ͑x , y͒.A boundary condition at the lateral surface of the vessel follows from Eq. ͑25͒: ‫ץ‬⌽ ‫ץ‬n = 0. ͑33͒ Equations ͑32͒ and ͑33͒ represent the Neumann eigenvalue problem for the Helmholtz equation.It is known ͓56͔ that all eigenvalues of this problem are real and of finite degree of degeneracy; a set of eigenvalues is countable and may be numbered in order of their increase with the eigenvalues growing without limit with the increase of the number.For example, in the case of a discharge vessel with a circular cross section, this problem can be solved in terms of the Bessel functions and its spectrum is k s = j ,s Ј / R, where =0,1,2,..., s =1,2,3,..., R is the radius of the tube, and j ,s Ј is the sth zero of the derivative of the Bessel function of the first kind of order ͑the initial four nontrivial zeros are j 1,1 Ј Ϸ 1.84, j 2,1 Ј Ϸ 3.05, j 0,2 Ј Ϸ 3.83, and j 3,1 Ј Ϸ 4.20͒.A graphic interpretation of the function ⌽͑x , y͒ for this case can be found, e.g., in ͓42͔.Solving Eq. ͑30͒ for g and substituting the result into Eq.͑29͒, one arrives at a fourth-order ordinary differential equation for the function f: where = E 0 2 ͑fЉ − k 2 f͒ and

͑35͒
Boundary conditions for Eq.͑34͒ at the cathode surface are

͑36͒
If the calculation domain includes the whole of the discharge between parallel electrodes, then further boundary conditions for the function f apply at the surface of the anode and are obtained from Eq. ͑26͒ with U 1 =0:

͑37͒
If the calculation domain includes only the near-cathode space-charge sheath, then the boundary conditions ͑27͒ must be employed.The LHS of the first boundary condition ͑27͒ is proportional to ⌽͑x , y͒; the quantity E ϱ1 is constant.This condition may be satisfied only with E ϱ1 = 0, in which case this condition assumes the form fЈ → 0. Given that U 1 = 0, the second boundary condition ͑27͒ assumes the form f → 0. Thus, a boundary condition for Eq.͑34͒ at large distances from the cathode reads The asymptotic behavior of a general solution to Eq. ͑34͒ at large z includes two exponentially decaying terms, an exponentially growing term, and a term that is either exponentially growing or constant; see Appendix C. In order that this asymptotic behavior be compatible with the boundary condition ͑38͒, exponentially growing and constant terms must be excluded.Hence, the boundary condition ͑38͒ allows one to determine two integration constants and the problem ͑34͒, ͑36͒, and ͑38͒ is complete.

C. Procedure for finding bifurcation points
In order to find bifurcation points, one should first determine the spectrum of the eigenvalue problem ͑32͒ and ͑33͒.The trivial eigenvalue that corresponds to a constant eigenfunction is discarded.The first nontrivial eigenvalue k = k 1 is substituted into the problem constituted by Eqs.͑34͒, ͑36͒, and ͑37͒ or ͑38͒.The problem obtained is solved for different j 0 until values of j 0 are found for which a nontrivial solution exists.͓In other words, the problem constituted by Eqs.͑34͒, ͑36͒, and ͑37͒ or ͑38͒ is solved for a given k = k 1 as an eigenvalue problem, j 0 being the eigenparameter.͔Value͑s͒ of j 0 determined in this way define͑s͒ bifurcation point͑s͒ corresponding to k = k 1 .Multidimensional solutions that branch off from ͑or join͒ the 1D solution at these points are described in the vicinity of these points by eigenfunction͑s͒ ⌽͑x , y͒ of the problem ͑32͒ and ͑33͒ associated with the eigenvalue k = k 1 .After that, the second eigenvalue k = k 2 is substituted into the problem constituted by Eqs.͑34͒, ͑36͒, and ͑37͒ or ͑38͒ and corresponding bifurcation points are found.Multidimensional solutions that branch off from ͑join͒ the 1D solution at these points are described in the vicinity of these points by eigenfunction͑s͒ of the problem ͑32͒ and ͑33͒ associated with the eigenvalue k = k 2 .Then the third value k = k 3 is substituted, etc.This procedure includes solving three problems: ͑1͒ The 1D problem governing the functions n i0 ͑z͒, n e0 ͑z͒, and 0 ͑z͒ ͑in the following, the index 0 attributed to the 1D solution will be dropped for brevity͒; ͑2͒ the eigenvalue problem ͑32͒ and ͑33͒; and ͑3͒ the eigenvalue problem constituted by Eqs.͑34͒, ͑36͒, and ͑37͒ or ͑38͒.The 1D problem may be conveniently solved as described in Appendix B. Note that while dealing with the 1D solution it is convenient to replace U 0 as a given parameter with j 0 the density of electric current from the plasma to the cathode.Then the 1D problem can be reduced to a boundary-value problem for a second-order ordinary differential equation, if the calculation domain includes the whole discharge from the cathode to the anode, and to an initial-value problem for a first-order ordinary differential equation, if the calculation domain includes only the nearcathode space-charge sheath.
A solution to the Neumann eigenvalue problem for the 2D Helmholtz equation, represented by Eqs.͑32͒ and ͑33͒, is known for discharge tubes of simple shapes, for example, in the above-mentioned case of a tube with a circular cross section.
The eigenvalue problem constituted by Eqs.͑34͒, ͑36͒, and ͑37͒ or ͑38͒ must be solved numerically.In this work, it has been solved with the use of the Petukhov method ͓57͔ ͑which is a method of numerical solution of boundary-value problems for linear ordinary differential equations of the second or third order, or for partial differential equations of parabolic type, based on a finite-difference scheme of the fourth order of accuracy͒ generalized for the case of a system of equations.

IV. COMPUTATION RESULTS
In this section, results are considered of calculations of bifurcation points at which modes with stationary spots branch off from the distributed mode.The calculations have been performed by means of the procedure described in Sec.III and refer to the conditions of experiments with microdischarges in high-pressure xenon plasma reported in Refs.͓4-6,9͔.The mobility of ions Xe 2 + in Xe was evaluated by means of the formula i = 2.1ϫ 10 21 m −1 V −1 s −1 / n ͑here n is the density of the neutral gas͒, which is an approximation of the measurements of Ref. ͓58͔.The mobility of the electrons and Townsend's ionization coefficient were calculated using the zero-dimensional Boltzmann solver BOLSIG ͓59͔.The calculation results for mobility were approximated by the formula e = 19 Torr m 2 V −1 s −1 / p, where p is the pressure of the plasma.The calculation results for Townsend's ionization coefficient are well approximated by Eq. ͑4.6͒ of Ref. ͓1͔ in the range of reduced electric fields of interest.The coefficient of dissociative recombination of molecular ions Xe 2 + was set equal to 2 ϫ 10 −13 m 3 / s, which is the value given in Refs.͓60,61͔ for the neutral gas temperature T 0 = 300 K and the electron temperature T e = 1 eV.The effective frequency of diffusion losses, dif , was estimated for a cylindrical discharge vessel by means of Eq. ͑8.18͒ of Ref.
͓1͔: dif = D a ͑2.4/R͒ 2 , where D a = T e i is the coefficient of ambipolar diffusion ͑which was evaluated at T e =1 eV͒.The effective secondary emission coefficient ␥ was set equal to 0.03.
Calculation results given in this work for a discharge between parallel electrodes refer to the interelectrode gap h = 0.5 mm.Diffusion losses have been found to produce a very small effect on the computation results for realistic values of the discharge vessel radius ͑down to R = 0.25 mm͒; we indicate for definiteness that the results presented in this work for a parallel-plane discharge have been calculated without diffusion losses.
In Fig. 2, values of the current density are shown at which a steady-state multidimensional solution F = F͑x , y , z͒ with a given k branches off from the 1D solution F = F͑z͒, or, in other words, at which a spot mode with a given k branches off from the distributed mode.We recall that k has the meaning of a wave number characterizing variation of the bifurcating spot mode along the cathode surface.The dashed lines in Fig. 2 represent values j min of the current density at which the discharge voltage in the distributed mode attains a minimum value.
Note that data shown in Fig. 2 are suitable for determination of bifurcation points for any geometry of the discharge vessel cross section: after the eigenvalue problem Eqs.͑32͒ and ͑33͒ for a particular cross section has been solved and its spectrum k 1 ,k 2 ,k 3 ,... found, the graphs shown in Fig. 2 allow one to identify values of j at which the first, second, third, etc. spot modes branch off from the distributed mode.For example, the lowest nontrivial value of k for a circular cathode of the radius R = 0.375 mm is k 1 Ϸ 1.84/ R Ϸ 4.9 ϫ 10 3 m −1 ; if the pressure is between 75 and 760 Torr, one can see from Fig. 2 that the first spot mode, which is the one with a spot at the edge of the cathode, branches off from the distributed mode very close to the point of the minimum of the CDVC.
There are two disconnected families of bifurcation points at p = 75 Torr ͓see Fig. 2͑a͔͒, one of these families being represented by full circles and the other by open circles.As the plasma pressure increases, the families come closer and become connected at p around 100 Torr; see Fig. 2͑b͒.As the plasma pressure increases further, two disconnected families appear again ͓see Fig. 2͑c͔͒; however, the topology of these families is different from that seen in Fig. 2͑a͒: one can say that the families have exchanged branches.Only one family can be seen on the graph for atmospheric pressure, Fig. 2͑d͒.
It is useful to normalize the wave number in order to better understand the computation results.It is convenient to use as a normalization multiplier the ratio of the discharge voltage to the electric field E w at the cathode surface, ⌬ = U / E w , evaluated for the distributed mode.At very low current densities ͑the limiting case of Townsend discharge͒, the electric field in the gap is not appreciably distorted and ⌬ tends to the gap width, ⌬ Ϸ h.At high current densities, where a dominating contribution to the discharge voltage is given by a thin space-charge sheath near the cathode, ⌬ can be viewed as the scale of thickness of the sheath.This is illustrated by Fig. 3, in which distributions of parameters across the discharge, obtained by numerical calculations, are depicted for p = 75 Torr and different values of the current density.Note that the value j = 1.51ϫ 10 3 A / m 2 , to which Fig. 3͑c͒ refers, corresponds to the point of the minimum of the CDVC.The arrows below the abscissa axis represent values of ⌬.One can see that the numerical results conform to the above-described physical meaning of the length scale ⌬.
Normalized values of the wave number of the spot mode as a function of the current density at which this mode branches off from the distributed mode are shown in Fig. 4. Also shown in this figure are the CDVCs of the distributed mode.Horizontal dashed lines in Fig. 4 represent the values of the discharge voltage corresponding to the limiting case of the Townsend discharge, which have been calculated with the use of Eq. ͑B22͒ of Appendix B. As expected, these lines represent asymptotes of the CDVCs at low j.For convenience, values of the current density j min at which the discharge voltage in the distributed mode attains a minimum value are depicted in Fig. 4 by the vertical dashed lines.
One of the conditions of applicability of the model of drift-dominated motion of the ions and electrons is that the characteristic time of drift of the electrons in the axial direction be much smaller than the time of diffusion of the electrons in the direction along the cathode.The former time is of the order of ⌬ 2 / e U; the latter is of the order of 1 / D e k 2 , where D e = T e e is the electron diffusion coefficient.The above validity condition amounts to k⌬Ӷ ͱ U / T e ϳ 10. ͑A similar condition should be imposed on the ions; however, it not more restrictive than that for the electrons given that the ion temperature does not exceed T e .͒This is why the range of normalized wave numbers shown in Fig. 4 is limited from above by 10.Only one family of bifurcation points exists in this range; the other family is positioned beyond the region of applicability of the model and therefore is not of particular interest.With account of this limitation, the dependence of k⌬ on j does not change dramatically with the plasma pressure, in contrast to what is seen in Fig. 2, and this is why data Any extreme point of a CDVC of the distributed mode represents a bifurcation point associated with k =0.͑This is a saddle-node bifurcation; see the discussion in Ref. ͓51͔.͒This explains why bifurcation points shown in Fig. 4 in the limiting case of small wave numbers ͑or long wavelengths͒, k⌬Ӷ1, approach the line j = j min , i.e., are positioned in the vicinity of the point of minimum of the CDVC.As k⌬ increases, the bifurcation points are shifted along the rising section of the CDVC in the direction of increasing j.As k⌬ grows further, the bifurcation points reach the current density value j = j A and then start moving in the direction of decreasing j; here A is the rightmost point of the curve k⌬͑j͒ as shown in Fig. 4. On reaching the current density value j = j B , the bifurcation points change the direction of their displacement once again.
The above-described bifurcations occur in the range of high enough current densities where the discharge possesses a well-developed structure composed of a ͑uniform͒ plasma column and near-electrode sheaths; see Figs. 3͑c͒ and 3͑d͒.A question arises: are the bifurcations generated by the nearcathode sheath, as one might think intuitively?And if so, do the bifurcations in a 0.5-mm-wide gap occur in the same way as in a near-cathode sheath bordering an ͑infinitely long͒ uniform plasma column?
In order to answer these questions, let us turn to the second model described in Sec.II, namely, to the model in which the calculation domain includes only the near-cathode space-charge sheath bordering the uniform plasma column.The dotted lines in Fig. 4 represent CDVCs of the distributed mode calculated by means of this model.͑Note that the quantity U has different meanings in the two models, the discharge voltage in the framework of the first model and the sheath voltage in the framework of the second model.Therefore, the product E ϱ h accounting for the voltage drop in the plasma column was added to the sheath voltage before plotting the dotted curves in Fig. 4.͒ One can see that the model of a near-cathode space-charge sheath bordering the uniform plasma column describes the distributed mode with good accuracy in the range of j of interest.The question is whether it can provide a good accuracy also for bifurcation points.
Diffusion losses have been found to produce a small effect for the conditions of interest in the framework of the sheath model, similarly to what has been found for the model of a parallel-plane discharge.We indicate for definiteness that the results given in this work for the sheath model have been calculated with the diffusion losses evaluated for the discharge vessel radius of 0.4 mm.
Bifurcation points in the plane ͑k⌬ , j͒, calculated by means of the sheath model, are shown in Fig. 5. Also shown are the CDVCs of the distributed mode, U͑j͒ ͑solid lines; we recall that in the framework of the sheath model U designates the sheath voltage͒.The values j min of the current density at which the sheath voltage in the distributed mode attains a minimum value are depicted in Fig. 5 by the vertical dashed lines.As in the case of the model of a discharge between parallel electrodes, values of k⌬ of the order of 10 or higher are beyond the range of applicability of the assumption of motion of charged particles dominated by drift; nevertheless, data in the range of k⌬ up to 10 2 are shown in Fig. 5 in order to give an idea of the topology of bifurcation points.
One can see two disconnected families of bifurcation points in each of Figs.5͑a͒ and 5͑b͒.One of these families is localized in the range k⌬Շ0.1 and represented by full circles; the second family is localized in the range k⌬տ0.1 and represented by open circles.At small values of the normalized wave number, k⌬Շ10 −2 , the first family approaches the dashed line j = j min , i.e., the bifurcation points are positioned in the vicinity of the point of minimum of the CDVC of the distributed mode.As k⌬ increases, the bifurcation points of the first family are shifted along the falling section of the CDVC of the distributed mode in the direction of decreasing j, and disappear in the region of very low current densities.This behavior is different from that discussed above for the model of parallel-plane discharge.Bifurcation points belonging to the second family emerge from the region of very small j.As k increases and the product k⌬ increases as well ͓at p = 75 Torr, the product k⌬ first slightly decreases and then starts to increase; see Fig. 5͑a͔͒, the bifurcation points of the second family are shifted in the direction of increasing j and eventually attain the region j Ͼ j min , i.e., the growing section of the CDVC.As k⌬ grows further, the bifurcation points change the direction of their displacement and start shifting in the direction of decreasing j.Then they leave the growing section of the CDVC and eventually once again disappear in the region of very low j.
Let us proceed to comparison of results for bifurcation points given by the model of parallel-plane discharge and by the sheath model.This comparison is meaningful if the length of attenuation of perturbations in the plasma column ͑i.e., the length of decay of a solution to the eigenvalue problem at large z͒ in the sheath model does not exceed by order of magnitude the gap width in the model of parallel-plane discharge.At small k, the length of attenuation of perturbations in the plasma column equals ͉ 2 ͉ −1 , where 2 is given by the second equation in Eq. ͑C8͒ of Appendix C, and is of the order of k −1 .Hence, the above comparison is meaningful provided that k տ h −1 .In addition, the current density must be high enough so that the thickness of the near-electrode  space-charge sheaths is much smaller than the gap width, which is the basic condition of applicability of the sheath model.Taking into account the comparison between CDVCs predicted by the two models, which is shown in Fig. 4, one can assume that j must exceed 10 3 -10 4 A / m 2 .
Bifurcation points in the plane ͑k , j͒ predicted by the two models in the domain k Ն 10 3 m −1 , j Ն 10 3 A / m 2 are shown in Fig. 6.The topology of the dependence j͑k͒ predicted by the model of parallel-plane discharge is more involved than that predicted by the sheath model; however, those sections of the dependence that are common for the two models are rather close to each other.The latter conforms to the assumption that the bifurcations are generated by the near-cathode sheath.On the other hand, bifurcations in a 0.5-mm-wide gap occur in not quite the same way as in a near-cathode sheath bordering an ͑infinitely long͒ uniform plasma column.

V. DISCUSSION
Information on the relationship between the wave number of a spot mode and the value of the current density at which this mode branches off from the distributed mode allows one to hypothesize about the spot patterns that can occur in a particular discharge.Let us consider as examples the depen-dences k͑j͒ schematically shown in Fig. 7.In contrast to Figs. 2 and 6, these dependences are shown in Fig. 7 by lines rather than by points; ͑solid͒ lines I and II.j min in Fig. 7 designates the value of the current density at which the discharge ͑or sheath͒ voltage in the distributed mode attains a minimum value; k 1 , k 2 , k 3 , and k 4 are the four smallest nontrivial eigenvalues of the problem ͑32͒ and ͑33͒.Let us for definiteness assume once again that the discharge tube has a circular cross section; then spot modes that branch off at bifurcation points associated with eigenvalues k 1 , k 2 , k 3 , and k 4 are those with, respectively, a spot at the edge of the cathode; two spots at the edge; a spot at the center of the cathode or a ring spot at the edge; and three spots at the edge.The modes with one, two, or three spots at the edge are 3D and branch off from the distributed mode through a pitchfork bifurcation, supercritical or subcritical, while the mode with a spot at the center of the cathode or a ring spot at the edge is axially symmetric and branches off through a transcritical bifurcation; see Ref. ͓49͔.j 1 ͑I͒ -j 4 ͑I͒ and j 1 ͑II͒ -j 4 ͑II͒ in Fig. 7 designate values of the current density at which the first to fourth bifurcations occur in the case I or, respectively, II.Both lines I and II start at the point ͑j = j min , k =0͒; a consequence of the above-mentioned fact that any extreme point of the CDVC of a distributed mode represents a bifurcation point associated with k =0.
One can expect that the transition from the distributed mode to a mode with one or more spots occurs as follows.The distributed mode is stable beyond all bifurcation points, i.e., at j high enough.As j decreases, the distributed mode loses stability as the first bifurcation point is encountered and the discharge switches into the spot mode that branches off at this bifurcation point.In the example depicted by line I, bifurcations occur on the falling section of the CDVC.The first bifurcation point to be encountered as j decreases from high values is j = j 1 ͑I͒ .Thus, the distributed mode is stable on the growing section of the CDVC, j Ն j min , and on an initial segment of the falling section, j 1 ͑I͒ Ͻ j Յ j min .At j = j 1 ͑I͒ , the distributed mode becomes unstable and the discharge switches into a mode with a spot at the edge of the cathode.Whether this switching can be realized in a quasistationary way or is necessarily accompanied by a jump of electric parameters depends on whether the bifurcation at j = j 1 ͑I͒ ͑which is pitchfork, as indicated above͒ is super-or, respectively, subcritical.
In the example depicted by line II, the first bifurcation to be encountered at decreasing j occurs at j = j 4 ͑II͒ , i.e., on the growing section of the CDVC.The distributed mode is stable at j Ͼ j 4 ͑II͒ ; at j = j 4 ͑II͒ it loses stability and the discharge switches into a mode with three spots at the edge of the cathode.Again, this switching can be quasistationary or is accompanied by a jump, depending on whether the ͑pitch-fork͒ bifurcation at j = j 4 ͑II͒ is super-or, respectively, subcritical.
These examples show that the distributed mode gives way with decreasing j to a mode with one spot, if all bifurcation points are positioned on the falling section of the CDVC, or to a mode with several spots, if there are bifurcation points on the growing section.According to calculations of this work, both the model of parallel-plane discharge and the sheath model predict the presence of bifurcation points on the growing section of the CDVC and, consequently, a transition from the abnormal mode to a mode with multiple spots.The fact that this happens on the growing section of the CDVC is somewhat unexpected, given that loss of stability of a homogenous stationary state and pattern appearance in gas discharges are usually associated with a negative differential resistance of the discharge.
One could think of an analogy with spot modes that were recently observed on the left-hand side of the Paschen curve, where the CDVC is monotonically growing; see Ref. ͓40͔.However, a subsequent investigation reported in Ref. ͓41͔ has revealed that the discharge is essentially nonstationary in such conditions.It is unclear therefore to what extent the above-mentioned analogy would be meaningful.One could mention also the work ͓62͔, in which the role of negative differential conductivity ͑or, equivalently, negative differential resistance͒ of the glow discharge was studied by means of a nonstationary one-dimensional numerical simulation.However, the results of this work are not quite clear: on one hand, no numerical counterexamples have been found where oscillations would occur while the CDVC of the discharge shows a positive differential conductivity; on the other hand, the authors reported, without giving details, the existence of a counterexample where a state with positive differential resistance is unstable and develops into a limit cycle.
The analysis of this work is based on Eqs.͑1͒-͑3͒ or, equivalently, on Eqs.͑3͒, ͑9͒, and ͑10͒.The only term with a second derivative in the latter set of equations is the one on the LHS of the Poisson equation ͑3͒.One can say that the only way in which spots may "interact" each with other and with the walls of the discharge tube is through electrostatic force.Hence, the only mechanism governing patterns in the framework of the present theory is electrostatic.According to results of this work, the electrostatic mechanism favors appearance of modes with multiple spots.
Other mechanisms are also possible, the most important being diffusion of the charged particles.One can expect that diffusion tends to smooth out perturbations, i.e., produces a stabilizing effect over the distributed mode, and this effect becomes stronger with decreasing wavelength of perturbations, i.e., with increasing k.Therefore, diffusion causes a shift of the bifurcation points, which represent points of neutral stability of the distributed mode, in the direction of lower j, and this shift increases with increasing k.In a situation where diffusion is the only mechanism governing patterns, one can expect that the bifurcation curve will be similar to the one depicted by line I in Fig. 7. ͓In fact, one can expect to encounter this kind of bifurcation curve in any situation where patterns are governed by dissipative mechanism͑s͒; for example, a bifurcation curve of this kind appears in the problem of current transfer to hot arc cathodes, where patterns are governed by thermal conduction in the cathode body.͔One can say that diffusion favors appearance of a mode with one spot, i.e., of the normal mode.
Thus, the two most important mechanisms governing the appearance of patterns on glow cathodes, the electrostatic mechanism and diffusion, produce competing effects: while the former favors appearance of modes with multiple spots, the latter favors appearance of a mode with one spot.This allows one to explain, at least in principle, different patterns that occur on glow cathodes: if the electrostatic mechanism prevails, then the abnormal mode gives way at decreasing current to multiple-spot patterns reported in Refs.͓4-9͔; if diffusion prevails, then the abnormal mode gives way at decreasing current to a mode with one spot, i.e., to the normal mode, which is observed in most cases.
The above hypothesis should be validated by a multidimensional simulation of spot modes in the whole range of their existence and of the stability of these modes.A bifurcation analysis is capable of providing valuable guidance in this work, as it has done in the theory of near-cathode phenomena in high-pressure arc discharges ͓38,39,49-51͔.Note that a solution describing the distributed mode in a discharge tube with a circular cross section, obtained with account of diffusion, is axially symmetric due to diffusion of the charged particles to the walls of the tube.Then one needs to deal with branching of 3D solutions from an axially symmetric solution, instead of branching of axially symmetric and 3D solutions from a 1D solution, which is dealt with in the present work.This can be done by means of a procedure similar to the one developed in Ref. ͓50͔.In addition to diffusion, the model of glow discharge used in the present work disregards a number of other effects that may play a role, such as the presence of multiple ion and/or neutral species and variations of the electron and heavy-particle temperatures.However, a similar approach based on the bifurcation theory may be applied also in the framework of more complex models, such as, e.g., the model of microhollow cathode discharges in xenon developed in Ref. ͓63͔.

VI. CONCLUSIONS
Bifurcation points have been calculated at which steadystate spot modes on dc glow cathodes branch off from the abnormal mode or the mode corresponding to the falling section of the CDVC.The simplest self-consistent model of dc glow discharge is employed, namely, the one in which the plasma contains a single ion species and motion of the ions and the electrons is dominated by drift.For both glowdischarge configurations treated, which are a discharge between parallel electrodes and a space-charge sheath separating a cathode from a uniform plasma column, bifurcations in the abnormal mode have been detected; an unexpected result given that loss of stability and pattern appearance in dc gas discharges are usually associated with a negative differential resistance of the discharge.
This result may be viewed as an indication that the electrostatic mechanism, which is the only mechanism governing patterns in the framework of the present theory, favors appearance of modes with multiple spots.Diffusion of the charged particles, which is another potentially important mechanism, favors appearance of modes with one spot.A competition of these two mechanisms can explain why the abnormal mode at decreasing current densities in certain cases gives way to modes with multiple spots and in other cases to the normal mode, i.e., to a mode with one spot.
The bifurcation approach developed in this work can be applied also to more complex models of glow discharges and is capable of providing valuable guidance for multidimensional simulation of spot modes in the whole range of their existence and of stability of these modes, as it has done in the theory of near-cathode phenomena in high-pressure arc discharges.

APPENDIX A: ASYMPTOTIC BOUNDARY CONDITIONS FOR THE SHEATH EQUATIONS
Equations ͑1͒-͑3͒ allow one to describe a smooth ͑asymptotic͒ transition from the sheath to the plasma column.In order to exploit this opportunity, let us first discuss properties of the plasma column described by these equations.
As usual, the electric field in the plasma column is spatially uniform and directed along the z axis.The plasma column is neutral, i.e., the ion and electron densities are equal here.Loss processes in the plasma column are balanced by ionization, i.e., the net production rate S vanishes and Eq.͑4͒ gives i − dif − ␤ dis n e = 0. ͑A1͒ It follows that the charged particle density in the plasma column is spatially uniform in the framework of the model considered.The electric field in the plasma column, E ϱ , is related to the current density in the plasma column ͑which also is spatially uniform͒.The distribution of potential in the plasma column is a linear function of z: = E ϱ z + U, where U is a constant.This constant represents the difference between the value of the potential, = U, which is obtained by extrapolation to the cathode surface of the potential distribution in the plasma column, and the actual value of the potential of the cathode surface, = 0.It is natural to call this difference the sheath voltage.In fact, this is the only logical definition of a sheath voltage; any other definition would involve a "sheath edge," something that does not exist in reality ͑since the sheath joins the plasma column smoothly͒ and cannot be introduced except arbitrarily.The sheath voltage U is related to the current density in the plasma column, as is the electric field in the plasma column, E ϱ .
Thus, the missing boundary conditions for Eqs.͑1͒-͑3͒ in the near-cathode space-charge sheath assume the form ͑11͒.
Since both E ϱ and U are related to the current density in the plasma column, these quantities are not independent; therefore one of them may be considered as a given parameter and the other should be found as a part of the solution.It is convenient for the purposes of this work to consider U as a given parameter and E ϱ as a parameter to be found as a part of the solution.
The above needs to be slightly modified in the case where the model of a diffusion-controlled discharge is considered.In the framework of this model, the last term on the LHS of Eq. ͑A1͒ is discarded and this equation becomes independent of n e .The electric field E ϱ in the plasma column can take only one value, namely, the one that ensures i being equal to dif .Hence, in this case one must consider U as a given parameter and E ϱ is a root of the equation i = dif .Note that the charged particle density and the current density in the plasma column are not spatially uniform in this case and may depend on x and y ͑but not on z͒; they are governed by the solution in the near-cathode sheath.Let us first consider the case where the calculation domain includes the whole discharge from the cathode to the anode.The current continuity equation ͑12͒ in the 1D case may be integrated to give e i n i E + j e = j.͑B1͒ Here j e = e e n e E is the density of electric current transported from the plasma to the cathode by the electrons, E = Ј ͑the prime denotes d / dz͒.Multiplying Eq. ͑3͒ by E and transforming the RHS with the use of Eq. ͑B1͒, one can write

͑B5͒
where j 1 = j / 0 i and q 4 is given by Eq. ͑22͒.It was assumed for brevity that i and e are constant while writing this problem.
Let us consider the current density in the discharge, j, as a given ͑control͒ parameter instead of U. Then the nonlinear boundary-value problem for a second-order ordinary differential equation constituted by Eqs.͑B4͒ and ͑B5͒ is complete and may be solved numerically.Finding solutions for various j, one can construct a complete description of the CDVC of the 1D mode of operation of the discharge.In this work, this problem has been solved by iterations with the use of Newton's linearization and the Petukhov method ͓57͔.Now let us consider the case where the calculation domain includes only the near-cathode space-charge sheath.In this case, one can go a step further and reduce the problem to an initial-value problem for a first-order ordinary differential equation.Eliminating n e in terms of j e and n i with the use of Eq. ͑B1͒, one can rewrite Eq. ͑2͒ as where L = ␤ dis / e i e .Equation ͑B2͒ may be solved for EЈ: Values j ew and j eϱ of the electron current density j e at the cathode surface and in the plasma column may be expressed in terms of j by, respectively, using the boundary condition Eq. ͑5͒ and setting EЈ = 0 in Eq. ͑B7͒: ͑The index ϱ here and below denotes values of corresponding quantities in the plasma column.͒Note that the asymptotic behavior of solutions to Eqs. ͑B6͒ and ͑B7͒ at large z can be found by means of linearizing the equations around j eϱ and E ϱ and these solutions tend to j eϱ and E ϱ exponentially; see the next section.
It is convenient to transform the problem to the independent variable j e and the dependent variable E. A differential equation for the new unknown function E͑j e ͒ is obtained by dividing Eq. ͑B7͒ by Eq. ͑B6͒: There is an uncertainty on the RHS of Eq. ͑B9͒ at j e = j eϱ .It follows from results of the next section that the ratio EЈ / j e Ј tends to E ϱ C 2 / j eϱ at large z, where C 2 and are given by Eqs.͑B14͒ and ͑B16͒.Hence, the RHS of Eq. ͑B9͒ at j e = j eϱ equals E ϱ C 2 / j eϱ .Note that this result can be derived also by a direct resolution of the uncertainty.
The procedure of solution is as follows.A value of the current density in the discharge, j, is chosen.The value of the electric field in the plasma column, E ϱ , is found by solving Eq. ͑A1͒ written in the form After this, the first-order differential equation ͑B9͒ is numerically solved on the interval j ew Յ j e Յ j eϱ with the initial condition E͑j eϱ ͒ = E ϱ .After this equation has been solved, one can find the function j e ͑z͒ ͓and, consequently, the function E͑z͔͒ by means of numerical evaluation of an integral that follows from Eq. ͑B6͒: The sheath voltage corresponding to the value j being considered is given by the formula

Asymptotic behavior of 1D solution to the sheath equations at large z
In order to find the asymptotic behavior of solutions to Eqs. ͑B6͒ and ͑B7͒ at large z, let us represent j e ͑z͒ = j eϱ + w͑z͒ and E͑z͒ = E ϱ + v͑z͒, where w and v are small quantities, ͉w͉ Ӷ j eϱ and ͉v͉ Ӷ E ϱ .Retaining in Eqs.͑B6͒ and ͑B7͒ terms of the first order in w and v, one gets the equations where the constant coefficients C 1 and C 2 are While deducing Eqs.͑B13͒, the function K = K͑E͒ was expanded as K = K ϱ + K 1ϱ ͑E − E ϱ ͒.Note that, since the function K͑E͒ is growing, K 1ϱ Ͼ 0. For brevity, the treatment is restricted to the case where the derivatives ͉͑d i / dE͉͒ E=E ϱ and ͉͑d e / dE͉͒ E=E ϱ are zero.Equations ͑B13͒, being a system of linear equations with constant coefficients, have an exponential solution: w ϳ e z , v ϳ e z , where is a constant quantity satisfying the quadratic equation

͑B15͒
The roots of this equation are of opposite signs.The negative root is Thus, solutions to Eqs. ͑B6͒ and ͑B7͒ tend to j eϱ and E ϱ exponentially.

1D solution at small current densities
In the case where the calculation domain includes the whole discharge from the cathode to the anode, the limit of low current densities corresponds to the Townsend discharge, in which a distortion of the applied electric field by plasma charges is negligible.Let us introduce a normalized current density where E ͑0͒ = U / h is the applied electric field.͑Note that the kinetic and transport coefficients of the plasma to a first approximation may be treated as constant in the gap since the electric field in the gap also is constant to a first approximation.͒The parameter ␦ may be interpreted as the ratio of the density of space charge in the gap, j / i E ͑0͒ , to a characteristic value of the charge density needed to distort the applied electric field, 0 E ͑0͒ / h.The limiting case of the Townsend discharge corresponds to small values of this parameter.The electric field in the gap is sought as a superposition of the applied electric field and of a ͑small͒ distortion of the applied field due to plasma charges,

͑B18͒
where the ellipsis represents terms of the order of ␦ 2 E ͑0͒ and of higher orders.Substituting the expansion ͑B18͒ into Eqs.
͑B4͒ and ͑B5͒, expanding in ␦, and retaining linear terms, one arrives at

͑B20͒
Here ␣ ˜= ␣ − dif / e E ͑0͒ .Note that the recombination term of Eq. ͑B4͒ is of the order of q 4 ␦ relative to the term on the LHS; if the coefficient q 4 is of the order unity or smaller, the recombination term may be dropped.͑Using data cited at the beginning of Sec.IV, one finds q 4 = 5.8ϫ 10 −7 p / Torr, so this coefficient indeed is small under the considered conditions.͒A solution to Eq. ͑B19͒ satisfying the first boundary condition ͑B20͒ is In order that this solution satisfy the second boundary condition ͑B20͒, one should set This equation represents a well-known condition for initiating a self-sustained discharge and governs E ͑0͒ .
Let us now consider the case where the calculation domain includes only the near-cathode space-charge sheath.As the current density tends to zero, the densities of charged particles both in the plasma column and in the sheath tend to zero.As in the case of a parallel-plane discharge, a distortion of the electric field in the sheath by plasma charges becomes negligible and the recombination becomes negligible as well, i.e., the discharge becomes controlled by diffusion.The latter means that E ϱ , the electric field in the plasma column, satisfies to a first approximation the equation K͑E ϱ ͒ = 0 and is independent of j; the former means that the electric field throughout the sheath is close to E ϱ .
A normalized current density that plays the role of a small parameter is introduced by an equation similar to Eq. ͑B17͒ where K 1ϱ −1 plays the role of a length scale.We seek a solution to Eq. ͑B9͒ in the form of a two-term asymptotic expansion: E͑j e ͒ = E ϱ + ͱ ␦E ͑1͒ ͑j e ͒ + ¯.

͑B24͒
Substituting this expansion into Eq.͑B9͒ and expanding in ␦, one arrives to a first approximation at ͱ ␦ dE ͑1͒ dj e = 1 0 i j e ͑1 + ␤͒j e − j K 1ϱ ͱ ␦E ͑1͒ .

͑B25͒
Note that the second term in the denominator of the second multiplier on the RHS of Eq. ͑B9͒ is of the order of q 4 ͱ ␦ relative to the first term and therefore can be dropped.
A deviation of E ϱ from the root of the equation K ϱ ͑E ϱ ͒ = 0 may be estimated from Eq. ͑B10͒ and is of the order of q 4 ␤␦E ϱ .Hence, a boundary condition for Eq.͑B25͒ at j e = j eϱ is E ͑1͒ = 0.A solution reads Thus, the sheath voltage tends to a finite value at small current densities.One can see from Eq. ͑B27͒ that the sheath thickness grows infinitely ͑proportionally to ␦ −1/2 ͒ at small current den-sities.This obviously restricts the physical applicability of the results obtained.However, the results, in particular, the conclusion that the sheath voltage tends to a finite value at small current densities, may be useful for understanding bifurcations of solutions in the model considered.

APPENDIX C: ASYMPTOTIC BEHAVIOR OF SHEATH PERTURBATIONS AT LARGE z
Let us start with the particular case of a diffusioncontrolled glow discharge, where one can set ␤ dis = 0.One finds for the coefficients q 2 , q 3 , and q 5 at large z q 2 → n eϱ , q 3 Ϸ K 1ϱ ͑E − E ϱ ͒, q 5 → − n eϱ K 1ϱ E ϱ .͑C1͒ As shown in Appendix B, the difference E − E ϱ at large z decays proportionally to e z , where =− ͱ C 2 K 1ϱ in the case of the diffusion-controlled discharge being considered.Equation ͑34͒ to a first approximation assumes the form ͕e −z ͓fٞ − ͑k 2 + 2 ͒fЈ͔͖Ј = 0. ͑C2͒ This equation has exponential solutions f ϳ e z , with satisfying an algebraic equation of the fourth order, The latter equation has four roots, which are, in the order of their increase, 1 = − ͱ 2 + k 2 , 2 = , 3 = 0, 4 = ͱ 2 + k 2 .

͑C4͒
Hence, the asymptotic behavior of a general solution to Eq. ͑34͒ at large z includes in the particular case of a diffusioncontrolled glow discharge two exponentially decaying terms, a constant term, and an exponentially growing term.Now let us consider the general case where the discharge is not controlled by diffusion, ␤ dis 0. All the coefficients involved in Eq. ͑34͒ tend at large z to finite limits in this case, in particular, q 3 → − K ϱ , q 5 → − n eϱ ͑K 1ϱ E ϱ + K ϱ ͒, q 6 → 0 K ϱ eE ϱ ͑1 + ␤͒ .

͑C5͒
͑ i and e are treated as constant in this appendix, as in Sec.III; hence there is no need to attribute the index ϱ to ␤.͒ Equation ͑34͒ to a first approximation assumes the form This is an equation with constant coefficients which has exponential solutions f ϳ e z , with satisfying the algebraic equation of the fourth order Let us first find approximate solutions to Eq. ͑C7͒ in the limiting cases k → 0 and k → ϱ.In the limiting case k → 0, this equation has two finite roots, which are found by setting k = 0, and two roots of the order k, which are found by dropping the first and second terms of this equation.In the order of their increase, the roots are 1 = , 2 = − kͱ ␤C 1 Note that 2 tends to zero proportionally to k at small k, which means that the length of attenuation of the function f͑z͒ at large z, being equal to ͉ 2 ͉ −1 , is much larger than the length of attenuation of functions j e ͑z͒ − j eϱ and E͑z͒ − E ϱ , which equals ͉͉ −1 .This complicates a numerical solution of the eigenvalue problem.
In order to treat the limiting case k → ϱ, it is convenient to rewrite Eq. ͑C7͒ as ͑C10͒ Dropping the term on the RHS, one finds the roots ͑again in the order of their increase͒ One can see from Eqs. ͑C8͒ and ͑C11͒ that in both limiting cases Eq. ͑C7͒ has two positive roots and two negative roots, which means that the asymptotic behavior of a general solution to Eq. ͑34͒ at large z includes two exponentially growing terms and two exponentially decaying terms.Let us prove that the latter is the case not only for small or large k, but for any k.Equation ͑C7͒ may be written as a standard quartic equation, a 0 4 + a 1 3 + a 2 2 + a 3 + a 4 = 0 ͑C12͒ with the coefficients With the aim of invoking the Routh-Hurwitz criterion, let us consider the sequence where Given that C 1 , C 2 , and K 1ϱ are positive, the terms of the sequence have signs as shown below each term.There are two sign changes in the sequence.According to the Routh-Hurwitz criterion, Eq. ͑C7͒ has two roots with positive real parts and two roots with negative real parts.Hence, the asymptotic behavior of solutions to Eq. ͑34͒ at large z includes two exponentially growing terms and two exponentially decaying terms.
FIG. 1. Schematic of the discharge.
FIG. 2. ͑Color online͒ Relationship between the wave number of a spot mode and the value of the current density at which this mode branches off from the distributed mode.Discharge between parallel electrodes.p = ͑a͒ 75, ͑b͒ 100, ͑c͒ 150, and ͑d͒ 760 Torr.The filled and open circles represent different families of bifurcation points.
FIG. 4. ͑Color online͒ Points: normalized wave number of the spot mode vs current density at which this mode branches off from the distributed mode.Curves: CDVCs.Points and solid curve, discharge in Xe between parallel plates; dotted curve, near-cathode space-charge sheath in Xe. p = ͑a͒ 75 and ͑b͒ 760 Torr.
FIG. 5. ͑Color online͒ Points: normalized wave number of the spot mode vs current density at which this mode branches off from the distributed mode.Curve: CDVC.Near-cathode space-charge sheath.p = ͑a͒ 75 and ͑b͒ 760 Torr.The filled and open circles represent different families of bifurcation points.
ln x͒ 1/2 , ͑B27͒ U = E ϱ K 1ϱ −1 ln 1 + ␥ ␥͑1 + ␤͒ .͑B28͒ , ͑16͒ where = ͑E͒ = d ln͑1+␤͒ / dE.͑We recall that E is the modulus of the local electric field: E = ٌ͉͉.͒Eliminating n i on the RHS of Eq. ͑4͒ by means of Eq. ͑3͒ and substituting the resulting expression into Eq.͑16͒, one arrives at e Kn e − ␤ dis n e ͩn e − 0 e ٌ 2 ͪ + n e e ٌ E • ٌ ͑2͒ and ͑5͒, and from the first equation in Eq. ͑8͒ n i and n e with the use of Eq. ͑B3͒, one arrives at a boundary-value problem for the function E. Introducing a new unknown function = E 2 , one can write this problem in the form