Yann Munschy
email: yann.munschy@cea.fr

Emily Bourne

Guilhem Dif-Pradalier

Peter Donnel

Philippe Ghendrih

Virginie Grandgirard

Yanick Sarazin

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.

I. INTRODUCTION

The prompt plasma recombination when in contact with a cold, dense and electrically neutral media drives a plasma flow towards 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 [START_REF] Loarte | Tsitrone, and the ITPA Scrape-off Layer and Divertor Group[END_REF]. The region where this physics take place has been named Scrape-Off Layer (SOL). Field lines in the Scrape-Off Layer 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 [START_REF] Loarte | Tsitrone, and the ITPA Scrape-off Layer and Divertor Group[END_REF]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 de-manding 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 cross-field 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 [START_REF] Bohm | The Characteristics of Electrical Discharges in Magnetic Fields (ch 3[END_REF]. In a fluid approach one can then use this boundary condition at the sheath entrance [START_REF] Bufferand | [END_REF]. 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 [START_REF] Joncquieres | A 10-moment fluid numerical solver of plasma with sheaths in a Hall Effect Thruster[END_REF]. 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 [START_REF] Boyd | Proceedings of the Physical Society. Section B[END_REF][START_REF] Persson | [END_REF]. 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][START_REF] Coulette | Plasma Physics and Controlled Fusion[END_REF][START_REF] Moritz | [END_REF][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 plasma-wall interaction, can be found in Ref. [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 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 [START_REF] Harrison | [END_REF] 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 Ref. ([11]).

We thus revisit these various approaches of the selforganized SOL/sheath/wall physics. We first address the SOL-wall interplay model in Section II: presenting some results from a reference kinetic simulation that motivate this effort (Section II A), the underlying kinetic equations (Section II B), highlighting the source terms (Section II B 2), and collision operators (Section II B 3). We then step to analyzing the SOL/sheath self-organization in the fluid framework Section III. 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 A, 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 d) and that of higher moment closures (Appendix A 2 e). 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 III A and III B. We discuss means of completing the fluid model in Appendix A 3 considering noncollisional closures (Appendix A 3 a), 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 A as well as the comparison between the fluid predic-tions and kinetic simulations are summarized in Section III C. Kinetic issues are discussed in Section IV. We revisit the previously published kinetic constraint in Section IV B and show that collisions are mandatory in the self-organization process (Section IV C). The steady state description of the distribution functions is presented in Section IV D and the role the source and collisional terms investigated. A Summary, discussion and conclusion Section closes the paper, Section V.

II. SOURCE TO SINK PLASMA ACCELERATION A. 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 [START_REF]ITER Physics Expert Group on Confinement and Transport, ITER Physics Expert Group on Confinement Modelling and Database, and ITER Physics Basis Editors[END_REF]. 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 Scrape-Off Layer 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 [START_REF] Stangeby | The Plasma Boundary of Magnetic Fusion Devices[END_REF]. 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 selforganization 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 Scrape-Off Layer 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 [START_REF] Gunn | [END_REF]. 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 FIG. (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 FIG. (1). In steady state, the particle source must be compensated by the build-up of a parallel particle flux Γ, FIG. (2). Because of the symmetry of the simulation set-up, only the right hand side of the simula- tion 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,[START_REF] Stangeby | The Plasma Boundary of Magnetic Fusion Devices[END_REF]. Accordingly, one expects that when the source term vanishes in such a quasineutral region, the particle flux becomes constant and plasma acceleration stops [4,[START_REF] Stangeby | The Plasma Boundary of Magnetic Fusion Devices[END_REF]. 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 ui = |du i /dx|/u i , L -1 ni = |dn i /dx|/n i .
Here n i stands for the ion density, hence the subscript i, FIG. (3). The length scale L ui is the blue curve with open triangles. It increases in the source region, rolling over to a range of values of order L ∥ . Outside the source region, L ui levels-off, then, left of the vertical dash-dot line, L ui exhibits a strong decrease towards values of the order of 10 Debye lengths. The length scale L ni for the ion density exhibits the same behavior as L ui 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 ui → 0 while L ni → +∞. In a narrow region close to the wall, on the right hand 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,[START_REF] Stangeby | The Plasma Boundary of Magnetic Fusion Devices[END_REF], 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 Γ, FIG. (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 sug- gests 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. In Section II A 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 Eq. (1a) that governs the evolution of the distribution function f a is then:

∂ t f a + v∂ s f a - e a m a ∂ s U ∂ v f a = C(f a ) + S(f a ) (1a)
where C(f a ) is the collision operator, standing for both self-collisions 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 ρ c .

-∂ 2 s U = ρ c /ε 0 ; ρ c = species a e a n a (1b)
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:

∂ t f a + v∂ s f a - T 0 m a ∂ s ϕ ∂ v f a = C(f a ) + S(f a ) (2a) ε 2 D ∂ 2 x ϕ = species a e a n a en 0 ; ε 2 D = λ 2 D0 L 2 0 = ε 0 T 0 n 0 e 2 L 2 0 (2b)
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 FIG.

(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 Eq.(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 Eq.(2a) is similar to the Gysela source terms [START_REF] Sarazin | Angelino[END_REF]. 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,[START_REF] Dif-Pradalier | 19th International Conference on Plasma-Surface Interactions in Controlled Fusion[END_REF], 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 [START_REF] Valade | [END_REF]. 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 [START_REF] Baschetti | Radiation Driven Bifurcations in Fusion Plasmas[END_REF]. This particular regime has implications when addressing the plasma collisionality, see Section II B 3. The particular self-organization 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 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 Eq.(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.

S(f

a ) = M s (x) s k 2πT s /m a exp - m a v 2 a 2T s (3)
The mask M s , plotted in black on FIG.

(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 SOL-plasma 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 [START_REF] Sarazin | [END_REF]. 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 [START_REF] Ghendrih | Proc. 20th Int. Conf. Vilamoura[END_REF]. 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 Eq. (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 IV C. 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 Eq. ( 4), similar to the standard ν ⋆ parameter [START_REF] Doyle | ITPA Transport Physics Topical Group, ITPA Confinement Database and Modelling Topical Group, and ITPA Pedestal and Edge Topical Group[END_REF] 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 [START_REF] Loarte | Tsitrone, and the ITPA Scrape-off Layer and Divertor Group[END_REF]. 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 [START_REF] Dif-Pradalier | 19th International Conference on Plasma-Surface Interactions in Controlled Fusion[END_REF] 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 the Landau and Lenard-Balescu non-linear operators [START_REF] Nicholson | Introduction to plasma theory[END_REF]. The inter-species collisions are set to drive the distribution functions towards 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

ν ⋆ 0 = ν ⋆ D0 L ∥ /λ D0 . ν ⋆ 0 = 4 √ π 3 e 2 4πε 0 2 n 0 T 2 0 log Λ L ∥ = ν ⋆ D0 L ∥ λ D0 = L ∥ L coll (4)
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.

III. 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 II A, 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.

A. Fluid properties in steady state plasma conditions

In steady state conditions discussed in Section A 1 d, one can combine the electron and ion conservation equations to obtain plasma balance equations for the number of particles, plasma momentum and plasma energy, Eq.(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 Eq.(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, Eq.( A6c) is simplified so that the steady state plasma balance equations take the form:

d s Γ = S n (5a) d s Π = 0 (5b) d s Q -T s Γ = 0 (5c)
As discussed in the Appendix, the kinetic simulation is found to exhibit the three conserved fluxes Eq.( 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 A 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:

Γ = s 0 ds S n ; Γ w = L ∥ 0 ds S n (6a) Π qn = constant (6b) Γ 3T + m i u 2 -2T s = 0 (6c)
Given the velocity V defined by

m i V 2 = T and M V = u/V , one can write Γ = n i M V V , Π qn = n i T (1 + M 2 V ) and Q c = ΓT (3 + M 2 V )
. According to Eq. ( 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 Eq.(6c), one can determine T in terms of M V .

T = 2T s 3 + M 2 V (7a)
Given the expression of Π qn , one finds p = Π qn /(1 -

+ M 2 V ), then dp/p = 2M V dM V /(1 + M 2 V ) 2 ,
4M V dM V /(3 + M 2 V )/(1 + M 2 V ).
Since the polytropic index γ p is defined by dp/p = γ p dn i /n i , one then finds:

γ p = 3 + M 2 V 2 (7b)
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 Eq.(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

Π qn by Γ/(M V V ), so that Π qn = T /V Γ(M V + 1/M V ).
One can then set Π qn V /(ΓT ) = 2/A, and the change of variable 2/A = M V + 1/M V [4]. However, T /V = m i V also depends on M V as found in Eq. (7a). It is then more relevant for this non-isothermal case to consider a further change of variable.

A 2 s = 16Γ 2 m i T s 9Π 2 qn = 8M 2 1 + M 2 1 + 3M 2 2 (7c)
One then finds that M increases from M = 0 at the stagnation point with Γ = 0 and therefore A 2 s = 0 to M = 1 for A 2 s = 1 where A 2 s is maximum. For M > 1, A 2 s decreases and at M → +∞, A 2 s → 8/9. In Eq. (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 Eq. (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 Eq. (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 Section A 4 is similar to that used in Ref. [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 A s = 1 when Γ = Γ w .

Π qn = 4 3 m i T s Γ w (7d)
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 Eq. (7c), and T = (T s /3) 2/(1 + M 2 ), Eq. (7a). At constant momentum flux Π qn , set by Eq. (7d), the density n i decreases from the value at the stagnation point

M 2 = 0, n M 0 = Π qn /(2(T s /3)) = n M 0 = 2 m i /T s Γ w , to n M 1 = n M 0 /2 at M 2 = 1, since n = n i = n M 0 (1 + M 2 )/(1 + 3M 2 ).
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 Eq.(7a):

T e = T i = (T s /3)/(1 + M 2 ) (8)
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:

d s Γ( 3 2 T i + 3 2 T M 2 -1 2 T s ) = 0 (9a) d s Γ( 3 2 T e -1 2 T s ) = 0 (9b)
This yields:

T e = (T s /3) (10a)
T i = (T s /3) 1 -M 2 1 + M 2 (10b)
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.

T e = (T s /3) 3 γ (11a)
Then setting γi = γQ ci /Q, one obtains:

T i = (T s /3) 2γ i -3M 2 γ 1 + M 2
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:

T i = (T s /3) γ -3M 2 γ 1 + M 2 (11b)
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.

B. Fluid interpretation of the kinetic simulation

The control parameters of the reference simulation used as evidence are listed in Table I. 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 IV B, Table III -criterion 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 II. 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 FIG.(4), but comparing u i to c eff = 1.505 V 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. The latter is a smooth transition out of quasineutrality conditions. The choice of c eff = 1.505 V would have to be slightly adapted to changes in the sheath entrance definition. The variation of A s , Eq.(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 Eq.(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 I. 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 , FIG. (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. Towards 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 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, FIG. [START_REF] Bufferand | [END_REF], 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 towards 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, FIG. [START_REF] Joncquieres | A 10-moment fluid numerical solver of plasma with sheaths in a Hall Effect Thruster[END_REF], 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 head-up 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. Towards 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 heat-flux 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 towards the wall and contribute to the energy exhaust. This point will be further discussed when addressing non-collisional closure in Section A 3. With such a closure, one finds that the heat-flux 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 heat-flux 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 towards T e /T s ≈ 0.148 at the sheath entrance, see FIG. [START_REF] Harrison | [END_REF]. Blue curve closed circles, T e . The equipartition value T eq is given in Eq.( 8) but taking into account that Q c /Q < 1 due to the non-vanishing contribution of the heat-flux,

T eq = T s /(3(1+M 2 ))Q c /Q.
Black curve 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 by a factor 1.7 at the wall boundary. Given Eq.(11a), and the value of T e /T s = 0.148 at the sheath entrance, one can determine the sheath transmission factor γ ≈ 6.8. Sheath constrained electron thermal energy T e /T s = 0.148, dashed horizontal line on FIG. [START_REF] Harrison | [END_REF]. One finds that the equipartition value overestimates T e while the sheath constraint underestimates T e . 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, FIG. (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 Eq.( 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 al- though 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 Eq.(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.

C. 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 FIG.

(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 Γ/Γ 0 = 0.992 ± 2. 10 -3 , Π/Π 0 = 1.1814 ± 3. 10 -4 , Q/Q 0 = 1.00075 ± 3.5 10 -4 , where

V 2 0 = T 0 /m i , Γ 0 = n 0 V 0 , Π 0 = n 0 T 0 and Q 0 = n 0 T 0 V 0 .
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 e, non-collisional closure in Appendix A 3 a. 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 c, A 2 d, A 2 e and A 3 b. 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 heat-flux 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 selforganization 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 towards the wall.

IV. KINETIC PROCESSES IN THE SOL-SHEATH SELF-ORGANIZATION

We now analyze the kinetic constraint that determines the sheath entrance. It was first published in ref. [START_REF] Harrison | [END_REF], then in later papers [13]. We first present the original derivation and then adapt it to the actual conditions met in the simulations.

2 ε 2 D dϕ dx 2 -1 2 ε 2 D dϕ dx 2 ∞ = G(ϕ) (12b)
G(ϕ) = - ϕ ϕ ∞ dϕ ′ ρ c (12c)
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:

1 2 ε 2 D dϕ dx 2 = G(ϕ ∞ ) + G ′ (ϕ ∞ ) ϕ -ϕ ∞ + 1 2 G ′′ (ϕ ∞ ) ϕ -ϕ ∞ 2 + . . . (13a) By definition G(ϕ ∞ ) = 0 and G ′ (ϕ ∞ ) = -ρ ∞ .
Since quasineutrality holds at the sheath horizon ρ ∞ → 0, one then finds that:

1 2 ε 2 D dϕ dx 2 ≈ 1 2 G ′′ (ϕ ∞ ) ϕ -ϕ ∞ 2 (13b)
The sheath horizon constraint is therefore that G ′′ (ϕ ∞ ) ≥ 0, so that:

dn e dϕ ϕ ∞ - dn i dϕ ϕ ∞ ≥ 0 (14a)
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 ref. [START_REF] Harrison | [END_REF] 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 H = K + ϕ where

K = 1 2 v 2 .
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 f ∞ = f i (K ∞ , ϕ ∞ ).
Let us consider a source distribution function initially generated at x = 0, the SOL symmetry point where the electric potential is maximum, FIG. (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 towards 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 FIG. 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 towards 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 ϕ ∞ .
potential and cannot reach the point x = 0 along the characteristics, an example is given by the dashed black lines on FIG. (13).

Let us now consider the neighborhood of the sheath horizon with electric potential ϕ ∞ and given a decreasing electric potential towards the right hand side where the wall is localized, dashed region on the right hand side. We first only consider positive velocities since no ions are assumed to flow back from the wall. From ϕ ∞ , black dashdot 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

K g ≥ ϕ ∞ -ϕ. n i (ϕ) = +∞ vg dv f i (v, ϕ) = +∞ 0 dK f i (K, ϕ) v(K, ϕ)
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

f i (K, ϕ) = f i (K ∞ , ϕ ∞ ) provided K +ϕ = K ∞ +ϕ ∞ . One then obtains the identity f i (K, ϕ)dK = f i (K ∞ , ϕ ∞ )dK ∞ .
The density n i (ϕ) can then be determined by the values of the distribution function at the sheath horizon f i (K ∞ , ϕ ∞ ).

n i (ϕ) = +∞ 0 dK ∞ f i (K ∞ , ϕ ∞ ) v(K, ϕ)
One can then use v(K, ϕ) = 2(K ∞ + ϕ ∞ -ϕ) to determine dn i /dϕ:

dn i dϕ ϕ = +∞ 0 dK ∞ f i (K ∞ , ϕ ∞ ) v(K, ϕ) 3
Then setting ϕ → ϕ ∞ one finally obtains:

dn i dϕ ϕ ∞ = +∞ 0 dv ∞ f i (K ∞ , ϕ ∞ ) v 2 ∞ ( 14b
)
We have therefore recovered the expression of Ref. ([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 f i (v ∞ )/v 2 ∞ → 0 as v ∞ → 0.
As will be shown in the following Section IV B, this proves to restrict the use of this elegant result.

B. 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 Eq. (A5). In 1-D, 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 16) 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 towards the fourth region, in the vicinity of the wall, where ρ c ≫ ε 2 D and ρ c exhibits a sharp exponential growth by 4 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 Eq.(13a) does not hold. One must then consider the complete expression of the electric field vari- one can read in this table the wall position, x w and the transition points for ϕ, E, ρ c , |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 param-TABLE III. 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 Eq.(15b). The position x sheath is defined as the transition point for dρc/dϕ, x sheath = 195.9. eters 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 Eq. (14a), using Eq. (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 Eq. ( 14b)

n K (ϕ) = +∞ 0 dv f i (K, ϕ) ; dn K dϕ = +∞ 0 dv f i (K, ϕ) v 2 (15a) R(ϕ) = dn e dϕ ϕ - dn K dϕ ϕ (15b)
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 Eq. ( 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 Eq.(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 FIG. 15, the slowing down of the electrons towards 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 Eq. (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 Eq. (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 Eq. ( 15) are difficult to use. 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 Eq. ( 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, towards 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ϕ, FIG. (18).

C. 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 to 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, towards 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 FIG. (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 FIG. 20. Electron phase space characteristics using Ksign(v) as velocity coordinate with K = 1 2 v 2 /v 2 T 0a 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.

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 FIG. [START_REF] Coulette | Plasma Physics and Controlled Fusion[END_REF], 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.

D. 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 FIG. (23), 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 point x = 0, ϕ = ϕ 0 towards the wall located on the right hand side. Conversely, particles with negative velocities flow towards x = 0 and beyond towards the wall on the left hand side. The 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.09 and 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, FIG. (24), 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, FIG. (24). 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 FIG. [START_REF] Moritz | [END_REF]. 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 IV C, 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 ), FIG. (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 towards 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, FIG. (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 left-hand 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.

V. SUMMARY, DISCUSSION AND CONCLUSION

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 Scrape-Off Layer plasma of magnetic confinement devices dedicated to fusion. We have also restricted the analysis to a 1-D 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 II A. 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 a 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 cross-field 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 II B. 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 heat-fluxes. 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 Eq.(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 III A. 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, Section 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 Section 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 1-D 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 Section 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 and the transport equations 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, Section A 3. The analysis of this projection in Fourier space indicates that the heatflux 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, Section 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, Section A 2 e. 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. How-ever, 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 [START_REF] Harrison | [END_REF]. 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 IV A 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 towards 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 towards the wall. The Harrison-Thompson criterion must therefore be adapted to actual plasma behavior as displayed in the kinetic simulation, see Section IV B. 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 IV C. 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 towards 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 towards 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 1-D 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 SOL-Sheath 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 2-point model [START_REF] Stangeby | The Plasma Boundary of Magnetic Fusion Devices[END_REF], the plasma total pressure is linked to the upstream particle and energy flux and to sheath constraints [START_REF] Baschetti | Radiation Driven Bifurcations in Fusion Plasmas[END_REF][START_REF] Capes | [END_REF]. 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 veloc-ity V div . Note that M div is not the Mach number but corresponds to M V (defined and used in Section III A) 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:

Π 2 = m i γT e,div /T div Γ up Q up M div + 1 M div 2 (16)
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 Π ⋆ = m i Γ up Q up , and on sheath constraints, γT e,div /T div and M div . Using results from the analytical Navier-Stokes model presented in Section III A, 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 ) and A s,div is defined in Eq.(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 self-organization 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 towards 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 heat-flux, 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. The Navier-Stokes fluid moments can be computed using the dimensional form of the Boltzmann equation Eq.(1a). For species a, one readily obtains:

∂ t n a + ∇ s Γ a = S a n (A1a) m a ∂ t Γ a + ∇ s Π a = -e a n a ∇ s U + C a mΓ (A1b) ∂ t E a + ∇ s Q a = -e a Γ a ∇ s U + C a E + S a E (A1c)
where C a mΓ and C a E stand for the momentum and energy collisional exchange and S a n , S a E for the particle and energy sources. With the standard definitions of the density n a , particle flux Γ a , momentum flux Π a , energy density E a , and energy flux Q a for species a:

n a = +∞ -∞ dv f a (A2a) Γ a = +∞ -∞ dv v f a = n a u a (A2b) Π a = +∞ -∞ dv m a v 2 f a = n a T a + m a n a u 2 a = 2E a (A2c) Q a = +∞ -∞ dv 1 2 m a v 2 v f a = q a + u a E a + p a (A2d)
These definitions also introduce the mean particle velocity u a , the thermal pressure p a = n a T a , the thermal energy T a , and the heat flux q a . The latter is the difference between the energy flux Q a and the convective energy flux Q ca = u a E a + p a . In the 1-D limit, the momentum flux Π a and energy density E a are related Π a = 2E a . One can thus remark that the Navier Stokes system Eq.(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.

b. 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:

∂ t n i + ∇ s Γ i = S n (A3a) ∂ t ρ c + ∇ s j = 0 (A3b)
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.

∂ t (m e Γ e ) + ∇ s Π e = en e ∇ s U + C e mΓ (A4a) ∂ t (m i Γ i + m e Γ e ) + ∇ s (Π i + Π e ) = -ρ c ∇ s U (A4b) ∂ t E e + ∇ s Q e = eΓ e ∇ s U + C e E + S E (A4c) ∂ t (E e + E i ) + ∇ s (Q e + Q i ) = -j∇ s U + 2S E (A4d)
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 Eq.(A4d).

c. 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 Eq.(2b) in the asymptotic limit ε 2 D → 0. The charge balance equation Eq.(A3b) then enforces that ∇ s j = 0, hence a 1-D 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:

∂ t n + ∇ s Γ = S n (A5a) m i ∂ t Γ + ∇ s Π = -ρ c ∇ s U (A5b) ∂ t E + ∇ s Q = 2S E (A5c)
where Π = Π e + Π i = p e + p i + m i n i u 2 i , E = E e + E i = Π/2, and Q = Q e + Q i . 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 Eq.(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 Eq.(A5).

d. Plasma balance equations in steady state

Let us now consider the steady state regime with respect to the ion evolution time. The various fields now only vary in space and one can replace the space derivative ∇ s by d s .

d s Γ = S n (A6a) d s Π = -ρ c d s U (A6b) d s Q = 2S E (A6c)
Since the energy source is convective, proportional to the particle source S n , and T s /2 energy per injected particle, this model yields 2S E = T s S n . Using this relation to modify Eq.(A6c) and the quasineutral regime ρ c → 0 to modify Eq.(A6b), the system takes the form:

d s Γ = S n (A7a) d s Π = 0 (A7b) d s Q -T s Γ = 0 (A7c)
Since both Γ and Q are anti-symmetric with respect to the mid-box s = 0 while Π is symmetric, one obtains:

Γ = s 0 ds S n (A8a) Π = constant (A8b) Q = ΓT s (A8c)
These equations define three invariants that describe the quasineutral plasma behavior.

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. 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

M (a) ℓ (s, t) = dv v ℓ f a (s, v, t) (A9a) M (a) 0 = n a ; M (a) 1 = M (a) 0 u a = Γ a (A9b)
Then defining w = v -u a one defines:

N (a) ℓ (s, t) = dv w ℓ f a (s, v, t) (A9c) N (a) 0 = M (a) 0 ; N (a) 1 = 0 (A9d)
The moment M ℓ is related to the lower moments N j , j ≤ ℓ, and the moment N ℓ to the lower moments M j ,

j ≤ ℓ. M (a) ℓ = dv (w + u a ) ℓ f a = ℓ j=0 C j ℓ (u a ) j N (a) ℓ-j (A10a) N (a) ℓ = dv (v -u a ) ℓ f a = ℓ j=0 C j ℓ (-u a ) j M (a) ℓ-j (A10b) C j ℓ = ℓ! (ℓ -j)!j! (A10c)
Let us consider the fluid hierarchy for species a:

∂ t M (a) ℓ + ∂ s M (a) ℓ+1 = X (a) ℓ + C (a) ℓ + S (a) ℓ (A11a) X (a) ℓ = e a m a ∂ s U ℓM (a) ℓ-1 (A11b) C (a) ℓ = dv v ℓ C a (f a ) ; S (a) ℓ = dv v ℓ S a (f a ) (A11c)
The mass center combination of the moments m i M

ℓ + m e M (e) ℓ = (m i + m e )M ℓ then yields:

∂ t M ℓ + ∂ s M ℓ+1 = e∂ s U ℓ M (i) ℓ-1 -M (e) ℓ-1 m i + m e + m i C (i) ℓ + m e C (e) ℓ m i + m e + m i S (i) ℓ + m e S (e) ℓ m i + m e (A12a)
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 Eq.(A12a) leads to Eq.(A5a).

∂ t M 0 + ∂ s M 1 = S n (A12b)
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. Eq.(A12a) therefore leads to:

∂ t M 1 + ∂ s M 2 = 0 (A12c)
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:

∂ t M 2 + ∂ s M 3 = 4 m i + m e S E (A12d)
For the higher moments, the left hand side of the moment equations keeps the same structure as found in Eq.(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:

∂ t n + ∇ s Γ = S n (A13a) ∂ t Γ + ∇ s Π/m i = 0 (A13b) ∂ t Π/m i + ∇ s 2Q/m i = 4S E /m i (A13c)
Let us now consider fluctuations of a steady state solution. We assume that these fluctuations are not driven by fluctuations of the source term, these are chosen fixed.

The fluctuation amplitudes are assumed to be small, allowing us to linearize the evolution equations, and use Fourier modes as convenient basis to determine the eigenmodes. The phase of the Fourier modes is chosen to be -ωt+ks, so that the phase velocity of interest is c = ω/k. The fluctuating part of any field F is noted F and its Fourier amplitude F . In the quasineutral limit, the system Eq. ( A13) is closed but for the coupling to higher moments governed by the heat flux q contributing to the total energy flux Q. To compute the possible values of c, we close the fluid hierarchy by assuming the linearized energy flux to be a linear combination of the lower moments, density n, particle flux Γ and momentum flux Π.

We then find that non-trivial solutions are obtained when the following determinant is null.

-

c 1 0 0 -c 1 a n a Γ a π -c = 0 (A14a)
Where the coefficients are define as: a n = ∂ n (2Q/m i ) at fixed Γ and Π, a Γ = ∂ Γ (2Q/m i ) at fixed n and Π and a π = ∂ Π (2Q) at fixed n and Γ. The dispersion equation is therefore cubic in c:

c 3 -a π c 2 -a Γ c -a n = 0 (A14b)
Finding the roots of this dispersion equation is eased when the chosen closure allows a class of density fluctuations to exhibit a phase velocity c equal to the mean velocity u.

b. Density fluctuations propagating at mean plasma velocity

We address here a particular behavior of the system that stems from the left hand side structure of the moment equation Eq.(A12a). Removing the drive of the fluctuations by the source terms, setting the electric field E = -∂ s U to zero and ignoring the contribution of the collisions. These assumptions ensure that the plasma moment equations are similar to that of a neutral fluid at local thermodynamic equilibrium, only the left hand side is retained.

∂ t M ℓ + ∂ s M ℓ+1 = 0 (A15a)
When stepping to the dispersion relation this linear equation readily leads to:

M ℓ+1 = c M ℓ (A15b)
the fluid prediction to the evidence of the kinetic simulation that this simplifying assumption does not hold in particular for the electrons. This shortfall emphasizes the kinetic nature of the self-organization process. However, this does not mean that a fluid description cannot address the heat-flux physics since one can increase the number of moments used in the fluid description, or address a better closure of the Navier-Stokes 3-moment fluid model by allowing the heat to be finite. In many models of plasma-wall interaction a collisional closure for the heat flux is retained, so that for a species a q a = -κ a ∂ s T a .

One finds that κ a does not depend on density, exhibits a non-linear dependence on the thermal energy T 5/2 a , and is typically proportional to the inverse of the mass square root 1/ √ m a . Such a local formulation assumes that L coll ∂ s T a /T a ≪ 1 where L coll is the collision mean free path. When addressing plasma-wall interaction this assumption is an issue so that non-local features are being addressed [54,55]. Regarding the comparison to the low collisionality reference simulation, ν ⋆ 0 = 0.1, see Table I, the collision mean free path is typically 10L ∥ so that the random walk conditions underlying the collisional closure are not fulfilled. Furthermore, we have observed that the electron heat-flux q e builds-up in the source region at nearly constant, or even slightly positive temperature gradient. It then appears more interesting to address a non-collisional closure such that the heat-flux q is a linear combination of the lower moments [56][57][58]. The starting point are the linearized time dependent fluid equations.

Setting Ω = ω/(kV e ) with m e V 2 e = T 0 , p 0 = n 0 T 0 , and linearizing n U ≈ n 0 U and Γ U ≈ Γ 0 U one then obtains in Fourier space:

Ω n n 0 = Γ n 0 V e (A37a) Ω Γ n 0 V e = p e p 0 - e U T 0 (A37b) Ω p e p 0 = 2 Q e p 0 V e - 2Γ 0 n 0 V e e U T 0 (A37c)
As done when determining the sound velocity in Section A 2 c, 2 Q e /(p 0 V e ) is then projected on the lower moments n, Γ and Π e ≈ p e , therefore:

2 Q e p 0 V e = a n n n 0 + a Γ Γ n 0 V e + a π p e p 0 (A38a)
So that Eq. (A37c) becomes

Ω -a π p e n 0 T 0 = a n + Ωa Γ n n 0 - 2Γ 0 n 0 V e ϕ (A38b)
and therefore Eq. (A37b) becomes n n 0 = ϕ Ω -a π + 2Γ 0 /(n 0 V e ) a n + Ωa Γ + Ω 2 a π -Ω 3 (A38c)

In Section A 2 c, the projection coefficients a n , a Γ and a π were determined by approximation Q e ≈ Q ce . In the proposed non-collisional closure we determine these coefficient by matching the relation between the electron density response to fluctuations of the electric potential of the fluid and kinetic frameworks. In that case, using the Poisson equation, one enforces the same dispersion relation for Langmuir waves, recovering Landau damping in the fluid framework. However, 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:

∂ t f e + v∂ s f e + e m e ∂ s U ∂ v f e = 0 (A39a) so that one obtains for the Fourier modes :

v - ω k f e = - T 0 m e ϕ ∂ v f e (A39b)
The equilibrium contribution is chosen to be a shifted Maxwellian:

f e = n 0 √ 2πV e exp(-w 2 )
where w 2 = (v -u) 2 /(2V 2 e ) and such that ∂ v f e = √ 2wf e /V e , therefore:

f e = 1 + y w -y ϕ f e (A40a) y = |ω/k -u| √ 2V e = √ 2 2 |Ω -ε u | (A40b)
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:

n n 0 = 1 √ π +∞ -∞ dw 1 + y w -y e -w 2 ϕ (A40c)
The Fried and Conte function Z is defined by:

Z(y) = 1 √ π +∞ -∞
dw e -w 2 w -y one obtains:

n n 0 = 1 + yZ(y) ϕ (A40d)
Eq.(A58c). Using the control parameter A s Eq. (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 Eq. (A51a) and given Π = constant the sign of dΠ qn /dϕ is identical to that of dρ c /dϕ.

dA s A s dϕ = - dΠ qn Π qn dϕ = dM M dϕ 1 -M 2 (1 + M 2 )(1 + 3M 2 ) dM M dϕ = - 1 1 + 2M 2 dn ndϕ < 0 dΠ qn Π qn dϕ = dρ c ρ 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.

Discontinuity of the steady-state solution

In this Section we use the formalism of Section A 2 e to investigate the transition into the sheath, and examine the criterion based on the discontinuity of the steady state solution. The starting point is therefore the system Eq.(A27) for steady state conditions and restricting the right hand side term Z ℓ to the source term as in Eq.(A28a) We therefore consider Z ℓ = S ℓ-1 completed by Eq. (A10b).

d s M ℓ = S ℓ-1 (A59a) N ℓ = ℓ j=0 C j ℓ (-u) j M ℓ-j (A59b)
We differentiate Eq.(A59b) and therefore:

d s N ℓ = ℓ j=0 C j ℓ (-u) j d s M ℓ-j -d s u ℓ j=1
jC j ℓ (-u) j-1 M ℓ-j

One can then replace d s M ℓ-j by the source S ℓ-j-1 according to Eq.(A59a), with some caution when handling S -1 that we set by construction identical to d s M 0 , S -1 = d s M 0 . One take into account that M 0 d s u + ud s M 0 = d s M 1 = S 0 and, recall the relation Eq.(A30) obtained in Section A 2 e: ℓ j=1 jC j ℓ (-u) j-1 M ℓ-j = ℓN ℓ-1

One then relates the derivatives of M 0 and N ℓ to the source terms:

-u ℓN ℓ-1 M 0 + (-u) ℓ d s M 0 + d s N ℓ = ℓ-2 j=0
C j ℓ (-u) j S ℓ-j-1 + ℓ (-u) ℓ-1 -

N ℓ-1 M 0 S 0 (A60)
Two types of closures have been addressed in Section A 2 e, for ℓ = 2, we have defined N 2 = V 2 M 0 with V 2 constant, while for ℓ ≥ 3 N ℓ = 0 has been considered. For either choice of the closure, one can then use Eq.(A60) to determine ud s M 0 . A singularity will appear if the coefficient of ud s M 0 on the left hand side of Eq.(A60) vanishes while the right hand side cannot be null. For ℓ = 2, we have d s N 2 = V 2 d s M 0 . Given N 1 = 0, this leads to:

V 2 -u 2 d s M 0 = S 1 -2uS 0 (A61a)
One finds that the derivative d s M 0 diverges for u 2 → V 2 whenever S 1 -2uS 0 ̸ = 0. At the divergence point u 2 = V 2 , we have shown that the phase velocity of the sound wave verifies (c -u) 2 = V 2 = c 2 s , Eq.(A33a), so that one finds (c-u) 2 = u 2 and therefore either c = 2u or c = 0. The latter condition leading to a standing shock at the divergence point. However, if S 1 -2uS 0 = 0 can be satisfied, then V 2 -u 2 = 0 where S 1 -2uS 0 = 0, so that there is no divergence and a smooth transition at u 2 = V 2 . When replacing Π with Π eq , Π = Π qn -ηT e Eq.(A51a), one modifies S 1 , S 1 ≈ T e d s η/m i . Here we have neglected the contribution proportional to the electric field E, eηE/m i , assuming that the condition T e d s η/m i ≫ eηE/m i previously verified when discussing the expansion at the sheath horizon, Section A 4. In the reference simulation the ratio between these two terms exceeds 100 in the region of interest, prior to the sheath entrance and nearly up to the wall. Close to the wall these two terms become comparable. A smooth transition at the Bohm condition u 2 = c 2 s occurs together with the constraint:

d s η = 2m i uS 0 T e (A61b)
If one lets S 0 small, one obtains a threshold in d s η ≥ 0. Since S 1 is assumed null for η = 0 in the quasineutral region, the transition at u 2 = V 2 , therefore the sheath entrance is characterized by the Bohm criterion M 2 = 1, together with d s η → 0 + . When relaxing the quasineutrality constraint, one removes the singularity at M 2 = 1 by enforcing the constraint Eq.(A61b) at the sheath entrance. Then both denominator and numerator change sign so that M 0 remains a monotonically decreasing function.

For ℓ = 3, one closes the hierarchy with N 3 = 0, so that Eq.(A60) with N 2 = V 2 M 0 reads:

-u d s M 0 3V 2 -u 2 = S 2 -3uS 1 + 3S 0 u 2 -V 2 = S 0 2T s m i + 3u 2 -3V 2 -3u T e m i d s η
When taking into account Eq.(6c) yielding 2T s /m i = 3V 2 + u 2 to obtain:

-u d s M 0 3V 2 -u 2 = 4u 2 S 0 -3u T e m i d s η (A62a)
One finds that there is no singularity for u = 0 and that the singularity for u 2 = 3V 2 is avoided provided:

d s η = 4m i uS 0 3T e (A62b)
This constraint is similar to that determined for ℓ = 2 and indicates that for S 0 → 0 the sheath constraint is set by d s η ≥ 0, the condition d s η = 0 being met when u 2 = 3V 2 , the Bohm condition M 2 = 1. The sheath entrance is therefore characterized by the continuity of the derivative d s M 0 with a value corresponding to that of the SOL plasma, M 0 monotonically decreasing on the characteristic length scales of the SOL, prior to a sharper decrease with variations on the Debye scale towards the wall.

Generalizing the result obtained for ℓ = 3 with the closure N ℓ defined to be null, therefore d s N ℓ = 0, one obtains d s M 0 :

-ud s M 0 = N ℓ D ℓ (A63a)

N ℓ = ℓ-1 j=0 C j ℓ (-u) j S ℓ-j-1 -ℓ N ℓ-1 M 0 S 0 (A63b) D ℓ = ℓN ℓ-1 M 0 -(-u) ℓ-1 (A63c)
One recovers the possible divergence when the denominator D ℓ is null. The latter condition occurs when (-u) ℓ-1 = c ℓ-1 s , therefore depending on the sound velocity, Eq.(A35). One recovers here sign issues for the even values of ℓ = 2ℓ ′ that can only be resolved if N 2ℓ ′ -1 ≤ 0. For the odd values of ℓ, D ℓ = 0 yields u = ±c s where c s verifies Eq.(A35b). For example, for ℓ = 4, D 4 = 0 leads to u 3 = -4N 3 /M 0 . The sign of the solution for u is the opposite of the sign of N 3 . This condition is not met in the reference kinetic simulation where the heat flux q ∝ N 3 exhibits the same sign as the mean velocity u. One can override such difficulties by prescribing that the closure N ℓ = 0 can only be considered for odd values of ℓ. The divergence of ud s M 0 governed by D ℓ = 0 can be removed by setting the numerator to zero N ℓ = 0. As discussed for S 1 , one must include possible contributions proportional to the electric field in the source terms S j with j > 2. One can still prescribe a condition on d s η, but the constraint will become more complicated.

The possible discontinuity of the derivative of the moment M 0 = n, and consequently of other fields such as the Mach number, is found to be closely related to the sound velocity and occurring via a vanishing denominator at Mach = 1. Analyzing the numerator, one finds that it vanishes when a constraint of the derivative of the charge density η = n i -n e is satisfied, namely d s η ∝ S 0 , S 0 being the particle source. A smooth transition is thus enforced where both denominator and numerator change sign. The plasma then switches from subsonic with d s η → 0 to supersonic when d s η becomes finite. We therefore recover the results obtained by asymptotic expansion in the vicinity of the sheath horizon Section A 4. In the asymptotic limit of the sheath, one assumes that S 0 → 0 so that the sheath criterion is d s η ≥ 0 and the sheath entrance the point where d s η departs from zero. However, the asymptotic matching in Ref. [4] is determined for a finite value of S 0 , and the sheath entrance at Mach = 1 is therefore shown to occur for d s η > 0. One must remark here that in Ref. [4], the asymptotic matching is performed by assuming that at constant thermal energy T the unique control parameter is A, therefore A s Eq.(7c) with varying thermal energy. This meaning that all quantities of interest, therefore both the Mach number and d s η only depend on A s . In the reference simulation, S 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 FIG. (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 FIG. (3).

FIG. 1 .

 1 FIG. 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.

FIG. 2 .

 2 FIG. 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.

FIG. 3 .

 3 FIG. 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.

B. Kinetic model to address open field line self-organization 1 .

 1 Kinetic equation and normalization

  FIG. 5. Mach number M versus As, Eq. (7c), maximum of As = 1 at M = 1.

FIG. 6 .

 6 FIG. 6. Profile of the mean plasma velocity ui ≈ 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. ui 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.

FIG. 7 .

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

FIG. 8 .

 8 FIG. 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.

FIG. 9 .

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

FIG. 10 .

 10 FIG. 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 qi, black line closed triangles.

FIG. 12 .

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

A 2 D d dx 2 ϕ

 22 . 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 Eq. (2b), slightly modified here to account for a possible change of the normalization length scale L. ε = -ρ c (12a) 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 Eq.(12a) by dϕ/dx and integrating over x.

1

 1

FIG. 13 .

 13 FIG.13. 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.

FIG. 15 .

 15 FIG. 15. 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 towards 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 Eq.(14b) as v → 0 since the distribution function has finite values at v = 0.

FIG. 16 .

 16 FIG. 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.

FIG. 19 .

 19 FIG. 19. Derivatives with respect to ϕ of the electron density ne, blue curve closed circles, and ion density ni, black curve open circles, compared to the kinetic calculation dnK /dϕ, red curve open triangles. While the derivatives of ne and ni converge in the quasineutral region ϕ > ϕ sheath , vertical dash dot line, the derivatives of ne and nK 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.

FIG. 21 .

 21 FIG. 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.

FIG. 22 .

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

FIG. 23 .

 23 FIG. 23. Phase space contour plot, plain blue curve open circles and open squares, of the electron distribution function log 10 (fe), 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 = K0 + ϕ -ϕ0 plotted for different values of the kinetic energy Ke0 at x = 0.

FIG. 24 .

 24 FIG. 24. Phase space contour plot, plain blue curve open circles and open squares, of the electron distribution function log 10 (fe), 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 = K0 + ϕ -ϕ0 plotted for different values of the kinetic energy Ke0 at x = 0.

FIG. 25 .

 25 FIG. 25. Phase space contour plot, plain blue curve open circles and open squares, of the ion distribution function log 10 (fi), coordinates ϕ -ϕ0 for the position and Ki sign(v) for the velocity, with Ki < 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 Ki = K0 -(ϕ -ϕ0) plotted for different values of the kinetic energy Ki0 at x = 0.

  Appendix A: Plasma-wall self-organization properties in the fluid framework 1. Fluid description of a 1-D plasma a. Navier-Stokes conservation equations

2 .

 2 Plasma sound waves a. Sound waves in quasineutral plasma As discussed in Section II A and illustrated by FIG.(

∂

  t n + ∇ s Γ = 0 (A36a) ∂ t Γ + ∇ s p e /m e = en U /m e (A36b) ∂ t Π e /m e + ∇ s 2( Q e )/m e = 2eΓ U /m e

TABLE I .

 I Reference simulation parameters.

	L ∥ /λD 0	Γw/Γ0	Ts/T0	mi/me	ν ⋆ 0
	≈ 202	1	1	400	0.1

TABLE II .

 II Reference locations in the simulation domain.

	≈ 95	195.9	202.6

source: xsce sheath: x sh wall: xw

Acknowledgments

This work has also benefited from discussions at the 12 th Festival de Théorie in Aix-en-Provence (2022, France). The work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No 101052200 EU-ROfusion) and training program and from the European Union's Horizon 2020 research and innovation program under Grant Agreement No. 824158 (EoCoE-II). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them. This work was supported by the EUROfusion Theory and Advanced Simulation Coordination (E-TASC) initiative under the TSVV-04 (Theory, Simulation, Verification and Validation) "Plasma Particle/Heat Exhaust: Gyrokinetic/Kinetic Edge Codes" project.

which yields the structure highlighted by the two first lines in the determinant Eq.(A14a), the third being determined by the chosen closure. Stepping back to Eq.(A10a), one finds that if the eigenmode is chosen such that N j = 0 for j > 0, then all the moment equations are equivalent to:

Further choosing the eigenmode such that u = 0, then leads to:

and the dispersion equation c = u for this particular eigenmode. Density fluctuations, with no other fluctuating field, are then found to exhibit a phase at rest in the frame moving at the plasma velocity u. From the kinetic point of view, the chosen fluid eigenmode corresponds to a perturbation of the distribution functions, identical for electrons and ions, of the form n g M (v, u, T M → 0) where g M is a Maxwellian with mean velocity u, thermal energy T M → 0 and v-space integral equal to 1. All moments of this perturbation are then proportional to n and to the mean velocity at some power, u ℓ for the moment of order ℓ. This property is readily recovered with the system Eq.(A13) when setting u = 0, T = 0 and q = 0, so that Γ = u n, Π/m i = u 2 n, and 2 Q/m i = u 3 n. For this eigenmode, one finds the dispersion relation c = u. This properly designed eigenmode enforces that one of the roots of the cubic dispersion equation is c = u. The procedure to specify the eigenmode generating a particular eigenvalue depends on the chosen fluid closure. For the dispersion relation Eq.(A14b), when requiring the system to exhibit the root c = u, one enforces following constraint on the closure coefficients.

This equation can be read as a constraint on one of the coefficients, a π , a Γ or a n , typically a n = u 3 -a π u 2 -a Γ u so that the dispersion relation allowing c = u as possible root should then also read:

For plasmas, unlike the case of standard neutral fluids, requiring this particular eigenmode and eigenvalue does not appear to be generic. For instance, for the moments ℓ > 3 both the electrical force and collisions will modify the dispersion relation. Furthermore, the closure that is chosen can rule out this particular eigenmode and eigenvalue. This issue will be further discussed in the following.

c. Sound wave with null heat-flux closure

For a quasineutral plasma, given T = T e +T i , one finds the following expressions for Π and the convective energy flux:

For a 1D geometry, with no energy exchange with other dimensions, γ c = 1 2 + 1. We generalize these expressions by introducing σ V and σ p , with γ p = σ p /σ V , analogous to the heat 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 Eq.(A20a) one then obtains:

With these coefficients the constraint Eq.(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 Eq.(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,

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, FIG.

(4), one finds that the sheath entrance is observed for u = 1.505 V , 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 Section A 2 d.

d. 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 Eq.(A12b) with zero source and Eq.(A12c) we then have:

According to Eq.(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:

One can then determine the expression of ∂ s M 2 .

Stepping to Fourier space and linearizing the latter expression, one obtains therefore:

Consequently:

Non trivial solutions with M 0 ̸ = 0 are obtained when the dispersion relation is fulfilled

By construction, the polytropic closure excludes the root c = u since N 2 is not null. However, the dispersion relation Eq.(A26b) allows one to recover the sound velocity

With no heat exchange γ p = σ p /σ n = 3 is the adiabatic index and one recovers the same sound velocity as obtained in Section A 2 c, Eq.(A22c). If one now identifies this result to the simulation result, 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.

e. Sound wave velocity with higher moment closure

The moment equation of order ℓ Eq.( A12a) is of the form:

For the sound waves we step to Fourier space, therefore:

As in the previous Sections, we shall consider here the case with null electric field and neglect the collisional contributions. If one sets Z ℓ = S ℓ-1 in Eq.(A27a), one can use this fluid hierarchy to determine the steady state solutions. Addressing Eq.(A27b) in the fluctuation framework, one assumes that there is no drive by source fluctuations, so that Z ℓ is identified to the time derivative of moment ℓ -1. With such assumption, the steady state and fluctuation equations, the drive term Z ℓ in Eq.(A27) is then written as:

when addressing the steady state solutions, while for the fluctuations in Fourier space one has:

As already highlighted, the fluid hierarchy in Fourier space exhibits a very simple structure given the phase velocity c = ω/k

For the purpose of the following calculation, it is useful to use the identity:

Linearizing Eq. (A10b) in Fourier space, one determines M ℓ-1 :

The second term can be simplified according to the identity Eq.(A30):

For ℓ ≥ 2 one then finds:

while for ℓ = 1:

From the latter equation one can determine u, which can be used in Eq.(A31a) to obtain for ℓ ≥ 2:

Closing the fluid hierarchy at moment ℓ, one finds the dispersion relation.

In this expression we have retained the term N ℓ to be able to address the specific closures for the moment ℓ = 2:

The polytropic closure could be used here leading to the results of Section A 2 d. As an alternative, we define

This corresponds to the isothermal closure which is in fact a particular case of the polytropic closure for γ p = 1.

Regarding the closure at moment ℓ = 3, setting N 3 = 0 corresponds to the closure q = 0 of Section A 2 c.

which yields the root c = u together with c 2 s = 3V 2 and the roots:

One recovers the result of Section A 2 c. If one now assumes N ℓ = 0 as rule for the closure at moment ℓ > 2, the existence of non-trivial solutions M 0 ̸ = 0 allows one to simplify the dispersion relation Eq.(A32b).

One can then generalize the result obtained for ℓ = 3, a root c = u together with the roots:

The sound velocity crucially depends on N ℓ-1 for which we have no particular insight for ℓ > 2. The roots will depend on the sign of N ℓ-1 as well as on the parity of ℓ. Odd values of ℓ will allow the roots ±c s , while even values require that c -u has the sign of the roots c s . One can note that applying Eq.(A35) to the case ℓ = 2, yields c s = 0. The results for ℓ = 2 have been obtained because we have used another closure than N ℓ = 0. Stepping to higher moments can be used to adjust the fluid prediction to the simulation result but one does not gain in predictive capability. Indeed, defining the transition into the sheath region with the Bohm criterion M 2 = 1, with M = u/c s , does not provide a boundary condition on u since c s depends on the chosen closure.

Given the result of Section A 2 c, recovered in Eq.(A33c), we shall consider c s = √ 3 V , unless specified differently, and, for convenience, we define The fluid closure we have chosen in Section III A is to set the heat-flux to zero. We have seen when comparing

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 Eq. ( 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 the electron convective energy flux 3 2 T 0 Γ. This indicates that the heat flux contributes to the energy flux and exhibits a contribution proportional to the particle flux typically of the form γ h ΓT e with γ h = 0.83. This result is consistent with the fact that the heat-flux contributes to the energy flux impinging onto the wall. This closure also modifies the calculation of the sound wave velocity as addressed in the following Section.

b. Modified sound velocity with non-collisional closure

With the non-collisional closure one has determined the linear dependence of the heat-flux on the lower order moments, namely n, Γ and Π. One can then compute the coefficients a n = ∂ n (2Q/m i ), a Γ = ∂ Γ (2Q/m i ) and a π = ∂ Π (2Q/m i ) stepping beyond the approximation Q = Q c as in Eq.(A21). The coefficients a n , a Γ and a π obtained with the non-collisional closure Eq.(A43), can then be used to determine the sound velocity and the departure from that previously obtained Eq. (A22). We concentrate on the real contribution since imaginary terms will govern damping or instability processes. We now consider Eq. (A14b)

where the three coefficients a π , a Γ and a n are modified by taking into account the non-collisional closure of the electron heat flux q e so that

are those previously used to determine the sound velocity Eq.(A21). The noncollisional closure yields coefficients that do not meet the constraint Eq.(A17a) and therefore exclude the root c = u.

Using Eq. (A44c) and defining δV 2 = 2γ h Te mi one then obtains:

The root such that c = u is excluded. We consider the root c 0 = u + δu and assume δu ≪ V so that the dispersion equation is then:

Assuming δu ≪ V , the dispersion equation is then:

At lowest order in δu, the constraint to ensure that c 0 is a root is therefore:

Then for c ̸ = c 0 one must solve:

The roots of this equation are then:

One finds that the sound wave with positive phase velocity c + is shifted to a sound velocity larger than √ 3V 2 while the sound wave with negative phase velocity, hence opposite to the fluid velocity for the right hand side sheath, is shifted to a modulus of the sound velocity smaller than

We therefore find that the fluid velocity that must be reached so that the sound waves propagating against the flow are at rest in the frame attached to the wall, is smaller than previously computed. The condition to generate a standing shock wave are therefore met for a mean velocity u smaller than √ 3 V . This trend agrees with the observation made in the reference simulation where the sheath horizon location suggests that u ≈ 1.5 V < √ 3 V , see FIG. 4 and discussion in Section II A and Section III B. Note that the downshift of the sound velocity for the sound waves propagating against the mean fluid flow computed here is only approximate since it has been determined in the limit δV 2 ≪ 3V 2 .

Asymptotic matching at the sheath horizon

We address here the plasma properties at the sheath horizon within the fluid framework. The horizon is assumed to be at z → +∞ where z is the curvilinear abscissa with origin at the wall and oriented towards the plasma. We then expand the various fluid moments in powers of 1/z [4]. Within the sheath, the total momentum flux of the plasma Π can be written as:

As done for the SOL plasma Q c stands for the convective energy flux, equal to the total energy flux when neglecting the heat flux. The asymptotic conditions are determined by assuming that Π is constant so that in Eq. (A51a), the variation of Π qn is balanced by the variation of the last term. One therefore expands Π qn in powers of 1/z and identify this expansion with the expansion of ηT e . One determines η = ρ c /e using the Poisson equation such that ρ c ∝ d 2 z δϕ. The term Π qn is shown to be a function of δϕ where δϕ is the variation of ϕ.

We choose as reference electric potential the asymptotic value at the sheath horizon ϕ ∞ . One therefore expects a shift by two orders of the expansion of η compared to that of δϕ and therefore of Π qn (δϕ). Such a shift together with the identification of the remaining leading order contributions, then yields consistency constraints at the sheath horizon. When the electron thermal energy T e is taken constant, the momentum balance for the electrons yields the standard adiabatic response for the electron density:

Choosing the space normalization defining z to be the Debye length at the sheath horizon, λ 2 D = ε 0 T e /(e 2 n ∞ ), the normalized charge density η = ρ c /n ∞ determined by the Poisson equation then allows one to determine the ion density n i :

Let us now define the δϕ functions that determine the ion density.

One can then write:

Due to scale separation, one can assume that the variation of both the particle flux Γ and energy flux Q c are small, consequently negligible at the relevant scale at the sheath horizon. One thus finds that the variation of Π qn only depends on that of the ion mean velocity u i = Γ/n i . Therefore, given Γ = Γ ∞ one can write:

Let us now define Π Q , Π Γ and Π η :

One can then write:

At the sheath horizon z → +∞ and δϕ → 0 so that one can expand the latter expression in powers of F (δϕ) -G(δϕ) given that by definition F (δϕ) → 0 and G(δϕ) → 0 so that F (δϕ) -G(δϕ) → 0.

and therefore:

If one now assumes the leading order term of F (δϕ) to be:

One can then determine the leading term in the expan

Since Π is assumed to be constant, Π = Π ∞ , one then obtains the two first constraints:

so that the remaining contributions to the expansion are then:

The constraint Π Q = Π Γ yields:

and consequently determines both u ∞ and Π:

Finally one must enforce that the chosen solution is not the trivial solution F (δϕ) = 0 and G(δϕ) = 0 of Eq. (A55c) by enforcing the lowest order constraint:

which yields:

To close this calculation one can note that F (δϕ) ≈ δϕ. Since A F K , δϕ ≤ 0, the sheath electric potential is smaller than the value at the sheath horizon, and similarly for the densities since n e ≈ n i ≈ δϕ. Since ρ c = -n ∞ G(δϕ), and given the constraint Eq. (A55c), one finds that:

One therefore recovers that the charge density within the sheath is positive n i ≥ n e and that the charge density derivative 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 , Eq. ( 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 towards 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 Eq.(A56a), and dn e /dϕ ≥ dn i /dϕ