Numerical stability of plasma sheath

We are interested in developing a numerical method for capturing stationary sheaths, that a plasma forms in contact with a metallic wall. This work is based on a bi-species (ion/electron) Vlasov-Amp{\`e}re model proposed in [3]. The main question addressed in this work is to know if classical numerical schemes can preserve stationary solutions with boundary conditions, since these solutions are not a priori conserved at the discrete level. In the context of high-order semi-Lagrangian method, due to their large stencil, interpolation near the boundary of the domain requires also a specific treatment.


Introduction
When a plasma is in contact with a metallic wall, stationary boundary layers, called sheaths, form. They are due to both the reflexion/absorption properties of the wall and the charge imbalance coming from the mass difference between electrons and ions: electrons leave the plasma faster than ions. A drop of the self-consistent potential at the wall is indeed built up to accelerate ions and decelerate electrons just as to ensure an equal flux of ions and electrons (a zero current) at the wall. In this work, the starting point is a plasma sheath model proposed by Badsi et al. [3], for which a so-called kinetic Bohm criterion on the incoming ion distribution has been found out. Numerically, the sheath potential is computed by a non-linear Poisson solver.
We here consider the corresponding non-stationary ion/electron Vlasov-Ampère model (note that the stability of its linearized version has been studied in [2]) and its numerical implementation through a semi-Lagrangian scheme [6]. Semi-Lagrangian schemes are transport solvers with no stability conditions on the time step: time step is only constrained by the physical dynamics that has to be captured. High-order schemes, using high degree interpolation and time splitting, can be devised and parallel implementation can be considered. We raise here the following questions: • Question 1: How do semi-Lagrangian scheme preserve this equilibrium in time ?
• Question 2: How can we treat boundary conditions with large stencil interpolation ?
We will see that the equilibrium is stable provided that the incoming current, not exactly zero at the discrete level, is made small. This will be ensured by taking a fine grid in space and velocity and a time step small enough. For the second question, fictitious values of the distributions functions will be set at the spatial boundaries. Finally, we emphasize that the time step is constrained by the electronic dynamics and thus is very small here.

.1 Stationary solution
Following [3], we consider a one-dimensional domain [0, 1]: particles enter the domain on the left (x = 0) and the metallic wall is located at abscissa x = 1. We consider in this work purely absorbing boundary conditions at the wall.
Electrons and ions are described by their distribution functions in phase space, denoted f sh se (x, v) and f sh s i (x, v) with spatial variable x ∈ [0, 1] and velocity variable v ∈ R, satisfy the dimensionless stationary Vlasov system: v ∂ x f sh where µ 1 is the mass ratio between electrons and ions, and E sh (x) denotes the stationary electric field. It satisfies the Poisson equation: where φ sh (x) is the electric potential, ε 1 is the dimensionless Debye length of the plasma and ρ sh (x) is the charge density. Boundary conditions. At the entrance of the domain, incoming ions and electrons have distributions: where n 0 is an electronic density parameter, that will be determined to have a given charge density ρ 0 at the entrance of the domain. Thus, electrons have a Maxwellian distribution. At the wall, ions and electrons are supposed here to be totally absorbed: The boundary conditions of the potential (3) are defined as follows. The potential values at x = 0 is set to 0 and the potential at the wall, φ w = φ sh (1), also called floating potential, is fixed such that the current vanishes at the wall: Note that it implies that the current vanishes in the whole domain, since both electron and ion momentum are constant in space.
Solutions. Using the characteristic lines of the transport equation and under the hypothesis that the potential is decreasing in the domain, the solution to problem (1)-(2)-(3) is given by: where the sheath potential is the solution to the following non-linear Poisson equation: (10) complemented by the boundary conditions (6)-(7).
Let consider a given charge density ρ 0 ∈ R at the entrance of the domain. The problem is solved in two steps: 1. a negative potential drop at the wall φ w ≤ 0 is uniquely defined satisfying (7) provided that the ionic distribution satisfies the relation: It satisfies the following equation: and the associated electronic density parameter is given by: 2. then, once this potential wall exists, the non-linear Poisson equation (10), complemented with the boundary conditions φ sh (0) = 0 and φ sh (1) = φ w , has a unique solution under the kinetic Bohm criterion: Parameters. The equilibria has been computed by the gradient algorithm described in [3] with the following parameters: • the dimensionless Debye length is set to ε = 0.01.
• the charge is zero at the entrance of the domain: ρ 0 = 0.
• the incoming ion distribution is given by (12) Z is the macroscopic ion velocity when entering the domain.
• the incoming electronic distribution is Maxwellian as given by (5). For the sake of completeness, we recall the expression: µ is the mass ratio of a deuterium plasma.
Given these parameters, we can compute n 0 , φ w with (11)-(1), and the sheath electric potential φ sh by solving the non-linear Poisson equation (10). In the numerical results, we take: In practice, for future use in the non-stationary model, we store the nonlinear potential φ sh on a uniform grid of N = 2048 cells, that is, at points j/N, j = 0, . . . , N .

Non-stationary model
In the non-stationary setting, electron and ion distribution functions, f se (t, x, v) and f s i (t, x, v), and the electric field, E(t, x), depend also on time t > 0. They satisfy the non-linear Vlasov-Ampère system: where J denotes the current density: Initial data. The initial data are given by the stationary solution: We note that the initial electric field E(x, 0) satisfies the following Poisson equation: where the densities n s i and n se are given by: and φ w is the floating potential at x = 1.

Numerical scheme
We solve system (14)-(15)-(16) using a semi-Lagrangian scheme. Due to the large mass ratio between electrons and ions, the computational domain of electrons will be larger than the ion one and thus velocity meshes have to be different for the two kinds of particles. Mesh notations. We consider uniform cartesian meshes for both spatial domain [0, 1] and velocity domain [v min , v max ]. Considering N x +1 points in the spatial direction and N v +1 points in the velocity one, mesh points are denoted ( Note that the definitions relative to the velocity domain implicitly depend on the particle s ∈ {s e , s i } under consideration. We denote by f n s,(i,j) the approximate value of the distribution function f s at point (x i , v j ) at time t n = n∆t, with ∆t > 0 and n ∈ N, and E n i the approximate value of the electric field at point x i .

Initialization of the electric field
The solution to system (17)-(18) is given by: The charge densities n s , s ∈ {s e , s i } are computed at grid points x i using the trapezoidal formula in the velocity direction from the distributions functions 1 . We then consider a reconstruction at any point where L k are the elementary Lagrange polynomials: Note that the densities have to be defined at points outside of the domain when d > 0. This will be done by expanding the definition of the distribution functions. This point will be detailed in Section 3.3. The discrete electric field at grid points x i is then obtained from (19) in which n s are replaced by their discrete counterparts n s,h . The involved integrals are computed exactly and the overall scheme has O(∆x 2d+2 ) accuracy.

Splitting
To solve the Vlasov-Ampère system, we use a splitting between space and velocity dynamics and write it as a succession of one-dimensional advections. More precisely, we consider the following dynamics: the kinetic transport system given by: and the electric transport system, given by: Since the electric field is constant in time in this second step, both dynamics T and U are advections at constant velocities. In practice, to obtain second order accuracy in time, we consider the Strang splitting which consists in computing for s ∈ {s e , s i } f n+1 where T h,τ and U h,τ are discrete approximations of T and U over a time interval τ . Each transport equation is solved using a semi-Lagrangian scheme with centered Lagrange interpolation of degree 2d + 1. We thus consider the discretized version of the kinetic and electric transport dynamics. Operator T h,τ : ((f s,(i,j) ), consists in the following: for any grid points (x i , v j ), we define the shifted index i * and α ∈ [0, 1] such that x i − v j τ = x i * + α∆x and then the distribution function is given by: where the (L k ) are defined in (20), and the electric field at point x i is given by: where the current density (J i ) i (resp. (J * i ) i ) is computed by trapezoidal rule in velocity from the discrete distribution function (f i,j ) (resp. (f * i,j )). We thus exactly solve in time the transport equations (21)-(22) at the grid points, starting from the interpolated distribution function in the spatial direction, while the Ampère equation (23) is computed using a second order Crank-Nicolson scheme. The interpolation requires values of the distribution function outside the domain: we explain in the next section how to extrapolate it.
Operator U h,τ : ((f s,(i,j) ), (E i ) i ) → ((f * s,(i,j) ), (E * i ) i ) consists in the following: for any grid points (x i , v j ), we define the shifted indexes j * s i , j * se and α s For this advection in velocity U h,τ , periodic boundary conditions are used; we have here made the presentation for Lagrange interpolation, as it is used for advection in space, but other advection scheme can be used; in particular, in the numerical results, we will use cubic splines.

Boundary conditions
In both the initial computation of the electric field and the advection in space T h,τ , the proposed numerical scheme requires to take values of the distribution function outside the physical domain.
For any x i = i∆x < 0 and v j , we consider the following extrapolation at the entry x = 0: For any x i+Nx = (i + N x )∆x > 1 and v j , we consider the following extrapolation at the wall x = 1: This corresponds to a purely Dirichlet condition for incoming velocities and an extension by imparity for leaving velocities (also called butterfly procedure).

Sheath test-case
We here used the following set of parameters:  Figure 2), we represent the distribution of electrons (resp. ions) at time t = 0 (top) and t = 4 (bottom). We clearly see that the equilibrium is well preserved. We see that the maximum principle is not exactly preserved, as no limiting procedure is introduced both in space advection (d = 8, i.e. Lagrange interpolation of degree 17) and in velocity advection (cubic splines). However, this seems to be not crucial here, as it is not far from being preserved and it can be used as a measure of accuracy of the simulation as other theoretically preserved quantities. We see that the distribution function of electrons presents a discontinuity; this has the effect that the quadrature in velocity for computing the current converges slowly, and thus a high number of points in velocity is needed.
For ions, we see that the distribution is not constant in space near the wall (x = 1) for a given velocity. At this boundary, the butterfly procedure does not destroy the C 1 property of the distribution function and seems to be more adapted than the prolongation by a constant value, which is used for incoming velocities. Note also that the fictitious boundary values are only used to interpolate inside the domain, as the sub time step τ is here always positive (negative time steps, could be however considered when going to higher order splitting schemes in time, but are not studied in this work).
The space discretization is rather fine; this is needed for the sharp gradient near the wall. Non-uniform grids could be useful here to save memory and computations [4], but are not tackled here for simplicity. We mention also that high order interpolation is used, which permits not to have to refine too much. High order schemes in a non uniform setting could be considered with a Semi-Lagrangian Discontinuous Galerkin method (SLDG) (see [5], for example), or with non uniform cubic splines (see for example [1]).
Finally, we note that the time step is chosen very small. In order to capture the oscillations due to the dynamic of the electrons, it should be indeed smaller than 2 · 10 −4 . However, as we are interested in studying the stability of the equilibrium for long times (at least long with respect to the electron dynamics), we have to reduce the time step to one order of magnitude in order to reduce the time error (due to the Strang splitting and the Crank-Nicolson scheme for the Ampère equation). Here again, higher order schemes in time might be useful, in order to use larger time steps (see for example [1]).
On Figure 3 (resp. 4), we represent the difference between the distribution function of electrons (resp. ions) at time 0 and at time t = 4 and t = 8. We see that the error for the electrons is mainly localized near the wall and at the discontinuity of the distribution. As regards ions, the error is localized near the wall. We see that the error near the wall increases in time, especially for the ions. We remark also an exceptional value for the electrons, that seems to be located at v = 0 and x = 1: this might be explained by the fact that the function is not modified by the space advection (at v = 0), and is thus not diffused.
On Figure 5, we represent 2d views of the distribution function of electrons and ions at time t = 8. We recognize here the pictures of [3].
Current density at the entry. On Figure 6, we represent the current density at the entry x = 0 (which is zero at the continuous level). We expect a better behavior on the this boundary than at the wall, where the convergence is more delicate. We remark an oscillatory behavior due to the electron dynamics. After a violent behavior at time around 0.01, where the value can have a peak around 0.02 (that value is guessed to be linked to the space discretization), it tends to decrease and stabilize, while still oscillating around 0, with the same frequency.
Other diagnostics. On Figure 7, we represent the total energy, whose continuous expression is We see that it is not exactly conserved. Except at the very beginning, it is decreasing and reaches a relative error of about 1.4% at final time. We then represent the time evolution of total density and L 1 norm of the distribution function 2 for ions and electrons on Figure 8 and the L 2 norm of the distribution function on Figure 9. We see that their time behavior is similar to the total energy one. Note that the L 1 norm and total density are indistinguishable (for both ions and electrons), which shows that the positivity of the distribution functions is rather well preserved.

Comparison with other parameters
On Figure  . We compare the current density at the entry x = 0 and the total energy. We observe that the change in velocity domain has no influence, while the change in time step ∆t leads to less accurate results.
On Figure 11, we present the results when taking less discretization points N x = 256, N v = 512 and velocity domain [−500, 500] for electrons and [−10, 10] for ions and ∆t = 10 −4 . This enables to consider longer time simulation, but we see that the results are degraded especially in large time.
On Figure 12, we consider low order interpolation d = 0. We consider ∆t = 5 · 10 −6 and N v = 8192 or N v = 65536. The velocity domain are [−500, 500] for electrons and [−10, 10] for ions. We see that the current density at the entry x = 0 is improved when using more points 2 Total density is given by:         in velocity. This is however not true as regards the other diagnostics: for the total energy, the case d = 0 is really more diffusive and can not compete with the high order scheme.
On Figures 13 and 14, we consider the case d = 0, with ∆t = 10 −4 and N x = 256 or N x = 2048 and N v = 8192 or N v = 512. Comparing with the previous case with d = 8 ( Figure  12), we see that there is no such oscillatory behavior in large time. However, for N x = 256 and N v = 512, we see that diffusion is really important and leads to different behavior: the current density at entry tends to a value around 0.1 instead of 0.

Conclusion
We have studied the behavior of the numerical solution of the Vlasov equation, initialized with a sheath equilibrium [3]. Thanks to high resolution in velocity, high order interpolation and very small time steps, we are able to recover the equilibrium accurately for relatively long time. This work is a first step as regards the numerical method. We mention here several directions of research: • to reduce the constraint on the time step, asymptotic preserving schemes could be designed, the fictitious boundary values may be improved and we could consider unsplit time integration, • to reduce the constraint on space/velocity grid, adaptive/Discontinuous-Galerkin method or delta-f method could be considered, • scheme ensuring a discrete Gauss law could be develop and its impact on the numerical results could be analyzed, • we could enhance mixed openmp/mpi parallelization to get full performance on current and future architectures.