A variational sheath model for gyrokinetic vlasov-poisson equations

. 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


Introduction
The mathematical description and numerical simulation of a plasma interacting with a metallic surface is an active topic of research in plasma physics [8,16,33,40,44,45].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/inafewlineswhich 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 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 wellestablished [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][4][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 gyrocenter positions in 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  < 0, their velocity component along the magnetic field is denoted  ∈ R and their magnetic moments is denoted  ∈ 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   and an ionic Larmor radius   .These radii are arbitrary in the final model.
For the mathematical analysis, it is convenient to choose equal Larmor radius   =   in order to minimize technical difficulties.This hypothesis, the same unique radius for ions and electrons, is similar of other monokinetic 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.4is 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.

Gyroaverage operator
In magnetized plasma physics, charged particles follow helicoidal trajectories.A key notion is the gyroaveraging on the Larmor radius.Let B ̸ = 0 be a constant non zero given magnetic field 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 = ( 1 ,  2 ,  3 ) involves the Larmor vector   () :=  (cos  v 1 + sin  v 2 ) ∈ R 3 with  ∈ [0, 2), and where v 1 = (− sin , cos , 0) and v 2 = (0, 0, 1) are two vectors orthonormal to B. With these notations, a trajectory with radius  > 0, parallel direction b = B |B| , parallel velocity  ‖ ̸ = 0 and cyclotron frequency   ̸ = 0 verifies the Newton's second law of motion where the convolution kernel  is defined almost everywhere by The convolution kernel  belongs to   (R) for 1 ≤  < 2 and its total mass is normalized, ∫︀ R ()d = 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   (R) to  + 1 2 (R) for any  ∈ R.

Gyrokinetic equations
Following [11,13,20], a plasma made of one species of ions and electrons subject to a given magnetic field B is modeled with gyrokinetic Vlasov-Poisson equations.Positions are denoted as  ∈ R. Velocities in the  direction are denoted as  ∈ R. Magnetic moments are denoted as  = ±|v ⊥ | ∈ R.
The unknowns are the electrostatic potential  :  ∈ R ↦ → () ∈ R, the particle density for ions   : (, , ) ∈ R × R × R ↦ →   (, , ) ∈ R + and the particle density for electrons They satisfy the gyrokinetic Vlasov-Poisson system (2.7)-(2.9) ) (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  , () =  , √︀ || > 0 where  , > 0. We refer to p. 5 of Reference [11] for the definition of the Debye length  > 0 in function of physical quantities in the context of fusion plasmas.
In particular small Debye length arise from the physical scaling One has  ≈ 6.4−5 in fusion plasmas ( [20], p. 126).The mass ratio between ions and electrons is a natural small parameter denoted as  =   ≈ 1/2000.The coefficient sin 2  in equation (2.9) comes from the fact that it is the projection along the first coordinate of the operator ∇ * ⊥ ∇ ⊥ where ∇ ⊥ is the gradient orthogonal to the magnetic field B. Since one also has a sin  in the gyroaveraging (2.5), it is immediate to observe that sin  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 ) By integration of the Vlasov equations (2.7) and (2.8), we find that the parallel current densities   and   are constant.

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 (, ) ∈  = (−∞, 0) × 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  = 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   (0,  < 0) = 0 and   (0,  < 0) = 0. (2.15) We will extend the functions inside the wall with (2.16) These closure relations express that the spatial domain (0, +∞) 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.

Electronic phase diagram
Starting from a given smooth potential 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     + 1   ′ ()    = 0 in .The oscillations of the gyroaveraged potential  are caused by similar oscillations of the physical potential .
The level curves are defined by 1 2  2 − 1  () = constant.Where the potential is regular enough ( ∈  1 loc for example), the trajectories correspond the characteristics curves which are the solutions of the characteristics equations  ′ =  and  ′ = 1   ′ ().The characteristics curves correspond to connected components of the level curves.Therefore, if two points (, ) ∈  and (, ) ∈  are connected by a characteristic curve then For  = 0, one gets the condition Making the assumption that the minimum of  is reached at the wall, that is then the inequality (2.18) delimits 3 open regions in , as described in Figure 2 }︁ .
The superior region  1 is connected to  = −∞, so particles coming from the infinity with a positive velocity travel through  1 and are absorbed at the wall.The internal region  2 contains closed loops because the function  may have some oscillations (we will prove this fact).Particles in this region do not reach the wall.The inferior region  3 is connected to the wall at  = 0. However no electron with negative velocity is emitted at the wall (2.15), so this region is empty of electrons. originates from the pure electrostatic setting (that is b = 0) for which the physics of the sheath is dominant [4].The case (0) ̸ = min R −  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.
We consider where  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   defined in (2.19) is a weak solution in  to the equation The "physical" electronic density  phys  = ∫︀   d, see (2.14), can be written in different ways Another possibility that will be used later is where the Gauss error function is erf The constant electronic current is equal to its limit Figure 3. Four different regions of the ion phase space.The dotted lines are general characteristics.The bold line is the characteristic 1 2  2 + () =  + .

Ionic phase diagram
Next we construct the phase diagram for ions for the potential The function  may also have oscillations as in Figure 1.With the method of characteristics  ′ =  and  ′ = − ′ (), if two points (, ) ∈  and (, ) ∈  are connected with a characteristic curve then 1  2  2 +() = is a barrier of the potential.For simplicity as illustrated in Figure 1, we assume that the maximum value  + > 0 is positive (consistently with the condition at infinity lim −∞  = 0) and that  + is unique.Indeed a point (, ) such that 1 2  2 + () <  + cannot be connected to ( + , ).This is the reason of the decomposition of the domain in four open regions Remark 2.3.Note that the assumption that  + is unique simplifies the phase diagram 3. The case of several points where the maximum value  + > 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     −  ′ ()    = 0 in  based on the diagram of Figure 3, we consider a function at infinity for incoming ions  ∞  ().Using the approach [4], we assume that  ∞  is continuous with compact support in { > √︀ 2 + + }.More precisely where 0 <  is a small number and  is taken large enough so that √︀ 2 + +  < .We consider the continuous function where the last equality is obtained using the fact that  ∞  vanishes for  < √︀ 2 + .The ionic density at infinity is The electronic current is constant in space, so it is equal to its value at infinity

Determination of the sheath potential 𝜓 wall and the reference density 𝑛 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 Considering equations (2.21) and (2.22), we denote Observe that the condition of neutral charge at −∞ and the condition of neutral flux take the form (2.26) In this system, the right hand side is given because it depends on the function  ∞  which is independent of  and .The unknown is the pair ( wall ,  ref ) with the sign conditions The function  is strictly monotone, with (0) = 1 and (+∞) = +∞.Therefore there exists a positive solution, necessarily unique, if and only if 1 Once  wall is determined, the reference density is given.

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 , 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.

Potential representation of the densities
From the general definition (2.13) and (2.19), one gets Let us define the function By definition one has that   (()) =   ⋆ ( ′  (  ⋆ )) ().In (2.28), the function   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 By definition the function of the claim vanishes at the origin   (0) = 0.So it is equal to (2.28).
We use a similar approach for the ions.Let us consider Define the function By definition   (0) = 0.
Proof.One has Plugging this identity in (2.30) yields the claim.

The energy
For mathematical correctness of the material presented below, it is necessary to have continuous extension of the function   () for  <  wall and continuous extension of the function   () for  >  + .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   ,   ∈  0 (R),   ,   are piecewise  1 (R) and   ,   ∈  2 (−, ) for some  > 0.
Definition 2.8.The energy is defined formally by )︂ d. (2.32) An affine functional space naturally adapted to this energy is where we remind that Ḣ1 (R) is the set of measurable functions  such that  ′ ∈  2 (R).The definition of  accounts for the first condition of (2.16) inside the wall and for the boundary condition (2.25) at the wall.Functions  ∈  are such that lim +∞ () = 0. We observe that  is a complete space for the norm of ) and functions in  are constant on R + , then one has then Remark 2.10.The need for  1 integrability is due to a mismatch between the   and   functionals at −∞, it will be relaxed in the mathematical analysis of Section 3.

A stability condition at infinity
Here we study a condition such that  is convex with respect to small perturbation at  ≈ −∞ of the null function.The condition is written as a convex condition around  = 0 Then the condition (2.34) of convexity at infinity holds.
Proof.Using the physical densities, one gets easily

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 (  ,   ,  ∞  , . . .).Note however that a general convexity result will be established in Section 3 but for   =   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 us define now the non linear operators where 1 R − () = 1 for  ≤ 0 and = 0 for  > 0 is the characteristic function of R − .Using the physical densities (2.20) and (2.24), we observe that so that ñ () and ñ () coincide with the gyroaveraged densities   =   ⋆  phys

𝑒
and   =   ⋆  phys  for  < −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  ′ ∈  0 (R − ), and the Euler-Lagrange relations where  wall < 0 is the (electronic) gyroaveraged wall potential defined by Proposition 2.5, and  ∈ R is a Lagrange multiplier associated to the boundary constraint.
Remark 2.15.For  < −  away from the wall,   vanishes and the first equation in (2.36) coincides with the initial gyrokinetic Poisson equation (2.13).For the same reason, the perturbed densities ñ , ñ vanish for  > max(  ,   ), 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 (, ) = () +  (  ⋆ (0) −  wall ) for which a convenient functional space is W0 =  0 ∩  1 (R − ) with Note that Ṽ ⊂ W0 and Ṽ ̸ = W0 .Observing that the energy functional  is continuously defined from W0 into R (and not only from Ṽ as in Lem.2.9), we see that the Lagrangian  is also continuously defined from W0 into R. Let us calculate the variation with respect to test functions ℎ ∈ W0 where we note that for ℎ(0) ̸ = 0 this relation is only justified provided  ′ ∈  0 (R − ), then for  ∈ {, }, using the symmetry of the kernels   , This yields (with the same reservation as above if ℎ(0 Taking an arbitrary ℎ ∈ W0 with ℎ(0) = 0 gives the first equation of (2.36).This shows that  ′′ ∈  1 loc (R − ), so that  ′ ∈  0 (R − ) and (2.39) holds indeed with ℎ(0) ̸ = 0.The second equation of (2.36) is then the remaining coefficient in front of ℎ(0).Finally the boundary condition in (2.36) follows from the variations in .

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 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   <   .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  > 0 with a positive wall potential  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: with an energy functional now defined as The functional space  is redefined as Here the parameters are the function  : R → R whose properties will be discussed below, the gyroaveraging kernel  defined in (2.6) and the wall potential  wall , which is now assumed positive for mathematical simplicity, as specified above.We also remind that the space Ḣ1 is defined as the set of measurable functions  which derivative is in  2 , and that functions  in the affine space  are such that lim +∞  = 0.The tangent space is Note that if the function  is replaced by a Dirac mass, we recover a standard nonlinear Poisson problem with a Dirichlet constraint at  = 0.With our simplification (3.1), the function  : R → R corresponds to the difference   −   from the modeling Section 2.8.1.Thus it should gather some main properties of the density potentials   and   .For the mathematical analysis we assume that it satisfies the following hypotheses: •  (0) = 0; (3.6) The condition (3.6) expresses the physical condition   (0) =   (0) = 0.The condition  ′ (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  ′ lays between the curves  =  and  =  as illustrated in Figure 4. We also draw in Figure 4 the fact that the function  ′′ may have a mild singularity at  wall since it is the case for the physical model, see equations (2.28)-(2.29).Since the singularity there scales like (3.10)

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.

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 ).The minimization problem given by (3.2), (3.3) and (3.4) admits a unique solution.This solution  ∈  satisfies  ′ ∈  0 (R + ) and it is characterized by the Euler-Lagrange equations where and  ⋆ (0) =  wall > 0. (3.13) The real number  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   =   = 1.Concerning the function , minimizer of  on  , we are able to extract the following elementary properties.Proposition 3.2 (Basic properties of the minimizer ).Let  ∈  be the solution to the minimizing problem (3.2), (3.3) and (3.4).Then the solution  satisfy the following properties.

Local monotonicity of 𝑤 ⋆ 𝜑 near the wall
The property of local monotonicity near the wall guarantees that, at least in the vicinity of the wall, the function  ⋆  takes values below  wall .This property is important because the physical modeling of Section 2.5 give a priori zero information for values above  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  + > 0 also helps).

Oscillatory nature of the solution for large 𝛼
The analysis of the solution near the wall is rather complicated but we expect it to be concave decreasingat least on a neighborhood of the wall.The difficulty comes from the presence of the Lagrange multiplier  in equation (3.11), which is active only for  < 1 since the support of  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  is large.We remind that, by Proposition 3.2-(), there exists  0 ∈ [1; +∞) such that ( 0 ) ̸ = 0. Proposition 3.4 (Oscillatory nature of the solution).Assume that  defined at (3.10) satisfies  ≥ 5 2 .Take  0 ∈ [1; +∞) such that ( 0 ) ̸ = 0. Then there exists  1 ∈ [ 0 ;  0 + 10] such that ( 0 ) ( 1 ) < 0. Remark 3.5.From this proposition, we conclude that the solution  oscillates infinitely many times around the value  = 0. Indeed, by using this proposition recursively we get an increasing sequence (  ) ∈N such that for all index  we have (  ) ( +1 ) < 0. Proposition 3.2-() 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  1 is lower than  0 + 10.
This phenomenon generates new open questions.What is the limit value of ?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.
One also has ℎ is non-negative.The property (3.9) combined with (3.18) yields On the other hand, an integration by part gives where again  ≥ 0 on [ 0 ;  0 + 10].One now considers  very large so that  >  and obtains Moreover, as a consequence of (3.17) the function  is non-negative on [ 0 ;  0 + 10] but it is not identically 0 since ( 0 ) > 0. Thus, the right-hand side of (3.• Step 3: Here, we prove that  > 5 2 .The idea here is to consider a particular function ℎ for the formula (3.21) and to provide an estimate on its associated coefficient .
• Firstly one defines a function : R → R + (this notation has nothing to do with the previous function  in the modeling section) as follows By definition, one has that  ∈  2,∞ (R) is supported on [0; 2], is non-negative and is such that One computes ∫︀∫︀  ()d d by splitting the computation on these 6 sub-domains and using directly (3.24).First, one has ∫︁ ∫︁ () d d = 0 One gets the 4 remaining terms by a direct computation using the Fubini theorem, that is in order Invoking now the fact that  is symmetrical with respect to  = 1, one can do the very same reasoning with  ∈ [3/2; 2].Therefore Gathering now the estimates (3.28) and (3.29), one concludes where  ′′ was computed at (3.25).

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  = 0 yields the linear problem.It is also a particular case of the general situation, since it is sufficient to take  =  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 Lemma 3.6.One has () =  0 () where  0 is the first modified Bessel function.
Proof.The proof is by intersection of graphs, as in the Figure 6.
Now we study equation (3.32) for  >  * , and for possible complex solutions  ∈ C.
Lemma 3.8.One can rewrite equation (3.32) under the form where  is an analytic series (with a certain radius of convergence) with real coefficients.
Proof.One rewrites equation (3.32) as A fixed point method yields two solutions in the vicinity of  * , but in the complex plane, that is  =  + i with ,  ∈ R and  ̸ = 0.
In summary of this discussion is that the solutions of (3.32) are either real or complex.
- -If  =  * , it is also possible to construct solutions that change sign under the form -Since  ⋆  = () and  ⋆  ⋆  = () 2 , then the changes of sign hold also for the convolved functions  ⋆  and  ⋆  ⋆ .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].

A physical interpretation of the oscillatory solutions
The two methods explained in Sections 3.1.3and 3.1.4exhibit oscillatory solutions for  > 0 large enough, which means for a function  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 , 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.

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  < 0, with negative wall potentials  wall < 0.

Description of the numerical method
The numerical method is made of two steps.
-The first step consists in assembling the density potentials   and   and computing the reference density  ref > 0 and the Dirichlet boundary condition on the wall potential  wall < 0 by solving (2.26)-(2.27).We will consider a function  ∞  that is piecewise continuous, so that standard quadrature formulas for velocity integrals can be used to assemble the potential   .The other potential poses no additional difficulty.Next,  ref > 0 and  wall < 0 are computed by means of a Newton-Raphson method applied to the function As for the Dirichlet boundary condition at the wall, we use a standard penalty method [23].We consider  > 0 small enough and the following penalized minimization problem on  ℎ 0 : find  ℎ ∈  ℎ 0 such that where the Lagrangian   is defined for all  ℎ ∈  ℎ 0 by with  the energy functional from (2.32).For any  > 0 the minimization problem (4.1) has a unique solution   ℎ since  ℎ 0 is a finite dimensional subspace of  0 and   is coercive and strictly convex.It is standard that ) where  ℎ is the unique minimizer of  on  ℎ 0 such that   ⋆  ℎ (0) =  wall .To compute   ℎ , we use a standard gradient method which consists in computing recursively the sequence ( , ℎ ) ∈N ⊂  ℎ 0 given by where  > 0 is a given parameter.By the Riesz representation theorem, the gradient ∇  ( , ℎ ) ∈  ℎ 0 is the unique solution of the variational problem 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  given by equation (2.6).with  th = 0.5,   = 1.0,   = 3/2.Provided the condition  > 5 For the finite element approximation we choose   = 1000.For the gradient algorithm we choose  = 1 × 10 −2 ,  = 10 −4 .The gradient algorithm is stopped when the  ∞ norm of the gradient of   is smaller than  = 10 −15 .We plot in Figures 7-9 the computed electrostatic potential, its gyroaverage, and the ionic and electronic densities.where  < 1 is the mass ratio between electrons and ions.In this case, we have modified the quadrature formulas used to calculate the potentials   and   .Indeed it is better to guarantee that the equilibrium at  ≈ −∞ is perfectly respected in order to reduce numerical oscillations and to accelerate the numerical convergence.Such For the finite element approximation we choose   = 1000.For the gradient algorithm we choose  = 1 × 10 −2 and we normalize the gradient of   with respect to the first iterate.We stop the gradient algorithm when the relative  ∞ norm of the gradient is smaller than  = 5 × 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  = 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.
where for the second inequality one uses Proposition 3.2-() in Appendix A. It is straightforward that the functional  defined at (3.3) is strictly convex on  0 as a consequence of the strict convexity of the function  given by Hypothesis (3.7).Recall that the affine space  (3.4) inside which the minimization problem is defined is an affine subspace of  0 .
To obtain existence and uniqueness of this constrained minimization problem, it remains to prove that the functional  is coercive on  1 (R + ).Note that one has separately where one uses the Fubini Theorem for the last inequality.Using now the Cauchy-Schwarz inequality and then using Property () of Lemma A.1 in appendix with parameters  =  = 4/3 and  = 2, one gets Then one has proved (B.10) and this complete the proof.

Figure 4 .
Figure 4. Condition (3.9) states that the function  ′′ may show integrable singularities, as it is the case for the physical model in the vicinity of the wall at the value  wall .
23) is positive.Plugging this back into (3.21)gives lim →0  ( +  ℎ) −  ()  > 0, which eventually contradicts the minimality of the function .It already proves the claim for  large enough.It remains to show the lower bound on .

2 ; 2 )Figure 5 .
Figure 5.An illustration of the set   ⊆ R 2 , subdivided into  ,1 • • •  ,6 represented by the different shades of grey.This subdivision intervenes in the computation of the integral (3.26).

Figure 7 .
Figure 7. Plot of the electrostatic potential and gyroaveraged electrostatic potential.

Figure 8 .
Figure 8. Plot in log-scale of the electrostatic potential and gyroaveraged electrostatic potential.The drop at  ≈ −15 is due to machine precision, because the boundary condition was numerically set to  = 0 on the left boundary to emulate the condition at infinity.

Figure 9 .
Figure 9. Plot of the ionic density and electronic density.

Figure 10 .
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  ref  vanish for low value of  ref  .

Figure 11 .
Figure 11.Electric potential  (linear-scale on top and log-scale on bottom) for four values of the mass ratio .One observes that the oscillations disappear for small value of .

4. 3 .
Numerical illustrations for   ̸ =   We provide some numerical illustrations when the Larmor radii   and   are different and scale as   = √

Figure 12 .
Figure 12.Ionic and electronic densities for different values of the mass ratio:  = 1 160 (top left),  = 1 320 (top right),  = 1 640 (bottom left),  = 1 1280 (bottom right).One observes that the charge separation close to  = 0 tends to increase for smaller values of the mass ratio.