# High order numerical methods for Vlasov-Poisson models of plasma sheaths

Valentin Ayot, Mehdi Badsi, Yann Barsamian, Anaïs Crestetto, Nicolas Crouseilles, Michel Mehrenberger, Averil Prost, Christian Tayou-Fotso

► To cite this version:

Valentin Ayot, Mehdi Badsi, Yann Barsamian, Anaïs Crestetto, Nicolas Crouseilles, et al. High order numerical methods for Vlasov-Poisson models of plasma sheaths. 2023. [hal-03926305](https://hal-03926305)

HAL Id: hal-03926305

https://inria.hal.science/hal-03926305

Preprint submitted on 6 Jan 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
High order numerical methods for Vlasov-Poisson models of
plasma sheaths

Valentin Ayot,\* Mehdi Badsi,† Yann Barsamian,‡ Anaïs Crestetto,§
Nicolas Crouseilles,¶ Michel Mehrenberger,∫ Averil Prost,\*\* Christian Tayou-Fotso.††

Abstract

This article is a report of the CEMRACS 2022 project, called HIVLASHEA, standing for "High order methods for Vlasov-Poisson models for sheaths". A two-species Vlasov-Poisson model is described together with some numerical simulations, permitting to exhibit the formation of a plasma sheath. The numerical simulations are performed with two different methods: a first order classical finite difference (FD) scheme and a high order semi-Lagrangian (SL) scheme with Strang splitting; for the latter one, the implementation of (non-periodic) boundary conditions is discussed. The codes are first evaluated on a one-species case, where an analytical solution is known. For the two-species case, cross comparisons and the influence of the numerical parameters for the SL method are performed in order to have an idea of a reference numerical simulation.

Acknowledgements

Centre de Calcul Intensif d'Aix-Marseille is acknowledged for granting access to its high performance computing resources. This project has been supported by the French Federation for Magnetic Fusion Studies (FR-FCM) and by the French ANR project MUFFIN ANR-19-CE46-0004. This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No 101052200 EUROfusion). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them.

Finally, authors would like to thanks the CEMRACS organizers, Pierre Navaro for the computational support and the CIRM's boards.

Introduction

Plasmas are neutral at the equilibrium in a sufficiently large domain. However, near a boundary, a charge imbalance may be observed in a thin layer called sheath. This phenomenon stems from the interaction of ions and electrons with the boundary media (a cold metallic wall, for instance). Both species will be absorbed by the wall, but with a rate proportional to their speed. Since the electrons are moving faster than the ions, a positively charged layer (the Debye sheath) forms near the boundary.
Plasma sheaths are particularly challenging to simulate for several reasons. In this region, the plasma parameters (temperature or density) develop steep gradients. Due to their different mass, electrons and ions

<footnote>\*Institut de Mathématiques, CNRS, UMR 5251, Université de Bordeaux, France. valentin.ayot@u-bordeaux.fr</footnote>

<footnote>†Nantes Université, CNRS, Laboratoire de Mathématiques Jean Leray, LMJL, UMR 6629, F-44000 Nantes, France.
mehdi.badsi@univ-nantes.fr</footnote>

<footnote>‡École européenne de Bruxelles 1, Belgique, yann.barsamian@teacher.eursc.eu, ⊕</footnote>

<footnote>§Nantes Université, CNRS, Laboratoire de Mathématiques Jean Leray, LMJL, UMR 6629, F-44000 Nantes, France.
anais.crestetto@univ-nantes.fr</footnote>

<footnote>¶Inria (Mingus team), IRMAR UMR 6625 and ENS Rennes, Univ Rennes, France. nicolas.crouseilles@inria.fr</footnote>

<footnote>||I2M, CMI, UMR 7373, Aix-Marseille Université, France. michel.mehrenberger@univ-amu.fr</footnote>

<footnote>\*\*LMI (EA 3226 - FR CNRS 3335), INSA de Rouen, France. averil.prost@insa-rouen.fr</footnote>

<footnote>††Labo. J. A. Dieudonné, UMR 6621, Université Nice-Sophia Antipolis, France. christian.tayou-fotso@unice.fr</footnote>
have a very different behaviour close to the boundary, and a two-species model is required to describe the formation of the sheath. Moreover, kinetic models are necessary to capture the velocity effect of the particles (see [Rie91]), and we have to deal with different scales to account for each species. Hence, we are considering in this work a two-species Vlasov model to run simulations of the plasma sheath.

Recent works are dedicated to the non homogeneous equilibrium sheath, both on the mathematical side and on the numerical side [Bad21, BBC21, DPB16]. Regarding the dynamical approaches, one refers to [CM14, BMG+22] or [BMN18] for which this study is a follow-up. In the latter work, we studied the behaviour of the numerical solution of the Vlasov equation, initialized with a sheath homogeneous equilibrium. In this work, our purpose is to investigate numerically the formation of a sheath, when we include ionization in the model, inspired by the recent work [ALPM+20] where a fluid model is proposed.

Then, we are concerned by the numerical approximation of the two-species Vlasov model including ionization and boundary. To do so, we propose a high order semi-Lagrangian method combined with a time splitting method. The presence of boundaries requires specific adaptation to define correctly the semi-Lagrangian method close to the boundary. We adapt the strategies developed in [CL20, BNS+21] to define inflow and outflow ghost points according to the order of the Lagrange interpolation used in the semi-Lagrangian method.

The paper is organized as follows: in Section 1, we introduce the two-species Vlasov-Poisson system with boundary, the numerical methods are described in Section 2 and Section 3 is devoted to numerical results. The different codes are available respectively on https://github.com/JuliaVlasov/DynamicElectricSheath.jl for the Finite Difference solver, and on https://github.com/averil-prost/HiVlaSheaSL for the Semi-Lagrangian solver.

# 1 Plasma sheaths

**The model** Let $t \in \mathbb{R}^+$ denote the time variable, $x \in [-1, 1]$ denote the spatial variable in a normalized one-dimensional domain, and $v \in \mathbb{R}$ denote the speed variable. Each species is described through their density in the phase space, denoted by $f_i : (t, x, v) \in \mathbb{R}^+ \times [-1, 1] \times \mathbb{R} \mapsto \mathbb{R}$ for the ions, and $f_e : (t, x, v) \in \mathbb{R}^+ \times [-1, 1] \times \mathbb{R} \mapsto \mathbb{R}$ for the electrons. Macroscopic quantities will be needed in the modelling, so that we consider the densities $\rho_{i,e}$ and currents $J_{i,e}$ defined by

$$ \rho_{i,e}(t, x) := \int_{v \in \mathbb{R}} f_{i,e}(t, x, v)dv, \quad \text{and} \quad J_{i,e}(t, x) := \int_{v \in \mathbb{R}} v f_{i,e}(t, x, v)dv. \qquad (1.1) $$

In the sequel, we will denote $\rho(t, x) := \rho_i(t, x) - \rho_e(t, x)$, and $J(t, x) := J_i(t, x) - J_e(t, x)$. Phase space densities are assumed to obey the Vlasov equation where the force term is the self-consistent electrostatic potential $\varphi$ which obeys the Poisson equation. In short, the system writes

$$ \left\{ \begin{array}{ll} \partial_t f_i + v \partial_x f_i - \partial_x \varphi \partial_v f_i = \nu f_e, & (t, x, v) \in \mathbb{R}_*^+ \times (-1, 1) \times \mathbb{R}, & (1.2a) \\ \partial_t f_e + v \partial_x f_e + \frac{\partial_x \varphi}{\mu} \partial_v f_e = 0, & (t, x, v) \in \mathbb{R}_*^+ \times (-1, 1) \times \mathbb{R}, & (1.2b) \\ -\lambda^2 \partial_{xx}^2 \varphi = \rho, & (t, x) \in \mathbb{R}^+ \times (-1, 1). & (1.2c) \end{array} \right. $$

The physical parameters $\nu$, $\mu$ and $\lambda$ have the following meaning:

- $\nu \ge 0$ is the ionization frequency. It describes the rate of creation of ions in presence of electrons.

- $\mu := m_e/m_i$ is the mass ratio between electrons and ions.

- $\lambda > 0$ is the Debye length.

In the sequel, we may use the electric field $E(t, x) := -\partial_x \varphi(t, x)$ instead of the potential. Then, the Poisson equation rewrites as

$$ \lambda^2 \partial_x E(t, x) = \rho(t, x), \quad (t, x) \in \mathbb{R}^+ \times (-1, 1). \qquad (1.3) $$
**Remark 1.1.** To reduce the notations, we will use $f_s$, $s \in \{i, e\}$ to denote either the ionic or the electronic distribution. The Vlasov equations (1.2a) and (1.2b) rewrite

$$ \partial_t f_s + v \partial_x f_s - c_s \partial_x \varphi \partial_v f_s = S_s, $$

with the coefficients $c_s$ and source terms $S_s$ defined as

$$ c_i := 1, \quad c_e := -\frac{1}{\mu}, \quad S_i := \nu f_e, \quad S_e := 0. $$

The densities $f_i$ and $f_e$ are subject to initial and boundary conditions, given by

$$ \begin{cases} f_s(0, x, v) := f_s^0(x, v), & (x, v) \in (-1, 1) \times \mathbb{R}, \\ f_s(t, x = 1, v < 0) = f_s(t, x = -1, v > 0) := 0, & t \in \mathbb{R}_+. \end{cases} \quad (1.4a) \quad (1.4b) $$

The homogeneous boundary condition (1.4b) stems from the non-emitting wall model: the boundary absorbs particles without any reflection. This loss of ionic particles is compensated by the ionization source term in the right hand side of (1.2a). For the electrons, it is expected that the most energetic ones are absorbed by the wall, while the others are confined in the core.

To completely describe the model, we still need to provide boundary conditions for the Poisson problem (1.2c). A first one is given by the choice of a reference potential in the plasma core (at $x = 0$)

$$ \varphi(t, 0) = 0, \quad \forall t \in \mathbb{R}^+. \quad (1.5) $$

To derive a second boundary condition, we introduce a symmetry assumption.

**Symmetries** We will look for a symmetric potential

$$ \varphi(t, x) = \varphi(t, -x), \quad (t, x) \in \mathbb{R}^+ \times [-1, 1]. \quad (1.6) $$

By derivation with respect to $x \in (-1, 1)$, we immediately obtain

$$ \partial_x \varphi(t, x) = -\partial_x \varphi(t, -x), \quad \text{i.e.} \quad E(t, x) = -E(t, -x). $$

In fact, assuming (1.6), the other quantities of the system also exhibit nice symmetric properties. Let us notice that the Vlasov equations (1.2a) and (1.2b) are driven by the vector fields

$$ (t, x, v) \rightarrow (1, v, E(t, x)) := V_i(t, x, v) \quad \text{and} \quad (t, x, v) \rightarrow (1, v, -E(t, x)/\mu) := V_e(t, x, v). $$

Both these fields satisfy the radial symmetry $V_s(t, x, v) = V_s(t, -x, -v)$. In consequence, if we assume that $f_s^0(x, v) = f_s^0(-x, -v)$, the solutions $f_s(t, x, v)$ will be radially symmetric around $(t, 0, 0)$, i.e.

$$ f_i(t, x, v) = f_i(t, -x, -v) \quad \text{and} \quad f_e(t, x, v) = f_e(t, -x, -v), \quad \forall (t, x, v) \in \mathbb{R}^+ \times [-1, 1] \times \mathbb{R}. $$

In particular, we have for the densities and currents

$$ \rho_s(t, x) = \int_{v \in \mathbb{R}} f_s(t, x, v) dv = \int_{w \in \mathbb{R}} f_s(t, x, -w) dw = \int_{w \in \mathbb{R}} f_s(t, -x, w) dw = \rho_s(t, -x), $$

$$ J_s(t, x) = \int_{v \in \mathbb{R}} v f_s(t, x, v) dv = - \int_{w \in \mathbb{R}} w f_s(t, x, -w) dw = - \int_{w \in \mathbb{R}} w f_s(t, -x, w) dw = - J_s(t, -x). $$

**Remark 1.2** (Additional symmetry of $f_e$). Notice that the function $f: (t, x, v) \mapsto f_e(t, x, v) - f_e(t, x, -v)$ satisfies the linear equation

$$ 0 = \partial_t f(t, x, v) + v \partial_x f(t, x, v) - \frac{E(t, x)}{\mu} \partial_v f(t, x, v). $$

The boundary condition (1.4b) gives $f(t, \pm 1, \pm v < 0) = 0$. If, in addition, we assume that the initial condition $f_e^0$ satisfies $f_e^0(x, v) - f_e^0(-x, -v) = 0$, then we obtain

$$ f_e(t, x, v) = f_e(t, x, -v), \quad \forall (t, x, v) \in \mathbb{R}^+ \times [-1, 1] \times \mathbb{R}. \quad (1.7) $$
**Boundary conditions** If the potential $\varphi$ is regular, the symmetry assumption gives us $E(t, 0) = -E(t, 0)$, and we may use the centered Neumann boundary condition

$$ \partial_x \varphi(t, 0) = 0 \quad \text{or equivalently} \quad E(t, 0) = 0 \qquad (1.8) $$

to close the Poisson equation (1.2c). However, this relies on the assumption that $E$ is continuous at $x = 0$, and it is not clear whether any initial condition leads to a stationary state with such a regular electric field. We may avoid this constraint by deriving another Neumann condition (at $x = \pm 1$) on the potential which seems to us more physical.

First, we derive with respect to time the Poisson equation (1.2c)

$$ -\partial_t(\lambda^2 \partial_{xx}^2 \varphi) = \partial_t \rho, $$

and considering the difference between the $v$-integration of the Vlasov equations (1.2a) and (1.2b) gives (recalling $\rho = \rho_i - \rho_e$ and $J = J_i - J_e$)

$$ \partial_t \rho = \nu \rho_e - \partial_x J, $$

leads to (using $E = -\partial_x \varphi$)

$$ \partial_x(\lambda^2 \partial_t E + J) = \nu \rho_e, \quad \forall x \in (-1, 1). $$

Integrating now over $x \in [-1, 1]$ leads to

$$ \lambda^2 \partial_t E(t, 1) + J(t, 1) = \lambda^2 \partial_t E(t, -1) + J(t, -1) + \nu \int_{-1}^{1} \rho_e(t, x) dx, \qquad (1.9) $$

and using the symmetries $E(t, 1) = -E(t, -1)$ and $J(t, 1) = -J(t, -1)$, it comes

$$ \lambda^2 \partial_t E(t, \pm 1) + J(t, \pm 1) = \pm \frac{\nu}{2} \int_{-1}^{1} \rho_e(t, x) dx. \qquad (1.10) $$

Time integration gives a condition of the form $E(t, \pm 1) = C_{\pm}(t)$, where

$$ E(t, \pm 1) = C_{\pm}(t) := E(0, \pm 1) - \frac{1}{\lambda^2} \int_0^t J(s, \pm 1) ds \pm \frac{\nu}{2\lambda^2} \int_0^t \int_{-1}^1 \rho_e(s, x) dx ds. \qquad (1.11) $$

This more elaborate boundary condition could allow for discontinuities at $x = 0$. Notice that the values $E(t, \pm 1)$ now depends on the initial value $E(0, \pm 1)$, whereas the centered boundary condition (1.8) does not rely on the initial state.

# 2 Numerical methods

In this section, the numerical methods used to approximate the system (1.2) are described. First, we focus on the Poisson equation and then, two numerical methods are presented for the Vlasov part.

## 2.1 Poisson equation

The Poisson problem (1.2c) is solved with integral representations of the variable $E$. First, we consider the centered Neumann boundary condition (1.8). Then, integrating the Poisson equation (1.3) over $[0, x]$ yields

$$ E(t, x) = 0 + \frac{1}{\lambda^2} \int_0^x \rho(t, y) dy = \frac{1}{\lambda^2} \int_0^x \int_{v \in \mathbb{R}} [f_i(t, y, v) - f_e(t, y, v)] dv dy. \qquad (2.1) $$

**Remark 2.1.** Instead of imposing $E(t, 0) = 0$, we can impose $\int_{-1}^1 E(t, x) dx = 0$ (coming from the symmetry of $\varphi$ (1.6)) or even $E(t, 1) = -E(t, -1)$. Then we compute (2.1) up to a constant and adjust the constant a posteriori to verify the chosen constraint.
Let us now consider the boundary condition (1.11). The spatial domain $[-1, 1]$ is split into its positive and negative parts, and integrating (1.3) gives

$$E(t, x) = \begin{cases} E(t, 1) - \frac{1}{\lambda^2} \int_x^1 \rho(t, y) dy = C_+(t) - \frac{1}{\lambda^2} \int_x^1 \int_{v \in \mathbb{R}} [f_i(t, y, v) - f_e(t, y, v)] dv dy, & x \in [0, 1], \\ E(t, -1) + \frac{1}{\lambda^2} \int_{-1}^x \rho(t, y) dy = C_-(t) + \frac{1}{\lambda^2} \int_{-1}^x \int_{v \in \mathbb{R}} [f_i(t, y, v) - f_e(t, y, v)] dv dy, & x \in [-1, 0). \end{cases} \quad (2.2)$$

Moreover, $C_\pm$ is computed at the required time by integrating (1.10), between two consecutive time steps. For the integrals under consideration (in space, velocity and time), we will approximate them by the trapezoidal method.

## 2.2 Finite Differences (FD)

Define a numerical computational domain $\Omega := [-1, 1] \times [-\bar{V}, \bar{V}]$, with a large enough maximum speed $\bar{V}$. Let $(x_j, v_k)_{k \in [0, K]}^{j \in [0, J]}$ be a cartesian grid of $\Omega$ of step $(\Delta x, \Delta v)$. We discretize the advection equations on the subgrid $(x_j, v_k)_{k \in [1, K-1]}^{j \in [1, J-1]}$ by an explicit Euler scheme in time, and the upwind scheme in space:

$$\frac{f_{s,j,k}^{n+1} - f_{s,j,k}^{n}}{\Delta t} + D_{j,k}^{-} f_{s}^{n} \left( c_{s} E_{j}^{n} \right)_{+} + D_{j,k}^{+} f_{s}^{n} \left( c_{s} E_{j}^{n} \right)_{-} = S_{s,j,k}^{n}, \quad (2.3)$$

where $a_+ = \max(a, 0)$ and $a_- = \min(a, 0)$ are respectively the coordwise positive and negative parts, and the uncentered finite differences are defined as

$$D_{j,k}^{\pm} f := \pm \left( \frac{f_{j \pm 1, k} - f_{j,k}}{\Delta x}, \frac{f_{j,k \pm 1} - f_{j,k}}{\Delta v} \right).$$

The values of $f_{s,j,k}^n$ on the boundary ($j = 0, J$ and $k = 0, K$) are taken as follows:

- the boundary condition (1.4b) yields $f_{s,j,k}^n = 0$ whenever $x_j = -1, v_k > 0$ or $x_j = 1, v_k < 0$,

- it is considered that $\bar{V}$ is large enough to take the values on the speed boundary $v_k = \pm \bar{V}$ equal to 0,

- the remaining values $f_{s,j,k}^n$, $x_j = -1, -\bar{V} < v_k \le 0$ or $x_j = 1, 0 \le v_k < \bar{V}$ may be computed using the scheme (2.3), since the sign of the speed allows to use only inner points.

The upwind scheme is known to be diffusive, and stable under the CFL condition

$$1 - \max_k |v_k| \frac{\Delta t}{\Delta x} - |c_s| \max_j |E_j^n| \frac{\Delta t}{\Delta v} \ge 0, \quad \forall s \in \{i, e\} \text{ and } n \in [1, N].$$

Given $\Delta x$ and $\Delta v$, we deduce a sufficiently small value of $\Delta t$ with the bound

$$\Delta t \le \min \left( \frac{\Delta x}{\bar{V}}, \min(1, \mu) \frac{\Delta v}{E_{\max}} \right), \quad E_{\max} > 0 \text{ postulated a priori.} \quad (2.4)$$

## 2.3 Semi-Lagrangian (SL)

The full model (1.2) nicely lends itself to approximation by time splitting. Indeed, consider the following Strang splitting decomposition
$$ \frac{\Delta t}{2} : \begin{cases} \partial_t f_s + v \partial_x f_s = 0 & \text{Linear advection along } x, \\ \lambda^2 \partial_x E = \rho_i - \rho_e & \text{Poisson problem,} \end{cases} $$

$$ \frac{\Delta t}{2} : \partial_t f_i = \nu f_e \qquad \text{Ionization,} $$

$$ \Delta t : \partial_t f_s + c_s E \partial_v f_s = 0 \qquad \text{Linear advection along } v, $$

$$ \frac{\Delta t}{2} : \partial_t f_i = \nu f_e \qquad \text{Ionization,} $$

$$ \frac{\Delta t}{2} : \begin{cases} \lambda^2 \partial_x E = \rho_i - \rho_e & \text{Poisson problem,} \\ \partial_t f_s + v \partial_x f_s = 0 & \text{Linear advection along } x. \end{cases} $$

Each of the splitting steps may be solved exactly in time. Indeed, the Poisson problems are solved by the integral representation (2.1). The ionization steps are pointwise ODE with time-independent source term, and are exactly solved by the explicit Euler scheme. Finally, notice that each advection is at constant speed with respect to the advection variable. This allows for the use of elementary 1D solvers.

The use of the decentered boundary condition requires a discretization of (2.2), and may penalize the time order of the splitting method. In the sequel, we will consider only the symmetric boundary condition (1.8).

**Numerical treatment of the boundaries** Let us focus on the elementary advection equation with constant speed $a > 0$

$$ \partial_t f(t, x) + a \partial_x f(t, x) = 0, \quad f(t, -1) = 0, \quad \forall (t, x) \in \mathbb{R}_*^+ \times (-1, 1). $$

Let $(x_j)_{j \in [0, J]}$ be a space mesh of step $\Delta x := 2/J$, and $(t_n)_{n \in [0, N]}$ be a time mesh of step $\Delta t := T/N$. We follow the work of [CL20, BNS$^{+21}$], and consider a semi-Lagrangian scheme defined as

$$ f_j^{n+1} = \text{Lagrange interpolation} (f^n, x_j - a\Delta t) := \sum_{k=-d}^{d+1} f_{j_0+k}^n L_k(\alpha), $$

with $(L_k)_{k \in [-d, d+1]}$ the Lagrange polynomials of degree $(2d+1)$ defined by $L_k(z) = \prod_{\ell=-d, \ell \neq k}^{d+1} \frac{z-\ell}{\ell-k}$ (which satisfy $L_k(\ell) = \delta_{k\ell}$ for $\ell \in [-d, d+1]$), and $x_j - a\Delta t = x_{j_0} + \alpha \Delta x$, $j_0 \in \mathbb{Z}$, $\alpha \in [0, 1[$. The boundaries are treated as follows:

- the *inflow* side, corresponding to $x = -1$, relies on the analytical solution $f(t, x) = 0 \quad \forall x \le at$. Whenever the scheme needs a value $f_j^n$ with $j < 0$, it may be exactly taken equal to 0.

- in the case $d > 0$, the Lagrange stencil may also need *outflow* values $f_j^n$ with $j > J$. Such values may be determined by polynomial extrapolation. Let $k_b \in \mathbb{N}$, and let $p$ be the unique polynomial of degree $k_b$ interpolating $(x_j, f_j^n)$ for $j \in [J - k_b, J]$. The *outflow ghost points* will be defined by $f_j^n := p(-1 + j\Delta x)$, $\forall j > J$. The couple $(d, k_b)$ characterizes the chosen scheme.

## 3 Numerical results

The semi-Lagrangian code is written in C; it uses subroutines (translated from Fortran) of the `selalib` library<fnref n="1" />. The code works in parallel using MPI. The finite difference scheme has been written in Julia.

<footnote n="1">¹https://selalib.github.io/selalib.html</footnote>
<figure><img src="image_1.png" /><figcaption>Figure 1: Malkov solutions (3.4) on [-1.5, 1.5] x [-2, 2]. Left: electric field E. Right: density f.</figcaption></figure>

The electric field is extended outside of [-1, 1] by a constant. The density f is represented in the domain [-1.5, 1.5] x [-2, 2], and truncated to 10.

## 3.1 One species validation test case

We rely on the work of [MK20] to provide an analytical solution in a one-species case. Consider the simplified model describing the density of particles $f = f(t, x, v)$, and the potential $\varphi = \varphi(t, x)$:

$$ egin{cases} 
\partial_t f + v \partial_x f - \partial_x \varphi \partial_v f = 0, & (t, x, v) \in \mathbb{R}_*^+ \times (-1, 1) \times \mathbb{R}, \\ 
\partial_{xx}^2 \varphi = \int_{v \in \mathbb{R}} f dv, & (t, x) \in \mathbb{R}^+ \times (-1, 1). 
\end{cases} 
\quad (3.1a) 
$$

The initial and boundary conditions are given by

$$ egin{cases} 
f(0, x, v) := f^0(x, v), & f(t, x = \pm 1, \pm v < 0) = 0, \\ 
\varphi(t, 0) = \partial_x \varphi(t, 0) = 0. 
\end{cases} 
\quad (3.2a) 
$$

This model may be seen as a particular case of the two-species Vlasov-Poisson (1.2), upon taking the following parameters:

$$ f_i^0 = 0, \quad \nu = 0, \quad \mu = -1, \quad \lambda = 1, \quad f_e^0 = f^0. $$

One may verify that (3.1) is solved in $\mathbb{R}^+ \times [-1, 1] \times \mathbb{R}$ by the following stationary couple:

$$ f(t, x, v) := \begin{cases} \frac{1}{\pi} (1 - x^2 - v^2)^{-1/2} & \text{if } x^2 + v^2 < 1, \\ 0 & \text{otherwise,} \end{cases} \quad \text{and} \quad \varphi(t, x) := \frac{x^2}{2}. \quad (3.3) $$

It is numerically relevant to extend the Malkov solution (3.3) to spatial domains $x \in [-1 - \varepsilon, 1 + \varepsilon]$ by

$$ f(t, x, v) := \begin{cases} \frac{1}{\pi} (1 - x^2 - v^2)^{-1/2} & \text{if } x^2 + v^2 < 1, \\ 0 & \text{otherwise} \end{cases}, \quad \text{and} \quad \varphi(t, x) := \begin{cases} x^2/2, & |x| < 1, \\ |x| - \frac{1}{2}, & |x| \ge 1. \end{cases} \quad (3.4) $$

Figure 1 illustrates the stationary solutions.

### Validation of the FD scheme
This stationary test first let us validate the Finite Difference scheme presented in Subsection 2.2 associated to the boundary condition (1.8) for the electric field (computed by (2.1)). Considering (x,v)-domain [-1, 1] x [-2, 2], we compute the L²-errors on electric field at time 0.1.
<figure><img src="image_2.png" /><figcaption>Figure 2: FD code for one species case. L²-errors on electric field at time 0.1, as a function of Δx (left) or Δv (right).</figcaption></figure>

Figure 2 presents the error as a function of Δx for different fixed values of $N_v = K + 1$ (on the left) and as a function of Δv for different fixed values of $N_x = J + 1$ (on the right). Time step Δt is computed thanks to the CFL condition (2.4), K (resp. J) being chosen such that Δt is fixed along each curve. In both plots (presented in log-log scale), we recover convergence of rate 1, as expected, until a saturation is reached. In Figure 3 we aim at comparing the three approaches for the computation of electric field presented in Subsection 2.1. Imposing $E(t, x = 0) = 0$ as proposed in (1.8) or $\int E(t, x)dx = 0$ as discussed in Remark 2.1 leads to two similar solutions, that are exact at $x = 0$ due to the properties of this case. Now, computing E from (1.11)-(2.2) leads to an exact solution on boundaries $x = -1$, $x = 1$ but a jump in $x = 0$. Even if the exact solution is continuous in $x = 0$, it is not always the case in physically relevant applications.

**Validation of the SL scheme** Results for the semi-Lagrangian code are given on Figure 4. Instead of evaluating the numerical solution on a point $(x, v)$, we use a quadrature to approximate its average on the cell of size $(\Delta x, \Delta v)$ centered in $(x, v)$. We evaluate the error on the electric field for different values of $N_x = N_v = N$. We see on the left figure that the scheme converges when $N$ increases at initial time. Then, on the right figure, we plot the error on the electric field at $T = 1$. Dividing the time step by two (resp. by four), we observe that the error is divided by around 4 (resp. 16), which illustrates the second order in time of the Strang splitting. The test is quite difficult (the function goes to infinity on one side and then is equal to zero on the other side) and large values of $N$ have to be used to clearly see the convergence.

## 3.2 Two species case

In this part, we focus on the two-species model (1.2) and use the following physical parameters:

$$ \lambda = \frac{1}{2}, \quad \mu = \frac{1}{100}, \quad \text{and} \quad \nu = 20. \qquad (3.5) $$

The initial conditions are chosen as the thermodynamic equilibrium in an infinite spatial domain, or in a domain with periodic condition. The densities are then given by

$$ f_i^0(x, v) := \frac{\exp\left(-\frac{v^2}{2}\right)}{\sqrt{2\pi}} \quad \text{and} \quad f_e^0(x, v) := \sqrt{\mu} \frac{\exp\left(-\mu\frac{v^2}{2}\right)}{\sqrt{2\pi}}. $$

In order to satisfy the boundary conditions, we multiply $f_{i,e}$ by a mask, defined as

$$ \text{mask}(x, v) := \frac{1}{2} \left( \tanh \left( \frac{x - (-0.8)}{0.1} \right) - \tanh \left( \frac{x - 0.8}{0.1} \right) \right). $$
<figure><img src="image_3.png" /><figcaption>Figure 3: FD code for one species case. Parameters $N_x = 2048$, $N_v = 2049$, time 0.1. Computed minus exact $E$ as a function of $x$ with three approaches (imposing $E(t, x=0) = 0$, imposing $\int E(t,x)dx = 0$, computing $E(t,x)$ from $J$ with (2.2)).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 4: SL code for one species case. Left: Electric field error for $x \in [-1, 1]$, at initial time. Right: at time $T = 1$, error of electric field for $\Delta t = 0.1$, 4 times error of electric field for $\Delta t = 0.05$, 16 times error of electric field for $\Delta t = 0.05$. We use $N \in \{16, 32, 64, 128, 256, 512, 1024\}$ cells in each direction. For each cell we compute the average value of $f$ at initial time using $N$ points in each direction.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 5: Initial conditions $f_i^0$ (left) and $f_e^0$ (right).</figcaption></figure>

If the boundary condition is not satisfied, the sharp profile at the boundary may induce spurious oscillations in the Lagrange interpolation during the first steps of the simulation. Figure 5 illustrates the resulting initial conditions.

The simulations run over the spatial domain $x \in [-1, 1]$. The semi-Lagrangian (SL) code computes the electron velocities on $v_e \in [-60, 60]$, and ion velocities on $v_i \in [-50, 50]$. The finite differences (FD) code uses the same mesh for ions and electrons, chosen as $v \in [-60, 60]$. To simplify the comparison, visualisations of $f_i$ are restricted to the coordinates $f_{i,j,k}$ such that $v_k \in [-5, 5]$. For the finite differences code, we use $N_x = 512$, $N_{v_e} = N_{v_i} = 513$; the time step is computed in order to satisfy the CFL condition; a further time step is used for terminating to the final time $T$. For the semi-Lagrangian code, we will use

<table><thead><tr><th>Run</th><th>N<sub>x</sub></th><th>N<sub>v<sub>i</sub></sub></th><th>N<sub>v<sub>e</sub></sub></th><th>Δt</th></tr></thead><tbody><tr><td>Run0</td><td>512</td><td>513</td><td>513</td><td>0.000250</td></tr><tr><td>Run1</td><td>256</td><td>2049</td><td>8193</td><td>0.000250</td></tr><tr><td>Run2</td><td>1024</td><td>2049</td><td>8193</td><td>0.000250</td></tr><tr><td>Run3</td><td>1024</td><td>2049</td><td>8793</td><td>0.000025</td></tr></tbody></table>

For the SL scheme, we always use $d = 8$ with periodic boundary conditions for the interpolation in velocity, and $k_b = 1$ together with $d = 2$ for the spatial interpolation. Lagrange interpolation of degree 3 (that is, $d = 1$) is used for passing from the ion velocity mesh to the electron velocity mesh (which is needed for the ionization step).

As diagnostics, we represent the electric field $E$, the charge density $\rho = \rho_i - \rho_e$, the ion distribution function $f_i$ and the electron distribution function $f_e$ at different times $T \in \{0.1, 0.2, 1, 2, 5, 20\}$.

**Short time** We first look at the results for short times: $T = 0.1$ (Figure 6) and $T = 0.2$ (Figure 7). For the semi-Lagrangian code, we use the parameters of Run2, which gives a reference solution. We see that the results are very similar with the FD code for $T = 0.1$, which permits to validate the results by cross comparisons. For $T = 0.2$, differences begin to appear; the sharp profile of $f_e$ is not well reproduced by the FD scheme which has also a coarser mesh. We can see also the differences on the charge density.

**Long time** Now, we turn to the long-term simulations, with $T = 20$, on Figure 8. The finite difference code suffers from numerical diffusion, and the approximations of $f_i$ and $f_e$ are very damped. Electric field and charge density are really different. On the contrary, the SL scheme gives much better results (in comparison with refined runs, as we will see later), while it uses the same number of points. We can distinguish a little degradation of the symmetry, by looking at the charge density.
<figure><img src="image_6.png" /><figcaption>Figure 6: Comparison between finite differences (left) and semi-Lagrangian (right) schemes at $T = 0.1$. The semi-Lagrangian code uses parameters of Run2.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 7: Comparison between finite differences (left) and semi-Lagrangian (right) schemes at $T = 0.2$. The semi-Lagrangian code uses parameters of Run2.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 8: Comparison between finite differences (left) and semi-Lagrangian (middle and right) schemes at $T = 20$. The semi-Lagrangian code uses parameters of Run0 (middle $E(t, 0) = 0$, right $\int_{-1}^{1} E(t, x)dx = 0$.)</figcaption></figure>
**Looking for a reference solution** It is not easy to find a reference solution, as we look for long time solution, and the ionization parameter needs to be adjusted to ensure the stationarity of the solution. The space scale must be able to capture thin structures, as the Vlasov equation is prone to filamentation. We have here also the problem that electrons and ions have different time scales, which implies that short time steps have to be used in order to be able to follow the dynamics (even if we use here a semi-Lagrangian scheme, whose time step is not restricted by the strong CFL condition of the finite difference scheme).
It is interesting to notice that we can have converged solutions in space for a fixed time step, taking for example $\Delta t = 0.0025$ (or even $\Delta t = 0.025$), but the corresponding solution is then not (at all) converged in time (results are shown on Figure 9; we can remark that the results on charge density can be quite different; looking at the electron distribution function, we see some shifting of the solution, and even more changes for $\Delta t = 0.025$).

So we have diminished the time step to $\Delta t = 0.00025$, and the comparison with the same simulations with $\Delta t = 0.000025$ on Figure 10 is now much similar, which indicates that the time step is now fine enough. Note that for semi-Lagrangian schemes, taking the smallest possible time step does not necessarily lead to a better result (due to accumulation of errors, as the number of interpolations increases).

We choose to use a very fine discretization in velocity for the electrons by taking $N_{v_e} = 8193$ (some results, not shown, with $N_{v_e} = 16385$ have lead to indistinguishable results). We then make vary $N_x$ and $N_{v_i}$; we remark that $N_x$ can be quite lower and $N_{v_i}$ a little lower also, as we have similar results taking $(N_x, N_{v_i}) = (256, 2049)$ and $(N_x, N_{v_i}) = (1024, 4097)$, as we can see on Figure 11 (the first one is more diffusive which is coherent). We remark also that the spatial interpolation does not lead to unstable results, which can occur sometimes when extrapolation is used (see [BMN18]). To our understanding, this is due to the fact that interpolation of nontrivial functions happens only on the outflow boundary, and any oscillation is evacuated immediately without propagation in the computational domain. We have used an odd number of points in order to prevent from having 0 as a mesh point, which would lead to increase of the value at constant rate, in the ionization step.

The study of the numerical equilibrium and the study of the behaviour of the ion density around (0,0), which is quite complex at equilibrium, are not tackled here and are out of the scope of this work.

In Figure 12, we compare the solutions given by the semi-Lagrangian code when three different boundary conditions for the Poisson solver are used. At the numerical level, we see that the first one $E(t, 0) = 0$ yields a loss of symmetry in the solution while the two others preserve it. At a physical level, we observe that the electric field is monotone increasing and consequently low energetic electrons (see Figure (d)) are essentially confined while ions are prone to be accelerated and lost at the wall (see Figure (c)).

# 4 Conclusion

This work discussed the development of a semi-Lagrangian scheme for the numerical simulation of a two-species plasma model involving sheath. Specific treatment of the boundary conditions have been proposed. For comparison, a 1st-order finite difference scheme has been implemented. Both schemes are first validated on a one-species test case for which an analytic solution is known. Then we obtained simulations for the full two-species case.

We obtain first promising results that encourage us to go further in the study of this semi-Lagrangian scheme coupled with boundary conditions. Long-time behaviour is a difficult question for which it would be very interesting to have a reference solution, given by the stationary problem. This is however out of the scope of this work.

# References

[ALPM$^{+}$20] A. Alvarez Laguna, T. Pichard, T. Magin, P. Chabert, A. Bourdon, and M. Massot. An asymptotic preserving well-balanced scheme for the isothermal fluid equations in low-temperature plasmas at low-pressure. *Journal of Computational Physics*, 419:109634, October 2020.
<figure><img src="image_9.png" /><figcaption>(a) $\rho$ for $\Delta t = 0.025$. $N_x = 256$, $N_{v_e} = N_{v_i} = 1023$ (left): $N_x = 4096$, $N_{v_e} = 8193$, $N_{v_i} = 16385$ (right)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(b) $f_e$ for $\Delta t = 0.025$. $N_x = 256$, $N_{v_e} = N_{v_i} = 1023$ (left): $N_x = 4096$, $N_{v_e} = 8193$, $N_{v_i} = 16385$ (right)</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(c) $\rho$ for $\Delta t = 0.0025$. $N_x = 256$, $N_{v_e} = N_{v_i} = 1023$ (left): $N_x = 512$, $N_{v_e} = N_{v_i} = 8193$ (right)</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(d) $\rho$ for $\Delta t = 0.0025$. $N_x = 256$, $N_{v_e} = N_{v_i} = 1023$ (left): $N_x = 512$, $N_{v_e} = N_{v_i} = 8193$ (right)</figcaption></figure>

<caption>Figure 9: Density $\rho$ and electron distribution function $f_e$ for time $T = 5$; with $\Delta t = 0.025$ and $\Delta t = 0.0025$
The semi-Lagrangian code is used; coarse mesh on the left and fine mesh on the right.</caption>
<figure><img src="image_13.png" /><figcaption>Figure 10: Density ρ (left) and electron distribution function (middle and right) for time $T \in \{1, 2, 5, 10\}$; The semi-Lagrangian code is used with parameters of Run2 (left, middle) and of Run3 (right).</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 11: Comparison of semi-Lagrangian code: parameters of Run1 (left column) and parameters of Run2 (right column).</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 12: Comparison of the numerical approximation of the Poisson equation in the semi-Lagrangian code: $E(t, 0) = 0$ (left column), $E(t, 1) = -E(t, -1)$ (middle column) and $\int_{-1}^{1} E(t, x)dx = 0$ (right column). $\Delta t = 0.025$. $N_x = 256$, $N_{ve} = N_{vi} = 1023$.</figcaption></figure>
[Bad21] Mehdi Badsi. Collisional sheath solutions of a bi-species Vlasov-Poisson-Boltzmann boundary value problem. *Kinetic & Related Models*, 14(1):149, 2021.

[BBC21] Mehdi Badsi, Christophe Berthon, and Anaïs Crestetto. A stable fixed point method for the numerical simulation of a kinetic collisional sheath. *Journal of Computational Physics*, 429:109990, March 2021.

[BMG+22] Emily Bourne, Yann Munschy, Virginie Grandgirard, Michel Mehrenberger, and Philippe Ghen-drih. Non-Uniform Splines for Semi-Lagrangian Kinetic Simulations of the Plasma Sheath. August 2022.

[BMN18] Mehdi Badsi, Michel Mehrenberger, and Laurent Navoret. Numerical stability of plasma sheath. *ESAIM: Proceedings and Surveys, SMAI 2017*, 64:17-36, 2018.

[BNS+21] B. Boutin, T.H.T. Nguyen, A. Sylla, S. Tran-Tien, and J.-F. Coulombel. High order numerical schemes for transport equations on bounded domains. *ESAIM: Proceedings and Surveys*, 70:84-106, 2021.

[CL20] Jean-François Coulombel and Frédéric Lagoutière. The Neumann numerical boundary condition for transport equations. *Kinetic & Related Models*, 13(1):1-32, 2020.

[CM14] David Coulette and Giovanni Manfredi. An Eulerian Vlasov code for plasma-wall interactions. *Journal of Physics: Conference Series*, 561(1), November 2014.

[DPB16] Bruno Després, Martin Campos Pinto, and Mehdi Badsi. A minimization formulation of a bi-kinetic sheath. *Kinetic and Related Models*, 9(4):621-656, September 2016.

[MK20] Evgeny A. Malkov and Alexey N. Kudryavtsev. Non-stationary Antonov self-gravitating layer: Analytics and numerics. *Monthly Notices of the Royal Astronomical Society*, 491:3952-3966, January 2020.

[Rie91] K.-U. Riemann. The Bohm criterion and sheath formation. *Journal of Physics D: Applied Physics*, 24(4):493, April 1991.