Yann Munschy
email: yann.munschy@cea.fr

Emily Bourne

Guilhem Dif-Pradalier

Peter Donnel

Philippe Ghendrih

Virginie Grandgirard

Yanick Sarazin

Kinetic plasma-wall interaction using immersed boundary conditions

The interaction between a plasma and a solid surface is studied in a (1D-1V) kinetic approach using immersed boundary conditions and penalization to model the wall. Two solutions for the penalized wall region are investigated that either allow currents to flow within the material boundary or not. Essential kinetic aspects of sheath physics are recovered in both cases and their parametric dependencies investigated. Importantly, we show how the two approaches can be reconciled when accounting for relevant kinetic effects. Non-Maxwellian features of the ion and electron distribution functions are essential to capture the value of the potential drop in the sheath. These features lead to a sheath heat transmission factor for ions 60% larger than usually predicted. The role of collisions is discussed and means of incorporating minimally-relevant kinetic sheath physics in the gyrokinetic framework are discussed.

I. INTRODUCTION

The physics of plasma/material boundary interfaces is akin to that of time-dependent boundary conditions for hyperbolic systems in which discontinuities and shock waves develop [1]. Langmuir and Tonks were the first to study the properties of a plasma in contact with a solid surface [2,3]. The discontinuity region at the plasma/boundary interface needs to be resolved through departure from quasineutrality across a narrow region [4] which scales with the Debye length, namely the reference scale for the electric field generated by charge separation.

The reference problem in plasmas is that of an open field line with a source of particles and energy, possibly of momentum, which extends over the whole domain (typically L ∥ ) and intercepts on both ends a material boundary. Characteristic distances in fusion plasmas are such that the connection length L ∥ ≳ 10 m whilst the Debye scale λ D is in the range of λ D ≲ 10 -4 m, hence leading to a typical small parameter ℓ = λ D /L ∥ ≲ 10 -5 .

With respect to the incoming plasma, a surrounding dense gas, a liquid or a solid are seen as dense and neutral "walls". Boundary conditions for the plasma are thus governed by the respectively very large density of these surroundings, which governs a prompt recombination of the impinging plasma ions, so that the boundary condition acts as a close to perfect particle, momentum and energy sink. Unlike neutral fluids, the plasma disappears into the boundary material, which drives a plasma flow into it. Furthermore, it is generally assumed that the boundary material is supposed to be at constant conditions, and grounded. It thus further acts as a thermal bath and as an infinite reservoir of particles.

From a fluid point of view, it can be shown that the sheath region where there is a departure from charge neutrality must be in a supersonic regime whilst the presheath [5], which exhibits a stagnation point due to the assumed antisymmetry of the plasma flow must remain subsonic. The sheath boundary layer can be understood as a region where the plasma self organizes such that an electric field develops in order to confine bulk plasma electrons, preventing macroscopic charge separation. This yields conditions for a dynamical steady state such that the volumetric source of particles and energy is balanced by the plasma particle and energy loss by recombination in the boundary material. In the sheath, the plasma becomes non-neutral. Dynamical polarization of the sheath, whilst ensuring confinement of mobile electrons contributes to accelerating ions so that at the plasma/sheath transition net zero plasma current is achieved statistically, in steady state.

The presheath is understood here as the plasma region where the source effects prevail and where quasineutrality conditions are met, in contrast to the sheath region where the sources, in particular the particle source become negligible and departure from quasineutrality occurs. Asymptotic matching [4,6] between the subsonic presheath and supersonic sheath is often used to yield the typical distance from the wall where the so-called Bohm condition is met [START_REF] Bohm | The Characteristics of Electrical Discharges in Magnetic Fields (ch 3[END_REF] -where the ion flow becomes supersonic. This location usually defines the sheath entry point. This criterion results from the strong scale separation between the plasma typical extent L ∥ and the sheath typical width λ D . In this problem, the fact that most electrons are reflected by the sheath potential whilst all the ions are absorbed by the wall material significantly distorts the distribution functions away from Maxwellian distributions. This indicates that the kinetic rather than the fluid framework must be considered to investigate the sheath boundary region, as first formulated by Harrison and Thompson [START_REF] Harrison | Proceedings of the Physical Society[END_REF] and then clarified by Riemann [6,[START_REF] Riemann | [END_REF][10][11]. In this context, taking into account collisional processes was early recognized to be of great importance [START_REF] Boyd | Proceedings of the Physical Society. Section B[END_REF][START_REF] Persson | [END_REF]. An extensive review of plasma models in the vicinity of the sheath, along with numerical methods that can be used to study plasma-wall interac-tion, can be found in Ref. [14].

Plasma-wall interaction has strong implications in magnetic fusion devices. The need to understand, predict and control heat and particle deposition patterns is indeed a major concern for Iter. Erosion, re-deposition, impurity control and the onset and stability of the edge transport barrier that develops in the vicinity of the separatrix in H-mode plasmas [15] all greatly depend on the Scrape-Off Layer dynamics [16][17][18][19] which in turn are strongly influenced by boundary layer physics and plasma-wall interaction processes. Fluid approaches have long been developed to take into account the key features governed by the kinetic physics without addressing the sheath physics. In this case, plasma-wall interaction is usually treated by imposing a supersonic outflow as boundary condition [20] along with specifying the outgoing plasma energy flux using the so-called sheath heat transmission factor. In many simulation settings, the latter condition prescribes a nonphysical heat diffusion across the boundary. Conveniently for fusion plasmas, the mesh size is much larger than the sheath; the constraints governed by the latter are taken into account but the sheath is not resolved. It has been shown that completing the standard fluid framework based on particle momentum and energy conservation equation, the Navier-Stokes equations, by considering higher fluid moments that capture the departure from Maxwellian distributions, allows one to recover key aspects of the sheath physics [START_REF] Joncquieres | Joint Propulsion Conference[END_REF]. Although interesting from the numerical point of view, this approach requires that one meshes at the Debye scale and does not alleviate the usual question regarding the fluid closure. Reference to the kinetic solution remains crucial to validate such results. Approaches that capture fine kinetic details of sheath physics exist and have allowed detailed studies of a variety of situations: regimes with multiple ions [START_REF] Tskhakaya | [END_REF], the effects of E × B drift [START_REF] Tskhakaya | [END_REF]23] or of grazing angle magnetic fields [24][START_REF] Coulette | Plasma Physics and Controlled Fusion[END_REF][START_REF] Moritz | [END_REF][27]. However, the self-consistent description within a kinetic or gyrokinetic framework of the plasma edge, the Scrape-Off layer and the sheath is beyond current numerical capabilities. The scope of this paper is to investigate the kinetic properties of the global sheath and presheath physics by implementing immersed boundary conditions that could be extended to the gyrokinetic framework, where the mesh step is typically 10 -2 L ∥ , i.e. 10 3 Debye length. Such immersed boundary conditions have been implemented in fluid models for plasma-wall interaction using a penalization technique [28][29][30][31], and the required spatial organization of particle and momentum sinks to recover the appropriate sheath behavior have been documented [4]. These results need to be extended to the kinetic framework.

Having these questions in mind, we thus restrict ourselves in the current paper to investigating selforganization of the plasma-wall interface under the combined outfluxes of electrons and ions, at normal incident angle and without any interaction with neutral particles. Note that the model considered in the present workone-dimensional, with a magnetic field line that is considered to impact a wall with normal incidence-is completely equivalent to studying a non-magnetized plasma in contact with a solid wall. The problem is therefore that of reaching steady state plasma conditions such that the plasma volumetric source balances absorption due to recombination at the surface of field-line intercepting materials. Material boundaries are assumed to remain neutral, i.e. to be totally unaffected by plasma recombination and heat deposition at their surface.

We show that plasma-wall interaction-including sheath physics-with normal incidence of magnetic field lines onto the wall can be adequately described within a kinetic framework by means of immersed boundary conditions, much in the spirit of Ref. [31]. The results obtained in such framework and presented in the following are complemented by two companion papers. Together with the present work these papers aim at providing a self-contained study of kinetic plasma-wall interaction from both numerical and physical standpoints. The first of the two companion papers adresses the comparison of the kinetic simulation results with the fluid theory and highlights in particular that the fluid framework cannot describe the parallel physics at play in the kinetic simulations [32]. It further questions the validity of Bohm's criterion-and of its kinetic equivalent, the Harrison-Thompson or Riemann criterion [6,33]-to define the entrance of the sheath. Collisional processes appear to be a key player in explaining the departure of the simulations results from the fluid theory. This departure becomes especially blatant when computing the heat flux from the simulation results. The second companion paper investigates how to deal with inhomogeneities of physical quantities that appear in the sheath region, in particular steep gradients-at or below the typically resolved scales of turbulence [34]. Specific numerical techniques are indeed needed to efficiently handle the strong scale separation that exists between the quasineutral plasma and the sheath. This second companion paper shows the use of non uniform cubic splines, along with GPU parallelization. This alleviates those numerical issues leading to fast and well resolved simulations. In particular, it is shown that non-uniform splines lead to a reduction of memory requirements by 89%, whilst keeping numerical oscillations due to the interpolation step of the numerical scheme to a low level. Even if non-uniform spline schemes are more costly than their uniform counterpart, parallelization on GPUs lead to fast simulation times. Overall, a non-uniform simulation runs 5.5 times faster than a uniform simulation resolved enough to provide equivalent physical results.

In the present paper we therefore discuss the features of the simulated plasma-wall interaction with the Voice code [START_REF] Stangeby | The Plasma Boundary of Magnetic Fusion Devices[END_REF]. This code solves a Boltzmann-Poisson system with a numerical scheme which is very close to the one implemented in the 5D gyrokinetic code Gysela [START_REF] Grandgirard | [END_REF].

Namely, the coupling between Boltzmann and Poisson equations is solved using a predictor corrector scheme. The Boltzmann equation solver further involves a splitting of Strang, with a semi-Lagrangian method to describe the advections. The code describes a 1D1V phase space with possibly non uniform meshes in both spatial and velocity directions. The Voice code has already been used as a test bench for non-uniform splines that are to be implemented in Gysela [34]. Using this code we present a detailed analysis of the robustness of key sheath outcomes-namely the potential drop in the sheath and the shape of the ion and electron distribution functions-to the impact of reduced electron to ion mass ratios, plasma source and system size. Establishing this robustness is especially helpful for reduced models (fluid or gyrokinetic) for which nominal parameters may not be readily accessible due to computational limitations. In Section II we present the model of immersed boundary conditions that we use to treat the problem of the interaction between a plasma and a solid surface. Two other possible solutions that allow the integration of the plasma-wall interaction into a global kinetic or gyrokinetic code are presented as well in Section V. General features of simulations of plasma-wall interaction are discussed. In Appendix C the code verification is presented, and the parameters used in the simulation are discussed. Section III is devoted to the thorough analysis of the kinetic results. They are put in perspective with the more familiar fluid-like approach where distribution functions are assumed to remain close to Maxwellians, either shifted of truncated. In particular, the issue of the sheath heat transmission factor is addressed. Then, in Section IV we assess the influence of non-nominal parameters (system size, electron to ion mass ratio, electron temperature) on key observables relevant to more complex and numerically-intensive codes (ion and electron distribution functions, potential drop). On the basis of the previous results, the discussion in Section V explores possible ways to implement sheath-like boundary conditions -without resolving the sheath-in flux-driven gyrokinetic codes such as Gysela [START_REF] Grandgirard | [END_REF] with the help of immersed boundaries.

II. KINETIC MODELING OF PLASMA-WALL INTERACTION USING IMMERSED BOUNDARY CONDITIONS

We consider a Scrape-Off layer magnetic field line represented on FIG. (1), left graph. This line is connected to the divertor wall at both of its ends. Curvature effects are neglected and the field line is flattened, giving the geometry of FIG. (1), top right graph. This flattened line impacts the divertor with a normal incidence, the treatment of a magnetic presheath in the case of oblique incidence being left for future work. The effect of perpendicular transport between field lines is restricted to generating particle, momentum and energy source terms standing for the divergence of the respective fluxes. One can thus restrict the analysis to a single field line. The model can be extended to handle a charge source as would result from the inhomogeneity of the magnetic field amplitude. With such simplifications we restrict the problem to onedimension in position space and one dimension in velocity space. We use the s variable to specify the considered position on the field line. The v variable expresses the parallel velocity of the considered particle. We furthermore introduce the distribution function f a (s, v, t), expressing the density of particles at time t and at point (s, v) of phase space. The subscript a denotes the species the quantity refers to: e for electrons, i for ions. The mass of any particle of species a is written m a , and its charge is e a . The electric field is written E, the magnetic field B. In the present work, the magnetic field is assumed to be constant and homogeneous along the considered field line. A 1D-1V Boltzmann Eq. ( 1) is used to describe the evolution of the distribution function f a over time in phase space.

∂ t f a + v∂ s f a - e a m a ∂ s ϕ ∂ v f a = C(f a ) + S(f a ). (1)
The C(f a ) term accounts for collisions. The S(f a ) operator represents sources and sinks. As stated in Section I, this Boltzmann equation can be used to describe either a one-dimensional plasma along a single magnetic field line that impacts the wall at normal incidence, or a one-dimensional non-magnetized plasma in contact with a solid wall. These two situations are formally equivalent. In the electrostatic limit the electric field is -∂ s ϕ. The latter is determined by the Poisson equation (2), the electrostatic limit of the Maxwell-Gauss equation, which relates the electric potential to the charge density ρ c .

-∂ 2 s ϕ = ρ c /ε 0 with ρ c = species e a n a , (2)
where we introduced the local density n a of species a, which is a function of time and space. It is defined as the integral over the velocity space n a = dv f a . The Boltzmann-Poisson system is normalized using a reference density n 0 and temperature T 0 . Time is normalized to the inverse of the electron plasma frequency ω pe0 = n 0 e 2 /m e ε 0 . The space variable s is normalized to the length scale relevant for plasma-wall interaction studies, i.e. the Debye length λ D0 = ε 0 T 0 /n 0 e 2 and we set x = s/λ D0 . The electrostatic potential is normalized to T 0 /e. The phase space velocity variable v is normalized to the thermal velocity of each species v T0a as v a = v/v T0a , with the thermal velocity of species a written as v T0a = T 0 /m a . Particle distribution functions are normalized to the reference particle density in phase space n 0 /v T0a . We use the notation A a = m e /m a to express the mass ratio between electrons and species a (in particular A e = 1). The normalized charge of species a is written as Z a = e a /e. The normalized Boltzmann- I.

Poisson system then reads as follows

∂ t f a + A a v a ∂ x f a -Z a ∂ x ϕ ∂ va f a = C(f a ) + S(f a ), (3a)
∂ 2 x ϕ = -ρ c with ρ c = species Z a n a . (3b)
In the remainder of the paper, each quantity of interest is normalized. If a reference to some dimensional quantities is needed, it will be explicitly pointed out. The right hand side of the considered Boltzmann equation is composed of a collision term C and S, which stands for both sources and sinks. The collision operator accounts for self-collisions C aa and inter species collisions C ab . The S operator is made of three different components. First, the particle sink term S sink,n that is the critical term for the immersed boundary condition technique. As in previous work using this technique [4,31], we do not impose explicit boundary conditions on the distribution function at the plasma-wall interface. Rather, we include the specific term S sink,n in Eq. (3a) that acts as particle sinks in the wall region. This term is discussed in detail in Section II A. Second, a source term S src that injects particles and energy in the plasma-and defines the source region. This term can readily be changed to inject heat or momentum. In the present work, the source region is localized, and does not extend along the whole field line. Lastly, an energy sink term S sink,T whose purpose is to cool the remaining electrons in the wall region, and ensures that these electrons exhibit a low temperature Maxwellian distribution. The latter term provides a useful damping mechanism of numerical artifacts in the wall region, in particular when we allow currents to exist there.

A. Immersed wall models

Because of the very large electron density available in the wall material, ions leaving the plasma and entering the wall promptly recombine. Both ions and electrons penetrating the wall are slowed down until their energy is reduced to the wall thermal energy. From the plasma point of view, the wall acts as a perfect particle, momentum and energy sink. Because of the inertia difference between electrons and ions, the electron flux leaving the plasma would be significantly larger than the ion flux when assuming comparable temperature for the two species and free particle motion. Consequently, a large parallel electric field develops at the plasma-wall interface to confine the electrons, reducing their effective mobility. In the simplest form, the steady state self-organization of the plasma in contact with the wall consists of a bulk quasineutral plasma surrounded by a boundary layer at the plasma wall interface. There, quasineutrality breaks down due to the large electric field repelling the electrons. This boundary layer is called the sheath. It regulates the ion and electron flux to the wall so that a steady state can be reached, given the plasma particle and energy source that must balance the losses to the wall. The charge imbalance, which defines the sheath region, also governs its characteristic scale as determined by the Poisson equation (2), namely the Debye scale λ D . On average, the sheath ensures ambipolarity at the plasma boundary, i.e. it prevents net charge build-up within the plasma region. The analysis based on matching sheath and presheath conditions indicates that plasma self-organization governs plasma acceleration due to the source terms, in particular the particle source, in the region of quasineutrality, and that the actual sheath extent depends on a balance between the electron repulsion by the electric field and plasma acceleration by the source [4]. It is important to keep in mind that there is no abrupt transition between the quasineutral plasma and the sheath region where the charge density becomes comparable to the particle density. The key difference is a matter of scales. From Poisson equation we see that a potential drop ∆ϕ spread across a characteristic gradient length ℓ relates to the charge density ρ c that gives rise to this potential drop as

e∆ϕ T 0 ≈ - ℓ 2 λ 2 D ρ c n 0 e ,
where all quantities are expressed in dimensional units and n 0 , T 0 are the plasma typical density and temperature. The potential drop in the presheath-or quasineutral region-is of the same order of magnitude as the potential drop in the sheath, as we have typically e∆ϕ presheath /T 0 ≈ -0.7 and e∆ϕ sheath /T 0 ≈ -3 [START_REF] Stangeby | The Plasma Boundary of Magnetic Fusion Devices[END_REF].

The typical gradient length in the sheath is given by the Debye length, i.e ℓ sheath ≈ λ D whereas in the presheath the gradient length is given by the length of the magnetic field line, i.e. ℓ presheath ≈ L ∥ . For a typical magnetic field line of the Scrape-Off layer one has L ∥ ≫ λ D , thus the charge density in the sheath is very large as compared to the charge density in the presheath as

ρ presheath c /ρ sheath c ≈ λ 2 D L 2 ∥ ≪ 1.
Therefore the electric field that develops in the sheath is also very large as compared to the presheath electric field since from Poisson equation we have

E presheath /E sheath ≈ λ D L ∥ ≪ 1.
It should be highlighted that the analysis conducted in the present section and more generally in the present paper is valid only when considering a magnetic field line that impacts the wall at normal incidence. When going towards more realistic model that takes into account tilted field lines with respect to the divertor/limiter, the distinction between a quasineutral presheath and a positively charge Debye sheath becomes irrelevant. Indeed, in this situation a magnetic pre-sheath develops-the socalled Chodura sheath-that modifies the properties of the plasma close to the material boundary [24]. Furthermore, under realistic divertor conditions where the magnetic field lines impact the divertor wall at shallow angles it has also been shown that the Debye sheath may simply disappear [37].

Taking into account plasma-wall interaction in simulations implies dealing with the physics of these two plasma regions. However, in well resolved plasma turbulence simulations the mesh size in the parallel direction is proportional to the ion Larmor radius. Moreover, this mesh size is usually larger than the Larmor radius because of the slow variation in the parallel direction compared to the transverse direction that governs the ion Larmor radius scaling. In simulations of interest, the ion Larmor radius is an intermediate scale between the Debye scale and the plasma scale L ∥ . The sheath physics, described above, and which account for plasma-wall interaction, are then unresolved. Within the fluid equation framework, the standard approach is to set the boundary of the simulation domain at the sheath entrance. One then prescribes a particular dependence on the local plasma parameters for the particle and energy fluxes as well as the electrical current. These depend on a few dimensionless parameters, independently determined by the sheath physics, namely the Mach number at the sheath entrance, the sheath energy transmission factor, and the potential drop within the sheath normalized by T e /e, where T e is the electron thermal energy. When contemplating turbulence simulations in the gyrokinetic framework one must set as many boundary conditions as mesh points used to describe the distribution function dependence on the parallel velocity. A generic formulation for such a task is not straightforward and alternative means are currently used, these are described in Sections V 1 and V 2 in Section V. We investigate here an alternative path, namely penalization, already successfully used for plasmas in the fluid description [31]. There, penalization based on a simplified model of the plasma-wall interaction allows one to recover the sheath constraints without resolving the sheath scale [4,31]. Before stepping directly to this approach in a 5D gyrokinetic framework, we investigate here the plasma and sheath self-organization in a 1D-1V kinetic framework using an immersed boundary and the same simplified wall model. We can then analyze consistently the main kinetic properties that must be recovered whenever one addresses kinetic boundary conditions without resolving the sheath.

The wall penalization by a particle sink

The immersed boundary condition is introduced in the Boltzmann Eq. (3a) using the particle sink term S sink,n . It is a Bhatnagar-Gross-Krook operator [30] that acts as a particle, momentum and energy sink. It is written as

S sink,n (f a ) = -M w ν wa (f a -g w ) . (4)
The structure of this operator is that of a restoring force that drives the distribution function f a (v a ) towards the low density and cold target distribution g w (v a ) For the sake of simplicity, this target function is a Maxwellian expressed in normalized units as

g w (v a ) = n w √ 2πT w exp - v 2 a 2T w . (5)
It is chosen to be species independent. The density n w and temperature T w are constant in space and time. The wall target density n w is orders of magnitude smaller than the plasma density, so the operator acts as a particle sink.

Typically we take n w ≈ 10 -9 in the simulations. The choice for the particular value taken for n w is arbitrary, but is motivated by physical and numerical considerations. The density of the wall should be as small as possible as compared to the plasma typical density, so that the wall effectively acts as a sink of particles. From a numerical point of view, it is easier to take n w ̸ = 0 to avoid negative values in the distribution function. Indeed, as the voice code uses a semi-Lagrangian scheme for solving the advective part of a Boltzmann equation, an interpolation method is required. Typically we use cubic spline interpolation to perform this procedure. Interpolating data points that are very close to zero may lead to negative values as the spline often shows small oscillation around the interpolated data points [34]. To avoid this, the chosen value for the wall density is non-zero. The authors have checked that the potential drop in the sheath region is independent of the particular value of n w as long as this parameter remains small as compared to the typical plasma density. The Maxwellian has a zero mean velocity so that the momentum of the distribution function f a is removed by the BGK sink. The target temperature T w is smaller than the source temperature that injects particles in the plasma, see Section II B 1. Because of the finite resolution of the velocity space in the simulations, this temperature is however much larger than the wall temperature achieved in experimental conditions. The particle sink operator is localized using the mask function M w (x). This mask function therefore defines the wall region in the simulation domain. It is expressed as

M w (x) = 1 - 1 2 tanh x -x ℓ w d w (6) -tanh x -x r w d w .
This function is represented on FIG.

(1), bottom left graph. The x ℓ w and x r w parameters define the spatial positions where the mask function is equal to 1/2 as it steps down from 1 in the wall to 0 in the plasma. The wall region is therefore defined by M w (x) = 1 while the plasma region is such that M w (x) = 0. Intermediate values characterize gray regions where the wall sink is active but the plasma density is still larger than the target value n w . The characteristic width of this intermediate region is d w , typically a fraction of the Debye length, see Table I. The plasma-wall interfaces are therefore blurred, which must be accounted for when analyzing the results. Their characteristic location are x ℓ w (left) and x r w (right), these being symmetric with respect to the middle of the simulation box. The third parameter in the chosen restoring force Eq. ( 4) is the species dependent amplitude ν wa . Here we fix the value ν wi for the ions, hence enforcing the recombination rate, and address two different particle loss amplitudes for the electrons. With the operator presented here, the absorbing wall boundary condition is therefore "immersed" in the simulation box. In other words, from a numerical point of view the wall appears in FIG. 2. Treatment of the material boundary. The mathematical boundary condition is relegated away from the region of interest where the plasma-wall transition takes place. The wall region acts as a particle sink. This region plays the role of an immersed boundary condition from the plasma point of view.

our simulations as a plasma with altered (i.e. penalized) dynamics. The boundary layer at the wall interface then develops independently of any numerical boundary conditions that are necessary to solve the Boltzmann-Poisson system given in Eq. (3a-3b). Actual boundary conditions are imposed at x = 0 and x = L x , far from the sheath region, see FIG. (2). The penalization method allows recovering a perfectly absorbing plasma boundary in the limit of an infinitely steep wall, a zero wall density, and an infinitely strong restoring force in the wall. That is to say in the limit where (i) d w → 0 in the definition of the wall mask Eq. ( 6), (ii) n w → 0 in the definition of the target Maxwellian Eq. ( 5) and (iii) ν w → +∞ in the expression of the BGK operator Eq. (4). Indeed in this case the wall mask takes a value of one in the plasma, zero in the wall and the transition between these two regions does not exist anymore. In the plasma the restoring force is therefore not active S sink,n = 0, hence the plasma is not affected by the restoring force. To be completely rigorous we need to add the constraint ν w M w → 0 in the plasma region. In the wall region the restoring force is infinitely strong, hence f a (x ∈ wall) = g w = 0. Thus, any particle (electron or ion) that comes into contact with the wall is absorbed.

Setting ν we = ν wi , therefore constant in time and space one finds that a current and charge separation is enforced by the penalization term. Indeed, considering the simplified kinetic problem governed by the equation v T a ∂ x n a (x) = -ν a n a (x), one finds an exponential fall-off where the e-folding length v T a /ν a is the wall stopping distance for particles of species a. We have chosen here the thermal velocity v T a as characteristic velocity. For the chosen absorption amplitude ν we = ν wi , one finds that the stopping distance is much longer for electrons compared to that of ions, which then effectively acts as a negative charge and electrical current source. In the following, we refer to the wall with equal sink amplitude as the conducting wall because it is found to exhibit an electrical current with reduced charge density, as will be shown in Section II A 4. This naming is only for convenience because neither the wall resistivity nor any electrical current properties of the wall are addressed in the chosen reduced wall model. An alternative choice is to ensure that the chosen wall particle sink does not generate any charge. For an electron and singly charged ion plasma, assuming that the BGK amplitude does not depend on velocity, then according to Eq. ( 4) the normalized charge source is S ρ = -M w {ν we (n e (x, t) -n w ) + ν wi (n i (x, t) -n w )}. Then setting:

ν we (x, t) = ν wi n i (x, t) -n w n e (x, t) -n w . ( 7
)
one enforces the charge source to be null. In this wall model the ion recombination leads to an electron absorption by the wall at the same location, the difference in location is assumed not to be resolved. This process then removes the electrical currents that must appear in the conducting case to connect the different ion and electron absorption locations, and which tend to recover neutral conditions in the wall. In simulations with this wall model, hence ν wi = constant and ν wi determined by Eq. ( 7), one observes a charge build-up in the wall region, while the electrical current remains small, Section II A 4. Given these properties and the chosen name for the first wall model, we shall refer to this case as the insulating wall although this actually makes little sense for a region that is understood to be an infinite and grounded charge reservoir. One can note that an intermediate model would be interesting combining constant amplitude of the restoring forces, ν wi = constant and ν we = constant, but such that ν we = ν wi √ m i / √ m e so that the ion and electron stopping distances are comparable. This would tend to reduce the charge source within the wall without using a restoring force amplitude depending on the particle densities, and which can generate singular properties whenever n e or n i tend to n w .

3. Electron heat sink in the insulating wall model.

One can notice in the constraint Eq. ( 7) that as n i (x, t) → n w , ν we (x, t) → 0. Therefore, once the bulk ions are absorbed so that n i (x, t) → n w , the absorption of the fast electrons stops, so that the latter then freely propagate further into the wall and are then constrained by the chosen boundary conditions imposed at the two ends of the simulation domain. A specific heat sink S sink,T is then used to ensure that these remaining electrons tend to a low temperature Maxwellian distribution function. This operator is a Bhatnagar-Gross-Krook 3. Total current normalized to the ion saturation current in the sheath vicinity, for the conducting wall in plain red, and the no-current wall in dashed blue. The spatial axis has been renormalized so that the sheath entrance and wall position coincide in the conducting and no-current cases. The wall mask defined in Eq. ( 6) is represented in dashed grey.

term of the form

S sink,T (f e ) = -M sink,T ν sink,T (f e -g sink,T ). (8)
The mask M sink,T is similar to the wall mask of Eq. ( 6) but restricted to the deeper half portion of each wall, see FIG. (1). This is to ensure that the heat sink term does not interplay with the plasma-wall transition region. This heat sink coexists with the particle sink Eq. ( 4). The amplitude ν sink,T of the restoring force is chosen constant. The target function g sink,T is a Maxwellian

g sink,T (v) = n sink,T 2πT sink,T exp - v 2 2T sink,T .
It is characterized by a local density n sink,T (x, t) equal to the electron density n sink,T (x, t) = f e (x, v, t) dv at each time step and position. This operator thus conserves density locally, it does not absorb particles. We set the temperature of the target distribution function T sink,T as equal to a fraction of the electron temperature at the wall position: T sink,T (t) = αT e (x w , t), where x w is the position where the wall mask defined by Eq. 6 steps from zero to one. Typically we take α = 1/3. This operator then acts as an energy and momentum sink for the remaining electrons in the wall. Furthermore, it restores the electron distribution function towards a Maxwellian. Note that this heat sink operator is necessary only for the insulating wall model. 4). The potential is represented in black, in the conducting wall case. ϕ sh is the electric potential taken at the sheath entrance in the insulating wall case.

Conducting and insulating wall behavior.

To compare the wall models that we have named conducting wall and insulating wall, we anticipate here on the results of Section III and IV dedicated to the simulations. The location of the sheath entrance and wall position are discussed in these latter Sections and specifically in Section III A. We use them for convenience without specifying their exact definition at this stage. We focus on the current and charge density in the wall region, mostly in the vicinity of the sheath, since the penalization terms are set to achieve evanescent plasma conditions deeper into the wall region. Restricting the plasma to electrons and one species of singly charged ions, the dimensional form of the total electrical current is by definition:

j tot = j e + j i = q e Γ e + q i Γ i ,
where Γ a is the particle flux of species a. On FIG.

(3) the total electrical currents are plotted in the sheath region at the right-hand side plasma-wall interface for the two wall models. A negative current can therefore be associated to a positive flux of electrons. The total current is normalized by the ion saturation current at the sheath entrance, j sh sat = q i Γ sh i . For the conducting wall model, the current, plain red curve, is non-vanishing in the wall region and is non negligible compared to the ion saturation current. The dash-dot blue curve corresponds to the current in the insulating wall model, where the ν w coefficients are adapted to prevent strong currents being directly driven by the particle sink term. In both cases we remark that the total current is negligible at the sheath entrance compared to the large particle flux. Furthermore, it is still small at the chosen wall position: we retrieve a fundamental feature that one expects to occur within the sheath: ambipolarity is ensured at the plasma boundary to prevent continuous build-up of the plasma charge, a necessary condition to sustain quasineutrality in the plasma bulk. In steady state the normalized charge balance equation is simply written as

∂ x j tot = species Z a dv S sink,n (f a ) = S ρ (9a) S ρ = -M w ν wi n i -n w -ν we n e -n w (9b)
The right hand side of Eq. ( 9a) is the charge source term S ρ . For the case of electrons and singly charged ions it is given in Eq. (9b). In the plasma region where M w = 0 the charge source is null and a constant current is solution of Eq. ( 9a). With such a solution a current flows from one end point of the field line to the other. This is physically possible, as with thermo-electric currents. In the present work we assume a symmetry with respect to the middle of the box so that the constant current is constrained to be null. In the conducting wall model with ν wi = ν we , one obtains S ρ = -M w ν wi (n i -n e ). Charge separation and non-vanishing wall mask then govern a variation of the electrical current. Given the finite steepness of the wall mask, the charge source becomes effective prior to the defined wall position. The maximum steepness of the wall mask is constrained on the one hand by the refinement of the spatial grid. The more points are used to describe the wall mask transition region, the steeper the wall mask can be. On the other hand, the steepness of the wall is also limited by the value of the timestep adopted in the simulation. The VOICE code used to perform the simulations operates with a splitting of Strang to solve the Boltzmann equation. Thus, the routine to solve one timestep can be written as follows:

1. Solve ∂ t f a = S(f a ) on dt/2; 2. Solve ∂ t f a = -v∂ x f a -q a E/m a ∂ v f a on dt; 3. Solve ∂ t f a = S(f a ) on dt/2;
Where we did not mention the collisions in the splitting scheme for the sake of simplicity. This splitting is no longer valid if the effect of the source term S(f a ) on the particles trajectories cannot be neglected when performing the advection step (b). This can happen if the timestep if too large, or if the source term S(f a ) presents a very steep spatial variation, as it is the case when the wall mask is too steep. Having these constraints in mind we take the maximum steepness which is compatible with the chosen discretization refinement and timestep. The steepness of the wall mask that we have considered has been checked to be large enough, that the observed sheath physics become independent of this parameter. The positive charge separation within the sheath then drives a negative current in the wall, FIG.

(3) plain red curve, until the difference between the electron and ion stopping distance generates a negative charge separation so that the wall current increases again and gradually tends to zero, FIG.

(3). In the insulating wall model the amplitude of the electron sink is defined in Eq. ( 7) so that S ρ = 0. As observed on FIG.

(3), dash dot blue curve, the electrical current remains close to zero. On FIG. ( 4) are compared the steady state charge density ρ c profiles in the sheath region for the two wall models conducting, plain red curve, and insulating, dash-dot blue curve. In both cases the region between the chosen sheath entrance and wall positions is characterized by a positive and increasing charge density as expected in the sheath region. In the wall region, the positive charge density is observed to decrease and changes sign. The region of negative charge in the wall, spans over a larger distance in the insulating wall model compared to the conducting wall model, see FIG. (4). In the conducting wall model, the charge density and electrical current behavior are consistent with the particle source properties, which govern a larger stopping distance for the electrons than for the ions. In the insulating wall model, charge conservation is enforced for the complete simulation domain. Therefore, the positive charge density that prevails in the plasma region must be balanced by a negative charge density in the wall region. The proposed wall models-conducting and insulating-are found to provide appropriate conditions to investigate the plasma response when interacting with the wall. However, these models are too crude to properly describe the wall behavior, in particular with respect to currents, charge and consequently electric potential variation.

B. Plasma model

Plasma source

The source/sink term S(f a ) of Eq. (3a) is composed of three distinct entities S = S sink,n + S sink,T + S src . The two sink terms, S sink,n and S sink,T characterize the wall region. They act as particle, momentum and energy sinks. To achieve steady state conditions a source is mandatory. In the SOL region of magnetically confined plasmas, cross-field turbulent transport is found to govern the heat source. Depending on the operating conditions, neutral particle ionization on the open field lines and cross-field turbulent transport determine the particle source. The chosen source term S src for particles and energy is similar to the Gysela source terms [38]. While in Gysela this term stands for actual particle, momentum and energy sources 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. Considering that turbulent transport is observed to be ballooned to the low field side, 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 as well as spurious external heating of SOL electrons [39]. For the sake of simplicity, we shall not address in this paper such cases nor the possibility of momentum and charge sources. We focus here on conditions which tend to 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. The SOL problem of interest is therefore that of the hot plasma regime where atomic processes are expected to have a relatively weak effect. This particular regime has implications when addressing the plasma collisionality, see Section II B 2. The particular self-organization problem is to determine the plasma total pressure and the plasma temperatures on the open field lines, given the wall losses monitored by the parallel transport and the sheath constraints on the fluxes. Together with the source conditions, the latter strongly depends on the cross-field transport properties. To model the source we use a BGK-like source but without the restoring contribution, similar to that implemented in the Gysela code [38]. This source target distribution chosen here is a Maxwellian, with zero mean velocity, and therefore no momentum transfer, and is written as

S src (f a ) = s k M src (x) Lx 0 M src (x) dx (10)
× 1 √ 2πT src exp - v 2 a 2T src .
The mask M src plotted on FIG.

(1) localizes the source region in the central part of the plasma. It has the same form as the wall mask, see Eq. ( 6). However, the width of the transition region being comparable to the width of the source region, there is no flat top in the source amplitude. The source temperature T src and the source magnitude s k determine the rate of energy transfer and particle transfer. These two control parameters are chosen here to be constant both in space and time and independent of SOL conditions. To ensure energy transfer from the source in the plasma bulk to the sink in the wall, the obvious constraint is to set the temperature of the source to be larger than the temperature of the target distribution function in the wall, see Eq. ( 5). The ratio chosen in the simulation verifies T src /T w > 1. However, to properly describe the distribution functions in the simulations with typically 10 2 points in velocity space, the chosen value for this ratio is not as large as would be consistent with experimental conditions.

Collisions

The collision operator C(f a ) in the simulations of the Boltzmann equation Eq. ( 1) is the sum of the linearized self and inter species collision operators. A detailed description of the chosen linearized operators is found in Appendix A. As will be discussed, it is critical to include a self scattering mechanism in our system so that a steady state can be reached, see also the companion paper [32]. This illustrates the fundamental difference between a collisionless plasma with ν ⋆ 0 = 0 and a weakly collisional plasma, for which ν ⋆ 0 → 0 + . The dimensionless control parameter ν ⋆ 0 is defined in Eq. ( 13). It is similar to the standard ν ⋆ parameter used in Gysela [START_REF] Grandgirard | [END_REF], slightly modified to account for properties of the SOL collisional transport. The collisional regularizing operator in velocity space must handle cases with large departure from Maxwellian distributions. Furthermore, to investigate the possible role of collisions on the sheath physics, it must address a regime characterized by variations on the Debye scale. While the non-linear Landau and Lenard-Balescu collision operators do not prescribe near Maxwellian distributions [START_REF] Nicholson | Introduction to plasma theory[END_REF], their derivation assumes spatial scale separability, i.e. that each distribution function only depends on velocities at the Debye scale. To our knowledge, there is no published closure of the kinetic hierarchy allowing the calculation of a collision operator that would remain valid in conditions with plasma inhomogeneity on the Debye scale. Further collisional issues are met when stepping to realistic SOL conditions where electron-impurity and ion-neutral collisions become particularly important in the high recycling and detached divertor regime. Our aim in the present work is to first address the so-called sheath limited regime in the gyrokinetic simulations of turbulent transport. We thus consider a low-collisionality SOL regime where such complex atomic processes can be ignored. Finally, in the simulations we take into account collisions with two linearized operators, one specific to inter-species electron-ion collisions, and the other for self-collisions, both electron-electron and ion-ion collisions. The latter has the key features to recover the properties of parallel transport, namely that the collision exchange decreases with increasing particle velocity. This is most important in determining the electron heat diffusivity, which is a key point to address the divertor regimes. The chosen expression is devised to satisfy the fundamental conservation properties of elastic binary collisions. The self-collision operator C aa we use is inspired from [START_REF] Dif-Pradalier | [END_REF]. It conserves density, momentum and energy and Maxwellian distribution functions belong to its kernel. A complete description is given in Appendix A. Given the chosen normalization it is written as

C aa (f a ) = ν ⋆ D0 ∂ va D v f Ma ∂ va f a f Ma . (11)
The f Ma function is an artificial distribution constructed to be a Maxwellian with mean velocity u Ma and a ther-mal energy T Ma computed at each time step and at each spatial position to constrain momentum and kinetic energy conservation. The diffusion coefficient in velocity space D v is defined in Eq. (A4). It is proportional to

ν ⋆ aa = √ A a Z 4 a n a /T 3/2
a , Eq. (A5). It 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 control parameter ν * D0 is the reference collision frequency normalized by the plasma frequency. It is further discussed below. The inter species collision operator is adapted from [42] and fully described in Appendix A. The inter-species operator is set to drive the distribution functions towards thermal equipartition and equal mean velocity. With the chosen normalization, it is written as

C ab (f a ) = ν ⋆ D0 2Q ab n a T a H 2 + R ab n a T 1/2 a H 1 F Ma (12a)
In this expression H 1 and H 2 stand for the Hermite polynomials of order 1 and 2 that depend on the normalized velocity

w a = (v -u a )/ √ T a H 2 = 1 2 (w 2 a -1) ; H 1 = w a and consistently F Ma = exp(-1 2 w 2 a )/
√ 2π is the normalized Maxwellian distribution. In these expressions T s and u s are the normalized thermal energy and mean velocity of the species s distribution function f s . It must be reminded here that while the normalization of the thermal energy is identical for all species, that of velocities is species dependent, which leads to √ A s factors when addressing species s. The term proportional to R ab then accounts for momentum exchange between species a and b, hence a friction force, while the term proportional to Q ab stands for the energy transfer, leading to equipartition. These two terms, expressed in normalized units, are

R ab = -n a ν ⋆ ab u a -u b A b /A a , (12b)
Q M ab = -3n a ν ⋆ ab (T a -T b ) A b A b + A a -u a R M ab (12c)
The species dependent contribution to the collision frequency is

ν ⋆ ab = √ 2A a (1 + A b /A a )n b Z 2 a Z 2 b /T 3/2
a , Eq. (A7). The collisionality is set using ν ⋆ 0 as control parameter, such that ν ⋆ 0 = ν ⋆ D0 L ∥ /λ D0 , see Eq. ( 13). It is more convenient than the standard ν * parameter for the SOL physics since the parallel transport process does not depend on trapped particles. As for ν ⋆ the reference parallel scale is the system size in the parallel direction, hence the difference with respect to ν D0 defined with λ D0 as normalizing scale. The control parameter ν ⋆ 0 can be read as 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. The expression of

ν ⋆ 0 is ν ⋆ 0 = 4 √ π 3 e 2 4πε 0 2 n 0 T 2 0 log Λ L ∥ = ν ⋆ D0 L ∥ λ D0 = L ∥ L coll (13
) where log Λ is the Coulomb logarithm. The control parameter ν ⋆ 0 does not depend on particle mass. Collisions in the plasma are a continuous random walk process and thus contribute even for ν ⋆ 0 ≤ 1. Furthermore, their effective amplitude depends on both the control parameter ν ⋆ 0 and on the values taken by the operator. For the self-collision operator Eq. 11, the latter vanishes for Maxwellian distributions and one can expect values larger than unity whenever the velocity space derivative becomes large. It is also important to keep in mind that the expression of ν ⋆ 0 is an average value where the low velocity particle have a strong weight. Supra thermal electrons exhibit lower collisionality, at twice the thermal velocity, Eq. 13 indicates a factor four reduction of collisionality for this class of particles. In a high performance and high power fusion plasma, SOL plasma density and temperature values, n 0 = 2.5 10 19 m -3 and T 0 = 100 eV, are reasonable in the sheath limited regime close to the transition to the high recycling regime. These values lead to λ D0 ≈ 1.5 10 -5 m. One finds therefore that L ∥ = 10 6 λ D0 for L ∥ ≈ 15 m. In such a sheath limited regime, the parallel gradients of the temperature are small and homogeneous conditions are expected along the field line. We have used this consideration to address smaller box sizes and consequently reduce the numerical resources required to perform the simulations. As indicated in Table I, the size of the simulation domain in units of λ D0 is L x = 700 and 2L ∥ = 406. The connection length 2L ∥ = 406 of the simulated magnetic field is very small as compared to the usual connection length of a SOL magnetic field line in a typical tokamak device. In the following we refer to the simulated plasma as a SOL plasma in the sense that (i) it connects two walls that represents limiter/divertor walls, and (ii) that the current model takes into account a particle source term that simulates cross-field transport from the core region. This abusive SOL plasma denomination is taken only for the sake of simplicity. The authors would like to point out that the connection length considered here is not realistic and that they do not model an entire SOL field line. One should keep in mind however that the connection length scan investigated in Section IV C did not demonstrate any dependence of the sheath potential drop on this parameter. This facts hints that the connection length considered in the present work is large enough to study some properties of plasma-wall interaction such as the sheath potential drop. More importantly the width of the transition region between plasma and wall conditions is set to 1/5 of λ D0 . We have checked that the transition is sharp enough that the simulation behavior then hardly depend on this parameter. In particular, it is this parameter that defines the size of the region where the penalization operates. On FIG 16 it can be seen that all of the particle, momentum and energy sink terms S na , S ua and S ha are non-zero only close to the wall position, and vanish in the wall after a few Debye lengths. It is these sink terms that carry the penalization process as they are responsible for absorbing particles, momentum and energy in the wall region. Moreover, the scan performed in Section IV C did not demonstrate a large dependence of the sheath potential drop on the connection length. The normalized target wall density n w is much smaller than unity so that evanescent plasma conditions are enforced in the wall. Regarding the thermal energies, we set the source value at the temperature T 0 so that T src = 1 is the upper value that is achieved when the plasma reaches the source temperature. In the wall, the smaller value T w is chosen. Given a fixed number of points in velocity, one must properly describe both temperature ranges of the distribution function, towards the high velocities at the "hot" end, and towards the small velocities at the "cold" end. It appears that the actual value of the "cold" wall plays a minor role in the penalization process compared to the key role of the wall target density n w . The appropriate resolution determined by the velocity mesh is found to be more important than the actual value of T w provided it is smaller than T src . The choice of collisionality is also made according to a similar constraint on the velocity mesh resolution. Before discussing this point let us first address the collisionality regime for the chosen plasma parameters. Given the definition of the Debye scale λ D0 and that of ν ⋆ 0 , one obtains

C. Plasma regime and simulation control parameters

L coll = λ 2 D0 /λ L0 (3 √ π/ log Λ)
, where the Landau scale λ L0 = (e 2 /(4πε 0 ))/T 0 is the distance of closest approach. It is also interesting to note that λ D0 /L coll = 3 √ π log Λ/Λ given Λ = λ D0 /λ L0 , so that for standard plasmas with Λ ≫ 1, one finds that the physics on a Debye scale is always weakly collisional L coll ≫ λ D0 . In practice, with λ L0 ≈ 10 -11 m and λ D0 ≈ 1.5 10 -5 m one obtains Λ = 1.5 10 6 , and therefore ν ⋆ D0 = λ D0 /L coll ≈ 5 10 -5 , so that for L ∥ = 10 6 λ D0 , ν ⋆ 0 ≈ 50. However, for the reduced box size chosen for the simulation, L ∥ ≈ 200λ D0 , one finds ν ⋆ 0 ≈ 10 -2 . The large value considered in the simulations ν ⋆ 0 = 50 (cf. Table I) is a compromise between realistic physics that require ν ⋆ 0 ≈ 10 -2 and numerical constraints. Indeed, the smaller ν ⋆ 0 , the smaller the grid step ∆v required to resolve steep gradients. The latter tend to generate negative values of the distribution functions. Balancing the time scale for collisional diffusion on the scale ∆v by a time scale that characterizes the convection, one finds that the number of mesh point in the velocity space will scale typically like 1/ ν ⋆ 0 . The cost to use the realistic value of ν ⋆ 0 is therefore a factor 70 increase in the number of points required to mesh the velocity space. The increase of collisionality by a factor 5 × 10 3 in the simulation has several implications. First, as discussed above, it avoids generating negative values of the distribution functions in the regions of strong gradients, as in the sheath and wall interface. Second it prevents the accumulation of slow electrons in the bulk plasma, which would be otherwise confined. Third, it modifies parallel transport such as the collisional heat diffusivity, which is found to scale like 1/ν ⋆ 0 . However, the choice made for ν ⋆ 0 is such that the collisionality in the bulk plasma is unchanged when reducing the system size from 10 6 λ D0 in the reference experimental conditions to 200λ D0 in the simulations. The bulk SOL plasma remains at fixed collisionality despite the reduced size used in the simulations. The collisionality change is important in the sheath where it is increased to ν ⋆ D0 = 0.25 in the simulations. Convergence tests have been performed, showing that the main results presented in this paper do not depend on the precise value of ν ⋆ 0 .

III. FLUID AND KINETIC PROPERTIES OF PLASMA-WALL SELF-ORGANIZATION

A. On the self-consistent definition of the sheath entrance and wall position

To be able to compute the potential drop in the sheath we need a systematic definition of the sheath entrance and the wall position. It turns out that one cannot derive a general expression of the sound speed in a kinetic plasma that would remain valid in every collisionality regime. Indeed, as shown in the companion paper [32], such an expression critically depends on the fluid closure that is retained: the definition of the sound speed is then model dependent. Fortunately, when considering divertor plasmas under realistic tokamak conditions, the polytropic closure generally yields a well-defined sound speed compatible with kinetic simulations [43]. In the companion paper we additionally observe that the kinetic equivalent of Bohm criterion-the so-called Riemann or Harrisson-Thomson criterion-fails when collisions are accounted for in the sheath region, and when the sheath is no longer considered as being an asymptotic limit for the plasma. Both these hypotheses that are necessary for the kinetic criterion to apply are not satisfied in the present work. Indeed, we include a Fokker-Plank operator to model collisional processes in the sheath region, see Section II B 2, and the Debye scale is resolved in the present simulations. Therefore, slow ions for which v ≈ 0 can be found at all positions in the simulated plasma, thus the velocity integral of the Harrisson-Thomson criterion-which is roughly of the form f i /v 2 -is not defined. In the companion paper [32] an intuitive definition is proposed instead, where the sheath entrance corresponds to the location where the derivative of the charge density ρ c with respect to the potential dρ c /dϕ transitions from a quasi-linear variation in the presheath to a rapid exponential growth in the sheath region. This is the definition we use in the present work. This definition is somewhat arbitrary, but since the Debye length is resolved in the present work, there is by definition no abrupt transition between a quasineutral plasma and a positively charged sheath region. Moreover, the results presented in this work mainly focus on the potential drop in the sheath region, which is not very sensitive to the precise choice of the definition of the entrance of the sheath. Indeed, as can be seen on FIG. 4, the variation of the potential at the entrance of the sheath is small as compared to its variation inside the positively charged region. Note that this definition cannot be used to define the entrance of the sheath when the wall is tilted with respect to the magnetic field line. In this situation the entrance of the sheath should be computed using the so-called Bohm-Chodura condition [24]. The definition taken in the present work is suitable only for our specific case of a magnetic field that impacts the wall with a perpendicular incidence. For the definition of the wall position x w we use the location where the charge density is maximum. In the classical sheath model the electron and ion densities diverge from each other in an increasing monotonous way from the sheath entrance to the wall position. We see on FIG.

(3) that this maximum corresponds approximately to the location where the wall mask starts its transition from zero to one. This position corresponds as well to the beginning of the region where particles are absorbed by the sink term defined in Eq. ( 4). Also, it appears that the energy of the plasma is conserved up to this position, as discussed later in Section III C. Indeed at this location the normalized total current j tot starts to deviate from zero, meaning that the sink term begins absorbing particles. A rigorous mathematical formulation for the sheath entrance and wall position is given below. We consider only the right half of the plasma defined by x ∈ [L x /2, L x ], so that the wall is located at x ≥ x w . The definition of x sh and x w for the left half of the plasma is similar since the simulation box is symmetric. The wall position is defined through the following implicit relationship

ρ c (x w ) = max x≥Lx/2 ρ c
The sheath entrance location is defined as the intersection between the two linear fits defined as

D 1 : L log dρ c dϕ with x ∈ [L x /2, x w ] and α dρ c dϕ ∞ ≤ dρ c dϕ ≤ dρ c dϕ ∞ ,
where D 1 is the label for the linear fit, L is a linear fit operator and α ∈ [0, 1] is an arbitrary coefficient that sets the sensitivity of the sheath entrance definition with the charge density fluctuation. The higher α, the deeper in the positively charged region the entrance of the sheath will be. Typically we take α = 0.1. The ∥ •∥ ∞ notation refers to the maximum of the considered quantity taken on the [L x /2, x w ] interval. The second linear fit is defined as

D 2 : L log dρ c dϕ with x ∈ [L x /2, x w ] and 0 ≤ dρ c dϕ ≤ α dρ c dϕ ∞ .
The first linear fit D 1 corresponds to performing a linear fit of log |dρ c /dϕ| in the positively charged region, whereas the second linear fit D 2 is performed in the quasineutral (presheath) region. Note that if the definition of the sheath entrance depends on the chosen value for the α parameter, the potential drop in the sheath does not depend much on this choice since the variation of the potential at the sheath entrance is small as compared to its variation within the sheath. For a graphical illustration of this linear fit process see the companion paper [32]. The sheath entrance and wall position defined above can be viewed on FIG. (4).

B. Distribution functions

We define the electric potential drop in the sheath region as ∆ϕ sh = ϕ(x w ) -ϕ(x sh ), with the sheath entrance and wall position defined in the previous section. This negative potential drop confines the plasma electrons that have a velocity at the sheath entrance smaller than a critical velocity v c . Using the conservation of energy in the sheath for electrons, we obtain for this critical velocity the expression in normalized units

v c = -2∆ϕ sh . (14)
On the top left and bottom left graphs of FIG.

(5) are shown a plot of the ion and electron distribution functions in the entire simulation box, against the (x, v) variables. Both distribution functions tend to vanish in the two wall regions: the penalized sink of particles from Eq. ( 4) absorbs particles there. The iso-contours of the Hamiltonian for both species are superimposed on each graph as red lines. In our case the Hamiltonian for species a is written in dimensional units as

H a = 1 2 m a v 2 + e a ϕ.
The iso-contours of the Hamiltonian gives the trajectories of particles in phase space in the absence of collisions, source and sink terms. For positive velocities, corresponding to the top half of each plot, the iso-contours are to be followed from left to right. For negative velocities (bottom half of each plot), it is the opposite. The two graphs on the right of FIG.

(5) show a zoom of phase space in the vicinity of the sheath near the right wall. The following discussion focuses on the physics of the right part of the plasma, for which x ≥ 0. The x < 0 region is the symmetrical of the x ≥ 0 region. Ions experience an acceleration in both plasma, sheath and wall regions.

The acceleration in the quasineutral plasma is negligible as compared to the one in the sheath and the wall, where the iso-contours bend towards large velocities. This is a consequence of the steep drop of the electric potential that exists there. Electrons on the other hand are slowed down by this potential drop. The passing electrons for which v ≥ v c at the sheath entrance overcome the potential barrier of the sheath and penetrate the absorbing wall region. On the contrary, electrons with v ≤ v c are reflected in the sheath and thus are "trapped" in the plasma. We see that some of the passing electrons are able to come back in the plasma region. This effect is reminiscent of two distinct causes. First, the penalized sink term from Eq. ( 4) does not absorb perfectly every electron that enters the wall region, as the magnitude ν w of this operator is not infinite, and because the transition between the plasma and the wall region is not perfectly abrupt given that the wall mask defined in Eq. ( 6) has a finite steepness. Second, the presence of fast electrons at the sheath entrance coming back to the plasma is also the consequence of collisional processes taking place in the sheath region. Indeed, as the overall effect of plasmawall interaction on the electron distribution function is to absorb fast electrons, the distribution function of this species departs from a Maxwellian. Collisions thus act on electrons to bring them back to thermodynamical equilibrium, that is to say to a Maxwellian distribution function. This effect can be viewed on the ion distribution function as well, since as all ions that reach the wall are absorbed in this region, then collisions are the only mechanism that allows for explaining the presence of ions with negative velocities at the entrance of the sheath. The impact of electron re-emission on the potential drop measured in the sheath region is discussed in detail in Section IV A. We see on the ion distribution graph that the iso-contours of the Hamiltonian alone are not enough to explain all the trajectories of the particles in phase space. In particular, in the area denoted by the dashed red line near the (1) marker we see that the ion distribution function presents an increasing number of particles with negative velocities as we go towards the plasma center. This observation is not consistent with the trajectories depicted by the Hamiltonian iso-contours, but is a direct consequence of self species collisions. Collisional processes tend to replenish this region of phase space to make the distribution function closer to a Maxwellian.

The trapped electrons for which v ≤ v c do not contribute to the electron flux at the sheath entrance. To illustrate this, on FIG. (6) we plot the partial electron flux Γ ⋆ e against the velocity variable. This partial flux is defined as

Γ ⋆ e (v) = 1 Γ sh e v -v dv ′ v ′ f sh e (v ′ ), (15)
where Γ sh e is the total electron flux at the sheath entrance. The partial electron flux is negligible for velocities lower than the critical velocity v ≤ v c and starts to become non zero when v ≥ v c . The electron flux at the sheath entrance is therefore essentially carried by the high velocity |v| > v c component of the electron distribution function.

On the same graph we also plot the negative velocity part of the electron distribution function at the sheath entrance f sh- e (|v|) ≡ f sh e (-|v|) on top of its positive velocity part f sh+ e (|v|) ≡ f sh e (|v|). There exists a clear imbalance between these two parts of the electron distribution function for velocities higher than the critical velocity v c . The majority of passing electrons directed towards the wall are absorbed there. We still find some fast electrons at the sheath entrance that are directed towards the plasma, i.e. f sh- e (|v| ≤ -v c ) ̸ = 0. Their origin is twofold. First, the wall is not perfectly absorbing: the magnitude ν w of the particle sink term of Eq. ( 4) is not infinite. Therefore, some of the passing electrons at the sheath entrance are able to penetrate the wall region and come back into the plasma. Second, the self species collision operator tends to make the distribution function converge towards a Maxwellian, replenishing the depleted v ≤ -v c velocities in the sheath region. In summary, this situation departs from the academic case of a collisionless plasma interacting with a perfectly absorbing wall, where f sh- e (|v| ≤ -v c ) would be strictly zero. The ion distribution function is shown in FIG. [START_REF] Bohm | The Characteristics of Electrical Discharges in Magnetic Fields (ch 3[END_REF]. It is shifted towards positive velocities. Let us notice that the underlying mechanism that results in a positive particle flux oriented towards the wall differs between the two species. The electron flux is essentially carried by the imbalance of the electron distribution function at high velocities. Conversely, the ion flux is mostly due to a global shift of the distribution function towards positive velocities. We further notice that the ion distribution function is depleted for velocities v ≤ 0. All the ions present at the sheath entrance have been accelerated in the quasineutral plasma, therefore they have acquired a positive velocity upshift. The potential drop of the sheath further accelerates the ions towards the wall where they are all absorbed. The remaining ions around v ≈ 0 originate once again from self species collisions that tend to replenish the low velocity part of the ion distribution function. There is in some sense a competition between this self scattering process and the acceleration of ions driven by the electric field which depletes the slow ions part of the distribution function. On FIG. [START_REF] Bohm | The Characteristics of Electrical Discharges in Magnetic Fields (ch 3[END_REF] we see that the electron distribution function appears to be composed of two truncated Maxwellian. The first one contributes to the overall distribution function at low velocities, whereas the second one, which has a higher temperature and lower density compared to the first one, is visible at larger velocities. This effect is better seen for the positive velocity part of the electron distribution function which has not been affected by the wall sink. This feature is essentially due to intra-species collisions. The central Maxwellian part of the electron distribution has the same characteristics as the kernel distribution function f Me of the intra species collision operator, see Eq. ( 11) and Appendix B. On the other hand, the high velocity part of the electrons and ions functions share a form which is close to the source distribution function, FIG. 8. Sheath heat transmission factor for ions (red triangles) and electrons (blue circles). The total (ions + electrons) sheath heat transmission factor is represented in black. The predictions for each of these quantities are written with the sh superscripts to indicate that they are computed at the sheath entrance.

plotted as the blue dotted curve on FIG. [START_REF] Bohm | The Characteristics of Electrical Discharges in Magnetic Fields (ch 3[END_REF]. The trajectory of these high velocities particles have not been much altered by the self scattering processes: the magnitude of the intra species collision operator decreases like |v| -3 at large velocities. The effect of the collisionality on the shape of the distribution functions is investigated more extensively in Appendix B.

C. Sheath heat transmission coefficient

One of the key outcomes of plasma-wall interaction studies is the estimation of the heat flux deposited by the plasma on its material boundaries. This is especially critical for the Iter project, where the power deposited on the divertor target must stay around 10 MW/m 2 in steady state [44]. The so-called sheath heat transmission coefficient γ is defined by the relation

Q = γΓT e ,
where Γ, Q and T e are respectively the particle flux, energy flux and electron temperature taken at the plasmawall boundary. In fusion devices the electron temperature along with the particle flux can be experimentally measured using Langmuir probes. Knowing the value of the sheath heat transmission coefficient thus allows the computation of the incident energy flux on the plasma facing components. Another way of defining the sheath heat transmission factor consists in separating the electron and ion energy fluxes and thus writing

Q e = γ e ΓT e , and
Q i = γ i ΓT i .
These relations prove especially usefull in fluid codes, where the value of the electron and ion temperatures along with the particle flux at the plasma-wall interface give access to boundary condition values for the energy flux. Using the results from the reference simulation we can compute the value of the heat transmission factors presented above, with the energy flux of species a defined as

Q a = ∞ -∞ dv 1 2 m a v 3 f a . (16)
Note that in our case the energy flux, along with the particle fluxes and the temperatures are functions of the position and time variables, and thus the same is true for the heat transmission factors. On FIG. 8 are shown the values for the heat transmission factors (electrons, ions and total) between the entrance of the sheath and the wall location. Let us first remark that the total heat factor γ is approximately constant in the whole sheath region. This is a direct consequence of the conservation of energy in the sheath-and further indication of the consistence for the choice of the wall location, see Section III A. At first order neither energy nor particle sources exist in the sheath, therefore the total energy flux is constant. The ion heat factor slightly increases in the far sheath while the electron factor decreases. As expected the sheath enables energy transfer between species: fast electrons slow down in the sheath, their energy is transferred to ions that therefore accelerate.

In Appendix E we detail a model developed in [START_REF] Stangeby | The Plasma Boundary of Magnetic Fusion Devices[END_REF] which is commonly used to predict a value for the sheath heat transmission coefficient. The prediction yields γ pred. e = 4 for electrons, and γ pred. i = 2.5 for ions. As explained in Appendix E, when considering a 3D geometry the predictions are different: γ 3D,pred. e = 5 and γ 3D,pred. i = 3.5. As shown on FIG. 8, the predicted heat factor value for electrons and ions are close to the one found in the simulation. The ion sheath heat transmission factor is approximately 60% higher than the predicted value. Thus, the total heat transmission factor obtained from the kinetic simulation is approximately 23% higher than the prediction. Lastly, the prediction suggests that ions deposit less power on the wall than electrons: γ pred. i = 2.5 < γ pred. e = 4. From the kinetic simulations we find that ions deposit 8% more power on the wall as compared to electrons. The discrepancy between the ion heat transmission factor and the prediction can be explained from kinetic effects. The prediction assumes the ion distribution function to be a Maxwellian shifted towards positive velocities. As shown on FIG. 7 we see that the actual ion distribution function is not Maxwellian. In particular, most of the ions at sheath entrance have velocities higher than the mean velocity of the ion distribution function. Said differently, the median of f sh i is higher than its mean value. Therefore, the ion heat flux found in the simulation is higher than its counterpart computed from a Maxwellian distribution.

IV. INFLUENCE OF SOME RELEVANT PLASMA PARAMETERS ON THE POTENTIAL DROP IN THE SHEATH REGION

A prediction for the electric potential drop in the sheath region can be found in [START_REF] Stangeby | The Plasma Boundary of Magnetic Fusion Devices[END_REF]. It is written in dimensional units as

∆ϕ pred. sh = T sh e e log   2π
m e m i M 2 1 + T sh i T sh e   , (17)
where the ion and electron temperatures are supposed to be taken at the sheath entrance, as indicated by the sh superscript. M is the Mach number defined as M = u i /c s where c s is the velocity of plasma sound waves. When considering an isothermal closure of the fluid equations, namely the conservation of density, momentum and energy this sound speed can be written as c s = (T i + T e )/m i . In the companion paper [32] we show that this expression does not necessarily match the velocity of sound waves of our kinetic simulation. However, the value of the velocity of sound waves does not matter for the following discussion. Indeed, from Eq. (F3) and by equating the electron and the ion fluxes at sheath entrance we see that the quantity the quantity inside the logarithm of Eq. ( 17) is simply equal to the ratio between the ion flux and the forward flux of electrons computed at the entrance of the sheath, that is to say ∆ϕ pred. Furthermore, the quantity Γ sh i /Γ sh+ e is independent of the choice of the definition for the speed of sound waves, as we have Γ sh i = n sh i M c s whatever the definition of c s . Thus, another choice for the sound wave velocity would lead to another expression for the potential drop prediction that would yield the same numerical value when computed. The derivation of the sheath potential prediction is well known and is is recalled in Appendix F 1. The prediction assumes a perfectly absorbing wall that does not emit any particles back into the plasma. From this expression we can construct an effective potential drop as FIG. 9. Sheath potential drop dependence on the value of the ion-to-electron mass ratio. The red squares and blue circles correspond to conducting and insulating walls respectively. The effective potential drop defined in Eq. ( 18) and the refined potential drop defined in Eq. ( 22) are both represented. The dashed grey curve corresponds to the prediction, which is the same in both effective and refined cases.

Looking at Eq. ( 17) we see that the predicted value for this effective potential drop only depends on the mass ratio m e /m i as

∆ϕ pred. sh = 1 2 log m e m i . (19)
A. Potential drop dependence on the ion-toelectron mass ratio value

The mass ratio between ions and electrons m i /m e is fixed by the respective electron and chosen ion masses. In global gyrokinetic codes it is common practice to use a different ratio, corresponding to heavy electrons, or equivalently lighter ions. In Gysela for instance, this ratio is often set to m i /m e ≈ 100. This allows for resolving the ion and electron timescales at a numerical cost that remains affordable. The main constraint lies in making sure that the physics observed with an adapted mass ratio value is similar to the one that would be obtained with a realistic value. In the perspective of an integration in Gysela of the plasma-wall interaction physics, we assess the dependence of the potential drop on the mass ratio between ions and electrons. We investigate the effect of mass ratios in the range {40, 100, 400, 1835, 3600} on the sheath potential drop. The m i /m e = 1836 value corresponds to a hydrogen plasma, while m i /m e = 3600 is similar to that of a deuterium plasma. Two types of simulations have been performed, either with an insulating or conducting wall -see Section II. As explained in Section III B the present penalized model leads to observing electrons at the entrance of the sheath with velocities such that they should have been absorbed in the wall. This effects stems from two different causes. First, the wall in our simulations is not perfectly absorbing and therefore some electrons that enter the wall region are able to come back into the plasma. Second, collisions taking place in the sheath region tend to replenish the part of the electron distribution function that corresponds to electrons that should have been absorbed in the wall region. The resulting flux from these fast electrons is not taken into account in the classical prediction for the potential drop in the sheath given by Eq. ( 17), which assumes a perfectly absorbing wall. We therefore adapt the prediction of the potential drop given by Eq. ( 17) to account for this electron re-emission effect. The derivation is detailed in Appendix F 2. This leads to the following formula for the refined potential drop prediction

∆ϕ ref. sh = T sh e e log 2π m e m i M 2 1 + T sh i T sh e × 1 + Γ - e Γ sh i , (20)
Where Γ - e is the flux at the sheath entrance of high velocity electrons coming back from the wall

Γ - e = - -vc -∞ dv vf sh e (v). (21)
Using this refined prediction we can define in a similar fashion as in Eq. ( 18) an effective potential drop as

∆ϕ ref. sh ≡ e∆ϕ sh T sh e -log 2πM 2 1 + T sh i T sh e × 1 + Γ - e Γ sh i . (22)
This refined potential drop is represented on FIG. [START_REF] Riemann | [END_REF] using plain squares for the conducting wall, and plain circles for the insulating wall. We see a good agreement between the refined potential drop and the prediction, which is given by Eq. ( 19). The potential drop increases in absolute value with the mass ratio, as expected. The higher the inertia difference between electrons and ions, the easier it is for the light species to enter the wall region as compared to the heavy one. Hence, the potential drop required to ensure ambipolarity at the plasma wall transition is larger-i.e. more negative. The effective potential drop without refinement given by Eq. ( 18) that does not take into account the electron flux coming back from the wall is represented as well on the same figure as the baseline scenario, using empty squares and circles for the conducting and insulating walls respectively. We see in this case a discrepancy between the computed potential drop and the predicted value that is larger than when taking into account the electron flux coming back from the wall. We observe additionally that this discrepancy increases when approaching realistic ion-to-electron mass ratios, for both insulating and conducting walls: this electron flux becomes essential when the electrons becomes more mobile. We also notice that the discrepancy is more significant in the insulating wall case as compared to the conducting case. In the insulating wall there exists much more electrons than in the conducting wall, as shown by the charge densities profiles of FIG. (4). Therefore, the wall emits more electrons back into the plasma in the insulating case, departing further from the model of a perfectly absorbing wall. In the present section we therefore have shown that when varying the electron-to-ion mass ratio and accounting for relevant kinetic effects, the dependence of the potential drop in the sheath region obtained from the simulation is compatible with the simple sheath analytical model presented in Appendix F. Namely, the main discrepancies between the analytical model and the kinetic simulations lies in the re-emission of electrons in the plasma that we observe as a consequence of the (i) the finite strength of the penalization term in the wall region and (ii) the presence of collisions in the sheath region. It is therefore reassuring from the perspective of model reduction-e.g. from the perspective of incorporating Debye (unmagnetised) sheath physics in gyrokinetics-to observe this robustness of the potential drop with reduced mass ratios in the unmagnetised sheath problem.

B. Potential drop dependence on the electron temperature

It is well admitted that the establishment of the Hmode transport barrier at the last closed flux surface in tokamak plasmas is closely linked to the shearing of the electric field at this location [45,46]. The electrostatic potential radial variation is governed by turbulent transport and neoclassical physics in the core plasma, whereas it is mainly driven by plasma-wall interaction in the Scrape-Off layer. In particular we expect the potential drop in this region to increase linearly in magnitude with the electron temperature. This is quite intuitive: if electrons have a large velocity on average at the sheath entrance, i.e. if the electron temperature is large, then the potential barrier needs to be larger to confine a fraction of these electrons in the plasma. This prediction is usually written as ∆ϕ pred. sh = ΛT e /e. Taking into account the electron flux coming back from the wall at the sheath entrance the Λ term is given by Eq. ( 20) as

Λ = log   2π m e m i M 2 1 + T sh i T sh e 1 + Γ - e Γ sh i   . ( 23
)
Note that this Λ term is usually considered constant because of the relatively slow variation of the logarithm with its argument. In the following we assess the dependence of the potential drop in the sheath with the electron temperature. To vary the electron temperature at the sheath entrance we perform several simulations with a different value for the electron temperature of the particle source, which is defined in Eq. ( 10). We use temperatures in the range T src,e ∈ {0.5, 0.75, 1, 1.25, 1.5}. The temperature of the ion source remains constant in all of these simulations, we use a mass ratio value of m i /m e = 400 and a plasma length of L ∥ = 300λ D0 . FIG. (10) shows the quantity e∆ϕ sh /Λ for both insulating and conducting wall cases, with the Λ term defined in Eq. ( 23). The dataset is referred to on the graph using the label refined. We observe the expected linear dependence of the normalized potential drop e∆ϕ sh /Λ on the electron temperature at the sheath entrance. The slope of the plotted data is in good agreement with the predicted slope of 1. A constant offset exists between the values obtained from the simulation and the prediction. It can be mainly explained by the departure of the electron distribution function from a Maxwellian at the sheath entrance. In Appendix F 1 to derive the prediction of Eq. ( 20) we write the equality of ion and electrons fluxes at the sheath entrance as

Γ sh i = Γ +pred.
Note here that f sh e is the electron distribution function at the sheath entrance obtained from the code output, i.e. it is a simulation-dependent quantity. In this framework the prediction for the potential drop is correct only if the predicted flux matches the actual forward fluxi.e. if Γ + e = Γ +pred. e . As shown on FIG. [START_REF] Bohm | The Characteristics of Electrical Discharges in Magnetic Fields (ch 3[END_REF] the electron distribution function is not a Maxwellian for velocities higher than v c but resembles more to a superposition of two Maxwellians of different temperatures, as discussed in detail in Section III B. Therefore, the forward flux of the simulations departs from the predicted forward flux: Γ + e ̸ = Γ +pred.

e

. This departure can be quantified by noticing that f sh e (v ≥ v c ) is approximately the sum of two Maxwellians of different temperatures and densities. We thus adapt the temperature used in the predicted flux expression of Eq. ( 24) so that it matches the forward flux of the simulation, i.e. so that Γ + e = Γ +pred.

e

. This corresponds to taking the temperature of the best Maxwellian fit for the v ≥ v c part of the electron distribution function at the sheath entrance. As discussed in Section III B, this temperature is close to the one of the particle source term S src -it would be equal to it in the absence of collisions. By doing this we obtain the appropriate set of points on FIG. (10), which does not show any offset as compared to the refined case. This adaptation of the temperature of the predicted forward flux further highlights one of the limitations of the analytic model we use. In particular, the shape of the distribution function we obtain results from purely kinetic effects, and of collisions in particular. Using fluid moments such as the temperature, obtained by integrating over the entire velocity space, does not allow one to explain the details of the plasma-wall interaction physics we observe.

C. Potential drop dependence on the parallel length

The L ∥ term represents the length of the field line separating the right and left wall. If this length is already large enough we do not expect that further increasing L ∥ will have a significant effect on the plasma-wall interaction physics. The prediction of Eq. (17) shows no dependence of the potential drop on L ∥ . Indeed, in the classical sheath theory the presheath is supposed to be long enough so that the small electric field that develops there is able to accelerate ions up to at least sonic velocities to ensure that ambipolarity is reached at the sheath entrance. In the following we investigate parallel length in the range L ∥ ∈ {300, 600, 1200, 2400}. The width of the walls does not change between these simulations: the fraction of the total length of the simulation domain L x FIG. 11. Effective potential drop dependence on magnetic field line length L ∥ , for conducting (red squares) and insulating (blue circles) walls. The effective potential drop definition is the one used in the mass ratio scan, i.e. Eq. (18).

that accounts for the wall region decreases when L ∥ increases. The particle source is kept constant in these simulations. Also, the collisionality of the plasma ν ⋆ 0 defined in Eq.13 has been taken constant in the simulations. On FIG. (11) we plot the effective potential drop defined in Eq. ( 18) both conductive and insulating walls. The same scale is used as in the mass ratio scan of FIG. [START_REF] Riemann | [END_REF]. We see that the values from the simulations stay close to the constant predicted value of potential drop, with a weak dependence in L ∥ . When using the refined potential drop that is adapted for the electron temperature as explained in Section IV B we see almost no variation of this quantity with the parallel length.

V. DISCUSSION: SHEATH-LIKE BOUNDARY CONDITIONS WITHOUT RESOLVING THE SHEATH

Simulations of plasma-wall interaction implies dealing with two kinds of boundary conditions. On the one hand the sheath can be viewed as a boundary condition for the plasma since it forms at the plasma boundary and influences the plasma behavior. On the other hand a mathematical boundary condition is also needed to solve numerically the considered model. The immersed boundary setting also allows one to investigate periodic simulation domains, in such a case the mean potential is not constrained and can be fixed independently. In our case for instance, we need to specify some explicit boundary conditions (Dirichlet, Neumann...) to solve the Boltzmann-Poisson system (1, 2). There exist different ways of dealing with these two boundary conditions. A first approach consists in specifying the sheath boundary condition and the numerical boundary condition at the same spatial location in the system. Examples of this approach are the logical and insulating sheath techniques, which are detailed below in Sections V 1, V 2. A second solution consists in decoupling spatially these two boundary conditions. This is the immersed boundary conditions approach, detailed in Section II. Interestingly, such a geometry is reminiscent of that met when the plasma interacts with gas as contemplated in the detached divertor regime. Different specification of the wall properties would then allow one to address this particular regime. For the more standard case of a hot plasma interacting with a dense material, characterized by prompt recombination, it will be shown how such treatment differs from the logical and conducting sheath approaches, and how changes in control parameters allow one to explore nevertheless the physics of both no-current and conducting sheaths. In the following, we briefly review the logical and conducting sheath models, before introducing the penalization technique that is the backbone of our approach.

Logical sheath

The logical sheath boundary condition has first been introduced in [47]. In this approach, the quasineutral plasma only is modeled. The positively charged region that would form at the plasma boundary is not resolved, see FIG. [START_REF] Boyd | Proceedings of the Physical Society. Section B[END_REF]. The logical sheath condition allows one to include the effects of the sheath on the plasma anyway. This condition ensures ambipolarity at the plasma boundary most of the time, see below. It therefore mimics the sheath equilibrium physics. The logical sheath also plays the role of a numerical boundary condition since it proposes a procedure to deal with outgoing particles at the simulation boundary. This condition is applied by the authors to a simulation that uses a particlein-cell method. This procedure can however be adapted to other numerical methods, such as a semi-Lagrangian simulation for instance, see below. The authors simulate a one-dimensional electrostatic and unmagnetized plasma (or equivalently the parallel dynamics of a magnetized plasma with normal incidence). Let us consider that the plasma spans on the -L ∥ /2, L ∥ /2 region. The z = ±L ∥ /2 coordinate denotes the plasma boundary. In a particle-in-cell code, a pusher method moves the particles from time t to t + ∆t. ∆t stands for the simulation time step. The logical sheath condition consists in the following sequential steps, see FIG. [START_REF] Boyd | Proceedings of the Physical Society. Section B[END_REF] as well:

1. Starting from the particle positions at time t, apply the pusher method to compute their new position at time t + ∆t. Let us remark that situation 3. is the most common one since electrons are more mobile than ions. In this case the net current to the wall is exactly zero. Therefore ambipolarity is ensured by the logical sheath condition. Situation 4. rarely happens, and may result from fluctuations in the system. When this case arises, at next time step the wall is positively charged. Therefore, it repels ions and attracts electrons, making situation 3. more probable again. Apart from this particular case the wall bears no charge. The plasma-wall interaction effect on the system arises only as a result of the logical sheath algorithm detailed above. No potential drop develops at the plasma boundary. Therefore, it is not needed to resolve the sheath typical time scales and spatial extents, allowing for larger ∆t and coarser spatial discretization in the simulation. The sheath potential drop ∆ϕ sh is never required in the course of the numerical treatment, but it can be retrieved during the post-processing phase of the simulation. This virtual potential drop is such that it confines in the plasma all electrons crossing the wall boundary that are slower than the cutoff velocity v c . Using conservation of energy, it is therefore expressed in dimensional units as

∆ϕ sh = - m e 2e v 2 c .
The logical sheath algorithm has been described in the framework of a particle-in-cell method, but can be adapted to continuum 5D gyrokinetic simulations [48,49] as well. In this particular case, let us consider that the z variable denotes the position along the considered magnetic field line. The other (x, y) spatial variables refer to the perpendicular plane. Let us finally introduce the parallel velocity v and the magnetic moment µ. Then, the equality of ion and electrons fluxes at the plasma boundary z = L ∥ /2 reads as follows

∞ 0 dv vf i (x, y, L ∥ /2, v, µ) = ∞ vc dv vf e (x, y, L ∥ /2, v, µ).
This equality determines a cutoff velocity v c when situation 4. is not encountered, see above. At the plasma boundary position z = L ∥ /2, all outgoing ions are removed from the simulation at the next time step. The fastest electrons that have velocities larger than v c are removed as well, while the slowest v < v c electrons are reflected back into the plasma. An extra step involving the rescaling of the distribution function is required in this case, since the cutoff velocity v c is not a velocity grid point in general. This process is described in detail in the conducting sheath framework, see below. A drawback of the logical sheath is that it tends to enforce a zero electrical current at the boundary at each time and boundary mesh point. Such a constraint fulfills the physical requirement that the plasma charge does not build-up and consequently that the plasma remains quasineutral but for the sheath boundary layer, it appears to be too stiff. Indeed, the global constraint of zero electrical current does not require zero electrical current at all times and all points. The situation appears to be over constrained when considering experimental evidence. It can therefore modify the turbulence behavior compared to cases where the electrical current locally fluctuates around zero net current [50].

Conducting sheath

Another way of including the plasma-wall physics at the plasma boundary without resolving the sheath scale is presented in [51]. We call this boundary condition conducting sheath in the sense that it allows finite currents to flow freely in and out of the wall, possibly forming current loops when dealing with a 2D plasma boundary. This approach differs from the logical sheath model where the electron flux is adapted to match exactly the ion flux at the wall-at least when situation 3. applies, see above. In the conducting sheath framework the potential at the plasma boundary ϕ sh determines a proportion of reflected electrons. In the meantime, outgoing ions leave the simulation freely. The authors of [51] apply this boundary condition to a 5D gyrokinetic simulation. Let us consider a three dimensional plasma bounded by a wall surface. We use the x, y, z coordinates and the plasma boundary is located at z = ±L ∥ /2. The v variable still refers to the velocity along the z axis. Let us focus on the upper plasma boundary located at z = L ∥ /2. The conducting sheath algorithm consists in the following:

1. The potential at the plasma boundary ϕ sh (x, y) = ϕ(x, y, L ∥ /2) is computed from the quasineutrality equation. This equation can be written as presented in [START_REF] Grandgirard | [END_REF] as

-∇ ⊥ m i n i eB 2 ∇ ⊥ ϕ + e 2 ⟨n e ⟩ FS ⟨T e ⟩ FS (ϕ -⟨ϕ⟩ FS ) = σ g , (26)
where n i (resp. n e ) is the ion (resp. electron) charge density and T e stands for the electron temperature. The bracket notation ⟨. . .⟩ FS stands for the average over a given magnetic surface. The guiding center charge density is defined as σ g = e(Z i n Gi -n e ), where n Gi is the ion gyrocenter density and Z i is the normalized ion charge (equal to +1 for hydrogen isotopes). ∇ ⊥ refers to the gradient in the plane perpendicular to the magnetic field.

2. When the potential ϕ sh (x, y) is positive it defines a cutoff velocity which is such that 1 2 m e v 2 c = eϕ sh . At the considered wall position, all slow electrons that have velocities v such that 0 ≤ v ≤ v c are reflected back into the plasma. The fastest electrons that are not reflected, along with all outgoing ions leave the simulation domain at next time step.

To clarify the reflection procedure let us introduce the velocity grid in the parallel direction (v i ). It is characterized by a discretization step ∆v. The value of the distribution function f e (x, y, z, v i , µ) accounts for the density of particles that have velocities spanning the [v i -∆v/2, v i + ∆v/2] interval. In general the cutoff velocity v c does not lie on a specific point of the velocity grid, see FIG. [START_REF] Persson | [END_REF]. Rather, there exists an index j which is such that the cutoff velocity is inside the v j cell. That is to say we have v c ∈ [v j -∆v/2, v j + ∆v/2]. Therefore, for all velocities 0 ≤ v i < v j the value of the distribution function f e (v i ) is copied to a ghost cell v ′ i which is such that v ′ i = -v i . This ghost cell stores the values of the reflected distribution function and is used as a boundary condition in the advection step, see below. The value of the distribution function at the particular cell f e (x, y, L ∥ /2, v j , µ) is copied to the corresponding ghost cell v ′ j = -v j and scaled by a factor c defined as

c(x, y) = vc vj -∆v/2 dv µmax 0 dµ f e (x, y, L ∥ /2, v, µ) vj +∆v/2 vj -∆v/2 dv µmax 0 dµ f e (x, y, L ∥ /2, v , µ)
.

It ensures that all outgoing electrons verifying 1 2 m e v 2 ≤ eϕ sh are indeed reflected to the plasma, independently of the velocity grid in the z direction. To illustrate the need for a ghost cell to store the values of the reflected electrons we take the case of a gyrokinetic code that uses a splitting scheme along with a forward Eulerian method to solve the advection steps of the gyrokinetic equation. In this framework, a time step can be decomposed into the following steps:

1. Solve the quasineutrality equation.

2. Apply the conducting sheath method.

3. Solve the advection steps, and all the other operators of the gyrokinetic equation.

In the advection step, some boundary conditions are needed at the wall location for particles with negative velocities. That is to say, we have to find a value for f a (x, y, L ∥ /2, v i , µ) for all v i ≤ 0. Since no ion comes back from the wall, we simply set

f i (x, y, L ∥ /2, v i , µ) = 0 ∀v i ≤ 0.
For the electrons we use the ghost cell values as a boundary condition, i.e.

f e (x, y, L ∥ /2, v i , µ) = f e (x, y, L ∥ /2, v ′ i , µ) ∀v i ≤ 0.
The 1 2 m e v 2 c = eϕ sh equation defines a cutoff velocity for electrons only if the potential at the wall position ϕ sh is positive. This is usually the case. Consider the first time step of the simulation: we start from a potential and a charge density of guiding centers σ g equal to zero everywhere. Since ϕ sh = 0 there is no reflection at the plasma boundary: all the outgoing electrons leave the simulation domain. Because of the inertia difference between ions and electrons, the latter move faster than the former. Thus, more electrons than ions are lost in this first advection step. Therefore, at the wall boundary we shall have a positive charge of guiding center σ g (L ∥ /2) > 0. In the quasineutrality equation [START_REF] Moritz | [END_REF], at first order we have -∇ ⊥ m i n i /(eB 2 )∇ ⊥ ϕ ≪ ϕ. This assumption is especially true when starting from an initial condition homogeneous in the perpendicular plane. Therefore the quasineutrality equation gives at the wall position (ϕ -⟨ϕ⟩ FS ) ≈ σ g ⟨T e ⟩ FS /(e 2 ⟨n e ⟩ FS ) > 0. Thus, the potential ϕ(L ∥ /2) becomes positive at the next timestep.

Thus the conducting sheath starts to reflect electrons back into the plasma. If the potential at the wall position were to become negative because of fluctuationsdue to turbulence for instance-no electrons would be reflected. Again, because of the large inertia difference between ions and electrons, the latter would tend to leave the plasma more easily than the ions. This would drive the charge density of guiding center to rise to positive values, causing the potential to become positive again. We have shown that immersed boundary conditions are efficient to recover the kinetic properties of plasma-wall interaction in the parallel direction, along the magnetic field lines. Voice simulations resolve both the quasineutral plasma and the non-neutral sheath. As already noted, sheath physics is sub-scale for gyrokinetics, which focuses on turbulent and transport scales. The question of the extrapolation of the Voice results to gyrokinetics is both timely and important, in particular for Gysela in which penalization has recently been used to study plasma-solid interaction [18] from the standpoint of the turbulent edge. However, techniques akin to those described here may not straightforwardly apply, in particular because the reduction from Poisson equation to gyrokinetic quasi-neutrality postulates a time and scale separation between fast and small-scale sheath equilibration and slower, larger-scale turbulence. The Poisson equation reads

λ D0 L ∥ 2 ∇2 ∥ φ + λ D0 L ⊥ 2 ∇2 ⊥ φ = n e n 0 - Z i n 0 (n Gi + n pol ), (27)
with λ D0 the Debye length at density n 0 and temperature T 0 , L ∥ and L ⊥ the characteristic gradient length of the electric potential in the parallel and transverse directions, respectively. Here, gradient and divergence operators are normalized to these characteristic gradient lengths, and φ = eϕ/T 0 is the dimensionless potential. n e ≈ f Ge d 3 v is the electron density and n Gi = J i .f Gi d 3 v the ion gyrocenter density where J i is the gyro-average operator. The gyro-average operator of a given function g written J i g is defined as the average of g over a cyclotron motion as

J i g = 2π 0 g dφ c 2π ,
where φ c is the cyclotron phase. The ion polarization density n pol reads, in the long wavelength limit L ⊥ ≫ ρ i (with ρ i the ion gyro-radius):

n pol ≈ ∇ ∇ ∇ ⊥ .[(m i n e /eB 2 )∇ ∇ ∇ ⊥ ϕ].
In any case, the scaling L ⊥ ≫ λ D0 is well fulfilled in drift wave interchange turbulence of tokamak plasmas. Turbulent eddies indeed develop at scales larger than the ion gyro-radius, so that the second term in the left hand side of eq.( 27) is neglected in usual orderings. For consistency, it readily appears that the first term can also be neglected to describe turbulence, since L ∥ ≈ qR ≫ ρ i ≫ λ D0 . Notable exceptions to this ordering is in the vicinity of the sheath, where L ∥ ≈ λ D0 , and of the wall when L ⊥ < ρ i . Outside of this thin layer, the left hand side of the gyrokinetic Poisson equation ( 27) is vanishingly small, leading to the quasi-neutrality constraint n e = Z i (n Gi + n pol ). The gyrokinetic quasi-neutrality equation thus only involves at leading order the transverse operator ∇ ∇ ∇ ⊥ , which enters through the expression of n pol and the parallel direction is ignored. Two correlated consequences follow, particularly important when discussing plasma-wall interaction issues in the context of gyrokinetics: (i) the electric potential can be calculated everywhere without any knowledge of its properties along the parallel direction, and (ii) the parallel electric field is a priori unconstrained unless specific boundary conditions are explicitly enforced on the distribution functions and possibly on the transverse electric field in order to mimic the presence of a material boundary somewhere along the open field line. A key question is thus how to use the kinetic knowledge from Voice to infer the correct electron and ion distribution functions as well as the sheath potential in gyrokinetics at the transition between plasma and material boundaries, without resolving the Debye scale of the sheath.

One can first notice important points regarding the distribution functions at the sheath entrance. Considering the case of the wall on the right hand side of the plasma, and up to collisional effects and the re-emission of a few electrons due to the finite value of the Krook coefficients ν wa (cf. Eq. ( 4)), the electron distribution function f e appears to be centered at v = 0 and asymmetric with a sharp cut-off at velocities smaller than the cut-off velocity -v c (cf. Figs. 6-7 and associated discussion). Also, the ion distribution function f i is shifted (maximum at v ̸ = 0) and vanishes in the negative half of velocity space, see Fig. 7.

Combining these observations with the method proposed for a conducting sheath (cf. section V 2 and reference [51]), the following strategy is contemplated for Gysela. Before stepping to the details of the method proposed below, let us first remark that the effect of a tilted magnetic field line with respect to the plasma solid boundary is not taken into account in the present study, since the considered model is one-dimensional in both spatial and velocity direction. In a gyrokinetic code such as Gysela, the magnetic geometry is more realistic and magnetic field lines impact the limiter with non-normal angles given by the value of the safety factor at the separatrix. Nevertheless, even in such more realistic magnetic topology, the first ingredient of plasma-wall interaction that should be included in a gyrokinetic code is a mechanism allowing for the electron and ion particle flux to compensate on average at the plasma boundary. This observation follows one of the main outcomes of the present study: the total current at the wall always vanishes, independently of the choice of the penalization operator in the wall region, the ion-to-electron mass ratio, the elec-tron temperature and the sheath entrance and the magnetic field line length. This vanishing current is therefore a robust property of plasma-wall interaction. When going to tilted magnetic field lines, a magnetic presheath appears that modifies the details of plasma-wall interaction. The combination of the magnetic presheath and of the Debye sheath however still ensures that in average, the current flowing into the solid wall vanishes so that the plasma stays quasineutral on scales larger than the Debye length. Therefore, similarly to the conducting and logical sheath approaches, the future implementation of realistic plasma-wall interaction in the Gysela code will focus firstly in adding a mechanism to limit the electron flux at the divertor/limiter surface, the magnetic presheath effects being left for a possible future study. Lastly, note that here, one considers the plasma-wall interaction such that positive velocities go from the plasma to the wall. The treatment of the other side of the limiter can be obtained by symmetry.

• Define the mask of an axi-symmetric limiter or divertor M w (r, θ), the subscript "w" standing for "wall": M w is equal to 0 in the plasma and 1 in the wall, with an hyperbolic tangent transition from plasma to wall. Here, the width of the transition region is a free parameter, typically a few units of the local gyro-radius.

• Penalize the ion distribution function f Gi using a Krook term similar to the one described in Eq. ( 4) of the form

S Gi = -M w ν f Gi -f T Gi ,
where the target gyrocenters distribution function

f T Gi is f T Gi (r, θ, v ∥ , µ, t) = f Gi (r, θ, v ∥ , µ, t) × H(v ∥ ),
and ν is a constant coefficient that defines the Krook term magnitude. Here H(v ∥ ) is the Heaviside function such that H(v ∥ < 0) = 0 and H(v ∥ ≥ 0) = 1. The fact that the target ion distribution function f T Gi is vanishing for negative velocities ensures that no ions are allowed to come back from the wall into the plasma.

• The electron gyrocenters distribution function f Ge is relaxed towards a target function truncated for velocities v ≤ -v ∥c , in the spirit of the results found in Section III B. To ensure that currents can freely enter and exit the limiter region (possibly forming current loops) the relaxation process is applied on the average of the electron distribution function rather than on the distribution function itself.

More specifically, we propose the following procedure:

1. Define the average of the electron distribution function over a portion Z of the limiter surface ⟨f Ge ⟩ Z . The average is not conducted over the entire limiter due to the parallelization of the code: such a global average would be very costly since it requires a large amount of communications between processors. Using this average, the electron distribution function can be written as f Ge = ⟨f Ge ⟩ Z + δf Ge , where δf Ge = f Ge -⟨f Ge ⟩ Z .

2. Define a critical velocity v ∥c such that on average the electron flux due to electrons for which v ≥ -v ∥c balances the parallel ion flux on the considered limiter portion:

+∞ -v ∥c v ∥ ⟨f Ge ⟩ Z d 3 v = Γ ∥i Z .
3. Using this critical velocity, the average of the electron distribution function is relaxed towards the target function

f T Ge (v ∥ , t) = ⟨f Ge ⟩ Z (1 -H(v ∥ + v ∥c ))
through the evolution equation

∂ t ⟨f Ge ⟩ Z = -M w ν ⟨f Ge ⟩ Z -f T Ge .
4. Solving the latter equation yields an updated value of the averaged distribution function ⟨f Ge ⟩ new Z that is used to reconstruct the electron distribution function for the next timestep:

f Ge (r, θ, v ∥ , µ, t + dt) = ⟨f Ge ⟩ new Z + C × δf Ge ,
Where the C(r, θ, φ) coefficient is adapted such that quasineutrality is ensured locally, i.e.

f Ge d 3 v = f Gi d 3 v.
The steps above should ensure relaxation of the distribution functions towards their expected shapes at the sheath entrance. While the logical and conducting sheath approaches (see Section V 2) impose these constraints on the last grid point of the simulation domain, the penalization technique brings in some flexibility, the transition taking place on several grid points governed by the sharpness of the mask function M. With these conditions, two constraints are satisfied: (i) the relaxed ion distribution function f T Gi is such that no ion comes back from the wall, and (ii) the relaxed electron distribution function f Ge ensures that on average over a portion of the limiter the ion and electron currents compensate. Since no constraint is applied on the local values of the current, the latter are therefore free to evolve. In particular, they are not forced to balance each other locally, so that nonvanishing parallel currents can freely go in and out of the limiter boundary. This degree of freedom is physically desirable as transverse currents should be allowed to exist in the sheath and in a skin depth of the solid material. The protocol described above further ensures that quasineutrality is fulfilled in the whole simulation domain, hence evacuating the problem of explicitly calculating the sheath potential drop at the limiter boundary. The ion flow is not a priori required to become supersonic. A proper definition of a Mach number is non trivial in kinetics but more importantly flow acceleration is an outcome of the modelling and should remain such. Whether the proposed approach (or any existing existing approach for that matter which does not resolve the detailed kinetic physics of the sheath) would provide an accurate estimate of actual current dynamics in the vicinity of the material boundaries is important, complex and remains largely unaddressed.

VI. CONCLUSION

We have discussed the interaction between constituents of a 1D unmagnetised plasma and a material boundary. This unmagnetized plasma is equivalent to a onedimensional magnetized plasma streaming along an open magnetic field that intercepts at normal incidence a wall. The adopted viewpoint is kinetic (1D-1V), for both ions and electrons. The material boundaries are treated using immersed boundary conditions, which allows decoupling of the physical boundary condition of interest-here the sheath-from the mathematical-numerical boundary conditions, often poorly representative of the physical complexity. Two different versions of penalization have been investigated in the Voice framework, either allowing for a conducting wall or preventing currents from flowing in the wall (effectively realizing an insulating plasma-surface interface). These two solutions are reminiscent of the logical and conducting sheath approaches which have been proposed to include plasma-wall interaction physics in kinetic or gyrokinetic codes. A linearized intra-and inter-species collision operator is included, which allows for momentum and energy transfer between ions and electrons. Collisions are found to play a crucial role in preventing velocity-space filamentation and numerical blow up at small scales.

In the results presented here, the system freely evolves towards an equilibrium where a non-neutral positively charged boundary layer a few Debye lengths wide is observed, which screens the quasineutral plasma from material boundaries. Penalization allows the recovery of the main kinetic features of one-dimensional plasma-wall interaction with a magnetic field line that impacts a wall at normal incidence. The sheath layer reflects low energy electrons back into the plasma and ensures ambipolarity at the plasma boundary: the ion particle flux compensates the electron particle flux there. These kinetic features have been retrieved for two different versions of the penalization operator, either allowing currents to exist in the wall region-as a result of the absorption of electrons and ions in the wall on different characteristics lengths-, or without. Note that these two cases-in addition to allowing one to test the robustness of the immersed boundary condition technique in two different settingsare also relevant from a physical standpoint, as they address different wall properties that can be encountered experimentally. The impact on plasma dynamics of such currents in the wall remains an open issue. The potential drop in the sheath has been computed for both conducting and insulating wall solutions, and its dependence on key simulation parameters, namely the electron-to-ion mass ratio, the electron temperature and the magnetic field line length have been investigated. The obtained values for the potential drop have been compared with an analytical prediction that derives from the simplest sheath model one could think of, in which the sheath reflects slow electrons in the plasma and lets fast ones be absorbed in the wall. The cutoff between the fast and slow populations is determined by the equality of ion and electron fluxes at the entrance of the sheath. Given this analytical model, the departure of the sheath potential drop computed from the simulation results in the conducting and insulating wall cases has been quantified and its origin has been detailed. The observed discrepancies mainly come from the re-emission of fast electrons in the plasma. This re-emission is not necessarily encountered in actual plasmas as it originates from (i) collisions happening in the Debye sheath but also comes from (ii) the finite strength of the penalization operator. This re-emission has been shown to be more important in the insulating wall case, as more electrons are present in the wall. Additionally, the conducting wall solution implies less arbitrary choices for fixing the numerical parameters of the simulation as compared to the insulating wall case. The conducting wall solution thus appears preferable to describe one-dimensional plasma-wall interaction using immersed boundary conditions. The insulating wall remains nevertheless of interest since the question of the distribution of wall currents in higher dimensions remains open, in particular when taking the cross-field transport into account. When taking the reemission phenomenon into account, the potential drop from the simulations can be reconciled with the analytical prediction. This gives confidence in the fact that the immersed boundary condition technique can be used to describe one-dimensional plasma-wall interaction. Furthermore, we would like to stress that techniques other than penalisation also have important shortcomings that are less present in our present approach. Such shortcomings mostly revolve around assumptions made on currents flowing to the wall or flowing across flux surfaces in higher dimensionality. This point is briefly touched upon whilst discussing the logical and conducting sheath methods. We believe that penalisation does allow in that respect for a more controlled and flexible way to deal with such current loops, which could be at least as important to the physics than the possible re-emission of unwanted electrons. We would also like to stress that to a large degree the amount of unwanted electrons that are re-emitted due to penalisation only is dependent on the numerical price one is willing to pay through strength of the Krook coefficient and resolution (steepness of the masks). We however agree with the referee that these questions (reemission of electrons and current loops through sheath and wall) and the trade-off that they imply between numerical cost and impact on the physics should be carefully assessed by comprehensive comparisons of various boundary techniques, especially in higher dimensionality. These questions are however left for future work.

Whilst computing the sheath heat transmission factor, we find that the ion heat transmission factor is 60% larger than usually assumed in fluid prescriptions. This result can only be captured through kinetic modeling as it is found to be primarily the consequence of the ion distribution function departing from a Maxwellian. The deposition of heat by the plasma on the wall (divertor) is found to be carried 8% more by ions than by electrons in the kinetic code, on the contrary to the common prediction assuming shifted (ions) or truncated (electrons) Maxwellians that states that electrons carry more energy to the wall than ions.

We have shown here that immersed boundary conditions, in the form of penalization BGK (Krook) terms in the Vlasov equation, are efficient to model and recover the kinetic physics of the interaction between the plasma and a solid material. The simplicity of the model allows one to understand the full extent of the kinetic physics at play and in particular both the potential drop in the sheath and the shapes of the distribution functions at the sheath entrance. In contrast to the logical and conducting sheath techniques, here the sheath boundary layer is actually resolved. However, this level of detail is impractical in higher dimensionality, e.g. in flux-driven global gyrokinetic turbulence codes where the required resolution would be orders of magnitudes greater than currently affordable. An important question is therefore whether-and whereby-the present method based on penalization of the plasma dynamics can be extended to gyrokinetics to include kinetic features of plasmasurface interaction without resolving the Debye sheath, or the magnetic presheath. To this end, we discussed a prospective solution that should allow the incorporation of sheath physics in gyrokinetics without resolving the sheath. The proposed solution bears similarity with the conducting sheath method, while using the flexibility of immersed boundary conditions. This prospective method is under active investigation in Gysela and may provide an alternative-and complementary-approach to the logical and conducting sheath methods proposed in the literature, whilst adding the versatility of penalization.

This question is timely given the central role that outer boundary conditions on the open field lines may have for heat or particle deposition loads or access to improved confinement. In the perspective of implementing plasmawall interaction physics in global gyrokinetics, we took special care to assess robustness of key observables (elec-tric potential drop and the distribution functions) with respect to relevant physical parameters such as ion to electron mass ratio, magnetic field length, abruptness of the plasma-wall transition or the possibility of net currents flowing across the sheath. The potential drop in the Debye sheath region is robustly proportional to the electron temperature. Interestingly for gyrokinetics, we find that many of the basic features of plasma-wall interaction are robust even at (non-realistic) low electron-to-ion mass ratios. Details of the sheath physics are also weakly dependent on the parallel plasma length (or connection length of the Scrape-Off layer) in the investigated range.

Important open questions remain, left for future work: what is the impact of non-normal (grazing) angles between open field lines and wall, i.e. what is the impact of the magnetic presheath on the the onset and the depth of the potential drop, as well as on the shapes of the distribution functions? With what implications with respect to energy deposition on the first wall? Could neutrals modify sheath equilibration to the extent that such corrections should be incorporated in models which do not resolve the sheath, such as gyrokinetics? All of the quantities described in the following are normalized unless the contrary is explicitly specified. The intra species collision operator C aa (f a ) is inspired from [START_REF] Dif-Pradalier | [END_REF]. It conserves density, momentum and energy. It is written as

C aa (f a ) = ν ⋆ D0 ∂ va D v f Ma ∂ va f a f Ma . (A1)
This operator automatically ensures particle conservation since dv C aa (f a ) = 0. To furthermore fulfill momentum and kinetic energy conservation, we adjust the fluid moments of the Maxwellian function f Ma . We set the density of f Ma as equal to the density of f a . Then, its two other moments u Ma and T Ma are adapted at each time step and spatial position. Namely the two equations dv a v a C aa (f a ) = 0 and dv a v 2 a C aa = 0 are solved for u Ma and T Ma . This process gives

u Ma = -⟨v a D v ⟩⟨∂ va (v a D v )⟩ + ⟨v 2 a D v ⟩⟨∂ va D v ⟩ ⟨v a D v ⟩⟨∂ va D v ⟩ -⟨D v ⟩⟨∂ va (vD v )⟩ , (A2) T Ma = ⟨v a D v ⟩ 2 -⟨D v ⟩⟨v 2 a D v ⟩ ⟨v a D v ⟩⟨∂ va D v ⟩ -⟨D v ⟩⟨∂ va (v a D v )⟩ , (A3)
where the bracket notation is defined as an integral against f a , for instance ⟨D v ⟩ = dv a f a D v . Note that these two moments depend on both space and time: u Ma (x, t) and T Ma (x, t). In the case where f a itself is a Maxwellian, then u M a and T M a are its actual mean velocity and temperature. The diffusion coefficient in velocity space D v of Eq. (A1) depends on the velocity and space variables. It is written as

D v = D 0 Φ -G y a , (A4)
with

D 0 = 3 √ 2π 4 ν ⋆ aa T a and y a = |v a | / 2T a .
Note that the y a variable is also a time and space dependent quantity since in this expression T a is the temperature of species a at the considered position x. The Φ and G functions are defined as

Φ(y a ) = 2 √ π ya 0 e -z 2 dz, G(y a ) = Φ -y a Φ ′ 2y 2 a .
The self-collision frequency ν ⋆ aa is derived from the classical expression of the two-species collisions frequency ν ⋆ ab defined in Eq. (A7). These frequencies are normalized by the reference plasma frequency ω pe0 . It is expressed as

ν ⋆ aa = A a Z 4 a n a T 3/2 a . (A5)
This collision frequency is also a function of space and time since it involves the density and temperature of species a at point x and time t. The ν ⋆ 0 parameter is used as an input of the code to control the plasma collisionality, and represents the number of collisions experienced by a particle when traversing the entire simulation box, see Eq. ( 13) for its precise definition. The two species collision operator is inspired from [42]. It takes the form

C ab (f a ) = ν ⋆ D0 2Q M ab n a T a m a 2T a (v a -V a ) 2 - 1 2 F Ma (A6) + ν ⋆ D0 R M ab n a T a (v a -V a )F Ma ,
where a and b are the two species colliding with each other. In the above expression the F Ma function is a

Maxwellian that has the same first three moments as the distribution function f a . Note that this Maxwellian F Ma is different from the one f Ma appearing in the definition of the self-collision operator, Eq. (A1). The friction force term R M ab and the thermal energy transfer term Q M ab are expressed as

R M ab = -n a ν ⋆ ab u a - A b A a u b , Q M ab = -3n a A b A b + A a ν ⋆ ab (T a -T b ) -u a R M ab .
Or equivalently, using dimensional quantities for all the variables involved

R M ab,dim. = -n a m a ν ab (u a -u b ), Q M ab,dim. = -3n a m a m a + m b ν ab (T a -T b ) -u a R M ab .
In these latter expression one can clearly check that the friction term describing the momentum exchange is a restoring force proportional to u a -u b , which consequently vanishes when u a = u b . Similarly, the energy equipartition term exhibits a thermal equipartition part acting as a restoring force proportional to T a -T a , hence vanishing when T a = T b together with a term describing the energy exchange governed by the momentum exchange. The inter species collision frequency is written in normalized units as:

ν ⋆ ab = 2A a n b Z 2 a Z 2 b T 3/2 a 1 + A b A a 1 + A b T b A a T a -3/2 (A7)
One thus recovers the collision frequency ν ab , in dimension units since ν ab = ν * 0 ν * ab (λ D0 /L ∥ )ω pe0 and therefore:

ν ab = 4 √ 2π 3 n b m a e a e b 4πε 0 2 × log Λ 1 (v 2 Ta + v 2 T b ) 3/2 m a + m b m a m b .
Where v Ta is the thermal velocity of species a at the considered spatial location: v T a (x, t) = T a (x, t)/m a , and log Λ is the Coulomb logarithm that also depends on the plasma properties. expected, the higher the collisionality, the closer to a Maxwellian the electron distribution function becomes, see the blue dotted line for ν ⋆ 0 = 1 and the blue dashdot line for ν ⋆ 0 = 0.1. As mentioned in Section III B, we see that the electrons of small velocities have a distribution that tends to the kernel Maxwellian whatever the value of the collisionality. This is not the case for the higher velocity part of the distribution function that becomes Maxwellian as the collisionality increases. This dichotomy can be understood by looking at the collisional diffusion coefficient of Eq. (A4). This coefficient, whose magnitude determines the strength of the self scattering processes, decreases with the norm of the velocity. Therefore, fast particles are less sensitive to collision compared to slow ones. This results in a distribution function that looks like a sum of two Maxwellians. Note that a similar double Maxwellian structure has already been observed in [43], that shares a common origin with the present case. The kinetic source term presented in Section II B 1 injects a Maxwellian contribution to the plasma at each timestep. This source models cross-field transport. Since the sheath acts as a filter that removes fast electrons from the plasma, this boundary layer has a cooling effect on electrons. Thus, high velocity particles injected by the source appear as supra-thermal electrons when looking at FIG. 14. These supra-thermal electrons travel towards the entrance of the sheath and are absorbed in the wall. solve the coupling between Boltzmann and Poisson equation, a midpoint method is used. The Poisson solver uses a finite elements method. Second order Strang splitting is used to solve the Boltzmann equation. The advection operators are handled with a backward semi-Lagrangian scheme with spline interpolation. The particle sink, heat sink and inter species collisions operators (4,8) are solved using a second order Runge-Kutta scheme. The source term ( 10) is solved analytically, and a Crank-Nicholson method is used for the intra species collisions. The x and v directions are discretized using a nonuniform mesh. More precisely, the mesh is split between several uniform zones. Each zone contains evenly distributed points, and the density of points varies between each zone. The companion paper [34] details the algorithm that generates this non-uniform mesh by refining the discretization in the regions of phase space where a physical quantity of interest-the distribution function in our case-presents large variations. It is also shown that GPU parallelization is more efficient for non uniform grids than for uniform grids. The efficiency in this case is defined as the ratio between simulation times for the serial parallelized versions of the code respectively. In the present case a total number of N x = 1500 points is used to discretize the spatial direction. The simulation box has a total length of L x = 700λ D0 in the reference case. This length includes the wall region. Most of the points are placed in the vicinity of the sheath, where steep gradients develop. Conversely the quasineutral plasma and the wall region are discretized with a lower density of points since the distribution function is more regular there. Overall the point density can be viewed on Fig. (15). The figure only represents the right half of the simulation box, and the discretization in its left part is symmetric. A total of 5 uniform zones is used.

N v = 700 points are used in the velocity direction. The velocity grid spans over the interval [-v max , v max ]. The maximum velocity v max is equal to six thermal velocities, i.e. v max = 6v T0 . The points are distributed between three uniform discretization zones. The first one is the [-3v T0 , 3v T0 ] region of the velocity interval, which is discretized using 600 points. The remaining 100 points are used to discretize the rest of the velocity interval.

A list of the most important parameters used in the ref-erence simulation can be found in Table I. L ∥ represents the size of the plasma in the simulation box. The d w parameter defines the steepness of the wall mask functionsee Section II A 2 and Eq. ( 6). The target density in the wall n w is defined in Eq. ( 5). The collisionality ν ⋆ 0 is defined in Eq. [START_REF] Persson | [END_REF]. Note that typical values of density and temperature in the Scrape-Off Layer of about 2.5 10 19 m -3 and 100 eV would lead to a collisionality ν ⋆ 0 of order unity. The large value considered in the simulations (cf. Table I) is a compromise between realistic physics and numerical constraints. Indeed, the smaller ν ⋆ 0 , the smaller the grid size ∆v in the velocity space (i.e. the larger the number of velocity grid points) to avoid negative values of the distribution function. These are likely to occur when the distribution functions develop velocity gradient lengths of the order of or less than the velocity step size. The collision operator helps maintaining this length at a finite value, all the bigger since ν ⋆ 0 is large. The relation between ∆v and the collisionality ν ⋆ 0 can be obtained as follows. The distribution functions evolve at the plasma frequency ω pa = v T a /λ D . Balancing this dynamics with the diffusive part of the collision operator -of characteristic transport coefficient

D v ≈ v 2 T a ν aa -leads to ∆v/v T a ∼ (ℓ mfp /λ D ) 1/2 ∼ (ν ⋆ 0 λ D0 /L ∥ ) 1/2
, where ℓ mfp stands for the mean free path. Taking L ∥ /λ D0 = 7.10 5 (cf. Section II B 2) and ν ⋆ 0 = 50 leads to ∆v/v T a ≈ 10 -2 , which is about the considered velocity grid step (∆v/v T0 = 6/600 in the interval [-3v T0 , 3v T0 ]). Convergence tests have been performed, showing that the main results presented in this paper do not depend on the precise value of ν ⋆ 0 .

Appendix D: Fluid equations for code verification For the sake of numerical efficiency the velocity normalization used in the kinetic framework is species dependent, while both time and space normalization are species independent. This particularity introduces an explicit dependence in √ A a in the reference kinetic equation Eq. (3a). This control parameter is the ratio of the reference kinetic convective frequency and the reference plasma frequency. It weighs all phase space convection terms, including collisions. It is furthermore assumed that the characteristic time governing the source is species independent. With this normalization, we now define the various moments in the following way, retaining the same notations for the normalized quantities, con-

T a = +∞ -∞
dv a v a -u a 2 f a /n a (D3)

Q a = +∞ -∞ dv a v 3 a f a (D4)
where Q a stands for the energy flux of species a. Note that the heat flux is consequently defined as

q a = +∞ -∞ dv a (v a -u a ) 3 f a . ( D5
)
We define then define the source terms as:

S na = +∞ -∞
= ∂ t n i + √ A i ∂ x (n i u i ) -S ni max x∈[0,Lx] |∂ t n i | , √ A i ∂ x (n i u i ) , |S ni | ,
and it is computed in a similar fashion for the other conservation equations for both species. We see on FIG. (16) that this error term never exceeds 10 -2 , hence the numerical error is under a percent.

Appendix E: Prediction of the sheath heat transmission factor

We recall in the following a derivation of the sheath heat transmission factor for both electrons and ions species. This derivation is well known in the literature and can be found for instance in [35, p. 92]. As stated in Section III C the ion and electron sheath heat transmission coefficients are defined respectively as

γ e = Q e T e Γ γ i = Q i T i Γ ,
where the energy flux Q a for species a is defined in normalized units in Eq. (D4) or equivalently in dimensional units in Eq. ( 16). Note that electron and ion particle fluxes are assumed to be equal at the plasma boundary, i.e. Γ i = Γ e = Γ. In the following of this section all quantities are taken at the entrance of the sheath. As explained in Appendix F, we use the conservation of energy for electrons to define the critical velocity v c . This quantity represents the minimal velocity an electron must have at sheath entrance to overcome the repulsing electric field of the sheath and hit the wall:

v c = -2e∆ϕ sh m e . ( E1
)
where ∆ϕ sh is the negative electric potential drop in the sheath region. We consider the wall as perfectly absorbent. Therefore no electron can escape the wall with a finite velocity, so that their velocity is at most -v c at the sheath entrance: there, electron velocities satisfy v ≥ -v c . It follows that the electron distribution function at sheath entrance is necessarily truncated for velocities v ≤ -v c . Further assuming a Maxwellian shape leads to And therefore the electron heat transmission factor is in this case

γ e = 1 - e∆ϕ sh T e ≈ 4.
The last equality is justified by the prediction of the potential drop in the sheath given in Eq. (F4) that yields for a Deuterium plasma -e∆ϕ sh /T e ≈ 3. One should keep in mind that the expression obtained for the electron energy flux in Eq. (E2) is only valid in a simplified 1D geometry. When going towards more realistic 3D models one has to also take into account the component of the velocity vector perpendicular to the magnetic field along with its parallel component. This leads to a factor 2 in front of the T e Γ e term of Eq. (E2). Consequently, the prediction for the electron heat transmission factor is modified so that γ 3D e ≈ 5.

Regarding the ion distribution function at sheath entrance, we consider a Maxwellian of temperature T i with a mean velocity u i = M c s , where c s is the plasma sound wave velocity. Without any loss of generality the sound speed can be defined as m i c s = χT , with T = T e + T i and χ is a coefficient that depends on the chosen closure of the fluid hierarchy, see the companion paper [32]. For instance, with an isothermal closure T e = T i and we find χ = 1. The energy flux associated to such distribution function reads

Q i = ΓT i 3 2 + χM 2 2
T e T i + 1 . ( E3
)
Therefore with we find the ion heat transmission factor at the sheath entrance

γ i = 3 2 + χM 2 2
T e T i + 1 .
Assuming an isothermal closure of the fluid hierarchy we further get χ = 1, and T e = T i , thus γ isoth. i = 5/2. In this case the total sheath heat transmission factor is γ = γ isoth. i + γ e ≈ 6.5.

Lastly, the same remark as the one formulated for the electrons applies here. The value found above for the ion sheath heat transmission factor depends on the dimensionality of the model. When going towards a 3D model, the ion energy flux expression is modified as

Q 3D i = ΓT i 5 2 + χM 2 2
T e T i + 1 .

Thus with the isothermal closure one would get γ 3D,isoth. i = 3.5 and then γ 3D ≈ 7.5.

Appendix F: Prediction of the potential drop in the sheath

Classical potential drop prediction

The sheath potential drop prediction detailed below is taken from [START_REF] Stangeby | The Plasma Boundary of Magnetic Fusion Devices[END_REF]. Let us consider a perfectly absorbing wall at position x w , bounding a semi infinite plasma located at x ≤ x w . We consider the wall to be grounded and to not emit any secondary particles. In front of the wall a sheath develops to ensure ambipolarity. As usual x sh refers to the sheath entrance position. The flux of ions at the sheath entrance can be written as

Γ sh i = n sh i M c s . (F1)
The Mach number M and sound speed velocity c s are taken at the sheath entrance as well. We ignore the effects of collisions in the sheath. Under these assumptions the Hamiltonian function for electrons is simply H e = 1 2 m e v 2 -eϕ. Since this quantity does not explicitly depends on time it is conserved along the electrons trajectories. On a particular trajectory connecting the sheath entrance to the wall we can therefore write that 1 2 m e v(x sh ) 2 -eϕ(x sh ) = 1 2 m e v(x w ) 2 -eϕ(x w ), or equivalently v(x w ) 2 = v(x sh ) 2 + 2e∆ϕ sh m e .

Here we see that a necessary condition for the electrons to reach the wall is v(x sh ) ≥ v c = -2e∆ϕ sh /m e . The electrons that are faster than the velocity v c at the sheath entrance are therefore absorbed in the wall. The rest of the electrons are reflected back into the plasma. Since the wall does no emit any secondary electron, the above calculation shows that there are no electrons with velocities v(x sh ) < -v c at the sheath entrance. We therefore assume that the electron distribution function at the sheath entrance is a Maxwellian truncated for velocities v ≤ -v c . This truncated Maxwellian is characterized by a density n e and a temperature T sh e . It is written as Note that the quantity n e is not the electron density at the sheath entrance n sh e . The density at sheath entrance n sh e can be found by integrating against velocity Eq. (F2). By doing so, we find The last equality holds because the v ≥ v c part of the Maxwellian at the sheath entrance is exactly the fast electrons that are found at the wall position. If the electron density at the wall is negligible compared to the density at the sheath entrance, we can make the assumption n e ≈ n sh e . It leads to the well-known prediction of the sheath potential drop, see below Eq. (F4). In our simulations this assumption is not always valid. We observed nevertheless that taking into account the correct density at the sheath entrance n sh e = n e -n e (x w ) did not change noticeably the computed value of the potential drop from the simulations. We therefore consider in the following that n e ≈ n sh e . The expression of the electron flux at the sheath entrance is therefore As explained in Section IV A we observed some discrepancies between the potential drop obtained from the simulations and the prediction detailed above. These discrepancies are mainly due to electrons coming back from the wall region. In the simple model of Section F 1, the electron distribution function at the sheath entrance vanishes for v ≤ v c , that is to say no fast electrons come back from the wall region. On FIG. (6) wee see that it is not the case in our simulations. To express the equality of the ion and electron flux at the sheath entrance, it is therefore more accurate to write In this expression f sh e is the distribution function at the sheath entrance obtained from the simulation. Note that Γ - e is a positive quantity, therefore the corrected electron flux at the sheath entrance is lower than the value obtained with a perfectly truncated Maxwellian. After some straightforward calculation we find that the prediction for the potential drop in the sheath is now ∆ϕ ref.

1 FIG. 1 .

 11 FIG. 1. (a) Poloidal projection of a tokamak magnetic geometry. (b) Simplified geometry with no curvature. The field line represented as a dashed blue line hits the divertor walls at both ends. (c) Masks used to define the wall and plasma area for a typical simulation whose parameters can be found in TableI.

  FIG.3. Total current normalized to the ion saturation current in the sheath vicinity, for the conducting wall in plain red, and the no-current wall in dashed blue. The spatial axis has been renormalized so that the sheath entrance and wall position coincide in the conducting and no-current cases. The wall mask defined in Eq. (6) is represented in dashed grey.

FIG. 4 .

 4 FIG.4. Charge density in the sheath vicinity, for the conducting wall in plain red, and the no-current wall in dashed blue. A renormalization of the spatial axis has been performed, as in FIG.(4). The potential is represented in black, in the conducting wall case. ϕ sh is the electric potential taken at the sheath entrance in the insulating wall case.

FIG. 5 .FIG. 6 .FIG. 7 .

 567 FIG. 5. Distribution functions in phase space. Top left: ion distribution function, bottom left: electron distribution function.The right graphs represent a zoom of both distribution functions in the vicinity of the sheath. For clarity a different color scale is used in the zoomed graphs as compared to the two left graphs. The gray boxes on the two left graphs remind the reader that the right graphs are a zoom of some portion of the initial left graphs. The boxes do not however cover the exact zoomed area, but are drawn bigger to make them more visible.

e

  is the predicted forward flux of electrons with velocities v ≥ v c . The actual forward flux of electrons with velocities v ≥ v c from the simulations is rather expressed as

3 . 4 .

 34 More electrons than ions have crossed the boundary, i.e. N e > N i . In this case: (a) Order all N e electrons from slowest to fastest. (b) Remove all N i ions from the simulation, and the fastest N i electrons. The velocity of the slowest electron that is absorbed in this step defines a cutoff velocity v c . (c) The remaining N e -N i electrons that have not been absorbed have velocities lower than the cutoff velocity. These electrons are reflected back into the plasma, by flipping the sign of their velocities for the next time step. More ions than electrons have crossed the boundary, i.e. N i > N e . In this case: (a) All these N i ions and N e electrons are absorbed by the wall, and are thus removed from the simulation. (b) To fulfill charge conservation in this case, a positive charge N i -N e is attached to the plasma boundary position at z = L ∥ /2 for the next time step only.

FIG. 13 .

 13 FIG.[START_REF] Persson | [END_REF]. Velocity grid in the z direction at the wall position L ∥ /2. Because the cutoff velocity falls inside the vj velocity cell, when reflecting the slow electrons to the plasma the distribution function in this cell has to be rescaled.

Appendix B : FIG. 14 .

 :14 FIG. 14. Electron distribution function at the sheath entrance for different values of the collisionality in solid lines, and kernel collisional Maxwellian in dashed lines.

FIG. 15 .

 15 FIG. 15. Points density in the right part of the simulation box. The blue curve is the wall mask function. The red curve corresponds to the density of points in each uniform discretization zone.

  FIG.16. Top, from left to right: conservation of mass for electrons, ions and conservation of momentum for electron. Bottom, from left to right: conservation of momentum for ions, conservation of energy for electrons and ions.

  dv a S(f a ) S ua = +∞ -∞ dv a v a -u a S(f a ) S ha = +∞ -∞ dv a v a -u a 2 S(f a )With these definition one obtains:∂ t n a + A a ∂ x (n a u a ) -S na = 0 (D6) ∂ t (n a u a ) + A a ∂ x (n a T a + n a u 2 a ) -Z a A a n a E -S ua -R V ab = 0 (D7) ∂ t (n a T a ) + A a u a ∂ x (n a T a ) + A a (∂ x Q a + 3n a T a ∂ x u a ) -S ha -2Q V ab = 0 (D8)where Eq. (D6) expresses the conservation of the plasma density, Eq. (D7) the conservation of fluid velocity and Eq. (D8) the conservation of thermal energy. These equations are obtained by integrating the Boltzmann Eq. (1) against powers of v. The intra species collision operator C aa defined in Eq. (A1) does not intervene in the fluid conservation equations since by construction it conserves density, momentum and energy, see Appendix A. Integrating the inter species collision operator C ab defined in Eq. (A6) against powers of v a simply yields the friction force term R M ab in the momentum conservation equation (D7), and the thermal energy transfer term Q M ab in the energy conservation equation (D8). All the terms involved in the fluid conservation equations are shown on FIG.(16). The accuracy of the simulation can be checked by looking at the value of the relative error term denoted by |lhs -rhs| rel. on the same figure, right axis. This term is computed by substracting the right-hand side part from the left-hand side part of each conservation equation. Then, this term is divided by the maximum value of all the terms of the equation. To give a concrete example, the |lhs -rhs| rel. term is computed for the equation expressing the conservation of density for ions as |lhs -rhs| rel.

  f e (x sh , v) = n e m e 2πT e e -m e v 2 2T e H(v + v c ), where H(v +v c ) is the Heaviside function that equals zero for v ≤ -v c and one otherwise. n e is the electron density at sheath entrance. After some straightforward calculations and using the definition of v c given by Eq. (E1) we find the energy flux of this distribution function Q e = T e Γ e -e∆ϕ sh n e e e∆ϕ sh T e T e 2πm e (E2) Using the expression of the electron flux at sheath entrance given in Eq. (F3) we rewrite this energy flux as Q e = T e Γ e -e∆ϕ sh Γ e .

  f e (x sh , v) = n e m e 2πT sh e e -m e v 2 2T sh e H(v + v c ), (F2) = f e H(v + v c ),where H(v +v c ) is the Heaviside function that has a value of one when v ≥ -v c , and zero otherwise. For the sake of compactness of the following computation we introduced the f e notation that stands for the non-truncated Maxwellianf e = n e

=

  n e -n e (x w ).

2 .

 2 ion flux Eq. (F1) and electron flux Eq. (F3) we find the classical expression of the potential drop in the sheath region ∆ϕ sh = On the undesired electron re-emission from the wall

  corresponds to the expression of the electron flux at the sheath entrance in the framework of the classical model developed in the last section. The Γ - e accounts for the flux of fast electrons coming back from the wall at the sheath entrance, i.e.

TABLE I .

 I Reference simulation parameters.

	Lx	2L ∥	nsrc	Tsrc	dw	nw	Tw	ν ⋆ 0
	700	406	1	1	0.2	10 -9	1/3	50

Conducting and insulating immersed boundary conditions.

Acknowledgments

The authors would like to acknowledge E. Serre for his early involvement in the development of immersed boundary conditions to address plasma-wall interaction in numerical simulations, E. Caschera for the first application to the Voice code of this technique, S. Brunner and A. Geraldini for sharing their insight on the adaptation of the sheath model to a gyrokinetic framework, and D. Told, D. Tskhakaya and all members of the TSVV-04 project for many fruitful discussions. We acknowledge as well Claudia Negulescu and Etienne Lehman for bringing valuable remarks from the mathematicians community into our work. 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 EUROfusion) 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. We acknowledge the Centre de Calcul Intensif d'Aix-Marseille University for granting us access to its high performance computing

Appendix C: Numerical scheme and simulation parameters

The semi-Lagrangian version of the Voice code [39] is used to solve numerically the Boltzmann-Poisson system (3a), in a very similar way to the one used in Gysela. To