<figure><img src="image_1.png" /></figure>

HAL
open science

# A variational sheath model for gyrokinetic vlasov-poisson equations

Mehdi Badsi, Martin Campos Pinto, Bruno Després, Ludovic Godard-Cadillac

► To cite this version:

Mehdi Badsi, Martin Campos Pinto, Bruno Després, Ludovic Godard-Cadillac. A variational sheath model for gyrokinetic vlasov-poisson equations. ESAIM: Mathematical Modelling and Numerical Analysis, 2021, 55 (6), pp.2609-2642. 10.1051/m2an/2021067. hal-03160211v2

HAL Id: hal-03160211

https://hal.sorbonne-universite.fr/hal-03160211v2

Submitted on 1 Feb 2022

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.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A VARIATIONAL SHEATH MODEL FOR STATIONARY GYROKINETIC VLASOV–POISSON EQUATIONS

MEHDI BADSI<fnref n="1" />*, MARTIN CAMPOS-PINTO<fnref n="2" /> ID,
BRUNO DESPRÉS<fnref n="3" /> AND LUDOVIC GODARD-CADILLAC<fnref n="4" /> ID

**Abstract.** We construct a stationary gyrokinetic variational model for sheaths close to the metallic wall of a magnetized plasma, following a physical extremalization principle for the natural energy. By considering a reduced set of parameters we show that our model has a unique minimal solution, and that the resulting electric potential has an infinite number of oscillations as it propagates towards the core of the plasma. We prove this result for the non linear problem and also provide a simpler analysis for a linearized problem, based on the construction of exact solutions. Some numerical illustrations show the well-posedness of the model after numerical discretization. They also exhibit the oscillating behavior.

**Mathematics Subject Classification.** 78A30, 49J05, 34C10.

Received March 11, 2021. Accepted October 5, 2021.

## 1. INTRODUCTION

The mathematical description and numerical simulation of a plasma interacting with a metallic surface is an active topic of research in plasma physics [1, 10, 15, 21]. Applications are wide spread over a large spectrum, from the simulation of solar wind to design of laboratory plasma devices [20] in particular for the ITER project https://www.iter.org/proj/inafewlines which is to produce electric energy from thermonuclear fusion. One of the main feature of isolated plasmas interacting with a metallic surface is the development near the surface of a thin positively charged layer of several Debye length in thickness. This layer is called a sheath and results from the relative mobility difference between ions and electrons. In the sheath, a significant electric field aims at repelling the electrons in the core plasma while ions are accelerated towards the surface so that an equal flux of ions and electrons leaving the plasma to the wall is reached. The mathematical description of this transition layer between the core plasma and the wall is fundamental for the understanding of the plasma properties in its globality. Significant efforts have been made in this direction at the physical level [14, 16, 33, 45]. At the

*Keywords and phrases.* Plasma-wall interaction, Vlasov–Poisson equation, gyrokinetic model, gyroaverage operator, Bohm criterion, floating potential, stationary solutions.

<footnote n="1">¹ Université de Nantes, CNRS UMR 6629,2 rue de la Houssinière, 44332 Nantes, France</footnote>

<footnote n="2">² Max Planck Institute for Plasma Physics, Boltzmannstr. 2, 85748 Garching, Germany</footnote>

<footnote n="3">³ Laboratoire Jacques Louis Lions, Sorbonne Université, 4 place Jussieu, 75005 Paris, France</footnote>

<footnote n="4">⁴ Università degli Studi di Torino, dipartimento di Matematica, Giuseppe Peano. Via Carlo Alberto 10, Torino, TO 10123, Italia</footnote>

<footnote>*Corresponding author: mehdi.badsi@univ-nantes.fr</footnote>
mathematical level, the kinetic description of plasmas has become standard and the mathematical theory of kinetic plasma models be it in the absence of boundaries or in the case of bounded plasma is by now well-established [1, 6, 7, 9, 26–29, 38, 39, 41]. However efficient, little works on the mathematical side have focused on the precise analysis of sheaths [3–5, 18]. In [4], the authors showed the existence and uniqueness of a solution for a two species Vlasov-Poisson system with boundaries under a moment condition on the incoming flow of ions that takes the form of a so called kinetic Bohm criterion [8, 40]. The extension of this work to the case where a constant magnetic field is imposed was intended in [2]. At the physical level, when a magnetic field is considered, numerical simulations [16, 33, 35] seem to show off different physical scenarios according to both the angle of incidence of the magnetic field with respect to the surface and the ratio of force between the electric field and the intensity of the magnetic field. The existence of physical stationary states seems nevertheless subjected to the validity of the so called Bohm-Chodura condition [16, 45] whose mathematical justification is up to our knowledge a condition for the existence of a monotonic electrostatic potential obtained via a linearized model.

As far as the mathematical modeling of magnetized plasma is concerned, gyrokinetic theory has emerged to reduce the computational cost of computing the particles distribution functions in perpendicular directions to the magnetic field. In the classical kinetic description, particles are usually described with their six coordinates of positions and velocities. In the gyrokinetic description however, particles are identified, through their gyro-center positions in $\mathbb{R}^3$ and two other components: their parallel component along the magnetic field and their magnetic moment. The magnetic moment enables an identification with the perpendicular velocity. It yields a representation of the particles distribution function in a five dimensional phase space. The mathematical and physical justifications can be found in [10, 13, 19, 22, 34, 43]. In practice [11], the transport of particles distribution functions is made through the gyroaveraged fields which are non local in space. This non locality poses the question of what are appropriate boundary conditions when material boundaries are considered.

By far, the construction of solutions in the presence of spatial boundaries for a two species gyrokinetic Vlasov-Poisson system is up to our knowledge an open problem. In this work, we study a magnetized plasma interacting with a metallic wall by considering a simplified two species gyrokinetic Vlasov-Poisson system. As far as gyrokinetic Vlasov-Poisson equations are concerned, it is known at the physical level that transient solutions of such equations are prone to small scales instabilities which are the cause of plasma turbulence [21, 24]. These instabilities are interesting in their own and generally occur on time scale much smaller than that of fusion for which relaxation towards a stationary state is expected. Our main concern in this work, is the construction of reference stationary solutions that are compatible with the presence of a material boundary. It is an important topic [32, 42] but rigorous mathematical contributions are very sparse. Using symmetries, particles distribution functions in our model live in a three dimensional phase space where particles gyrocenter positions are spotted by the variable $x < 0$, their velocity component along the magnetic field is denoted $v \in \mathbb{R}$ and their magnetic moments is denoted $\mu \in \mathbb{R}$. The angle of the prescribed magnetic field is arbitrary. Incoming flow of particles is considered at the infinity and partial absorption is considered at the wall. In our model, the infinity is in fact the core of a plasma assumed to be neutral and at rest. The electrostatic potential is assumed to vanish at the infinity, while at the wall its gyroaveraged value is computed so as to ensure the neutrality of the current. These natural equations are supplemented by what we call the closure relations. These closure relations are the simplest ones, have a clear physical meaning even if they are arbitrary, and yield a well-posed global problem, as we show in this work.

To construct the model, our methodology follows a standard approach: we integrate the gyrokinetic Vlasov equations with respect to the gyroaveraged electrostatic potential, and we then consider a non linear Poisson equation that we treat as an extremalization problem. The density functions are peaked around an electronic Larmor radius $r_e$ and an ionic Larmor radius $r_i$. These radii are arbitrary in the final model.

For the mathematical analysis, it is convenient to choose equal Larmor radius $r_e = r_i$ in order to minimize technical difficulties. This hypothesis, the same unique radius for ions and electrons, is similar of other mono-kinetic hypotheses in the mathematical literature [12]. This simplification has the advantage that well posedness (convexity, existence and uniqueness of a minimum, etc.) is proved. A new result is a rigorous justification of the oscillating behavior of the electric potential, a phenomenon that has been observed numerically [16] but, to
the best of our knowledge, never analyzed mathematically. Specifically, we are able to prove that under certain conditions the exact electric potential admits an infinite number of oscillations of comparable length. This is first shown with a method adapted to the non linearity of the problem, and also with a linear method which calculates exact solutions to the linear problem far from the wall. Some illustrations are eventually provided in the numerical section. These illustrations fully confirm our original theoretical findings.

The organization is as follows. The variational model is constructed in Section 2, by means of exact integration of the kinetic Vlasov equations in the phase space. The value of the wall potential is obtained by solving a pair of non linear equations which express the neutrality at infinity and the ambipolarity principle. The variational principle is then formulated with a stability condition at infinity, corresponding to a kinetic Bohm–Chodura inequality. In Section 3, the model is simplified so as to establish its main properties. The mathematical solution is based on the minimization of a strictly convex functional. The monotonicity of the solution at the wall is proved in Section 3.1.2. The oscillating behavior of the solution is then proved by two different methods. The first one is described in Section 3.1.3. It is non linear and proceeds by means of inequalities and comparison with convenient test functions. The second one in Section 3.1.4 is linear in nature and is based on the calculation of exact complex exponential solutions. Finally Section 4 is dedicated to simple numerical tests which illustrate the theoretical findings. Several technical proofs and results are gathered in the Appendix.

## 2. CONSTRUCTION OF A VARIATIONAL SHEATH MODEL

### 2.1. Gyroaverage operator

In magnetized plasma physics, charged particles follow helicoidal trajectories. A key notion is the gyroaveraging on the Larmor radius. Let $\mathbf{B} \neq 0$ be a constant non zero given magnetic field

$$ \mathbf{B} = |\mathbf{B}| (\cos \theta \mathbf{e}_x + \sin \theta \mathbf{e}_y) \in \mathbb{R}^3 \quad \text{with} \quad \theta \in (0, \frac{\pi}{2}). \qquad (2.1) $$

Following [11,13], the dynamics of charged particles is described by averaging along helicoidal trajectories which oscillate around guiding-center trajectories. A typical helicoidal trajectory starting at $X = (x_1, x_2, x_3)$ involves the Larmor vector $\rho_r(\alpha) := r(\cos \alpha \mathbf{v}_1 + \sin \alpha \mathbf{v}_2) \in \mathbb{R}^3$ with $\alpha \in [0, 2\pi)$, and where $\mathbf{v}_1 = (-\sin \theta, \cos \theta, 0)$ and $\mathbf{v}_2 = (0, 0, 1)$ are two vectors orthonormal to $\mathbf{B}$. With these notations, a trajectory with radius $r > 0$, parallel direction $\mathbf{b} = \frac{\mathbf{B}}{|\mathbf{B}|}$, parallel velocity $v_\parallel \neq 0$ and cyclotron frequency $\omega_c \neq 0$ verifies the Newton's second law of motion

$$ X'(t) = V(t), \qquad (2.2) $$

$$ V'(t) = \omega_c V(t) \wedge \mathbf{b} \qquad (2.3) $$

whose general solution is

$$ X(t) = X + \rho_r(\omega_c t) + v_\parallel \mathbf{b}t := \tilde{X}(t) + v_\parallel \mathbf{b}t. \qquad (2.4) $$

Given a characteristic length $L$, gyrokinetic models are based on the approximation $u(X) \leftarrow \frac{\omega_c}{2\pi} \int_0^{2\pi/\omega_c} u(\tilde{X}(t))dt$ in the regime $\omega_c \gg v_\parallel/L$. Assume now that $u$ is an univariate function of the coordinate $x = x_1$. By averaging on the Larmor circle, one gets the approximation $u(x) \leftarrow \langle u \rangle_r(x) = \frac{1}{2\pi} \int_0^{2\pi} u(x + \rho_r(\alpha) \cdot \mathbf{e}_x)d\alpha$. Using the change of variable $\tilde{x} = \rho_r(\alpha) \cdot \mathbf{e}_x = -r \sin \theta \cos \alpha$, the local gyroaverage operator takes the form of a convolution operator. Namely for all $u \in C_c^\infty(\mathbb{R})$,

$$ \langle u \rangle_r(x) := \int_{\mathbb{R}} w \left( \frac{x-y}{r \sin \theta} \right) u(y) \frac{dy}{r \sin \theta}, \qquad (2.5) $$

where the convolution kernel $w$ is defined almost everywhere by

$$ w(x) = \frac{1}{\pi\sqrt{1-x^2}} \quad \text{for } |x| < 1 \quad \text{and} \quad w(x) = 0 \quad \text{for } |x| > 1. \qquad (2.6) $$
The convolution kernel $w$ belongs to $L^p(\mathbb{R})$ for $1 \le p < 2$ and its total mass is normalized, $\int_{\mathbb{R}} w(x)dx = 1$. Some properties of the associated gyroaverage operator are stated in the Appendix. In particular it can be extended to a linear continuous operator from $H^s(\mathbb{R})$ to $H^{s+\frac{1}{2}}(\mathbb{R})$ for any $s \in \mathbb{R}$.

## 2.2. Gyrokinetic equations

Following [11, 13, 20], a plasma made of one species of ions and electrons subject to a given magnetic field $\mathbf{B}$ is modeled with gyrokinetic Vlasov-Poisson equations. Positions are denoted as $x \in \mathbb{R}$. Velocities in the $x$ direction are denoted as $v \in \mathbb{R}$. Magnetic moments are denoted as $\mu = \pm|\mathbf{v}_{\perp}| \in \mathbb{R}$.

The unknowns are the electrostatic potential $\phi: x \in \mathbb{R} \mapsto \phi(x) \in \mathbb{R}$, the particle density for ions $f_i: (x, v, \mu) \in \mathbb{R} \times \mathbb{R} \times \mathbb{R} \mapsto f_i(x, v, \mu) \in \mathbb{R}^+$ and the particle density for electrons $f_e: (x, v, \mu) \in \mathbb{R} \times \mathbb{R} \times \mathbb{R} \mapsto f_e(x, v, \mu) \in \mathbb{R}^+$. They satisfy the gyrokinetic Vlasov-Poisson system (2.7)-(2.9)

$$v \partial_x f_i - \langle \partial_x \phi \rangle_{R_i(\mu)} \partial_v f_i = 0, \quad (x, v, \mu) \in \mathbb{R}^3, \qquad (2.7)$$

$$v \partial_x f_e + \frac{1}{\eta} \langle \partial_x \phi \rangle_{R_e(\mu)} \partial_v f_e = 0, \quad (x, v, \mu) \in \mathbb{R}^3, \qquad (2.8)$$

$$-\lambda^2 \sin^2 \theta \partial_{xx} \phi = n_i(x) - n_e(x), \quad x \in \mathbb{R}. \qquad (2.9)$$

In the equations, the notations are as follows. The Larmor radius is proportional to the square root of the magnetic moment for both species, that is $R_{i,e}(\mu) = \alpha_{i,e} \sqrt{|\mu|} > 0$ where $\alpha_{i,e} > 0$. We refer to p. 5 of Reference [11] for the definition of the Debye length $\lambda > 0$ in function of physical quantities in the context of fusion plasmas.

In particular small Debye length arise from the physical scaling

$$\frac{\rho_s^2}{\lambda^2} = \frac{4\pi nmc^2}{|\mathbf{B}|^2} \gg 1. \qquad (2.10)$$

One has $\lambda \approx 6.4e-5$ in fusion plasmas ([20], p. 126). The mass ratio between ions and electrons is a natural small parameter denoted as $\eta = \frac{m_e}{m_i} \approx 1/2000$. The coefficient $\sin^2 \theta$ in equation (2.9) comes from the fact that it is the projection along the first coordinate of the operator $\nabla_{\perp}^* \nabla_{\perp}$ where $\nabla_{\perp}$ is the gradient orthogonal to the magnetic field $\mathbf{B}$. Since one also has a $\sin \theta$ in the gyroaveraging (2.5), it is immediate to observe that $\sin \theta$ can be conveniently eliminated after a rescaling of the space (it will be performed in the next section).

The macroscopic gyroaveraged densities and parallel current densities are

$$n_i(x) = \int_{\mathbb{R}^2} \langle f_i \rangle_{R_i(\mu)}(x, v, \mu) \mathrm{d}v \mathrm{d}\mu, \quad n_e(x) = \int_{\mathbb{R}^2} \langle f_e \rangle_{R_e(\mu)}(x, v, \mu) \mathrm{d}v \mathrm{d}\mu, \qquad (2.11)$$

$$J_i = \int_{\mathbb{R}^2} \langle f_i \rangle_{R_i(\mu)}(x, v, \mu) v \mathrm{d}v \mathrm{d}\mu, \quad J_e = \int_{\mathbb{R}^2} \langle f_e \rangle_{R_e(\mu)}(x, v, \mu) v \mathrm{d}v \mathrm{d}\mu.$$

By integration of the Vlasov equations (2.7) and (2.8), we find that the parallel current densities $J_i$ and $J_e$ are constant.

Using the order of magnitude of physical quantities in fusion plasma ([20], p. 146) the mean Larmor radius for ions is $r_i \approx 5e-3m$ and the mean Larmor radius for electrons is $r_e \approx 8.3e-5m$. The ratio comes from $r_e/r_i = \sqrt{\eta} = \sqrt{m_e/m_i}$. Therefore, for the ease of mathematical simplifications, we will consider particles density which are peaked at their corresponding (positive) magnetic moment, $\mu_i$ such that $R_i(\mu_i) = r_i$ and $\mu_e$ such that $R_e(\mu_e) = r_e$. That is we will consider

$$f_i(x, v, \mu) \leftarrow f_i(x, v) \otimes \delta(\mu - \mu_i) \quad \text{and} \quad f_e(x, v, \mu) \leftarrow f_e(x, v) \otimes \delta(\mu - \mu_e).$$

## 2.3. Dimensionless equations

There are three characteristic lengths in the model which are $\lambda$, $r_i$ and $r_e$. It is convenient to perform a rescaling of the space variable

$$\hat{x} = \frac{x}{\lambda \sin \theta}, \quad \hat{r}_i = \frac{r_i}{\lambda}, \quad \hat{r}_e = \frac{r_e}{\lambda}.$$
The convolution kernels given by (2.6) are rescaled as well, that is $\hat{w}_i(\hat{x}) = \frac{1}{\hat{r}_i} w(\frac{\hat{x}}{\hat{r}_i})$ and $\hat{w}_e(\hat{x}) = \frac{1}{\hat{r}_e} w(\frac{\hat{x}}{\hat{r}_e})$. The rescaled density functions are $\hat{f}_i(\hat{x}, v) = f_i(x, v)$ and $\hat{f}_e(\hat{x}, v) = f_e(x, v)$. The rescaled electric potential is $\hat{\phi}(\hat{x}) = \phi(x)$. One gets a dimensionless Vlasov equations

$$v \partial_{\hat{x}} \hat{f}_i(\hat{x}, v) - (\hat{w}_i \star \hat{\phi})'(\hat{x}) \partial_v \hat{f}_i(\hat{x}, v) = 0 \quad \text{and} \quad v \partial_{\hat{x}} \hat{f}_e(\hat{x}, v) + \frac{1}{\eta} (\hat{w}_e \star \hat{\phi})'(\hat{x}) \partial_v \hat{f}_e(\hat{x}, v) = 0 \quad (2.12)$$

and a dimensionless Poisson equation

$$-\partial_{\hat{x}\hat{x}} \hat{\phi}(\hat{x}) = \hat{n}_i(\hat{x}) - \hat{n}_e(\hat{x}). \quad (2.13)$$

Here the dimensionless densities correspond to the gyroaveraged ones $\hat{n}_i(x) = n_i(x)$ and $\hat{n}_e(x) = n_e(x)$, namely

$$\hat{n}_i(x) = \int (\hat{w}_i \star \hat{f}_i)(\hat{x}, v) dv = \hat{w}_i \star \hat{n}_i^{\text{phys}} \quad \text{and} \quad \hat{n}_e(x) = \int (\hat{w}_e \star \hat{f}_e)(\hat{x}, v) dv = \hat{w}_e \star \hat{n}_e^{\text{phys}}, \quad (2.14)$$

where we have denoted by $\hat{n}_i^{\text{phys}} := \int \hat{f}_i(\cdot, v)dv$ and $\hat{n}_e^{\text{phys}} := \int \hat{f}_e(\cdot, v)dv$ the “physical” densities in dimensionless variables. In the rest of this work, the non dimensional notation $\cdot$ will be discarded.

## 2.4. Modeling sheath and boundary conditions

A sheath is the boundary layer at a metallic wall observed in real devices. We write (2.12) and (2.13) in the domain

$$(x, v) \in \mathcal{D} = (-\infty, 0) \times \mathbb{R}$$

and we must complement the equations with boundary conditions. However there is still a major difficulty near the wall, which is that some fundamental assumptions behind the gyroaveraging procedure are not satisfied near the wall at $x = 0$: indeed the validity of helicoidal trajectories (2.4) in the vicinity of the wall is quite questionable. This problem is well known in physical literature [8, 16, 33, 40, 44, 45].

In our approach, we will firstly use natural boundary conditions and we will secondly close the model with the energy extremalization principle. The natural boundary conditions are generalization of a previous work [4]. The energy extremalization principle can be seen as a way to recover some compatibility with the Hamiltonian description of Vlasov-Maxwell equations [36, 46].

The zero reflection law for particles at the wall writes

$$f_i(0, v < 0) = 0 \quad \text{and} \quad f_e(0, v < 0) = 0. \quad (2.15)$$

We will extend the functions inside the wall with

$$\left\{ \begin{array}{ll} \partial_x \phi = 0 & \text{in } (0, +\infty), \\ f_i = f_e = 0 & \text{in } (0, +\infty) \times \mathbb{R}^2. \end{array} \right. \quad (2.16)$$

These closure relations express that the spatial domain $(0, +\infty)$ corresponds to the inside of a perfectly conducting wall. Consequently there is neither an electric field inside, nor plasma particles that travel inside. These closure relations are arbitrary in a sense, however they are simple and have a clear physical meaning.

## 2.5. Electronic phase diagram

Starting from a given smooth potential

$$\psi = w_e \star \phi, \quad \psi \in C^1(-\infty, 0], \quad (2.17)$$

which may have oscillations as described in Figure 1, the electronic phase diagram constructs a physically and mathematically admissible solution of the gyroaveraged kinetic equation $v \partial_x f_e + \frac{1}{\eta} \psi'(x) \partial_v f_e = 0$ in $\mathcal{D}$.
<figure><img src="image_3.png" /><figcaption>FIGURE 1. Graph of potentials $\phi, \psi = w_e \star \phi$ and $\varphi = w_i \star \phi$ with oscillations. A priori $\phi \neq \psi \neq \varphi$, but the general shape is similar.</figcaption></figure>

The oscillations of the gyroaveraged potential $\psi$ are caused by similar oscillations of the physical potential $\phi$. The level curves are defined by $\frac{1}{2}v^2 - \frac{1}{\eta}\psi(x) = \text{constant}$. Where the potential is regular enough ($\psi \in C_{\text{loc}}^1$ for example), the trajectories correspond the characteristics curves which are the solutions of the characteristics equations $x' = v$ and $v' = \frac{1}{\eta}\psi'(x)$. The characteristics curves correspond to connected components of the level curves. Therefore, if two points $(x, v) \in \mathcal{D}$ and $(y, w) \in \mathcal{D}$ are connected by a characteristic curve then

$$ \frac{1}{2}v^2 - \frac{1}{\eta}\psi(x) = \frac{1}{2}w^2 - \frac{1}{\eta}\psi(y). $$

For $y = 0$, one gets the condition

$$ \frac{1}{2}v^2 - \frac{1}{\eta}\psi(x) \geq -\frac{1}{\eta}\psi(0). \qquad (2.18) $$

Making the assumption that the minimum of $\psi$ is reached at the wall, that is

$$ \psi(0) = \min_{\mathbb{R}^{-}} \psi, $$

then the inequality (2.18) delimits 3 open regions in $\mathcal{D}$, as described in Figure 2

$$ \begin{aligned} \mathcal{D}_1 &= \{ (x, v) \in \mathcal{D}, \quad v > \sqrt{\frac{2}{\eta}(\psi(x) - \psi(0))} \}, \\ \mathcal{D}_2 &= \{ (x, v) \in \mathcal{D}, \quad -\sqrt{\frac{2}{\eta}(\psi(x) - \psi(0))} < v < \sqrt{\frac{2}{\eta}(\psi(x) - \psi(0))} \}, \\ \mathcal{D}_3 &= \{ (x, v) \in \mathcal{D}, \quad v < -\sqrt{\frac{2}{\eta}(\psi(x) - \psi(0))} \}. \end{aligned} $$

The superior region $\mathcal{D}_1$ is connected to $x = -\infty$, so particles coming from the infinity with a positive velocity travel through $\mathcal{D}_1$ and are absorbed at the wall. The internal region $\mathcal{D}_2$ contains closed loops because the function $\psi$ may have some oscillations (we will prove this fact). Particles in this region do not reach the wall. The inferior region $\mathcal{D}_3$ is connected to the wall at $x = 0$. However no electron with negative velocity is emitted at the wall (2.15), so this region is empty of electrons.

**Remark 2.1.** The assumption $\psi(0) = \min_{\mathbb{R}^{-}} \psi$ originates from the pure electrostatic setting (that is $b = 0$) for which the physics of the sheath is dominant [4]. The case $\psi(0) \neq \min_{\mathbb{R}^{-}} \psi$ would yield a different phase diagram from Figure 1. Especially, one would obtain the existence of trapping sets, that is, sets which contain closed trajectories that do not intersect the boundary. On these sets, the density function may take arbitrary values, which makes the analysis more technical if they are non zero. Although such sets are possible in theory, our numerical experiments in Section 4 seem to sustain the fact that our simplifying assumption remains valid when a magnetic field is considered.
<figure><img src="image_4.png" /><figcaption>FIGURE 2. Three different regions of the electron phase space. The dotted lines are general characteristics. The bold line is the characteristic $\frac{1}{2}v^2 - \frac{1}{\eta}\psi(x) = -\frac{1}{\eta}\psi(0)$.</figcaption></figure>

We consider

$$ 
\begin{aligned}
f_e(x, v) &= n_{\text{ref}} \exp\left(-\frac{\eta}{2}v^2 + \psi(x)\right) & \text{for } (x, v) \in \mathcal{D}_1 \cup \mathcal{D}_2, \\
f_e(x, v) &= 0 & \text{for } (x, v) \in \mathcal{D}_3,
\end{aligned} 
\qquad (2.19) $$

where $n_{\text{ref}}$ is a reference density that will be determined later in Section 2.7. This representation is compatible with standard assumptions in plasma physics [15].

**Lemma 2.2.** The function $f_e$ defined in (2.19) is a weak solution in $\mathcal{D}$ to the equation $v\partial_x f_e + \frac{1}{\eta}\psi'(x)\partial_v f_e = 0$.

The “physical” electronic density $n_e^{\text{phys}} = \int f_e dv$, see (2.14), can be written in different ways

$$ 
\begin{aligned}
n_e^{\text{phys}}(x) &= n_{\text{ref}} \int \sqrt{\frac{2}{\eta}(\psi(x) - \psi(0))} \exp\left(-\frac{\eta}{2}v^2 + \psi(x)\right) dv \\
&= n_{\text{ref}} \exp\psi(x) \sqrt{\frac{2}{\eta}} \int \sqrt{\psi(x) - \psi(0)} e^{-t^2} dt.
\end{aligned} 
$$

Another possibility that will be used later is

$$ n_e^{\text{phys}}(x) = n_{\text{ref}} \exp \psi(x) \sqrt{\frac{\pi}{2\eta}} \left( 1 + \text{erf} \left( \sqrt{\psi(x) - \psi(0)} \right) \right), \qquad (2.20) $$

where the Gauss error function is $\text{erf}(x) = \frac{2}{\sqrt{\pi}} \int_{0}^{x} \exp(-t^2) dt$. The density at infinity $n_e^{\infty} = \lim_{-\infty} \int f_e(x, v) dv$ is

$$ n_e^{\infty} = n_{\text{ref}} \sqrt{\frac{\pi}{2\eta}} \left( 1 + \text{erf} \left( \sqrt{-\psi(0)} \right) \right). \qquad (2.21) $$

The constant electronic current is equal to its limit $J_e^{\infty} = \lim_{-\infty} \int f_e(x, v) v dv$, that is

$$ J_e^{\infty} = \frac{n_{\text{ref}}}{\eta} \exp \psi(0). \qquad (2.22) $$
<figure><img src="image_5.png" /><figcaption>FIGURE 3. Four different regions of the ion phase space. The dotted lines are general characteristics. The bold line is the characteristic $\frac{1}{2}v^2 + \psi(x) = \psi^+.</figcaption></figure>

## 2.6. Ionic phase diagram

Next we construct the phase diagram for ions for the potential

$$ \varphi = w_i \star \phi, \quad \varphi \in C^1(-\infty, 0]. $$

The function $\varphi$ may also have oscillations as in Figure 1. With the method of characteristics $x' = v$ and $v' = -\psi'(x)$, if two points $(x, v) \in \mathcal{D}$ and $(y, w) \in \mathcal{D}$ are connected with a characteristic curve then $\frac{1}{2}v^2 + \varphi(x) = \frac{1}{2}w^2 + \varphi(y)$. The maximal value

$$ \varphi^+ = \sup_{x \in \mathbb{R}^-} \varphi(x) = \varphi(x_+) $$

is a barrier of the potential. For simplicity as illustrated in Figure 1, we assume that the maximum value $\varphi^+ > 0$ is positive (consistently with the condition at infinity $\lim_{-\infty} \varphi = 0$) and that $x_+$ is unique. Indeed a point $(x, v)$ such that $\frac{1}{2}v^2 + \varphi(x) < \varphi^+$ cannot be connected to $(x_+, w)$. This is the reason of the decomposition of the domain in four open regions

$$ \begin{aligned} \mathcal{D}_4 &= \{ (x, v) \in \mathcal{D}, \quad v > \sqrt{2(\varphi^+ - \varphi(x))} \}, \\ \mathcal{D}_5 &= \{ (x, v) \in \mathcal{D}, \quad x_+ < x < 0, \quad -\sqrt{2(\varphi^+ - \varphi(x))} < v < \sqrt{2(\varphi^+ - \varphi(x))} \}, \\ \mathcal{D}_6 &= \{ (x, v) \in \mathcal{D}, \quad v < -\sqrt{2(\varphi^+ - \varphi(x))} \}, \\ \mathcal{D}_7 &= \{ (x, v) \in \mathcal{D}, \quad x < x_+, \quad -\sqrt{2(\varphi^+ - \varphi(x))} < v < \sqrt{2(\varphi^+ - \varphi(x))} \}. \end{aligned} $$

**Remark 2.3.** Note that the assumption that $x_+$ is unique simplifies the phase diagram 3. The case of several points where the maximum value $\varphi^+ > 0$ is reached would yield the existence of trapping sets, that is, set which contain closed trajectories that do not intersect the boundary. On these sets, the density may take arbitrary values, which makes the analysis more technical if they are non zero.

To construct a weak solution of $v\partial_x f_i - \varphi'(x)\partial_v f_i = 0$ in $\mathcal{D}$ based on the diagram of Figure 3, we consider a function at infinity for incoming ions $f_i^\infty(v)$. Using the approach [4], we assume that $f_i^\infty$ is continuous with compact support in $\{v > \sqrt{2\varphi^+ + \varepsilon}\}$. More precisely

$$ f_i^\infty \in C^0(\mathbb{R}), \quad f_i^\infty(v) = 0 \quad \text{for } v \in (-\infty, \sqrt{2\varphi^+} + \varepsilon] \cup [A, +\infty), $$
where $0 < \varepsilon$ is a small number and A is taken large enough so that $\sqrt{2\varphi^+} + \varepsilon < A$.
We consider the continuous function $f_i$

$$ 
\begin{aligned}
f_i(x, v) &= f_i^\infty \left( \sqrt{v^2 + 2\varphi(x)} \right) && \text{for } (x, v) \in \mathcal{D}_4, \\
f_i(x, v) &= 0 && \text{for } (x, v) \in \mathcal{D}_5 \cup \mathcal{D}_6 \cup \mathcal{D}_7.
\end{aligned} 
\qquad (2.23) $$

**Lemma 2.4.** The function $f_i$ defined in (2.23) is a weak solution in $\mathcal{D}$ to the equation $v\partial_x f_i - \varphi'(x)\partial_v f_i = 0$.

The “physical” ionic density $n_i^{\text{phys}} = \int f_i dv$, see (2.14), is

$$ 
\begin{aligned}
n_i^{\text{phys}}(x) &= \int_{\sqrt{2(\varphi^+ - \varphi(x))}}^{+\infty} f_i^\infty(\sqrt{v^2 + 2\varphi(x)}) dv \\
&= \int_{\sqrt{2\varphi^+}}^{+\infty} f_i^\infty(w) \frac{w}{\sqrt{w^2 - 2\varphi(x)}} dw \\
&= \int_0^{+\infty} f_i^\infty(w) \frac{w}{\sqrt{w^2 - 2\varphi(x)}} dw,
\end{aligned} 
\qquad (2.24) $$

where the last equality is obtained using the fact that $f_i^\infty$ vanishes for $v < \sqrt{2\varphi^+}$. The ionic density at infinity is

$$ n_i^\infty = \int_0^{+\infty} f_i^\infty(w) dw. $$

The electronic current is constant in space, so it is equal to its value at infinity

$$ J_i^\infty = \int_0^\infty f_i^\infty(w) w dw. $$

## 2.7. Determination of the sheath potential $\psi_{\text{wall}}$ and the reference density $n_{\text{ref}}$

In the physics of sheaths, the value of the potential at the wall is of critical importance. Here since the current densities are constant we may use the method from [4] which is to consider the equations obtained from $n_e^\infty = n_i^\infty$ and $J_e^\infty = J_i^\infty$. Considering equations (2.21) and (2.22), we denote

$$ \psi_{\text{wall}} = \psi(0). \qquad (2.25) $$

Observe that the condition of neutral charge at $-\infty$ and the condition of neutral flux take the form

$$ 
\left\{ 
\begin{array}{l}
n_{\text{ref}} \sqrt{\frac{\pi}{2\eta}} (1 + \text{erf}(\sqrt{-\psi_{\text{wall}}})) = n_i^\infty, \\
\frac{n_{\text{ref}}}{\eta} \exp(\psi_{\text{wall}}) = J_i^\infty.
\end{array}
\right. 
\qquad (2.26) $$

In this system, the right hand side is given because it depends on the function $f_i^\infty$ which is independent of $\psi$ and $\phi$. The unknown is the pair $(\psi_{\text{wall}}, n_{\text{ref}})$ with the sign conditions

$$ \psi_{\text{wall}} < 0 \quad \text{and} \quad n_{\text{ref}} > 0. \qquad (2.27) $$

**Proposition 2.5.** There exists a unique solution $(\psi_{\text{wall}}, n_{\text{ref}})$ to (2.26) and (2.27) if and only if $\sqrt{\pi\eta J_i^\infty} < \sqrt{2n_i^\infty}$.

*Proof.* From (2.26) one gets by division the reduced equation

$$ H(u) = \frac{\sqrt{2n_i^\infty}}{\sqrt{\pi\eta J_i^\infty}}, \quad H(u) = (1 + \text{erf}(\sqrt{u})) \exp(u), \quad u = -\psi_{\text{wall}} > 0. $$

The function $H$ is strictly monotone, with $H(0) = 1$ and $H(+\infty) = +\infty$. Therefore there exists a positive solution, necessarily unique, if and only if $1 < \frac{\sqrt{2n_i^\infty}}{\sqrt{\pi\eta J_i^\infty}}$. Once $\psi_{\text{wall}}$ is determined, the reference density is given. $\square$
## 2.8. A principle of energy extremalization

Energy minimization, or more generally energy extremalization, is a general principle in mathematical physics. It allows to identify those solutions which realize a minimum, or more generally a critical point, of some potential referred to as an energy. In the sequel, we show how to use this principle in order to complete the construction of our model.

Firstly we express the densities in function of the electric potential $\phi$, secondly we give the energy and thirdly we write the Euler-Lagrange equations for critical points. A stability condition is expressed in the form of a Bohm-Chodura condition.

### 2.8.1. Potential representation of the densities

From the general definition (2.13) and (2.19), one gets

$$n_e(\phi)(x) = w_e \star \left( n_{\text{ref}} \sqrt{\frac{\pi}{2\eta}} \exp(w_e \star \phi(\cdot)) \left( 1 + \text{erf}\left(\sqrt{w_e \star \phi(\cdot) - \psi_{\text{wall}}}\right) \right) \right) (x).$$

Let us define the function

$$N_e(\psi) = n_{\text{ref}} \sqrt{\frac{\pi}{2\eta}} \int_0^\psi \exp(z) \left(1 + \text{erf}\left(\sqrt{z-\psi_{\text{wall}}}\right)\right) dz, \quad \psi \ge \psi_{\text{wall}}. \qquad (2.28)$$

By definition one has that $n_e(\phi(x)) = w_e \star (N'_e(w_e \star \phi))(x)$. In (2.28), the function $N_e$ is defined as a double integral because the Gauss error function erf is itself an integral. Another definition is possible with just one integral.

**Lemma 2.6.** One has another formula

$$N_e(\psi) = 2n_{\text{ref}} \sqrt{\frac{\pi}{2\eta}} \left( \exp \psi - 1 - \frac{2}{\sqrt{\pi}} \int_{-\psi_{\text{wall}}}^{\infty} \exp(-w^2) w \left( \sqrt{w^2 + \psi} - w \right) dw \right), \quad \psi \ge \psi_{\text{wall}}. \qquad (2.29)$$

*Proof.* Differentiate the claim $N'_e(\psi) = n_{\text{ref}} \sqrt{\frac{\pi}{2\eta}} \left( 2 \exp \psi - \frac{2}{\sqrt{\pi}} \int_{-\psi_{\text{wall}}}^{\infty} \exp(-w^2) \frac{w}{\sqrt{w^2+\psi}} dw \right)$. Make a change of variable $z^2 = w^2 + \psi$ with $\sqrt{(\psi - \psi_{\text{wall}})} \le z < \infty$. One obtains

$$\begin{aligned}
N'_e(\psi) &= n_{\text{ref}} \sqrt{\frac{\pi}{2\eta}} \left( 2 \exp \psi - \frac{2}{\sqrt{\pi}} \int_{\sqrt{\psi - \psi_{\text{wall}}}}^{\infty} \exp(-z^2 + \psi) dz \right) \\
&= n_{\text{ref}} \sqrt{\frac{\pi}{2\eta}} \exp \psi \left( 2 - \frac{2}{\sqrt{\pi}} \int_{\sqrt{\psi - \psi_{\text{wall}}}}^{\infty} \exp(-z^2) dz \right) \\
&= n_{\text{ref}} \sqrt{\frac{\pi}{2\eta}} \exp \psi \left( 1 + \frac{2}{\sqrt{\pi}} \int_{0}^{\sqrt{\psi - \psi_{\text{wall}}}} \exp(-z^2) dz \right) \\
&= n_{\text{ref}} \sqrt{\frac{\pi}{2\eta}} \exp \psi \left( 1 + \text{erf}\left(\sqrt{\psi - \psi_{\text{wall}}}\right) \right)
\end{aligned}$$

By definition the function of the claim vanishes at the origin $N_e(0) = 0$. So it is equal to (2.28). $\square$

We use a similar approach for the ions. Let us consider

$$n_i(\phi)(x) = w_i \star \left( \int_{\sqrt{2(\phi^+ - w_i \star \phi(\cdot))}} \int_{\sqrt{v^2 + 2w_i \star \phi(\cdot)}}^{\infty} f_i^{\infty}(\sqrt{v^2 + 2w_i \star \phi(\cdot)}) dv \right) (x). \qquad (2.30)$$

Define the function

$$N_i(\varphi) = -\sqrt{2} \int_{\sqrt{2\varphi^+}}^{\infty} f_i^{\infty}(z) z \left( \sqrt{z^2 - 2\varphi} - z \right) dz, \quad \varphi \le \varphi^+ \qquad (2.31)$$

By definition $N_i(0) = 0$.
**Lemma 2.7.** One has the identity $n_i(\phi)(x) = w_i \star (N'_i(w_i \star \phi))(x)$.

*Proof.* One has $N'_i(\varphi) = \int_{\sqrt{2\varphi^+}}^{\infty} f_i^{\infty}(z) \frac{z}{\sqrt{z^2-2\varphi}} dz$. Make the change of variable $z = \sqrt{v^2+2\varphi}$ or equivalently $v = \sqrt{z^2-2\varphi}$ which satisfies $zdz = vdv$. So $N'_i(\varphi) = \int_{\sqrt{2(\varphi^+-\varphi)}}^{\infty} f_i^{\infty}(\sqrt{v^2+2(\varphi-\varphi^+)} ) dv$. Plugging this identity in (2.30) yields the claim. □

### 2.8.2. The energy

For mathematical correctness of the material presented below, it is necessary to have continuous extension of the function $N_e(\psi)$ for $\psi < \psi_{wall}$ and continuous extension of the function $N_i(\varphi)$ for $\varphi > \varphi^+$. It can be made on many ways, even if we expect that physical solutions might not depend on these continuous extensions, as it will be visible in the numerical simulations at the end of this work. So, from now on, we will assume that

$N_e, N_i \in C^0(\mathbb{R})$, $N_e, N_i$ are piecewise $C^1(\mathbb{R})$ and $N_e, N_i \in C^2(-\varepsilon, \varepsilon)$ for some $\varepsilon > 0$.

**Definition 2.8.** The energy is defined formally by

$$J(\phi) = \int_{-\infty}^{0} \left( \frac{1}{2} |\phi'(x)|^2 + N_e(w_e \star \phi(x)) - N_i(w_i \star \phi(x)) \right) dx. \quad (2.32)$$

An affine functional space naturally adapted to this energy is

$$V = \{ \phi \in \dot{H}^1(\mathbb{R}) \cap L^2(\mathbb{R}^-) \mid \phi' \in L^2(\mathbb{R}), \phi(x) = \phi(0) \text{ for } x > 0, w_e \star \phi(0) = \psi_{wall} \},$$

where we remind that $\dot{H}^1(\mathbb{R})$ is the set of measurable functions $\phi$ such that $\phi' \in L^2(\mathbb{R})$. The definition of $V$ accounts for the first condition of (2.16) inside the wall and for the boundary condition (2.25) at the wall. Functions $\phi \in V$ are such that $\lim_{+\infty} \phi(x) = 0$. We observe that $V$ is a complete space for the norm of $H^1(\mathbb{R}^-)$, $\|\phi\| = \sqrt{\|\phi\|_{L^2(\mathbb{R}^-)}^2 + \|\phi'\|_{L^2(\mathbb{R}^-)}^2}$. Since $\dot{H}^1(\mathbb{R})$ is continuously embedded in $L^\infty(\mathbb{R}^-)$ and functions in $V$ are constant on $\mathbb{R}^+$, then one has then

$$V \subset C^0(\mathbb{R}) \cap L^\infty(\mathbb{R}).$$

**Lemma 2.9.** The function $J$ is continuously defined from $\tilde{V} := V \cap L^1(\mathbb{R}^-)$ into $\mathbb{R}$.

**Remark 2.10.** The need for $L^1$ integrability is due to a mismatch between the $N_i$ and $N_e$ functionals at $-\infty$, it will be relaxed in the mathematical analysis of Section 3.

*Proof.* We need to prove that the integral (2.32) is convergent for functions in $\tilde{V}$ and focus on the only problem which is at $-\infty$. Since any $\phi \in V$ tends to zero at $-\infty$, we have $|\psi(x)|, |\varphi(x)| \to 0$ as $x \to -\infty$ for $\psi = w_e \star \phi$ and $\varphi = w_i \star \phi$. Now, by construction $N_e(0) - N_i(0) = 0$ and

$$N'_e(0) - N'_i(0) = n_e^{\infty} - n_i^{\infty} = 0, \quad (2.33)$$

so that one has a structural bound $|N_e(\psi) - N_i(\psi)| \le C\psi^2$ for $-x$ large enough. In particular, one has

$$|N_e(\psi) - N_i(\varphi)| \le |N_e(\psi) - N_i(\psi)| + |N_i(\psi) - N_i(\varphi)| \le C\psi^2 + D|\psi - \varphi|$$

for $-x$ large enough. As $\|\psi\|_{L^p}, \|\varphi\|_{L^p} \le \|\phi\|_{L^p}$ for $p \ge 1$ (see Appendix), this a priori estimate shows the convergence of the integral for $\phi \in \tilde{V}$. Refining this arguments shows that $J$ is continuous from $\tilde{V}$ into $\mathbb{R}$. □
### 2.8.3. A stability condition at infinity

Here we study a condition such that $J$ is convex with respect to small perturbation at $x \approx -\infty$ of the null function. The condition is written as a convex condition around $\phi = 0$

$$N''_e(0) - N''_i(0) > 0. \qquad (2.34)$$

**Lemma 2.11.** Assume $\int_0^{+\infty} f_i^\infty(w) \frac{1}{w^2} dw \le \int_0^{+\infty} f_i^\infty(w) dw$. Then the condition (2.34) of convexity at infinity holds.

**Remark 2.12.** The condition of this Lemma takes the form of a kinetic Bohm-Chodura condition [16,33,44,45].

*Proof.* Using the physical densities, one gets easily

$$ \begin{aligned} N''_e(0) - N''_i(0) &= \frac{d}{d\psi} \left( n_{\text{ref}} \sqrt{\frac{2}{\eta}} e^{\psi} \int_{-\sqrt{-\psi-\psi_{\text{wall}}}} e^{-t^2} dt \right) (0) - \frac{d}{d\varphi} \left( \int_{\sqrt{2\varphi^+}}^{+\infty} f_i^\infty(w) \frac{w}{\sqrt{w^2-2\varphi}} dw \right) (0) \\ &= n_{\text{ref}} \sqrt{\frac{2}{\eta}} \int_{-\sqrt{-\psi-\psi_{\text{wall}}}} e^{-t^2} dt + n_{\text{ref}} \sqrt{\frac{1}{2\eta}} e^{\phi_{\text{wall}}} - \int_{\sqrt{2\varphi^+}}^{+\infty} f_i^\infty(w) \frac{1}{w^2} dw. \end{aligned} $$

Using $\varphi^+ \ge 0$, one gets $N''_e(0) - N''_i(0) > n_i^\infty - \int_0^{+\infty} f_i^\infty(w) \frac{1}{w^2} dw \ge 0$ which ends the proof. $\square$

### 2.8.4. Extremalization

Even if the functional is convex at infinity under the condition of Lemma 2.11, a general guarantee of convexity is not known so far for arbitrary values of the physical coefficients ($r_i, r_e, f_i^\infty, \dots$). Note however that a general convexity result will be established in Section 3 but for $r_i = r_e$ which is usually not true in real plasmas. This is the reason we focus on the weaker notion of extremal solutions in this section.

Let $\tilde{V}_0 = V_0 \cap L^1(\mathbb{R}^-)$ be the function space tangent to $V$, with

$$V_0 = \{h \in \dot{H}^1(\mathbb{R}) \cap L^2(\mathbb{R}^-) \mid h(x) = h(0) \quad \text{for } x > 0, w_e \star h(0) = 0\}.$$

For $\phi \in \tilde{V}$, the differential $dJ(\phi) \in \tilde{V}_0'$ is defined weakly by

$$\langle dJ(\phi), h \rangle = \int_{-\infty}^{0} [\phi'(x)h'(x) + N'_{e}(w_{e} \star \phi(x))w_{e} \star h(x) - N'_{i}(w_{i} \star \phi(x))w_{i} \star h(x)] dx, \quad \forall h \in \tilde{V}_{0}.$$

**Definition 2.13.** We call $\phi \in \tilde{V}$ an extremal solution of the gyrokinetic Poisson equation if $dJ(\phi) = 0$.

Let us define now the non linear operators

$$\tilde{n}_e(\phi) = w_e \star (1_{\mathbb{R}^-} - N'_e(w_e \star \phi)) \quad \text{and} \quad \tilde{n}_i(\phi) = w_i \star (1_{\mathbb{R}^-} - N'_i(w_i \star \phi)) \qquad (2.35)$$

where $1_{\mathbb{R}^-}(x) = 1$ for $x \le 0$ and $= 0$ for $x > 0$ is the characteristic function of $\mathbb{R}^-$. Using the physical densities (2.20) and (2.24), we observe that

$$\tilde{n}_e(\phi)(x) = w_e \star (1_{\mathbb{R}^-} - n_e^{\text{phys}})(x) \quad \text{and} \quad \tilde{n}_i(\phi)(x) = w_i \star (1_{\mathbb{R}^-} - n_i^{\text{phys}})(x)$$

so that $\tilde{n}_e(\phi)$ and $\tilde{n}_i(\phi)$ coincide with the gyroaveraged densities $n_e = w_e \star n_e^{\text{phys}}$ and $n_i = w_i \star n_i^{\text{phys}}$ for $x < -1$, see (2.11). Close to the wall they can be seen as a perturbation of the latter.
**Proposition 2.14.** An extremal solution satisfies $\phi' \in C^0(\mathbb{R}^-)$, and the Euler-Lagrange relations

$$ egin{cases} -\phi''(x) + \tilde{n}_e(\phi)(x) - \tilde{n}_i(\phi)(x) + \sigma w_e(x) = 0, & x < 0, \\ \sigma = -2[\phi'(0^-) + \int_{\mathbb{R}^+} (\tilde{n}_e(\phi) - \tilde{n}_i(\phi))(x)dx] & \\ w_e \star \phi(0) = \psi_{\text{wall}} & \end{cases} \qquad (2.36) $$

where $\psi_{\text{wall}} < 0$ is the (electronic) gyroaveraged wall potential defined by Proposition 2.5, and $\sigma \in \mathbb{R}$ is a Lagrange multiplier associated to the boundary constraint.

**Remark 2.15.** For $x < -r_e$ away from the wall, $w_e$ vanishes and the first equation in (2.36) coincides with the initial gyrokinetic Poisson equation (2.13). For the same reason, the perturbed densities $\tilde{n}_e, \tilde{n}_i$ vanish for $x > \max(r_e, r_i)$, so that in the second equation the integral only involves values of the potential close to the wall.

*Proof.* Such formulations for extremal solutions with linear constraints are very classical. The shortest path to the result is achieved with the Lagrangian $L(\phi, \sigma) = J(\phi) + \sigma(w_e \star \phi(0) - \psi_{\text{wall}})$ for which a convenient functional space is $\tilde{W}_0 = W_0 \cap L^1(\mathbb{R}^-)$ with

$$ W_0 = \{h \in \dot{H}^1(\mathbb{R}) \cap L^2(\mathbb{R}^-) \mid h(x) = h(0) \text{ for } x > 0\}. \qquad (2.37) $$

Note that $\tilde{V} \subset \tilde{W}_0$ and $\tilde{V} \neq \tilde{W}_0$. Observing that the energy functional $J$ is continuously defined from $\tilde{W}_0$ into $\mathbb{R}$ (and not only from $\tilde{V}$ as in Lem. 2.9), we see that the Lagrangian $L$ is also continuously defined from $\tilde{W}_0$ into $\mathbb{R}$. Let us calculate the variation with respect to test functions $h \in \tilde{W}_0$

$$ \int_{-\infty}^{0} [\phi'(x)h'(x) + N'_{e}(w_{e} \star \phi(x))w_{e} \star h(x) - N'_{i}(w_{i} \star \phi(x))w_{i} \star h(x)] dx + \sigma \int_{\mathbb{R}} w_{e}(x)h(x)dx = 0. \qquad (2.38) $$

One has

$$ \int_{\mathbb{R}^{-}} \phi'(x)h'(x)dx = -\int_{\mathbb{R}^{-}} \phi''(x)h(x)dx + \phi'(0^{-})h(0), \qquad (2.39) $$

where we note that for $h(0) \neq 0$ this relation is only justified provided $\phi' \in C^0(\mathbb{R}^-)$, then for $s \in \{e, i\}$, using the symmetry of the kernels $w_s$,

$$ \begin{aligned} \int_{\mathbb{R}^{-}} N'_{s}(w_{s} \star \phi(x))w_{s} \star h(x)dx &= \int_{\mathbb{R}} 1_{\mathbb{R}^{-}}(x)N'_{s}(w_{s} \star \phi(x))w_{s} \star h(x)dx = \int_{\mathbb{R}} \tilde{n}_{s}(\phi)(x)h(x)dx \\ &= \int_{\mathbb{R}^{-}} \tilde{n}_{s}(\phi)(x)h(x)dx + \left( \int_{\mathbb{R}^{+}} \tilde{n}_{s}(\phi)(x)dx \right) h(0) \end{aligned} $$

and finally

$$ \int_{\mathbb{R}} w_e(x)h(x)dx = \int_{\mathbb{R}^{-}} w_e(x)h(x)dx + \frac{1}{2}h(0). $$

This yields (with the same reservation as above if $h(0) \neq 0$)

$$ \int_{\mathbb{R}^{-}} [-\phi'' + \tilde{n}_e(\phi) - \tilde{n}_i(\phi) + \sigma w_e](x)h(x)dx + [\phi'(0^{-}) + \int_{\mathbb{R}^{+}} (\tilde{n}_e(\phi) - \tilde{n}_i(\phi))(x)dx + \frac{1}{2}\sigma]h(0) = 0. $$

Taking an arbitrary $h \in \tilde{W}_0$ with $h(0) = 0$ gives the first equation of (2.36). This shows that $\phi'' \in L_{\text{loc}}^1(\mathbb{R}^-)$, so that $\phi' \in C^0(\mathbb{R}^-)$ and (2.39) holds indeed with $h(0) \neq 0$. The second equation of (2.36) is then the remaining coefficient in front of $h(0)$. Finally the boundary condition in (2.36) follows from the variations in $\sigma$. $\square$
### 3. MATHEMATICAL STUDY

To have clearer insight into the existence, uniqueness and qualitative properties of the solution of the physical extremalization problem (2.13), we simplify the problem so as to concentrate on the main features. The main simplification is that the physical characteristic lengths discussed in Section 2.2 are taken all three of them equal to one

$$ \lambda = r_e = r_i = 1. \qquad (3.1) $$

This condition may not be justified on physical grounds, but it allows us to prove that the above model is well-posed. It will be relaxed in Section 4 where numerical experiments will be performed in the regime $r_e < r_i$. We also modify the problem for mathematical simplicity: instead of working in the physical domain consisting of a negative half-axis and a physical electric potential which is negative at the wall, we will now work in the positive half axis $x > 0$ with a positive wall potential $\psi_{\text{wall}} > 0$. This modification is applied also to the functional spaces which are (re)defined accordingly.

Thus, in this section we consider the following simplified mathematical minimization problem:

$$ \min_{\phi \in V} \mathcal{J}(\phi) \qquad (3.2) $$

with an energy functional now defined as

$$ \mathcal{J}(\phi) := \int_{0}^{\infty} \left[ \frac{1}{2} \left| \frac{d}{dx} \phi(x) \right|^2 + F(w \star \phi(x)) \right] dx. \qquad (3.3) $$

The functional space $V$ is redefined as

$$ V := \left\{ \phi \in \dot{H}^1(\mathbb{R}) \cap L^2(\mathbb{R}^+) : w \star \phi(0) = \psi_{\text{wall}} > 0 \quad \text{and} \quad \phi(x) = \phi(0) \quad \text{for} \ x < 0 \right\}. \qquad (3.4) $$

Here the parameters are the function $F: \mathbb{R} \to \mathbb{R}$ whose properties will be discussed below, the gyroaveraging kernel $w$ defined in (2.6) and the wall potential $\psi_{\text{wall}}$, which is now assumed positive for mathematical simplicity, as specified above. We also remind that the space $\dot{H}^1$ is defined as the set of measurable functions $\phi$ which derivative is in $L^2$, and that functions $\phi$ in the affine space $V$ are such that $\lim_{+\infty} \phi = 0$. The tangent space is

$$ V_0 := \left\{ \phi \in \dot{H}^1(\mathbb{R}) \cap L^2(\mathbb{R}^+) : w \star \phi(0) = 0 \quad \text{and} \quad \phi(x) = \phi(0) \quad \text{for} \ x < 0 \right\}. \qquad (3.5) $$

Note that if the function $w$ is replaced by a Dirac mass, we recover a standard nonlinear Poisson problem with a Dirichlet constraint at $x = 0$.

With our simplification (3.1), the function $F: \mathbb{R} \to \mathbb{R}$ corresponds to the difference $N_e - N_i$ from the modeling Section 2.8.1. Thus it should gather some main properties of the density potentials $N_e$ and $N_i$. For the mathematical analysis we assume that it satisfies the following hypotheses:

* $F(0) = 0$; \qquad (3.6)

* $F'$ is increasing; \qquad (3.7)

* $F''$ is in $L_{\text{loc}}^1(\mathbb{R})$; \qquad (3.8)

* $\exists 0 < \alpha < \beta, \forall s \in \mathbb{R}, \min\{\alpha s; \beta s\} \le F'(s) \le \max\{\alpha s; \beta s\}$. \qquad (3.9)

The condition (3.6) expresses the physical condition $N_i(0) = N_e(0) = 0$. The condition $F'(0) = 0$ expresses the charge neutrality at infinity (2.33). The monotonicity (3.7) is a generalization of the physical stability condition at infinity (2.34). Condition (3.9) means that the derivative $F'$ lays between the curves $y = \alpha x$ and $y = \beta x$ as illustrated in Figure 4. We also draw in Figure 4 the fact that the function $F''$ may have a mild singularity at $\psi_{\text{wall}}$ since it is the case for the physical model, see equations (2.28)–(2.29). Since the singularity there scales like
<figure><img src="image_6.png" /><figcaption>FIGURE 4. Condition (3.9) states that the function $F''$ may show integrable singularities, as it is the case for the physical model in the vicinity of the wall at the value $\phi_{\text{wall}}$.</figcaption></figure>

the inverse of a square root, it is locally integrable which gives a physical justification for the requirement (3.8). The condition (3.9) is justified in the vicinity of the null electric potential $s = \phi = 0$ by the local smoothness of the functions $N_i$ and $N_e$ and the local convexity of $N_e - N_i$. This insures the coercivity of the function $\mathcal{J}$. Since the physical problem satisfies the four conditions (3.6)–(3.9), they are not a restriction. Finally all conditions (3.6)–(3.9) imply the lower boundedness and upper boundedness inequalities for all $s \in \mathbb{R}$

$$ \alpha s^2 \le sF'(s) \le \beta s^2 \quad \text{and} \quad \frac{\alpha}{2}s^2 \le F(s) \le \frac{\beta}{2}s^2. \qquad (3.10) $$

## 3.1. Main results

In this section we state our main results concerning the minimization problem defined by (3.2), (3.3), and (3.4). For the clarity of the presentation, the proofs of the latter are postponed to Appendix B. Essentially we will show that the minimization problem has a unique solution, and also prove that for certain parameters the solution admits an infinite number of oscillations as it propagates towards the core of the plasma. This oscillating behavior has already been observed in numerical simulations [16], it is related with the gyroaverage operator.

### 3.1.1. Existence, uniqueness, first properties

For the sake of readability, the proofs of the two next propositions are postponed in the appendix. They are rather standard for such a problem.

**Proposition 3.1** (Existence and uniqueness of the minimizer $\phi$). The minimization problem given by (3.2), (3.3) and (3.4) admits a unique solution. This solution $\phi \in V$ satisfies $\phi' \in C^0(\mathbb{R}^+)$ and it is characterized by the Euler-Lagrange equations

$$ \forall x > 0, \quad -\phi''(x) + w \star (1_{\mathbb{R}^+} F'(w \star \phi))(x) = \sigma w(x), \qquad (3.11) $$

$$ \text{where} \quad \frac{\sigma}{2} = -\phi'(0^+) + \int_{-1}^{0} w \star (1_{\mathbb{R}^+} F'(w \star \phi))(x) dx, \qquad (3.12) $$

$$ \text{and} \quad w \star \phi(0) = \psi_{\text{wall}} > 0. \qquad (3.13) $$

The real number $\sigma$ is the Lagrange multiplier associated to the constraint at the wall (3.13).

These equations correspond to (2.36) rewritten in the setting of Section 3. In particular Remark 2.15 also applies here, with the simplification (3.1) namely $r_e = r_i = 1$. Concerning the function $\phi$, minimizer of $\mathcal{J}$ on $V$, we are able to extract the following elementary properties.
**Proposition 3.2** (Basic properties of the minimizer $\phi$). Let $\phi \in V$ be the solution to the minimizing problem (3.2), (3.3) and (3.4). Then the solution $\phi$ satisfy the following properties.

(i) One has an estimate of the $H^1$ semi-norm of $\phi$:

$$ \alpha \int_{-1}^{\infty} |w \star \phi(x)|^2 dx \leq \int_{-1}^{\infty} |\phi'(x)|^2 dx - \sigma \int_{-1}^{1} w(x) \phi(x) dx \leq \beta \int_{-1}^{\infty} |w \star \phi(x)|^2 dx. $$

(ii) With $\beta$ defined at (3.9), the energy satisfies $\mathcal{J}(\phi) \leq 2 \psi_{\text{wall}}^2 \sqrt{\beta}(1 + \beta)$.

(iii) The following identity holds $\phi'(0^+) = \frac{1}{2} \int_0^1 g(-x) - g(x) dx$ where $g = w \star (1_{\mathbb{R}^+} F'(w \star \phi))$.

(iv) The Lagrange multiplier is positive $\sigma > 0$. For $x \to 1^-$, one has the equivalent $\frac{\phi''(x)}{\sigma w(x)} \to 1$.

(v) The function $\phi$ is not identically equal to 0 on the interval $[1; +\infty)$.

### 3.1.2. Local monotonicity of $w \star \phi$ near the wall

The property of local monotonicity near the wall guarantees that, at least in the vicinity of the wall, the function $w \star \phi$ takes values below $\psi_{\text{wall}}$. This property is important because the physical modeling of Section 2.5 give a priori zero information for values above $\psi_{\text{wall}}$. In other words, the local monotonicity property states that the solution takes only physical values (in the vicinity of the wall). At infinity, the solution converges to zero with oscillation, as explained in the next two sections, but the oscillations are asymptotically in bounds with the admissible domain (for the physical model in Sect. 2.6, the fact that $\psi^+ > 0$ also helps).

**Proposition 3.3.** At $x = 0^+$ the gyroaverage of $\phi$ satisfies $\frac{d}{dx}(w \star \phi)(0^+) < 0$.

*Proof.* One considers the following test function:

$$ \chi := \gamma\phi + \phi'(x), \qquad (3.14) $$

where the coefficient $\gamma$ is chosen such that

$$ \gamma w \star \phi(0) + w \star (\phi')(0^+) = 0. \qquad (3.15) $$

One must take care that $\chi \notin \dot{H}^1$ because the right hand side in (3.11) is not square integrable (indeed $\sigma > 0$ and $w \in L^p(\mathbb{R})$ for $1 \leq p < 2$). Nevertheless one performs the calculation using a preliminary regularization argument (not detailed). Thus,

$$ \begin{aligned} 0 &= \int_{0}^{\infty} \phi'(x)\chi'(x) dx + \int_{0}^{\infty} F'(w \star \phi(x)) w \star \chi(x) dx \\ &= \gamma \int_{0}^{\infty} |\phi'(x)|^2 dx + \int_{0}^{\infty} \phi'(x)\phi''(x) dx \\ &\quad + \gamma \int_{0}^{\infty} F'(w \star \phi(x)) w \star \phi(x) dx + \int_{0}^{\infty} F'(w \star \phi(x)) w \star \phi'(x) dx, \end{aligned} \qquad (3.16) $$

where for the second equality one uses (3.14). Using now the definition of $\gamma$ given at (3.15), equation (3.16) becomes

$$ \begin{aligned} 0 = &-\frac{w \star \phi'(0^+)}{\psi_{\text{wall}}} \left( \int_{0}^{\infty} |\phi'(x)|^2 dx + \int_{0}^{\infty} F'(w \star \phi(x)) w \star \phi(x) dx \right) \\ &+ \int_{0}^{\infty} \phi'(x) \phi''(x) dx + \int_{0}^{\infty} F'(w \star \phi(x)) w \star \phi'(x) dx. \end{aligned} $$

Then one is led to

$$ 0 = -\frac{w \star \phi'(0^+)}{\psi_{\text{wall}}} \left( \int_{0}^{\infty} |\phi'(x)|^2 dx + \int_{0}^{\infty} F'(w \star \phi(x)) w \star \phi(x) dx \right) - \frac{|\phi'(0^+)|^2}{2} - F(\psi_{\text{wall}}). $$

Observing now that assumption (3.9) implies $F'(s)s \geq 0$ for all $s \in \mathbb{R}$ and also $F(s) \geq 0$ with (3.6), one concludes with the above equality that $w \star \phi'(0^+) < 0$. $\square$
### 3.1.3. Oscillatory nature of the solution for large $\alpha$

The analysis of the solution near the wall is rather complicated but we expect it to be concave decreasing – at least on a neighborhood of the wall. The difficulty comes from the presence of the Lagrange multiplier $\sigma$ in equation (3.11), which is active only for $x < 1$ since the support of $w$ is $[-1, 1]$. Concerning the behavior of the solution when the distance to the wall is larger than 1, we are able to extract the following property when $\alpha$ is large. We remind that, by Proposition 3.2-(v), there exists $x_0 \in [1; +\infty)$ such that $\phi(x_0) \neq 0$.

**Proposition 3.4** (Oscillatory nature of the solution). Assume that $\alpha$ defined at (3.10) satisfies $\alpha \ge 5\pi^2$. Take $x_0 \in [1; +\infty)$ such that $\phi(x_0) \neq 0$. Then there exists $x_1 \in [x_0; x_0 + 10]$ such that $\phi(x_0) \phi(x_1) < 0$.

**Remark 3.5.** From this proposition, we conclude that the solution $\phi$ oscillates infinitely many times around the value $\phi = 0$. Indeed, by using this proposition recursively we get an increasing sequence $(x_n)_{n \in \mathbb{N}}$ such that for all index $n$ we have $\phi(x_n) \phi(x_{n+1}) < 0$. Proposition 3.2-(v) is important to initialize this recursive argument. We also emphasize on the fact that this proposition gives an upper bound on the wave length of the oscillations. Indeed, it is required that $x_1$ is lower than $x_0 + 10$.

This phenomenon generates new open questions. What is the limit value of $\alpha$? Is it possible to characterize the amplitude and wave length of the oscillations? Some answers are possible with the method of exact solutions of next Section 3.1.4.

*Proof.** • Step 1: Let $x_0 \in [1; +\infty)$ such that $\phi(x_0) \neq 0$. Without loss of generality, we make the assumption that $\phi(x_0) > 0$ because the following reasoning works the same in the case $\phi(x_0) < 0$, by simply changing all the signs. Suppose for the sake of contradiction that

$$ \forall x \in [x_0; x_0 + 10], \quad \phi(x_0) \ge 0, \qquad (3.17) $$

We now construct a particular admissible direction $h \in V_0$ (with $V_0$ defined at (3.5)). We ask the function $h \in V_0 \cap C^1(\mathbb{R})$ to be non-negative, bounded, with support $\text{supp} h = [x_0 + 4; x_0 + 6]$ and such that $h'' \in L^\infty$. There exists an infinite choice of such functions. Since $w$ is also non negative and has its support is a segment of size 2, one has

$$ \text{supp } w \star h = [x_0 + 2; x_0 + 8], \quad \text{and} \quad \text{supp } w \star w \star h = [x_0; x_0 + 10]. \qquad (3.18) $$

One now studies the variations of the functional $\mathcal{J}$ in the direction $h$. One gets by a direct computation

$$ \lim_{s \to 0} \frac{\mathcal{J}(\phi + sh) - \mathcal{J}(\phi)}{s} = \int_{x_0+4}^{x_0+6} \phi'(x) h'(x) dx + \int_{x_0+2}^{x_0+8} F'(w \star \phi)(x) (w \star h)(x) dx. \qquad (3.19) $$

Since the function $\phi$ is supposed to be non-negative on $[x_0; x_0+10]$, then $w \star \phi$ is non-negative on $[x_0+2; x_0+8]$. One also has $h$ is non-negative. The property (3.9) combined with (3.18) yields

$$ \begin{aligned} \int_{x_0+2}^{x_0+8} F'(w \star \phi)(x) (w \star h)(x) dx &\ge \alpha \int_{x_0+2}^{x_0+8} (w \star \phi)(x) (w \star h)(x) dx \\ &= \alpha \int_{\mathbb{R}} (w \star \phi)(x) (w \star h)(x) dx \\ &= \alpha \int_{\mathbb{R}} \phi(x) (w \star w \star h)(x) dx \\ &= \alpha \int_{x_0}^{x_0+10} \phi(x) (w \star w \star h)(x) dx. \end{aligned} \qquad (3.20) $$

On the other hand, an integration by part gives

$$ \int_{x_0+4}^{x_0+6} \phi'(x) h'(x) dx = - \int_{x_0+4}^{x_0+6} \phi(x) h''(x) dx. $$

Therefore (3.19) implies

$$ \lim_{s \to 0} \frac{\mathcal{J}(\phi + sh) - \mathcal{J}(\phi)}{s} \ge \int_{x_0}^{x_0+10} \phi(x) \left[ \alpha(w \star w \star h)(x) - h''(x) \right] dx. \qquad (3.21) $$
• **Step 2:** Given the function $h$, one defines the quantity

$$ \mu := \sup_{x \in [x_0; x_0+10[} \frac{h''(x)}{w \star w \star h(x)} = \sup_{x \in [x_0+4; x_0+6]} \frac{h''(x)}{w \star w \star h(x)} \quad (3.22) $$

because $w \star w \star h$ is positive on $[x_0; x_0 + 6]$ as stated in (3.18), whereas $h''$ has its support contained in the support of $h$ which is $[x_0 + 4; x_0 + 6]$. These two facts imply that the quantity $\mu$ is finite. On the other hand, one has $0 = h'(x_0 + 6) - h'(x_0 + 4) = \int_{x_0+4}^{x_0+6} h''(x) dx$. Since $h''$ is not identically 0 on the support of $h$, then the equation above gives that

meas $\{x \in [x_0 + 4; x_0 + 6]: h''(x) > 0\}$ is positive.

One concludes that $\mu > 0$. One now uses the definition of $\mu$ to study the right-hand side of equation (3.21). One has $h''(x) \le \mu(w \star w \star h)(x)$ for all $x \in [x_0; x_0 + 10]$. Therefore,

$$ \int_{x_0}^{x_0+10} \phi(x) [\alpha(w \star w \star h)(x) - h''(x)] dx \ge (\alpha - \mu) \int_{x_0}^{x_0+10} \phi(x)(w \star w \star h)(x) dx, \quad (3.23) $$

where again $\phi \ge 0$ on $[x_0; x_0 + 10]$. One now considers $\alpha$ very large so that $\alpha > \mu$ and obtains $(\alpha - \mu)(w \star (w \star h))(x) > 0$ for all $x \in [x_0; x_0 + 10]$. Moreover, as a consequence of (3.17) the function $\phi$ is non-negative on $[x_0; x_0 + 10]$ but it is not identically 0 since $\phi(x_0) > 0$. Thus, the right-hand side of (3.23) is positive. Plugging this back into (3.21) gives

$$ \lim_{s \to 0} \frac{\mathcal{J}(\phi + sh) - \mathcal{J}(\phi)}{s} > 0, $$

which eventually contradicts the minimality of the function $\phi$. It already proves the claim for $\alpha$ large enough. It remains to show the lower bound on $\alpha$.

• **Step 3:** Here, we prove that $\alpha > 5\pi^2$. The idea here is to consider a particular function $h$ for the formula (3.21) and to provide an estimate on its associated coefficient $\mu$.

• Firstly one defines a function $\psi: \mathbb{R} \to \mathbb{R}^+$ (this notation has nothing to do with the previous function $\psi$ in the modeling section) as follows

$$ \forall x \in \mathbb{R}, \quad \psi(x) = \begin{cases} 0 & \text{if } x \le 0, \\ x^2/2 & \text{if } 0 \le x \le 1/2, \\ -(x-1)^2/2 + 1/4 & \text{if } 1/2 \le x \le 3/2, \\ (x-2)^2/2 & \text{if } 3/2 \le x \le 2, \\ 0 & \text{if } 2 \le x. \end{cases} \quad (3.24) $$

By definition, one has that $\psi \in W^{2,\infty}(\mathbb{R})$ is supported on $[0; 2]$, is non-negative and is such that

$$ \psi''(x) = 1_{(0; \frac{1}{2})}(x) - 1_{(\frac{1}{2}; \frac{3}{2})}(x) + 1_{(\frac{3}{2}; 2)}(x) \text{ almost everywhere.} \quad (3.25) $$

Since $w$ is lower bounded by $\frac{1}{\pi}$ on the interval $[-1, 1]$, one has

$$ w \star w \star \psi(x) = \int_{x-1}^{x+1} \int_{y-1}^{y+1} w(x-y) w(y-z) \psi(z) dy dz \ge \frac{1}{\pi^2} \int_{x-1}^{x+1} \int_{y-1}^{y+1} \psi(z) dy dz $$

because $|x - y| \le 1$ and $|y - z| \le 1$. One is now going to compute explicitly this last integral. One writes it as

$$ \frac{1}{\pi^2} \iint_{D_x} \psi(z) dz dy \quad (3.26) $$
<figure><img src="image_7.png" /><figcaption>FIGURE 5. An illustration of the set $D_x \subseteq \mathbb{R}^2$, subdivided into $D_{x,1} \cdots D_{x,6}$ represented by the different shades of grey. This subdivision intervenes in the computation of the integral (3.26).</figcaption></figure>

where $D_x := \{(y, z) \in \mathbb{R}^2 : x - 1 \le y \le x + 1 \text{ and } y - 1 \le z \le y + 1\}$. Subdivide the set $D_x$ into the six following subsets

$$ 
\begin{align} 
D_{x,1} &:= D_x \cap \{(y,z) \in \mathbb{R}^2 : z \le 0\}, & D_{x,2} &:= D_x \cap \{(y,z) \in \mathbb{R}^2 : 0 \le z \le x\}, \nonumber \\
D_{x,3} &:= D_x \cap \{(y,z) \in \mathbb{R}^2 : x \le z \le 1/2\}, & D_{x,4} &:= D_x \cap \{(y,z) \in \mathbb{R}^2 : 1/2 \le z \le 3/2\}, \tag{3.27} \\
D_{x,5} &:= D_x \cap \{(y,z) \in \mathbb{R}^2 : 3/2 \le z \le 2\}, & D_{x,6} &:= D_x \cap \{(y,z) \in \mathbb{R}^2 : 2 \le z\}. \nonumber
\end{align} 
$$

One computes $\iint_{D_x} \psi(z)dz dy$ by splitting the computation on these 6 sub-domains and using directly (3.24). First, one has

$$ \iint_{D_{x,1}} \psi(z) dz dy = \iint_{D_{x,6}} \psi(z) dz dy = 0 $$

One gets the 4 remaining terms by a direct computation using the Fubini theorem, that is in order

$$ 
\begin{align*}
\iint_{D_{x,2}} \psi(z) dz dy &= \int_0^x \int_{x-1}^{z+1} \psi(z) dy dz = \int_0^x \int_{x-1}^{z+1} \frac{z^2}{2} dy dz = \frac{x^3}{3} - \frac{x^4}{24}, \\
\iint_{D_{x,3}} \psi(z) dz dy &= \int_x^{1/2} \int_{z-1}^{x+1} \psi(z) dy dz = \int_x^{1/2} \int_{z-1}^{x+1} \frac{z^2}{2} dy dz = -\frac{x^4}{24} - \frac{x^3}{3} + \frac{x}{48} + \left(\frac{1}{24} - \frac{1}{128}\right), \\
\iint_{D_{x,4}} \psi(z) dz dy &= \int_{1/2}^{3/2} \int_{z-1}^{x+1} \psi(z) dy dz = \int_{1/2}^{3/2} \int_{z-1}^{x+1} \left[\frac{1}{4} - \frac{(z-1)^2}{2}\right] dy dz = \frac{5}{24}(x+1)
\end{align*} 
$$
and

$$ \iint_{D_{x,5}} \psi(z) dz dy = \int_{3/2}^{2} \int_{z-1}^{x+1} \psi(z) dy dz = \int_{3/2}^{2} \int_{z-1}^{x+1} \frac{(z-2)^2}{2} dy dz = \frac{x}{48} - \frac{1}{128}. $$

Thus one concludes that $\iint_{D_x} \psi(z) dz dy = -\frac{x^4}{12} + \frac{x}{4} + \frac{15}{64}$. An analysis of the function $x \in [0, 1/2] \mapsto -x^4/12+x/4$ shows that it is non-negative and thus $\iint_{D_x} \psi(z) dz dy \ge \frac{15}{64} > \frac{1}{5}$. Plugging this back into (3.26) gives

$$ \forall x \in [0; \frac{1}{2}], \quad (w \star w \star \psi)(x) > \frac{1}{5\pi^2} \qquad (3.28) $$

Invoking now the fact that $\psi$ is symmetrical with respect to $x = 1$, one can do the very same reasoning with $x \in [3/2; 2]$. Therefore

$$ \forall x \in [\frac{3}{2}; 2], \quad (w \star w \star \psi)(x) > \frac{1}{5\pi^2} \qquad (3.29) $$

Gathering now the estimates (3.28) and (3.29), one concludes

$$ \forall x \in \mathbb{R}, \quad (w \star w \star \psi)(x) > \frac{1}{5\pi^2} \psi''(x), \qquad (3.30) $$

where $\psi''$ was computed at (3.25).

* Finally one defines $h(x) := \psi(x - x_0 - 4)$. It is straight-forward to check that $h$ is in $V_0$, non-negative, supported on $[x_0 + 4; x_0 + 6]$ and satisfy $h'' \in L^\infty(\mathbb{R})$. As a consequence of (3.30), the coefficient $\mu$ defined at (3.22) associated to this function $h_0$ is smaller than $5\pi^2$. It shows that $\alpha \ge 5\pi^2$. □

### 3.1.4. Linearization of the equation far from the wall

The equation (3.11) linearized at infinity (that is far from the wall) around the asymptotic value $\phi = 0$ yields the linear problem. It is also a particular case of the general situation, since it is sufficient to take $\alpha = \beta$ in (3.10). It appears that the linear problem helps to explain the oscillating behavior of the solution by means of the construction of exact solutions.

The homogeneous linear equation that we consider is written as

$$ -u''(y) + \gamma w \star w \star u(y) = 0, \quad y > 0. \qquad (3.31) $$

By compatibility with (3.10), the coefficient $\gamma$ is positive with $\beta \ge \gamma \ge \alpha > 0$. If the convolution operator would be atomic (a Dirac mass), then the linearized equation would be the classical Laplace equation $-u'' + \gamma u = 0$ with classical particular real solutions $u_r(y) = \exp(\pm\sqrt{\gamma}y)$. But the situation is now very different due to the convolution operator.

Let us look for a general complex exponential solution $u_r(y) = \exp(-ry)$ with $r \in \mathbb{C}$. One has by a direct calculation that $w \star u_r(y) = \int_{\mathbb{R}} w(x) \exp(-r(y-x))dx = \int_{\mathbb{R}} w(x) \exp(rx)dx \exp(-ry) = C(r)u_r(y)$ with $C(r) = w \star u_r(0)$. Plugging in (3.31), one obtains the compatibility relation $r^2 = \gamma C(r)^2$, that is $r = \pm\sqrt{\gamma} C(r)$.

**Lemma 3.6.** One has $C(r) = I_0(r)$ where $I_0$ is the first modified Bessel function.

*Proof.* One starts from $C(r) = \frac{1}{\pi} \int_{-1}^{1} \exp(rx) \frac{1}{\sqrt{1-x^2}} dx$. By derivation and integration by parts, one gets $C'(r) = \frac{1}{\pi} \int_{-1}^{1} \exp(rx) \frac{x}{\sqrt{1-x^2}} dx = \frac{r}{\pi} \int_{-1}^{1} \exp(rx) \sqrt{1-x^2} dx$. Another derivation yields $C''(r) = \frac{1}{\pi} \int_{-1}^{1} \exp(rx) \frac{x^2}{\sqrt{1-x^2}} dx$. One can also write $C(r) - C''(r) = \frac{1}{\pi} \int_{-1}^{1} \exp(rx) \frac{1-x^2}{\sqrt{1-x^2}} dx = \frac{1}{\pi} \int_{-1}^{1} \exp(rx) \sqrt{1-x^2} dx$ so one obtains the differential equation

$$ C(r) - C''(r) = \frac{1}{r} C'(r) \iff r^2 C''(r) + rC'(r) - r^2 C(r) = 0. $$

This is the modified Bessel equation, see Formula 10.25.1 of p. 248 from [37]. Therefore $C = aI_0 + bK_0$ where $a, b \in \mathbb{R}$. Since $C(0) = 1$ and $C'(0) = 0$, one gets $a = 1$ and $b = 0$ and the proof is ended. □
<figure><img src="image_8.png" /><figcaption>FIGURE 6. Intersection of $x \mapsto x$ with $x \mapsto tI_0(x)$ for three different values of $t$. The correspondence is: $f$ is for $t > t_*$, $g$ is for $t = t_*$ and $h$ is for $t < t_*$.</figcaption></figure>

The first modified Bessel function $I_0$ is even, so it is sufficient to study positive solutions $z > 0$

$$z = tI_0(z), \quad t = \sqrt{\gamma} > 0. \tag{3.32}$$

The functions $u_r(y) = \exp(-ry)$ solution to (3.31) that we study correspond to $r = \pm z$.

For real $x \in \mathbb{R}$, the graph of the function $x \mapsto I_0(x)$ is strictly convex and lower bounded by a quadratic, that is $I_0(x) \ge 1 + cx^2$ with $c = \frac{1}{2}I_0''(0) > 0$. Therefore there exists $t_* > 0$ such that the graph of $x \mapsto t_*I_0(x)$ is tangent to the graph $x \mapsto x$. It is characterized by the transcendental equations

$$\begin{cases} x_* = t_* I_0(x_*), \\ 1 = t_* I'_0(x_*). \end{cases}$$

**Lemma 3.7.** For $0 < t < t_*$, there exists two real solutions $0 < x_-(t) < x_* < x_+(t)$ of the equation $x = tI_0(x)$.

*Proof.* The proof is by intersection of graphs, as in the Figure 6. $\square$

Now we study equation (3.32) for $t > t_*$, and for possible complex solutions $z \in \mathbb{C}$.

**Lemma 3.8.** One can rewrite equation (3.32) under the form

$$z - x_* = \pm \sqrt{\frac{t_* - t}{t_*}} H(z - x_*) \tag{3.33}$$

where $H$ is an analytic series (with a certain radius of convergence) with real coefficients.

*Proof.* One rewrites equation (3.32) as

$$z - x_* = t_* (I_0(z) - I_0(x_*)) + (t - t_*) I_0(z),$$

then $(t_* \frac{I_0(z) - I_0(x_*)}{z - x_*} - 1)(z - x_*) = (t_* - t)I_0(z)$, then $(t_* \frac{I_0(z) - I_0(x_*)}{z - x_*} - t_* I'_0(x_*))(z - x_*) = (t_* - t)I_0(z)$, then

$$(z - x_*)^2 = \frac{t_* - t}{t_*} \times \frac{I_0(z)}{\frac{I_0(z) - I_0(x_*) - I'_0(x_*)(z - x_*)}{(z - x_*)^2}} \tag{3.34}$$
By continuity, the function $G(z) = \frac{I_0(z)}{\frac{I_0(z)-I_0(x_*)-I_0'(x_*)(z-x_*)}{(z-x_*)^2}}$ admit a finite value at $z = x_*$, namely $G(x_*) = 2\frac{I_0(x_*)}{I_0'(x_*)} > 0$. The function $I_0$ is an analytic series in the whole complex plane with real coefficients so $G(z)$ is also an analytic series in the whole complex plane with real coefficients. Therefore one can take the square root of $G$

$$H(z) = \sqrt{G(z)} = h_0 \sqrt{1 + \frac{G(z) - G(x_*)}{h_0^2}} \quad \text{where} \quad h_0 = \sqrt{G(x_*)} > 0.$$

This is an analytic series $H(z) = \sum_{n \ge 0} h_n (z - x_*)^n$ in function of $z - x_*$ with real coefficients and with a certain radius of convergence (because of the square root). Therefore taking the square root of (3.34) yields the claim. $\square$

It is now evident to study the solution of equation (3.32), at least for $|t - t_*|$ small enough.

First case $0 < t < t_*$ The square root $\sqrt{t_* - t} \in \mathbb{R}$ takes real values. Then the equation can be solved with a fixed point method for $z \in \mathbb{R}$. In this case one recovers the two real solutions of Lemma 3.7.

Second case $t_* < t$ The square root $\sqrt{t_* - t} \in \mathbb{R}$ takes imaginary values. Then it is better to write the equation as

$$z - x_* = \pm i \sqrt{\frac{|t_* - t|}{t_*}} H(z - x_*).$$

A fixed point method yields two solutions in the vicinity of $x_*$, but in the complex plane, that is $z = a + ib$ with $a, b \in \mathbb{R}$ and $b \neq 0$.

In summary of this discussion is that the solutions of (3.32) are either real or complex.

- If $t > t_*$, one can rewrite complex solutions of the linearized equation (3.31) as real solutions under the form

$$u(y) = \exp(-ay) \cos(by), \quad \text{where } z = a + ib \text{ solves (3.33) with } a > 0, b \neq 0, \text{ and } y > 0.$$

This solution is clearly oscillating.

- If $t < t_*$, it is also possible to construct solutions that change sign under the form

$$u(y) = a \exp(-x_-(t)y) + b \exp(-x_+(t)y), \quad b > -a > 0, \quad y > 0.$$

- If $t = t_*$, it is also possible to construct solutions that change sign under the form

$$u(y) = at \exp(-x_*y) + b \exp(-x_*y), \quad b, -a > 0, \quad y > 0.$$

- Since $w \star u = C(r)u$ and $w \star w \star u = C(r)^2u$, then the changes of sign hold also for the convolved functions $w \star u$ and $w \star w \star u$. So the convolution operator has no smoothing effect for these functions.

The construction of this section and the result of Proposition 3.4 exemplify the fact that the behavior of the mathematical solutions is deeply affected by the convolution operator. To our understanding, it seems that this phenomenon has not been studied in the mathematical literature [19, 22, 25, 30].
### 3.1.5. A physical interpretation of the oscillatory solutions

The two methods explained in Sections 3.1.3 and 3.1.4 exhibit oscillatory solutions for $\alpha > 0$ large enough, which means for a function $F$ with a magnitude large enough. Going back to the physical equations (2.7)-(2.9) and to the physical scaling (2.10) of the Debye length $\lambda$, a natural physical interpretation is that the physical Debye length is small enough in relative units. We arrive at the conclusion that the gyroaveraging procedure has the ability to generate oscillatory mathematical solutions for the range of parameters considered in physical papers such as [11]. We also note that oscillatory potentials have already been observed in numerical simulations [16], where they are associated to the so-called magnetic pre-sheath.

## 4. NUMERICAL STUDY

The objective hereafter is, on the one hand to explain that it is possible to discretize the minimization/extremalization problem with standard numerical methods, and the other hand to show that the oscillatory/non oscillatory behavior of the solutions predicted in the theoretical section is observed. In this section we use the sign convention of the physical modeling Section 2: we work in the negative half line $x < 0$, with negative wall potentials $\psi_{\text{wall}} < 0$.

### 4.1. Description of the numerical method

The numerical method is made of two steps.

- The first step consists in assembling the density potentials $N_i$ and $N_e$ and computing the reference density $n_{\text{ref}} > 0$ and the Dirichlet boundary condition on the wall potential $\psi_{\text{wall}} < 0$ by solving (2.26)-(2.27). We will consider a function $f_i^\infty$ that is piecewise continuous, so that standard quadrature formulas for velocity integrals can be used to assemble the potential $N_i$. The other potential poses no additional difficulty. Next, $n_{\text{ref}} > 0$ and $\psi_{\text{wall}} < 0$ are computed by means of a Newton-Raphson method applied to the function $W: 0 < u \mapsto H(u) - \frac{\sqrt{2n_i^\infty}}{\sqrt{\pi\eta}J_i^\infty}$, see Proposition 2.5. Provided $\sqrt{\pi\eta}J_i^\infty < \sqrt{2n_i^\infty}$, we know that this function has a unique zero, $W(-\psi_{\text{wall}}) = 0$.

- The second step consists in solving the minimization problem (3.2), (3.3) and (3.4) by means of a gradient method based on a finite element approximation. Since the solution $\phi$ to the minimization problem (3.2), (3.3) and (3.4) belongs to $H^1(\mathbb{R}^-)$, it has a zero limit at $-\infty$, that is $\lim_{x \to -\infty} \phi = 0$. As a consequence, for the numerical method we take $L > 0$ large enough such that for $x < -L$, then $|\phi(x)| < \varepsilon_{\text{machine}}$ where $\varepsilon_{\text{machine}}$ corresponds to the relative machine error. Given $N_x \in \mathbb{N}$, we subdivide the computational domain $[-L, 0]$ in $N_x + 1$ intervals $I_k = [x_k, x_{k+1}]$ of uniform size $h = \frac{L}{N_x+1}$, We then consider the finite element space $W_0^h \subset W_0$ made of globally continuous and piecewise affine functions which are constant outside the domain $[-L, 0]$. That is, we set

$$W_0^h = \{v_h \in C^0(\mathbb{R}) : v_h|_{I_k} \in \mathbb{P}_1(I_k) \forall k = 0, \dots, N_x \text{ and } v_h(x) = 0 \text{ for } x \le -L, v_h(x) = v_h(0) \text{ for } x > 0\}.$$

As for the Dirichlet boundary condition at the wall, we use a standard penalty method [23]. We consider $\varepsilon > 0$ small enough and the following penalized minimization problem on $W_0^h$: find $\phi_h \in W_0^h$ such that

$$\phi_h^\varepsilon = \arg \min_{v_h \in W_0^h} L_\varepsilon(v_h) \quad (4.1)$$

where the Lagrangian $L_\varepsilon$ is defined for all $v_h \in W_0^h$ by

$$L_\varepsilon(v_h) = J(v_h) + \frac{1}{2\varepsilon}(w_e \star v_h(0) - \psi_{\text{wall}})^2$$

with $J$ the energy functional from (2.32). For any $\varepsilon > 0$ the minimization problem (4.1) has a unique solution $\phi_h^\varepsilon$ since $W_0^h$ is a finite dimensional subspace of $W_0$ and $L_\varepsilon$ is coercive and strictly convex. It is standard that
<figure><img src="image_9.png" /><figcaption>FIGURE 7. Plot of the electrostatic potential and gyroaveraged electrostatic potential.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIGURE 8. Plot in log-scale of the electrostatic potential and gyroaveraged electrostatic potential. The drop at $x \approx -15$ is due to machine precision, because the boundary condition was numerically set to $\phi = 0$ on the left boundary to emulate the condition at infinity.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>FIGURE 9. Plot of the ionic density and electronic density.</figcaption></figure>

$\phi_h^\epsilon \to \phi_h$ in $H^1(\mathbb{R}^-)$ as $\epsilon \to 0^+$ (see Thm. 7.1 of [23]) where $\phi_h$ is the unique minimizer of $J$ on $W_0^h$ such that $w_e \star \phi_h(0) = \psi_{\text{wall}}$. To compute $\phi_h^\epsilon$, we use a standard gradient method which consists in computing recursively the sequence $(\phi_h^{\epsilon,n})_{n \in \mathbb{N}} \subset W_0^h$ given by

$$ \begin{cases} \phi_h^{\epsilon,0} \in W_0^h \\ \phi_h^{\epsilon,n+1} = \phi_h^{\epsilon,n} - \rho \nabla L_{\epsilon}(\phi_h^{\epsilon,n}) \quad \forall n \in \mathbb{N}, \end{cases} $$

where $\rho > 0$ is a given parameter.

By the Riesz representation theorem, the gradient $\nabla L_{\epsilon}(\phi_h^{\epsilon,n}) \in W_0^h$ is the unique solution of the variational problem

$$ (\nabla L_{\epsilon}(\phi_h^{\epsilon,n}), v_h)_{H^1(\mathbb{R}^-)} = dL_{\epsilon}(\phi_h^{\epsilon,n})(v_h) \quad \forall v_h \in W_0^h, $$

where $(\cdot, \cdot)_{H^1(\mathbb{R}^-)}$ is the usual inner product on $H^1(\mathbb{R}^-)$ and $dL_{\epsilon}(\phi_h^{\epsilon,n})$ is the Fréchet derivative of $L_{\epsilon}$ at $\phi_h^{\epsilon,n}$,

$$ dL_{\epsilon}(\phi_h^{\epsilon,n}) : v_h \mapsto dJ(\phi_h^{\epsilon,n})(v_h) + \frac{1}{\epsilon}(w_e \star \phi_h^{\epsilon,n}(0) - \psi_{\text{wall}})(w_e \star v_h(0)). $$

To apply the gyroaverage operator (2.5) we use a Gauss-Chebyshev quadrature which is well suited for the numerical treatment of the convolution with the singular kernel $w$ given by equation (2.6).

## 4.2. Numerical illustrations for $r_i = r_e$

### 4.2.1. Oscillatory solutions

We desire to exhibit the oscillatory nature of the solution $\phi$ to the minimization problem (3.2), (3.3) and (3.4). The values of the parameters are: $r_i = r_e = 1$, $\lambda = 1$, $\theta = \frac{\pi}{4}$, $\eta = 1/3000$, $L = 15$. We consider the ion boundary condition

$$ f_i^{\infty}(v) = \frac{n_i^{\text{ref}}}{\sqrt{2\pi v_{\text{th}}}} (v - v_c)^2 e^{-\frac{(v-v_d)^2}{2v_{\text{th}}^2}} \mathbf{1}_{\{v>v_c\}}(v), \quad (4.2) $$
<figure><img src="image_12.png" /><figcaption>FIGURE 10. Plot of the electric potential, linear-scale on top and log-scale on bottom, for three different values of the reference ionic density at infinity. One observes that the oscillations for high value of $n_i^{\text{ref}}$ vanish for low value of $n_i^{\text{ref}}$.</figcaption></figure>

with $v_{\text{th}} = 0.5$, $v_c = 1.0$, $v_d = 3/2$. Provided the condition $\alpha > 5\pi^2$ is satisfied, our theory predicts that at least oscillatory solutions exist. In our case, a rough estimation $\alpha \approx \int_0^{+\infty} f_i^{\infty}(v)dv - \int_0^{+\infty} \frac{f_i^{\infty}(v)}{v^2}dv$ shows that $\alpha$ is proportional to $n_i^{\text{ref}}$. Therefore choosing $n_i^{\text{ref}}$ large enough is a simple way to enforce $\alpha > 5\pi^2$. In the numerical illustrations reported in Figures 7, 8 and 9, we set $n_i^{\text{ref}} = 100$. For the finite element approximation we choose $N_x = 1000$. For the gradient algorithm we choose $\varepsilon = 1 \times 10^{-2}$, $\rho = 10^{-4}$. The gradient algorithm is stopped when the $L^\infty$ norm of the gradient of $L_\varepsilon$ is smaller than $\delta = 10^{-15}$. We plot in Figures 7–9 the computed electrostatic potential, its gyroaverage, and the ionic and electronic densities.
<figure><img src="image_13.png" /><figcaption>FIGURE 11. Electric potential $\phi$ (linear-scale on top and log-scale on bottom) for four values of the mass ratio $\eta$. One observes that the oscillations disappear for small value of $\eta$.</figcaption></figure>

### 4.2.2. Transition from oscillatory solutions to non oscillatory solutions

Here we recompute the problem of the previous section, except that $n_i^{\text{ref}}$ in (4.2) is set to three different values

$$n_i^{\text{ref}} = 100, 40 \text{ and } 10.$$

We plot the three profiles of the electric potential $\phi$ (in linear-scale and log-scale) in Figure 10. The oscillations for $n_i^{\text{ref}} = 100$ are mitigated with a larger wavelength for $n_i^{\text{ref}} = 40$, and eventually vanish for $n_i^{\text{ref}} = 10$. We believe that this behavior is in accordance with the results of the theoretical section.

## 4.3. Numerical illustrations for $r_i \neq r_e$

We provide some numerical illustrations when the Larmor radii $r_e$ and $r_i$ are different and scale as $r_e = \sqrt{\eta} r_i$ where $\eta < 1$ is the mass ratio between electrons and ions. In this case, we have modified the quadrature formulas used to calculate the potentials $N_i$ and $N_e$. Indeed it is better to guarantee that the equilibrium at $x \approx -\infty$ is perfectly respected in order to reduce numerical oscillations and to accelerate the numerical convergence. Such
<figure><img src="image_14.png" /><figcaption>FIGURE 12. Ionic and electronic densities for different values of the mass ratio $\eta = \frac{1}{160}$ (top left), $\eta = \frac{1}{320}$ (top right), $\eta = \frac{1}{640}$ (bottom left), $\eta = \frac{1}{1280}$ (bottom right). One observes that the charge separation close to $x = 0$ tends to increase for smaller values of the mass ratio.</figcaption></figure>

modified quadrature techniques are standard and are not described.

The values of the parameters are: $r_i = 1$, $\lambda = 1$, $\theta = \frac{\pi}{4}$, $L = 15$. We consider the ion boundary condition (4.2) with $v_{th} = 0.5$, $v_c = 1.0$, $v_d = 3/2$ and $n_i^{\text{ref}} = 100$. We solve the minimization problem (4.1) for four different values

$$ \eta = \frac{1}{160}, \frac{1}{320}, \frac{1}{640} \text{ and } \frac{1}{1280}. $$

For the finite element approximation we choose $N_x = 1000$. For the gradient algorithm we choose $\epsilon = 1 \times 10^{-2}$ and we normalize the gradient of $L_\epsilon$ with respect to the first iterate. We stop the gradient algorithm when the relative $L^\infty$ norm of the gradient is smaller than $\delta = 5 \times 10^{-4}$. We plot in Figures 11 and 12 the computed electrostatic potential, and the ionic and electronic densities for different values of the mass ratio. We essentially observe that the smaller the mass ratio, the weaker the oscillations on the electric potential and the larger the charge separation approaching $x = 0$.

We believe that it is in accordance with the mathematical structure of the Poisson equation (2.36) for which decreasing the mass ratio tends to make the electronic density dominant with respect to the ionic density. At the limit, the effect of the convolution operator vanishes so potential oscillations are less prone to manifest.
# APPENDIX A. PROPERTIES OF THE GYROAVERAGE OPERATOR

We state the main mathematical properties of the gyroaverage operator, which is the convolution with the kernel $w$ defined by (2.6).

**Lemma A.1 (Properties of the gyroaverage operator).**

One has:

(i) The function $w$ belong to $L^p(\mathbb{R})$ for all $p \in [1, 2[$.

(ii) For all $q, r \in [1, +\infty]$ and $p \in [1, 2[$ such that $\frac{1}{p} + \frac{1}{q} = 1 + \frac{1}{r}$, one has $\|w \star f\|_{L^r} \le \|w\|_{L^p} \|f\|_{L^q}$. In particular, by taking $p = 1$, one gets $\|w \star f\|_{L^q} \le \|w\|_{L^1} \|f\|_{L^q} = \|f\|_{L^q}$.

(iii) Define $q$ and $r$ as above. If $f \in L_{\text{loc}}^q(\mathbb{R})$, then $w \star f \in L_{\text{loc}}^r(\mathbb{R})$.

(iv) Let $p \in [1, +\infty[$ and let $s \in \mathbb{R}$. Suppose that $f \in \dot{W}^{s,p}(\mathbb{R})$. Then $w \star f \in \dot{W}^{s+\frac{1}{2},p}$.

*Proof. Proof of point (i):* The function $w$ is even. Then, it is enough to study its integrability near $x = 1$. Let $x \ge 0$, one has $0 \le \pi w(x) = \frac{1}{\sqrt{1-x^2}} \le \frac{1}{\sqrt{1-x}} = |1-x|^{-1/2}$. This last function belongs to $L^p([0, 1])$ for all $p \in [1, 2[$.

*Proof of point (ii):* This is the Young inequality for the convolution [31].

*Proof of point (iii):* Let $K := [a, b]$ be a measurable compact interval of $\mathbb{R}$ and let $f \in L_{\text{loc}}^q(\mathbb{R})$. The support of the function $w \star (f.1_K)$ is contained in $[a - 1, a + 1]$ and coincide with $w \star f$ on the set $[a + 1, b - 1]$ as a consequence of supp($w$) = $[-1, 1]$. By applying the result given by the previous point to the two functions $w$ and $f.1_K$, one gets that $w \star f \in L^r([a + 1, b - 1])$. This implies that $w \star f \in L_{\text{loc}}^r(\mathbb{R})$.

*Proof of point (iv):* As a consequence of the standard theory on fractional Sobolev spaces [17], it is enough to prove that the Fourier transform of $w$ satisfies $|\hat{w}(\xi)| \le \frac{C}{1+|\xi|^{1/2}}$ for all $\xi \in \mathbb{R}$. Using the parity of the function $w$ and a change of variable, one gets

$$|\hat{w}(\xi)| = \left| \int_{-1}^{1} \frac{e^{i\xi x}}{\sqrt{1-x^2}} dx \right| = 2 \left| \int_{0}^{1} \frac{e^{i\xi x}}{\sqrt{(1-x)(1+x)}} dx \right| \le 2 \left| \int_{0}^{1} \frac{e^{i\xi x}}{\sqrt{x}} dx \right| = \frac{2}{|\xi|^{1/2}} \left| \int_{0}^{\xi} \frac{e^{iy}}{\sqrt{y}} dy \right|$$

Since $\left| \int_{0}^{\xi} \frac{e^{iy}}{\sqrt{y}} dy \right| \le \min(c, |\xi|^{1/2})$ holds for some constant $c > 0$, one gets $|\hat{w}(\xi)| \le \frac{C}{1+|\xi|^{1/2}}$ for all $\xi \in \mathbb{R}$. The rate of decrease of $W$ at infinity yields the gain of regularity of the claim. $\square$

# APPENDIX B. BASIC PROPERTIES OF THE MINIMIZER

## B.1. Proof of Proposition 3.1

Equations (2.36) correspond to (3.11)–(3.13) (see below), so that we only need to prove the existence and uniqueness of a solution to the minimization problem (3.2)–(3.4). For this we first check that the functional $J$ is well-defined on the vector space

$$W_0 := \{ \psi \in \dot{H}^1(\mathbb{R}) \cap L^2(\mathbb{R}^+) : \psi(x) = \psi(0) \text{ for } x < 0 \}, \quad (B.1)$$

where the homogeneous Sobolev space $\dot{H}^1$ is the set of functions that admit a weak derivative in $L^2$. Up to the redefinition of the physical space for the $x$ variable, this space is identical to the one in (2.37). One has that

$$\forall \psi \in W_0, \quad \int_0^\infty |\psi'(x)|^2 dx < +\infty.$$

On the other hand, as a consequence of (3.10),

$$\forall \psi \in W_0, \quad \int_0^\infty F(w \star \psi)(x) dx \le \frac{\beta}{2} \int_0^\infty |w \star \psi(x)|^2 dx \le \frac{\beta}{2} \int_{-1}^\infty |\psi(x)|^2 dx < +\infty,$$
where for the second inequality one uses Proposition 3.2-(ii) in Appendix A. It is straightforward that the functional $\mathcal{J}$ defined at (3.3) is strictly convex on $W_0$ as a consequence of the strict convexity of the function $F$ given by Hypothesis (3.7). Recall that the affine space $V$ (3.4) inside which the minimization problem is defined is an affine subspace of $W_0$.

To obtain existence and uniqueness of this constrained minimization problem, it remains to prove that the functional $\mathcal{J}$ is coercive on $H^1(\mathbb{R}^+)$. Note that one has separately

$$ \mathcal{J}(\phi) \geq \frac{1}{2} \int_{0}^{\infty} \left| \frac{d}{dx} \phi(x) \right|^2 dx \quad \text{and} \quad \mathcal{J}(\phi) \geq \frac{\alpha}{2} \int_{0}^{\infty} |w \star \phi(x)|^2 dx \quad (B.2) $$

so one needs to obtain a good lower bound on the second term. One has $w \star \phi(x) = \phi(x) + R(x)$ where

$$ R(x) = \int_{-1}^{1} w(y)(\phi(x-y) - \phi(x))dy = \int_{y=-1}^{1} \int_{z=x}^{x-y} w(y)\phi'(z)dzdy. $$

One has directly (Fubini Theorem then Cauchy-Schwarz inequality)

$$ |R(x)| \leq \int_{y=-1}^{1} \int_{z=x-1}^{x+1} w(y)|\phi'(z)|dzdy = \int_{z=x-1}^{x+1} |\phi'(z)|dz \leq \sqrt{2} \left( \int_{z=x-1}^{x+1} |\phi'(z)|^2 dz \right)^{\frac{1}{2}}. $$

One has the basic inequality $|w \star \phi(x)|^2 = |\phi(x) + R(x)|^2 \geq \frac{1}{2} |\phi(x)|^2 - |R(x)|^2$. Therefore one can write

$$ \int_{0}^{\infty} |w \star \phi(x)|^2 dx \geq \frac{1}{2} \int_{0}^{\infty} |\phi(x)|^2 dx - 2 \int_{x=0}^{\infty} \int_{z=x-1}^{x+1} |\phi'(z)|^2 dzdx $$

that is

$$ \int_{0}^{\infty} |w \star \phi(x)|^2 dx \geq \frac{1}{2} \int_{0}^{\infty} |\phi(x)|^2 dx - 4 \int_{0}^{\infty} \left| \frac{d}{dx} \phi(x) \right|^2 dx. \quad (B.3) $$

Plugging (B.3) in (B.2), one gets

$$ \begin{aligned} (6\alpha + 1)\mathcal{J}(\phi) &\geq 3\alpha \int_{0}^{\infty} \left|\frac{d}{dx}\phi(x)\right|^2 dx + \frac{\alpha}{4} \int_{0}^{\infty} |\phi(x)|^2 dx - 2\alpha \int_{0}^{\infty} \left|\frac{d}{dx}\phi(x)\right|^2 dx \\ &\geq \frac{\alpha}{4} \int_{0}^{\infty} |\phi(x)|^2 dx + \alpha \int_{0}^{\infty} \left|\frac{d}{dx}\phi(x)\right|^2 dx. \end{aligned} $$

This is the coercivity estimate, so the proof of the existence is ended. The functional being strictly convex, the minimum is unique.

Equations (3.11)-(3.13) then follow by rewriting Proposition 2.14 in the setting of Section 3. Indeed for the mathematical analysis we have replaced the physical variable and potential (denoted here $\tilde{x} < 0$ and $\tilde{\phi}$ to avoid confusion) by $x = -\tilde{x} > 0$ and $\phi(x) = -\tilde{\phi}(\tilde{x})$. With the simplification (3.1) we have $w = w_e = w_i$, and the function $F$ involved in the model energy (3.3) corresponds to $F(\psi) = (N_e - N_i)(-\psi)$, in particular (2.35) gives

$$ (\tilde{n}_e(\tilde{\phi}) - \tilde{n}_i(\tilde{\phi}))(\tilde{x}) = -w \star (1_{\mathbb{R}^+} F'(w \star \phi))(\tilde{x}). $$

Together with $\tilde{\phi}''(\tilde{x}) = -\phi''(x)$ and $\tilde{\phi}'(0^{-}) = \phi'(0^{+})$, these relations allow to rewrite (2.36) into (3.11)-(3.13). In passing we note that the variational principle (2.38) reads here

$$ \int_{0}^{\infty} \phi'(x) h'(x) dx + \int_{0}^{\infty} F'(w \star \phi(x)) (w \star h(x)) dx - \sigma w \star h(0) = 0 \quad (B.4) $$

for any test function $h \in W_0$ defined by (B.1).
## B.2. Proofs for Proposition 3.2

### B.2.1. Proof of Proposition 3.2-(i)

Using the variational principle (B.4) with the test function $h = \phi$, one ends up with

$$ \int_{0}^{\infty} |\phi'(x)|^2 dx + \int_{0}^{\infty} F'(w \star \phi(x)) (w \star \phi(x)) dx = \sigma \int_{-1}^{1} w(x) \phi(x) dx \quad (B.5) $$

One now uses assumption (3.9) to write for all $x \in \mathbb{R}^+$,

$$ \min \{\alpha w \star \phi(x); \beta w \star \phi(x)\} \le F'(w \star \phi(x)) \le \max \{\alpha w \star \phi(x); \beta w \star \phi(x)\}. $$

The combination of these two facts with $\alpha \le \beta$ gives

$$ \alpha \int_{0}^{\infty} |w \star \phi(x)|^2 dx \le \int_{0}^{\infty} |\phi'(x)|^2 dx - \sigma \int_{-1}^{1} w(x) \phi(x) dx \le \beta \int_{0}^{\infty} |w \star \phi(x)|^2 dx. $$

□

### B.2.2. Proof of Proposition 3.2-(ii)

Take $\nu > 0$. One defines

$$ \phi_{\nu}(x) := \begin{cases} \nu e^{-\sqrt{\beta}(x+1)} & \text{if } x > 0, \\ \nu e^{-\sqrt{\beta}} & \text{otherwise.} \end{cases} $$

A direct computation gives

$$ \frac{1}{2} \int_{0}^{\infty} |\phi'_{\nu}(x)|^2 dx = \frac{\nu^2 \sqrt{\beta}}{4} e^{-2\sqrt{\beta}}. \quad (B.6) $$

On the other hand, using (B.10) one has

$$ \int_{0}^{\infty} F(w \star \phi_{\nu}(x)) dx \le \frac{\beta}{2} \int_{0}^{\infty} |w \star \phi_{\nu}(x)|^2 dx \quad (B.7) $$

Since the function $x \mapsto \phi_{\nu}(x)$ is non-increasing and since $\int w = 1$, one has that $0 \le w \star \phi_{\nu}(x) \le \phi_{\nu}(x-1)$. Therefore (B.7) implies $\int_{0}^{\infty} F(w \star \phi_{\nu}(x)) dx \le \frac{\beta}{2} \int_{0}^{\infty} |\phi_{\nu}(x-1)|^2 dx$. One infers by a direct computation that

$$ \int_{0}^{\infty} F(w \star \phi_{\nu}(x)) dx \le \frac{\nu^2 \sqrt{\beta}}{4} e^{-2\sqrt{\beta}}(1 + 2\beta). \quad (B.8) $$

Gathering (B.6) and (B.8) gives

$$ \mathcal{J}(\phi_{\nu}) \le \frac{\nu^2 \sqrt{\beta}}{2} e^{2\sqrt{\beta}}(1 + \beta) \quad (B.9) $$

One sets the value of $\nu$ such that $w \star \phi_{\nu}(0) = \psi_{\text{wall}}$ so that $\phi_{\nu} \in V$. This means

$$ \psi_{\text{wall}} = \nu \int_{-1}^{0} w(x) e^{-\sqrt{\beta}} dx + \nu \int_{0}^{1} w(x) e^{-\sqrt{\beta}(x+1)} dx $$

Thus,

$$ \nu = \frac{\psi_{\text{wall}}}{\int_{-1}^{0} w(x) e^{-\sqrt{\beta}} dx + \int_{0}^{1} w(x) e^{-\sqrt{\beta}(x+1)} dx} \le \frac{2}{e^{-\sqrt{\beta}} + e^{-2\sqrt{\beta}}} $$

where for the last inequality one used $e^{-\sqrt{\beta}(x+1)} \ge e^{-2\sqrt{\beta}}$ when $x \in [0; 1]$. Therefore (B.9) gives $\mathcal{J}(\phi_{\nu}) \le 2^{2}\sqrt{\beta}(1 + \beta) \frac{e^{-2\sqrt{\beta}}}{(e^{-\sqrt{\beta}} + e^{-2\sqrt{\beta}})^2}$. Since one has $e^{-2\sqrt{\beta}} \le (e^{-\sqrt{\beta}} + e^{-2\sqrt{\beta}})^2$, one concludes $\mathcal{J}(\phi) \le \mathcal{J}(\phi_{\nu}) \le 2 \psi_{\text{wall}}^2 \sqrt{\beta} (1 + \beta)$. □
### B.2.3. Proof of Proposition 3.2-(iii)

One starts by integrating equation (3.11) and get

$$ \int_{0}^{\infty} -\phi''(x) \, dx + \int_{0}^{\infty} w \star (1_{\mathbb{R}^{+}} F'(w \star \phi))(x) \, dx = \sigma \int_{0}^{\infty} w(x) \, dx. $$

Therefore, $\phi'(0^+) + \int_0^\infty w \star (1_{\mathbb{R}^+} F'(w \star \phi))(x) \, dx = \frac{\sigma}{2}$. One plugs this in (3.12) to get the desired equation. □

### B.2.4. Proof of Proposition 3.2-(iv)

From (B.5), one has

$$ \int_{0}^{\infty} |\phi'(x)|^2 dx + \int_{0}^{\infty} F'(w \star \phi(x)) (w \star \phi(x)) dx = \sigma \psi_{\text{wall}}. $$

The left hand side is positive and $\psi_{\text{wall}}$ is also positive. So $\sigma$ is positive.

As a consequence of equation (3.11) and since $w \to +\infty$ at +1, it is enough to prove that

$$ x \mapsto w \star (1_{\mathbb{R}^+} F'(w \star \phi))(x) \in L^{\infty}([0, 1]). \quad (B.10) $$

For $x \in [0, 1]$, one has by triangular inequality $|w \star (1_{\mathbb{R}^+} F'(w \star \phi))(x)| \le \int_{-1}^2 w(x-y) |F'(w \star \phi)(y)| dy$. One uses Hypothesis (3.9) and gets

$$ \begin{aligned} |w \star (1_{\mathbb{R}^+} F'(w \star \phi))(x)| & \le \beta \int_{-1}^{2} w(x-y) |w \star \phi(y)| dy \\ & \le \beta \int_{-1}^{2} \int_{-2}^{3} w(x-y) w(y-z) |\phi(z)| dz dy \quad (B.11) \\ & \le \beta \int_{-2}^{3} |\phi(z)| (w \star w)(x+z) dz, \end{aligned} $$

where one uses the Fubini Theorem for the last inequality. Using now the Cauchy-Schwarz inequality and then using Property (ii) of Lemma A.1 in appendix with parameters $p = q = 4/3$ and $r = 2$, one gets

$$ |w \star (1_{\mathbb{R}^+} F'(w \star \phi))(x)| \le \|\phi\|_{L^2} \|w \star w\|_{L^2} \le \|\phi\|_{L^2} \|w\|_{L^{4/3}}^2. $$

Then one has proved (B.10) and this complete the proof. □

### B.2.5. Proof of Proposition 3.2-(v)

As a consequence of Proposition 3.2-(v), the function $\phi$ cannot change sign infinitely many times as $x \to 1$ with $x < 1$. In other words, there exists $\varepsilon$ with $0 < \varepsilon < 1$ such that

$$ \forall x \in ]1 - \varepsilon, 1[, \quad \phi(x) \neq 0. \quad (B.12) $$

Suppose now by the absurd that $\phi(x) = 0$ for all $x > 1$. One is going to prove that in this case $\phi''(3 - \varepsilon) \neq 0$, which yields a contradiction. For that purpose, one computes

$$ w \star (1_{\mathbb{R}^+} F'(w \star \phi))(3 - \varepsilon) = \int_{2-\varepsilon}^{4-\varepsilon} w(3 - \varepsilon - y) F'\left(\int_{y-1}^{y+1} w(y-z) \phi(z) dz\right) dy $$

Using the fact that $\phi(x) = 0$ for all $x > 1$ with $F'(0) = 0$, consequence of (3.9), the equality above becomes

$$ w \star (1_{\mathbb{R}^+} F'(w \star \phi))(3 - \varepsilon) = \int_{2-\varepsilon}^{2} w(3 - \varepsilon - y) F'\left(\int_{y-1}^{1} w(y-z) \phi(z) dz\right) dy \quad (B.13) $$
Nevertheless, as stated by (B.12), the function $\phi$ does not change sign on $]-1; -1 + \varepsilon[$. Suppose for instance and without loss of generality that $\phi(x) > 0$ for all $x \in ]1 - \varepsilon; 1[$. Then $\int_{y-1}^{1} w(y-z) \phi(z) dz > 0$ for all $y \in ]2 - \varepsilon; 2[$. Since $F' > 0$ on $\mathbb{R}_{*}^{+}$ as stated by (3.9) and since $w$ is positive on its support, one gets $\int_{2-\varepsilon}^{2} w(3 - \varepsilon - y) F' \left( \int_{y-1}^{1} w(y-z) \phi(z) dz \right) dy > 0$. Plugging this back into (B.13) and using equation (3.11) leads to $\phi''(3 - \varepsilon) > 0$. Since $\varepsilon < 1$, this is in contradiction with $\phi(x) = 0$ for all $x > 1$. $\square$

*Acknowledgements.* Mehdi Badsi and Bruno Després acknowledge the support of ANR project MUFFIN under contract ANR-19-CE46-0004. Ludovic Godard-Cadillac acknowledges the support of ANR project ODA under contract ANR-18-CE40-0020.

## REFERENCES

[1] A. Arsenev, Existence in the large of a weak solution of Vlasov's system of equations. *Mat. Mat. Fiz.* **15** (1975) 136-147.

[2] M. Badsi, *Etude mathématique et simulation numérique de modèles de gaines bi-cinétiques*. Ph.D. thesis, Université Pierre et Marie Curie, Paris (2016).

[3] M. Badsi, Linear electron stability for a bi-kinetic sheath model. *J. Math. Anal. Appl.* **453** (2017) 954-872.

[4] M. Badsi, M. Campos Pinto and B. Després, A minimization formulation of a bikinetic sheath. *Kinet. Relat. Models* **9** (2016) 621-656.

[5] M. Badsi, M. Merhemberger and L. Navoret, Numerical stability of plasma sheath. *ESAIM: Proc.* **64** (2018) 17-36.

[6] C. Bardos and P. Degond, Global existence for the Vlasov-Poisson equation in 3 space variables with small initial data. *Ann. Inst. Henri Poincaré* **2** (1985) 101-118.

[7] N. Ben Abdallah, Weak solutions of the initial-boundary value problem for the Vlasov-Poisson system. *M2AS* **17** (1994) 451-476.

[8] D. Bohm, The characteristics of electrical discharges in magnetic fields. McGraw Hill, New York (1949).

[9] M. Bostan, Existence and uniqueness of the mild solution for the 1D Vlasov-Poisson initial-boundary value problem. *SIAM J. Math. Anal.* **37** (2005) 156-188.

[10] M. Bostan, Gyrokinetic models for strongly magnetized plasmas with general magnetic shape. *Discrete Contin. Dyn. Syst* **5** (2012) 257-269.

[11] A. Bottino and E. Sonnendrucker, Monte Carlo particle-in-cell methods for the simulations of the Vlasov-Maxwell gyrokinetic equations. *J. Plasma Phys.* **81** (2015) 435810501.

[12] Y. Brenier, Convergence of the Vlasov-Poisson system to the incompressible Euler equations. *Comm. Partial Diff. Equ.* **25** (2000) 737-754.

[13] A.J. Brizard and J.W. Burby, Gauge-free electromagnetic gyrokinetic theory. *Phys. Lett.* **18** (2019) 2172-2175.

[14] R. Chalise and R. Khanal, A kinetic trajectory simulation model for magnetized plasma sheath. *Plasma Phys. Control. Fusion* **54** (2012) 095006.

[15] F.F. Chen, Introduction to Plasma Physics and controlled fusion. Springer (1984).

[16] R. Chodura, Plasma-wall transition in an oblique magnetic field. AIP Publishing (1982).

[17] E. Di Nezza, G. Palatucci and E. Valdinoci, Hitchhiker's guide to the fractional Sobolev spaces. *Bull. Sci. math.* **136** (2012) 521-573.

[18] M. Feldman, S.Y. Ha and M. Slemrod, A geometric level-set formulation of a plasma sheat interface. *Arch. Rat. Mech. Anal.* **178** (2005) 81-123.

[19] E. Frénod and E. Sonnendrücker, The finite Larmor radius approximation. *SIAM J. Math. Anal.* **32** (2001) 1227-1247.

[20] J.P. Friedberg, Plasma Physics and Fusion Energy. Cambridge University Press (2007).

[21] X. Garvet, P. Mantica, C. Angioni and E. Asp, Physics of transport in tokamaks. *Plasma Phys. Control. Fusion* **46** (2004) B557-B574.

[22] P. Ghendrih, M. Hauray and A. Nouri, Derivation of a gyrokinetic model. Existence and uniqueness of specific stationary solutions. *Kinet. Relat. Models* **2** (2009) 707-725.

[23] R. Glowinski, *Numerical methods for nonlinear variational problems*. Springer-Verlag (1984).

[24] V. Grandgirard and Y. Sarazin, Gyrokinetic simulations of magnetic fusion plasmas. In: *Numerical models for fusion* **39** (2013) 91-176.

[25] E. Grenier, Oscillatory limits with varying spectrum. In: *Congrès National de Mathématiques Appliquées et Industrielles, Volume 35 of ESAIM Proc.* EDP Sciences, Les Ulis (2011) 46-58.

[26] Y. Guo, Regularity for the Vlasov equations in a half space. *Indiana Univ. Math. J.* **43** (1994) 255-320.

[27] Y. Guo, The dynamics of a plane diode. *SIAM J. Math. Anal.* **35** (2004) 1617-1635.

[28] J.-H. Hwang and J. Schaeffer, Uniqueness for weak solutions of a one-dimensional boundary value problem for the Vlasov-Poisson system. *J. Diff. Equ.* **244** (2008) 2665-2691.
[29] H.-J. Hwang and J.J.L. Velazquez, Global existence for the Vlasov-Poisson system in bounded domains. *Arch. Rat. Mech. Anal.* **195** (2010) 763–796.

[30] J. Kromes, Dielectric response and thermal fluctuations in gyrokinetic plasma. *Phys. Fluids B Plasma Phys.* **5** (1993) 1066.

[31] E. Lieb and M. Loss, Analysis (2nd version). Graduate Studies in Mathematics, Vol. 14. American Mathematical Society (2001).

[32] N. Mandell, A. Hakim, G. Hamett and M. Francisquez, Electromagnetic full-f gyrokinetics in the tokamak edge with discontinuous galerkin methods. *J. Plasma Phys.* **86** (2020) 905860109.

[33] G. Manfredi and D. Coulette, Kinetic simulations of the Chodura and Debye sheaths for magnetic fields with grazing incidence. *Plasma Phys. Control. Fusion* **58** (2016) 025008.

[34] E. Miot, The gyrokinetic limit for the Vlasov-Poisson system with a point charge. *Nonlinearity* **32** (2019) 654–677.

[35] J. Moritz, E. Faudot, S. Devaux and S. Heuraux, Plasma sheath properties in a magnetic field parallel to the wall. *Phys. Plasmas* (2016).

[36] P.J. Morrison, The Maxwell-Vlasov equations as a continuous Hamiltonian system. *Phys. Lett. A* **80** (1980) 383–386.

[37] F.W. Olver, D.W. Lozier, R.F. Boisvert and C.W. Clark, NIST Handbook of mathematical functions. Cambridge University Press (2010).

[38] F. Poupaud, Boundary value problems for the stationary Vlasov-Maxwell system. *Forum Math.* **4** (1992) 499–527.

[39] P.A. Raviart and C. Greengard, A boundary-value problem for the stationary Vlasov-Poisson equations: the plane diode. *Commun. Pure Appl. Math.* **43** (1990) 473–507.

[40] K.U. Riemann, The Bohm criterion and sheath formation. *J. Phys. D: Appl. Phys.* **24** (1991) 493.

[41] J. Schaeffer, Global existence of smooth solutions to the Vlasov-Poisson system in three dimensions. *Commun. Partial Diff. Equ.* **16** (2009) 1313–1335.

[42] E. Shi, A. Hakim and G.W. Hammet, A gyrokinetic 1D scrape-off layer model of an elm heat pulse. *Phys. Plasmas* **22** (2015) 022504.

[43] L. St-Raymond, Control of large velocities in the two-dimensional gyrokinetic approximation. *J. Math. Pure. Appl.* **81** (2002) 379–399.

[44] P. Stangeby, The plasma boundary of magnetic fusion devices. Institute of Physics Publishing (2000).

[45] P. Stangeby, The Chodura sheath for angles of a few degrees between the magnetic field and the surface of divertor targets and limiters. *Nucl. Fusion* **52** (2012) 083012.

[46] A. Weinstein and P.J. Morrison, Comments on the Maxwell-Vlasov equations as a continuous Hamiltonian System. *Phys. Lett. A* **86** (1981) 383–386.

Subscribe to Open (S2O)
A fair and sustainable open access model

<figure><img src="image_15.png" /></figure>

This journal is currently published in open access under a **Subscribe-to-Open model (S2O)**. S2O is a transformative model that aims to move subscription journals to open access. Open access is the free, immediate, online availability of research articles combined with the rights to use these articles fully in the digital environment. We are thankful to our subscribers and sponsors for making it possible to publish this journal in open access, free of charge for authors.

**Please help to maintain this journal in open access!**

Check that your library subscribes to the journal, or make a personal donation to the S2O programme, by contacting
**subscribers@edpsciences.org**

More information, including a list of sponsors and a financial transparency report, available at: https://www.edpsciences.org/en/maths-s2o-programme