Kinetic plasma-sheath self-organization

The interaction between a plasma and a solid surface is studied in a (1D–1V) kinetic framework using a localized particle and convective energy source. Matching the quasineutral plasma region and sheath horizon is addressed in the fluid framework with a zero heat flux closure. It highlights non-polytropic nature of the physics of parallel transport. Shortfalls of this approach compared to a reference kinetic simulation highlight the importance of the heat flux as the measure of kinetic effects. Non-collisional closure and higher moment closure are used to determine the sound velocity. Within these frameworks, no gain in the fluid predictive capability is obtained. The kinetic constraint at the sheath horizon is discussed and modified to account for conditions that are actually met in simulations, namely quasineutrality with a small but finite charge density. Analyzing the distribution functions shows that collisional transfer is mandatory to achieve steady-state self-organization on the open field lines.


Introduction
The prompt plasma recombination when in contact with a cold, dense and electrically neutral media drives a plasma flow toward this medium, which then appears as void for the plasma.Such a flow characterizes the plasma-wall interaction.It is akin to that observed in time-dependent boundary conditions for hyperbolic systems in which discontinuities and shock waves develop [1].A vast literature addresses these two classes of problems.For plasma physics, it goes back to the pioneering work of Langmuir and Tonks describing the properties of a plasma in contact with a solid surface [2,3].Unlike the neutral fluids where the shock waves are not resolved in the fluid framework, the shock waves in plasmassuch as the plasma/boundary interface-can be resolved when taking into account a departure from quasineutrality across a narrow region [4], the sheath region at the plasma/boundary.The width of the latter region depends on the characteristic scale of the Poisson equation relating the divergence of the electric field to the charge density.The key issue is then asymptotic matching between disparate scales, that of the quasineutral plasma 2L ∥ with L ∥ ≳ 10 m, and the Debye scale λ D ≲ 10 −4 m.The small parameter ε D = λ D /L ∥ ≲ 10 −5 is the control parameter of the self-organization process of the plasma interacting with a wall, a solid in most cases although a liquid or a dense gas could be considered.
The physics of plasma-wall interaction in fusion devices is a growing concern for burning plasma operation [5].The region where this physics take place has been named scrapeoff layer (SOL).Field lines in the SOL are connected to the wall at both ends.Three major issues are being investigated.First, parallel transport along the field lines and onto the wall is large so that particle and energy deposition patterns are narrow and become challenging problems for the present technology of heat and particle extraction [5,6].Second, wall erosion and impurity transport from the wall component into the confined plasma play a major role in the aging of the plasma facing components and on plasma performance.Third, the variation of the electric field in the SOL modifies turbulence instabilities and very likely plays a role in the onset and stability of the edge transport barrier [7].Control of such issues mostly depends on the plasma pressure, electron thermal energy and electric potential.These fields are shown to result from the balance between the sources that sustain the plasma on the open field lines and the losses governed by the sheath properties.
Our aim here is to analyze the self-organization along a field line between the quasineutral plasma, that we refer to as the SOL plasma, and the sheath boundary layer where charge separation occurs.The idea addressed in the companion paper is the means to incorporate the key results of this physics without having to resolve the Debye scale.This is especially important for the simulation effort of turbulent transport in magnetic confinement devices.To optimize the numerical codes, this physics is addressed for quasineutral plasmas so that the property of near constant fields in the parallel direction allows contemplating scales of order 100 or more Debye lengths when meshing in the parallel direction.The sheath physics cannot be solved and boundary conditions incorporating the sheath effects must be introduced, both in fluid codes stemming from the three coupled Navier-Stokes conservation equations, and more recently gyrokinetic codes.The specific issues in the gyrokinetic framework are discussed in the companion paper using kinetic simulations of the SOL and sheath plasma regions as a test bed for the physics.However, when analyzing the latter kinetic simulations, novel theoretical considerations have been developed.The scope of this paper is to present this part of the work, using one of these kinetic simulations as a reference to illustrate the key findings.
Addressing the full 6D kinetic problem is by far too demanding but the problem can conveniently be simplified to 1D-1V assuming that the chosen direction is parallel to the magnetic field line.The particle, momentum and energy sources on such a field line are mostly governed by crossfield transport due to turbulence and collisions.In the parallel direction, classical transport, in first approximation of free streaming particles, governs prompt losses onto the wall components.The reference problem is then reduced to a source region for particles and energy, possibly of momentum, which extends over a large part of the parallel domain (typically L ∥ ), connected by parallel transport to the narrow sheath region, which in turn controls the deposition onto the wall.The difference in parallel transport between electrons and ions generates an electric field.In steady state, quasineutrality is sustained over most of the plasma, in the SOL or presheath regions, and charge separation is localized in the sheath.To capture this self-organization process a two species plasma of electrons and singly charged ions is a minimum plasma model.An alternative to the kinetic description of the electron-ion plasma consists in using a fluid projection, which then depends on the number of moments that are selected and on the choice made to close the system.
The most familiar analysis of the SOL/sheath selforganization is addressed using the Mach number, ratio of the ion mean velocity and sound velocity.Given that the particles must flow outward at both ends of the field line, a stagnation point with zero Mach number must be found.When symmetry is enforced, the stagnation point is localized in the middle of the plasma domain.Because of the source terms, plasma acceleration drives an increase of the Mach number, which remains subsonic in the SOL region.The transition to a supersonic flow defines the sheath entrance.This is the well known Bohm condition at the sheath boundary [8].In a fluid approach one can then use this boundary condition at the sheath entrance [9].In the standard fluid framework based on particle, momentum, and energy conservation equations, the Navier-Stokes equations, the boundary conditions of the even moments, particle and momentum fluxes, are set by the source terms.Conversely, the Mach = 1 condition defined at the sheath location is the boundary constraint for the odd moment, the momentum flux.When considering higher fluid moments, which capture the departure from Maxwellian distributions, qualitative aspects of the sheath physics are recovered [10].However, reference to the kinetic solution remains crucial to validate such results.
Besides the fluid approach, the distortion of the distribution functions away from Maxwellians leads one to address this problem in the kinetic framework where the Mach number has little meaning.A kinetic criterion has been derived by Harrison and Thompson [11], then recalled in several papers [12][13][14][15].The importance of collisional processes has also been been discussed [16,17].Some kinetic details of the sheath physics have been addressed: regimes with multiple ions [18], the effects of E × B drift [18,19] or of grazing angle magnetic fields [20][21][22][23], wall parallel to the magnetic field [24].An extensive review of plasma models in the vicinity of the sheath, along with numerical methods used to study plasmawall interaction, can be found in [25].
We address here the self-organization of the SOL/sheath/wall system in a kinetic framework.The companion paper [26] aims at developing immersed boundary conditions [27][28][29][30] that would be appropriate in the gyrokinetic framework [31,32].The Voice code [33] has been used to perform the simulations, see companion paper [26].Analyzing the simulation data has shown that open issues arise whenever one addresses a problem with ε D small but finite.This is particularly true when determining the sheath entrance.While the Bohm criterion appears to be a robust condition since it is consistent with the picture of a standing shock wave, the kinetic simulations indicate that the relevant sound wave is not readily determined.The Bohm criterion is then model dependent with consequently a reduced predictive capability.In the course of the present paper, this issue appears to be coupled to the number of degrees of freedom used in the model, ranging from a small number in usual fluid descriptions to infinite in the kinetic descriptions.Furthermore, using the kinetic formulation published by Harrison and Thompson [11] raised new issues.Applying this formulation to the simulation output yields spurious sign changes of a quantity that is defined positive according to the assumptions made in reference [11].
We thus revisit these various approaches of the selforganized SOL/sheath/wall physics.We first address the SOLwall interplay model in section 2: presenting some results from a reference kinetic simulation that motivate this effort (section 2.1), the underlying kinetic equations (section 2.2), highlighting the source terms (section 2.2.2), and collision operators (section 2.2.3).We then step to analyzing the SOL/sheath self-organization in the fluid framework section 3. To facilitate the reading of the paper, and follow the main line dedicated to the self-organization properties, important contributions to the analysis within the fluid framework are presented in appendix, in particular the fluid conservation equation that stem from the kinetic two species equations (appendix A.1).The important issue of determining the sound wave velocity (appendix A.2) includes the discussion of the polytropic closure (appendix A.2.4) and that of higher moment closures (appendix A.2.5).The fluid predictions and their shortfalls in providing an accurate description of the steady state plasma self-organization, as determined by the kinetic simulations, are addressed in sections 3.1 and 3.2.We discuss means of completing the fluid model in appendix A.3 considering non-collisional closures (appendix A.3.1), matching conditions at the sheath horizon (appendix A.4) and the discontinuity of the steady state quasineutral fluid solution (appendix A.5).The findings presented in appendix as well as the comparison between the fluid predictions and kinetic simulations are summarized in section 3.3.Kinetic issues are discussed in section 4. We revisit the previously published kinetic constraint in section 4.2 and show that collisions are mandatory in the self-organization process (section 4.3).The steady state description of the distribution functions is presented in section 4.4 and the role the source and collisional terms investigated.A summary, discussion and conclusion section closes the paper, section 5.

First look at open field line self-organization
In this paper and its companion paper [26] we aim at determining appropriate boundary conditions when addressing turbulent plasma transport in magnetically confined devices.For a strongly magnetized plasma, a large asymmetry is enforced between the fast transport in the parallel direction and the reduced transport transverse to the magnetic field.In the chosen geometry of the magnetic field, the parallel transport is free at first approximation while it is hindered in the transverse directions and restricted to the Larmor gyration.This confinement in the directions transverse to the magnetic field is reduced by collisional and turbulent cross-field transport, both scaling typically with the Larmor gyration radius.The control parameter that accounts for this asymmetry is the parameter ρ ⋆ ≪ 10 −2 , ratio of a reference ion Larmor gyration radius and the plasma minor radius [34].Because the electron Larmor radius is smaller than the ion Larmor radius, one expects the ion dynamics to govern the properties of transverse transport.The ion Larmor radius then characterizes the size of the turbulent eddies.The ion Larmor radius being larger than the Debye scale, quasineutral plasma conditions prevail when addressing turbulent transport.We address here the plasma region named SOL at the bulk plasma outer boundary where field lines are connected to a wall component.In this region, under steady state conditions, the transverse transport is balanced by parallel transport.For the latter, the electron mobility is larger than the ion mobility and determines the smallest time scale to be addressed.Furthermore, for a field line connected to a wall component, the Debye scale cannot be ignored [35].Based on these simple arguments one finds that the physics of parallel transport to a wall departs from the conditions that prevail when investigating plasma transverse transport.The means of recovering boundary conditions mostly governed by parallel transport to be used in a quasineutral plasma framework is discussed in the companion paper [26].We address here the physics of parallel transport to a wall and the self-organization in the parallel direction between a quasineutral region and the sheath region where the charge density drives a large electric field, which in turn governs the losses to the wall.
For this purpose we consider a 1D-1V kinetic model describing the parallel physics along a magnetic field line in the SOL of a magnetically confined plasma.This line is connected to a wall at both ends.Simulations of this model have been performed and a reference simulation is used in this paper as a guideline.The simulation model and equations are presented in the companion paper [26].In the present work, magnetic field curvature drifts are neglected because these are small compared to the turbulent transport taken into account by the source terms.The role of these drifts can be addressed and modeled by charge sources and sinks, which will modify the symmetry.As reference case, one can have a divertor geometry in mind.However, the model is more generic and can be used to address other configurations of interest.This field line is then connected to the divertor low and high field target plates where plasma promptly recombines.One can note that in the present model, one only requires prompt recombination so that the wall component can be a solid, liquid or dense gas.The key property for the present purpose is that plasma particles, momentum and energy are absorbed, as if destroyed, when in contact with the wall, meaning on a time and space scale that is smaller than any scale relevant to plasma modeling.This condition is akin to that of a gas opening onto vacuum.For convenience, we shall consider a fixed location for the wall boundary as for a solid.Achieving steady state with non vanishing plasma conditions then requires a source to balance the outflux to the wall at both field line ends.In the regime we have chosen to address, the main source term on the open field lines stems from the divergence of the cross-field turbulent fluxes.The field line is assumed to intercept the divertor with a normal incidence, the treatment of a magnetic presheath in the case of oblique incidence requiring one to address the Larmor gyration, therefore stepping to at least a 2D-2V model.Furthermore, in the latter geometries the wall can be intercepted during the gyration motion at a distance from the intersection point of the field line with the wall [36].This introduces a dependence on the particle gyro-angle which questions the key assumption made to address plasma turbulence in the gyrokinetic framework [37] or the drift expansion for the fluid framework [38].To avoid such issues we address here the rather standard simplified geometry illustrated on figure 1.Of course such issues are important and 6D kinetic simulations in simplified geometries of plasma-wall interaction will enable one to identify the gyrokinetic shortfalls and means to circumvent them when addressing realistic geometries of plasma-wall interaction.
Before stepping to the description of the kinetic model, let us recall the expected behavior for the plasma between the volumetric source and the sink boundary condition as depicted on figure 1.In steady state, the particle source must be compensated by the build-up of a parallel particle flux Γ, figure 2. Because of the symmetry of the simulation set-up, only the right hand side of the simulation region is shown, from the middle of the box at x = 0 to the wall on the right hand side.As will be recalled in the following, the growth of the particle flux in the subsonic regime, hence in the plasma region prior to the sheath, governs plasma acceleration [4,35].Accordingly, one expects that when the source term vanishes in such a quasineutral region, the particle flux becomes constant and plasma acceleration stops [4,35].Without particular conditions, this flow with mean velocity u i remains subsonic up to the wall where a standing shock wave appears to develop.The latter is in fact characterized by a sharp increase of the gradients so that very small scales are generated, typically at the Debye scale λ D .In plasmas, the shock wave discontinuity can be resolved by addressing a departure from quasineutrality, typically on the Debye scale [4].The reference simulation allows one to illustrate this point.We define the characteristic scales L −1 Here n i stands for the ion density, hence the subscript i, figure 3   Outside the source region, L u i levels-off, then, left of the vertical dash-dot line, L u i exhibits a strong decrease toward values of the order of 10 Debye lengths.The length scale L n i for the ion density exhibits the same behavior as L u i between the end of the source region and the wall.The behavior is different in the source region, in particular close to x = 0 where for symmetry reasons L u i → 0 while L n i → +∞.In a narrow region close to the wall, on the right hand side of the vertical dashdot line on figure 3, one observes the expected sharp increase of the gradients.This effect is more pronounced for the density than for the mean velocity.Indeed L n i is reduced to the order of 1 Debye length.However, since the Debye scale is resolved, no discontinuity is observed.
If one now considers the mean velocity of the plasma u i , figure 4, one can observe first a rapid increase in the source region, followed by a slower, close to monotonic increase typically from the source region to the wall, steepening in the wall vicinity, as expected from the behavior of L u i figure 3. On figure 4, the plot of u i is the blue curve open circles, and is compared to that of the thermal velocity V, black curve headdown open triangles, and an expression of the sound velocity c s obtained in the fluid framework, black curve head-up open triangles.The thermal velocity V depends on the thermal energy T = T e + T i , sum of the electron T e and ion T i thermal velocities: m i V 2 = T where m i is the ion mass.The sound velocity is defined by c s = √ 3 V.One can readily notice that u i exceeds the thermal velocity V before the end of the source region, then further increases but does not appear to reach the sound velocity.With respect to the standard plasma description on an open field line [4,35], as sketched above, the reference simulation exhibits unexpected features.First, both fields n i and u i vary in the quasineutral region with vanishing source, hence constant particle flux Γ, figure 2. We stress here that these are steady state features of the simulation.Second, the mean velocity does not appear to reach the sound velocity at the sheath entrance, indicated by the vertical dash-dot line.The latter observation suggests that either the Bohm criterion does not hold, or the effective sound velocity is not properly determined by the usual fluid framework.Since the Bohm criterion is based on physics arguments, it is tempting to state that the Bohm criterion is fulfilled but that the sound velocity is not accurately determined.This problematic discrepancy for plasma turbulence modeling in the fluid framework is further addressed in the following.Alternatively one could question the location of the sheath entrance.Determining the sheath entrance is also addressed in the following and various criteria are analyzed and compared using the reference kinetic simulation.

Kinetic model to address open field line self-organization
2.2.1.Kinetic equation and normalization.In section 2.1 we have presented some of the results obtained with the reference kinetic simulation.We now describe the model used in this paper and for the simulations.We address the kinetic behavior of a two species plasma.For simplicity we shall consider singly charged ions.The effect of cross-field transport, mostly turbulent transport, is taken into account by a source term.With this simplification, one can restrict the model to a single field line with one-dimension in position space s and one dimension in velocity space v.One now introduces the distribution function f a (s, v, t), expressing the density of particles of species a, e for electrons and i for ions, at time t and at point (s, v) of phase space.Mass and charge of any particle of species a are m a and e a .The ratio m e /m i is a free parameter, usually chosen small.Varying this parameter gives access to the role of mobility on the physics at hand.The differential mobility between positive and negative charge plays a strong role in the sheath physics, for a plasma the negative charge mobility is smaller than the positive charge mobility, while the inverse is met in dusty plasmas.This parameter also controls the simulation cost since steady state conditions for the ions is typically m i /m e longer than for electrons.For the reference simulation we have chosen m e /m i = 1/400, which is small and yields a factor 3 gain on the computation cost compared to that required for a realistic mass ratio of a deuterium plasma.In the companion paper [26] simulations scanning the mass ratio and including realistic mass ratios, are presented and analyzed.The electric field E in the electrostatic limit is determined by the gradient of the electric potential U, E = −∂ s U.The 1D-1V Boltzmann equation (1a) that governs the evolution of the distribution function f a is then: where C(f a ) is the collision operator, standing for both selfcollisions with particles of the same species a and inter-species collisions, and where S(f a ) is the source term.The electric potential is determined by the Poisson equation (1b), the electrostatic limit of the Maxwell-Gauss equation, which relates the electric potential to the charge density Here we have introduced the local density n a of species a, which is a function of time and space.It is the velocity space integral of the distribution function n a = ´dv f a .In the present paper we only normalize the Poisson equation introducing the reference thermal energy T 0 , density n 0 , and length scale L 0 .We define ϕ = eU/T 0 and x = s/L 0 , so that: We have introduced the Debye length λ D0 to define ε 2 D , the control parameter in the Poisson equation.In this paper the source term S is the SOL-plasma source mainly governed by cross-field transport.Specific sink terms are addressed in the companion paper to implement the immersed boundary conditions.All simulation data are normalized.A normalization is used on the various plots presenting the simulation output, and recalled on each plot for clarity.On the plots, densities and thermal energies are normalized by n 0 and T 0 respectively, and all velocities are normalized by V 0 = T 0 /m i .Consequently, the particle flux normalization on figure 2 is Γ 0 = n 0 V 0 .Note that λ D , the scale normalization in the plots, stands for λ D0 .This normalization choice on the various plots differs from the normalization made in the companion paper.Let us recall that the only normalized variables in the equations of the present paper are the electric potential ϕ and position x, see equation (2).

The forcing problem: the plasma source.
Because of the plasma losses onto the wall, the source term is a crucial aspect of the self-organization process along the SOL field line.In the SOL region of magnetically confined plasmas, cross-field turbulent transport is found to govern the heat source.The latter, and, depending on the operating conditions, neutral particle ionization on the open field lines determine the particle and convected energy sources.The chosen source term S in equation (2a) is similar to the Gysela source terms [39].While in Gysela this term stands for actual particle, momentum and energy source governed by particle ionization, torque injection and plasma heating, the source term in the present 1D model is understood to be the divergence of the cross-field fluxes from neighboring field lines.Considering that turbulent transport is observed to be ballooned to the low field side [40,41], it is reasonable to address a localized source term as done here, which then differs from the homogeneous source addressed in some works [4].One could extend the model and consider a particle source by electron impact ionization of neutrals, or a spurious external heating of SOL electrons as that investigated for lower hybrid launchers [42].One could also address problems with charge or momentum sources, which would induce parallel electric currents and break the symmetry of the plasma region.We focus here on conditions that enforce the symmetry of the plasma and that minimize the role of atomic processes as suitable to perform global gyrokinetic turbulent transport simulations with somewhat reduced complexity.Neglecting the SOL ionization source in the present work is consistent with the SOL behavior in standard limiter configurations and in the so-called sheath limited divertor regime.Some features of the transition into the high recycling divertor regime could also be addressed insofar that atomic processes can be ignored to explain the plasma properties.Estimating the particle source by impact ionization, coupling the gyrokinetic plasma to kinetic neutrals, or accounting for impurity line radiation in a kinetic framework can be developed.It will require appropriate testing of the physics and of the computation cost, typically in simplified setting such as that provided by Voice [33].The SOL problem of interest is therefore that of the hot plasma regime where atomic processes are expected to have a relatively weak effect [43].This particular regime has implications when addressing the plasma collisionality, see section 2.2.3.The particular selforganization problem is to determine the plasma total pressure and the plasma thermal energies on the open field lines, given the wall losses monitored by the parallel transport, the sheath constraints on the fluxes as well as the source properties.A particular choice of the source term is made in this paper and the companion paper.While one can expect E × B drift convection to act as a source of particles and heat from neighboring field lines and a sink of particles and heat from the chosen field line, suggesting a Bhatnagar-Gross-Krook (BGK-like) source term, we have chosen a particle and convected energy source independent of the plasma condition on the field line.The source in the kinetic equation (2a) is proportional to a target distribution function, chosen to be a Maxwellian, with zero mean velocity, and therefore no momentum transfer.This source term determines a particle and convected energy source with zero heat source-energy source with zero particle source, therefore exchanging a cold by a hot particle The mask M s , plotted in black on figure 2, right hand side scale, localizes the source region in the central part of the plasma.The shape is based on hyperbolic tangent functions.The width of the transition region is comparable to the width of the source region so that there is no flat top in the source amplitude.This mask function is defined so that its integral in the plasma volume is equal to 1.The parameter s k then defines the source amplitude as the number of particles of species a injected in the plasma per unit time.The source thermal energy T s determines the energy of the injected particles, and therefore the convected energy source.These two control parameters are chosen to be species independent and constant both in space and time.This source does not depend on the SOLplasma conditions and is symmetric with respect to the center of the simulation domain at x = 0. Together with the absence of momentum and charge source, this enforces the symmetry of the solution with respect to x = 0 that defines consequently the stagnation point where the plasma mean velocity is null.
A notable simplification is made here compared to intermittent transport that is understood to prevail in the SOL region [44].With intermittent cross-field transport events, avalanches or blobs, the source term exhibits a time dependent pattern with a short burst of plasma convected into the SOL, with hot and dense plasma conditions, followed by a quiescent phase prior to a new burst [45].Both the duration of the bursts and that of the quiescent phase play a role on transport properties and therefore on the source properties.One can then use the statistics or time traces of the divergence of turbulent fluxes to determine the time dependence of the source term in Voice simulations [26,33].The interesting issues raised by such dynamics of the source term and interplay with the parallel transport are left for future work.The constant source approximation allows one to address the steady state properties without having to perform statistical averages on the distribution of burst properties.

Collisions.
In the simulations, the collision operator C(f a ) in the simulations of the Boltzmann equation (2a) is the sum of the linearized self and inter species collision operators.A detailed description of the chosen linearized operators is given in the companion paper [26].
Collisions are important as regularizing operators smoothing the velocity space variations of the distribution functions.They are also observed to contribute to the transport process and are found to be crucial to achieve steady state conditions, see section 4.3.In the self-organization physics on open field lines, one recovers the fundamental difference between a collisionless plasma with ν ⋆ 0 = 0 and a weakly collisional plasma, ν ⋆ 0 → 0 + .The dimensionless control parameter ν ⋆ 0 is defined in equation ( 4), similar to the standard ν ⋆ parameter [46] which is used in Gysela [47].However, it is slightly modified to account for the fact that we do not emphasize trapped particles in the SOL collisional transport, see also [5].To address the so-called sheath limited regime in gyrokinetic simulations of turbulent transport, hence a low-collisionality SOL regime without complex atomic collisions, simplified collision operators can be considered.In the Voice simulations [26,33], the collision operators stand for self-collisions as in [41] and inter-species collisions adapted from [48].They conserve density, momentum and energy.Maxwellian distribution functions belong to the kernel of the self-collision operator.The latter also exhibits the 1/|v| velocity dependence of the Landau and Lenard-Balescu non-linear operators [49].The inter-species collisions are set to drive the distribution functions toward thermal equipartition and equal mean velocity.For the present purpose, we assume that the details of the collisional process are not crucial so that the chosen features of the collision operators are sufficient.The amplitude of the collisional term is determined by the chosen collisionality ν ⋆ 0 , such that where log Λ is the Coulomb logarithm.The control parameter ν ⋆ 0 does not depend on particle mass.It is more convenient for SOL physics than the standard ν * parameter since the parallel transport process does not depend on trapped particles.The control parameter ν ⋆ 0 can be read as the inverse of the collisional mean free path divided by the reference field line length in the SOL: L ∥ ≈ π qR, R stands for the major radius and q is a relevant value of the safety factor, typically the value taken one energy e-folding length within the separatrix.

Analysis of SOL-sheath self-organization within the fluid framework
The fluid framework is most often understood to be restricted to the Navier-Stokes equations, namely the three first moments of the Boltzmann equation.A closure is introduced to bound the infinite fluid hierarchy to these three first moments.In the literature of fusion plasma turbulence this closure is usually a collisional closure that determines the heat fluxes in terms of the temperature gradients.The Navier-Stokes equations are also interesting when analyzing the kinetic results since these provide physical insight into the ongoing processes.However, given the mismatch between the expected behavior and simulation results presented in section 2.1, one can expect a shortfall of the standard Navier-Stokes closure.In this section we compare the fluid predictions to the kinetic reference simulation.Non-collisional closure and higher moment closure are addressed.It turns out that neither closure is able but to provide a predictive guideline to analyze the simulation evidence.

Fluid properties in steady state plasma conditions
In steady state conditions discussed in appendix A.1.4,one can combine the electron and ion conservation equations to obtain plasma balance equations for the number of particles, plasma momentum and plasma energy, equation (A6).The derivative d s is the derivative along the curvilinear abscissa along the field line that stems from the divergence of the fluxes, particle flux Γ, total plasma momentum flux Π and energy flux Q.When enforcing quasineutrality, hence a vanishing charge density ρ c in equation (A6b), these balance equations are a closed system.Since the chosen source term enforces a convective energy source, proportional to the particle source S n , with T s /2 characteristic energy of the injected particles, equation (A6c) is simplified so that the steady state plasma balance equations take the form: As discussed in the appendix, the kinetic simulation is found to exhibit the three conserved fluxes equation ( 5) while the distribution functions vary with no direct constraint enforcing such conservation.Since no heat is injected by transverse transport, a reasonable assumption to close this system is to set the parallel heat flux to zero q = 0, hence Q = Q c (to simplify the expressions we also set γ c used in appendix to 2γ c = 3).In the quasineutral limit, we also replace Π by Π qn since Π qn = n i T + n i m i u 2 = Π + ηT e , where T = T i + T e and η = ρ c /e.With these simplifications, and taking into account that with the chosen source both Γ and Q are anti-symmetric with respect to the mid-box position s = 0 while Π is symmetric, one obtains: Given the velocity V defined by According to equation ( 6), Γ increases from 0 at x = 0 to a maximum value Γ w in the source region and is then constant up to the wall.When η = 0 in the quasineutral regime, Π ≈ Π qn is constant.Finally, given equation (6c), one can determine T in terms of Given the expression of Π qn , one finds Since the polytropic index γ p is defined by dp/p = γ p dn i /n i , one then finds: This index is not constant and varies monotonically with M 2 V .In particular, it ranges in the plasma region from a minimum value of 3/2 for M 2 V = 0, to the value of the adiabatic index γ p = 3 for M 2 V = 3.One finds therefore that the present model, with the chosen fluid closure q = 0, does not support the polytropic closure since γ p is not a constant.However, both the calculation of the sound velocity based on three fluid moments equation (A33c) and the analysis in terms of the polytropic index yield identical expressions for the sound velocity.
Let us now replace the density in the plasma pressure [4].However, T/V = m i V also depends on M V as found in equation (7a).It is then more relevant for this nonisothermal case to consider a further change of variable One then finds that M increases from M = 0 at the stagnation point with Γ = 0 and therefore decreases and at M → +∞, A 2 s → 8/9.In equation (7c), we have replaced M V by the actual Mach number M = M V / √ 3. The numerical coefficients that appear in the definition of A s are chosen to enforce that A 2 s = 1 is the maximum value of A s .One can also use equation ( 7c) to determine M 2 as a solution of a second order equation.When varying A s , one finds a first regime with a single positive solution for 0 ⩽ A 2 s ⩽ 8/9.For A 2 s = 8/9, the system is singular since equation (7c) reduces to a first order equation in M 2 with positive solution M 2 = 1/3.When A 2 s > 8/9, two positive solutions of the quadratic equation are obtained for M 2 .At constant momentum flux Π qn and positive particle source S n ⩾ 0, the parameter A s in the source region increases monotonically from 0 at s = 0.It leads to M 2 ⩽ 1 in the quasineutral plasma.The matching constraint presented in appendix A.4 is similar to that used in [4] and yields M ⩾ 1 in the sheath region, with therefore M 2 = 1 at the sheath entrance.The self-organization process therefore enforces that the momentum flux is adjusted so that The plasma total momentum flux in the SOL Π is identical to Π qn when quasineutrality holds.It appears to be determined by the upstream source, as highlighted by the dependence on Γ w and T s , and to the sheath constraints via the matching condition.
With such steady state conditions, one can then write A s = Γ/Γ w so that A s is fully determined by the particle source, which in turn yields M given equation (7c), and T = (T s /3) 2/(1 + M 2 ), equation (7a).At constant momentum flux Π qn , set by equation (7d), the density n i decreases from the value at the stagnation point With further assumptions, one can step to determining T i and T e .
A first regime is found at high collisionality, such that the equipartition transfer is the leading contribution in the energy balance equations.In this asymptotic regime one then finds T e = T i = T/2 and therefore from equation (7a): In a second regime characterized by weak energy coupling between ions and electrons, one splits the energy balance equation into the electron and ion energy balance equations assumed to be independent.One thus neglects the collisional equipartition as well as the Joule energy transfer eΓ|d s U| from the electrons to the ions.With these assumptions, one obtains: This yields: The possible change of sign of T i can be read as the signature that the Joule heating of the ions by the electrons cannot be ignored.This energy transfer corresponds to the ion acceleration by the expanding electrons.One can remark that with these assumptions T e ⩾ T i so that the collisional energy transfer would also occur from the electrons to the ions.This result can therefore be understood as yielding an upper bound for T e and a lower bound for T i .
In the previous calculations, the sheath constraints are not taken into account.At the wall one should recover Q = γΓT sh e where γ is the sheath energy transmission factor, γ > 3 and T sh e is the electron thermal energy at the sheath entrance.A third calculation of the thermal energies T e and T i is achieved by enforcing the sheath conditions.Balancing the energy exhaust Q = γΓT sh e with the energy source ΓT s then yields the electron thermal energy T e ≈ T sh e when assuming the electron thermal energy to be constant in the parallel direction Then setting γi = γQ ci /Q, one obtains: The ratio Q ci /Q must still be estimated.Neglecting the heat fluxes, one can assume the same energy on the electron and ion channel so that Q ci /Q ≈ 1 2 , therefore: This constraint on the energy exhaust can be seen as giving a lower bound for T e and consequently an upper bound for T i .One can note that in such a regime one finds that T i > T e for M 2 small enough.Collisional energy exchange then transfers energy from the ions to the electrons and is opposite to the Joule energy transfer.

Fluid interpretation of the kinetic simulation
The control parameters of the reference simulation used as evidence are listed in table 1.The mass ratio m i /m e is a free parameter in the simulations and is scanned in the companion paper [26] including realistic values.The size of the simulation domain L ∥ /λ D0 is also a free parameter.For the chosen reference simulation, both mass ratio and domain size are chosen to reduce the simulation cost in the appropriate asymptotic regime m i ≫ m e , L ∥ ≫ λ D0 .The collisionality ν ⋆ 0 is chosen to be consistent with the effect of collisions in the sheath limited divertor regime.The source, proportional to a Maxwellian with zero mean velocity, is characterized by two control parameters, the injected particle flux Γ w and the injected energy flux Q w = Γ w T s , where T s is the source thermal energy.The latter defines T 0 .Given T 0 , L ∥ /λ D0 and ν ⋆ 0 are two functions of the remaining unknowns, L ∥ and n 0 .The normalization parameters n 0 and T 0 , and consequently λ D0 , are therefore defined by the control parameters.
The simulation is run until steady state conditions are achieved.The various moments of the distribution that define the fluid quantities are then computed by integrating the electron and ion distribution functions.The electric potential is an output of the simulation.
Before analyzing the simulation evidence, we first introduce key positions, the position characterizing the end of the source region x sce , the sheath entrance x sh and the wall location x w .The sheath entrance is determined in section 4.2, table 3criterion dρ c /dϕ.The wall location, with penalized wall conditions, is addressed in the companion paper [26].The three reference positions are presented in table 2. In most figures, the end of the source region is indicated by a vertical dashed line and the sheath entrance by a thin vertical dash-dot line, the wall position is the maximum value on the x-axis.
We first address the variation of u i ≈ u e as on figure 4, but comparing u i to c eff = 1.505V for x ranging from the stagnation point at x = 0 to the sheath entrance.This sound velocity c eff is determined from the simulation data to intersect u i at the sheath entrance as shown on figure 6.The latter is a smooth transition out of quasineutrality conditions.The choice of c eff = 1.505V would have to be slightly adapted to changes in the sheath entrance definition.The variation of A s , equation (7c), is also plotted, right hand side axis.The control parameter is computed using the simulation output for Γ and Π qn .One can observe that A s becomes constant outside of the source region, to the left of the vertical dashed black line.Furthermore, A s exceeds unity for x ⩾ 38, indicated by the vertical black dash-dot line.Since the value of Γ is close to that predicted, A s is overestimated by ≈12% because Π qn is observed to be ≈12% smaller in the kinetic output than predicted by the fluid model equation (7d).
In the region where A s levels off at the end of the source region, one finds that V < u i < √ 3 V, however u i is still increasing while V is decreasing so that the ratio u i /V, proportional to the Mach number, is increasing with A s constant.In the kinetic simulation, the Mach number cannot only depend on A s .In the fluid framework, the number of degrees of freedom is fixed by the number of chosen moments.We find here that with three degrees of freedom, the Navier-Stokes system does not capture all the feature of the kinetic physics.
Let us now address the energy flux.The energy source Q s = ΓT s /2 is identical for both species and determined by the particle source and the thermal energy of the injected particles T s .The latter is used to define T 0 normalizing energies of the simulation output, see table 1.Within the source region, one finds that the electron and ion energy flux, respectively Q e and Q i , increase with the energy source Q s , figure 7. Q e is the blue curve with closed circles, Q i the black curve with open circles, and Q s the black dashed curve.Toward the end of the source region, Q e tends to level off Q e < Q s , before decreasing.Conversely, Q i steadily increases, Q i > Q s .One finds that energy is transferred from the electrons to the ions at conserved total energy, Q = 2Q s = Q e + Q i .At the sheath entrance, the ion channel has increased to 60% of the total energy flux Q, and the electron channel is reduced to 40% of Q.One finds therefore that about 20% of Q s , the energy coupled to the electron channel, has been transferred to the ion channel.As will be shown in the following T i ⩾ T e .The energy transfer is governed by the term eΓE from the electrons to the ions, which must also balance the collisional equipartition transfer from the ions to the electrons.One can split the total energy flux Q = Q c + q into convective flux Q c up to 75% of Q and heat-flux q, typically 25% of Q.Both Q c and q mostly build-up in the source region, to the left of the vertical dashed line.Toward the end of the source region and up to the wall, there is a slight energy transfer from q to Q c , figure 8. Black curve head-up open triangles for Q, blue curve open triangles for Q c and black curve head-down closed triangles for q.One can notice that the heat-flux is sustained at a level of about 20% of Q at the wall, where Q must be absorbed in steady state conditions, the plasma energy sink at the wall surface must compensate the plasma volumetric energy source.This implies that the heat-flux q as well as the convective energy flux must contribute to the energy flux leaving the plasma via the sheath and deposited onto the wall.
For the ions, figure 9, one finds that most of the energy flux Q i is convective Q ci .The ion heat flux q i levels off at less than 10% of Q s and slowly decreases toward the wall.Black curve head-up open triangles for Q i , blue curve open squares for Q ci , black curve head-down closed triangles for q i , and dash-dot black curve for Q s .
For the electrons, figure 10, the energy flux Q e , the convected energy flux Q ce and the heat-flux q e increase in the source region and then gradually decreases.Black curve open headup triangle Q e , blue curve open squares Q ce , black curve closed head-down triangles q e and Q s dash-dot curve.One also finds that the convective flux account for slightly more than 50% of the electron energy flux.It is therefore of the same magnitude as the heat flux q e .Both components of the electron energy flux tend to decrease downstream from the source.Toward the end of the source region, when the source tends to zero, one can observe that the convective energy flux Q ce is roughly constant, while the total energy flux Q e decreases with the heatflux q e .The combined effect of collisions and energy transfer via the electric field eΓE govern a rather complex reorganization process of the energy transport.As for the total heat-flux q, both the ion and electron heat flux do not vanish toward the wall and contribute to the energy exhaust.This point will be further discussed when addressing non-collisional closure in appendix A. 3.With such a closure, one finds that the heatflux is a function of the particle flux-it exhibits a non-zero projection on the particle flux.This result does not agree with the understanding of the heat-flux as being an energy transfer with zero net particle transfer.The latter picture of the heatflux suggests an exchange of hot against cold particles with equal particle flux in opposite directions.At the wall contact, a non-vanishing heat-flux then requires a particle flux moving out of the wall, contradicting the complete plasma absorption by the wall.One is then led to argue that the wording heat-flux for the flux proportional to the fluid moment N 3 might be inappropriate when kinetic features cannot be ignored.Coupling to higher moments must then be addressed and consequently the understanding of the moment N 3 should be adapted.
The electron thermal energy T e determined by the reference simulation is characterized by a near constant value.It very slightly increases in the source region with T e /T s ≈ 0.163.The variation of T e is small x = 0 to x ≈ 100, then slowly decreasing toward T e /T s ≈ 0.148 at the sheath entrance, see figure 11.Blue curve closed circles, T e .The equipartition value T eq is given in equation ( 8) but taking into account that Q c /Q < 1 due to the non-vanishing contribution of the heat-flux, open triangles T eq /T s .Note that T eq the equipartition prediction of T e yields T e /T s ≈ 0.15 at the sheath entrance.The equipartition condition appears to be met close to the sheath entrance.This could be a coincidence.The fluid description using q e = q i = 0 as closure, leads to T e constant T e /T s = 1/3, therefore too large by roughly a factor 2. Reducing the effective source of convective energy flux to account for the energy transferred to the electron heat-flux q e yields a reasonable order of magnitude at x = 0, T e /T s ≈ 0.154, but increasing with x and overestimating T e /T s The ion thermal energy T i exhibits a variation in the parallel direction in qualitative agreement with the expected behavior when neglecting the heat flux, figure 12. Blue curve closed circles T i , black curve open triangles T eq , black curve open diamonds T iγ .One finds however that T i is smaller than the equipartition value equation ( 8) but exhibits a comparable profile.Correcting this value to account for the heat-flux yields T eq , which is found to underestimate T i and to give a similar profile although with reduced amplitude.One can expect that T eq gives a lower bound of T i because of the energy transfer from the electrons to the ions via the parallel electric field.Finally one determines T i according to equation (11b), T iγ = T s (γ − 3M 2 )/(3γ e (1 + M 2 )).The latter expression gives a fair agreement with T i overestimating T i close to x = 0 and then underestimating T i for x ≳ 100.One can notice than T eq is close to the values of T i and consequently T e at the sheath entrance.Equipartition is observed T i ≈ T e at this location.Otherwise, from x = 0 to the sheath entrance, one observes that T i > T e .Despite the same energy source for the ion and electron channel, the fact that the ion energy flux is larger than the electron energy flux and that the electron heat-flux is larger than the ion heat-flux q e > q i enforces that 3T i + m i u 2 i is sufficiently larger than 3T e to yield T i > T e despite the large value taken by u i in most of the plasma region.One can also determine the sheath transmission factor for each species γ i = Q i /(ΓT i ) and γ e = Q e /(ΓT e ), which yields γ i ≈ 4.1 and γ e = 2.7 both values being larger than 1.5.The large value for γ i reflects the importance of the energy proportional to u 2 i , of order 3(T i + T e ), in the ion energy flux.

Fluid prediction capability of SOL-sheath self-organization
Let us first recall that when referring to the fluid description, we mostly consider a model based on the three Navier-Stokes equations for each plasma species of interest.

The plasma fluid description.
We summarize here the results presented in appendix A.1.In the quasineutral limit and neglecting the terms depending on the electron inertia in the momentum and energy fluxes, one finds that three invariants constrain the steady state plasma variation along the field line.For an open system, with volumetric particle source, hence depending on the curvilinear abscissa s, and boundary particle sink, the known particle source determines the particle flux.For a symmetric source with respect to the middle of the box at s = 0, the particle flux is antisymmetric and null at the stagnation point s = 0.By definition, the particle flux is proportional to the plasma density n and ion mean velocity u.In the region with vanishing particle source, the particle flux is invariant.
The momentum balance is enforced by the fact that no external source is applied to the plasma, and that binary collisions ensure momentum conservation, both in the collision process and in the mean field momentum transfer via the internal electric field.The plasma momentum flux is therefore invariant.The simplification brought to the expression of the momentum flux in the quasineutral limit and when neglecting the electron inertia then allows one to express the momentum flux in terms of the plasma density n, ion mean velocity u, equal to the plasma mean density, and plasma thermal energy T.
Finally energy conservation, with volumetric energy source and energy sink at the boundary, yields the third invariant.For a convective energy source, the energy source being proportional to the particle source, the symmetry of the source region enforces that the energy flux is antisymmetric and null at the stagnation point s = 0. Neglecting the heat flux, the energy flux is then proportional to the plasma density n, plasma mean velocity u and plasma thermal energy T. In the region with vanishing particle source, and hence vanishing energy source, the energy flux is invariant.
When restricting the fluid model to the Navier-Stokes equation for the plasma, one finds that one can replace the three fluid unknowns n, u and T by the three fluxes Γ, Π and Q, which are invariant when the volumetric sources are null.This property can be used to analyze the accuracy of the kinetic simulations.
In the region with vanishing source, left hand side of the vertical dashed line on figure 2, Γ = Γ w is constant and equal to the particle flux impinging onto the wall.Consequently Q = Q w , the energy flux impinging onto the wall.In the reference simulation these three conservation properties are observed Since all normalized values are of order 1, the relative and absolute error are identical.One finds that the largest error is made on the particle flux, with a systematic error with respect to the expected value Γ/Γ 0 = 1.Since the particle flux is a moment of the distribution functions, determined without requiring any approximation, one can use this comparison to estimate the numerical error of the simulation, typically smaller than 10 −2 .The kinetic simulation is therefore found to exhibit three conserved fluxes of the fluid description while the distribution functions vary with no direct constraint enforcing such conservation.

Model dependent plasma sound velocity.
Fluid models and closures at higher moments are discussed in appendix A.2.5, non-collisional closure in appendix A.3.1.As discussed in appendix A.5, the sheath entrance in the vicinity of the wall can be understood as a discontinuity or bifurcation point of the quasineutral plasma variation along the field line.As such, it corresponds to sound waves emitted in the counter direction to the mean plasma flow and generating a shock wave when the mean plasma flow reaches the sound velocity, namely when the sound waves are emitted at zero phase velocity in the wall framework.This criterion is the well known Mach 1 Bohm criterion for the sheath entrance.Although this argument is quite robust, one finds that its practical use crucially depends on the value of the sound velocity, which is not readily determined as addressed in appendix sections A.2.3, A.2.4, A.2.5 and A.3.2.
The Bohm condition emphasizes the particular role of the sound velocity c s on the plasma behavior in the very vicinity of the wall.However, the sound velocity itself is not readily determined and is found to depend on the particular closure of the fluid system.The reference dimensionless parameter should be the Mach number u/c s , however, the uncertainty on the appropriate value of c s leads one to use M V = u/V, with the previously determined variable V = T/m i .In the isothermal limit one finds c s = V while for convective energy flux and convective energy source one obtains c 2 s = 3 V 2 : the sound velocity is model dependent.For higher moment closures, the sound velocity is found to depend on the closure that is chosen, since the latter is arbitrary, this does not give access to predictive capability.
The polytropic closure yields c 2 s = γ p V 2 .With no heat exchange γ p = σ p /σ n = 3 is the adiabatic index and one recovers the same sound velocity as obtained for convective energy transport.If one identifies this result to the simulation result, where c s ≈ 1.5 V is observed, one finds γ p = 2.265, which is 25% smaller than the adiabatic index.This result suggests that heat transport plays a role so that the temperature variation is less important than would be expected with the adiabatic closure.
The non-collisional closure also leads to a result emphasizing effects governed by the heat flux and yielding c s < √ 3 V.Although this appears to be the proper trend, there is no agreement with the behavior observed in the simulation.These model dependent results do not provide a prediction of the plasma mean velocity at the sheath entrance.

Qualitative predictions of fluid models.
The fluid prediction of the plasma behavior, compared to the reference kinetic simulation, appears to give a qualitative description of the plasma variation from the source region to the wall.Conserved fluxes are recovered.The fluid mean velocity and the thermal energies of each species determined with the fluid model exhibit the appropriate trends but the agreement is rather poor.Another shortfall of the fluid description is that outside the source region, where the plasma particle, momentum and energy fluxes are constant, one can observe variations of the thermal energies, mean fluid velocity, heat flux, etc.Such variations are not expected in the fluid framework.These can only be understood as kinetic effects, including collisions, that modify the distribution functions.In the kinetic response one finds therefore that the number of degrees of freedom of the system is larger than the three retained within the fluid Navier-Stokes framework.
One can also observe that the actual value of the momentum flux Π is smaller in the kinetic simulation than predicted.The self-organization properties are different.Given these differences, as stated above, there is no accurate prediction regarding the sound velocity.Consequently the Bohm criterion does not provide a predictive capability to determine the sheath entrance.However, knowing the location of the sheath entrance, and therefore the value of the fluid mean velocity at that location, one can modify the closure properties to recover consistency.
Another difference, is the contribution of heat-fluxes to the energy flux.The latter is modest for the ion channel but amounts to nearly 50% of the electron channel.The heat-fluxes are also found to contribute to the energy fluxes to the wall, a property that does not agree with the picture of the heatflux being an energy exchange at zero particle flux.(Since no plasma outflux from the wall is possible, with only neutral outflux after plasma perfect recombination, therefore zero plasma particle outflux, the zero net particle flux constraint for heat exchange would then enforce zero heat-flux to the wall.) Alternative closures, polytropic, non-collisional, and at higher fluid moments do not improve the predictive capability but remain a tool for interpreting the evidence from kinetic simulations.They provide in particular a means supporting the role of the heat-flux in the self-organization process between the SOL plasma and the sheath constraints.
Finally, the analysis of the transition into the sheath indicates that one expects a change in the characteristic scale of the charge density gradient d s η > 0 at the sheath entrance.Resolving the Debye scale then removes the shock wave like discontinuity at the sheath entrance, driving a smooth transition into the sheath and steepening of the gradients on the Debye scale toward the wall.

Kinetic processes in the SOL-sheath self-organization
We now analyze the kinetic constraint that determines the sheath entrance.It was first published in [11], then in later papers [13].We first present the original derivation and then adapt it to the actual conditions met in the simulations.

Free streaming plasma at the sheath horizon
The constraint is addressed in steady state at the sheath horizon assumed to be standing at x → +∞ from the wall.The starting point is the normalized Poisson equation (2b), slightly modified here to account for a possible change of the normalization length scale L This equation takes the form of a Newton equation with mass ε 2 D = λ 2 D0 /L 2 , position ϕ, time x and applied force −ρ c = n e − n i .The sheath entrance stands at the horizon x → +∞ and is labeled ∞.In a first step we compute the first integral, multiplying equation (12a) by dϕ /dx and integrating over x At the sheath horizon x → +∞, the potential ϕ ∞ is assumed to have a finite value, which enforces a vanishing derivative dϕ/dx.The function G(ϕ), the opposite of the potential energy of the dynamical system, is expanded in the neighborhood of ϕ ∞ , which yields: Since quasineutrality holds at the sheath horizon ρ ∞ → 0, one then finds that: The sheath horizon constraint is therefore that This constraint is not specific of either the fluid or kinetic framework and indeed one can readily show that this constraint is equivalent to the Bohm constraint on the Mach number M 2 ⩾ 1.The kinetic formulation of [11] is determined by computing dn i /dϕ and obtaining an expression depending on the ion distribution function f i .To achieve this calculation we consider the neighborhood of the sheath horizon and the energy conservation of the ions with potential energy ϕ, hence Here H is the total particle energy normalized by the thermal energy T 0 so that the velocity appearing in K is the phase space velocity normalized by the ion thermal velocity.For electrons, v will stand for the phase space velocity normalized by the electron thermal velocity.For both cases, ϕ is the previously normalized electric potential.When neglecting the collisional interactions H is conserved along the characteristics so that one can relate the distribution function f i (K, ϕ) to that at the sheath horizon Let us consider a source distribution function initially generated at x = 0, the SOL symmetry point where the electric potential is maximum, figure 13.Positive velocity particles evolve to the right, x > 0, for the ions along the characteristics K = K 0 + ϕ 0 − ϕ, the subscript 0 refers to the position x = 0. Conversely, the negative velocities evolve to the left, x < 0. In terms of the kinetic energy, this evolution opens a gap in the neighborhood of the axis K = 0.In the simulations, this gap in the distribution function is filled by the source term and by the collisions that act as a restoring force toward a Maxwellian distribution, therefore with non vanishing values for all values of K.One can note that the particles born within the gap, such that K ⩽ ϕ 0 − ϕ, are trapped by the electrostatic potential and cannot reach the point x = 0 along the characteristics, an example is given by the dashed black lines on figure 13.
Let us now consider the neighborhood of the sheath horizon with electric potential ϕ ∞ and given a decreasing electric potential toward the right hand side where the wall is localized, dashed region on the right hand side, see figure 14.We first only consider positive velocities since no ions are  assumed to flow back from the wall.From ϕ ∞ , black dash-dot vertical line, to ϕ < ϕ ∞ , black dashed vertical line, the distribution indicated by the shaded region is convected along the characteristics, plain and dashed blue lines.Since there is no source term in this region, the gap that is generated between K = 0 and K ϕ = ϕ ∞ − ϕ can only be fed by collisions.On a small distance, comparable to the Debye scale, we shall first assume that collisions are negligible and that no particles are generated within the gap.At position ϕ, the density can be determined by the velocity integration of the distribution function f i (v, ϕ) from v g to +∞.Here v g is such that for v < v g , the distribution f i (v, ϕ) = 0.Here v g is taken positive because we assume that the ion population is outgoing and no particles stream back from the wall.Given K g = 1 2 v 2 g , one must also have To perform the next steps we have conveniently set the lower bound of the integral to zero, which leaves the integral unchanged because of the gap without particles, and use the kinetic energy as integration variable.In the calculation, one can then take into account that the characteristics conserve δK, δK = δK ∞ and that the distribution is constant along the characteristics One then obtains the identity The density n i (ϕ) can then be determined by the values of the distribution function at the sheath horizon One can then use v(K, ϕ) = 2(K ∞ + ϕ ∞ − ϕ) to determine dn i /dϕ: Then setting ϕ → ϕ ∞ one finally obtains: We have therefore recovered the expression of reference ( [11]) and the intrinsic difficulties in handling this expression [13,14].Indeed such a result only holds if the ion population at small or negative velocity is fully depleted so that As will be shown in the following section 4.2, this proves to restrict the use of this elegant result.

Kinetic constraint on the ion outflux
The simulation results highlight several shortfalls of the constraints proposed in the previous sections to determine the location of the sheath entrance.Regarding the fluid approach, a key issue is the closure of the Navier-Stokes system equation (A5).In 1D, one finds that both the particle and momentum balance equation are exact since it only relies on the definitions of the fluid moments.The only issue is therefore the closure of the third equation and the importance of the heat flux.As discussed earlier, modifying this closure will change the value of the sound velocity, so that the Bohm criterion M = 1 is not applicable in the kinetic framework because the sound velocity is not readily determined.
In the analysis of the sheath constraint at the sheath horizon, an important issue is that of charge neutrality, quasineutrality with ε D → 0. In realistic cases, this asymptotic limit is not met since ε D is small but finite.Furthermore, when defining as ε D as ε D = λ D0 /L ∥ , we have implicitly assumed that a unique scale characterizes the SOL behavior.However, as the distance to the wall decreases, the relevant scale in the parallel direction can be seen to decrease.The effective value of ε D , which controls the variation of the electric field E = −dϕ /dx, therefore gradually increases.The picture of an abrupt change from scale L ∥ to scale λ D0 does not hold, in agreement with the idea of a smooth transition into the sheath conditions.
When plotting |ρ c | obtained in the reference simulation, figure 15 one finds four regions.In the two first regions |ρ c | ≈ ε 2 D and changes sign with the sign of the curvature of the electric potential, hence the dips where ρ c changes sign, figure 15 blue curve.In the third region, the electric potential exhibits a weak curvature: ρ c is positive and ρ c ≪ ε 2 D .It is also observed to increase exponentially toward the fourth region, in the vicinity of the wall, where ρ c ≫ ε 2 D and ρ c exhibits a sharp exponential growth by four orders of magnitude.One can note that the transition point where ρ c = ε 2 D occurs at the position x = 195.1,therefore at a distance ≈8λ D0 from the wall.The presheath exponential growth is characterized by an e-folding length of 21λ D0 , while the sharp exponential growth is characterized by an e-folding length of 0.43λ D0 .
The SOL charge density being small but finite, the key assumption ε 2 D = 0, consequently ρ c = 0, used to describe the sheath horizon and the expansion at the sheath horizon in Table 3. Location in x where the various fields switch from a SOL behavior to sheath variation on the Debye scale.The position labeled Qn corresponds to the point where ρc ⩾ ε D .The position of the transition point increases with the order of integration from that for dρc/dϕ to that of ϕ.The wall position is recalled as well as the transition point determined when using R defined in equation (15b).The position x sheath is defined as the transition point for dρc/dϕ, x sheath = 195.9.|dρ c /dϕ| and finally the point where ρ c = ε 2 D labeled Qn.As can be expected, one finds that the higher the derivative order, the further upstream from the wall one observes the sheath entrance.The sheath transition is gradual, and, when defining the transition according to the change of slope of |dρ c /dϕ|, hence x sheath ≈ 196, one finds that most plasma parameters exhibit the same behavior, the large increase of the derivatives being mostly located close to the wall.

Qn
In view of the latter discussion, one could expect that the kinetic criterion equation (14a), using equation (14b) to determine dn i /dϕ, should provide a means to determine the sheath horizon.For convenience, let us define the ion density n K determined with equation (14b) In the previous calculation, |dρ c /dϕ|, both dn e /dϕ and dn i /dϕ were computed with the derivatives of n e and n i with respect to ϕ using the moment 0 of the distribution functions to determine the densities.The computation of R(x) in equation ( 15) is therefore an alternative calculation of |dρ c /dϕ| where the derivative of the ion density with respect to ϕ is determined directly in the kinetic framework.Provided the latter calculation is correct, one expects that R ⩾ 0, to be small for x < x sheath , and to increase when x > x sheath .
One can note that the 'kinetic' criterion equation (14a ) does not address the ion and electron population on the same footing although the calculation leading to the expression of dn i /dϕ can be performed for any species.However, as sketched on figure 18, the slowing down of the electrons toward the wall generates a population of trapped electrons so that one cannot identify the distribution function at ϕ with that at ϕ ∞ .Furthermore, the electron population in the trapped region and in particular for K = 0 will govern a singularity of the integral equation (14b).In simulations, the same issue arises for the ions since collisions will tend to populate the gap region, with both positive and negative velocities so that the calculation of the ion density at ϕ when restricting the integral to positive velocities with will tend to underestimate the ion density.The kinetic formula will therefore yield n K such that n K < n i where n i is the exact value.Furthermore, the singularity of the integral for K → 0 will tend to yield dn K /dϕ > dn i /dϕ because of the finite value of the distribution function for v → 0 + .Depending on the values of the distribution function f i , the fraction f i /v 2 for v → 0 + can dominate the integral in equation (15a).This discrepancy will tend to increase as the potential difference |ϕ − ϕ ∞ | increases because of the growing possibility of generating ions in the trapped region by collisions.Consistently, as one approaches the wall, the probability of ion collisional transfer from the bulk into the gap region tends to decrease with the distance to the wall.The kinetic calculation should become more accurate.
The results of the reference simulation confirm that the kinetic formulation of the sheath boundary conditions equation ( 15) are difficult to use.First dn K /dϕ, red curve open triangles on figure 19 is found to seriously overestimate dn i /dϕ, black curve open circles for ϕ − ϕ 0 > −0.15.For ϕ − ϕ 0 < −0.15 the two curves get closer but one now finds dn i /dϕ ⩾ dn K /dϕ.Note that n i and dn i /dϕ are exact values determined using the simulation evidence.Conversely, the kinetic formulation equation ( 15) is an approximate expression, which is therefore found to be not at all accurate and yields misleading values.In the quasineutral region ϕ − ϕ 0 > −0.13, one finds R(ϕ) ≪ −1 while the direct calculation yields −dρ c /dϕ → 0 + .Indeed, the derivative of the two densities dn i /dϕ and dn e /dϕ, blue curve open circles, are found to converge for ϕ ⩾ ϕ sheath , vertical black dash-dot line.In that regime, the calculation leading to R(ϕ) does not hold.Closer to the wall, toward decreasing values of ϕ, the difference between the direct calculation dn i /dϕ and the approximate formulation dn K /dϕ tends to reduce.One can then define ϕ K at the intersection point R(ϕ K ) = 0.The potential ϕ K is found to be reached at x K ≈ 201.5, much closer to the wall than the sheath entrance x sheath , with ϕ sheath = ϕ(x sheath ), determined using the criterion on dρ c /dϕ, figure 18.

Mandatory collisional process in steady state
Let us consider the again the phase space characteristics using the variables Ksign(v) in velocity and ϕ 0 − ϕ in position.The energy conservation for the electrons leads to K = K 0 − (ϕ 0 − ϕ).We choose here the stagnation point at x = 0 as reference phase space location K 0 , ϕ 0 .We then consider the source at this location Sδ t .The chosen source for the model does not depend on the distribution function and is assumed to be proportional to a Maxwellian.We shall consider this case to illustrate the need for phase space collisional transport to achieve steady state conditions.For simplicity we shall consider the source contribution to the distribution function at x = 0 and follow the characteristics from x = 0. Given ϕ w the electric potential of the plasma at x = x w , namely the wall position, one then obtains K w = K 0 − (ϕ 0 − ϕ w ).When K 0 ⩾ (ϕ 0 − ϕ w ) one then finds K w ⩾ 0, these electrons reach the wall and are absorbed by the wall.Conversely, for K 0 < (ϕ 0 − ϕ w ) the electrons do not reach the wall and the velocity changes sign at position, ϕ such that K 0 = ϕ 0 − ϕ.The latter class of electrons are trapped by the variation of the electric potential between x 0 and x w .The variation of the electric potential is twofold, in the quasineutral SOL plasma, the electric potential drop is governed by the density drop, and consequently the plasma acceleration, and, toward the wall, the electric potential drops to ensure charge balance by enforcing equal ion and electron flux to the wall.Based on this argument the typical potential drop is ϕ 0 − ϕ w ≈ log(2 m i /m e ).In the reference simulation m i /m e = 400 so that ϕ 0 − ϕ w ≈ 3.5.For the present discussion the exact value of the potential drop is not needed, see companion paper [26].On figure 20 the characteristics starting from x = 0, ϕ = ϕ 0 are plotted.For convenience the velocity space coordinate is chosen to be Ksign(v), which allows one to track the velocity reversal.On this plot the potential drop ϕ 0 − ϕ w is set to unity.The shadowed area corresponds to the region of trapped electrons, subscript t.At the stagnation point, the critical electron kinetic energy is K t = ϕ 0 − ϕ w : electrons generated by the source at ϕ = ϕ 0 with K < K t are trapped and the electrons generated with K > K t are lost to the wall.Plotting the phase space characteristics for the electrons illustrates the kinetic behavior of the electrons generated by the source in steady state conditions.The blue line for K = K t at ϕ = ϕ 0 splits this population between the trapped and loss regions.One thus finds that all electrons generated with K < K t are confined and cannot be lost.With this physics, the particle source is always larger than the particle sink, and a steady state regime with source and sink balance cannot be achieved.In the framework of the present model, collisions are the only available physical mechanism of phase space transport that allows steady state conditions to be achieved.Trapped electrons accumulate due to the source term until collisional detrapping is efficient enough to transfer all the electrons generated with K < K t into the phase space region K > K t , as required to achieve steady state conditions.As illustrated on figure 21, for K t = 3.5 less than 1 % of the particle flux governed by a Maxwellian source is directly injected in the loss region.All the particle flux injected in the trapped region must then be balanced by a collisional flux from K < K t to K > K t .While the particle source is localized, the collisional particle transfer from the trapped to the loss regions can take place in the whole plasma.One thus finds that the electron population is split in two populations coupled by the collisions, both electron-electron and electron-ion.This particular kinetic feature cannot be captured within the fluid framework unless two electron fluids are used together with ad hoc transfer mechanisms.Furthermore, the constant evaporation of the fast particles governs a specific cooling of the trapped electron population so that one readily expects T e < T s where T s is the thermal energy of the particle source.This cooling is enhanced by the energy transfer from the electrons to the ions governed by the ion acceleration mechanism via the electric field confining the electrons.All these kinetic effects depend on the plasma collisionality and underline the fundamental difference between weak collisionality and non-collisional regimes.

Distribution functions
Two keys properties of the distribution functions have been addressed in the previous sections.First we have used the conservation of the distribution function along the characteristics at constant energy E a = K a + Z a ϕ for species a, where K a is the normalized kinetic energy K a = 1 2 v 2 and Z e = −1, Z i = 1 are the normalized charges.These characteristics account for the kinetic equation but for the collision and source terms.We can therefore use them to investigate the specific effects on the distribution functions governed by these two terms.
Let us first examine the distribution function at x = 0.These are symmetric f a (v) = f a (−v) because of the symmetry of the simulation domain.One finds that the distribution departs from the Maxwellian built with the same, density n i (x = 0) = 2.47, mean velocity and thermal energy T i (x = 0) = 0.31, blue line open circles.It is better recovered with the sum of two Maxwellians with a cold component with density n ci and thermal energy T ci ≈ 0.11 and a hot component with density n hi and thermal energy T hi ≈ 0.85 and n hi /n ci ≈ 0.11.The hot component can be related to the source thermal energy T s = 1, while the cold component that is dominant at small normalized velocity |v| ⩽ 1 exhibits a thermal energy that is closer to that of the electrons T e (x = 0) = 0.16.A similar analysis can be performed for the electrons, but the hot component T he ≈ 0.87 has a much smaller density n he /n ce ≈ 5. 10 −3 than the cold component that is consequently close to the fit by a single Maxwellian.
We investigate the properties of the distribution functions in phase space using ϕ as position coordinate and K sign(v) as velocity coordinate.The constant energy characteristics, conserved energy K e − ϕ for the electrons and K e + ϕ for the ions, are then straight lines.For the electrons, the contour plots of log 10 (f e ) are compared to the constant energy characteristics K e = K e0 + (ϕ − ϕ 0 ).On figure 22, the phase space is restricted to the small velocities K e ⩽ 1.Note that the electric potential is inverted for convenience so that the particles with positive velocities head from the symmetry   contour lines for the values −2.5, −1.5, and −1 are the blue curves with open circles.For the electron trapped in the electric potential variation, the contour levels are −0.09and 0.13 and the contours are the blue curves with open squares.The constant energy characteristics correspond to the black lines with closed circles.They are computed to intersect the contour curves at K e = 0 for the trapped electrons and otherwise at ϕ − ϕ 0 = −0.2.One finds little departure between the contour curves and the constant energy characteristics for K e < 1 except for ϕ → ϕ 0 and K e sign(v) < −0.5.This effect is clear when stepping to K e sign(v) < −1, figure 23, where the value of the contour line range from −6 to −3 with steps of 0.5, blues lines open circles.The departure between the contour curves and the constant energy characteristics can be attributed to the effect of the source terms that lies to the left of the vertical dash-dot line and is maximum for ϕ → ϕ 0 .This source is found to induce a broadening of the distribution function by adding particle with thermal energy T s = 1.This effect is observed to be much more important than the changes governed by the energy conservation, figure 23.One can also note that a distortion with respect to the constant energy characteristics occurs for −0.075 > ϕ − ϕ 0 > −0.13.In this range of values the source is null and only collisions can explain this departure.One can also remark that the source should govern the development of a hot Maxwellian contribution to the distribution functions, with comparable properties for both electrons and ions.However, this hot component at x = 0 is observed to be much smaller for the electrons than for the ions with a ratio ≈5. 10 −2 figure 25.One also finds that the observed density of the hot ion component is typically a factor 0.63 smaller than would be obtained with only the source term at work.Finally, for both species one observes that the thermal energy of the hot component at x = 0 is ≈0.85, therefore smaller than T s .These features highlight the role of collisions, and as discussed in section 4.3, the mandatory evaporation of the hot electrons sustained by collisional transfer.
For K e sign(v) > 1, one can also observe the broadening of the distribution function governed by the source term, and, prior to the source region, by collisions.However, in this part of the distribution function, these changes are larger but comparable in magnitude to the variation along the constant energy characteristics.We now step to the properties of the ion distribution function.We first consider the small kinetic energy region of the phase space with K i < 0.2 and the contour of log 10 (f i ), figure 25.The contour levels with blue curves open circles range from −2 to 0.25 with steps of 0.25.The nearly closed contour line is governed by the variation of the ion density.The constant energy characteristics for K i sign(v) < 0, black line closed circles, is very different from the contour lines that are consistent with the effect of the source term for a vanishing ion distribution with negative velocity and ϕ < −0.075, between the source region and the wall.The contour lines also appear to pinch toward the K i = 0 line, then opening a region of very small values of the distribution function for K i < 0.05 as ϕ − ϕ 0 → −0.2 one finds that the contour lines tend to be parallel to the constant energy characteristics for the smaller values of log 10 (f i ), with an important mismatch at the larger values of log 10 (f i ).For the larger values of K i sign(v), K i sign(v) > 0.2, figure 26, one finds that the contours are aligned on the constant energy characteristics for typically ϕ − ϕ 0 < −0.13.For ϕ − ϕ 0 > −0.13, one recovers a behavior that is comparable to that reported for the electrons, a broadening governed by the source term on the lefthand side of the vertical dash-dot black line, ϕ − ϕ 0 > −0.075, and a broadening that can only be attributed to collisions for −0.075 > ϕ − ϕ 0 > −0.13.

Motivation
In this paper we have addressed the self-organization of a plasma with perfectly absorbing boundary conditions.The problem we have in mind is the physics of the SOL plasma of magnetic confinement devices dedicated to fusion.We have also restricted the analysis to a 1D geometry, typically in the direction parallel to the magnetic field, and assumed the location of the absorbing boundary condition to be fixed, as when the plasma interacts with a solid.However, the problem at hand is not specific of fusion plasma conditions and our results are therefore more general.Under such conditions a boundary layer develops in the vicinity of the wall.The so-called sheath then stands between a bulk plasma that remains quasineutral and the wall where the plasma promptly recombines.The sheath is well known to be a region with positive charge density extending on several Debye lengths and where the electric field confining the electrons becomes large.Our specific interest is the self-organization problem of the SOL plasma sustained by particle momentum and energy source terms and in contact via the sheath boundary layer with a perfect particle, momentum and energy sink, section 2.1.The sheath transmission properties together with the sources then determine the SOL plasma properties, in particular the plasma momentum flux (also called total plasma pressure) and the plasma electric potential.Our approach does not aim at addressing the sheath physics as a stand alone problem but rather the particular balance between the plasma source and the sheath transmission.A specific interest of this analysis is to define appropriate boundary conditions for simulations dedicated to plasma turbulent transport where quasineutrality is enforced and the Debye scale not resolved.The simulation effort is twofold, a large body of SOL simulation codes use the fluid representation, and, more recently, the gyrokinetic framework has been used for SOL simulations.The problem is then to define appropriate boundary conditions such that the SOL properties are recovered without having to address the physics on the Debye scale.At the present stage, atomic processes are ignored since our first goal is to address the so-called sheath limited regime, hence an SOL regime with low density and high temperature plasma that tends to minimize the role of atomic processes.
To illustrate our theoretical analysis we have used the results of a reference kinetic simulation with the 1D-1V VOICE code [33] evolving a two species plasma, electrons and singly charged ions, and where penalized wall conditions are used, see the companion paper [26].In the present work, we have considered that the plasma source, mostly governed by crossfield turbulent transport is localized and symmetric, so that the stagnation point with zero mean velocity for either species stands in the middle of the plasma region, section 2.2.In this work, we furthermore assume that the source injects hot particles, hence a particle and energy source, with no heat, charge or momentum source.With this specific setting, the control parameters are the particle source, its shape and magnitude, the thermal energy of the injected particles, and the plasma collisionality.The latter defines in fact the density chosen for the normalization, while the thermal energy of the source determines the energy normalization.A Maxwellian particle source is chosen, which minimizes the number of free parameters, and is more convenient when stepping to the fluid framework.

Predictions in the fluid framework
We consider a quasineutral plasma with zero electric current.The fluid model we first address is the standard Navier-Stokes model for particle, momentum and energy conservation, either for each species or summing the latter to obtain the plasma conservation equations used instead of the ion conservation equations.The plasma equations, considering the density, particle flux, the total momentum flux with null charge density and electric current are found to be independent and a closed system when neglecting both electron and ion heatfluxes.We use this system with three degrees of freedom, the three fields plasma density, particle flux and momentum flux, to determine analytic steady-state solutions.One finds that the mean ion and electron velocity as well as the plasma thermal energy T = T e + T i vary with the particle flux, itself fully determined by the particle source.In steady-state, the system can be recast to exhibit three invariants, the plasma momentum flux Π, the exchange between energy source and energy flux and the exchange between particle source and particle flux.Since the energy source is governed by the particle flux, the latter is the only effective control parameter of the chosen model in the fluid framework.Consistently, the analysis of the present work determines the Mach number in terms of a single control parameter, A s equation (7c), which is bounded and its absolute value chosen to vary between 0, where the Mach number is null, and 1, where the Mach number is equal to 1.This control parameter increases with the particle flux and decreases with Π + ηT e , where η = n i − n e describes the departure from quasineutrality, section 3.1.This calculation follows similar steps as previously published for the simpler isothermal closure [4].
The plasma-sheath transition is found to exhibit a twofold behavior.In the quasineutral SOL, the source terms governs the increase of the control parameter A s , so that the subsonic Mach number increases.At the sheath horizon, the control parameter goes through a maximum and decreases into the sheath, since Π + ηT e is no longer constant and varies with charge separation η > 0. The Mach number then further increases into the supersonic regime.This bifurcation from a subsonic to a supersonic regime, with no discontinuity, is made possible because of the departure from quasineutrality.This result for the transition at the sheath horizon is supported by the matching technique [4] between the sheath and presheath variation, appendix A.4.
However, the above discussion assumes the sound velocity to be known so that the Mach number can be defined.In this paper, the sound velocity is determined by analyzing the sound wave dispersion relation in appendix A.2, using as independent variables the density, the particle flux and the total plasma momentum flux-the latter being akin to the energy density for a 1D model.This calculation is completed by analyzing the possible discontinuity of the fluid equations of the quasineutral plasma for the variables, density, mean plasma velocity and plasma thermal energy in appendix A.5.One shows in particular that the two approaches lead to the same critical mean plasma velocity.The sound velocity for the non-isothermal solution is then found to be c s = √ 3 V with purely convected energy, compared to c s = V with the isothermal model, V = T/m i .Both predictions depart from the result of the kinetic simulation where the transition into the sheath is estimated to occur at 1.5 V.The latter specific result can be used to determine an effective polytropic index, which would then imply that heat transport must be accounted for.However, the kinetic simulation is not consistent with a particular polytropic index since the simulation output yields a varying polytropic index throughout the simulation domain.This shortfall can be expected since the use of the polytropic index is based on the assumption of prompt thermodynamic equilibrium induced by small scale turbulence, a transport regime that clearly does not make sense for parallel plasma transport.
In the chosen fluid framework, heat-fluxes are assumed to be null.The transport equations then enforce that all fields, including the mean plasma velocity and Mach number, vary within the source region and then remain constant up to the sheath entrance where the departure from quasineutrality governs specific small scale variations.The mean plasma velocity must therefore reach the sound velocity at the end of the particle and energy source region, and then remain at Mach 1 up to the sheath region.This prediction and other fluid predictions are compared to the reference kinetic simulation and qualitative agreement is found.However, key differences are observed.First, the heat-flux for the electrons is comparable to the convected energy flux and contributes to the total energy flux impinging onto the wall.Second, the variation of the steady state solution is not restricted to the source region.The energy exchange between the species and between the convected and heat-flux components continues from the end of the source region up to the wall.Similarly, the mean particle velocity continues to increase monotonically between the end of the source region and the wall.It remains smaller than √ 3V but larger than V so that a Bohm criterion defining the sheath entrance with the sound velocities obtained with usual fluid closures is not operational.These differences suggests that the kinetic simulation involves more degrees of freedom than selected when handling the Navier-Stokes system.Furthermore, since the density, mean fluid velocity and thermal energy are exactly defined by the three first moments of the distribution function-namely the density, particle flux and momentum flux-the link to other degrees of freedom within the distribution function must involve the heat-flux.Kinetic effects, such as collisional exchange between convected energy flux and heat-flux, can only be addressed within the fluid framework with an appropriate model for heat-flux transport.
Given these shortfalls, a non-collisional closure is investigated to determine the heat-flux in terms of the lower moments, namely the density, the particle flux and the total plasma momentum flux, appendix A.3.The analysis of this projection in Fourier space indicates that the heat-flux exhibits a contribution along the particle flux with a real proportionality coefficient.Conversely, the proportionality coefficients for the density and total momentum flux are imaginary, driving a damping behavior and consistency with Landau damping.This non-collisional closure is interesting insofar that the heat-flux contribution proportional to the particle flux provides an explanation for the heat-flux contribution to the total plasma energy outflux.Such a contribution also modifies the sound velocity, introducing an asymmetry between the co and counter sound waves with respect to the mean plasma velocity.The counter propagating waves are found to be slower, so that a shock wave generated when the flow velocity cancels the sound velocity, hence with standing waves in the wall reference frame, occurs at a lowered flow velocity, appendix A.3.This trend is consistent with the observations from the reference kinetic simulation.
Closures using higher moments of the fluid hierarchy are also discussed, appendix A.2.5.One finds that the sound velocity and the discontinuity of the steady state balance equation are in agreement but depend entirely on the choice made for the closure.This arbitrary feature of the result do not allow one to determine a particular Bohm criterion defining the sheath entrance.Predictions of such fluid models using a closure at higher moments can be made more accurate provided one uses evidence from kinetic simulations to constrain the closure.However, they do not offer a predictive capability that is fully independent from the kinetic computation.The kinetic simulation guideline therefore remains crucial in determining the transition from the quasineutral SOL plasma into the sheath region, and consequently to define the sheath entrance.To put things bluntly, one can expect that the Mach 1 Bohm criterion holds at the sheath horizon, but one has no means to specify the sound velocity, so that this fluid criterion is not operational.

Kinetic criterion at the transition into the sheath
In the kinetic framework, the Mach 1 discontinuity or bifurcation, and therefore the Bohm criterion, does not emerge naturally as underlying a particular behavior of the distribution function.A specific kinetic criterion to define the sheath horizon has been proposed by Harrison and Thompson [11].The constraint is based on an integral involving the ion distribution function.It has been regarded as the most relevant criterion to determine the sheath horizon in the kinetic framework [13][14][15].We have revisited the derivation of this criterion in section 4.1 and recovered that the chosen ion distribution integral stands for the derivative of the ion density with respect to the electric potential.The criterion is obtained when expanding the electric field energy, typically E 2 in terms of an electric potential difference.At the sheath horizon, it is assumed that the plasma is neutral n e → n i so that the constraint E 2 ⩾ 0 enforces that the derivative of n e with respect to the electric potential is larger than that of n i .At this stage, the criterion involves the densities and is not really specific of either the kinetic or the fluid framework (although handling the density is more a fluid point of view).
To determine the proposed criterion, the derivative of the electron density is analytically determined, usually for adiabatic electrons, while the derivative of the ion density is computed using the ion distribution function and the characteristics at constant energy E = K + ϕ.The transfer of kinetic energy K to electric potential energy ϕ along the characteristic then allows one to compute the ion density derivative with respect to ϕ.However, this calculation exhibits a divergence for K → 0 that is only alleviated if the ion distribution function tends to zero fast enough when K → 0. Furthermore, the use of the constant energy characteristics to determine the ion distribution function is only appropriate when both collisions and source/sink terms are not taken into account.In the kinetic simulation, one finds that the Harrison-Thompson formula to compute dn i /dϕ departs from the actual value in most of the plasma region and only becomes similar, but not identical, to the latter toward the wall.One can also note that applying this formalism to the kinetic simulation is not appropriate because collisions, however weak, will play a role when the distribution function exhibits strong variations in velocity space, which is the case for the ion distribution function toward the wall.
The Harrison-Thompson criterion must therefore be adapted to actual plasma behavior as displayed in the kinetic simulation, see section 4.2.In particular, one must consider that quasineutrality does not enforce a null charge density.We must then identify a change in reference scale from that of the domain size L ∥ to the characteristic Debye scale λ D0 .Monitoring the rapid change in variation scales then allows one to identify turning points from one asymptotic behavior to the other.However, depending on the chosen field this transition occurs as different distances from the wall.As can be expected, this sensitivity increases when stepping from the electric potential to the electric field and then to the charge density and the derivative of the charge density with respect to the potential.The gradual transition from quasineutral to Debye scale dependence does not provide a unique definition of the sheath entrance.However, a suitable criterion for the sheath horizon appears to be the change of variation scale of dn e /dϕ − dn i /dϕ.It combines the fact that it is a precursor, responding when the other field have not yet changed behavior, and that it is rather straightforward to compute.For the reference kinetic simulation, dn e /dϕ − dn i /dϕ exhibits an exponential growth in the quasineutral region with e-folding length ≊20λ D0 , which decreases to ≊0.5λ D0 within the sheath.From the sheath entrance to the wall, therefore over a distance of ≊7λ D0 , the normalized charge density then increases from the small value expected in the quasineutral region to a value of order unity.

Role of collisions and distribution function remodeling
An important feature revealed by the model problem used in this paper is the key role of collisions, and the coupling it drives between the various degrees of freedom involved by the kinetic solution, section 4.3.We have shown that the electron population at the symmetry point, the stagnation point, can be split into trapped electrons with kinetic energy smaller than the electric potential drop from the stagnation point to the wall and electrons with large enough kinetic energy to be able to reach the wall.The kinetic source term injects hot particles into the plasma.The source distribution function will also be split into injected electrons that are confined by the electric potential drop and those that can stream to the wall, and lost to the wall.In the model, collisions are the only mechanism that can transfer electrons from trapped to streaming.These are found to be mandatory to achieve steady state.The evaporation mechanism driven by the collisions governs a cooling down of the electrons to low thermal energy.Following the constant energy characteristics from the stagnation point toward the wall, and comparing these to the contour lines of the distribution functions, allows one to identify the role of the source and collision operators.One finds relative agreement toward the wall where the electric field becomes large and can be expected to be the drive for the changes in the distribution functions.For the electrons, dominated by a cold component at the stagnation point, the source effect is mostly observed at large kinetic energy, yielding a small density for the hot component.This source effect appears to be balanced by collisional cooling in the source region and close to the source region.Regarding the ions, collisional cooling is also observed since fast ions are lost faster than the slow ions, but this effect is less pronounced than for electrons.Consequently, the density of the hot ion component is larger.The departure between the contour lines and the constant energy characteristics also indicates that collisions play a role everywhere.One also finds that they tend to fill the energy gap K i → 0 so that the distribution function is small but tends to a finite value as K i → 0.
The collisional transfer that remodels the distribution function is shown to be mandatory to balance the cooling by evaporation and achieve steady-state.It is a key process in the kinetic simulations with the present model.When considering actual plasmas in fusion devices, therefore higher dimensions in space, it can be argued that the distribution function remodeling can also be achieved by cross-field turbulent transport.This physics is however partly accounted for by the chosen source term.Turbulent cross-field convection from neighboring field lines within the SOL should yield a source with similarly distorted distribution functions.Choosing localized Maxwellian source terms characterizes an exchange with the core plasma in the midplane region, where turbulent transport is ballooned.Although this simplified setting of our 1D model cannot capture all the complexity of cross-field transport, one can argue that turbulent remodeling of the distribution function is at least partly taken into account.The collisional remodeling we report should therefore remain a key process of the SOLsheath self-organization in actual experiments.

SOL-sheath self-organization
A crucial feature of the SOL/sheath self-organization is the plasma pressure that is achieved given the particle and energy source on the one hand, and the sheath constraints on the other hand.In fact, in models for plasma-wall interaction, such as the two-point model [35], the plasma total pressure is linked to the upstream particle and energy flux and to sheath constraints [43,50].Let us use the subscript div to localize the sheath entrance, as standard when addressing plasma-wall interaction.The expressions for the energy flux Q div and particle flux Γ div at the sheath entrance are then: Q div = γΓ div T e,div with Γ div = n div V div M div .We define here m i V 2 div = T div , T div being the sum of the electron and ion thermal energies at the sheath entrance.Given V div one defines M div , the mean plasma velocity divided by the thermal velocity V div .Note that M div is not the Mach number but corresponds to M V (defined and used in section 3.1) at the sheath entrance; it does not depend on the definition of the sound velocity.Without losses along the field line, Q div = Q up and Γ div = Γ up , where Q up and Γ up are the injected energy and particle flux, the subscript up standing for upstream.Given these expressions, one then finds: To determine this relation we have used the fact that Π is constant along the field line, hence Π div = Π.The plasma total momentum flux Π is therefore proportional to the injected particle and energy fluxes via ) and A s,div is defined in equation (7c).The sheath constraint enforces that |A div | or |A s,div | take their maximum value equal to 1 at the sheath entrance and consequently that Π takes a minimum value, typically Π ⋆ .The result of these two closures suggests that the SOL-sheath self-organization process enforces the minimum possible value for the plasma pressure.In that respect, the present result is akin to selforganization processes determining turbulent states [51,52] or transport [53] as the extremum of a given functional.In the reference simulation, the plasma total momentum flux Π is observed to be constant along the field line up to the sheath entrance.The self-organization feature is recovered, however the value of Π is 10% smaller than predicted by the Navier-Stokes model.

Kinetic features of SOL-sheath self-organization
The importance of kinetic effects is highlighted by the distortion of the distribution functions away from Maxwellians.As discussed above, the collisional transfer is also an important kinetic effect.As also recalled when considering the Navier-Stokes fluid equations, one finds that the kinetic effects are described by the heat-flux and the inter-species collisional equipartition toward identical mean velocity and thermal energy.All the other terms appearing in the equations are definitions, therefore not specific of the kinetic framework.In the reference kinetic simulation, we find that the heat-flux amounts to a fourth of the plasma energy flux and close to half the electron energy flux.
The heat-flux is found to build-up in the source region and then to be sustained, although it slightly decays, up to the wall, where the heat-flux therefore contributes to the flux impinging onto the wall.As discussed above, the non-collisional closure indicates that this heat-flux exhibits a component that is proportional to the particle flux.Unlike the result of the collisional closure that describes the heat-flux as a diffusive process exchanging energy with zero net particle flux, we find in this work that the heat-flux contributes to the energy flux convected with the particles out of the plasma.Furthermore, it appears that this contribution must be coupled to higher moments of the fluid hierarchy to recover a Bohm criterion.We find therefore that defining the heat-flux as an energy transfer with no net particle transfer, which is consistent with the result of the collisional closure, does not account for our observations in the kinetic regime.A better understanding of the physics of the heatflux, the skewness of the probability distribution functions, its dependence on collisionality and connection to the distortion of the distribution functions, appears to be important to assess the kinetic features of parallel transport.This issue is further discussed in the companion paper dedicated to the simulations [26], together with the heat transmission factor.are related Π a = 2E a .One can thus remark that the Navier-Stokes system equation (A1) is closed but for the term q a and eventually the collisional exchange terms.Possible closures are either setting these terms to zero, or reducing these terms to functions of the three first moments n a , Γ a and Π a , and setting the other contributions to zero.Conversely, any specific kinetic feature, by definition not taken into account by this fluid description, must translate into a behavior of q a and the collisional exchange terms that cannot be approximated by zero or the chosen functions of the first three moments.
A.1.2.Two species plasma conservation equations.For a two species plasma with electrons, charge e e = −e and mass m e , in singly charged ions, charge e i = e and mass m i , one can conveniently split the system into equations that are related to charge and equations addressing the mass center.The latter system will be named plasma since it combines the electrons and the ions.Because of the large mass ratio m i ≫ m e , charge and electron balance equations evolve on the electron time scale, while the plasma balance equations evolve on the ion time scale.Regarding particle balance, we shall consider the ion particle balance and the charge balance, defining the charge density ρ c = e(n i − n e ) = eη, η is the density difference η = n i − n e , and the total electrical current j = e(Γ i − Γ e ), therefore: We assume here that there is no charge source, hence S i n = S e n = S n .The total plasma momentum and energy conservation equations are completed by those for the electrons.
In the chosen model we have furthermore assumed that the energy source is identical for the ion and electron channel S e E = S i E = S E , hence the factor 2 in equation (A4d).
A.1.3.Quasineutral and large mass ratio limits.In the quasineutral limit one assumes that the charge density ρ c can be set to zero, therefore ignoring the difference between electron and ion densities.This regime is enforced by the Poisson equation (2b) in the asymptotic limit ε 2 D → 0. The charge balance equation (A3b) then enforces that ∇ s j = 0, hence a 1D current that is constant along the field line.Although a non vanishing current can exist, it must be sustained by a potential difference.Assuming the wall to be grounded, so that no potential difference is generated within the wall and between the two ends of the field line, enforces the electrical current to be null, j = 0.For the two species plasma, the small departure from neutrality end forces n e ≈ n i = n, η ≪ n i , so that the j = 0 yields Γ e ≈ Γ i = Γ.The quasi-identical density and particle flux enforces u e ≈ u i .Given the large mass ratio m i /m e ≫ 1, we assume the ordering m i n i u 2 i + m e n e u 2 e ≈ m i n i u 2 i and Π e ≈ n e T e .The plasma balance equations are then: where As written here, we have not enforced ρ c → 0. One is led to consider this term when one addresses steady state conditions on the time scale of the electron evolution.This also enforces ∇ s j = 0, hence j = 0 and equality of the fluxes.Provided m i /m e ≫ n i /n e ≈ 1, one recovers the ordering m i n i u 2 i + m e n e u 2 e ≈ m i n i u 2 i and Π e ≈ n e T e and the plasma balance equations equation (A5).However, because of this coupling to ρ c , this system must be completed by the electron conservation equations and the Poisson equation relating the electric potential ϕ to the charge density ρ c .More generally, one can state that whenever j ≈ 0, meaning |j| ≪ e|Γ i | and provided |ρ c | ≪ en e m i /m e , one recovers equation (A5).For an open system, with volumetric particle source, hence depending on the curvilinear abscissa s, and boundary particle sink, the known particle source determines the particle flux.For a symmetric source with respect to the middle of the box at s = 0, the particle flux is antisymmetric and null at the stagnation point s = 0.By definition the particle flux is proportional to the plasma density n and ion mean velocity u.In the region with vanishing particle source, the particle flux is invariant.
The momentum balance is enforced by the fact that no external source is applied to the plasma, and that binary collisions ensure momentum conservation, both in the collision process and in the mean field momentum transfer via the internal electric field.The plasma momentum flux is therefore invariant.The simplification brought to the expression of the momentum flux in the quasineutral limit and when neglecting the electron inertia then allows one to express the momentum flux in terms of the plasma density n, ion mean velocity u, equal to the plasma mean density, and plasma thermal energy T.
Finally energy conservation, with volumetric energy source and energy sink at the boundary, yields the third invariant.For a convective energy source, the energy source being proportional to the particle source, the symmetry of the source region enforces that the energy flux is antisymmetric and null at the stagnation point s = 0. Neglecting the heat flux, the energy flux is then proportional to the plasma density n, plasma mean velocity u and plasma thermal energy T. In the region with vanishing particle source, and hence vanishing energy source, the energy flux is invariant.
When restricting the fluid model to the Navier-Stokes equation for the plasma, one finds that one can replace the three fluid unknowns n, u and T by the three fluxes Γ, Π and Q, which are invariant when the volumetric sources are null.This property can be used to analyze the accuracy of the kinetic simulations.

A.2. Plasma sound waves
A.2.1.Sound waves in quasineutral plasma.As discussed in section 2.1 and illustrated by figure 4, an important issue when using the fluid framework to discuss the physics of the SOL plasma self-organization is the sound velocity.To address this point we consider fluctuations of various moments of the distribution functions Then defining w = v − u a one defines: The moment M ℓ is related to the lower moments N j , j ⩽ ℓ, and the moment N ℓ to the lower moments M j , j ⩽ ℓ Let us consider the fluid hierarchy for species a: The mass center combination of the moments m i M For particle conservation ℓ = 0, the order 0 moment M 0 is the density, the order 1 M 1 is the particle flux Γ.For ρ c = 0 and j = 0 these are identical for the two species.Furthermore the particle source term is identical for both species to conserve charge, and collisions conserve particles, therefore equation (A12a) leads to equation (A5a).
For momentum conservation, moment ℓ = 1, collisions conserve the total momentum, the electrical force is proportional to the charge density ρ c = 0, and we have assumed no momentum source.Equation (A12a) therefore leads to: Similarly for kinetic energy conservation ℓ = 2, the Joule heating is proportional to j = 0, collisions conserve the kinetic energy, and the energy source is identical for each species, therefore: For the higher moments, the left hand side of the moment equations keeps the same structure as found in equation (A12a).However, on the right hand side the contribution proportional to the electric field, the collisional exchange and the source contribution will not be null.In the large mass ratio limit m i ≫ m e , the system for the three first moments then takes the form: capacities of thermodynamics.For an ideal mono-atomic gas in a space of dimension d, σ V = d/2 and σ p = σ V + 1.The modified expressions of Π and Q c are then: With the chosen zero heat flux closure q = 0 we must express Ignoring the contribution of the heat-flux q to the plasma energy flux Q, one can determine the three coefficients a π , a Γ and a n : Here the prime notation indicates that the total energy flux is replaced by the convective energy flux in the definition of the closure coefficients.Using equation (A20a) one then obtains: With these coefficients the constraint equation (A17a) is satisfied and c = u is therefore one of the roots of the dispersion equation.The two other roots are then determined by a quadratic equation, see equation (A17b).
The phase velocity of the sound wave is computed to be the sum of a Doppler shift proportional to the fluid mean velocity (γ p − 1)u/2, with γ p = σ p /σ V , plus or minus the sound velocity c 2 s = 2σ p V 2 + u 2 (γ p − 3) 2 /4.If one enforces the Doppler shift velocity to be the mean velocity u, then one must have γ p = 3 identical to the perfect gas adiabatic index in 1D.One then finds that the sound velocity does not depend on the mean velocity u, c − u = ± 2σ p V.
(A22b) Therefore c 2 s = 2σ p V 2 .The condition to achieve a standing shock wave c = 0 is met for u = c s .For the standard ideal gas values σ p = 3/2 and σ V = 1/2, hence γ p = 3, one obtains: Referring to the simulation results, figure 4, one finds that the sheath entrance is observed for u = 1.505V, which then requires 2σ p = 2.265 < 3, approximately 25% smaller than the ideal gas value.Accordingly one must have 2σ V = 0.755 < 1.
The ideal gas values of the coefficients σ V and σ p , namely σ V = 1/2, σ p = 3/2 are those obtained when computing the moments of the kinetic equation.The departure of σ V and σ p from these values indicates that the chosen closure is not appropriate to accurately determine the sound velocity.An alternative closure is addressed in the following, appendix A.2.4.
A.2.4.Sound wave velocity with polytropic closure.In many papers the sound velocity is determined in the polytropic framework.The latter assumes a closure of the form: where γ p is the polytropic index, which can be different from the adiabatic index.This state equation is used to close the first moment equation.Using equation (A12b) with zero source and equation (A12c) we then have: According to equation (A10a) one can identify M 2 to: where we have taken into account N 1 = 0, N 0 = M 0 and u = M 1 /M 0 by definition.According to the polytropic closure, one has: a perfect match between the electron density response of the fluid and kinetic framework is not possible and we restrict here the match to the first terms of an expansion in Ω, hence for Ω → 0 and c = ω/k ≪ V e .
In the kinetic framework one considers the linearized Vlasov equation for the electrons, the distribution function f e being split into fluctuations f e and equilibrium f e contributions: e m e ∂ s U ∂ v f e = 0 (A39a) so that one obtains for the Fourier modes: The equilibrium contribution is chosen to be a shifted Maxwellian: where w 2 = (v − u) 2 /(2V 2 e ) and such that ∂ v f e = √ 2wf e /V e , therefore: with ε u = u/V e .The definition of y proportional to the absolute value of Ω − ε u anticipates on the integration over w in the complex plane and the continuity requirement when Ω − ε u changes sign.Integrating over w, one then finds: The Fried and Conte function Z is defined by: w − y one obtains: The Taylor expansion of the Fried and Conte function to first order is Z(y) ≈ i √ π − 2y so that the kinetic relation between the density and electric fluctuations is: Linearizing with respect to ε u one then obtains: To identify the fluid and kinetic result one must now expand equation (A38c) Identifying the kinetic and the fluid expansions then yields the following relations: Given 2R 2 i = −π these expressions can readily be simplified: At order 0 in ε u , a n = −a π and both are complex while the order 1 contributions to both a n and a π are real.Conversely, at order 0 a Γ is real while the order 1 is complex.The complex contributions account for the Landau damping.More interesting for our purpose is the contribution A Γ that leads to Q e = T 0 Γ2/(4 − π) ≈ 2.33T 0 Γ this contribution is larger than One therefore recovers that the charge density within the sheath is positive n i ⩾ n e and that the charge density ative with respect to δϕ has the same sign as δϕ and is therefore negative.Taking into account that ϕ = ϕ ∞ + δϕ, one can recast the latter property as: This constraint holds within the sheath and the equality being satisfied at the sheath horizon z → +∞.
To perform the expansion we have assumed that Π is constant in the neighborhood of the sheath horizon, therefore neglecting the electric force −ρ c ∇ϕ compared to the variation of ∇Π qn and T e ∇ρ c , T e being assumed constant.The constraint of balancing the variation of Π qn by that of ρ c T e , equation ( A57a) is satisfied at order 1/z 4 so that we have used Π constant up to order 1/z 5 with no constraints on the higher order terms of the expansion, terms of order 1/z 6 and higher.Since ρ c is of order 1/z 4 and ∇ϕ of order 1/z 3 , one finds that −ρ c ∇ϕ of order 1/z 7 .This term only contributes to the variation of Π at order 1/z 6 which is not addressed by the analysis.One finds therefore that the variation of Π governed by the electric force −ρ c ∇ϕ can be neglected at the sheath horizon so that Π is constant within the required precision of the calculation.It is also interesting to underline that taking into account the departure from quasineutrality on the Debye scale, there is no discontinuity at the sheath horizon.Furthermore, the high order expansion in 1/z suggests that the sheath horizon is smooth with little variation.The strong gradients are found to develop for z → 0 toward the wall location.
With this analysis of the plasma variation at the sheath horizon we have found two constraints: u 2 ∞ = 3V 2 ∞ , which corresponds to the Bohm criterion given the sound velocity c 2 s = 3V 2 equation (A56a), and dn e /dϕ ⩾ dn i /dϕ equation (A58c).Using the control parameter A s equation (7c), one can readily show that these two constraints are equivalent in the fluid framework.At the sheath horizon, the particle flux Γ is constant, therefore dΓ/dϕ = 0.One then finds that the sign of dA s /dϕ is opposite to the sign of dΠ qn /dϕ.According to equation (A51a) and given Π = constant the sign of dΠ qn /dϕ is identical to that of dρ c /dϕ Combining these expressions one finds that the sign of dρ c /dϕ is determined by the sign of 1 − M 2 .The sheath condition M 2 ⩾ 1 is therefore equivalent to the condition dρ c /dϕ ⩽ 0, the equality being met at the sheath horizon.As expected and discussed in the literature, these two fluid constraints are found to be equivalent.The understanding of the transition into the sheath regime that stems from the matching of the SOL and sheath behavior is governed by the variation of the control parameter A s .In the SOL, the region of quasineutral plasma A s increases with the particle flux Γ with Π qn constant.When A s reaches its maximum value A s = 1, M 2 = 1, at the sheath horizon, the variation of A s become negligible with respect to the increase of Π s governed by the departure from quasineutrality.The control parameter A s then decreases and the fluid can access the supersonic regime without generating a discontinuity.At the sheath horizon, for z → +∞ where z is linear in s with origin at the wall, no discontinuity is found because of the typical dependence on z −2 .According to the behavior at the sheath entrance, the expansion in 1/z for large values of z suggests that the discontinuity, if any, will take place at the wall location as z → 0. By construction of the sheath entrance, the behavior at this point is in the continuity of that of the SOL plasma, therefore similar to that in the SOL plasma.It is only the starting point of a gradual increase of the gradients that become maximum at the wall position.all quantities of interest, therefore both the Mach number and d s η only depend on A s .In the reference simulation, 0 → 0 so that in the fluid framework one expects the sheath entrance to correspond to the point where d s η departs from zero.However, the simulation evidence indicates that other control parameters than A s are at play since the mean velocity varies while A s is constant, see figure 6.This suggests that the actual sheath entrance criterion is a change in the scale dependence of d s η rather than the point where d s η departs from zero.One must stress here that the actual quasineutrality limit η → 0 departs from the condition η = 0 that is enforced, the latter requiring d s η = 0 while the former only requires |d s η| → 0. A criterion base on the change in the scale dependence of d s η is therefore in better agreement with the physics, and, quite consistently, much easier to implement when analyzing the simulation evidence.
We have found that a possible discontinuity of the derivative of the moment M 0 = n at Mach = 1 is resolved when enforcing a criterion on d s η, the gradient of the charge density.This criterion governs a smooth transition from subsonic to supersonic flow regime.The understanding of the sheath as a standing shock wave corresponds to the asymptotic limit L ∥ /λ D0 → +∞.When addressing charge separation on the Debye scale λ D0 , hence for finite values of L ∥ /λ D0 , the discontinuity is resolved.However, the scale of the gradients changes from typically L ∥ to λ D0 ≪ L ∥ within the sheath.One predicts therefore a smooth transition into the sheath region, the gradients only steepening as one gets closer to the wall limit, in agreement with the simulation evidence, see figure 3.

Figure 1 .
Figure 1.Sketch of the geometry: shape and location of the source terms, plain blue line, versus the curvilinear abscissa s, here normalized by the Debye length.The vertical dashed lines indicate the limit of the source region centered on x = 0, vertical dash-dot line.The wall locations, hatched regions, bound the plasma.
. The length scale L u i is the blue curve with open triangles.It increases in the source region, rolling over to a range of values of order L ∥ .

Figure 2 .
Figure 2. Growth of the particle flux Γ, normalized by Γ 0 , blue line closed blue circles, left hand side scale.The shape of the source term is identical to that of the particle source Sn, black line open black circles, right hand side scale.The vertical dashed line indicates the limit of the source region and the vertical dash-dot line the transition into the sheath region.

Figure 3 .
Figure 3. Variation of the characteristic scales Ln i for the ion density, black curve open triangles and Lu i for the mean ion velocity, blue curve open circles.The vertical dashed black line is the right hand side limit of the source region and the black vertical dash-dot line identifies the sheath region to its right hand side.For the sake of comparison, the source scale Ls ≲ 100, horizontal dash-dot line, and the size of the simulation domain L ∥ , horizontal dashed line, are also indicated.

Figure 4 .
Figure 4. Profile of the mean plasma velocity u i ≈ ue, blue line open circle symbols, and of the reference velocities V = √ T/m i , black line open head-down triangle symbols, and cs = √ 3 V, black line open head-up triangle symbols.

Figure 6 .
Figure 6.Profile of the mean plasma velocity u i ≈ ue, and of the effective sound velocity 1.505 V, left hand side axis.Profile of As, right hand side axis, As > 1 beyond x ≈ 38, vertical black dash-dot-dot line.u i blue line open circle symbols.1.505 V black line open squares.As blue line closed circles.On this figure, the maximum of x is the sheath entrance.

Figure 7 .
Figure 7. Plasma energy flux Q = Qe + Q i , where Qe and Q i are the electron and ion energy flux respectively.Energy source for either species Qs.Qe blue line closed circles.Q i black line open circles.Qs dashed black line.Vertical dashed line source boundary, vertical dash-dot line sheath entrance.

Figure 8 .
Figure 8. Profile of the plasma energy fluxes, convective energy flux Qc, heat-flux q and total energy flux Q = Qc + q.Q black curve open triangles.q black curve closed triangles.Qc blue line squares.Vertical dashed line source boundary, vertical dash-dot line sheath entrance.

Figure 9 .
Figure 9. Profiles of the ion energy flux Q i , sum of the convective energy flux Q ci and of the heat-flux q i , compared to the energy source on the ion channel Qs.Q i black curve open head-up triangles.Q ci blue curve open squares.q i black line closed head-down triangles.Qs black dash-dot curve.Vertical dashed line source boundary, vertical dash-dot line sheath entrance.

Figure 10 .
Figure 10.Profiles of the electron energy flux Qe, blue line open circles, of the convective contribution Qce, blue line closed squares, and of the heat flux q i , black line closed triangles.

Figure 11 .
Figure 11.Profile of the electron thermal energy Te of the equipartition thermal energy taking into account the non vanishing heat-flux Teq.Blue line closed circles Te.Black line open triangles Teq, dashed horizontal black lines, sheath value of Te leading to the sheath transmission factor γ ≈ 6.8.

Figure 12 .
Figure 12.Profile of the ion thermal energy T i , of the equipartition thermal energy equipartition thermal energy taking into account the non vanishing heat-flux Teq, and T iγ .T i blue line closed circles, Teq black line open triangles, T iγ black curve open diamonds.

Figure 13 .
Figure13.Change of the ion kinetic energy K when varying the electric potential ϕ 0 − ϕ difference, where ϕ 0 is a maximum at the symmetry point x = 0. Given the sign of the velocity sign(v), ions evolve to the right for positive velocities and to the left for negative velocities according to the blue arrows.The distribution indicated by the shaded region splits accordingly, generating in particular a gap in the neighborhood of the axis K = 0.

Figure 14 .
Figure 14.Change of the ion kinetic energy K with the variation electric potential difference ϕ ∞ − ϕ.Ions for positive velocities, plain blue line and dashed blue line accelerate toward the wall on the right hand side.When neglecting the possible occurrence of trapped ions, plain and dashed black lines, the ion distribution function at ϕ is identical to that at ϕ ∞ as sketched by the shaded region.In such a framework, there is a gap with no ions for K ⩽ Kg at ϕ and for K ⩽ Kg − (ϕ ∞ − ϕ) at ϕ ∞ .

Figure 15 .
Figure 15.Profile of |ρc|, blue curve in log-scale.The reference value of the small control parameter in the Poisson equation ε 2 D is indicated by the horizontal dash-dot line.One can observe that ρc shoots-up above ε 2 D in the vicinity of the wall.

Figure 16 .Figure 17 .
Figure 16.Profile of ρc, blue curve in log-scale.The reference value of the small control parameter in the Poisson equation ε 2 D is indicated by the horizontal dash-dot line.One can observe that ρc shoots-up above ε 2 D in the vicinity of the wall.

Figure 18 .
Figure 18.Change of the electron kinetic energy K when varying the electric potential ϕ ∞ − ϕ close the vicinity of the sheath horizon.Electrons with positive velocity are slowed down toward the wall.Those with energy K ⩽ ϕ ∞ − ϕ are trapped, and reflected with negative velocity before reaching the wall.The passing electrons with K > ϕ ∞ − ϕ are absorbed by the wall.For the electrons the existence of trapped particles does not allow one to identify the electron density at ϕ given the distribution function at ϕ ∞ .An indicator of this issue is the singularity of the integral equation (14b) as v → 0 since the distribution function has finite values at v = 0.

Figure 19 .
Figure 19.Derivatives with respect to ϕ of the electron density ne, blue curve closed circles, and ion density n i , black curve open circles, compared to the kinetic calculation dn K /dϕ, red curve open triangles.While the derivatives of ne and n i converge in the quasineutral region ϕ > ϕ sheath , vertical dash dot line, the derivatives of ne and n K intersect at ϕ K , vertical black dashed line.Note that the bottom scale for ϕ is reversed so that the wall still stands on the right hand, typically at ϕ − ϕ 0 = −0.2.

Figure 20 .
Figure 20.Electron phase space characteristics using Ksign(v) as velocity coordinate with K = 1 2 v 2 /v 2 T0a and ϕ 0 − ϕ as position coordinate.Trapped electrons are confined to the shaded region, bounded by the characteristics issued from K = Kt at ϕ = ϕ 0 .For a source localized at ϕ = ϕ 0 , all electrons generated with K < Kt cannot be lost unless collisional transfer from K < Kt to K > Kt is taken into account.

Figure 21 .
Figure 21.Fraction of the source population with a Maxwellian distribution generated in the loss region with K > Kt.For a typical value of the potential drop ϕ 0 − ϕw = Kt ≈ 3.5, one finds that less than 1% of the source in injected directly in the loss region.Particle conservation then requires that the collisional flux from K < Kt to K > Kt balances the source flux to the trapped region, ensuring source sink balance.

Figure 22 .
Figure 22.Phase space contour plot, plain blue curve open circles and open squares, of the electron distribution function log 10 ( f e ), coordinates ϕ − ϕ 0 for the position and Ke sign(v) for the velocity, with Ke < 1.The wall is located on the right hand side in the direction of decreasing ϕ − ϕ 0 .The black lines with closed circles are the constant energy characteristics Ke = K 0 + ϕ − ϕ 0 plotted for different values of the kinetic energy K e0 at x = 0.

Figure 23 .
Figure 23.Phase space contour plot, plain blue curve open circles and open squares, of the electron distribution function log 10 ( f e ), coordinates ϕ − ϕ 0 for the position and Ke sign(v) for the velocity, with Ke sign(v) < −1.The wall is located on the right hand side in the direction of decreasing ϕ − ϕ 0 .The black lines with closed circles are the constant energy characteristics Ke = K 0 + ϕ − ϕ 0 plotted for different values of the kinetic energy K e0 at x = 0.

Figure 24 .
Figure 24.Ion distribution function at x = 0, plain blue curve closed circles.It departs from the Maxwellian with identical fluid moments f M i , thin blue curve open circles, but is well approximated by f fit , black curve open diamonds, the sum of a cold Maxwellian, black dashed line open triangles, and a hot Maxwellian, black dashed line open squares.

Figure 25 .
Figure 25.Phase space contour plot, plain blue curve open circles and open squares, of the ion distribution function log 10 ( f i ), coordinates ϕ − ϕ 0 for the position and K i sign(v) for the velocity, with K i < 0.2.The wall is located on the right hand side in the direction of decreasing ϕ − ϕ 0 .The black lines with closed circles are the constant energy characteristics K i = K 0 − (ϕ − ϕ 0 ) plotted for different values of the kinetic energy K i0 at x = 0.

Figure 26 .
Figure 26.Phase space contour plot, plain blue curve open circles and open squares, of the ion distribution function log 10 ( f i ), coordinates ϕ − ϕ 0 for the position and K i sign(v) for the velocity, with K i sign(v) > 0.2.The wall is located on the right hand side in the direction of decreasing ϕ − ϕ 0 .The black lines with closed circles are the constant energy characteristics K i = K 0 − (ϕ − ϕ 0 ) plotted for different values of the kinetic energy K i0 at x = 0.

Table 2 .
Reference locations in the simulation domain.
and on sheath constraints, γT e,div /T div M div .Using results from the analytical Navier-Stokes model presented in section 3.1, one can then show that Π ∝ Π ⋆ /|A div | with the isothermal closure and Π ∝ Π ⋆ /|A s,div | with the null heat-flux closure, where |A div | = 2|M div |/(1 + M 2 div