One-dimensional multi-fluid model of the plasma wall transition

A one-dimensional multi-fluid model of the plasma wall transition is presented. The plasma is composed of singly charged positive ions, electrons and negative ions. Positive ions are assumed to obey the continuity, momentum transfer and energy transport equation, while the electrons and negative ions are treated by the continuity and momentum transfer equation. For the positive ions the closure is done by assuming that heat flux in the plasma is zero, while both negative species are treated as isothermal, with the pressure gradient term being proportional to the density gradient. The potential drop in the pre-sheath can be determined either by the electrons or the negative ions – depending on their densities in the unperturbed plasma. It is shown how the boundary density of both particle species can be in some cases determined very precisely.


Introduction
The plasma sheath problem is one of the oldest in plasma physics [1], but it still attracts the attention of many researchers.Plasmas containing negative ions [2] are important in many areas of plasma physics [3].In such plasmas the pre-sheath potential drop can be determined by either the electrons or the negative ions [2].Which of the two predominates, depends on the density of the two species.With the presented model the limiting density can be determined very accurately and very easily.

Model
Plasma is assumed to be composed of three groups of charged particles: singly charged positive ions, electrons, and singly charged negative particles -called particles 2. The electrons and particles 2 are treated by the continuity and momentum transfer equations, while the positive ions are treated by the continuity, momentum transfer, and pressure equations.For simplicity, all collisions in which momentum or energy is transferred between different types of particles are neglected.The system of model equations reads as follows: The potential profile (x) is determined by the Poisson equation: The meaning of the symbols is standard: ni, ne, n2 are particle densities, ui, ue, u2 are fluid velocities of the respective particle species, mi, me, m2 are particle masses, Ti, Te, T2 are particle temperatures, i, e, 2 are density part of the source terms with the unit m -3 , i, e, 2 are temporal part of the source terms, kB is the Boltzmann constant and x is the space coordinate.
The system of equations ( 1) -( 6) can only be solved numerically and for that purpose it must be transformed in a dimensionless form.The following variables are introduced: Here n0 is an arbitrary density used for normalization and X is the normalized space coordinate -see below formulas (15) -(18).Using (7) the system (1) -( 6) is written in the following form: Alternatively, the Poisson equation ( 13) can be replaced by the neutrality condition:  Both systems of equations ( 8) -( 13) and ( 8) -( 12), ( 14) are systems of eight ordinary differential equations with eight unknown functions of X: Ni(X), Ne(X), N2(X), Vi(X), Ve(X), V2(X), (X) and (X).They are both solved as initial value problem.Initial values (at X = 0) of the unknown functions are selected and then the system (8) -( 13) or ( 8) -( 12), ( 14) is integrated numerically in the positive X direction.For a unique solution of the system (8) -( 12), ( 14) eight initial conditions must be selected: For the system (8) -( 13) on the other hand nine initial conditions must be specified, since the Poisson equation ( 13) is of the second order: .
The additional initial condition is the initial value of the electric field 0.Quantities that are selected independently, as given data are: the neutrality parameter , the source terms si, se and s2, the masses  and  as well as  -the temperature of particles 2. In addition the normalization of X must be selected, this means that either (16) or (18) must be selected for A1 and A2.
Before some results of the model are presented a few more remarks about the model should be added.
First let us define the so-called polytrophic function (X) [4].The ion pressure gradient, which appears in ( 9) and ( 12) can be written in the following form: ( ) From (21) the polytrophic function [4] can be read: So, Eq. ( 9) can be written in an alternative form: Next the ion density gradient term is eliminated from the first of the continuity equation (8): Then ( 24) is inserted into equation ( 23): After some rearrangement (25) can be written in the following form: The first term on the right hand side of ( 26) is positive, since the derivative of the potential is negative.
As long as Vi is large enough the second term cannot prevail and the right hand side of (26) remains positive.But the left hand side of ( 26) is positive only as long as If Vi drops below this value, a singularity appears.This puts the limit to the selection of Vi0.Even if 0 = 0 is selected, Vi0 must be larger than zero.
Next electron density gradient is eliminated from equation ( 8): When ( 28) is inserted into equation (10), the following equation is obtained: 1 .
The right hand side of (30) is negative, while the left hand side is negative only as long as Ve is smaller than If Ve exceeds this value a singularity appears.Similar exercise can be performed using the last equation of ( 8) and (11).When negative ion density gradient is eliminated from equation (8), equation ( 11) can be, after some rearrangement, written as: The right hand side of (32) is negative, while the left hand side of (32) is negative only as long as: If V2 exceeds this value a singularity appears.To resume: the system (8) -( 12) becomes singular whenever one of the conditions ( 27), (31) or (33) is violated.

Results
For the results shown in figure 1 the following parameters are selected:  = 1/3670.482(deuterium ions),  = 1 (negative ions with the same mass as positive ions),  = 0.02, si = se = 0.1, s2 = 0.The space coordinate is normalized to the ionization length L, so A1 = 1 is selected.The following initial conditions are selected: 0 = 0, Ni0 = 1, 0 = 0, Vi0 = 10 -5 , Ve0 = V20 = 0.In addition, two sets of initial densities of negative particles are selected.For the top plots (a) -(e) Ne0 = 0.2413 and N20 = 0.7587 is taken, while for the bottom graphs (f) -(j) Ne0 = 0.2412 and N20 = 0.7588 is chosen.The difference in the initial conditions is very small, but the difference in the resulting profiles is very obvious.The system (8) -( 12), ( 14) is integrated in the positive X direction.The integration breaks down, when the neutrality condition ( 14  In figure 2 the solutions of the system (8) -( 13) is shown for the same parameters and initial conditions as in figure 1.Additional parameters that must be selected, are  = 310 -5 , 0 = 310 -6 , A2 =  2 .The initial electric field 0 is found from figure 1  the integration is stopped at Xc = 1.2539 and the potential at this point is (Xc) = -1.45993.Both values are selected arbitrarily, but in both cases the integration goes well into the space charge dominated sheath region.In both cases the potential drop is determined by the electrons -the negative particles with larger energy.Another thing is observed.The profiles of the solutions Ni(X), Ne(X), N2(X), Vi(X), Ve(X), V2(X), (X), (X) and (X) exhibit oscillations, which can be best observed in the electric field profiles.These oscillations appear approximately in the range between both peaks of the electric field observed in figure 1 (b).Similar was observed also by Yasserian and Aslaninejad [5].
In plots (d) and (i) velocity profiles Vi(X) and Ve(X) are shown.They both increase and at some X the electron velocity Ve would reach the singularity (31).For the initial conditions Ne0 = 0.2413 and N20 = 0.7587 this would occur approximately at Xc = 1.2575815.The potential there would be (Xc) = -4.2548.For the initial conditions Ne0 = 0.2412 and N20 = 0.7588 the same would occur at Xc = 1.25711 with the corresponding potential (Xc) = -4.25075.But in both cases the integration is stopped earlier, since there is no need to go to so highly negative potentials.Let us return to figure 1, where it is shown, how the solution of the system (8) -( 12), ( 14) jumps from high solution to the low solution, when the initial density of the negative ions N20 is increased above some critical value.This value depends on the temperature of the negative ions .If  is increased the corresponding N20 decreases.This is shown in figure 3 where the following parameters are selected:  = 1/3670.482, = 1, si = se = 0.1, s2 = 0, 0 = 0, Ni0 = 1, 0 = 0, Vi0 = 10 -5 , Ve0 = V20 = 0.
The temperature  of the negative ions is gradually increased and then the initial density N20 of the negative ions is found, where the transition between the low and the high solution of the system (8) -( 12), (14) occurs.As  increases, the corresponding N20 decreases.But when  becomes larger than approximately 0.09 -this means that the electron temperature is less than eleven times larger than the negative ion temperature -a sharp jump between the low and the high solution cannot be observed anymore.
This is illustrated additionally in figure 4. In this figure the sheath edge position XS (plots (a) and (f), the corresponding potentials (XS) (graphs (b) and (g)), particle densities Ni(XS), Ne(XS), N2(XS) (figures The selected parameters are:  = 1/3670.482, = 1, si = se = 0.1, s2 = 0, 0 = 0, Ni0 = 1, 0 = 0, Vi0 = 10 -5 , Ve0 = V20 = 0.For the plots (a) -(e) in the top row  = 0.04 is selected, while for the graphs (f) -(j) in the bottom row  = 0.1 is taken.In the top plots (a) -(e) a very sharp jump of all the quantities at the transition from N20 = 0.739 to N20 = 0.740 can be observed.On the other hand, no sharp transition can be observed in graphs (f) -(j), but it can still be noticed that around N20 = 0.65 the system goes from the solution dominated by the electrons to the solution dominated by the negative ions.

Conclusions
A one-dimensional multi-fluid model was developed and used for the analysis of the plasma wall transition in a plasma containing electrons, singly charged positive ions, and singly charged negative ions with the same mass as the positive ions.Electrons and negative ions are described by continuity and momentum transfer equations, while the positive ions are described by the continuity, momentum transfer, and pressure equations.The collisions where momentum or energy is transferred between different types of particles are neglected.If the problem is analysed at the pre-sheath scale, so that the Poisson equation ( 13) is replaced by the neutrality condition (14), two types of solutions can be clearly distinguished.When the density of low-energy negative ions is below a certain limit, the potential gradient in front of the sheath is determined by the electrons (which have a higher energy).This is the high solution.When the density of the low energy negative ions is above this limit, the pre-sheath potential drop is determined by the negative ions (which have lower energy).This is the low solution.The transition from the low to the high solution is very sharp, provided that the temperature of the negative ions is small enough.If the temperature of the negative ions is too large, a sharp jump between the high and low solution is replaced by a smooth transition between the two regimes.However, if the problem is treated at the full scale (pre-sheath and sheath scale together), so that Poisson's equation ( 13) is used instead of the neutrality condition (14), the separation between the low and high solutions disappears.The profiles of the solutions Ni(X), Ne(X), N2(X), Vi(X), Ve(X), V2(X), (X), (X) and (X) show oscillations and resemble the high rather than the low solution of the system (8) -( 12), (14).Some details of the model were presented.The model equations have singularities that restrict the choice of initial conditions.
At the end some comments about the selection of the source terms si, se, s2 and negative ion temperature  are in order.The negative ion temperature  is in an experiment determined by the method of the plasma production and consequently in the model it is selected as an independent parameter.In atypical low pressure plasma discharges electron temperature is around 1 -2 eV, while the ions have temperature close to the room temperature or a little higher, which corresponds to the values between 0.02 eV up to max 0.1 eV and this corresponds to the selected values of  in our work.
As for the source terms the argumentation is the following.In experimental situations ionization mechanisms can be various.If the main reason are collisions between electrons and neutral particles, the source terms are proportional to the local electron density and therefore space dependent.But in order to keep the model reasonably simple we decided to select only constant source terms.In practice this would mean that the main ionization mechanism is some uniform cosmic or some other radiation which ionizes the gas in our vessel.Since any ionization mechanism must be such that equal number of positive and negative particles is created/annihilated per unit time and per unit volume the source terms must fulfill the relation 2 .
ie s s s =+ The simplest option would be si = se = s2 = 0.But in this case the solutions of the system (8)-( 12), ( 14) and (8)-(13) would not move from initial values.So at least one of the source terms se or s2 must be positive and consequently also si must then be positive.The larger are the source terms, the faster is the increase of the particle densities and fluid velocities and the smaller becomes the value of X, where the electrons and/or negative ions reach the singularity (31) or (33).It turns out (but in the 8 page paper there is no room to discuss this) that the ion source term si has the strongest impact to the length of the system.The larger is si, the shorter is the system.So, the value si = 0.1 is selected in such a way that the length of the system is "appropriate" for the graphs and calculations.With the negative ion mass  = 1 and negative ion temperature  = 0.02 the singularity other hand x is normalized to the ionization length L, ) is not fulfilled anymore.For the results shown in plots (a) -(e) this happens at XS = 1.25293.This point is identified as the sheath edge.At this point the potential is (XS) = -0.71162.At the plots (f) -(g) the neutrality and therefore also the integration of the equations breaks down at XS = 0.859959.Here the potential is (XS) = -0.04154.because of this difference between the absolute values of the sheath edge potentials the solution shown in plots (a) -(e) is called the high solution and the results shown in (f) -(j) are called the low solution.In the first case the pre-sheath potential drop is determined by the electrons and in the second case it is determined by the negative ions, which have much smaller temperature and therefore also the potential drop is considerably smaller.
(b).For the top graphs (a) -(e) the initial values Ne0 = 0.2413 and N20 = 0.7587 are selected and for the bottom plots (f) -(j) the initial values Ne0 = 0.2412 and N20 = 0.7588 are selected.This time the profiles displayed in the top and in the bottom plots are very similar and can hardly be distinguished.For the profiles shown in plots (a) -(e) the integration is stopped at Xc = 1.2543 and the potential at this point is (Xc) = -1.22652,for the graphs shown in plots (f) -(j)
The coefficients A1 and A2 depend on the normalization of the space coordinate x.If x is normalized to the Debye length D,