Derivation of a gyrokinetic model. Existence and uniqueness of specific stationary solutions

A finite Larmor radius approximation is derived from the classical Vlasov equation, in the limit of large (and uniform) external magnetic field. We also provide an heuristic derivation of the electroneutrality equation in the finite Larmor radius setting. Existence and uniqueness of a solution is proven in the stationary frame for solutions depending only on the direction parallel to the magnetic field and factorizing in the velocity variables.

Introduction. The ITER project is a challenge to the growing need of new sources of energy. It aims at producing energy by nuclear fusion. Nuclear reactions take place in a tokamak, where a high temperature plasma is confined. So far, confined plasmas are performed with relatively short energy confinement times due to microscale instabilities that generate turbulent transport [III99]. It is observed that the characteristic frequencies of these instabilities is several orders of magnitude smaller than the ion Larmor gyration frequency governed by the strong magnetic field. Studies of nuclear fusion in tokamaks are in full expansion, both experimentally and theoretically. Kinetic models are appropriate for studying the core of the plasma since the collisions have a very weak effect in these hot and low density plasmas. Physicists use gyrokinetic models and especially the finite Larmor radius approximation to model the core of the plasma [GIVT09]. Taking into account the fast Larmor gyration of the charged particles that characterizes magnetic confinement, these models allow one to average over that fast gyration and reduce the 6D kinetic problem to a 5D gyrokinetic one. That property is especially interesting for numerical simulation. In this paper the finite Larmor radius approximation is derived from the Vlasov equation, in the limit of large uniform magnetic field and with an external electric field. Because of the homogenization on the fast Larmor gyration, the limit equation (1.10) is written in 5D gyro-coordinates (x g , v , |v ⊥ |) defined in (1.11). These coordinates are the position of the so-called particle guiding center, x g in the 3D space together with the parallel velocity v and the amplitude of transverse velocity |v ⊥ | that is proportional to the magnetic moment, an adiabatic invariant of the particle motion in the strong magnetic field limit (given a constant magnetic field) . Its mathematical structure is a combination of the Vlasov equation in the direction parallel to the magnetic field and of the Euler equation in the perpendicular direction, where the original fields are replaced by the corresponding gyro-average fields. To close the system, physicists use the electroneutrality equation n e = Zn i where n e stands for the electron density and n i the ion density, each ion having a charge Z. In the following Z = 1 will be considered with no loss of generality. It can be shown that the electroneutrality equation is in fact the the Poisson equation solved on scales that are significantly larger than the Debye length. Since the latter governs the Laplacian term of the Poisson equation, the electroneutrality equation is an appropriate approximation on scales of the order or larger than the Debye length. Moreover, taking into account the difference between the density of particles and that of guiding centers leads one to introduce a polarization correction due to the non uniform distribution of particles on gyro-circles. In the second part of this paper, the electroneutrality equation (1.25) is carefully written. The coupling of the 5D gyrokinetic Vlasov equation (1.10) to the electroneutrality equation (1.25)is the model used for instance in the GYSELA code, a project that aims at modeling the turbulent transport in fusion plasmas [GSG + 06b, GSG + 06a]. While in the GYSELA code there is the possibility to use a collision operator we concentrate here on the actual Vlasov equation with no collisions.
A difficulty raised by the electroneutrality equation taken as such is the lack of an explicit regularization term for the electric potential. Consequently, the regularity of the latter is not sufficient to ensure a mathematical solution of the Vlasov transport equation. The analysis of the well-posedness of limit model based on the electroneutrality equation thus seems impossible today. For that reason we restrict ourselves to the 1D Vlasov equation (without gyro-kinetic effects), coupled to the electroneutrality condition. That model is a particular case of the 3D model, provided the solutions do not depend on the direction perpendicular to the magnetic field. Even for that simple kinetic model, not much is known about solutions. The difficulty lies in the fact that the force term is proportional to the derivative of the density in the field direction. Indeed, while the Vlasov equation ensures that all L p -norms of the distribution function can be bounded, there is no control on the norms of its derivatives. In this paper, we prove the existence and uniqueness of a steady state solution in a slab geometry, therefore between two boundary conditions, provided one fulfills some conditions, in particular that they are no particles trapped between the two boundaries.

Derivation of the finite Larmor radius approximation
The ion distribution in low density plasmas submitted to a magnetic field is well described by the Vlasov equation. The latter is valid provided one can neglect the two-particle distribution function altogether [Nic83] so that the evolution of the standard distribution function is governed by the Liouville equation for the conservation of particles where f is the phase space density depending on time, position (in the domain D of the plasma) and velocity (in R 3 ) of the ions, and where Ze and m i are the ion charge and mass respectively. As stated in the introduction, we shall consider Z = 1 in the following with no loss of generality. A priori, the electric and magnetic fields are governed by the Maxwell equations, but a scale analysis allows one to approximate, and thus simplify, these laws. This will be made clear in the following. For a strong external magnetic field, the charged particles exhibit a fast rotation motion around the magnetic field lines. The frequency of that gyration, the Larmor frequency, is several orders of magnitude larger than the observed frequency range of turbulence. Furthermore, in the quasineutral limit this frequency is larger than all other frequencies of the particle dynamics, thus providing the means for an efficient scale separation. It is noteworthy that this is the basis of magnetic confinement that is presently realized in devices such as the tokamaks. In this framework, one separates two parts in the particle motion, on the one hand the slow motion of the center of the Larmor gyration, and, on the other hand, the fast gyro-motion. The phase space reduction achieved by only considering the slow motion, thus ignoring the fast gyro-motion, leads one to the 5D gyrokinetic model. These are usually derived by physicists either by averaging the single particle dynamics ( [Cat78] and [CTB81]) or by a Lie-transform perturbative approach ( [Hah88], [BH07] and [GLV08]). Mathematically, this is achieved by looking at the limit of equation (1.1) when the modulus of the external magnetic field |B ext | tends to infinity. However, different models can be obtained, depending on the way the other control parameters vary when |B ext | → ∞. For a magnetized plasma, two relevant scales characterize the limit regimes that have been discussed. These are the Debye length λ D and the thermal Larmor radius of the ions ρ th .
The Debye length weights the Laplacian in the Poisson equation and thus defines the transition from the non-neutral description of the plasma required on the sub-Debye scales from the quasineutral plasma description for scales larger than the Debye length. The thermal Larmor radius is the radius of the fast Larmor gyration for ions with average speed v th . The averaging over the small time scales governed by the ion Larmor frequency then tends to translate into an averaging over the ion Larmor scale. The finite Larmor radius terms are therefore introduced to take into account this rather weak cut-off effect. The radius of the electrons gyration is generally much smaller for comparable ion and electron temperatures and will be neglected. The two scales discussed above are defined by where the electron thermal energy T e is introduced rather than that of the ions since it is more appropriate to characterize the large electron mobility and therefore strong response to the electric field. Let us introduce the characteristic scale of the system L, for instance introduced by the boundary conditions. The large magnetic field limit then corresponds to the vanishing ion Larmor radius limit, namely the gyro-center approximation where ρ th /L → 0. Within this framework, E. Grenier uses pseudodifferential calculus to prove [Gre97] the convergence of a 2D fluid model towards an incompressible fluid model when ε goes to zero with ρ th /L = ε and (λ D /L) 2 = ε. For the same scaling and cold initial distributions, i.e. the approximation of a Dirac measure at velocity zero for the distribution function, Y. Brenier proved in [Bre00] with modulated energy technics the convergence of solutions of the Vlasov-Poisson system towards dissipative solutions (introduced by P.-L. Lions) of the Euler equation, for weelprepared initial data. F. Golse and L. St-Raymond used the same technics in [GSR03] to prove the convergence of a Vlasov-Poisson model on the torus towards the 2D1/2 Euler equation. This is done for a vanishing ε with ρ th /L = ε and using another scaling for the Debye length, namely λ D /L = ε. This property is restricted to a torus size such that there is no resonant oscillation. This particular case appears to be relevant for electrons in a region close to the tokamak boundary.
Derivations have also been performed in the gyrocenter approximation, ρ th /L ⊥ small but finite. Here L ⊥ is a carcteristic length in the perpendicular direction, which is usually choosen smaller thant the one in the parallel direction in finite Larmor radius approximation. In [FS00a] (see also [FS00b] for a short version) E. Frénod and E. Sonnendrücker studied the convergence of a linear Vlasov equation (external magnetic field), in this finite Larmor radius limit (large B but finite Larmor radius), using two scale convergence methods. M. Bostan studied the 2D strong magnetic field limit, using Hilbert expansion technics, in a setting of finite Debye length and Larmor radius, λ D /(εL) = ρ th /(εL) = 1, where ε is small but finite [Bos09] and obtained strong convergence results for regular initial conditions. Here, we further examine the linear case studied by Frénod and Sonnendrücker, and obtain a simpler and more complete derivation. The finite Larmor radius approximation is correct for fusion plasmas, including the core and most of the edge plasma. However, the ion Larmor radius is much larger than the Debye length throughout the plasma. An interesting derivation should be done in the framework of present gyrokinetic calculations, namely ρ th /(εL) = 1, λ D /(εL) → 0, yet nothing has been done in this very demanding limit (see the beginning of 1.2 for more details).
The resolution of equation (1.1) is known in the case where the fields E and B are external and C 1 or at least BV . In the BV -case, one can use the DiPerna Lions theory of transport equation [DL89b] with its latest developments ( [Amb04]). In the case where the electric potential is given by self-interaction with the usual Poisson law, and the magnetic field B is still considered as external, results of existence and uniqueness obtained for the Vlasov equation without magnetic field may be used with the appropriate modifications. We refer to [LP91], [Hor93]. In the case where the self-induced magnetic field is not negligible any more, we refer to the work [DL89a] of DiPerna and Lions about the Vlasov-Maxwell equation. In the case of an electric field given by an electroneutrality equation, nothing is known from the mathematical point of view.

The finite Larmor radius gyro-kinetic approximation. Rigorous approach.
For the sake of simplicity, we neglect the variation of the external magnetic field B ext . In this cylindrical approximation, the curvature of the magnetic field is not considered, and one neglects the exploration of the magnetic field variation by the particles during the fast cyclotron motion. With respect to the relevant scales this effect if of the order of the Larmor radius ρ th divided by the characteristic radius of the magnetic field curvature, namely the tokamak major radius R. Although ρ th /R is a small parameter, this approximation is too strong with respect to many aspects of the tokamak physics, in particular since curvature effects are considered by physicists as the cause of a large class of micro-instabilities. The strong impact of such a small parameter can be traced back to a symmetry breaking of the system since the curvature governs the leading term that produces a charge separation. The present simplification must be regarded as a first step that is only valid when one addresses issues that do not lead one to symmetry breaking. Moreover, in order to avoid problems with boundary conditions, which are complex to handle, a periodic setting is used. In other words, we work on the torus T 3 in space, and R 3 in speed, with a constant magnetic field B = |B|(0, 0, 1).
Next, we define dimensionless variables with characteristic time τ , parallel length L , perpendicular length L ⊥ , velocity v th and electric field E 0 .
In this new set of dimensionless variables the Vlasov equation (1.1) may be rewritten as: where the subscript ⊥ (resp. ) stands for the projection on the perpendicular (resp. parallel) direction, and the superscript ⊥ stands for the projection on the plane orthogonal to the field line and the rotation by −π/2 around the field line direction: The time scale τ can conveniently be defined to reduce the number of control parameters in Eq.(1.3) by setting: In a similar fashion, the normalization of the electric field can be used to define the third control parameter, hence: The two remaining control parameters are then L /L ⊥ for the third term in Eq.(1.3) and L /ρ th for the last term. It will be assumed here that both parameters L ⊥ and ρ th exhibit the same asymptotic behavior L ⊥ ∝ ρ th . Let us define the perpendicular scale as L ⊥ = ρ th 2π/(k ⊥ ρ th ). In this expression, k ⊥ is the typical wave vector of the cross-field fluctuations. For a sufficiently small Larmor radius, one can then assume that the turbulence follows the so-called gyro-Bohm scaling such that k ⊥ ρ th is a constant [LEHT02]. In particular it does not depend on the magnitude of the magnetic field. One then finds the following relationship between the two control parameters that can be grasped by the following relation whenever one drops the proportionality constant In practice, one could also define the normalization scale as L ⊥ = ρ th so that the the two control parameters would be exactly the same. In this discussion, we have considered a transverse scale related to the fluctuation properties that allows one to characterize one of the terms in the Vlasov equation. However, other control parameters must be considered to account for the boundary conditions. The socalled ρ * parameter used in magnetic fusion is such a parameter since ρ * = ρ th /a where a is the plasma minor radius, the scale related to the boundary conditions while ρ th is the fluctuation scale. Relating the parameter ε to ρ * introduces aspect ratios that stem from the periodic boundary conditions.
The ratio a/L depends of the safety factor and the tokamak aspect ratio in the actual tokamak geometry but represents the ratio of the sizes of the domain in the radial and parallel direction with the present setting. The parameter ε can also be expressed in terms of the slow and fast characteristic times of the particle motion. The slow time τ introduced in the normalization of Eq.(1.3), is the characteristic time to explore the tokamak geometry along the parallel direction and the fast time is due to the gyration motion, hence: In this last step one thus finds that the ion Larmor frequency Ω = eB/m i is of order (τ ε) −1 . In other words, one thus assumes that the ion Larmor frequency is much larger than the parallel connection frequency τ −1 . That is a valid assumption in all regions of a tokamak. Indeed, for ITER conditions, the ion Larmor frequency is of the order of 2 10 8 Hz and the connection time τ ranges from 10 −3 s to 10 −4 s. These parameters mainly depend on the magnitude of the magnetic field and on the size of the device so that ε characterizes a given fusion device, ε ≈ 10 −5 for ITER parameters.
Dropping the primes, we obtain the following rescaled version of the Vlasov equation (1.1), When ε → 0, the largest terms are those proportional to the factor 1/ε. Retaining the latter, we obtain a transport equation associated to the following system of ODE in the plane transverse to the magnetic field,ẋ As B is homogeneous, the trajectories are circles of center x g = x ⊥ + v ⊥ ⊥ , and radius |v ⊥ | covered with frequency ε −1 . The global motion is the sum of this very quick motion of gyration and a slower and more complicated motion (with velocity of order one). In the limit of large |B|, particles are assumed to be evenly distributed on the gyro-circles, and their motion is the sum of a drift in the perpendicular direction, and a classical acceleration in the parallel direction. Heuristically, the electric drift v E may be obtained from the Newton law in normalized variables, Given the assumption that the particles have a fast motion of gyration, one can integrate the previous equation over a period of gyration, at lowest order, v E is given by the averaged velocity such that: (1.7) In the latter expression one readily recognizes the usual form of the electric drift velocity v E = E × B/B 2 , however where the electric field is averaged over a period E . One can show that this average of the field E corresponds to an average over a circle of radius |v ⊥ | in the perpendicular plane. Provided the only dependence on the gyrophase stems from the particle motion, this average translates into a Bessel operator defined as: where ρ L is the ion Larmor radius and e iϕc = (cos ϕ c , sin ϕ c , 0). We also introduce the operatorJ 0 ρ L which stands for a position-velocity version of this average, where e = (0, 0, 1).
Given the average on the gyrophase ϕ c , the motion is reduced to a 5D space. The gyration radius ρ L = |v ⊥ | (given the chosen normalization) is related to the magnetic moment µ = m i ρ 2 L /(2eB), a quantity that is an adiabatic invariant, i.e. a constant of motion in the large B limit. One can show that this invariant is the conjugate variable of the gyrophase. Here, as B is homogeneous, ρ L = |v ⊥ | will remain constant in the limit |B ext | → +∞. We thus obtain a 4D + 1D model, i.e. a 5D model with no dynamics in the variable derived from the magnetic moment. This reduction of the phase space is very interesting for numerical simulations. The following theorem states this reduction precisely.
Theorem 1.1 Let us assume that f 0 ε is uniformly bounded in L q , q > 1, that it weakly converges toward f 0 ∈ L q , and that E is a gradient and belongs to L p t,x where p −1 + q −1 = 1. For each ε > 0, let f ε be a solution of (1.4) with the initial condition f 0 ε . Then, up to a subsequence (ε n ), (f ε (t, x g − v ⊥ , v)) weakly converges tof in L q , wheref only depends on (t, x g , v , |v ⊥ | = ρ L ) and is solution of Remark 1.2 The operatorJ 0 is important to perform the adaptation of the 6D initial condition to the 5D limit model. Indeed, the very fast Larmor gyration creates an initial layer that instantaneously adapts the initial condition to the limit model.
Proof of Theorem 1.1. The phase space in position-velocity coordinates is not well adapted to perform the fast gyration averaging. To handle it more easily, it is convenient to change the system of coordinates and consider the gyro-coordinates defined by (1.11) x g is the position of the gyro-center and ρ L = |v ⊥ | is the ion Larmor radius. To express the gradient in x and v in this new system of coordinates, one can conveniently remark that: for any smooth function h, with the functionh defined byh(x g , v g ) = h(x, v). Then ∇ x = ∇ xg and ∇ v = ∇ vg − ∇ ⊥ xg , where ∇ ⊥ xg = (∂ xg,2 , −∂ xg,1 , 0). Note that ∇ ⊥ stands for the gradient vector rotated by −π/2 and not π/2. The gradients are taken at the corresponding points. For instance, the first equality reads ∇ x h(x, v) = ∇ xgh (x g , v g ). The anti-symmetry of ⊥, a · b ⊥ = −a ⊥ · b, has been used (and we shall make a wide use of it in the sequel). Given these relations, (1.4) can be modified leading one to an equation satisfied by the functionf ε (t, with initial conditionf 0 . Here the subscript ⊥ stands for the perpendicular components to the magnetic field, for instance E ⊥ = (E 1 , E 2 , 0). Barred quantities are functions of the gyro-coordinates. Since (1.6) is a conservative transport equation, i.e. it may be written as ∂ t f + div(. . . ) = 0, the L q -norms of f ε and then off ε are conserved. Thus, f ε L ∞ t (L q x,v ) is uniformly bounded. Then, up to the extraction of a subsequence, we may assume thatf ε weakly converges towards somef ∈ L ∞ t (L q x,v ). Upon multiplying equation (1.13) by ε, in the limit ε → 0, we obtain v ⊥ g · ∇ vgf = 0 , All the other terms from Eq.(1.13) are bounded in the sense of distributions and their product with ε therefore vanishes in the limit ε → 0. This reduced form of the Vlasov equation implies that the only dependence off on v g,⊥ is on |v g,⊥ |. Hence with no dependence on the gyrophase. Let us now consider equation (1.13) forf ε , when integrated against a smooth test fonction φ with support in time avoiding t = 0 (we will handle the initial conditions later). Let us further assume that the dependence on v g,⊥ is restricted to a dependence on |v g,⊥ |. For such a function, one readily finds that 1 ε f ε v ⊥ g ∇φdx g dv g = 0 6 hal-00477454, version 1 -29 Apr 2010 for symmetry reasons. As a consequence, the projection off ε takes the following form keeping in mind that E is calculated at the point x g − v ⊥ g . We have also used the relation a ⊥ · b = −a · b ⊥ as well as the fact that equation (1.13) may be written in a conservative form because ( if J is the matrix since E is a gradient. At this stage, we can take the limit and obtain thatf also satisfies equation (1.14). This is not yet a proper equation in the sense of distributions, even though f depends on |v g,⊥ |. Indeed, the electric field still exhibits the dependence on the particle position E ≡ E(t, x g − v ⊥ g ). To obtain an equation only depending on |v g,⊥ |, we use polar coordinates for v g,⊥ , i.e. v g = (ρ L cos ϕ c , ρ L sin ϕ c , v ) and Fubini's theorem to integrate the previous integral first in ϕ c , then in the other variables. This leads one to the following equation: where the the operator J 0 ρ L is defined in (1.8), and the term R is equal to ) · e iϕc dϕ c , with the previously introduced notation e iϕc = (cos ϕ c , sin ϕ c , 0). It can be shown that R is null since it is the circulation of the electric field E along a gyrocircle. The latter vanishes in the electrostatic case since E is a gradient. With R = 0, (1.15) is identical to the equation (1.10) written distribution wise for 2πρ Lf . As there is no dynamics in the ρ L direction, the factor ρ L is only a multiplicative constant. One can then introduce Φ/ √ ρ 2 L +η as test function, and then obtain thatfρ L / √ ρ 2 L +η satisfies the equation. Letting η → 0 in that linear equation, we obtain the appropriate result for the equation governing the evolution off .
Regarding the initial conditions we use a similar projection technique with a test function Φ depending on |v g, |, but not vanishing at t = 0. We then obtain equation (1.14), with the right-hand side replaced by (1.16) The asymptotic limit then leads one to: (1.17) Changing coordinate to recover the position-velocity coordinates in the ϕ c integral, we exactly obtain the initial conditionJ ρ L f 0 that is expected forf .
This proof is valid for external electric and magnetic fields. The case where the equation (1.13), invariant in the parallel direction, is coupled to the Poisson equation given a Debye length of the same order as the Larmor radius has been treated by Frénod and Sonnendrücker in [FS01]. It corresponds to a Debye length of order √ ε. Technically, we may handle that case with our technic to obtain weak solution as in Arsenev's work [Ars75]. The main point is to obtain L p estimates on the density n i = f dv. They can be obtained by classical estimates using upper bounds for the kinetic energy. The case where the Debye length is taken much smaller than the Larmor radius will be of greater interest, but the difficult problem is there to average the strong oscillations appearing at the scale of the Debye length. We refer to [Gre95], [Gre96] and [CG00] for more details on quasi-neutral plasma without magnetic field, and to [GSR03] for results in the guiding-center approximation with ρ L ∝ ε, λ D ∝ ε and ε → 0.

The electroneutrality equation. Heuristic approach.
In this section we address the self-consistent problem when linking the electric field to the charge distribution. The relevant equation is the Maxwell-Gauss equation that relates the divergence of the electric field to to the local charge governed by the particle density of charged particles. The latter must be determined using the ion distribution function for the guiding centers that is solution of the gyrokinetic Vlasov equation. A similar treatment for the electrons must be done. We will follow heuristic arguments together with assumptions that are not justified rigorously. However, this heuristic derivation bares some interest. It is an alternative to the physicists' presentation of that equation based on the Fourier transform. Furthermore, it allows one to recover the electroneutrality equation (1.25) used in the GYSELA code to close the gyrokinetic equation (1.10).
For quasineutral plasmas, the electric field response to any charge separation governs a restoring force. Should the charge separation extend on a scale larger than the Debye length, the restoring force would be to strong to allow any significant charge build-up. The plasma can thus be considered to be everywhere with near zero charge, hence quasineutral. As a consequence, the density of negative charges en e (n e being the density of electrons) equals the density of positive charges en i (n i being the density of ions). In the electrostatic limit, neglecting the time dependence on the vector potential, one can recover this physics based argument as an asymptotic limit of the Poisson equation for the electric potential.
The right hand side is dimensionless and so is eΦ/T e . The dimensionless control parameter on the left hand side operator thus appears as the square of the ratio of the Debye scale divided by the characteristic scale of the charge separation. In the limit λ 2 D → 0, one readily recovers the quasineutrality equation, namely: n e = n i . Let us consider the ions density. Given the Vlasov equation and its dependence on the electric field, it is likely that the ion distribution function (and thus the ion density) is an implicit function of the electrostatic potential. The same applies for the electron density. Then, in the limit λ 2 D → 0, there is no analytical dependence on the electric potential. The latter then becomes a Lagrangian multiplier associated to the quasineutrality equation. On such issues, we acknowledge the work of E. Grenier [Gre95], [Gre96] (the only work, at our knowledge, on that subject where the limiting model is kinetic and not fluid), Y. Brenier [BG94], S. Cordier [CG00] and N. Masmoudi [Mas01].
However, if we assume that the system is close to equilibrium, i.e. that the departure from a constant electric potential is small, then the electroneutrality equation n e = n i provides an explicit dependence on a mean-field electric potential. For an electron-ion plasma, the mass ratio is such that the electrons are far more mobile than the ions. One can then assume that their response to an electrostatic perturbation is adiabatic on magnetic field lines or surfaces. In that case, assuming that the equilibrium density of electrons n e,0 and that of ions n i,0 are equal to n 0 , one may write where Φ is the average of Φ on a closed magnetic field line or surface. The expansion on the right hand side holds if eΦ << T e . In the very simple geometry that is considered here (with B = (0, 0, B)), Φ is the average in the x direction: Φ = Φ(x) dx . The approximation of adiabatic electrons thus reintroduces an explicit dependence on the electrostatic potential. For the perturbation of the ion distribution, the first difficulty is to obtain the distribution of ions in physical space fromf written in gyro-coordinates. If the distribution in physical space is assumed to be constant on gyrocircles as shown in the last section for the limit model, the following formula is obtained, (1.19) Taking the integral in v leads to where J 0 ρ L only acts on the x variable (since ρ L = |v ⊥ | with our conventions). The quasineutrality equation, n e = n i , then leads one to (1.21) However, the assumption that the distribution f is constant on the gyrocircles is unrealistic whenever the electric potential is not constant. In order to express the inhomogeneity of the density on gyrocircles, we add a perturbation tof . It can be expressed as an adiabatic perturbation on the gyrocircles Φ is the average of Φ over the gyrocircle of a given particle, and f i (v) is the equilibrium distribution of ions. Note that this adiabatic perturbation depends on v ⊥ through the choice of the gyrocircle used in determining the averageΦ.
where e iϕc = (cos ϕ c , sin ϕ c , 0) and h i (ρ L ) = 2πρ L f i (ρ L , v ) dv . When the equilibrium distribution is a maxwellian, f i (v) = 1/(T i π) 3/2 e −|v| 2 /Ti , so that h i (ρ L ) = 2ρ L /T i e −ρ 2 L /Ti . Adding the adiabatic perturbation of the gyrocircles, we finally obtain the following electroneutrality equation, (1.24) Multiplying by T e /e, then yields is the convolution in the perpendicular plane with the radial function H T i (r) defined by (1.27)

See Appendix A for a detailed proof.
Since ρ L is a parameter in the equation of motion (1.10), the equilibrium distribution and the initial perturbation remain concentrated on a unique value of ρ L at all time provided it is the case for the initial conditions. A first step to solve the system might be to start under this assumption of single value of ρ L . This would likely be the most difficult step, since the general case is a superposition of such cases. The gyrokinetic Vlasov equation (1.10) is unchanged when considering a single value of ρ L , however the electroneutrality equation (1.25) may be simplified and becomes: wheren i (t, x) = 2πρ L f (t, x, v ) dv is the density of ions in gyro-coordinates. Still, after this further simplification, the gyrokinetic Vlasov equation (1.10) coupled to the electroneutrality equation (1.28) remains is a very difficult mathematical problem. The main difficulties are the following.
• The lack of regularity in the parallel direction. The potential Φ has the regularity of f in the parallel direction. This is to be compared to the Vlasov-Poisson case where D 2 Φ (∆Φ in the Poisson equation given above) has the regularity of f . We may overcome this problem by adding some viscosity in that direction, in other words by adding a term −λ∂ 2 v f in equation (1.10).
• A less important lack of regularity lies in the perpendicular direction. In fact, only the gyro-average of Φ appears in (1.10). Moreover, the term Φ is more regular than Φ, so that by (1.28) Φ has the regularity of J 0 ρ L (ρ). Hence J 0 ρ L (Φ) has the regularity of (J 0 ρ L ) 2 (ρ). The operator (J 0 ρ L ) 2 sends L 2 into H 1 . That is "almost" enough, in the sense that we could reach a sufficient regularity provided (J 0 ρ L ) 2 were compact from L 2 in H 1 .
In view of these difficulties, we focus in the next section on the lack of regularity in the parallel direction, and consider time-independent solutions depending only on x and v, of the form f (x , v )f ⊥ (|v ⊥ |).
Then the term f i has no incidence in the evolution equation and can thus be factorized in (1.10). Taking T e = 1 for the sake of simplicity, taking again into account the equilibrium density n 0 , not always constant in that section and therefore denoted by n, and replacing n i by ρ, the electroneutrality equation (1.25) then reads: Note that g ′ is the derivative of g with respect to z. Furthermore, we consider the problem in a slab geometry z ∈ [−1, 1], and therefore require given f ± as boundary conditions, Lemma 2.1 Given n positive, then any solution f to (2.1)-(2.2), such that ρ n is non decreasing, must satisfy: Hence V 2 (s) + 2 ρ n (Z(s)) = v 2 + 2 ρ n (z).

Remark 2.4
The assumption (H2) ensures that ρ is a perturbation of the equilibrium n, so that with assumption (H1) the force-field will always be oriented rightward (no possibilities of trapped particles).
The assumption (H3) ensures that there are two beams of ions (one coming from the right and the other from the left), and that all the ions coming from one side reach the other side (no turn-back). The (H4) assumption is more technical and enforces that the derivative of the density is bounded.