Kinetic plasma-wall interaction using immersed boundary conditions

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 $\lambda_\mathrm{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 $\Delta \phi$ spread across a characteristic gradient length $\ell$ relates to the charge density $\rho_\mathrm{c}$ that gives rise to this potential drop as

$$\begin{align*}\dfrac{e\Delta \phi}{T_0} \approx -\dfrac{\ell^2}{\lambda_\textrm{D}^2} \dfrac{\rho_\mathrm{c}}{n_0 e},\end{align*}$$

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\Delta \phi^\mathrm{presheath} / T_0 \approx -0.7$ and $e\Delta \phi^\mathrm{sheath} / T_0 \approx -3$ [34]. The typical gradient length in the sheath is given by the Debye length, i.e $\ell^\mathrm{sheath} \approx \lambda_\textrm{D}$ whereas in the presheath the gradient length is given by the length of the magnetic field line, i.e. $\ell^\mathrm{presheath} \approx L_\parallel$. For a typical magnetic field line of the Scrape-Off layer one has $L_\parallel \gg \lambda_\textrm{D}$, thus the charge density in the sheath is very large as compared to the charge density in the presheath as

$$\begin{align*}\rho_\mathrm{c}^\mathrm{presheath}/\rho_\mathrm{c}^\mathrm{sheath} \approx \dfrac{\lambda_\textrm{D}^2}{L_\parallel^2} \ll 1.\end{align*}$$

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

$$\begin{align*}E^\mathrm{presheath}/E^\mathrm{sheath} \approx \dfrac{\lambda_\textrm{D}}{L_\parallel} \ll 1.\end{align*}$$

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 so-called 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 [36].

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_\parallel$. 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 5.1 and 5.2 in section 5. 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.

2.1.1. The wall penalization by a particle sink.

The immersed boundary condition is introduced in the Boltzmann equation (3a ) using the particle sink term $\mathcal{S}_{\mathrm{sink,n}}$. It is a Bhatnagar–Gross–Krook operator [30] that acts as a particle, momentum and energy sink. It is written as

$$\begin{equation} \mathcal{S}_{\mathrm{sink,n}}\left(f_a\right) = - \mathcal{M}_\mathrm{w} \nu_{\mathrm{w}a} \left( {\,f_a - g_\mathrm{w}} \right). \end{equation} \tag{ 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

$$\begin{equation} g_{\mathrm{w}}\left(v_a\right) = \dfrac{n_\mathrm{w}}{\sqrt{2\pi T_\mathrm{w}} } \, \operatorname{exp} \left( -\dfrac{v_a^2}{2T_\mathrm{w}} \right). \end{equation} \tag{ 5 }$$

It is chosen to be species independent. The density $n_\mathrm{w}$ and temperature $T_\mathrm{w}$ are constant in space and time. The wall target density $n_\mathrm{w}$ is orders of magnitude smaller than the plasma density, so the operator acts as a particle sink. Typically we take $n_\mathrm{w} \approx 10^{-9}$ in the simulations. The choice for the particular value taken for $n_\mathrm{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_\mathrm{w} \ne 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 [33]. 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_\mathrm{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_\mathrm{w}$ is smaller than the source temperature that injects particles in the plasma, see section 2.2.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 $\mathcal{M}_\mathrm{w}(x)$. This mask function therefore defines the wall region in the simulation domain. It is expressed as

$$\begin{align} \mathcal{M}_\mathrm{w}\left(x\right) = 1 -\dfrac{1}{2} \left\{\operatorname{tanh} \left( \dfrac{x-x_{\mathrm{w}}^\mathrm{\ell}}{d_\mathrm{w}} \right) - \operatorname{tanh} \left( \dfrac{x-x_{\mathrm{w}}^{r}}{d_\mathrm{w}} \right) \right\}. \nonumber \end{align} \tag{ 6 }$$

This function is represented on figure 1, bottom left graph. The $x_{\mathrm{w}}^\mathrm{\ell}$ and $x_{\mathrm{w}}^\mathrm{r}$ 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 $\mathcal{M}_\mathrm{w}(x) = 1$ while the plasma region is such that $\mathcal{M}_\mathrm{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_\mathrm{w}$ . The characteristic width of this intermediate region is $d_\mathrm{w}$, typically a fraction of the Debye length, see table 1. The plasma-wall interfaces are therefore blurred, which must be accounted for when analyzing the results. Their characteristic location are $x_{\mathrm{w}}^\mathrm{\ell}$ (left) and $x_{\mathrm{w}}^\mathrm{r}$ (right), these being symmetric with respect to the middle of the simulation box. The third parameter in the chosen restoring force equation (4) is the species dependent amplitude $\nu_{\mathrm{w}a}$. Here we fix the value $\nu_{\mathrm{w}i}$ 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 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 equations (3a ) and (3b ). Actual boundary conditions are imposed at $x=0$ and $x = L_x$, far from the sheath region, see figure 2.

Zoom In Zoom Out Reset image size

Table 1. Reference simulation parameters.

Lx	$2L_\parallel$	$n_\mathrm{src}$	$T_\mathrm{src}$	$d_\mathrm{w}$	$n_\mathrm{w}$	$T_\mathrm{w}$	$\nu_0^{\star}$
700	406	1	1	0.2	10−9	1/3	50

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_\mathrm{w} \to 0$ in the definition of the wall mask equation (6), (ii) $n_\mathrm{w} \to 0$ in the definition of the target Maxwellian equation (5) and (iii) $\nu_\mathrm{w} \to +\infty$ in the expression of the BGK operator equation (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 $\mathcal{S}_{\mathrm{sink,n}} = 0$, hence the plasma is not affected by the restoring force. To be completely rigorous we need to add the constraint $\nu_\mathrm{w}\mathcal{M}_\mathrm{w} \to 0$ in the plasma region. In the wall region the restoring force is infinitely strong, hence $f_a(x\in \, \text{wall}) = g_\mathrm{w} = 0$. Thus, any particle (electron or ion) that comes into contact with the wall is absorbed.

2.1.2. Conducting and insulating immersed boundary conditions.

Setting $\nu_{\mathrm{w}e} = \nu_{\mathrm{w}i}$, 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_{Ta}\partial_{x}n_a(x) = -\nu_a n_a(x)$, one finds an exponential fall-off where the e-folding length $v_{Ta} / \nu_a$ is the wall stopping distance for particles of species a. We have chosen here the thermal velocity $v_{Ta}$ as characteristic velocity. For the chosen absorption amplitude $\nu_{\mathrm{w}e} = \nu_{\mathrm{w}i}$, 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 2.1.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 equation (4) the normalized charge source is $\mathcal{S}_\rho = -\mathcal{M}_\mathrm{w} \left\{\nu_{\mathrm{we}}\left( n_\mathrm{e}(x,t) - n_\mathrm{w} \right) + \nu_\mathrm{wi}\left( n_\mathrm{i}(x,t) - n_\mathrm{w} \right) \right\}$. Then setting:

$$\begin{equation} \nu_\mathrm{we}\left(x,t\right) = \nu_\mathrm{wi}\, \dfrac{n_\mathrm{i}\left(x,t\right) - n_\mathrm{w}}{n_\mathrm{e}\left(x,t\right) - n_\mathrm{w}}. \end{equation} \tag{ 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 $\nu_\mathrm{wi} = \text{constant}$ and $\nu_\mathrm{wi}$ determined by equation (7), one observes a charge build-up in the wall region, while the electrical current remains small, section 2.1.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, $\nu_\mathrm{wi} = \text{constant}$ and $\nu_\mathrm{we} = \text{constant}$, but such that $\nu_\mathrm{we} = \nu_\mathrm{wi} \sqrt{m_\textrm{i}}/\sqrt{m_\textrm{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_\mathrm{e}$ or $n_\mathrm{i}$ tend to $n_\mathrm{w}$.

2.1.3. Electron heat sink in the insulating wall model.

One can notice in the constraint equation (7) that as $n_\mathrm{i}(x,t) \to n_\mathrm{w}$, $\nu_\mathrm{we}(x,t) \to 0$. Therefore, once the bulk ions are absorbed so that $n_\mathrm{i}(x,t) \to n_\mathrm{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 $\mathcal{S}_{\mathrm{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 term of the form

$$\begin{equation} \mathcal{S}_{\mathrm{sink,T}}\left(f_e\right) = - \mathcal{M}_{\mathrm{sink,T}}\nu_{\mathrm{sink,T}}\, \left(f_e - g_{\mathrm{sink,T}}\right). \end{equation} \tag{ 8 }$$

The mask $\mathcal{M}_{\mathrm{sink,T}}$ is similar to the wall mask of equation (6) but restricted to the deeper half portion of each wall, see figure 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 equation (4). The amplitude $\nu_{\mathrm{sink,T}}$ of the restoring force is chosen constant. The target function $g_{\mathrm{sink,T}}$ is a Maxwellian

$$\begin{equation*} g_\mathrm{sink,T}\left(v\right) = \dfrac{n_\mathrm{sink,T}}{\sqrt{2\pi T_\mathrm{sink,T}} } \, \operatorname{exp} \left( -\dfrac{v^2}{2T_\mathrm{sink,T}} \right). \end{equation*}$$

It is characterized by a local density $n_\mathrm{sink,T}(x,t)$ equal to the electron density $n_{\mathrm{sink,T}}(x,t) = \int_{}^{} f_\mathrm{e}(x,v,t) \, \textrm{d}v$ 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_{\mathrm{sink,T}}$ as equal to a fraction of the electron temperature at the wall position: $T_{\mathrm{sink,T}}(t) = \alpha T_\mathrm{e}(x_{\mathrm{w}}, t)$, where $x_{\mathrm{w}}$ is the position where the wall mask defined by equation (6) steps from zero to one. Typically we take $\alpha = 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.

2.1.4. 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 sections 3 and 4 dedicated to the simulations. The location of the sheath entrance and wall position are discussed in these latter sections and specifically in section 3.1. 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:

$$\begin{align} j_\mathrm{tot} & = j_\mathrm{e} + j_\mathrm{i} = q_\mathrm{e} \Gamma_\mathrm{e} + q_\mathrm{i} \Gamma_\mathrm{i}, \nonumber \end{align}$$

where $\Gamma_a$ is the particle flux of species a. On figure 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_\mathrm{sat}^{\mathrm{sh}} = q_\mathrm{i}\Gamma_\mathrm{i}^{\mathrm{sh}}$. 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 $\nu_\mathrm{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.

Zoom In Zoom Out Reset image size

In steady state the normalized charge balance equation is simply written as

$$\begin{align} \partial_{x}\, j_{\text{tot}} & = \sum\limits_{\text{species}}^{} Z_a \int \textrm{d}v \, \mathcal{S}_{\mathrm{sink,n}}\left(f_a\right) = \mathcal{S}_\rho \end{align} \tag{ 9a }$$

$$\begin{align} \mathcal{S}_\rho & = - \mathcal{M}_\mathrm{w} \left(\nu_{\mathrm{w}i} \left( n_i - n_\mathrm{w} \right) - \nu_{\mathrm{w}e} \left( n_e - n_\mathrm{w} \right)\right). \end{align} \tag{ 9b }$$

The right hand side of equation (9a ) is the charge source term $\mathcal{S}_{\rho}$. For the case of electrons and singly charged ions it is given in equation (9b ). In the plasma region where $\mathcal{M}_\mathrm{w} = 0$ the charge source is null and a constant current is solution of equation (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 $\nu_{\mathrm{wi}} = \nu_{\mathrm{we}}$, one obtains $\mathcal{S}_\rho = - \mathcal{M}_\mathrm{w}\nu_{\mathrm{wi}} (n_\mathrm{i} - n_\mathrm{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 $\partial_t f_a = \mathcal{S}(f_a)$ on $\textrm{d}t/2$;
• 2.  
  Solve $\partial_t f_a = -v \partial_x f_a - q_aE/m_a \partial_v f_a$ on $\textrm{d}t$;
• 3.  
  Solve $\partial_t f_a = \mathcal{S}(f_a)$ on $\textrm{d}t/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 $\mathcal{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 $\mathcal{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, figure 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, figure 3. In the insulating wall model the amplitude of the electron sink is defined in equation (7) so that $\mathcal{S}_\rho = 0$. As observed on figure 3, dash dot blue curve, the electrical current remains close to zero.

On figure 4 are compared the steady state charge density $\rho_\mathrm{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 figure 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.

Zoom In Zoom Out Reset image size

2.2.1. Plasma source.

The source/sink term $\mathcal{S}(f_a)$ of equation (3a ) is composed of three distinct entities $\mathcal{S} = \mathcal{S}_{\mathrm{sink, n}} + \mathcal{S}_{\mathrm{sink,T}} + \mathcal{S}_{\mathrm{src}}$. The two sink terms, $\mathcal{S}_{\mathrm{sink,n}}$ and $\mathcal{S}_{\mathrm{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 $\mathcal{S}_{\mathrm{src}}$ for particles and energy is similar to the Gysela source terms [37]. 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 [38]. 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 2.2.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 [37]. This source target distribution chosen here is a Maxwellian, with zero mean velocity, and therefore no momentum transfer, and is written as

$$\begin{align} \mathcal{S}_{\text{src}}\left(f_a\right) & = s_\mathrm{k} \dfrac{\mathcal{M}_{\text{src}}\left(x\right)}{\int_{0}^{L_x} \mathcal{M}_{\text{src}}\left(x\right) \, \textrm{d} x }\times \dfrac{1}{\sqrt{2\pi T_{\text{src}}} } \, \operatorname{exp} \left( -\dfrac{v_a^2}{2T_{\text{src}}} \right). \nonumber \end{align} \tag{ 10 }$$

The mask $\mathcal{M}_{\text{src}}$ plotted on figure 1 localizes the source region in the central part of the plasma. It has the same form as the wall mask, see equation (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_\mathrm{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 equation (5). The ratio chosen in the simulation verifies $T_\text{src} / T_\mathrm{w} \gt 1$. However, to properly describe the distribution functions in the simulations with typically 10² points in velocity space, the chosen value for this ratio is not as large as would be consistent with experimental conditions: we typically take $T_\text{src} / T_\mathrm{w} = 3$ in the simulations.

2.2.2. Collisions.

The collision operator $C(f_a)$ in the simulations of the Boltzmann equation equation (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 $\nu_0^{\star} = 0$ and a weakly collisional plasma, for which $\nu_0^{\star} \to 0^{+}$. The dimensionless control parameter $\nu_0^{\star}$ is defined in equation (13). It is similar to the standard $\nu^{\star}$ parameter used in Gysela [35], 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 [39], 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. Let us remark that the collision operator presented in the following does not account for collisional processes with other velocity dimensions that the one solved by the 1D-1V model. The self-collision operator $\mathcal{C}_{aa}$ we use is inspired from [40]. 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

$$\begin{equation} \mathcal{C}_{aa}\left(f_a\right) = \nu^{\star}_{\mathrm{D}0}\partial_{v_a} \left[ D_v f_{\mathrm{M}a} \, \partial_{v_a} \left( \dfrac{f_a}{f_{\mathrm{M}a}} \right) \right]. \end{equation} \tag{ 11 }$$

The $f_{\mathrm{M}a}$ function is an artificial distribution constructed to be a Maxwellian with mean velocity $u_{\mathrm{M}a}$ and a thermal energy $T_{\mathrm{M}a}$ 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 equation (A4). It is proportional to $\nu^{\star}_{aa} = \sqrt{A_a} Z_a^4 n_a / T_a^{3/2}$, equation (A5). It exhibits the $1/|v|$ velocity dependence of the the Landau and Lenard–Balescu non-linear operators [39]. The control parameter $\nu^*_{\mathrm{D}0}$ is the reference collision frequency normalized by the plasma frequency. It is further discussed below.

The inter species collision operator is adapted from [41] 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

$$\begin{align} \mathcal{C}_{ab}\left(f_a\right) & = \nu^{\star}_{\mathrm{D}0} \left(\frac{2\mathcal{Q}_{ab}}{n_aT_a} H_2 + \frac{\mathcal{R}_{ab}}{n_aT_a^{1/2}} H_1\right) F_{\mathrm{M}a}. \end{align} \tag{ 12a }$$

In this expression H₁ and H₂ stand for the Hermite polynomials of order 1 and 2 that depend on the normalized velocity $w_a = (v-u_a)/\sqrt{T_a}$

$$\begin{align} H_2 & = \tfrac{1}{2} \left(w_a^2-1\right) \hspace{1em};\hspace{1em} H_1 = w_a \nonumber \end{align}$$

and consistently $F_{\mathrm{M}a} = \exp(-\tfrac{1}{2}w_a^2)/\sqrt{2\pi}$ is the normalized Maxwellian distribution. In these expressions $T_a$ and $u_a$ are the normalized thermal energy and mean velocity of the species $a$ distribution function $f_a$ . 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 $\sqrt{A_a}$ factors when addressing species $a$ . The term proportional to $\mathcal{R}_{ab}$ then accounts for momentum exchange between species a and b, hence a friction force, while the term proportional to $\mathcal{Q}_{ab}$ stands for the energy transfer, leading to equipartition. These two terms, expressed in normalized units, are

$$\begin{align} \mathcal{R}_{ab} & = - {n}_a \nu_{ab}^{\star} \left( u_a - u_b \sqrt{{A_b}/{A_a}} \right), \end{align} \tag{ 12b }$$

$$\begin{align} \mathcal{Q}_{ab}^\mathrm{M} & = - 3 n_a \nu_{ab}^{\star} \left( T_a - T_b \right) \dfrac{A_b}{A_b + A_a}- u_a \mathcal{R}_{ab}^\mathrm{M}. \end{align} \tag{ 12c }$$

The species dependent contribution to the collision frequency is $\nu_{ab}^{\star} = \sqrt{2A_a} (1+A_b/A_a) n_b Z_a^2Z_b^2/T_a^{3/2}$, equation (A9). The collisionality is set using $\nu_0^{\star}$ as control parameter, such that $\nu_0^{\star} = \nu^{\star}_{\mathrm{D}0} L_\parallel / \lambda_{\mathrm{D}0}$, see equation (13). It is more convenient than the standard $\nu^*$ parameter for the SOL physics since the parallel transport process does not depend on trapped particles. As for $\nu^{\star}$ the reference parallel scale is the system size in the parallel direction, hence the difference with respect to $\nu_{\mathrm{D}0}$ defined with $\lambda_{\mathrm{D}0}$ as normalizing scale. The control parameter $\nu_0^{\star}$ can be read as the collisional mean free path divided by the reference field line length in the SOL: $L_\parallel \approx \pi 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 $\nu_0^{\star}$ is

$$\begin{align} \nu_0^{\star} = \dfrac{4\sqrt{\pi}}{3}\left(\dfrac{e^2}{4\pi\varepsilon_0}\right)^2 \dfrac{n_0}{T_0^2} \operatorname{log}\Lambda ~L_\parallel = \nu^{\star}_{\mathrm{D}0} \frac{L_\parallel}{\lambda_{\mathrm{D}0}} = \frac{L_\parallel}{L_\textrm{coll}} \end{align} \tag{ 13 }$$

where $\operatorname{log} \Lambda$ is the Coulomb logarithm. The control parameter $\nu_0^{\star}$ does not depend on particle mass.

Collisions in the plasma are a continuous random walk process and thus contribute even for $\nu_0^{\star} \unicode{x2A7D} 1$. Furthermore, their effective amplitude depends on both the control parameter $\nu_0^{\star}$ and on the values taken by the operator. For the self-collision operator equation (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 $\nu_0^{\star}$ is an average value where the low velocity particle have a strong weight. Supra thermal electrons exhibit lower collisionality, at twice the thermal velocity, equation (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⁻³ 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 $\lambda_{\mathrm{D}_0} \approx 1.5\,10^{-5}\,$m. One finds therefore that $L_\parallel = 10^6 \lambda_{\mathrm{D}_0}$ for $L_\parallel \approx {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 1, the size of the simulation domain in units of $\lambda_{\mathrm{D}_0}$ is $L_x = 700$ and $2L_\parallel = 406$. The connection length $2L_\parallel = 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 4.3 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 $\lambda_{\mathrm{D}_0}$. 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 figure 16 it can be seen that all of the particle, momentum and energy sink terms $S_{\mathrm{n}a}$, $S_{\mathrm{u}a}$ and $S_{\mathrm{h}a}$ 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 4.3 did not demonstrate a large dependence of the sheath potential drop on the connection length. The normalized target wall density $n_\mathrm{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_\mathrm{src} = 1$ is the upper value that is achieved when the plasma reaches the source temperature. In the wall, the smaller value $T_\mathrm{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_\mathrm{w}$. The appropriate resolution determined by the velocity mesh is found to be more important than the actual value of $T_\mathrm{w}$ provided it is smaller than $T_\mathrm{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 $\lambda_{\mathrm{D}_0}$ and that of $\nu_0^{\star}$, one obtains $L_\textrm{coll} = \lambda_{\mathrm{D}_0}^2/\lambda_{\mathrm{L}_0} (3\sqrt{\pi}/\operatorname{log} \Lambda)$, where the Landau scale $\lambda_{\mathrm{L}_0} = (e^2/(4\pi\varepsilon_0))/T_0$ is the distance of closest approach. It is also interesting to note that $\lambda_{\mathrm{D}_0}/L_\textrm{coll} = 3\sqrt{\pi}\operatorname{log} \Lambda / \Lambda$ given $\Lambda = \lambda_{\mathrm{D}_0} / \lambda_{\mathrm{L}_0}$, so that for standard plasmas with $\Lambda \gg 1$, one finds that the physics on a Debye scale is always weakly collisional $L_\textrm{coll} \gg\lambda_{\mathrm{D}_0}$. In practice, with $\lambda_{\mathrm{L}_0} \approx \,10^{-11}\,$m and $\lambda_{\mathrm{D}_0} \approx 1.5\,10^{-5}\,$m one obtains $\Lambda = 1.5\,10^{6}$, and therefore $\nu^{\star}_{\mathrm{D}0} = \lambda_{\mathrm{D}_0}/L_\textrm{coll} \approx 5\,10^{-5}$, so that for $L_\parallel = 10^6 \lambda_{\mathrm{D}_0}$, $\nu_0^{\star} \approx 50$. However, for the reduced box size chosen for the simulation, $L_\parallel \approx 200 \lambda_{\mathrm{D}_0}$, one finds $\nu_0^{\star} \approx 10^{-2}$. The large value considered in the simulations $\nu_0^{\star} = 50$ (cf table 1) is a compromise between realistic physics that require $\nu_0^{\star} \approx 10^{-2}$ and numerical constraints. Indeed, the smaller $\nu_0^{\star}$, 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 $\Delta 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/\sqrt{\nu_0^{\star}}$. The cost to use the realistic value of $\nu_0^{\star}$ 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\times 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/\nu_0^{\star}$. However, the choice made for $\nu_0^{\star}$ is such that the collisionality in the bulk plasma is unchanged when reducing the system size from $10^6 \lambda_{\mathrm{D}_0}$ in the reference experimental conditions to $200 \lambda_{\mathrm{D}_0}$ 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 $\nu^{\star}_{\mathrm{D}0} = 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 $\nu_0^{\star}$.

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 [42]. 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 2.2.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 $\int f_\mathrm{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 $\rho_\mathrm{c}$ with respect to the potential $\textrm{d} \rho_\mathrm{c} / \textrm{d}\phi$ 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 figure 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_\mathrm{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 figure 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 equation (4). Also, it appears that the energy of the plasma is conserved up to this position, as discussed later in section 3.3. Indeed at this location the normalized total current $j_\mathrm{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\in [L_x/2, L_x]$, so that the wall is located at $x \unicode{x2A7E} x_\mathrm{w}$. The definition of $x_\mathrm{sh}$ and $x_\mathrm{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

$$\begin{align*}\rho_\mathrm{c}\left(x_\mathrm{w}\right) = \max\limits_{x \unicode{x2A7E} L_x/2} \rho_\mathrm{c}.\end{align*}$$

The sheath entrance location is defined as the intersection between the two linear fits defined as

$$\begin{align*} \mathcal{D}_1:\ \ & \mathcal{L}\left(\operatorname{log}\left|\dfrac{\textrm{d}\rho_\mathrm{c}}{\textrm{d}\phi}\right|\right) \quad \text{with} \quad x\in\left[L_x/2,\, x_\mathrm{w}\right]\, \\ & \text{and}\quad \alpha \bigg\lVert\dfrac{\textrm{d}\rho_\mathrm{c}}{\textrm{d}\phi}\bigg\rVert_\infty\! \unicode{x2A7D} \left|\dfrac{\textrm{d}\rho_\mathrm{c}}{\textrm{d}\phi}\right| \unicode{x2A7D} \bigg\lVert\dfrac{\textrm{d}\rho_\mathrm{c}}{\textrm{d}\phi}\bigg\rVert_\infty, \end{align*}$$

where $\mathcal{D}_1$ is the label for the linear fit, $\mathcal{L}$ is a linear fit operator and $\alpha \in [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 $\lVert \, \cdot \rVert_\infty$ notation refers to the maximum of the considered quantity taken on the $[L_x/2,\, x_\mathrm{w}]$ interval. The second linear fit is defined as

$$\begin{align*} \mathcal{D}_2:\ \ & \mathcal{L}\left(\operatorname{log}\left|\dfrac{\textrm{d}\rho_\mathrm{c}}{\textrm{d}\phi}\right|\right) \quad \text{with} \quad x\in\left[L_x/2,\, x_\mathrm{w}\right]\, \\ & \text{and}\quad 0 \unicode{x2A7D} \left|\dfrac{\textrm{d}\rho_\mathrm{c}}{\textrm{d}\phi}\right| \unicode{x2A7D} \alpha \bigg\lVert\dfrac{\textrm{d}\rho_\mathrm{c}}{\textrm{d}\phi}\bigg\rVert_\infty. \end{align*}$$

The first linear fit $\mathcal{D}_1$ corresponds to performing a linear fit of $\operatorname{log}|\textrm{d}\rho_\mathrm{c}/\textrm{d}\phi|$ in the positively charged region, whereas the second linear fit $\mathcal{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 figure 4.

We define the electric potential drop in the sheath region as $\Delta \phi_\mathrm{sh} = \phi(x_\mathrm{w}) - \phi(x_\mathrm{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_\mathrm{c}$. Using the conservation of energy in the sheath for electrons, we obtain for this critical velocity the expression in normalized units

$$\begin{equation} v_\mathrm{c} = \sqrt{-2\Delta \phi_\mathrm{sh}}. \end{equation} \tag{ 14 }$$

On the top left and bottom left graphs of figure 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 equation (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

$$\begin{align*}H_a = \dfrac{1}{2} m_a v^2 + e_a \phi.\end{align*}$$

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 figure 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 \unicode{x2A7E} 0$. The $x<0$ region is the symmetrical of the $x \unicode{x2A7E} 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 \unicode{x2A7E} v_\mathrm{c}$ at the sheath entrance overcome the potential barrier of the sheath and penetrate the absorbing wall region. On the contrary, electrons with $v\unicode{x2A7D} v_\mathrm{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 equation (4) does not absorb perfectly every electron that enters the wall region, as the magnitude $\nu_\mathrm{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 equation (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 plasma-wall 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 4.1. 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 $\left( 1 \right)$ 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.

Zoom In Zoom Out Reset image size

The trapped electrons for which $v \unicode{x2A7D} v_\mathrm{c}$ do not contribute to the electron flux at the sheath entrance. To illustrate this, on figure 6 we plot the partial electron flux $\Gamma_\mathrm{e}^{\star}$ against the velocity variable. This partial flux is defined as

$$\begin{equation} \Gamma_\mathrm{e}^{\star}\left(v\right) = \dfrac{1}{\Gamma_\mathrm{e}^{\,{\mathrm{sh}}}}\int_{-v}^{v} \textrm{d}v^{^{\prime}} \, v^{^{\prime}} f_\mathrm{e}^{\mathrm{sh}}\left(v^{^{\prime}}\right), \end{equation} \tag{ 15 }$$

where $\Gamma_\mathrm{e}^{\mathrm{sh}}$ is the total electron flux at the sheath entrance. The partial electron flux is negligible for velocities lower than the critical velocity $v \unicode{x2A7D} v_\mathrm{c}$ and starts to become non zero when $v \unicode{x2A7E} v_\mathrm{c}$. The electron flux at the sheath entrance is therefore essentially carried by the high velocity $\left| v \right| \gt v_\mathrm{c}$ component of the electron distribution function.

Zoom In Zoom Out Reset image size

On the same graph we also plot the negative velocity part of the electron distribution function at the sheath entrance $f_\mathrm{e}^{\,\,{\mathrm{sh-}}}(\left| v \right| ) \equiv f_\mathrm{e}^{\,\,{\mathrm{sh}}}(-\left| v \right| )$ on top of its positive velocity part $f_\mathrm{e}^{\,\,{\mathrm{sh+}}}(\left| v \right| ) \equiv f_\mathrm{e}^{\,\,{\mathrm{sh}}}(\left| v \right| )$. There exists a clear imbalance between these two parts of the electron distribution function for velocities higher than the critical velocity $v_\mathrm{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_\mathrm{e}^{\,\,{\mathrm{sh-}}}(\left| v \right| \unicode{x2A7D} -v_\mathrm{c}) \neq 0$. Their origin is twofold. First, the wall is not perfectly absorbing: the magnitude $\nu_\mathrm{w}$ of the particle sink term of equation (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 \unicode{x2A7D} -v_\mathrm{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_\mathrm{e}^{\,\mathrm{sh-}}(\left| v \right| \unicode{x2A7D} -v_\mathrm{c})$ would be strictly zero.

The ion distribution function is shown in figure 7. 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 \unicode{x2A7D} 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\approx 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.

Zoom In Zoom Out Reset image size

On figure 7 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_{\mathrm{Me}}$ of the intra species collision operator, see equation (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, plotted as the blue dotted curve on figure 7. 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 $\left| v \right|^{-3}$ at large velocities. The effect of the collisionality on the shape of the distribution functions is investigated more extensively in appendix B.

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⁻² in steady state [43]. The so-called sheath heat transmission coefficient $\gamma$ is defined by the relation

$$\begin{align*}Q = \gamma \Gamma T_e, \end{align*}$$

where Γ, Q and $T_\mathrm{e}$ are respectively the particle flux, energy flux and electron temperature taken at the plasma-wall 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

$$\begin{align*}Q_\mathrm{e} = \gamma_\mathrm{e} \Gamma T_\mathrm{e},\end{align*}$$

and

$$\begin{align*}Q_\mathrm{i} = \gamma_\mathrm{i} \Gamma T_\mathrm{i}.\end{align*}$$

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

$$\begin{equation} Q_a = \int_{-\infty}^{\infty} \textrm{d}v \, \dfrac{1}{2} m_a v^3 f_a. \end{equation} \tag{ 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 figure 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 $\gamma$ 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 3.1. 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.

Zoom In Zoom Out Reset image size

In appendix E we detail a model developed in [34] which is commonly used to predict a value for the sheath heat transmission coefficient. The prediction yields $\gamma_\mathrm{e}^\mathrm{pred.} = 2.7$ for electrons, and $\gamma_\mathrm{i}^\mathrm{pred.} = 2.5$ for ions. As explained in appendix E, when considering a 3D geometry and a realistic ion-to-electron mass ratio the predictions are different: $\gamma_\mathrm{e}^\mathrm{3D,pred.} = 5$ and $\gamma_\mathrm{i}^\mathrm{3D,pred.} = 3.5$. As shown on figure 8, the ion sheath heat transmission factor is approximately 60% higher than the predicted value and the electron sheath heat transmission factor is approximately 35% higher than the predicted value. Thus, the total heat transmission factor obtained from the kinetic simulation is approximately 50% higher than the prediction. The discrepancy between the heat transmission factor from the simulation 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 figure 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_\mathrm{i}^\mathrm{sh}$ 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. The discrepancy between the electron sheath heat transmission factor found in the simulation and its prediction can be explained by the following point. The prediction assumes that the electron distribution function is a truncated Maxwellian. On figure 7 it can be seen that the electron distribution function resembles a sum of two Maxwellian, with the positive high velocity tail of the distribution function having a temperature close to the one of the source, higher than the plasma temperature. Therefore, the high velocity electrons deposit more power on the wall than predicted.

The mass ratio between ions and electrons $m_\mathrm{i} / m_\mathrm{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_\mathrm{i} / m_\mathrm{e} \approx 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 $\left\{40, 100, 400, 1835, 3600 \right\}$ on the sheath potential drop. The $m_\textrm{i} / m_\textrm{e} = 1836$ value corresponds to a hydrogen plasma, while $m_\textrm{i} / m_\textrm{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 2.

As explained in section 3.2 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 equation (17), which assumes a perfectly absorbing wall. We therefore adapt the prediction of the potential drop given by equation (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

$$\begin{eqnarray} &&\Delta \phi_{\,{\mathrm{sh}}}^\mathrm{ref.} = \dfrac{T_\mathrm{e}^{\,{\mathrm{sh}}}}{e} \operatorname{log} \left( \sqrt{2\pi \dfrac{m_\mathrm{e}}{m_\mathrm{i}} M^2\left( 1 + \dfrac{T_\mathrm{i}^{\,{\mathrm{sh}}}}{T_\mathrm{e}^{\,{\mathrm{sh}}}} \right)} \times \left( 1 + \dfrac{\Gamma_\mathrm{e}^{-}}{\Gamma_\mathrm{i}^{\mathrm{sh}}} \right) \right), \end{eqnarray} \tag{ 20 }$$

where $\Gamma_\mathrm{e}^{-}$ is the flux at the sheath entrance of high velocity electrons coming back from the wall

$$\begin{equation} \Gamma_\mathrm{e}^{-} = -\int_{-\infty}^{-v_\mathrm{c}} \textrm{d}v\, vf_\mathrm{e}^{\,\,{\mathrm{sh}}}\left(v\right). \end{equation} \tag{ 21 }$$

Using this refined prediction we can define in a similar fashion as in equation (18) an effective potential drop as

$$\begin{align} \overline{\Delta }\phi _{\mathrm{sh}}^{\mathrm{ref.}} \equiv \dfrac{e\Delta \phi _\mathrm{sh}}{T_\mathrm{e}^{\,{\mathrm{sh}}}}- \operatorname{log} \left( \sqrt{2\pi M^{\,2}\left( 1 + \dfrac{T_\mathrm{i}^{\,{\mathrm{sh}}}}{T_\mathrm{e}^{\,{\mathrm{sh}}}}\right)} \times \left( 1 + \dfrac{\Gamma_\mathrm{e}^{-}}{\Gamma_\mathrm{i}^{\mathrm{sh}}} \right) \right). \end{align} \tag{ 22 }$$

This refined potential drop is represented on figure 9 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 equation (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.

Zoom In Zoom Out Reset image size

The effective potential drop without refinement given by equation (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 figure 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.

It is well admitted that the establishment of the H-mode transport barrier at the last closed flux surface in tokamak plasmas is closely linked to the shearing of the electric field at this location [15, 44]. 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 $\Delta \phi _\mathrm{sh}^{\mathrm{pred.}} = \Lambda T_\mathrm{e}/e$. Taking into account the electron flux coming back from the wall at the sheath entrance the Λ term is given by equation (20) as

$$\begin{equation} \Lambda = \operatorname{log} \left( \sqrt{2\pi \dfrac{m_\mathrm{e}}{m_\mathrm{i}} M^2\left( 1 + \dfrac{T_\mathrm{i}^{\,{\mathrm{sh}}}}{T_\mathrm{e}^{\,{\mathrm{sh}}}} \right) } \left( 1 + \dfrac{\Gamma_\mathrm{e}^{-}}{\Gamma_\mathrm{i}^{\mathrm{sh}}} \right) \right). \end{equation} \tag{ 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 equation (10). We use temperatures in the range $T_{\mathrm{src},\mathrm{e}} \, \in \, \left\{0.5,\, 0.75,\, 1,\, 1.25,\, 1.5 \right\}$. The temperature of the ion source remains constant in all of these simulations, we use a mass ratio value of $m_\mathrm{i} / m_\mathrm{e} = 400$ and a plasma length of $L_\parallel = 300\lambda_{\mathrm{D}_0}$. Figure 10 shows the quantity $e\Delta \phi_\mathrm{sh} /\Lambda$ for both insulating and conducting wall cases, with the Λ term defined in equation (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\Delta \phi_\mathrm{sh} / \Lambda$ 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 equation (20) we write the equality of ion and electrons fluxes at the sheath entrance as

$$\begin{align*}\Gamma_\mathrm{i}^{\mathrm{sh}} = \Gamma_\mathrm{e}^{+\mathrm{pred.}} - \Gamma_\mathrm{e}^{-},\end{align*}$$

where

$$\begin{equation} \Gamma_\mathrm{e}^{+\mathrm{pred.}} = \int_{v_\mathrm{c}}^{+\infty} \textrm{d}v\, v n_e^{\mathrm{sh}}\sqrt{\dfrac{m_\mathrm{e}}{2\pi T^{\,{\mathrm{sh}}}_\mathrm{e}}} e^{- \dfrac{m_\mathrm{e}v^2}{2T^{\,{\mathrm{sh}}}_\mathrm{e}}} \end{equation} \tag{ 24 }$$

is the predicted forward flux of electrons with velocities $v \unicode{x2A7E} v_\mathrm{c}$. The actual forward flux of electrons with velocities $v \unicode{x2A7E} v_\mathrm{c}$ from the simulations is rather expressed as

$$\begin{equation} \Gamma_\mathrm{e}^{+} = \int_{v_\mathrm{c}}^{+\infty} \textrm{d}v \, vf_\mathrm{e}^{\,{\mathrm{sh}}}\left(v\right). \end{equation} \tag{ 25 }$$

Note here that $f_\mathrm{e}^{\,{\mathrm{sh}}}$ 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 flux—i.e. if $\Gamma_\mathrm{e}^{+} = \Gamma_\mathrm{e}^{+\mathrm{pred.}}$.

Zoom In Zoom Out Reset image size

As shown on figure 7 the electron distribution function is not a Maxwellian for velocities higher than $v_\mathrm{c}$ but resembles more to a superposition of two Maxwellians of different temperatures, as discussed in detail in section 3.2. Therefore, the forward flux of the simulations departs from the predicted forward flux: $\Gamma_\mathrm{e}^{+} \neq \Gamma_\mathrm{e}^{+\mathrm{pred.}}$. This departure can be quantified by noticing that $f_\mathrm{e}^{\,{\mathrm{sh}}}(v \unicode{x2A7E} v_\mathrm{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 equation (24) so that it matches the forward flux of the simulation, i.e. so that $\Gamma_\mathrm{e}^{+} = \Gamma_\mathrm{e}^{+\mathrm{pred.}}$. This corresponds to taking the temperature of the best Maxwellian fit for the $v \unicode{x2A7E} v_\mathrm{c}$ part of the electron distribution function at the sheath entrance. As discussed in section 3.2, this temperature is close to the one of the particle source term $S_\mathrm{src}$—it would be equal to it in the absence of collisions. By doing this we obtain the appropriate set of points on figure 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.

The $L_\parallel$ 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_\parallel$ will have a significant effect on the plasma-wall interaction physics. The prediction of equation (17) shows no dependence of the potential drop on $L_\parallel$. 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_\parallel \, \in \, \left\{300, 600, 1200, 2400 \right\}$. The width of the walls does not change between these simulations: the fraction of the total length of the simulation domain $L_x$ that accounts for the wall region decreases when $L_\parallel$ increases. The particle source is kept constant in these simulations. Also, the collisionality of the plasma $\nu_0^{\star}$ defined in equation (13) has been taken constant in the simulations. On figure 11 we plot the effective potential drop defined in equation (18) both conductive and insulating walls. The same scale is used as in the mass ratio scan of figure 9. We see that the values from the simulations stay close to the constant predicted value of potential drop, with a weak dependence in $L_\parallel$. When using the refined potential drop that is adapted for the electron temperature as explained in section 4.2 we see almost no variation of this quantity with the parallel length.

Zoom In Zoom Out Reset image size

The logical sheath boundary condition has first been introduced in [45]. In this approach, the quasineutral plasma only is modeled. The positively charged region that would form at the plasma boundary is not resolved, see figure 12. 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 particle-in-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 $\left[ -L_\parallel/2,\, L_\parallel/2 \right]$ region. The $z = \pm L_\parallel/2$ coordinate denotes the plasma boundary. In a particle-in-cell code, a pusher method moves the particles from time t to $t+\Delta t$. Δt stands for the simulation time step. The logical sheath condition consists in the following sequential steps, see figure 12 as well:

• 1.  
  Starting from the particle positions at time t, apply the pusher method to compute their new position at time $t+ \Delta t$.
• 2.  
  As a result of step 1. particles may have crossed the simulation boundary at $z = L_\parallel/2$. Count the number of outgoing electrons $N_\mathrm{e}$ and ions $N_\mathrm{i}$. Then, two situations can occur:
• 3.  
  More electrons than ions have crossed the boundary, i.e. $N_\mathrm{e} \gt N_\mathrm{i}$. In this case:

  • (a)  
    Order all $N_\mathrm{e}$ electrons from slowest to fastest.
  • (b)  
    Remove all $N_\mathrm{i}$ ions from the simulation, and the fastest $N_\mathrm{i}$ electrons. The velocity of the slowest electron that is absorbed in this step defines a cutoff velocity $v_\mathrm{c}$.
  • (c)  
    The remaining $N_\mathrm{e} - N_\mathrm{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.
• 4.  
  More ions than electrons have crossed the boundary, i.e. $N_\mathrm{i} \gt N_\mathrm{e}$. In this case:

  • (a)  
    All these $N_\mathrm{i}$ ions and $N_\mathrm{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_\mathrm{i}-N_\mathrm{e}$ is attached to the plasma boundary position at $z = L_\parallel/2$ for the next time step only.

Zoom In Zoom Out Reset image size

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 $\Delta \phi_\mathrm{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_\mathrm{c}$. Using conservation of energy, it is therefore expressed in dimensional units as

$$\begin{equation*} \Delta \phi_{\mathrm{sh}} = - \dfrac{m_\mathrm{e}}{2e}v_\mathrm{c}^2. \end{equation*}$$

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 [46, 47] 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_\parallel/2$ reads as follows

$$\begin{equation*} \int_{0}^{\infty} \textrm{d}v\, vf_\mathrm{i}\left(x, y, L_\parallel/2, v, \mu\right) = \int_{v_\mathrm{c}}^{\infty} \textrm{d}v\, vf_\mathrm{e}\left(x, y,L_\parallel/2, v, \mu\right). \end{equation*}$$

This equality determines a cutoff velocity $v_\mathrm{c}$ when situation 4. is not encountered, see above. At the plasma boundary position $z = L_\parallel/2$, all outgoing ions are removed from the simulation at the next time step. The fastest electrons that have velocities larger than $v_\mathrm{c}$ are removed as well, while the slowest $v\lt v_\mathrm{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_\mathrm{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 [48].

Another way of including the plasma-wall physics at the plasma boundary without resolving the sheath scale is presented in [49]. 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 $\phi_\mathrm{sh}$ determines a proportion of reflected electrons. In the meantime, outgoing ions leave the simulation freely. The authors of [49] 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 = \pm L_\parallel/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_\parallel/2$. The conducting sheath algorithm consists in the following:

• 1.  
  The potential at the plasma boundary $\phi_\mathrm{sh}(x,y) = \phi(x,y,L_\parallel/2)$ is computed from the quasineutrality equation. This equation can be written as presented in [35] as

  $$\begin{equation} - \nabla_{\perp} \left(\dfrac{m_\mathrm{i} n_\mathrm{i}}{eB^2}\nabla_\perp \phi\right) + e^2\dfrac{\left\langle n_\mathrm{e}\right\rangle _\mathrm{FS}}{\left\langle T_\mathrm{e}\right\rangle _\mathrm{FS}}\left(\phi - \left\langle \phi \right\rangle _\mathrm{FS}\right) = \sigma_g, \end{equation} \tag{ 26 }$$

  where $n_\mathrm{i}$ (resp. $n_\mathrm{e}$) is the ion (resp. electron) charge density and $T_\mathrm{e}$ stands for the electron temperature. The bracket notation $\left\langle\ldots\right\rangle_\mathrm{FS}$ stands for the average over a given magnetic surface. The guiding center charge density is defined as $\sigma_g = e(Z_\mathrm{i} n_\mathrm{Gi} - n_\mathrm{e})$, where $n_\mathrm{Gi}$ is the ion gyrocenter density and $Z_\mathrm{i}$ is the normalized ion charge (equal to +1 for hydrogen isotopes). $\nabla_{\perp}$ refers to the gradient in the plane perpendicular to the magnetic field.
• 2.  
  When the potential $\phi_\mathrm{sh}(x,y)$ is positive it defines a cutoff velocity which is such that $\tfrac{1}{2}m_\mathrm{e} v_\mathrm{c}^2 = e \phi_{\mathrm{sh}}$. At the considered wall position, all slow electrons that have velocities v such that $0 \unicode{x2A7D} v \unicode{x2A7D} v_\mathrm{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 $\Delta v$. The value of the distribution function $f_\mathrm{e}(x,y,z,v_i, \mu)$ accounts for the density of particles that have velocities spanning the $\left[ v_i - \Delta v / 2, \, v_i + \Delta v /2 \right]$ interval. In general the cutoff velocity $v_\mathrm{c}$ does not lie on a specific point of the velocity grid, see figure 13. 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_\mathrm{c} \, \in \, \left[ v_j - \Delta v / 2, \, v_j + \Delta v /2 \right]$. Therefore, for all velocities $0\unicode{x2A7D} v_i \lt v_j$ the value of the distribution function $f_\mathrm{e}(v_i)$ is copied to a ghost cell $v_i^{^{\prime}}$ which is such that $v_i^{^{\prime}} = -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_\mathrm{e}(x,y,L_\parallel/2,v_j, \mu)$ is copied to the corresponding ghost cell $v_j^{^{\prime}} = -v_j$ and scaled by a factor c defined as

$$\begin{equation*} c\left(x,y\right) = \dfrac{\int_{v_{j} - \Delta v/2}^{v_\mathrm{c}}\textrm{d}v \int_{0}^{\mu_\text{max}} \textrm{d} \mu \, f_\mathrm{e}\left(x,y, L_\parallel/2, v, \mu\right) }{\int_{v_j - \Delta v /2}^{v_{j} + \Delta v /2} \textrm{d}v \int_{0}^{\mu_\text{max}} \textrm{d} \mu \,f_\mathrm{e}\left(x,y,L_\parallel/2,v, \mu\right)}. \end{equation*}$$

It ensures that all outgoing electrons verifying $\tfrac{1}{2} m_\mathrm{e} v^2 \unicode{x2A7D} e \phi_\mathrm{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.

Zoom In Zoom Out Reset image size

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_\parallel/2, v_i, \mu)$ for all $v_i\unicode{x2A7D} 0$. Since no ion comes back from the wall, we simply set

$$\begin{align*}f_\mathrm{i}\left(x,y,L_\parallel/2, v_i, \mu\right) = 0 \quad \forall v_i \unicode{x2A7D} 0.\end{align*}$$

For the electrons we use the ghost cell values as a boundary condition, i.e.

$$\begin{align*}f_\mathrm{e}\left(x,y,L_\parallel/2, v_i, \mu\right) = f_\mathrm{e}\left(x,y,L_\parallel/2, v_i^{^{\prime}}, \mu\right) \quad \forall v_i \unicode{x2A7D} 0.\end{align*}$$

The $\tfrac{1}{2}m_\mathrm{e} v_c^2 = e \phi_\mathrm{sh}$ equation defines a cutoff velocity for electrons only if the potential at the wall position $\phi_\mathrm{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 $\sigma_\mathrm{g}$ equal to zero everywhere. Since $\phi_\mathrm{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 $\sigma_\mathrm{g}(L_\parallel/2) \gt 0$. In the quasineutrality equation (26), at first order we have $- \nabla_{\perp} \left(m_\mathrm{i} n_\mathrm{i}/(eB^2)\nabla_\perp \phi\right) \ll \phi$. 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 $\left(\phi - \left\langle \phi \right\rangle_\mathrm{FS}\right) \approx \sigma_\mathrm{g} \left\langle T_\mathrm{e}\right\rangle_\mathrm{FS}/(e^2\left\langle n_\mathrm{e}\right\rangle_\mathrm{FS}) \gt 0$. Thus, the potential $\phi(L_\parallel/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 fluctuations—due 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 quasi-neutral 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

$$\begin{align} \left(\frac{\lambda_{\mathrm{D}_0}}{L_\parallel}\right)^2 \tilde\nabla_\parallel^2 \tilde\phi + \left(\frac{\lambda_{\mathrm{D}_0}}{L_\perp}\right)^2 \tilde\nabla_\perp^2 \tilde\phi = \frac{n_\mathrm{e}}{n_0} - \frac{Z_\mathrm{i}}{n_0}\left(n_\mathrm{Gi} + n_\mathrm{pol}\right), \end{align} \tag{ 27 }$$

with $\lambda_{\mathrm{D}_0}$ the Debye length at density $n_0$ and temperature $T_0$ , $L_\parallel$ and $L_\perp$ 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 $\tilde\phi = e\phi/T_0$ is the dimensionless potential. $n_\mathrm{e} \approx \int f_\mathrm{Ge}\ \textrm{d}^3v$ is the electron density and $n_\mathrm{Gi} = \int \mathcal{J}_{\mathrm{i}}.f_\mathrm{Gi}\ \textrm{d}^3v$ the ion gyrocenter density where $\mathcal{J}_{\mathrm{i}}$ is the gyro-average operator. The gyro-average operator of a given function g written $\mathcal{J}_{\mathrm{i}} g$ is defined as the average of g over a cyclotron motion as

$$\begin{align*}\mathcal{J}_{\mathrm{i}}\, g = \oint_0^{2\pi} g \, \dfrac{\textrm{d}\varphi_\mathrm{c}}{2\pi},\end{align*}$$

where $\varphi_\mathrm{c}$ is the cyclotron phase. The ion polarization density $n_\mathrm{pol}$ reads, in the long wavelength limit $L_\perp\gg \rho_\mathrm{i}$ (with $\rho_\mathrm{i}$ the ion gyro-radius): $n_\mathrm{pol} \approx {\pmb\nabla}_\perp. [(m_\mathrm{i}n_\mathrm{e}/eB^2) {\pmb\nabla}_\perp\phi]$. In any case, the scaling $L_\perp \gg \lambda_{\mathrm{D}_0}$ 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 equation (27) is neglected in usual orderings. For consistency, it readily appears that the first term can also be neglected to describe turbulence, since $L_\parallel \approx qR \gg \rho_\mathrm{i} \gg \lambda_{\mathrm{D}_0}$. Notable exceptions to this ordering is in the vicinity of the sheath, where $L_\parallel \approx \lambda_{\mathrm{D}_0}$, and of the wall when $L_\perp \lt \rho_\mathrm{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_\mathrm{e} = Z_\mathrm{i}(n_\mathrm{Gi} + n_\mathrm{pol})$. The gyrokinetic quasi-neutrality equation thus only involves at leading order the transverse operator ${\pmb\nabla}_\perp$, which enters through the expression of $n_\mathrm{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 $\nu_{\mathrm{w}a}$ (cf equation (4)), the electron distribution function $f_\mathrm{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 figures 6 and 7 and associated discussion). Also, the ion distribution function $f_\mathrm{i}$ is shifted (maximum at $v\neq0$ ) and vanishes in the negative half of velocity space, see figure 7.

Combining these observations with the method proposed for a conducting sheath (cf section 5.2 and [49]), 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 electron 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 $\mathcal{M}_\textrm{w}(r,\theta)$, the subscript '$\textrm{w}$' standing for 'wall': $\mathcal{M}_\textrm{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_\mathrm{Gi}$ using a Krook term similar to the one described in equation (4) of the form
  $$\begin{align*}\mathcal{S}_\mathrm{Gi} = -\mathcal{M}_\mathrm{w}\nu \left(f_\mathrm{Gi} - f_\mathrm{Gi}^\mathrm{T} \right),\end{align*}$$
  where the target gyrocenters distribution function $f_\mathrm{Gi}^\mathrm{T}$ is
  $$\begin{align*}f_\mathrm{Gi}^\mathrm{T}\left(r, \theta, v_\parallel, \mu, t\right) = f_\mathrm{Gi}\left(r, \theta, v_\parallel, \mu, t\right)\times\operatorname{H}\left(v_\parallel\right),\end{align*}$$
  and ν is a constant coefficient that defines the Krook term magnitude. Here $\operatorname{H}(v_\parallel)$ is the Heaviside function such that $\operatorname{H}(v_\parallel \lt 0) = 0$ and $\operatorname{H}(v_\parallel \unicode{x2A7E} 0) = 1$. The fact that the target ion distribution function $f_\mathrm{Gi}^\mathrm{T}$ 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_\mathrm{Ge}$ is relaxed towards a target function truncated for velocities $v \unicode{x2A7D} -v_{\parallel \mathrm{c}}$, in the spirit of the results found in section 3.2. 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 $\left\langle f_\mathrm{Ge}\right\rangle_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_\mathrm{Ge} = \left\langle\, f_\mathrm{Ge}\right\rangle_Z + \delta f_\mathrm{Ge}$, where $\delta f_\mathrm{Ge} = f_\mathrm{Ge} - \left\langle\, f_\mathrm{Ge}\right\rangle_Z$.
  • 2.  
    Define a critical velocity $v_{\parallel \mathrm{c}}$ such that on average the electron flux due to electrons for which $v \unicode{x2A7E} -v_{\parallel \mathrm{c}}$ balances the parallel ion flux on the considered limiter portion:

    $$\begin{align*}\hspace{1.5cm} \int_{-v_{\parallel \mathrm{c}}}^{+\infty} v_\parallel \left\langle\, f_\mathrm{Ge}\right\rangle _Z \textrm{d}^3v = \left\langle \Gamma_{\parallel \mathrm{i}}\right\rangle _Z.\end{align*}$$
  • 3.  
    Using this critical velocity, the average of the electron distribution function is relaxed towards the target function

    $$\begin{align*}\hspace{1.5cm} f_\mathrm{Ge}^{\,\mathrm{T}}\left(v_\parallel, t\right) = \left\langle \,f_\mathrm{Ge}\right\rangle_Z \left(1-\operatorname{H}\left(v_\parallel + v_{\parallel \mathrm{c}}\right)\right)\end{align*}$$

    through the evolution equation

    $$\begin{align*}\hspace{1.5cm} \partial_t \left\langle\, f_\mathrm{Ge}\right\rangle _Z = -\mathcal{M}_\mathrm{w}\nu \left( \left\langle f_\mathrm{Ge}\right\rangle_Z - f_\mathrm{Ge}^\mathrm{T}\right).\end{align*}$$
  • 4.  
    Solving the latter equation yields an updated value of the averaged distribution function $\left\langle\, f_\mathrm{Ge}\right\rangle _Z^\mathrm{new}$ that is used to reconstruct the electron distribution function for the next timestep:

    $$\begin{align*}\hspace{1.5cm} f_\mathrm{Ge}\left(r, \theta, v_\parallel, \mu, t+\textrm{d}t\right) = \left\langle\, f_\mathrm{Ge}\right\rangle_Z^\mathrm{new} + C\times \delta f_\mathrm{Ge},\end{align*}$$

    where the $C(r, \theta, \varphi)$ coefficient is adapted such that quasineutrality is ensured locally, i.e.

    $$\begin{align*}\hspace{1.5cm} \int f_\mathrm{Ge} \, \textrm{d}^3v = \int f_\mathrm{Gi} \, \textrm{d}^3v.\end{align*}$$

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 5.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 $\mathcal{M}$. With these conditions, two constraints are satisfied: (i) the relaxed ion distribution function $f_\mathrm{Gi}^\mathrm{T}$ 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 non-vanishing 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.

The sheath potential drop prediction detailed below is taken from [34]. Let us consider a perfectly absorbing wall at position $x_\mathrm{w}$, bounding a semi infinite plasma located at $x\unicode{x2A7D} x_\mathrm{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_\mathrm{sh}$ refers to the sheath entrance position. The flux of ions at the sheath entrance can be written as

$$\begin{equation} \Gamma_{\mathrm{i}}^{\mathrm{sh}} = n_\mathrm{i}^{\mathrm{sh}} M c_\mathrm{s}. \end{equation} \tag{ F1 }$$

The Mach number M and sound speed velocity $c_\mathrm{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_\mathrm{e} = \tfrac{1}{2} m_\mathrm{e} v^2 - e \phi$. 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

$$\begin{align*}\dfrac{1}{2}m_\mathrm{e} v\left(x_\mathrm{sh}\right)^2 - e\phi\left(x_\mathrm{sh}\right) = \dfrac{1}{2}m_\mathrm{e} v\left(x_\mathrm{w}\right)^2 - e\phi\left(x_\mathrm{w}\right),\end{align*}$$

or equivalently

$$\begin{align*}v\left(x_\mathrm{w}\right)^2 = v\left(x_\mathrm{sh}\right)^2 + \dfrac{2e\Delta\phi_\mathrm{sh}}{m_\mathrm{e}}.\end{align*}$$

Here we see that a necessary condition for the electrons to reach the wall is $v(x_\mathrm{sh}) \unicode{x2A7E} v_\mathrm{c} = \sqrt{-2e\Delta\phi_\mathrm{sh}/m_\mathrm{e}}$. The electrons that are faster than the velocity $v_\mathrm{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_{\mathrm{sh}})\lt -v_\mathrm{c}$ at the sheath entrance. We therefore assume that the electron distribution function at the sheath entrance is a Maxwellian truncated for velocities $v\unicode{x2A7D} -v_\mathrm{c}$. This truncated Maxwellian is characterized by a density $\overline{n}_\mathrm{e}$ and a temperature $T_\mathrm{e}^{\mathrm{sh}}$. It is written as

$$\begin{align} f_\mathrm{e}\left(x_{\mathrm{sh}},v\right) & = \overline{n}_\mathrm{e} \sqrt{\dfrac{m_\mathrm{e}}{2\pi T_\mathrm{e}^{\,{\mathrm{sh}}}}} \, e^{-\dfrac{m_\mathrm{e} v^2}{2 T_\mathrm{e}^{\,{\mathrm{sh}}}}} \operatorname{H}\left(v +v_\mathrm{c}\right),\nonumber\\ & = \overline{f_\mathrm{e}} \operatorname{H}\left(v +v_\mathrm{c}\right), \end{align} \tag{ F2 }$$

where $\operatorname{H}(v +v_\mathrm{c})$ is the Heaviside function that has a value of one when $v\unicode{x2A7E} -v_\mathrm{c}$, and zero otherwise. For the sake of compactness of the following computation we introduced the $\overline{f_\mathrm{e}}$ notation that stands for the non-truncated Maxwellian

$$\begin{align*}\overline{f_\mathrm{e}} = \overline{n}_\mathrm{e} \sqrt{\dfrac{m_\mathrm{e}}{2\pi T_\mathrm{e}^{\,{\mathrm{sh}}}}} \, e^{-\dfrac{m_\mathrm{e} v^2}{2 T_\mathrm{e}^{\,{\mathrm{sh}}}}}.\end{align*}$$

Note that the quantity $\overline{n}_\mathrm{e}$ is not the electron density at the sheath entrance $n_\mathrm{e}^{\mathrm{sh}}$. The density at sheath entrance $n_\mathrm{e}^{\mathrm{sh}}$ can be found by integrating against velocity equation (F2). By doing so, we find

$$\begin{align} n_\mathrm{e}^\mathrm{sh} & = \int_{-v_\mathrm{c}}^{+\infty} \textrm{d}v\, \overline{f_\mathrm{e}} \nonumber \\ & = \int_{-\infty}^{+\infty} \textrm{d}v\, \overline{f_\mathrm{e}} - \int_{-v_\mathrm{c}}^{+\infty} \textrm{d}v\, \overline{f_\mathrm{e}} \nonumber \\ & = \overline{n}_\mathrm{e} - \int_{v_\mathrm{c}}^{+\infty} \textrm{d}v\, \overline{f_\mathrm{e}}. \end{align}$$

In the following we neglect the contribution of the density of the tail of the electron distribution in the overall electron density at the entrance of the sheath, i.e. we assume that $\overline{n}_e \approx n_e^\mathrm{sh}$. It leads to the well-known prediction of the sheath potential drop, see below equation (F4).The expression of the electron flux at the sheath entrance is therefore

$$\begin{align} \Gamma_\mathrm{e}^{\mathrm{sh}} & = \int_{-v_\mathrm{c}}^{+\infty} \textrm{d}v\, v n_e^{\mathrm{sh}}\sqrt{\dfrac{m_\mathrm{e}}{2\pi T^{\mathrm{sh}}_\mathrm{e}}} e^{- \dfrac{m_\mathrm{e}v^2}{2T^{\mathrm{sh}}_\mathrm{e}}} \nonumber \\ & = n_\mathrm{e}^\mathrm{sh} e^{\dfrac{e\Delta \phi_\mathrm{sh}}{T^{\mathrm{sh}}_\mathrm{e}}} \sqrt{\dfrac{T^\mathrm{sh}_\mathrm{e}}{2\pi m_\mathrm{e}}}. \end{align} \tag{ F3 }$$

By equating the ion flux equation (F1) and electron flux equation (F3) we find the classical expression of the potential drop in the sheath region

$$\begin{align} \Delta \phi_{\mathrm{sh}} = \dfrac{T_\mathrm{e}^{\,{\mathrm{sh}}}}{e} \operatorname{log} \left( \sqrt{2\pi \dfrac{m_\mathrm{e}}{m_\mathrm{i}} M^2\left( 1 + \dfrac{T_\mathrm{i}^{\mathrm{sh}}}{T_\mathrm{e}^{\,{\mathrm{sh}}}}\right)}\right). \end{align} \tag{ F4 }$$

As explained in section 4.1 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 \unicode{x2A7D} v_\mathrm{c}$, that is to say no fast electrons come back from the wall region. On figure 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

$$\begin{align*}\Gamma_\mathrm{i}^\mathrm{sh} = \int_{-v_\mathrm{c}}^{+\infty} \textrm{d}v\, v n_e^{\mathrm{sh}}\sqrt{\dfrac{m_\mathrm{e}}{2\pi T^{\,{\mathrm{sh}}}_\mathrm{e}}} e^{- \dfrac{m_\mathrm{e}v^2}{2T^{\,{\mathrm{sh}}}_\mathrm{e}}} - \Gamma_\mathrm{e}^{-},\end{align*}$$

where the integral 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 $\Gamma_\mathrm{e}^{-}$ accounts for the flux of fast electrons coming back from the wall at the sheath entrance, i.e.

$$\begin{align} \Gamma_\mathrm{e}^{-} = -\int_{-\infty}^{-v_\mathrm{c}} \textrm{d}v\, vf_\mathrm{e}^{\,{\mathrm{sh}}}. \end{align} \tag{ F5 }$$

In this expression $f_\mathrm{e}^\mathrm{sh}$ is the distribution function at the sheath entrance obtained from the simulation. Note that $\Gamma_\mathrm{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

$$\begin{align} \Delta \phi_{\mathrm{sh}}^{\mathrm{ref.}} = \dfrac{T_\mathrm{e}^{\,\mathrm{sh}}}{e} \operatorname{log} \left( \sqrt{2\pi \dfrac{m_\mathrm{e}}{m_\mathrm{i}} M^2\left( 1 + \dfrac{T_\mathrm{i}^{\,{\mathrm{sh}}}}{T_\mathrm{e}^{\,{\mathrm{sh}}}}\right) }\times \left( 1 + \dfrac{\Gamma_\mathrm{e}^{-}}{\Gamma_\mathrm{i}^{\mathrm{sh}}} \right) \right). \end{align} \tag{ F6 }$$