Modeling and simulations of plasma and sheath edges in warm-ion collision-free discharges

It has been shown recently by Kos et al. [Phys. Plasmas 25, 043509 (2018)] that the common plasma-sheath boundary is characterized by three well defined characteristic points, namely the plasma edge (PE), the sheath edge (SE) and the sonic point. Moreover, it has been shown that the sheath profiles, when properly normalized at the SE, as well as the potential drop in the plasma–sheath transition region (PST), (region between between PE and SE) in collision-free (CF) discharges are rather independent of discharge parameters, such as the plasma source profile, ion temperature and plasma density, providing that the sheath thickness is kept well bellow the plasma length. While these findings were obtained by theoretical means under idealized discharge conditions, the question arises whether and to which extent they are relevant under more complex physical scenarios. As a first step toward answering this question the CF discharge with warm ions is examined in this work via kinetic simulation method in which some of the model assumptions, such as independence of time and the Boltzmann distribution of electrons can hardly be ensured. Special attention is payed to effects of ion creation inside the sheath. It is found that only with considerably increased sheath thickness the sonic point always shifts from SE towards the wall. Whether the absolute value of ion directional velocity at the sonic point will increase or decrease depends on the ion temperature and the source strength inside the sheath. In addition preliminary comparison of results obtained under CF assumption with the representative ones obtained with strongly enhanced Coulomb collisions (CC), indicate the relevancy of hypothesis that the VDF of B&J can be considered as a universal one in future reliable kinetic modeling and solving the plasma boundary and sheath problem in both collisional and collision-free plasmas.It has been shown recently by Kos et al. [Phys. Plasmas 25, 043509 (2018)] that the common plasma-sheath boundary is characterized by three well defined characteristic points, namely the plasma edge (PE), the sheath edge (SE) and the sonic point. Moreover, it has been shown that the sheath profiles, when properly normalized at the SE, as well as the potential drop in the plasma–sheath transition region (PST), (region between between PE and SE) in collision-free (CF) discharges are rather independent of discharge parameters, such as the plasma source profile, ion temperature and plasma density, providing that the sheath thickness is kept well bellow the plasma length. While these findings were obtained by theoretical means under idealized discharge conditions, the question arises whether and to which extent they are relevant under more complex physical scenarios. As a first step toward answering this question the CF discharge with warm ions is examined in this work via kinetic simulation method in which so...


I. INTRODUCTION
Ion velocity distribution functions (VDFs), which are characterized by well defined moments and have second moment comparable to the thermal pressure of electrons or even considerably larger, are of a particular interest in laboratory, fusion and space plasmas.In the scrape-off-layer (SOL) plasmas of Tokamak fusion devices [see, e.g., Refs. 1 and 2) such ion VDFs originate from supposedly Maxwellian ions which penetrate from the core plasma into SOL across the last closed magnetic flux surface (LCFS) during disruptive events, such as edge localized modes (ELMs), see e.g., Refs. 2 and 3].
a Electronic mail: leon.kos@lecad.fs.uni-lj.siWithin the SOL, where the magnetic flux tubes are usually terminated by electrically grounded electrodes, (i.e., limiters or divertors) a self-consistent electric field establishes in both (upstream and downstream) directions aligned with the magnetic field lines (so-called parallel directions).Direction of electric field is towards the terminating surfaces.Electric field has the largest value in the sheath, which is formed between the electrode and the quasi-neutral plasma.This electric field repels electrons away from the boundary electrodes.While even in collision-free (CF) plasmas, where the mean free path for binary processes is much longer than the distance between the plates, the bulk electrons can be modeled as Maxwellian, the motion of ions coming from perpendicular directions into the SOLregion is determined by their initial velocities, by the SOL electric field and, especially by the losses in parallel directions.Because of that in a steady-state the ion VDF strongly deviates from the initial distribution function, which is believed to be close to Maxwellian.For investigating plasma properties it is often convenient to define 2 the temperature and related quantities, such as pressure, ion-sound velocity, heat, energy and viscosity fluxes in terms of fluid quantities, i.e., in analogy to systems with VDFs in thermodynamic equilibrium.With such a definition different ion VDFs that have equal "temperatures" (hereinafter without quotation) can have different direct or feedback effects to plasma parameters, sheath properties and the plasma surface interactions at both microscopic and macroscopic levels.However, it is convenient to study the significance of such possible effects related to particular VDFs with the increased complexity of their shapes, starting from elementary ones, such as monoenergetic beams and (three dimensional in -space) velocity-shells, water-bag, kappa and Maxwellian distributions (see e.g., Refs.4-6), and their combinations, which can be expressed analytically or in a simple numerical form.For example a mono-energetic ion beam (released from e.g., a purposely created local fireball 7 ) around a positive electrode in so-called large-volume plasmas with multi-pole magnetic surface confinement, and accelerated by a strong double layer, after being reflected many times at the magnetized surface back to the main chamber, can establish a mono-energetic threedimensional velocity distribution function.In magnetic electrostatic plasma confinement (MEPC) devices 8 ions with high temperatures, comparable to electron temperatures are produced, while in low pressure double plasma devices, even without surface confinement (e.g., in experiments where they originate primarily from a high potential plasma 9 ) the ion temperature much higher than that of electrons can be realized.
However, in modeling the interaction of ions with a solid material (e.g., a diagnostic probe, electrode or a part of the wall) the parallel and perpendicular part of velocity distribution will both be characterized by the whole spectrum of velocities affecting the surface structure similarly as a three dimensional Maxwellian VDF in both parallel and perpendicular direction with respect to the flow direction.Furthermore, regarding the sheath formation, e.g., the critical (Bohm) velocity, only component of the VDF parallel to the flow and the corresponding parallel ion temperature play a role, while perpendicular components, certainly, enter for calculating higher moments, such as heat flux.For bounded plasmas confined with e.g., permanent magnets distributed over the chamber-wall surface such mono-energetic (three-dimensional) "velocity in shell" distribution function can be a perfect source for warm plasmas.However modeling the ion flow to the wall in collision-free plasmas must be done with the source VDF decomposed into the parallel and perpendicular part with respect to the flow direction, i.e., with the one-dimensional water-bag (parallel and perpendicular) components (see e.g., Ref. 4), each having its own temperature.
By keeping in mind above considerations, it is clear that in one-dimensional collision-free plasma models and kinetic simulations the actual ion VDF will deviate from the source VDF only in the flow direction.For each source VDF a different and unique ion-VDF is expected to be found, such that the complete plasma and sheath equation is satisfied everywhere.To our knowledge solving this problem with the water-bag source VDF has been just tackled in the past in Ref. 10, while most efforts on solving collision-free discharges with warm ions were done with a Maxwellian source [11][12][13][14][15][16][17][18][19][20] by assuming its strength to be either proportional to the electron density or independent of position (flat ion source).In this context it is important to mention that in Ref. 21 an artificial source has been constructed with intention to obtain Maxwellian ion VDF in the plasma center, resulting in the fluid quantities calculated at the plasma boundary which are not far from those obtained in works cited above. 2 On the other hand, in Ref. 22 it has been found that the fluid quantities at the plasma boundary, obtained with the cold ion source (Maxwellian with zero temperature, i.e., Dirac δ-function) in the famous CF Tonks-Langmuir 23 model and collision dominated charge exchange (CX) model also yield similar results, in spite of the fact that ion VDFs from the respective models exhibit none apparent similarities that might be identified.Above findings give hope that the ion VDF emerging from the original Bissel and Johnson (B&J) model 11 with a Maxwellian source, might serve as a reference in all relevant discharges with warm ion sources, yielding the moments, i.e., fluid quantities of a "universal" relevance which are, with a high degree of confidentiality independent of the source VDF.
However, the original B&J model implies intrinsic employment of mathematical two-scale approach 24 constrained to only one free/external parameter of the problem, i.e., to the ion-source temperature T n , with the Debye length disregarded (λ D = 0) and, moreover, numerical solution has been obtained for a few temperatures only, with the particular ion source profile s i proportional to electron (Boltzmann distributed) density n e (s i ∼ n e ).In this context it should be noted at least that the vanishing Debye length implies infinite electric field at the plasma edge, i.e., infinitely thin charged sheath (separating the neutral plasma from equally charged terminating planar surfaces) and, even worse, infinitely high plasma density -apparently contradicts the basic assumption of the model, i.e., about negligible cross sections for particle-particle interactions in a collision-free plasma.On the other hand collision-free plasmas (with very long mean-free paths for collisions) under realistic conditions, e.g., in laboratory and numerical experiments/simulations are characterized by intrinsically non-vanishing and/or externally variable free parameters, such as particle temperatures, density and source (particle production) profiles, as well as by spontaneously establishing ones, such as non-Maxwellian (e.g., truncated) electrons in vicinity of boundaries and time-dependent collective processes.
For above reasons a generalized B&J theoretical model, which takes into account many theoretically/computationally feasible free parameters like e.g., the temperatures, densities, source profiles, etc., such as examined in a series of investigations mentioned above [11][12][13][14][15][16][17][18][19][20]25 must be employed, however, in a manner such that the explicit functional dependencies on relevant plasma quantities on these parameters can be written down. One such functional dependencies are established, one can decide which of the above mentioned parameters is relevant for a particular physical scenario.Unfortunately, as emerges from almost one century long history of investigations of the T&L model (being just a particular, analytically manageable case of the B&J model) even such an apparently trivial discharge still requires considerable efforts for the problem to be closed in its basic aspects, i.e., concerning the most conceivable criterion for identifying the common plasma-sheath boundary and the right physical quantity for properly characterizing it in realistic plasmas, which are characterized by a finite Debye length and a non-negligible ion-source within the sheath.The B&J model with a non-vanishing ionsource temperature is mathematically for an order of magnitude more complex than the T&L one, and for the next order of magnitude more demanding from point of view of parametric dependence of a solution, i.e., a quantity of interest as a function of the ion-source temperature.In present work we update the results of our previous investigations reported in Refs.14, 16-20, and 25, however, here in both tabular and semi-analytic forms accompanied by new representative graphical results, which are obtained for previously unavailable combinations of parameters, primarily those concerning presence of both weak and strong ion source.According to Ref. 20 the main quantities of interest in identifying and characterizing the plasma and sheath edges and finding the correct sonic point, are the moments of ion VDF and the electrostatic pressure and their pseudo-gradients (derivatives over potential) as functions of potential, rather than as functions of position.The basic advantage of this approach is employment of the product of electric field and the Debye length, which turns out to be finite even in two-scale approach, i.e., at points of singularity of electric field alone, while its profile as well as profiles of its pseudo gradients in both plasma and sheath region turn out to be rather independent of the ion-source temperature and the Debye length, providing that the discharge length is properly normalized and the sheath thickness is kept well bellow the plasma length.The particular advantage of this approach is a new natural definition of the plasma and sheath edges.Namely, it is found that common plasma-sheath boundary is characterized by three well defined characteristic points, named the plasma edge (PE), the sheath edge (SE) and the sonic point, the last one formulated in terms of the differential ion polytropic coefficient function (DPCF) 22 in the form of unified Bohm criterion.A remarkable revealing should be pointed out: starting from SE in the wall direction the sheath profiles, as well as the potential drop within the plasma-sheath transition (PST), i.e., between PE and SE, appear quite insensitive in a wide range of above mentioned parameters for both types of ion sources considered in this work.The first one is the so-called exponential ion source, where the source strength is exponential function of potential.This means that it is proportional to the electron density.The second is the constant ion source, where the source strength is a given constant.In the second case ions are created also within the sheath.
The location of the sonic-point and corresponding ion-sound velocity, in our generalized B&J model however, appears to be dependent on all the three free parameters of interest, i.e., the Debye length, source temperature and the ion source profile.The ion source profile has been modeled as proportional to powers of the electron density (∼ n β e ), with β ≥ 0. The unified Bohm criterion has been found based on recent theory of the intermediate plasma-sheath solution for warm-ion plasmas from Ref. 19.In this work numerical results have been presented only for the exponential source ( β = 1), while the case with other profiles, such as constant ionization rate ( β = 0) has been rather unattended.
As the next step towards resolving possible effects of ion-creation inside the sheath to the sonic point location and the sheath profile the generalized B&J discharge is further examined here via kinetic particle-in-cell simulations in which some of the assumptions of the model, such as independence of time and the Boltzmann distribution of electrons can hardly be ensured.In addition it turns out that the essential features and quantities (e.g., related to derivatives of the moments of the ion distribution function and the field-pressure) from theoretical model can be reproduced by kinetic simulations with even better resolution than the resolution, which can be achieved by numerical solutions of the theoretical model.
Besides above mentioned updated theoretical considerations and results and new important quasianalytic expressions, it is found in this work that only at considerably increased sheath thickness the sonic point shifts from SE towards the wall.Absolute value of the ion directional velocity at the sonic point can either decrease or increase, depending on the ion temperature and the source strength inside the sheath.Physical reasons for this are discussed in detail.In addition, it turns out that simulated profiles and derivatives do not indicate any special role of deviation of electron VDF from Maxwellian ones.Maxwellian electron VDFs are traditionally employed in theoretical models, which is mathematically convenient but sometimes unrealistic.
The paper is organized as follows.In Section II A the theoretical model from Ref. 20 is again briefly presented and updated with new theoretical considerations and numerical data, especially those concerning enhanced ion-source in the sheath region.In particular the definitions of plasma edge (PE), sheath edge (SE) and sonic point are re-stated again, however with updated notation and from updated points of view.In Section II B particle in cell (PIC) simulations are described.In Section III theoretical and PIC simulation results are systematically compared under various discharge parameters, with special attention payed to possible effects of ion creation within the sheath, as well as to hypothesis about possible applicability of present model to plasma boundary problems under collision dominated discharge scenarios as well.Summary and discussion of results is given in Section IV.The relevant numerical values of plasma parameters not shown previously in literature are obtained and given in the Table I in the Appendix.

A. Theoretical approach and results
The basic equations of the model are the one-dimensional time-independent kinetic equations for the ion and electron velocity distribution functions (VDFs) f i,e (x, ) and the Poisson equation: to be solved under symmetric boundary conditions between two perfectly absorbing co-planar plates characterized by the electric potential Φ(±L) ≡ Φ W and located at positions x = ±L, under assumption that starting from the symmetry plane (x = 0, Φ = 0) the electrostatic potential Φ(x) monotonically decreases in directions x ≷ 0. Both normalized and unnormalized quantities will be used with the same notation of symbols (with rare exceptions), so that, e.g., the ion and electron densities n i,e = ∫ f i,e d and quantities related to higher velocity moments m i,e = ∫ f i,e m d , such as the directional velocities u i,e = with n 0 = n i (0) = n e (0), T e,i,0 = T e,i (0), e the positive elementary charge, k the Boltzmann constant, c se ≡ (kT e0 /m i )  14,25 or ascribed to an external ion-source originated from perpendicular direction, e.g., in the cases when the model is applied to scrape-off-layer 2 (SOL) plasma in contact with plasma-core of tokamak devices.It should be emphasized that, unlike Bissell and Johnson 11 we strictly distinguish the source temperature T n (with the common subscript "n") from self-consistently established ion temperature T i (ϕ).By introducing ε ≡ λ D /L where λ D = ( 0 kT e0 /n 0 e 2 ) 1/2 is the Debye length, 0 is the "vacuum permeability" and p E (Φ) ≡ ε 2 E 2 /2 is the abbreviation for the electrostatic pressure, the system (1) takes the form where prime denotes the derivative over the potential.Under these conditions the virial V(Φ) , introduced below, is constant. 27,28One should keep in mind that the electric field is finite everywhere even in the limit n i − n e = p E = 0.The normalized form of the total pressure balance therefore reads 2T − 0 E 2 /2 ≡ V(Φ), with 2T ≡ (n i,e kT i,e + n i,e m i,e u 2 i,e ).More precisely V(Φ) turns out to be constant at arbitrary closed surface of a box or cylinder (in present one dimensional geometry) having bases coplanar with the end plates/walls. 20If Boltzmann distribution of electrons is assumed, the non-dimensional virial takes the form and with normalized T e = T e0 ≡ 1 left explicitly for convenience.However, it should be noted that assumption about pure Boltzmann distributed, i.e., Maxwellian electrons is intrinsic idealization which is not consistent with the conservation of particle out-fluxes.Since the wall potential decreases with increased ion temperature, this might affect the virial conservation as well, so considering its behavior needs a detailed inspection, as will be done bellow.Since no further information can be extracted from the Vlasov equation for electrons, the system of Eqs.(1) (i.e., (3)) reduces to: and where K 0 in Eq. ( 5) is the Bessel function of zeroth order, y = 2 /2 and ϕ b stands for boundaries of integration ϕ PE or ϕ W , depending on whether ε is neglected or not.For easier comparison of this work with the previous paper 20 it should be noted that the subscript "PE" replaces the subscript "S" used in Ref. 20 and it refers to plasma edge.Also note that the abbreviation for the eigenvalue B of integral equation, introduced by B&J 11 for expressing equality of the ion and electron outfluxes, is in above equations simply replaced by equivalent quantity B ≡ (L/L i )/ √ 2πT n , (note that T n /T e0 ↔ T n /) i.e., with the ionization length taking the role of a physical eigenvalue. 25In numerical approach after Eq. ( 5) is solved iteratively (with E(ϕ) and formal eigenvalue value B) the wall potential ϕ W is found from particle flux balance, i.e., from B, for a particular set of parameters T n , ε and β (see, e.g., Refs.11, 12, 14, and 25), while L/L i is calculated from: where the integral in parentheses is the electron density n e,av averaged over the entire discharge from x = −L to x = L resulting in a value of n e,a which is always slightly smaller than unity.Note that E(ϕ) depends on parameters T n and β but this is not emphasized explicitly in Eq. (7).With E(ϕ) known, the ion velocity distribution function is calculated in a straightforward manner from Eq. ( 6).Then its moments, such as ion density, directional ion velocity, flux, etc. can be found easily.However, for comparing the theoretical results with those obtained in a numerical, experimental or simulation domain of a normalized physical length L = 1, it is important to remind here that the ion velocity distribution in CF discharges is a function of total particle energy rather than of position.This means that it is independent of the potential profile shape, which, on the other hand, has been shown to depend on source profile (see e.g., Refs. 25 and 29).It can be seen from Eq. ( 6) that a natural normalization of ion VDF, which is independent of L i ( β), must be of the form f i ( β)L i /L.This fact is illustrated for the first time in Fig. 1 where the ion VDFs, obtained numerically in the domain L = 1 for T n = 1 and ε = 0, for values β = 1 (red color) and β = 0 (blue color) at several potentials within plasma and sheath are plotted after being normalized each with L i ( β)f i ( β).It is evident that they are identical to each other at any potential in plasma and sheath regions.
The typical quantities, obtained from numerical solutions of Eqs. ( 5) -( 7) for ε = 0 and values β = 0 and β = 1, are tabulated in the Appendix and presented in Fig. 2 versus temperature T n .In plot (a) ion temperature T i0 at the center of the discharge and ion temperature at the plasma edge T iPE are shown.In graph (b) the potential at the plasma edge Φ PE is presented.Open circles show values obtained from numerical solutions of the system ( 5) - (7), while the solid line shows result obtained from approximate fitting formula ( 8 of exact solutions found from numerical solutions of Eqs. ( 5) -( 7) for ε = 0 is compared with obtained by fitting formula (9), given below.For quick estimation of the plasma edge and sheath edge potentials ϕ PE and ϕ SE the following approximate formulas can be used: 20 The ionization length can be expressed in the form while estimations for the differential ion polytropic coefficient function (DPCF) = 1 + (d ln T )/(d ln n), ion temperature and ion directional velocity are rewritten from Ref. 20 (note that here we replace the symbol γ, introduced in Ref. 22 and used by coauthors throughout their subsequent works, with for avoiding possible mixing with various other quantities used in fusion related texts, such as the sheath heat transmission coefficients 2 ).Approximate formulas therefore read: 105311-8 The value of ion temperature in the center of the discharge T i0 has to be read for a given source temperature T n from the Appendix.The quantity defines the "local ion-sound speed", with known also as the "screening temperature" and the ion DPCF usually calculated via i = 1 + n i T i dT i /dϕ dn i /dϕ .Above considerations are strictly valid only in the original B&J model (ε = 0), i.e., until the sheath thickness is negligible, so that ions produced in it can safely be neglected.According to Ref. 20 a general relation for the ion directional velocity: . The general behavior of these terms is illustrated in Fig. 3 where we show the characteristic points obtained numerically in the generalized T&L model for T n = 0, β = 1 with several finite values of ε.For comparing the cases obtained with β = 1 and with β = 0 we refer to Fig. 2 in Ref. 20.From the profiles in Fig. 3 it appears that (i) up to the point of inflection ϕ PE of the pseudo-gradient d(n i − n e )/dϕ the plasma quasineutrality condition (n i − n e 1) holds with a high degree of reliability, while (ii) up to the point of the pseudo-gradient maximum ϕ SE the electrostatic pressure can be neglected.Since this is a universal rule obeyed independently of external plasma parameters, such as the ion-source temperature T n , smallness parameter ε ≡ λ D /L ∼ λ D /L i and the source profile S i ∼ e −βϕ , the points ϕ PE and ϕ SE have been recognized in Ref. 20 as well formulated, natural definitions of the plasma and sheath edges, while the region ϕ PE − ϕ SE between them has been anticipated as the (weekly nonneutral) transition region (PST) between the quasineutral plasma and the electric field-dominated sheath.The magnitudes of field-pressure related quantities inside the PST in generalized (finite ε) T&L, B&J and fluid models, can be nowadays best estimated quantitatively at the inflection point via adopting the following scalings (see Ref. 20 and references therein): ε 10/9 ε 6/9 ε 2/9 B&J: ε 6/7 ε 4/7 ε −2/7 ε 10/7 ε 6/7 ε 2/7 fluid: ε 4/5 ε 2/5 ε −2/5 ε 6/5 ε 4/5 ε 2/5 (12)   from which it follows, that the conditions for validity of the general expression for the ion directional velocity: E within PST-region always holds for sufficiently small Debye lengths (1 ε ≥ 0).Explicitly, this relation states: FIG. 3. Characteristic points obtained for T n = 0, β = 1 with several finite ε values.Note that d 2 (n i − n e )/dϕ 2 is strongly reduced, i.e., divided by a proper constant factor.
where the term satisfying condition K > 0 20 describes the contribution of the ions originating from the symmetric part of the ion VDF, i.e., those ions which are created between the point of observation and the wall, with zero-velocity ions presented separately.In the case of cold ion sources this expression reduces to contribution of zero velocity ions only. 20In Eq. ( 13) the approximations based on estimations concerning the relevancy of the field-pressure terms are taken into account, so the strict equality sign stands there for convenience rather than the approximate one.
For the profiles illustrated for cold ion-source model in Fig. 3 it turns out, that in the case of non-vanishing n i − n e 1 (quasineutral plasma) the inflection point (PE) of d(n i − n e )/dϕ is rather insensitive to particular values of ε and T n , i.e., coincides with ϕ PE as obtained for strictly neutral plasma.Increase of ε e.g., above 10 −3 , however causes a shift of the inflection point towards the plasma center.The shift is of the order of one tenth of the electron temperature but, according to Fig. 2 in Ref. 20 with increased ion production within the sheath ( β = 0) this shift might appear even at smaller ε.For finite but small ε this requirement might be sensitive to the source profile (S i e −βϕ ), i.e., to value β, while with increased ion production within the sheath this insensitivity reduces, so that, e.g., replacing the exponentially decreasing source ( β = 1) with the flat one ( β = 0) would cause a shift of the inflection point towards the plasma center for about one tenth of the electron temperature already for ε > 10 −3 .For ϕ = ϕ PE the last term in expression E is negligible even when sheath thickness is considerable, while the kinetic term u 2 i,K is essential.This term decreases strongly from center towards the wall, so that in cases when the sheath thickness decreases, or the source strength either vanishes or its strength drops strongly in the region between ϕ = ϕ PE and ϕ = ϕ W , the directional velocity reduces to With increased temperature (c.f., Ref. 20) it appears that the width of the transition region between Φ PE and Φ SE does not depend on either temperature as well (Φ SE − Φ PE ≈ 1/3), as illustrated also in Fig. 3.This behavior, obviously, holds well for sufficiently small Debye length, e.g., ε < 10 −3 .
Definition of the sonic point u 2 i = T * e + i T i ≡ c 2 s requires vanishing of last two terms in Eq. ( 13), i.e., /dϕ, it is clear, that in the limit of infinitely thin sheath, i.e., ε = 0, both these terms vanish independently of each other at the inflection point so that the sonic point coincides with it.This means that ϕ B,ε =0 = ϕ PE , u 2 iB,ε=0 = 1+ iPE T iB ≡ c 2 sB .While all the quantities entering the limit ε = 0 can now be considered as completely elaborated, i.e., known this work and Ref. 20 for arbitrary ion source temperature and profile, situation is mathematically and physically much more complicated with non-vanishing ε.For taking a closer insight into this scenario we analyze the unified Bohm criterion in the form of three equations, i.e., originate from the definition of the ion sound Eq. ( 13) applied at the sonic point (ϕ = ϕ B , u 2 iB = c 2 sB ), while the third one is the expression Eq. ( 4) rewritten under the condition p E 0 with T e = 1.
One can see immediately that the second equation of the unified Bohm criterion Eqs.(15) simply gives the value of d(n i − n e )/dϕ at the sonic point, providing that the ion VDF, i.e., the kinetic term K(ϕ ϕ B ), as well as ϕ B and the directional velocity Γ eB /n eB = Γ iB /n iB = u iB are known.This equation is more important qualitatively than quantitatively.Namely, it states that, when the ionization within the sheath is increased, either by increased ε for a chosen ion source within the sheath ( β = 1), or when ε is kept constant but the ion source within the sheath is increased, ( β → 0).The quantity d(n i − n e )/dϕ > 0 at sonic point, in any case, must shift towards higher values as well.Quantitatively, however, one is interested in possible determination of ϕ B and u iB from first and third relation in Eqs.(15), but this task, obviously, can hardly be accomplished without knowledge of iB and/or iB T iB .
The method of bypassing the above problem is to rely on the scalings (12) for estimating the shift of the potential in the form ∆ϕ = ϕ B − ϕ PE = C ϕ δϕ where explicit form of the factor ϕ PE , which is of the order of unity, has been presented in Ref. 20. Substitution of this expression into third of Eqs.(15) expanded in vicinity of ϕ PE readily yields where ∆T iB is change in ion temperature inside the region ∆ϕ and terms in brackets can be replaced with the approximate expressions Eqs.(10) resulting in: It turns out that for sufficiently high T n 3 the coefficient C B = [1.266Ti0 + 2]C ϕ depends quite weekly on T n , i.e., C B ≈ 3 for β = 1 and drops for about 20% for β = 0, due to factor e 2β 7 ϕ PE .With decreased T n , however, C B shifts towards a doubled value.Moreover, it has been found that the numerically obtained ∆u 2 iB and ∆ϕ in Ref. 20 can be well fitted with theoretical predictions provided that C B ≈ 3 is multiplied by a correction factor of the order of 2. This quantitative discrepancy is not surprising with having in mind that presented method of "bypassing" the problem of exact determination of all relevant quantities at the sonic point is intrinsic estimation, and that, moreover, the derivation of C ϕ (T n , β) has been made on the intermediate scale theory from Ref. 19, where analytic considerations have been shown to hold for rather high ion-source temperatures (coincidentally T n 3) which still has not only to be upgraded to lower ones but also better justified.Secondly, the pressure balance equation employed in present "bypassing" of the problem (third of Eqs. ( 15)), as well as some other assumptions related to self-establishing ion VDF in the generalized B&J model, is strictly valid until the ion VDF is one-dimensional in phase-space.Last, but not least, it does not seem plausible to expect that characteristic plasma-sheath points, i.e., ϕ PE ϕ PS and ϕ B , and the relevant quantities therein (such as u 2 iB , ∆T iB and i remain insensitive to increased ion production, i.e., the ion VDF shape therein, at least not for any ion source temperature.For investigating these effects we apply a more powerful and more realistic method, i.e., PIC simulations, as follows.

B. Simulation approach
Similarly as in some previous works of co-authors addressed to particular problem of determining the ion DPFCs in T&L and B&J models 16,30 for particle-in-cell (PIC) simulation [31][32][33][34] we use the one-dimensional (1D3v) BIT1 35 code.This code was designed primarily for fusion oriented simulations [35][36][37] in the SOL region with possibilities to add/subtract a variety of atomic and plasmawall interaction microscopic processes with the basic capability to maintain a desired population of particles Maxwellian and cut-off Maxwellian, as measured in experiments (and as such assumed also in theoretical models) even when not expected, i.e., in collision-free (low pressure) plasmas.The unexpected (and still not fully understood) electron local thermodynamic equilibrium, known as the Langmuir paradox (see e.g., Ref. 38) is in BIT1 achieved with an artificial Coulomb collision mechanism, which turns out to have excellent performances (see Refs. 16, 30, and 39).In many other PIC codes under such conditions the source Maxwellian electrons behave non-locally i.e., they quickly leave the system so that it is very difficult to reach steady state at all.Next advantage of the present code is, of course, its capability of performing interactive simulations and obtaining a large number of results during relatively short times when the code is run in the parallel mode.However, a special care should be taken in preparing the simulation parameters and justifying their relevancy for particular physical scenarios, before the results from final runs can be considered as definite and further processed and interpreted, as follows.
In the code the distributions functions of the particles that are either created by the volume source or injected from the walls are given in terms of parallel T || and perpendicular T ⊥ temperature with respect to the external magnetic field.In present simulations a small magnetic field (B ∼ 10 −4 ) Tesla normal to both plates (in direction of x axis) is introduced, which does not affect either of distributions and fluxes in the x direction, The length of the system is 5 cm and the system is divided into 12000 cells.The length of one cell is therefore 4.167 × 10 −6 m.The surface of the boundary electrodes is 10 −4 m 2 , so the volume of the system is 5 × 10 −6 m 3 .Because we wish to get simulation results, that correspond to different Debye lengths and consequently to different values of ε, the source strength is varied.But typical order of magnitude is 10 22 electron ion pairs produced per second and per m 3 .With such source the plasma density is typically of the order between 10 15 and 10 16 m −3 .The parallel and the perpendicular electron temperatures are both set to kT e = 1 eV.So the electron Debye lengths are between 2.3 × 10 −4 and 7.4 × 10 −5 m.The length of the system is therefore between 210 and more than 670 Debye lengths.So the number of cells per Debye length is between 17 and 56.The electron plasma frequency is between f pe = ω pe 2π = 2.8 × 10 8 s −1 and f pe = 8.9 × 10 8 s −1 .Time step ∆t = 10 −11 s is selected.This results in more than 100 time steps per electron plasma period even for the largest plasma densities obtained in the simulation.Since we had rather powerful computational resources at our disposal, a rather small number of physical particles per computer superparticle was always selected.This number never exceeded 10 4 so we usually had up to 10 7 super-particles in the system.Since the potential, density, temperature and other profiles are obtained as an average over a large number of time steps (2 19 = 524288), the number of physical particles computer particle is not a crucial parameter and very smooth profiles of density, potential, temperature etc. are always obtained.
If number of physical particles per computer particle is very large, plasma potential oscillates with very large amplitude, which sometimes have frequency that corresponds to the local electron plasma frequency (see e.g., Ref. 30).If superparticles are increased even further, oscillations with rather strange spectra are sometimes observed.The time-averaged potential profiles, nevertheless, can still be identical to theoretically predicted ones, but on the other hand fine details of simulated VDFs are lost.As will be demonstrated bellow in present work the simulated VDFs perfectly fit theoretical ones (unlike referred to simulations from Ref. 30) and this is achieved by decreasing the size/charge of super-particles with a trial and error method for each density, until such oscillations become insignificant.The price that is paid for usage of smaller superparticles are increased simulation run-times and data storage resources.
It is interesting to note that the wavelength of the basic mode of plasma potential oscillations is 2(L − l sh ) rather than 2L.Here l sh is the apparent sheath thickness.This has been confirmed by observing the instant potential profiles from one time-step to another.As already mentioned these oscillations appear, when large superparticles are used.Another way to excite these oscillations is to inject positive and negative particles into the system with slightly unbalanced rates.The sheath region turns out to be time-independent even when the amplitude of the potential oscillations at the center and elsewhere in plasma region is extremely high, even larger than the sheath potential drop.The "sheath potential drop" is potential drop between the wall and the apparent sheath edge at a distance l sh from the wall, which does not change its position and value of the potential, as the simulation progresses.

A. Non-isothermal hot ion source T n > T e
The simulation method and results for cold and warm ion sources with β = 1 and β = 0 have already been presented in much detail elsewhere, e.g., for T n = 0 in Ref. 30 and for T n = T e = 1 in Ref. 16.For qualitative and quantitative comparisons with theoretical results obtained for hot ion sources and presented extensively in Ref. 20 we start here with the "familiar" temperature T n = 7, which has been exhaustively examined in Ref. 20.In Fig. 4 the virial V(ϕ) , given by formula (4) is plotted for 2 values of ε.Solid lines refer to ε = 3.1 × 10 −3 and dotted lines refer to ε = 10 −3 .Virial found from the numerical solution of the theoretical model is plotted with violet solid and dotted lines.It can be seen that it is constant.From the formula (4), using T e0 = 1 and T i0 = 2.18, one gets V = T i0 + T e0 ≈ 3.18.The value T i0 = 2.18 is read from the Table I in the Appendix at T n = 7.The virial found from the PIC simulation is shown by black solid and dotted lines.Also in this case the virial is constant.It should be emphasized that the lines obtained from theoretical model on one hand and from PIC simulation on the other can hardly be distinguished in spite of the fact that theoretical and simulation curves have been obtained with different source profiles ( β = 1 and β = 0, respectively).Small shift appearing between two types of curves just indicates that ion temperatures obtained by two methods with various ε can hardly be expected to perfectly match each other.In the same figure theoretical and simulated kinetic energy densities are compared also.Theoretical results are plotted with blue solid and dotted lines, while simulation results are shown with black lines.Lines corresponding to the virial and to kinetic energy densities can be seen in the upper part of the plot.In Fig. 4 also a more detailed decomposition of V into relevant contributing terms is illustrated, but of course only for theoretical V.
In Fig. 5 the properties of ∫ ϕ 0 (n i − n e )dϕ = ε 2 E(ϕ) 2 /2 and its derivatives, found from PIC simulation (green and blue curves) and from theoretical model 20 (red curve), are inspected.Red vertical dotted arrows mark the positions of PE, SE and the wall as found from the theoretical model for ε = 0.The plasma edge and the sheath edge are located at ϕ PE = 0.42 and ϕ SE = 0.7, respectively, while the wall is located at ϕ W = 2.77.The PE corresponds to the inflection point of d(n i − n e )/dϕ, while SE corresponds to the maximum of d(n i − n e )/dϕ.Green curves show results obtained from PIC simulation for ε = 10 −3 , while blue curves correspond to ε = 3.5 × 10 −3 The width of the PST region is approximately ∼ 0.28.As it can be seen from comparison of red, green and blue curves that this width is rather insensitive to ε values.Note that it is possible to calculate theoretical curves even for potentials that exceed the wall potential, ϕ W = 2.77.This is due to the fact that (starting from PE) the ion density profile in PST and sheath region can be calculated from the asymmetric (collision-free) part of ion energy distribution function simply shifted along the sheath potential.One can see that the simulated electrostatic pressure ε 2 E(ϕ) 2 /2 = ∫ ϕ 0 (n i − n e )dϕ and the charge density n i − n e = d(ε 2 E(ϕ) 2 /2)/dϕ in PST and the sheath region (i.e., right from ϕ PE ) slightly deviate from their theoretical counterparts.This is not surprising having in mind that the theoretical assumption about Boltzmann electrons is not fulfilled in PIC simulations.This is mainly due to two reasons.The first is that the electron VDF in the PIC simulation has a cutoff.The second reason is imperfect electron thermalisation in the kinetic code.In any case equality of the simulated ion and electron fluxes must be satisfied at any point, while in present theoretical model this is not the case.The reason is that in the model exact Boltzmann distribution of electrons is assumed, rather than the one corresponding to cutoff Maxwellian, see e.g., Ref. 40.In order to estimate whether this deviation of electrons from the Boltzmann distribution is a critical issue or not, several quantities are analyzed in Fig. 6.The electron and ion densities are shown in logarithmic scale in Fig. 6(a).From this figure it can be concluded that the deviation of electrons from Maxwellian is not dramatic even when one approaches very close to the wall.In plot (b) the ion temperature and the differential ion polytropic coefficient function are presented.It is obvious that effects to ion temperature can not be identified at all, while the effects to the ion polytropic coefficient can be considerable only close to the wall.Effects of the deviation of electrons from the Boltzmann distribution are more pronounced for smaller simulated particle densities but the effects, of whatever origin they are, to the sonic point seem to be irrelevant.This is illustrated in Fig. 6(c) It can be seen clearly that when ε is increased, the sonic point (marked at cross-sections of squared ion directional velocities and ion sound speeds obtained for two plasma densities) shifts towards the wall, while the value of the ion directional velocity increases, as predicted by theory.Expressions (10) for calculating u 2 iB and ∆ϕ in this case yield nice quantitative agreement with simulation results.With decreased source temperature T n the behavior described above still holds.In Fig. 7. ion VDFs obtained from PIC simulation for T n = 3 are plotted (blue curves) and compared with the theoretical ones (red lines) obtained from system Eqs.(5)-( 6) with ε = 0.It is obvious that, unlike examples from Ref. 16 the present ion VDFs, which are obtained with a large numbers of cells and a high number of super-particles, perfectly fit the theoretical ones.For convenience we insert the theoretical and simulation potential profiles with marked points of observations of VDFs.In simulations these points have been chosen at various positions and the corresponding potentials have been found after simulations were finished and the time averaged potential profiles have been plotted.After that the theoretical VDFs have been found at these potentials (rather than at the same positions) so that comparison of simulation results (obtained with arbitrary finite ε) is consistent with the theoretical reference obtained for ε = 0.In spite of the fact that the potential profiles are intrinsically different for ε 0 and ε = 0 such comparison of VDFs taken at the same potentials is then consistent.
Another purpose of the inserted graph with potential profile in Fig. 7 is to discuss the effects of ion source in the sheath.It is obvious that the sheath in Fig. 7 for ε = 6.9 × 10 −3 is relatively thick and that simulated VDFs are always characterized by slow and negative-velocity ions while these ions are completely absent in theoretical VDFs obtained with ε = 0.Although the contribution of the slow and negative-velocity ions to the total ion density in the simulations for T n = 3 is non-negligible, there is no evidence that such (gradually decreasing) contribution to ion density can affect the locations of the PE, SE and sonic point substantially.By recalling in mind that under common model and simulation assumptions about geometric symmetry of the discharge and symmetry of the ion source in -space, it is clear that the contribution of ions originated from the sheath (symmetric part of VDF) can be safely subtracted in evaluating d(n i − n e )/dϕ, so that one is left with the asymmetric part of the VDF there.This is, in fact, equivalent to considering the "idealized" VDFs i.e., the ones obtained theoretically for ε = 0 (red curves in Fig. 7).
Note that presented simulations are performed with the source, which is constant everywhere, including the sheath region ( β = 0).If a source that decreases strongly in the sheath region (e.g., β = 1) was used, the presence of symmetric part of the VDF for finite ε in PST and sheath region, would be even less significant, as it has been shown in Ref. 20.This means that, until the ion-source temperature T n is considerably above the electron temperature, the simulations and theoretical results will always fit nicely to one another.

B. Decreased source temperature -Case T n = 1
With source temperature T n decreasing, the plasma edge potential shifts towards the classical T&L value (ϕ PE = 0, 854. ..), so the plasma density there decreases in accordance to exp(−ϕ PE ).In discharges where the ion source profile decreases exponentially in the direction towards the wall ( β = 1), the contribution of ions with negative and zero velocities to ion total density should not depend on ϕ PE (and thus on T n ) at all, while their possible influence to characteristic plasma-sheath points, i.e., ϕ PE ϕ PS and ϕ B and the relevant quantities therein (such as u 2 iB , ∆ϕ, ∆T iB and i ), can be considered as negligible, as soon as the sheath thickness, is thin, i.e., ε sufficiently small.For above reasons it is, in general, hard to expect that in discharges characterized with the same ε but with different ion source profiles the effects of ions originated from the sheath will be the same.As indicated in Fig. 2, obtained for a high source temperature (T n = 7) with the "flat" ionization profile ( β = 0) the shifts of u 2 iB and ∆ϕ appear similar but a little smaller than found in Ref. 20 for the same temperature with β = 1, however, without any considerable effect to locations of inflection and maximum points of charge pseudo-gradients, i.e., ϕ PE and ϕ SE .On the other hand, in the limit of "cold" ion-source scenarios the value of plasma edge ϕ PE turns out to be considerably less insensitive to the source profile (c.f., Fig. 2 in Ref. 20), while the ion-polytropic coefficients i in the transition region and its vicinity (in plasma and the sheath) drop much faster with increased ε for β = 0 (and even become negative) than for β = 1 (see Ref. 30).While in cold ion-source model these effects to the ion sound velocity u 2 iB = + iB T iB ≡ c 2 sB might not be of special importance due to smallness of T iB ≈ 0.04 this might be relevant for the scenarios with the intermediate ion temperatures, such as those resulted with the source temperature T n = 1 that in the limit ε = 0 yields some "intermediate" ion temperature profile, i.e., T i0 = 0.421 and T iPE = 0.149 at plasma center and edges, respectively, with the plasma potential drop ϕ PE = 0.625 (c.f., Table I of Appendix A in Ref. 30).
Simulations T n = T e = 1 have been performed under conditions similar to those in Ref. 16 for several plasma densities, however, here with a much lower number of particles per super-particle (from 10 2 to 10 4 , depending on the density) to decrease the level of oscillations and for ensuring a high resolution in phase space.As an illustration potential profiles for several ε are presented in Fig. 8.In addition in the inserted graph 2 examples of the ion VDFs are presented but only for the largest and the smallest plasma density, which correspond to ε 1 = 7.58 × 10 −4 and to ε 4 = 2.67 × 10 −2 , respectively.For each selected density the distribution functions are shown at 2 positions.These positions are marked by numbers.It can be seen immediately that the ion VDF, even when obtained very deep in the sheath (e.g., number 4), is characterized with the long "tail" of ions which penetrate into the plasma region.On the other hand no considerable tail can be detected in the VDF for ε 1 located in vicinity of ϕ PE = 0.625, i.e., at position 2 (ϕ = 0.656).But on the other hand for ε 4 , even at position 3 (ϕ = 0.936) the VDF exhibits a tail, which contains a considerable number of ions moving to the "left" -i.e. away from the wall.So the contribution of the symmetric part of the VDF can have considerable effect in decreasing the directional ion velocity, and thus the sonic point as well.
In order to investigate described effects of constant source profile in more details additional profiles of the relevant quantities obtained from simulations with T n = 1 are plotted in Fig. 9 example obtained with extremely small density (n 0 = 3.1 × 10 13 m −3 ) is included as well, just to illustrate that it is out of range of the common rules that hold for small ε.Namely, from the profiles of temperatures, polytropic coefficients directional and ion sound velocities plotted for ε 1 and ε 2 it is clear that ϕ PE and ϕ SE are rather insensitive to ε, temperature and the source profile.This holds relatively well even for a rather high value of ε, like ε 3 = 7.86 × 10 −3 .However, while the sonic point still shifts towards the wall, as in the case of high ion temperatures and β = 1, this shift is here rather insignificant.More remarkable is the fact that the value of the ion directional velocity a that point, on contrary to discharges with high source temperatures and β = 1, decreases with increased ε.On one hand this means that the usage of expressions (10) for calculating ∆u 2 iB and ∆ϕ in this case is unjustified.But one should note that, unlike discharges with higher source-temperatures, where the square of the ion-velocity profile is almost independent of ε, the u 2 i profiles in Fig. 9 are quite dependent on ε.Physically, that means that the ions originated from the sheath obstructs the "cooling" of ions in the wall direction.
Nevertheless, it appears from above considerations that obtained shifts in quantities of interest are so small that they may both be disregarded even for relatively high values of ε, such as ε 3 = 7.86 × 10 −3 .The insignificance of these shifts, unexpectedly, appears to be even stronger for flat source ( β = 0) than for exponentially decreasing one ( β = 1) and indicate that values ϕ PE ≈ 0.625 and u 2 iB = c 2 s ≈ 1.5 obtained from the theoretical model with ε = 0, can be regarded as universal ones for any realistic collision-free discharge (ε sufficiently small, e.g., bellow 10 −3 ) with ion source temperatures being not far from T n .

C. Effects of departure of ion VDF from theoretical one
Finally we quickly check out the hypothesis about possible employment of the theoretical, i.e., B&J, ion VDF to realistic physical scenarios, when the experimental VDFs considerably departure from B&J one.For this purpose we perform the preliminary investigations of discharge with T n = 3, with the VDFs already presented in Fig. 7 for collision-free (CF) scenario.
In PIC Fig. 10(a) the results from CC simulations are replotted together with results obtained with strongly enhanced Coulomb collisions, and with the isotropic ion-source (perpendicular temperatures T ⊥ equal to the parallel one T ).The relevant temperature profiles corresponding to reference (theoretical) theoretical model (ε = 0) and the corresponding electrostatic terms are presented in Fig. 10(b), together with results obtained in simulations that are performed for a single, rather low density (ε = 2.7 × 10 −3 ).
There are several observations that we find important to be mentioned here as follows.Namely, in spite of intentionally strong Coulomb collisions no isotropy of the ion VDF and the temperature has been achieved in the ion-flow direction.Instead of a possible temperature expected in each direction to be, at least in the center of discharge close to (2T ⊥ + T )/3 ≈ (6 + 1.08)/3 = 2.36, its value in the flow direction is there just slightly above the one obtained under CF scenario.Inspecting the ion VDF shape (not presented here) shows that in direction of the flow it is just bell-shaped but far from any equilibrium (Maxwellian) one.Finally there is no relevant indication that either the relevant fluid quantities nor plasma and sheath edges and the sheath profiles are sensitive to a strong increase of collisions.These results, shortly speaking, at least initially, justify our hypothesis that the VDF of B&J can be considered as archetypal one, such that is appropriate for future detailed modeling and solving the plasma boundary and sheath problem, rather than employing apparently more physically plausible ones such as a shifted Maxwellian (see e.g., Ref. 41).In fact, the present finding just reflects the fact that the ion VDF in vicinity of a "perfectly absorbing" boundary establishes as a non-local, rather than the local one.

IV. SUMMARY AND DISCUSSION
According to our experience solving the collision-free discharge numerically with a high accuracy is an extremely stiff and CPU expensive task.As it can be seen from Ref. 20 that even when done very carefully, the differentiated profiles obtained on discrete grids are far from being perfect.The kinetic PIC simulation which runs in parallel on a large number of processors, on the other hand, turned out to be faster than available codes for numerical solution of the CF model with warm ions.Post-processing of the results of the parallel PIC code turned out to be simpler and even more reliable than the results obtained from earlier single processor version.In this work the parallel PIC code has been used to reproduce previous results (see Refs. 16 and 30) but with much better resolution in phase space due to largely increased number of cells and superparticles in the system.Also new results have been obtained, with special attention dedicated to the isothermal source (T n = T e ) with a constant profile (so-called β = 0, or flat source).In addition in this work the basic theoretical and numerical results of Bissel and Johnson model have been updated, especially the ones concerning the microscopic and fluid quantities and relations between them in the direction parallel to the flow and in the limit of vanishing ε.In kinetic simulations it is ensured, via employing an artificial electron thermalisation, that electron distribution function remains isotropic Maxwellian with Boltzmann distributed density in the longitudinal direction.This is in perfect agreement with the model assumptions, that the electron temperature is equal in all directions -parallel and perpendicular to the flow.The ions, which in both, model and simulations, are fully collision-free, are created by a Maxwellian source.In the direction parallel to the flow they move freely.Their temperature in perpendicular direction remains equal to initial one and it is independent of position.On the other hand the ion VDF, which is established in the parallel direction deviates strongly from Maxwellian at each point of the discharge and exhibits a strong "cooling" with a considerable temperature drop, with respect to the corresponding parallel source temperature.This is demonstrated by the typical values at the center of the discharge and at the plasma edge, presented in Table I.Perpendicular temperature does not enter the problem of identifying and characterizing the planar discharge and sheath at all.In fact, the choice of perpendicular temperature depends on physical scenario of interest, for which no data are provided here.In simulations we have first checked the results (and the code) via simulation runs performed with sources having zero and finite (equal to parallel) perpendicular temperatures.In fact knowledge of perpendicular temperature is relevant for calculating higher moments of the VDF and related quantities of interest, such as heat and energy fluxes.The ion VDF in the present model turns out to be a function of their energy.In the limit of vanishing ε it does not depend on ionization profile at all, provided it is normalized to the ionization length (rather than to the system length).That means that the moments of VDFs, their pseudo-gradients and related quantities, such as ion DPCF and the ion-sound speed are universally relevant (as reference ones) for each β of interest.By pseudo-gradient it is meant that the derivation is performed over the potential ϕ and not over the space coordinate x.But in the rest of this section the word gradient will be used instead of pseudo-gradient and it always refers to derivation over ϕ.The closure of fluid equations describing the ion motion in longitudinal direction is provided by the fact that ion temperature and its gradient are known at any point of the discharge.The plasma edge appears to coincide with the inflection point of the ion temperature and related quantities, such as DPCF and the corresponding local ion sound speed.In the limit ε → 0 this point coincides with the electric field singularity at ϕ PE as well as with the inflection point of charge density gradient which, for ε → 0 turns out to be continuous but non-differentiable function at plasma-side of ϕ PE .For finite but small values of ε, such as they are in all quasi-neutral plasmas, the charge density gradient exhibits several remarkable properties.The first is that for a given ion temperature the potential at which the inflection point of charge density gradient occurs is insensitive to particular ε.In fact the inflection point is located at ϕ PE (ε = 0), while the charge density itself, at the same time still, remains negligible there.The second property is appearance of the maximum of charge density derivative at the potential ϕ SE , which is such that ϕ SE − ϕ PE ∼ 1/3.In addition the value of charge density gradient at ϕ SE is quite insensitive not only to ε but also to ion-source temperature.At ϕ SE the violation of quasineutrality is already considerable, but the electrostatic pressure (the essential boundary condition for solving the sheath equation) appears to be still negligible there.Above properties suggest that separate definitions of plasma edge at the potential ϕ PE , sheath edge at the potential ϕ SE and the plasma-sheath transition region (PST) between them in the potential range ϕ SE − ϕ PE are physically very well founded and make perfect sense.
The third remarkable property emerges from apparent independence of charge density derivatives between ϕ SE and ϕ W of all the "external" discharge parameters, i.e., T n , ε and β.This has already been demonstrated in Ref. 20 for β = 1.Eq. ( 30) in Ref. 20 is referred to as the universal collision-free analytic sheath-wall asymptote.Figs. 5 and 9-d presented here for β = 0 as well as Fig 2 in Ref. 20 (the latter obtained for both β = 0 and β = 1) with temperatures T n = 7, 1 and 0, respectively, suggest that the hypothesis of the universal sheath-wall asymptote is a plausible one.Perhaps it will be possible to prove this hypothesis in the near future providing that a more reliable and universally applicable analytic or semi-analytic solution than available at present 17,19,42 is found for the ion VDF.
In context of sheath properties time independence of the sheath profile, observed in the PIC simulations, during large-amplitude potential oscillations in the quasi-neutral plasma region, irrespectively of nature of their triggering, is the fourth important property of the sheath boundary ϕ SE and the electrostatic pressure in the sheath.This issue deserves particular attention in forthcoming investigations, especially the ones related to sudden bursts of plasma particles, such as ELMs and related parallel plasma transport phenomena.
The ion directional velocity and the ion temperature at the plasma edge, appear to be the most important quantities to be determined for calculating related quantities of interest, such as heat fluxes.In this work it was demonstrated that in high density plasmas (negligible ε) the sonic point coincides with ϕ PE , only if the ion sound speed is expressed in terms of local ion temperature and its gradients (note that in case of non-Maxwellian electrons the local electron temperature and respective gradients must be known as well).The effects of non-vanishing ε, which are formulated with the unified Bohm criterion Eq. ( 13), i.e., the one that takes into account the electrostatic pressure and kinetic effects related to those ions which originate from the region between the point of observation and the wall.The effects of these ions is manifested as a slight shift of sonic point from ϕ PE towards the wall as ε is increased from zero to finite values, while the corresponding shift in the value of the Bohm velocity and ion directional velocity is observed towards either higher or lower values -depending on the ion-source temperature, but manifesting themselves through the corresponding local ion temperature and the source strength and profile within the sheath.In cases of ion sources with temperatures several times above the electron temperature (both parallel) the Bohm velocity increases independently of the source profile, as predicted by Eqs.(16), while in cases when the ion source temperature is comparable to or even smaller than the electron temperature, this behavior is reversed, however, only for "flat" ( β = 0) ion sources.These effects might manifest themselves through the term κ i T i .However, in PST region this term appears to become less significant as the ion temperature decreases.One can confirm this assumption by comparing the ion-sound speed profiles, presented for T n = 7 (T PE = 1.165) and T n = 1 (T PE = 0.146) in Fig. 6-(c) and Fig. 9-(c), respectively.As opposed to the ion sound speed, the directional ion velocity for small ion temperatures is rather sensitive to ε.This can be seen in the same figures.Ions originating from the sheath with negative velocities decrease the directional ion velocity in the whole plasma region.This effect is illustrated by simulations performed with very low plasma density (3.1 × 10 13 m −3 ).It is interesting to note that only in the case of a thick sheath such as the one corresponding to this density (c.f., Fig. 8 For plasmas with more realistic ε-values, it appears that discussed shifts are small enough to be considered as numerically and experimentally insignificant.In e.g., fusion-relevant and many other experimental plasmas, where ε is finite but small, the shifts of ϕ and u i relatively to the reference model (obtained with ε = 0) might be considered to be of the same order as the diagnostic errors or even smaller.Besides the question of the Bohm velocity, even if it is exactly known, the question arises about determining its location ϕ B .To do this one can simply use theoretical data, such as the data presented in Table I, but it must be remembered that the intrinsic information relating any directional ion velocity with the potential is already contained in conservation of the total (kinetic and electrostatic) energy density Eq. ( 4), named the virial.Unlike the conservation of energy for ions Eq. ( 13) (the unified Bohm relation), which is intrinsically approximate but well satisfied for small ε up to ϕ SE , the Eq. ( 4) is demonstrated here to be strictly satisfied at any location and for any ε.Moreover, with neglected electrostatic pressure and a spatially constant electron temperature (T e = T e0 = 1) Eq. ( 4) reduces to a very simple relation connecting the ion directional velocity and the plasma potential, that holds with a high reliability also right from the sheath edge (e.g., up to ϕ SE + 1/2).
It is a long-standing 26 practice to employ the total kinetic energy density Eq. ( 4) for determining the plasma-edge or the sonic-point potential.Alternatively, this can be done if the ion temperature drop with respect to plasma center and the ion-directional velocity at the point of observation are both known (see e.g., Ref. 20 and references therein).On the other hand, finding the Bohm velocity (e.g., in fluid simulations and experimental plasmas) requires knowledge about that location/potential, usually based on definitions of the (intermediate) electric field value/scaling and, moreover, by employing the standard isothermal/adiabatic constants, to be determined via various kinetic approaches and methods (see e.g., Refs.43-46).The history of plasma physics and diagnostics teaches us that looking for values and/or expressions such as unified and/or generalized Bohm criteria (see e.g., Ref. 20 and references therein) aiming for precise quantitative characterization of the common plasma-sheath boundary as the sonic point is not only demanding but physically highly disputable matter.However, based on present considerations one may take another point of view, i.e., to first identify a convenient plasma potential and with this (apparently arbitrary) input to simply employ the pressure balance Eq. ( 4) for calculating the directional velocity at that point.The term "convenient", however, does not give too much freedom since, according to present investigation any arbitrariness is restricted to a rather narrow PST-region, or eventually deeper into the sheath, but starting from ϕ PE , which is already a well tabulated quantity in the present model.Namely, the central question of interest in many experiments, theory and fluid simulations is, in fact, to identify and characterize a proper location where the ion directional velocity and/or kinetic energy has to be calculated, such that, between that point and the wall, the moments of ion and electron VDFs are unaffected by possible microscopic volume processes with accompanied gains and losses.So, e.g., in plasmas with finite ε the reference potential for calculating such a proper u 2 i , in discharges with β = 1 should be closer to ϕ PE than in discharges with β = 0, but the error caused by a particular choice of ϕ (and corresponding u 2 i ) is expected to be negligible until ε is reasonably small.In such cases the ion directional velocity found by this approach still remains close to the sonic one, calculated in terms of DPCF as tabulated, e.g., in Ref. 20.In the cases of extremely thick sheaths (which are more of academic than of practical interest), with ionization present in the entire sheath (such as that one simulated with very low plasma density illustrated in Figs. 8 and 9) it is hard to say, which location is the proper one, in spite of the fact that the sonic point and the Bohm velocity have been found exactly.Here it must be noted that even in this case the charge density gradient (Fig. 9) does not deviate from the regular ones (obtained for small ε), even relatively deep in the sheath, meaning that the charge-densities themselves differ from each other only for an additive constant that has been lost during normalization of electron and ion densities to unity.
Finally, we consider the hypothesis that the VDF of B&J can be considered as a universal one, such that is appropriate for reliable kinetic modeling and solving the plasma boundary and sheath problem with well formulated boundary conditions in collision dominated plasmas as well.However this task requires the analytic expression for the ion VDF that would be more reliable than currently available ones.

APPENDIX: NUMERICAL PLASMA PARAMETERS
In this Appendix numerical values of various plasma parameters, as obtained from the theoretical model are presented.Table I can be used for quick estimation of relevant plasma parameters using the formulas ( 8) -( 16) or for comparison with PIC simulations.The presented quantities are: ion source temperature T n , plasma edge potential ϕ PE , ion temperature in center of the discharge T i0 , ion temperature at the plasma edge T iPE , wall potential for exponential source profile, ϕ W ( β = 1), TABLE I. Source temperature T n , plasma edge potential ϕ PE , ion temperature in center of the discharge T i0 and at the plasma edge T iPE , wall potential ϕ W (β = 1), ionization length L i, β =1 , wall potential ϕ W (β = 0) and ionization length L i, β =0 .Note that here the ionization length is formulated in accordance to Riemann and it is by a factor √ 2 larger than the ionization length defined according to Harrison and Thompson.See Ref. 25
simplified for convenience (without essential loss of generality) with the assumption of Boltzmann distributed electrons n e = exp(−ϕ) and T * e = − 1 n e dn e dϕ −1 = 1.The first two of Eqs.(15), obviously,

FIG. 4 .
FIG.4.The exact (numerical) virial (solid and dashed violet lines) and the double kinetic energy density (solid and dashed blue lines) in comparison with corresponding quantities obtained from kinetic simulations (solid and dotted black lines) for T n = 7.Note that the decomposition of V(ϕ) is done only for numerical results.

FIG. 5 .
FIG. 5. Simulated electrostatic pressure∫ ϕ 0 (n i − n e )dϕ = ε 2 E(ϕ)2 /2 and its derivatives obtained for two densities for T n = 7, β = 0 in comparison with theoretical results obtained for ε = 0.The PE and SE (inflection and maximum points of d(n i − n e )/dϕ, marked with red vertical dotted lines) are found from theoretical derivatives at ϕ PE = 0.42 and ϕ SE = 0.7, respectively.The PST region (∼ 0.28) is, obviously quite insensitive to ε.The theoretical position of the wall (ϕ W = 2.77) is marked with dot-dashed red vertical line.

FIG. 6 .FIG. 7 .
FIG. 6.In plot (a) electron and ion densities n e (ϕ) and n i (ϕ) are shown in logarithmic scale.Red lines show theoretical results, while green and blue lines present the results of simulations.In figure (b) the ion temperature T i (ϕ) and DPCF i (ϕ) are presented.In the bottom plot (c) squares of the ion sound speed c 2 s (ϕ) and of the ion directional velocity u 2 i (ϕ) are presented.All the results are shown for T n = 7. Red dotted vertical lines indicate the locations of PE and SE.

FIG. 9 .
FIG. 8. Potential profiles obtained from simulations for several values of ε, i.e., densities with T n = T e = 1.Inserted graph represents the ion VDFs obtained for two of those densities -the largest and the smallest.The potentials at which VDFs have been monitored are indicated with open circles and numbers.

FIG. 10 .
FIG. 10.(a) Ion velocity distribution function at several points of discharge in collisional regime, dominated by Coulomb collisions (CC) in comparison with VDFs obtained in collision-free scenarios with other conditions kept unchanged, with the reference ion VDF (ε = 0) shown as well.(b) Ion temperature in the flow direction and profiles of the electrostatic pressure-related terms corresponding to discharge scenarios referred to in (a).
122, and m i the ion mass.E = −dΦ/dx ↔ dϕ/dx is the electric field and the ion source is modeled accordingly to Bissel and Johnson assumptions11in the form where the source strength S i = Rn n n e0 e βΦ/kT e e −m i 2 /2kT n /(2πkT n ) 1/2 , is modeled according to Harrison and Thompson26as a function of potential, i.e., ∼ n β e = e βΦ/kT e rather than of position.Factor β can take arbitrary values but here we employ only values β = 1 (after Bissel and Johnson11) and β = 0 (after Scheuer and Emmert12).The term Rn n 11 can be regarded either as the frequency of volume ionization ) DPFC drops to a value close to unity, otherwise it is considerable higher up to i,PE (T n = 1, ε = 0) 3.34 (see e.g., TableIin Ref. 20).
for more details.