Multiple solutions in the theory of direct current glow discharges: Effect of plasma chemistry and nonlocality, different plasma-producing gases, and 3D modelling

The work is aimed at advancing the multiple steady-state solutions that have been found recently in the theory of direct current (DC) glow discharges. It is shown that an account of detailed plasma chemistry and non-locality of electron transport and kinetic coefﬁcients results in an increase of the number of multiple solutions but does not change their pattern. Multiple solutions are shown to exist for discharges in argon and helium provided that discharge pressure is high enough. This result indicates that self-organization in DC glow microdischarges can be observed not only in xenon, which has been the case until recently, but also in other plasma-producing gases; a conclusion that has been conﬁrmed by recent experiments. Existence of secondary bifurcations can explain why patterns of spots grouped in concentric rings, observed in the experiment, possess in many cases higher number of spots in outer rings than in inner ones. V C 2013 AIP Publishing LLC . [http://dx.doi.org/10.1063/1.4826184]


I. INTRODUCTION
][6][7][8][9][10][11] The solutions [1][2][3] were found in the framework of the simplest self-consistent model of glow discharge, which accounts for a single ion species produced via a single effective ionization process and employs the local-field approximation.
The finding of these solutions represents the first step towards understanding and self-consistent modelling of selforganization in near-cathode regions of DC glow discharges.Many questions still remain unanswered.One of such questions is whether solutions similar to those reported in Refs.1-3 exist in the framework of more realistic models of glow discharges, which take into account multiple ionic species, various ionization channels, atomic excited states, excimers, and non-locality of electron transport and kinetic coefficients.A very interesting and important question is why the self-organization was observed in xenon but not in discharges in other plasma-producing gases, such as argon. 5,12nother important question is as follows.In patterns with many spots observed in the experiment, spots are grouped in concentric rings, and in many cases, the number of spots in outer rings is higher than that in inner ones.However, the number of spots in different rings in the modelling 3 is equal.It is very interesting to try to obtain in the modelling patterns with unequal number of spots in different rings similar to those observed in the experiment.
The aim of this work is to answer these three questions.The modelling is performed for a parallel-plane discharge tube in 1D, 2D (axial symmetry), or in 3D.A detailed model of glow discharge was employed for xenon.A basic model of glow discharge was employed for argon and helium.
The outline of the paper is as follows.The basic model employed for argon and helium and a detailed model employed for xenon are briefly discussed in Sec.II.The detailed model for xenon is compiled on the basis of data available in the literature and takes into account plasma chemistry and non-locality of the kinetic and transport coefficients of electrons (by means of differential equation of conservation of electron energy).The effect of plasma chemistry and non-locality on the multiple solutions is studied in Sec.III.In Sec.IV, multiple solutions for argon and helium are given.Patterns with unequal number of spots in different rings are given in Sec.V.The effect of absorbing walls of the discharge tube over 3D modes is studied in Sec.VI.Finally, conclusions are summarized in Sec.VII.

II. THE MODELS
The basic model employed for argon and helium is the well-known simplest self-consistent model of glow discharge, which accounts for a single ion species produced via a single effective ionization process and employs the driftdiffusion and local-field approximations.The model comprises equations of conservation of a single ion species and the electrons, transport equations for the ions and the electrons written in the drift-diffusion approximation (i.e., neglecting inertia of the charged particles; e.g., Ref. 13), and the Poisson equation.The (only) ionic species considered is Ar þ 2 in the case of argon discharge and He þ 2 in the case of helium discharge.The mobility of these ions in their parent gases was evaluated by means of the formula that has been obtained by fitting the measurements 14 and reads l i ¼ a i =n n , where n n is the density of the neutral gas, a i ¼ 7:1 Â 10 21 m À1 V À1 s À1 for argon and a i ¼ 5:5 Â 10 22 m À1 V À1 s À1 for helium.The mobility of the electrons was evaluated as l e ¼ a e =n n , where a e ¼ 10 24 m À1 V À1 s À1 for argon and a e ¼ 2:25 Â 10 24 m À1 V À1 s À1 for helium.Townsend's ionization coefficient was evaluated as a ¼ Cpexp½ÀDðp=EÞ 1=2 , where E is the electric field strength, p is the gas pressure, C ¼ 2:92 Â 10 3 m À1 Torr À1 and D ¼ 2:66 Â 10 2 V 1=2 ðmTorrÞ À1=2 for argon, and C ¼ 4:4 Â 10 2 m À1 Torr À1 and D ¼ 1:4 Â 10 2 V 1=2 ðmTorrÞ À1=2 for helium. 15(Note that this way of evaluation of l e and a yields results in good agreement with those given by the zero-dimensional Boltzmann solver BOLSIGþ. 16) The diffusion coefficients were evaluated by means of the Einstein relation with T i and T e the temperatures of the heavy particles and electrons equal to, respectively, 300 K and 1 eV.The coefficient of dissociative recombination of molecular ions was set equal to 5:38 Â 10 À14 m 3 s À1 for argon 17 and 3:75 Â 10 À17 m 3 s À1 for helium. 18he modelling reported in this paper refers to a cylindrical discharge tube with parallel plane electrodes.Boundary conditions at the cathode and anode are written in the conventional form: diffusion fluxes of the attracted particles are neglected as compared to drift (i.e., the normal derivative of the particle density vanishes); the flux of the electrons emitted by the cathode is related to the flux of incident ions in terms of the effective secondary emission coefficient c, which is assumed to characterize all mechanisms of electron emission (due to ion, photon, and excited species bombardment); 15 density of ions vanishes at the anode; electrostatic potentials of both electrodes are given.The effective secondary emission coefficient typically varies between 10 À3 and 10 À1 . 15In the present work, the effective secondary emission coefficient was set equal to 0.03 for both argon and helium.One boundary condition at the wall of the discharge tube is the conventional condition of zero electric current density, which with diffusion fluxes being neglected as compared to drift amounts to vanishing normal electric field.As far as the ions and the electrons are concerned, we will first assume that the wall is reflecting, or, in other words, will neglect neutralization of charged particles at the wall; then normal components of the charged particle flux densities vanish.The effect of neutralization of charged particles at the wall will be studied in Sec.VI.
The detailed model for xenon takes into account the following species: singly charged atomic ions, molecular ions, electrons, an effective species of excited atoms combining all of the excited states in the 6s manifold (6s½3=2 2 ; 6s½3=2 1 ; 6s 0 ½1=2 0 , and 6s 0 ½1=2 1 ), and excimers.In addition to direct ionization, stepwise ionization and ionization of excimers are taken into account.The detailed model comprises conservation equations for all species, transport equations for all species written in the drift-diffusion approximation for the charged species, and the diffusion approximation for excited atoms and excimers, the Poisson equation, and the equation of conservation of electron energy, which is written in the form recommended in Ref. 16.
The system of reactions taken into account in the framework of the detailed model, data necessary for calculations and corresponding references are summarized in Table I.The rate coefficients taken from Ref. 17 are functions of the electron temperature T e related to the average electron energy e in the usual way, T e ¼ 2 e=3. 17Reaction rates for processes of electron impact ionization and excitation from the ground state were calculated in terms of e using BOLSIGþ 16 and electron cross sections recommended in Ref. 19.Excited states of energy higher than 6s were assumed to decay instantly into the representative state Xe Ã .In other words, the total rate of excitation into higher excited states, which was computed with BOLSIGþ, was included in the source term of state Xe Ã .
The rate coefficient of ionization from state Xe Ã was estimated as 20 where M1I ¼ E M1I =T e with E M1I ¼ E I À E M1 being the energy difference between the ionized state Xe þ and the representative excited state Xe Ã .The mobility and diffusion coefficient of atomic ions Xe þ were evaluated as functions of the reduced electric field E=n n by means of the two-temperature displaced-distribution theory 22 (here, n n is the density of the neutral gas).Mobility of molecular ions Xe þ 2 was evaluated by means of the formula l i2 ¼ b=n n , where b takes the value 2:1 Â 10 21 m À1 V À1 s À1 , which is an approximation of the measurements. 14iffusion coefficient of molecular ions Xe þ 2 was found by means of the Einstein relation with T i ¼ 300 K. Mobility and diffusion coefficient of electrons as well as electron energy mobility and electron energy diffusion coefficient were evaluated in terms of e using BOLSIGþ 16 and cross sections recommended in Ref. 19.Diffusion coefficients of atoms in excited states and excimers were set equal to 10 À2 m 2 s À1 following Ref. 21. (Note that this value has little effect on results, which is in accordance with Ref. 21.)  Boundary conditions for electrons, ions, and the electrostatic potential are the same as in the basic model.Densities of excited atoms and excimers vanish at both electrodes.The effective secondary emission coefficient c was set equal to 0.03; no separate coefficients characterizing secondary emission due to bombardment by different ionic and excited species have been introduced.Deexcitation at the wall is neglected.
Trial computations involving equation of heavy particle energy have revealed that heating of the neutral gas is negligible under conditions of this work.Therefore, most of the results, including those reported in this paper, have been computed without account of heating of the neutral gas.
The problems describing both basic and detailed models have been solved by means of the commercial software COMSOL Multiphysics.The control parameter in calculations was discharge current I or discharge voltage U, depending on the slope of the current-voltage characteristics (CVC) at the state being computed.The key feature of the modelling is the use of a stationary solver.It is important to stress that stationary solvers allow one to decouple questions of numerical and physical stability, which is of critical importance for a systematic investigation of multiple solutions. 23The first step is to compute the 1D mode in a wide range of discharge currents.The second step is to compute bifurcation points by means of the eigenvalue solver (see Ref. 1 for more details).The third step is to construct the initial approximation for a given value of the discharge current as the 1D mode plus a small perturbation in the transversal directions.Multiple solutions should be sought for values of discharge current or voltage in the vicinity of a bifurcation point.The modelling was 2D (more precisely, axially symmetric) or 3D in the framework of the basic model, and 1D or 2D in the framework of the detailed model.Some further details of the numerical procedure are given in Sec.V.

III. EFFECT OF DETAILED KINETICS
Multiple 2D solutions computed in the framework of the detailed model described in Sec.II are given in this section.Calculations were performed for discharges in xenon under the pressure of 30 Torr; the radius of the discharge tube and the interelectrode gap are 0.5 mm.
We recall the pattern of axially symmetric multiple solutions found in the framework of the basic model without account of neutralization of charged particles at the wall of the discharge tube, cf.Fig. 1 of Ref. 2. The 1D solution, which is well known and goes back to von Engel and Steenbeck, exists at all discharge currents.Each 2D mode exists in a limited range of discharge currents and is constituted by two branches separated by two bifurcation points at which the 2D mode joins that 1D mode.One of the branches is associated with a spot pattern containing one central spot and possibly ring spot(s), while the other branch is associated with a pattern comprising ring spot(s) and no central spot.
The results obtained in the framework of the detailed model fit the same pattern.The current density-voltage characteristic (CDVC) of the 1D mode is shown in Fig. 1.This mode exists at all discharge currents and is generally similar to that calculated in the framework of the basic model.However, there is a very interesting difference: the 1D mode computed in the framework of the detailed model manifests a retrograde behavior in the current density range of 200 Am À2 Շ j Շ 300 Am À2 ; cf. the S-shape manifested by the line in Fig. 1.The existence of a retrograde section of the 1D mode is a surprising result, which is unrelated to multiple solutions and, as far as we know, has not been reported previously.The reason of its appearance is unclear.We note only that a similar retrograde behavior has been found in the framework of a similar detailed model developed for argon; the retrograde behavior in xenon disappears if stepwise ionization is neglected; 23 however, retrograde behavior in argon subsists without stepwise ionization.
The current-voltage characteristics of the 2D modes are very difficult to distinguish from the CDVC of the 1D mode.These modes are more conveniently represented in the coordinates (hji; j centre ), where hji is the average current density over a cross section of the discharge and j centre is the current density at centre of the cathode surface; see Fig.  an efficient way; however, the Plasma Module is meant to use time-dependent solvers and seems not to work with steady-state solvers.The question to what extent timedependent solvers can be useful in computing multiple solutions is extremely interesting; however, it is beyond the scope of this work and will be addressed elsewhere.Although 2D modes have been computed only in a part of their region of existence, enough information has been obtained to characterize these modes in some detail.Circles in Fig. 2 represent bifurcation points.Note that the bifurcation points were not found independently (in contrast to Ref. 2, where they have been found on the basis of bifurcation theory with the use of the eigenvalue solver of COMSOL Multiphysics); rather, they have been determined as states belonging to 2D modes where variations of discharge parameters in the radial direction become negligible.Following notation, 2 let us designate by a i and b i bifurcation points at which the i-th 2D mode branches off from or joins the 1D mode, a i being the bifurcation point positioned at higher discharge currents and b i being positioned at lower currents.Of the first 2D mode, two sections were found: b 1 A and b 1 B. Of the second 2D mode, two sections were found: b 2 C and b 2 D. A complete branch b 3 a 3 of the third 2D mode was computed in the entire region of its existence.
Schematics in Fig. 2 illustrate distributions of current density over the cathode surface corresponding to each section computed.One can see that the section b 1 B is associated with a spot at the centre of the cathode; the section b 1 A is associated with a ring spot at the periphery of the cathode; the section b 2 D is associated with a spot at the centre of the cathode and a (less bright) ring spot at the periphery; the section b 2 C is associated with a ring spot inside the cathode; the branch b 3 a 3 is associated with a ring spot inside the cathode and a ring spot at the periphery.The switching between the two patterns associated with the first mode takes place at the bifurcation point (b 1 ).The same is true for the second mode.For both modes, the branch that is associated with a pattern comprising a central spot is the one that bifurcates (from the points b 1 or, respectively, b 2 ) in the direction of lower discharge currents, the branch that is associated with a pattern without a central spot bifurcates in the direction of higher currents.The latter applies also to the direction of branching of the third mode from the point a 3 .
All the above features, including the type of pattern associated with each branch and the direction of bifurcations of all branches, coincide with those found in the framework of the basic model. 2 There is, however, a very interesting difference.Only two axially symmetric modes exist for the same conditions in the framework of the basic model; it was necessary to increase the discharge radius beyond 0.52 mm in order that the third 2D mode appears; cf.Fig. 4 of Ref. 2. On the other hand, at least three axially symmetric modes exist in the framework of the detailed model.One can conclude that the presence of different ion and excited species, various ionization channels, and non-locality of electron kinetic and transport coefficients favours appearance of multiple modes.
Examples of distributions of densities of atomic ions, molecular ions, electrons, atoms in excited states, excimers, and mean electron energy are shown in Figures 3-7.The maximum of mean electron energy is localized inside the cathode layer, as it should; results not shown here reveal the presence of field reversal at the edge of the cathode layer for all modes.
Patterns seen in distributions of densities of atomic ions, electrons, excited atoms, and excimers in the states belonging to the first 2D mode (Figs. 3 and 4) are similar, although the patterns seen in the distribution of excimers are somewhat more diffuse than those of the other species.The patterns seen in distributions of molecular ions are much less localized.The same is generally true for the states belonging to the second and third 2D modes (Figs.5-7); however, the degree of similarity of distributions of densities of atomic ions, electrons, excited atoms, and excimers progressively decreases.

IV. EFFECT OF DIFFERENT PLASMA-PRODUCING GASES
The discharges in argon and helium computed for the plasma pressure of 30 Torr and the radius of the discharge tube and the interelectrode gap of 0.5 mm in the framework of the basic model are obstructed (that is, the current density-voltage characteristic, CDVC, is rising for all current densities), 24 in contrast to that in xenon.The reason is the difference between the cross sections of elastic collisions of electrons with atoms.Unsurprisingly, the bifurcation analysis has shown the absence of bifurcations for argon and helium under these conditions, 24 which means that no multiple solutions exist and no self-organization is present.On the other hand, the bifurcation analysis 24 has revealed bifurcations for these gases if the pressure is sufficiently increased, and this suggests that multiple solutions and selforganization should exist for such conditions.
Accordingly, in this work, modelling for the same geometry as above was performed for the pressure of 75 Torr in argon and 530 Torr in helium.The CDVCs for these pressures are shown in Fig. 8 by the solid lines and reveal falling sections comparable to that in xenon for 30 Torr.Multiple solutions have been detected indeed and the pattern of selforganization is the same as in xenon.As an example, two 3D solutions are shown by the dashed lines c 8 f 8 and c 12 f 12 .(Here c i and f i designate bifurcation points at which the i-th 3D mode branches off from or joins the 1D mode, c i being the bifurcation point positioned at lower discharge currents and f i being positioned at higher currents.Modes c i f i with i 6 ¼ 8; 12 have been computed but are not shown for brevity.)Note that both modes have been computed in the whole range of their existence.
The patterns associated with the modes c 8 f 8 and c 12 f 12 are similar to those found in xenon for the pressure 30 Torr, 3 except that they are better pronounced, and are discussed here only in brief.In the vicinity of the bifurcation points, the mode c 8 f 8 is associated with a pattern comprising 6 spots at the periphery of the cathode.Away from bifurcation points, a pattern develops comprising 6 spots at the periphery and 6 spots inside the cathode and it is this pattern that is represented by the schematic in Fig. 8.In the vicinity of the bifurcation points, the mode c 12 f 12 is associated with a pattern comprising 4 spots at the periphery and 4 spots inside the cathode.Away from bifurcation points, a pattern develops comprising 4 spots at the periphery, 4 spots inside the cathode, and a spot at the centre, and it is this pattern that is represented by the schematic in Fig. 8.
Until recently, self-organized spot patterns have been observed in experiments with discharges in xenon but not in other plasma-producing gases, such as argon. 5,12The above numerical results indicate that such patterns in gases other than xenon should exist provided that the conditions, in the first place, the pressure, are right.For the particular case of krypton, modelling results (not shown in the present work) suggest that the minimum pressure required for the onset of self-organization is about twice as high as the required pressure for xenon under the same conditions.This prediction has been confirmed by recent experiments: 25 self-organization has been observed in krypton for plasma pressures higher than roughly twice the pressure required for the onset of selforganization under comparable conditions in xenon.

V. PATTERNS WITH UNEQUAL NUMBER OF SPOTS IN DIFFERENT RINGS
In addition to the 3D modes c i f i ði ¼ 1; 2; 3; …Þ that branch off from the 1D mode, the modelling has revealed the existence of further branching (secondary bifurcations): there are 3D modes that branch off from the modes c i f i .These modes may be conveniently represented in the coordinates (hji; j edge ), where j edge is the current density at a fixed point at the edge of the cathode surface.Examples are shown in Fig. 9. Here, the solid and dashed lines represent, respectively, the 1D solution and the 3D solution c 8 f 8 shown in Fig. 8; the dotted line g 1 d 1 and the dashed-dotted line g 2 d 2 represent 3D solutions that branch off from c 8 f 8 .
Evolution of distributions of current density over the cathode surface along the modes g 1 d 1 and g 2 d 2 is shown in Fig. 10.Note that the fixed point at the edge of the cathode surface to which the ordinate of Fig. 9 refers is the bottommost point of these schematics.Both modes possess azimuthal period of 2p=3, which is twice the period of the mode c 8 f 8 from which they branch off.Hence, the secondary bifurcations that occur at the points g i and d i represent pitchfork bifurcations usually called period-doubling.Such bifurcations have been detected previously in numerical modelling of plasma-cathode interaction in arc discharges 1 but not in the modelling of glow discharges.
There are 6 spots on the periphery of the cathode in the state d 1 .This state belongs also to the mode c 8 f 8 and has the azimuthal period of p=3.(Note that since the state d 1 is positioned in the vicinity of the bifurcation point f 8 , the associated spot pattern comprises 6 spots at the periphery rather than 6 spots at the periphery and 6 spots inside the cathode as shown by the schematic in Fig. 8.) On the section d 1 A, every second spot gradually moves from the periphery towards the centre of the cathode; this is how the period doubling occurs.A central spot is formed on the section ABC.The spots become less bright on the section CD.On the section Dg 1 , the central spot disappears and the spots inside the cathode gradually move towards the periphery.At the state g 1 , the spot pattern comprises 6 (hardly visible) spots at the periphery; the same pattern which is associated with the mode c 8 f 8 in the vicinity of the bifurcation point c 8 .
The period doubling on the mode g 2 d 2 happens differently.On the section d 2 EF, every second spot of the inner ring gradually moves towards the periphery and merges with a spot of the outer ring.The resulting spots are located close to the periphery: a pattern with 3 spots in the inner ring and 6 spots in the outer ring appears.
][6][7][8][9][10][11] It has been unclear up to now if the modelling is able to describe such modes, because all 3D modes computed up to now bifurcate from the 1D mode and the number of spots in each ring in these solutions is the same.(This is a characteristic feature of the spectrum of the eigenvalue problem which is associated with bifurcations of 3D modes from the 1D mode and has a solution in terms of the Bessel functions. 26) In this work, a mode with less spots in the inner ring than in the outer ring has been computed for the first time (the above-described mode g 2 d 2 ) and has been found to appear through secondary bifurcations.
Note that period-doubling bifurcations g 1 ; d 1 ; g 2 , and d 2 have been found in this work by means of the following procedure.The mode c 8 f 8 possesses azimuthal period of p=3.Since this mode possesses also plane symmetry, it is sufficient to compute only half the period of the solution.
Therefore, it is natural that the computation domain is a sector with azimuthal angle p=6, and such computations are indeed fast and efficient.On the other hand, such computations do not reveal the period-doubling bifurcations.If, however, the computation domain is a sector with azimuthal angle p=3, then in the vicinity of bifurcation points g 1 ; d 1 ; g 2 , and d 2 the iterations switch from the period of a solution describing the mode c 8 f 8 to half the period of a solution describing one of the modes g 1 d 1 or g 2 d 2 .For this to happen, the step on discharge current or voltage normally should be small.Period-doubling bifurcations will thus be detected.

VI. EFFECT OF ABSORBING WALLS
The above results have been computed for a discharge tube with a reflective wall.This section is concerned with the effect of neutralization of charged particles at the wall of the discharge tube.
The effect over 2D solutions for xenon discharges was reported in Refs. 1 and 2. It has been found that the effect over 2D solutions for argon and helium discharges is the same.Skipping figures for brevity, we note only that the number of solutions is reduced; a mode appears that exists at all discharge currents and comprises the Townsend, subnormal, normal, and abnormal discharges (the fundamental mode); 2D modes become disconnected from each other due to a destruction of bifurcations accompanied by exchange of branches.
The effect of neutralization at the wall of the discharge tube over 3D modes has not been investigated previously.In this work, this effect was described in the same way as in Refs. 1 and 2 for the case of 2D modes: the boundary condition for electrons and ions at the wall is written in the form ðs À 1Þn e;i =R þ s @n e;i =@r ¼ 0, where s is a given parameter that varies between 1 and 0 and R is the radius of the discharge tube.s ¼ 1 corresponds to a (totally) reflecting wall, or, in other words, to losses of the charged particles due to neutralization at the wall being neglected.s ¼ 0 corresponds to an absorbing wall, or, in other words, to neutralization at FIG. 7. Distributions of parameters in the interelectrode gap.Xenon, p ¼ 30 Torr; R ¼ h ¼ 0:5 mm.State E of the branch of the third 2D mode comprising a ring spot inside the cathode and a ring spot at the periphery.the wall being taken into account in full.The computations started from s ¼ 1 and then s was gradually reduced.Note that computations turn increasingly slow with reduction of s due to the thin diffusion layer that appears near the wall of the discharge tube and needs to be resolved.For this reason, parameter s in this work was not lowered down to s ¼ 0: calculations were stopped as densities of both electrons and ions at the wall were at least one order of magnitude lower than maxima of densities inside the discharge filaments attached to each of the cathodic spots.
As an example, the effect of neutralization at the wall of the discharge tube over the 3D mode associated with a normal spot is shown in Fig. 11.(Detailed numerical results on this mode without neutralization can be found in Ref. 3.) Neutralization at the wall forces the spot to move away from the wall while the maximum current density in the spot does not change, an effect which could be expected.
An example of a more complicated effect is seen in Fig. 12.Without neutralization (Fig. 12(a)), there are three spots at the periphery, three spots inside the cathode, and a central spot.Neutralization at the wall forces the spots at the periphery to move away from the wall, while the spots inside the cathode move towards the wall; the size of the central spot is increased.For s ¼ 0:66, there is only one ring of spots which comprises six spots near the periphery of the cathode; the azimuthal period of the pattern is now p=3 rather than 2p=3.This is an indication that a pitchfork bifurcation has just been crossed in the direction of periodhalving.It is interesting that the appearing pattern (six spots near the periphery and a central spot) is similar to the one occurring on the mode c 8 f 8 in the vicinity of the state f 8 except that the latter does not reveal the central spot; compare Fig. 12(b) and the distribution corresponding to state d 1 in Fig. 10(a).

VII. CONCLUSIONS
The inclusion in the model of detailed kinetics (multiple ionic species, atomic excited states, and excimers) and nonlocality of electron transport and kinetic coefficients favours existence of multiple solutions describing different modes of current transfer to cathodes of DC glow discharges: three 2D modes were found, in addition to the fundamental (1D) mode, for xenon in the framework of the detailed model, where only two 2D modes exist in the framework of the basic model under the same conditions.The multiple solutions fit the same pattern that those found in the framework of the basic model.One can conclude that the existence of multiple solutions is an inherent feature of self-consistent models of DC glow discharges.
Numerical results for argon and helium give a clear indication that the conditions of argon and helium microdischarges are less favorable for self-organization than those in xenon.This is a consequence of the difference between cross sections of elastic collisions between electrons and atoms of these gases.However, multiple solutions for discharges in argon and helium do exist at higher pressures.These results suggest the possibility of observation of self-organization in DC glow microdischarges in experiments with gases other than xenon provided that the experimental conditions, such as pressure, are right.This prediction has been confirmed by recent experiments. 25he existence of secondary bifurcations on 3D modes is demonstrated; 3D modes branching off through these bifurcations have their azimuthal period doubled compared to that of the parent 3D mode.In this way, patterns with unequal number of spots in different rings, similar to those observed in the experiment, [4][5][6][7][8][9][10][11] have been obtained in the modelling for the first time.
Neutralization at the wall of the discharge tube forces a normal spot to move away from the wall while the maximum current in the spot does not change.Another effect caused by neutralization at the wall, which is very interesting theoretically, is period-halving.
Important questions still remain open.In particular, it would be very interesting to try to observe and investigate in the experiment self-organization in DC glow microdischarges in different plasma-producing gases moving in the direction of higher pressures as predicted by the modelling.The first step in this direction has already been taken: 25 selforganized patterns of spots have indeed been observed in krypton for discharge pressures higher than 120 Torr, which is roughly twice as high as the pressure required for the onset of self-organization under comparable conditions in xenon.Also very interesting would be to try to numerically calculate the multiple solutions by means of a non-stationary solver. 2 FIG.1.CDVC of the 1D mode calculated in the framework of the detailed model.Xenon, R ¼ h ¼ 0:5 mm; p ¼ 30 Torr.

FIG. 3 .
FIG. 3. Distributions of parameters in the interelectrode gap.Xenon, p ¼ 30 Torr; R ¼ h ¼ 0:5 mm.State B of the branch of the first 2D mode comprising a spot at the centre of the cathode.

FIG. 4 .
FIG. 4. Distributions of parameters in the interelectrode gap.Xenon, p ¼ 30 Torr; R ¼ h ¼ 0:5 mm.State A of the branch of the first 2D mode comprising a ring spot at the periphery of the cathode.

FIG. 5 .
FIG. 5. Distributions of parameters in the interelectrode gap.Xenon, p ¼ 30 Torr; R ¼ h ¼ 0:5 mm.State D of the branch of the second 2D mode comprising a spot at the centre and a ring spot at the periphery of the cathode.

FIG. 6 .
FIG. 6. Distributions of parameters in the interelectrode gap.Xenon, p ¼ 30 Torr; R ¼ h ¼ 0:5 mm.State C of the branch of the second 2D mode comprising a ring spot inside the cathode.

ACKNOWLEDGMENTS
This work is dedicated to the memory of Lev Tsendin, who was a good friend of one of us since the early 1980s.Modelling presented in Sec.III was triggered by discussions with Professor Tsendin of the effect of non-locality of electron transport and kinetic coefficients on the multiple solutions.This work was supported by FCT of Portugal through projects PTDC/FIS-PLA/2708/2012 Modelling, understanding, and controlling self-organization phenomena in plasmaelectrode interaction in gas discharges: from first principles to applications and PEst-OE/MAT/UI0219/2011.

TABLE I .
Kinetic scheme for xenon and data used to describe each process.