A kinetic model for the charged triple layer in low pressure arc discharges

A Kinetic Model for the Charged Triple Layer in Low Pressure Arc Discharges

A. Tomimura

Instituto de Física, Universidade Federal Fluminense,

Gragoatá, CEP 24210.340, Niterói, RJ, Brazil

e-mail: [e-mail] and [e-mail]

and

H.S. Maciel

Instituto Tecnológico da Aeronáutica, CTA,

Departamento de Física,

CEP12228-900, S.José dos Campos, SP, Brazil

e-mail: [e-mail]

Received 27 February, 1998

A one dimensional model of a peculiar configuration of charged layers in equilibrium composed by one electron rich layer surrounded by two ion rich layers adjacent to plasmas at distinct potentials and which is formed in a low pressure arc discharge (usually know as a triple layer) has been constructed using the BGK method ^[1] viz., with the help of the Poisson-Vlasov's system of equations applied to the free and reflected populations of electrons and ions in a supposedly existing electrostatic potential with the free populations assumed to be monoenergetic beams and the reflected ones obeying the Maxwell-Boltzmann distribution. Sagdeev potentials derived for the charged region and matched by appropriate plasma boundary conditions are numerically integrated to obtain the electrostatic potential for some set of free input parameters, compatible with those of a specific group of experiments. Limitations of the model are addressed to.

I.Introduction

Space configurations of multiple charged layers, from single sheaths to quadruple ones, have long been reported in gas discharge experiments, either as the intermediate media separating a plasma from electrodes and probes ^[2-5] or as a region well inside the plasma where charge neutrality is violated and potential jumps are created ^[6,7,8] . They have also been held responsible for the anomalous resistivity in current carrying plasmas when the current density exceeds a certain critical value ^[9] and for particle acceleration in the ionization and excitation processes of neutrals in auroras ^[10] . Distinct physical reasons have been given for its onset but in many cases this is still an open question. Our concern in this work however is with the conditions which sustain the layer in its steady state. Hence, a one-dimensional model of a triple layer in steady state was constructed to numerically simulate the experimental potential profiles observed in some experiments on low pressure arc discharges.

The numerical solutions yield by this model are checked against a specific group of experimental results ^[6] , with most of the free input parameters of the model compatible with the corresponding ones in the experiment. However,some of these are left undetermined (there are more free parameters than boundary conditions to fix them all) so that the solutions are found in limited regions of the input parameter space (windows) rather than a unique solution for a specific set of input parameters. Also, direct comparisons with experiment were not possible in view of the fact that two of these free parameters were reported in the experiments as estimates (not measured values) of some crude model of their own.

II. Basic assumptions

The triple layer commonly observed in practice is described in this work as the result of five distinct populations of charged particles coexisting in equilibrium in some region along the discharge column, viz., an ion beam (free ions) coming from the anode region (the right side of the potential trough), thermal ions which are reflected by the potential barrier at the anode region, an electron beam (free electrons) directed from left to right and thermal electrons on the left and on the right of the potential trough (see Fig.1.a). The resulting eletrostatic potential profile is determined by the Poisson-Vlasov's system of equations with the free populations assumed to be monoenergetic and the thermal ones obeying the Maxwell-Boltzmann distribution. Phase space for both electron populations is illustrated in Fig.1.b. The model presents an ion "hole" in phase space,i.e., there is no thermal ions in the potential trough (see Fig.1.c). The equations are written down in two separate regions along the distance x of the column with the origin (separation plane) at the bottom of the potential trough f_m , the cathode side on the left with the eletrostatic potential f® 0 and the anode side on the right with f® f₀ , the last quantity standing for the flat top of the potential barrier; these values are supposed to be known from the experiments. The other relevant quantities in this model are the electron density of the beam directed left to right and near the cathode (at x®-¥ ) n_e0 , the ion density of the beam from right to left and near the anode (at x®+¥ ) n_i0 , the electron density and temperature of the thermal electrons near the cathode, n_ec and T_ec , and near the anode, n_ea and T_ea , the ion density and temperature of the thermal ions n_i and T_i , and finally, the electron and ion beams initial kinetic energies ef_e0º m_ev_e0²/2 and ef_i0º m_iv_i0²/2 , respectively.

Figure 1. (a) A sketch of a triple layer potential profile with its five populations indicated by the arrows; (b) phase space for the electron populations; (c) phase space for the ion populations.

Finally, with typical values for temperature and density around 2.eV and 10¹⁵m^-3, mean free path is about 10⁴cm ,i.e., much larger than the typical layer width observed in the experiments ( ~ cm); also, special laboratory conditions can provide ionization rate below 1% ^[8] .For these reasons, the layer is assumed to be free of collisions and ionization processes.

III. The model equations

Having established the basic assumptions,we now write down the various elements leading to Poisson's equation and the Sagdeev potential in analytical forms in two adjacent regions inside the charged region separately viz., to the left and to the right of the potential minimum.

Firstly, steady-state assumption and continuity equation lead to n_ebv_eb = n_e0v_e0 and n_ibv_ib = n_i0v_i0 for the electron and ion beam respectively and energy conservation applied on each of them leads to and where the index (0) for the electrons indicates values on the far cathode side (left on the Fig.1a. diagram) and for the ions on the far anode side (right on the same diagram).

and

where the index (0) for the electrons indicates values on the far cathode side (left on the Fig. 1a. diagram) and for the ions on the far anode side (right on the same diagram).

According to assumptions, the electron and ion number densities of reflected populations are obtained from their velocity distribution functions,

where c_e and c_i are normalization factors to be determined according to the number densities of reflected electrons and ions, n_ec and n_i , both defined on the far cathode side. The distribution function for the electrons reflected to the right have the same form as in eq.3 but the normalization factor is determined by its number density on the far anode side, n_ea.

The integrations in velocity space in order to get the number densities are restricted to well defined intervals which depend on the potential profile: for both electron populations ( the one reflected to the left and the other to the right of the potential trough), the limits are ± (2e/m_e)^1/2(f-f_m)^1/2 but for the ion population, the limits are ± [2e(f₀-f)/m_i]^1/2 for f > 0 and -[2e(f₀-f)/m_i]^1/2 , -(-2ef/m_i)^1/2 together with (-2ef/m_i)^1/2 , [2e(f₀-f)/m_i]^1/2 for f < 0 , the gap between the two subintervals accounting for the void in phase space of the reflected ions (there is no trapped ions in the potential trough). Now, having in mind the even parity of the distribution functions in velocity space and the respective adjustment of the normalization factors, the integrations yields

where f_m£ f £ 0 (cathode side, left on Fig.1b diagram),

where f_m£ f £ f₀ (anode side, right on Fig.1b diagram),

valid for the cathode side (left on Fig.1c diagram) and

valid for the anode side (right on the same diagram).In the formulae above,the square roots must be read as applied on the argument of the error function and H is the usual Heaviside function, used here to describe two distinct expressions for N_ia , viz., the one for f > 0 (H = 0) and the other for f < 0 (H = 1) ; this is a consequence of the " ïon hole" in phase space (see Fig.1c).

Poisson's equation can then be written as

and

valid for x < 0 and x > 0 , respectively.

At this point, the following non-dimensional quantities are defined in order to simplify the notation used hereinafter:the normalized potentials y = ef/kT_ec , y₀ = ef₀/kT_ec and y_m = ef_m/kT_ec , the electron and ion Mach numbers M_e = (2ef_e0/kT_ec)^1/2 , M_i = (2ef_i0/kT_ec)^1/2, the temperature ratios r_i = T_ec/T_i and r_a = T_ec/T_ea , the number density ratios n_e0 = n_e0/n_ec , n_i0 = n_i0/n_ec , n_ea = n_ea/n_ec and n_i = n_i/n_ec and the normalized distance x = x/l_D , where l_D = (e₀ kT_ec/n_ec e²)^1/2 is the Debye lenght and x is the distance measured from the point where f = f_m . Pre-sheat boundary conditions^[11] have been circumvented by solving the problem in a Debye's length scale where the infinities are located at both ends of the layer.

With notation above^[12] , Poisson's equations are integrated from y_m to y making use of the identity 2d²y/dx² = (d/dy)(dy/dx)² and the boundary condition [dy/dx]y_m = 0 . The double integrations implicit in the integration of an error function (erf) are reduced to a single one (another erf) by inversion in the order of integration and the definition of error function (see any recommended textbook of differential and integral calculus on the topic reduction of multiple integrals). The results for both intervals read,

valid for y_m£ y £ 0 ,x £ 0 and

valid for y_m£ y £ y₀ ,x ³ 0 .

Now, in order to integrate equations (12) and (14) numerically, boundary conditions (plasma conditions) have to be imposed over the pair of eqns.(11), (12) and (13),(14), viz., d²y/dx² = 0 and = 0 , the former pair at y(-¥) = 0 and the latter at y(+¥) = y₀ . The resulting set of equations relating the free parameters reads

An additional relation between the electron density of the thermal electrons near the cathode, n_ec ,and normalized parameters associated with the beams, can be obtained if the current density J is supposed to be known in the experiment. In this case, we write J/e = n_e0v_e0+n_i0v_i0 , divide by n_ec (kT_ec)^1/2 , and make use of the normalized parameters defined before; from the resulting relation one gets

where j and t_ec are the numerical values of the (uniform) current density J in Amps/m² and the temperature T_ec in eV respectively.

At this point, the model is now ready for at least a partial feasibility proof. Partial because experimental data available are not so proliferous and suficiently precise so as to give it (or not) support within its limited validity. Qualitative comparisons however are always possible. This checking has been made by feeding a computing code prepared for solving the equations above with a specific set of experimental data parameters available at present. In doing so,we have gone far beyond the scope of the initial objective by searching for larger parametric spaces where other solutions can be found.

IV.Numerical results and discussion

Now, the boundary conditions together with the additional constraint (19) above, make up five equations relating the ten free parameters n_i, n_e0, n_i0, n_ea , M_e (or, alternatively, f_e0 ), M_i (or, alternatively,f_i0 ), r_i (or,alternatively,T_i ), r_a, n_ec and T_ec .

Experimental conditions ^[8] suggest r_i values between 10. and 50. , and the known value of the current density fixes the value of n_ec through relation (19),so that one ends up with four equations having to be solved for four of these in terms of the other four. We chose f_e0, f_i0, r_a and T_ec as the independent free parameters and ran the computing code, inputing values of the last two for each pair of the first two. The input values of r_a were such as to exclude negative values of densities in the process of solving the set of boundary conditions. The resulting set of conditions, in its turn, feeds a numerical integration scheme for equations (12)and (14) by putting S^1/2(y) = dy/dx and integrating in y , i.e.,

where the upper limit in the integration is always kept under y₀ for x > 0 and under y = 0 for x < 0 . Invariably in the middle of the numerical process, further changes on the values of r_a were needed so as to exclude negative values of S.

Fig. 2 illustrates one tipical case, with y,r and S standing for the potential, charge density and Sagdeev's potential normalized to kT_ec , en_ec and (kT_ec/l_D)² respectively, with the values (see Fig.2) of the input parameters compatible with the experimental conditions leading to the formation of charged layers in low pressure mercury-arc discharges(compare with Fig. 4 and 5 in ref.[6]). For all cases run we kept m_e/m_i = 2.716·10^-6 , J = 1.414·10⁴ Amps/m² , f₀ = 24.Volts and f_m = -7.5 Volts , the first two used and the last two measured on the axis of the discharge tube ^[6] .

Figure 3. A typical solution with initial input values f_i0 = 16.Volts,f_e0 = 30. Volts, T_ec = 2.eV , iterated r_a = .11396 , and calculated n_ec = 4.88·10¹⁵m^-3 ,l_D = .150 mm .

Figure 4.
Windows a and b for solutions of the triple layer in the parametric space T_ecvs.T_ea with fixed values of (f_e0,f_i0) = (30.,16.)eV and (31.4,46.)eV respectively, both with a common value of r_i = 36. . Upper bound values for (T_ec,T_ea) are, for (a), (4.00,19.19)eV, and for (b), (8.32, 36.18)eV; lower bound values are for (a), (.10, 15.67)eV, and for (b), (.10, 21.76)eV.

Figure 5.
Average widths of the triple layer for the windows a and b in Fig.4, w₁ standing for the length of the region where r < 0 and w₂ the length joining the two points of maximum r , each one averaged over their values in the window for each T_ec .

Searching for solution in the parametric space f_e0 vs.f_i0 by a double iteration procedure in f_e0 or f_i0 and r_a and with a fixed value of r_i (e.g.,r_i = 36. ) one finds a region like the one shown in Fig.3 and described by the curve which limits the minimum values of f_e0 and f_i0 for solution. The smallest minimum for f_e0 , the electron beam energy, is limited by the value of -f_m (the potential energy barrier for electrons associated with the dip of the potential profile) but the smallest value of f_i0 depends on the (fixed) values of r_i and T_ec , the last remark also applied to the asymptotes to the curve. For relevant values of r_i and T_ec other than the one used in the example of Fig.3 however, the new curves constructed on the same basis do not show appreciable departures from it so that one can grossly speaking, take that curve as the representative one which delimits the region in the parametric space f_e0 vs.f_i0 for solution. Next, for each point (f_e0,f_i0) in such a region, the solutions can only be found in its complementary parametric space T_ecvs.T_ea over a limited and bounded region (window). The windows were constructed with specific sets of free input parameters borrowed from the ones used in an incomplete model fitting the experimental results in ref.[6] so as to allow comparisons. Just two tipical cases are shown in the plots of Fig.4 , viz., window (a) corresponding to the values of (f_e0,f_i0) = (30.,16.)eV and window (b) to (31.4,46.)eV. Each of these windows has its boundary represented by two curves having common points at its closure at its lower and upper bound values. Incidentally, the inequality f_i0 > T_ea/2 is always satisfied in this model, which is reminiscent of a Bohm's criteria in its inequality form for the onset of a plasma-sheat edge near a negative wall ^[4] . This imposes an upper limit for T_ea which, as a consequence, drives an upper bound for T_ec . Notice that by fixing the value of r_i for each window, one is imposing one further constraint on the model. Comparison of windows (a) and (b) of Fig.4 hint at the way the parameters T_ec and T_ea evolve with increasing values of f_e0 and f_i0 , viz., the increase in the range of their values. The widths of the layer for these cases are illustrated in Fig.5, where w₁ stands for the length of the region where r < 0 and w₂ the length joining the two points of maximum r , each one averaged over their values in the window for each Tec . Comparison of widths in both windows shows that a significant increase in widths is obtained the larger the range of definition of T_ec , i.e., the larger the extension of the window along the T_ec values. One can also see that the upper bound for T_ec is more sensitive to the values of f_i0 , viz., greater will be T_ec as f_i0 increases. Unfortunately, the complexity of the system of equations turns it impossible to define a scale length, representative of the layer, in terms of the relevant free parameters, in analitical form. Besides, the entire (graphical) picture would demand a great number of extra computations which might only be appropriate in a more compreensive version of this work.

V. Limitation of the model

A close look at the sequence of experimental plots in Fig.10, ref.[7], shows that the potential profiles along the direction of the discharge depend strongly on the radial distance from the axis where they were measured. One can see for example that the nearer to the axis the potential were measured the thinner and higher became the profile in the discharge direction. The potential profiles exhibited show clearly that charge distribution is at least two-dimensional in character and a corresponding 2D model should be more appropriate to describe the layer. Some radial profiles would have to be prescribed for the distribution functions of the thermal populations through their parallel and perpendicular temperatures, along and across the discharge direction as well as for the number densities on the far cathode and anode sites. Besides,the maximum and minimum potential should be left to be determined self consistently by the resulting system of equations describing the model. In our 1D model the radial dependency of the charge distribution was circumvented to extend our 1D code to any line parallel to the axis by the prescription of f_m and f₀ with a corresponding set of measured values on that line,at the expense of self consistency. This is clearly one limitation of the model. To turn it entirely self consistent however, one would need two more restrictions on the free parameters (for exemple,one restriction could be to get solutions with the potential profile satisfying a certain value at a typical inversion point of these profiles, like in ref.[6]), which can only be done at the expense of a reduction on the class of solutions we are looking for, since all the boundary conditions on the free parameters have been employed.

A quick comparison of the width of the potential profile exhibitted in Fig.2 ( ~ 2 mm) with the width of the experimental profile presented at ref.[6] ( ~ 1 cm) might show at a first sight another limitation of the model. However, the experimental profile has not been presented together with measured parameters associated with all of our free parameters, e.g., f_e0 and f_i0 ; for instance, thicker layer can be found at higher values of f_i0 as can be seen from Fig.5. Comparison is made even more difficult for the lack of further experimental data and corresponding potential profiles. In particular, the ones mentioned in Fig.10, ref.[7], cast some doubts about the precision of the experiments for measurements on the axis (the potential dip in the last plot seems more like a spike). This leaves somewhat undetermined such a comparison and gives no grounds to discard the present model.

As one can infer, the least necessary number of boundary conditions and constraints were used in the present model. Additional boundary condition like = 0 at plasma-sheath boundary (-¥,+¥ in the present model) ^[13] or some additional constraint restricting the potential profiles to some form ^[6] would reduce the windows regions to lines or even to points, which are done in these works in order to solve their respective model equations in a closed and self consistent way. Whatever the physical and mathematical reasons to justify the use of the additional boundary condition or constraint (in particular, the mentioned one in ref.[13] does not apply in the present model since the quantity above is singular at those points, unless a non-maxwellian velocity distribution is used), it is not easy to reproduce them in the laboratory, and one might conclude that perhaps windows rather than lines or points would be more like a rule than an exception in practice.

VI. Conclusion

One can say in conclusion that although the results above may seem pertinent to a particular laboratory experiment since most of the input parameters used in the exhibitted graphs in this work have been borrowed from the results of an incomplete model whose aim was to try to fit them with their experimental results ^[6] , they are nevertheless representative enough to validate the model for other aplications. All the exhibitted cases will show similar features for any point in the region delimited by the curve shown in Fig.3. By studying the windows boundaries one can draw some conclusions about the relationships of the relevant parameters used in the experiments which would sustain the layer in its steady-state. Even though these relationships are, most of them, of qualitative character, since even the more obvious simplification of the complete equations are analitically intractable and there is no Bohm's-like criterion ^[4] for such a model, the findings in this work may be useful in the sense that they point out the way of establishing the relevant input parameters to be chosen in similar experiments.

Further results and discussion will be published elsewhere.

Acknowledments

The authors are indebted to Dr.J.E.Allen from the Eng.Dept., Oxford University, UK, for valuable discussions.

References

1 I. B. Bernstein, J. M. Greene, M. D. Kruskal, Phys. Rev. 108, 546 (1957).

2 W. Schottky, Phys. Zeits. 25, 342 (1924).

3 L. Tonks, I. Langmuir, Phys. Rev. 34, 876 (1929).

4 D. Bohm D., The characteristics of Electrical Discharges in Magnetic Fields, ed A. Guthry and R. K.Wakerling, NY, MacGraw-Hill, 1949.

5 F. W. Crawford, A. B. Canara, J. Appl. Phys. 36, 3135 (1965).

6 H. S. Maciel, J. E. Allen, Proc. 2nd Symp. on Plasma Double layers and Related Topics, Insbruck, Austria, pg 218 (1984).

7 H. S. Maciel, J. E. Allen, J. Plasma Physics 42, 321 (1989).

8 H. S. Maciel, PhD thesis, Oxford University, Dept. of engineering Science, Oxford, England, 1985.

9 DeGroot et al., Phys. Rev. Lett. 38, 22, 1283 (1977).

10 M. Sanduloviciu, E. Lozneanu, Int. Conf. on Phenomena in Ionized Gases, Hoboken, NJ, USA, 1995, ed. Stevens Inst.of Tech., pgs.21-22,172.

11 K-U. Riemann, J. Phys. D: Appl. Phys. 24, 493 (1991).

12 H. Yamada, Z. Yoshida, J. Plasma Phys.48, 229 (1992).

13 J. G. Andrews, J. E. Allen, Proc. Roy. Soc. Lond.,A, 320, 459 (1971).

Publication Dates

• Publication in this collection
  13 May 1999
• Date of issue
  Sept 1998

History

• Received
  27 Feb 1998