Numerical study of an anisotropic Vlasov equation arising in plasma physics

Goal of this paper is to investigate several numerical schemes for the resolution of two anisotropic Vlasov equations. These two toy-models arise from a kinetic description of a tokamak plasma confined by strong magnetic fields. The simplicity of our toy-models permits to better understand the features of each scheme, in particular to investigate their asymptotic-preserving properties, in the aim to choose then the most adequate numerical scheme for upcoming, more realistic simulations.


Introduction
The present paper addresses a new approach for an efficient numerical resolution of anisotropic transport models, which simplified are of the type   subject to appropriate boundary conditions (here periodic ones). The unknown f ǫ stands for the quantity (distribution function) which is advected along the given (or self-consistently computed) field u in the domain Ω := [0, L x ] × [0, L y ] and the small scaling parameter ǫ ≪ 1 indicates that we have to deal with very strong advection fields u or equivalently with the long-time asymptotics of f ǫ . Such anisotropic transport models arise very often in physics, as simplifications of more complex systems. In Section 2 we detail some examples coming from plasma physics, as the Vlasov equation for the ion dynamics in the gyrokinetic regime. There are however several other examples arising in physics and leading to a simplified transport equation as (1.1), for example when one studies the long-time asymptotics of the incompressible Euler 2D equations, (1.1) representing then the vorticity equation, which has to be coupled (via u) with a Poisson equation for the stream-function computation [19].
A numerical resolution of problems of the type (1.1) is rather challenging in the regime ǫ ≪ 1, due to the singularity of the mathematical problem as ǫ → 0. Certainly, the exact solution of the simple transport-case (1.1) is known for ǫ > 0, however not in general situations, when u is self-consistently computed via f ǫ and when other (notstiff) terms are present. These general situations require then an efficient numerical treatment of (1.1). From a physical point of view we can say that we have to cope with a multiscale problem, the parameter ǫ being the stiffness parameter. Standard schemes (explicit hyperbolic approaches) require very restrictive CFL-conditions (dependent on ǫ) in order to accurately account for the microscopic ǫ-scales. Very often in such situations people are impliciting the stiff term [8], in order to avoid these too restrictive CFL-conditions. This can work in some situations, for example when the grid is aligned with the anisotropy, and only for a certain range of ǫ-values. However in more general configurations, not-aligned grids and ǫ-values covering all the interval [0, 1], impliciting the stiff term is no more sufficient, as shall be seen in this paper. We propose thus in this work a new numerical procedure, based on Asymptotic-Preserving arguments, being able to solve (1.1) in an efficient manner, uniformly accurate and stable in ǫ, and this on a simple, Cartesian grid. Asymptotic-Preserving methods are efficient, as they are designed in order to mimic on the discrete level the asymptotic behavior of the singularly perturbed problem solutions (see [15,22] for a detailed introduction). This paper was initiated by the repetitive remarks/questions one of the authors got during conferences, meaning that impliciting the stiff term in (1.1) is enough to get an efficient AP-scheme, which behaves well even in the limit ǫ → 0. The aim of this paper is to prove the contrary, AP-schemes are more than impliciting the stiff term. In order to understand in detail the main features of the here proposed AP-scheme, we preferred to keep the investigated model as simple as possible, so that a detailed numerical analysis is possible, permitting to perceive the differences of our scheme when compared to standard (implicit) schemes. We hope that doing so, we are able to resolve some of the confusion that surround AP-schemes. However, even if the here presented results are based on a simplified model as (1.1), the same Asymptotic-Preserving approach can be used for more involved anisotropic transport problems, such as those presented in Section 2 and which shall be the objective of an upcoming work. The AP-procedure we propose here was employed in other contexts by the authors (elliptic [6,7], parabolic [20]). The present setting is more stimulating, as we have to cope with highly oscillating problems when ǫ ≪ 1 and no more dissipative ones. In the present oscillating case, the limit (weak) ǫ → 0 is more challenging, and has to be defined adequately. We refer the reader to [3,4,16] for other AP-scheme references. This paper is laid as follows. Section 2 deals with the presentation of a physical situation leading, after scaling and simplification, to the anisotropic transport equation (1.1). Two simplified models which will be studied in the following, are presented. Section 3 reviews the mathematical framework necessary to study the first toy model, and investigates the asymptotic limit ǫ → 0. Section 4 introduces several numerical schemes that we shall apply for the resolution of the first toy model. Then, we present the numerical results obtained with these schemes in Section 5 and the numerical analysis in Section 6. The last section is dedicated to the mathematical and numerical study of the second toy model which considers variable coefficients. A conclusion gives some hints for our upcoming work, concerning the more realistic Vlasov equation (2.4).

Physical motivation and toy models
Let us shortly say here some words about the physical motivation of the present work and introduce the two simplified models we shall investigate numerically in the next sections. These simplified models are caricatures of typical asymptotic regimes encountered in plasma physics, as for example the gyro-kinetic regime, and contain all the numerical difficulties arising in the more complex real physical systems.
The core tokamak plasma can be considered as collisionless, such that the most appropriate model for the description of its dynamics is the Vlasov equation for each particle species (α = e for electrons and α = i for ions), i.e.
where e α = ±e resp. m α are the particle elementary charge resp. mass and E(t, x) resp. B(t, x) are the electric respectively magnetic fields, determined self-consistently from Maxwell's equations. In the electrostatic case (given field B), Maxwell's equations have to be replaced by Poisson's equation where Φ is the electrostatic potential, related to the electric field E by E(t, x) = −∇Φ(t, x). For more details about the modelling of magnetically confined fusion plasmas, we refer the interested reader to the textbooks [2,10,13].
From a numerical point of view, solving the system (2.2)-(2.3) is rather arduous, due among others to its high dimensionality (6 dimensional in the phase space (x, v)) and to the presence of several time and space scales in the dynamics, introduced for ex. by the strong magnetic field B which confines the plasma in the tokamak. We shall be concerned in the present work with the multi-scale aspects of the kinetic problem, difficulties which are described mathematically by the following rescaling of the Vlasov equation for the ions (see [1,9,11,12,21] for the gyrokinetic scaling) where ǫ stands for the ratio of the particle cyclotron period to the observation time. The electrons experience the appearance of a second small parameter, related to the small electron to ion mass ratio m e /m i , leading to additional numerical burden, we shall not consider here (see [5]). The effect of the intense magnetic field on the particle dynamics is that it introduces a strong anisotropy, the motion of the charged particles being splitted into a fast gyration around the magnetic field lines and a slow dynamics along these lines, separation which necessarily causes numerical complications.
Let us introduce now two simplified toy models, which contain all the numerical difficulties of the initial model. In the rest of this paper we shall consider a homogeneous magnetic field B = b b with fixed direction b := e z and constant magnitude |B| = b ≡ 1. Furthermore, let us also introduce the following notation Sometimes it is more convenient to shift in (2.4) from Cartesian coordinates to polar coordinates for the velocity, i.e.
The Vlasov equation (2.4), written in polar coordinates, has then the form The two formulations, (2.4) resp. (2.5), corresponding to a Cartesian (not fieldaligned) resp. polar (field-aligned) configuration, are different from a numerical point of view, and different numerical schemes are usually employed for their resolution. To understand this difference better, we shall investigate in the present work in detail some numerical schemes for simplified versions of (2.4) and (2.5). We deliberately simplified these equations in order to be able to do a complete numerical analysis and to understand in all details the features of the here introduced AP-schemes.
2.1. First toy model -Polar, field-aligned configuration. Let us start from the Vlasov equation (2.4), assume here that E ≡ 0, B = e z and consider furthermore only the dynamics in the perpendicular plane (x, y), i.e.
where ǫ ≪ 1 accounts as usual for very strong magnetic fields. In order to simplify the computations, one often shifts to polar coordinates for the velocity, leading to where the unknown now is F (t, x, y, r, θ). We recognize thus a simple 3D anisotropic transport equation, the variable r being considered as a parameter in (2.7).
Choosing an initial condition F in independent on the variable y, would even lead to a more simpler 2D transport model This problem represents the simplest example of an anisotropic advection equation, to be understood in detail before designing an efficient scheme for the resolution of the Vlasov equation in the gyrokinetic regime (2.4). It is sufficiently difficult in order to study the behavior of the various schemes we shall introduce, and shall be the starting point of Section 3.

2.2.
Second toy model -Cartesian, not field-aligned configuration. In this second part, we shall differently simplify our Vlasov equation in order to study a different behavior. In particular, setting E ≡ 0, B ≡ e z and taking an initial condition independent on the space variable, yields the following 2D equation, in Cartesian coordinates or equivalently The difference of this model to the previous one is that this time the characteristics are no more straight lines but curves, such that the numerical schemes will behave differently. As mentioned earlier, these two models correspond to simplified versions of a field-aligned, polar coordinate framework , as well as a not field-aligned, Cartesian framework, both associated to the Vlasov equation (2.4) in the gyro-kinetic regime.
2.3. Aim of the present paper. The main points we are interested in within this study are the following: • design of AP-schemes for an efficient numerical resolution of anisotropic Vlasov equations of type (2.8), (2.9). Important properties we are asking from the schemes are: (a) stability independent on ǫ; (b) numerical diffusion/accuracy independent on ǫ; (c) discretization of the limit model as ǫ → 0; • show that taking the stiff term 1 implicitly is not sufficient for having an AP-scheme, meaning that AP-schemes are more than taking "implicitly" the suitable terms. AP-schemes have to mimic at the discrete level the precise asymptotic behavior of the solution in the limit ǫ → 0; • perform a detailed numerical analysis of the presented schemes in the framework of the two simplified toy-models (2.8), (2.9) and identify exactly which are the particularities of each scheme and each equation; • understand the difference between a field-aligned framework (2.8) and a Cartesian framework (2.9), and this from a numerical point of view; • prepare the foundation for a future, more realistic work, dealing with the resolution of the initial Vlasov equation (2.4) in the gyro-kinetic regime. Finally, let us say some words about Asymptotic-Preserving schemes. In general, inaccuracy in numerical simulations can result from applying unstable algorithms to wellconditioned problems or stable algorithms to ill-conditioned problems. Dealing with singularly-perturbed problems is a hard task, as they are ill-conditioned from the beginning. A standard, stable discretization (implicit in this case) often results in inaccurate results. The essence of AP-procedures is to replace singularly-perturbed problems by equivalent problems, which are regularly perturbed, well-conditioned problems, leading to uniformly accurate results, if stable algorithms are used (AP-approach).

First anisotropic Vlasov toy model and its mathematical study
Let us investigate now in detail the following simplified toy model, corresponding to a field-aligned anisotropic Vlasov equation where f in is a given initial condition, a > 0 and b > 0 are for the moment constants and 0 < ǫ ≪ 1 is a parameter representing the strong anisotropy/stiffness of the problem. Our computational domain is a doubly periodic box Ω : We shall review here some standard numerical schemes as well as introduce some new ones for the resolution of such a singularly perturbed problem and discuss finally their advantages and disadvantages. In particular, one is interested in numerical schemes capable to solve (3.11) uniformly accurate in ǫ, so-called "Asymptotic-Preserving" schemes. Let us however start with a detailed mathematical study of the behavior of (3.11).
3.1. Singularly perturbed problem. Equation (3.11) is a simple advection problem, whose exact solution is given by the characteristic method, i.e.
Remark that this function is L x -periodic in the variable x, L y -periodic in the variable y. Concerning the time-variable, two time-scales are present in the problem, a slow time-scale t and a rapid one t/ǫ.
The term b ǫ ∂ y f ǫ in (3.11) is the dominant term in the case where ǫ ≪ 1, such that passing formally to the limit ǫ → 0, yields This system, called "reduced system", is ill-posed. Depending on the initial condition, it can admit or an infinite number of solutions, namely if ∂ y f in = 0, or no regular solution (if ∂ y f in = 0). From a numerical point of view, this ill-posedness in the limit is translated into the singularity of the matrix of the linear system obtained by discretization of this problem. In particular, trying to solve (3.11) in a standard manner will necessarily lead to a linear system which degenerates in the limit ǫ → 0. This shall induce sever numerical problems. More adequate schemes are hence necessary for an efficient resolution of (3.11), as for example "Asymptotic-Preserving" schemes which are uniformly stable and accurate independently on the small parameter ǫ, and are additionally able to capture the limit model as ǫ → 0.
3.2. Limit model. For a better comprehension of our singularly-perturbed problem as well as for the construction of efficient "Asymptotic-Preserving" schemes, we have to identify the limit problem (V ) 0 of (3.11) and its solution denoted by f 0 . The information we get from the reduced model is that the limit-function f 0 has to be y-independent. With this information we introduce now the average of the function f ǫ with respect to the direction yf Integration of the equation (3.11) with respect to y yields ∂ tf ǫ + a∂ xf ǫ = 0, which is an ǫ-independent problem. Passing then to the limit ǫ → 0 leads to the advection equation with solution The system (V ) 0 is what we call "limit-system" of the anisotropic Vlasov equation (V ) ǫ , as shall be proved in the next section.
3.3. Weak convergence. So far, we proved the existence of a unique solution f ǫ for the system (V ) ǫ resp. f 0 for the limit system (V ) 0 . The next step is now to show the weak-convergence of f ǫ towards f 0 as ǫ → 0, and this in a certain sense. To define this sense, we have to introduce the right mathematical framework. In the sequel the symbol ♯ shall underline the periodicity of the considered space.
Proof. To prove (3.15), which signifies It follows that the function g belongs to H 1 ). The L y -periodicity of g is seen by the simple computation Taking now an arbitrary test function φ ∈ C 1 0 (0, T ; L 2 ♯ (Ω)) and introducing for simplicity for each f, g ∈ L 2 ♯ (Ω) the bracket f, g := Ω f g dx dy, we have where C > 0 is a constant independent on ǫ. Therefore, which concludes the proof due to the dense injection C 1 0 (0, T ; L 2 ♯ (Ω)) ⊂ L 1 (0, T ; L 2 ♯ (Ω)).

Numerical schemes for the anisotropic Vlasov equation
In this section we shall now introduce several numerical schemes for the resolution of (3.11) and examine them in more details. Firstly, different time semi-discretizations will be presented and then some words mentioned about a standard upwind spacediscretization. The time-discretization is the most important step in the construction of AP-schemes. For this, let us first introduce the following homogeneous discretizations of our time interval [0, T ] as well as of our simulation domain Ω = [0, L x ] × [0, L y ] : We shall denote further by f ǫ,n resp. f ǫ,n ij the numerical approximation of f ǫ (t n , x, y) resp. f ǫ (t n , x i , y j ). Recall also that we consider a doubly-periodic framework, such that The first time semi-discretization we shall study will be the implicitexplicit (IMEX) Euler method, where the stiff term is taken implicitly, i.e To study the behavior of this scheme, as ǫ becomes smaller, let us formally let ǫ go to zero in (4.17) and get This equation admits an infinite amount of solutions, namely all periodic functions dependent only on x. This formal analysis permits hence to conclude that the IMEX scheme can not be an AP-scheme, as it does not capture correctly the asymptotic behavior of the problem, which is rather given by the limit problem (V ) 0 . This property shall be tested numerically in Section 5.

Fourier method/Micro-Macro method.
A different way to solve (3.11) is to use a partial Fourier transform in the variable y, which is possible here, as we are in a simplified periodic context with constant coefficients. Denoting indeed the Fourier coefficients bŷ where the Fourier coefficients are solutions of the system

(4.19)
A simple discretization of this problem can be Solving this system and using the inverse Fourier transform (4.18) permits to get the desired result, i.e. the values of the unknowns f ǫ,n ij , solution of (3.11).
Let us investigate now the behavior of this system when ǫ → 0. Formally we get Therefore, we find a discretized version of the Vlasov limit problem (V ) 0 , signifying that this method will be "Asymptotic-Preserving".
The Fourier method is very nice, however it can be applied only in a simplified periodic framework with constant coefficients. As a sort of generalization one can think at the micro-macro method [3], which is based on the decomposition of each quantity in its mean part over the variable y, denoted by H ǫ or simplyf ǫ , and the fluctuation part h ǫ or simply (f ǫ ) ′ , defined as follows Taking now the average of the advection equation (3.11) over y and subtracting the resulting equation then from the initial one, yields a system to be solved for the unknowns (H ǫ , h ǫ ), i.e. (4.20) Let us study now the behavior of this system when ǫ → 0. We have formally (4.21) The two last equations establish that h 0 ≡ 0. Hence the system (MM) 0 is nothing else than the Vlasov limit system (V ) 0 . Again, we have created a scheme which is a regular perturbation of the asymptotic limit model, and shall be hence "Asymptotic-Preserving". This method is rather similar to Fourier method, however more general, as it can be applied in rather broad contexts. To understand this similitude, remark that H ǫ is nothing else than the first Fourier coefficient f ǫ 0 and the fluctuation h ǫ regroups the remaining Fourier modes. However, there is still a disadvantage or difficulty, namely the implementation of the constrainth ǫ = 0, which is crucial for the passage to the limit ǫ → 0. It is this constraint which permits in the limit to get a unique h 0 and to have thus a well-posed limit problem (MM) 0 . But averaging along the anisotropy lines can be very difficult in more general contexts, for ex. when these lines are not aligned with the axes.

4.1.3.
Lagrange-multiplier method. The Lagrange-multiplier method is based on the idea to replace the stiff, dominant term b ǫ ∂ y f by a smoother one ∂ y q, yielding the system where the inflow boundary is defined as Γ in := {(x, y) ∈ ∂Ω / y = 0}. In the limit ǫ → 0 one remarks that q ǫ is a sort of Lagrange multiplier corresponding to the constraint ∂ y f 0 = 0, where the name of the method.
First, we will prove the equivalence between the system (La) ǫ and the Vlasov equation (V ) ǫ , proving thus the well-posedness of the reformulation (La) ǫ . For this, let us first consider the unique solution f ǫ of (V ǫ ) and prove the existence of a function q ǫ such Then we shall decompose f ǫ in the following manner, which is somehow similar to a Hilbert Ansatz : To have a unique decomposition, we have to single out the G-part of q ǫ , by fixing for example q ǫ on the inflow boundary Γ in , choosing q ǫ ∈ Q with Obviously, we have G ∩ Q = {0 V }, implying the uniqueness of the decomposition (4.23).
Replacing now this decomposition in the system (V ) ǫ , we obtain directly the system (La) ǫ , which proves the existence of a solution to (La) ǫ . The converse is trivial, meaning that for (f ǫ , q ǫ ) ∈ V × Q solution to (La) ǫ , f ǫ solves (V ) ǫ . Altogether, we have proved the equivalence between both systems. Now let us consider the limit problem of (La) ǫ , obtained by letting formally ǫ → 0 in (4.22) The second equation leads to f 0 =f 0 . Then, averaging the first equation of (4.24) in the y-variable, yields where we used that q 0 is L y -periodic. This equation permits the determination of the limit function f 0 . Furthermore, the remaining well-posed system can be solved to assure finally the existence of the unique solution (f 0 , q 0 ) for the limit problem (La) 0 . The Lagrangian scheme seems to be the most "far-reaching" AP-scheme . The only disadvantage of this method is that we have now two unknowns and hence two equations to be solved, meaning longer simulation times. However, we are no more forced to follow the anisotropy lines and can choose coarse Cartesian, not-field aligned grids.

4.2.
Space discretization for the IMEX scheme. For any numerical scheme presented above, we decided to consider the standard upwind method to discretize the transport terms in the equation (3.11). The idea behind this choice is that the spacediscretization is not the important step in the construction of an AP-scheme, such that we opted for a simple discretization, in order not to embroil the further numerical analysis as well as the understanding of the main ideas of our methods. The same arguments incited us to select only first order discretizations in time. A Runge-Kutta coupled to a second-order space-discretization would be naturally more accurate, changes however nothing in the essential concept of our AP-strategies. As mentioned earlier, in a forthcoming paper we shall be concerned with a realistic, fusion plasma situation, such that we shall adapt the most adequate of the here presented schemes to more accurate second order techniques, to gain in accuracy. Now, let us recall the first-order upwind forms We have analogous formulae for the partial derivative in the y-variable. Denoting now α := a∆t ∆x > 0 and β := b∆t ∆y > 0 and using the periodicity, i.e.
. the completely discretized IMEX scheme writes finally : . We remark that we can rewrite this scheme like a system of N x − 1 equations : where : . . .
At each time step, we resolve this system ∀i ∈ [1, N x − 1], to get the unknowns f ǫ,n+1 i,j . Remark that A = ǫ Id + C β is a regular perturbation of a singular, cyclic matrix C β .

Numerical simulations
In this part, we shall test numerically every scheme introduced in the previous Section for the resolution of the anisotropic Vlasov equation (3.11). The homogeneous time and phase-space discretization was previously introduced in (4.16) and we choose in the sequel the following parameters: T = 1, L x = 2π, L y = 2π, N t = 101, N x = N y = 201, a = 0.1 and b = 1. Changes in these parameters shall be explicitly mentioned. The initial condition we adopt is given by : We recall that the exact solution of (3.11) is known and reads, for each ǫ > 0: In Figure 1, we reveal two graphics which contain on the one hand f in and on the other hand f ǫ ex at the final time T = 1. Furthermore, in order to better figure out our problem, we plotted in Figure 2 the exact solution of the limiting Vlasov system (3.14) at the final time T , i.e. f 0 ex (T, x) = f in (x − aT ). Remark that this solution is homogeneous in the y-variable.
Finally, we show in Figure 3 the time-evolution of the exact solution f ǫ ex at one point only, i.e. (x Nx−1 , y Ny−1 ). We distinguish easily on the left plot (A) of Fig. 3 the two periods, one linked with the x-variable, and the other one corresponding to the y-variable. This last one is ǫ-dependent and we see that more ǫ is small, more the frequency of the timeoscillations becomes important. As the 2D situation is not so eloquent, we eliminate the x-variable in the problem and considered also a 1D problem, keeping only the term containing the parameter ǫ (i.e. a = 0). The time-evolution of the exact solution at the point y Ny−1 is now plotted in Fig. 3 (B). One observes here more easily that with smaller becoming ǫ, the frequency of the time-oscillations is increasing. In the limit ǫ → 0, f ǫ (t, y Ny−1 ) converges weakly towards the average, which is here the constant 1.

5.1.
Some results obtained with our schemes. Now we examine how the different numerical schemes introduced above cope with such an asymptotic behavior. For ǫ = 1, we recognize an approximation of the exact solution (see Figure 1) and for ǫ = 10 −10 , the limit solution is clearly obtained (see Figure 2). Briefly one can  say that the numerical solution follows the weak-⋆ convergence f ǫ ⋆ ⇀ ǫ→0 f 0 as ǫ becomes smaller and smaller. But, one can remark a numerical diffusion which leads to a loss of amplitude, especially visible in the non-limit case ǫ = 1 or ǫ = 10 −1 . To observe better this numerical diffusion, we show in the right plot of Fig. 5 the time-evolution of just one point of the numerical solution, corresponding again to a 1D situation as the one plotted on the right of Fig. 3, and this for several values of ǫ. As one can observe,  the damping is more and more pronounced if ǫ → 0. For small ǫ-values the numerical solution recovers quasi immediately the weak limit solution, here the constant 1. This damping phenomenon will be understood from the numerical analysis we shall fulfill in Section 6. 5.1.2. Fourier, Micro-Macro and Lagrange-multiplier schemes. Let us now present analogous results for the remaining schemes, namely the Fourier, Micro-Macro and Lagrangemultiplier schemes. The 2D plots are rather similar to the ones presented for the IMEXscheme (see Fig. 4-5). To examine the difference between these methods, we preferred to plot in Fig. 6 only the time-evolution of the numerical solution in the 1D-context again. We remark that the damping of the Fourier method is more slowly than the ones of the IMEX-scheme as well as Micro-Macro and Lagrange-multiplier scheme (which are completely overlapping). But, once again we observe that in the limit ǫ/t → 0, the fluctuations are completely damped out and we recover the weak limit solution.

5.2.
Convergence of the schemes for fixed ǫ > 0. Let us study now the convergence of the here presented schemes with respect to time and space, and this for fixed ǫ > 0, permitting to show their validity in the large ǫ-regime. For this, fix ǫ > 0 and consider the error between exact and numerical solutions as a function of the mesh-size, at the final time T. Firstly, concerning the convergence with respect to ∆t, we choose small space steps (N x = N y = 501) such that the space errors are much smaller than the time error and vary then the time step. Equally we apply the same strategy for the convergence with respect to ∆x and ∆y, by fixing a time step of N t = 501. In all cases, the parameter ǫ is fixed to 1. In Figure 7, we have plotted curves in log-log scale, showing the evolution of the errors as a function of ∆x, ∆y and ∆t, respectively.
As expected, we observe that all schemes are first order in time and space. Some comments are however necessary to understand Figure 7. First, the slop of the curves gets smaller than 1 in the small-grid ranges. This is due to the fact that the error to be investigated (for ex. in ∆t) becomes as small as the fixed error term (in ∆x, ∆y) and saturates. Secondly, the slope of the curves becomes also smaller in the large-grid ranges. This is usual, as for large discretization steps, the rest-terms in the Taylor series for the error analysis can no longer be neglected. Finally, we would like to draw the attention of the reader to the Fourier error curve, which has a constant slope in (B). This is completely natural, as the Fourier method has spectral accuracy.

5.3.
Asymptotic behavior as ǫ → 0. To begin the study of the asymptotic behavior, we define the following two errors where η ǫ (t) represents the L ∞ -error between the exact and the numerical solution at instant t, for fixed ǫ > 0, whereas γ ǫ (t) denotes the L ∞ -error at instant t between the numerical solution f ǫ num and the exact limit solution f 0 ex .
We are interested in the evolution of these two errors at the final time T as functions of ǫ. The curves corresponding to the different schemes are plotted in Figure 8. As expected, we observe a decrease of η ǫ (T ) and an increase of γ ǫ (T ) when ǫ → 1. For ǫ → 0 the converse behavior is observed. This plot shows that each scheme approximates well either the exact solution f ǫ ex for large ǫ, or the exact limit solution f 0 ex for small ǫ.
What can be said as a conclusion, is that all schemes seem to have the right asymptotic behavior in this simple test case. Indeed, for fixed ǫ > 0, each numerical solution f ǫ num converges to the expected solution f ǫ ex as long as the grid is refined (Fig. 7). For fixed discretization steps, all numerical solutions f ǫ num converge towards the limit solution f 0 when ǫ becomes smaller and smaller, underlying the AP property of our methods. It is worth mentioning however that the IMEX-scheme is no more working for ǫ smaller than 10 −14 , the matrix A of the IMEX linear-system (4.27), namely is becoming numerically singular in the limit ǫ → 0. This is not the case for the Micro-Macro as well as Lagrange-multiplier schemes, which give accurate results even for ǫ = 0. This difference in the behavior can be observed also from the study of the conditionnumber of the discretization matrices, paying attention especially on the ǫ-dependence. Remark here that an "Asymptotic-Preserving scheme" must have an ǫ-independent condition number, depending merely on the discretization parameters ∆x, ∆y.
In Fig. 9 we plotted thus the matrix condition-number cond(A) := ||A −1 || 2 ||A|| 2 corresponding to the three schemes (IMEX, Micro-Macro and Lagrange-multiplier) as a function of ǫ. What can be observed is that for the Micro-Macro and Lagrange-multiplier scheme, the condition-number is ǫ-independent (for ǫ ≤ 10 −2 ), which is a hint of the well-posedness of these two problems in the limit ǫ → 0, namely of (MM) 0 resp. (La) 0 . On the other hand, for the IMEX-scheme cond(A) is proportional to 1/ǫ (slope of the curve is approx. −1). This circumstance is the translation on the discrete level of the fact that the reduced model (3.13), obtained on the continuous level by letting formally ǫ → 0 in the IMEX time-discretization, is ill-posed, admitting an infinite amount of solutions.
However, even if these arguments show clearly that the IMEX-method should behave badly for very small ǫ-values, it is not the case in our simplified toy model, in particular it does not seem to be affected by the bad condition number. This will no more be the case in our second toy-model. To understand in detail what happens, a more refined error study could be profitable and shall be done in the next section. The final interpretation is postponed to Section 6.3 after having estimated the truncation error. One can only say here that the functioning of the IMEX-scheme is due to the fact that the investigated problem is very simple and specifically the anisotropy is aligned with the Cartesian mesh.

Numerical analysis
Let us now perform a numerical analysis study of our schemes introduced for the resolution of (3.11), permitting to understand in detail the behavior observed in the last section. In particular we shall detail only the error-analysis of the standard IMEXscheme and the Asymptotic-Preserving Lagrange-multiplier scheme. The error study of the other schemes is very similar. See [14,17] for more details on this analysis part.
6.1. IMEX scheme. We begin by recalling the full discretized form of the IMEX scheme : Finally, we observe that the IMEX scheme (6.28) is a second-order scheme for the modified Vlasov equation Proof: The local truncation error of the method (6.28) is defined by ∆y .
Supposing that f ǫ is sufficiently smooth in order to apply a Taylor development, we find where f ǫ is taken in (t n , x i , y j ). Since f ǫ satisfies the Vlasov equation (3.11), the O(1) terms drop out. Moreover, by differentiating the Vlasov equation along t, y and x, we express the partial derivatives ∂ tt f and ∂ ty f as functions of ∂ xx f and ∂ yy f . We find thus The local truncation error writes finally Remark 6.2. The modified equation (6.29) is an advection/diffusion equation. Note that the diffusion is stronger in the y-direction due to the term 1/ǫ. These diffusion terms are responsible for the damping that we observed in the numerical simulations (see Fig.  5 (B)), damping which tends towards infinity in the y-direction, as ǫ → 0. Note also that the diffusion coefficient is positive if α ≤ 1. This is precisely the stability condition of the upwind scheme, as we will see afterwards. If this condition is not respected, the diffusion becomes negative, leading to an ill-posed problem with exponentially growing solutions. Proof: To study the stability of our scheme, let us inject in (6.28) for fixed n ∈ N a plane wave of the form with k, l ∈ Z two arbitrary modes, and look how it evolves from one time-step to the other. Let us denote by ξ I the amplification factor for this passage t n → t n+1 , meaning Inserting now these terms in the discretized equation (6.28), yields, after simplification .
A scheme is said to be stable in the Von Neumann sense, if the amplification factor satisfies |ξ I | ≤ 1, such that the modes are not amplified from one time-step to the other. Straightforward computations yield Then, a necessary and sufficient condition to have the Von Neumann stability is : a∆t ∆x 1.
Remark 6.4. Note that in the case l = 0, when ǫ tends towards 0, the amplification factor converges towards 0. This means that for injected waves with mode l = 0, the scheme becomes more and more diffusive and attenuates completely the oscillations.
6.2. Lagrange-multiplier scheme. We do now the same work for the Lagrangemultiplier scheme, i.e. Proof: In order to prove the result, we write the local truncation error of the first equation. We find that Since the first equation of (6.30) is verified by (f ǫ , q ǫ ), we have Then, Since the second equation of (6.30) is verified, we have such that we find the same expression as for the IMEX scheme, i.e.
The just proved result confirms what we have seen on the numerical plots. Indeed, the IMEX and Lagrange-multiplier schemes have the same behavior when regarding the convergence and the asymptotic behavior.
Theorem 6.6. The Lagrange-multiplier scheme is stable in the Von Neumann sense if and only if the CFL condition a∆t ∆x 1 is satisfied.
Proof: Here, we have two unknown functions f ǫ and q ǫ . To study the Von Neumann stability, we write q ǫ,n+1 = ξ f f ǫ,n i,j , with the two amplification factors ξ q and ξ f . As usual, we insert these expressions in the discretized Lagrange-multiplier equations. We obtain a linear system where the unknowns are ξ q and ξ f . This system writes and is easy to invert. Computing the amplification factor ξ f , we remark that it is identical to the one calculated for the IMEX scheme.
6.3. AP-properties. We are now able to explain the numerical results obtained in Section 5, in particular to explain why the IMEX-scheme, even if being not an APscheme, gives in this simple field-aligned test case, good results up to a value of ǫ = 10 −14 .
For this, let us recall that two types of errors arise during a numerical resolution of the Vlasov equation (3.11). First of all we have the truncation errors, estimated in the last subsections, and secondly one has also the round-off errors, arising at each elementary computation. To be more precise, one has to consider the three linear systems, corresponding to (4.27): where to simplify notation we omitted all the time and space indices. Here we denoted by F ex the exact solution of the Vlasov equation (3.11), satisfying the linear system (4.27) up to a truncation error T , F is the exact solution of the linear system (4.27), supposing exact arithmetics, and finally F num is the solution to the linear system (4.27) obtained via a computer, hence contaminated with round-off errors. The error we are interested in, can be estimated as follows Stability and consistency permit to show that the first error term is of the order of the truncation error. For the estimate of the second error term, we have to take into account the condition number of the matrix, in particular one has the estimate [23] Performing our computations in double precision (machine accuracy of 10 −16 ), and as long as the condition number is not exceeding a value of 10 12 (see Fig. 9), the second error term is not so dangerous. For larger condition numbers, this term can give rise to erroneous results. In our test case, it is however rather the first error-term which leads to trouble, as the truncation error is 1/ǫ-dependent. In the first toy-model (3.11), the large truncation error impacts only the y-direction, leading to a large diffusion along the axes-aligned anisotropy and hence to the limit-model. We shall see a drastic difference in the second, not-field aligned toy-model.

Second Vlasov toy-model with variable coefficients
Finally, let us come now in this section to the second Vlasov toy model, given by :  We can write this system under matrix form : Denoting the rotation matrix by R ǫ (y) := e A y ǫ , one has We can easily verify that the characteristic curve passing through the point (x, y) is a spiral, whose projection on the (x, y)-plane is the circle with radius R := x 2 + y 2 and center (0, 0). All characteristics are 2 π ǫ-periodic (in t).
The exact solution f ǫ of (7.32) is now simply the advection of the initial condition along these characteristic curves, such that f ǫ (t, x, y) = f in (X(0, t, x, y), Y (0, t, x, y)) = f in cos t ǫ x−sin t ǫ y, sin t ǫ x+cos t ǫ y .

7.2.
Limit solution of the problem. The next step is to obtain the limit solution of the problem (7.32), as ǫ → 0. Keeping in mind that f ǫ is constant along the characteristic curves, we integrate (7.32) along C x,y ǫ , to get which comes from the periodicity of the characteristics, and denoting the average along a curve by f ǫ := 1 |C x,y ǫ | C x,y ǫ f ǫ dσ, with |C x,y ǫ | = 2 π ǫ, we have : Letting now formally ǫ → 0, we obtain the following limit problem associated to (7.32): (G) 0 f 0 = f in .
(7.35) The term (∆x∆y) γ q ǫ,n in (7.35) is a stabilization term permitting to have the uniqueness of the solution (f ǫ , q ǫ ). In the former "field-aligned" example, we fixed q ǫ on the anisotropy lines by setting q ǫ |Γ in = 0, but here it is more arduous from a numerical point of view. The stabilization aims equally to fix q ǫ , however in a different manner. It is very delicate to choose the magnitude of this term, in order not to destroy the problem, in particular we took here γ = 0.91. First it is a small perturbation of the equation, of the order of the truncation error. Secondly, averaging the second equation of the Lagrange-multiplier scheme along the anisotropy lines, permits to obtain (∆x∆y) γ q ǫ,n+1 = 0, which means that q ǫ is unique, by having zero average along the field lines. A more detailed study of this stabilization technique was performed in [18] for the elliptic framework. For the spatial discretization, we use again an upwind scheme, observing that this time the equation has no more constant coefficients. Thus, we define : The full discretization of the IMP scheme writes now with r x = ∆t ∆x and r y = ∆t ∆y . And for the Lagrange-multiplier scheme, we have :

Numerical simulations.
Here we present our simulations corresponding to both numerical schemes. We consider Ω = [−3, 3] 2 , T = 1 and the discretization parameters N t = 64 and N x = N y = 160. The initial data is defined by a Gaussian function : As we showed before, the exact solution is known thanks to the characteristic method.
In the present simple test case, one can easily prove that in other words, the exact solution is a stationary solution, independent of ǫ, the initial condition being constant along the anisotropy field. This simple test case permits in a very simple way to compare both methods with respect to the ǫ-dependence of the results, in particular to show that the IMP-scheme is not an Asymptotic-Preserving scheme. We shall investigate in a future paper a more involved, physical test-case, where we shall adapt the here introduced Lagrange-multiplier-method, which seems to be the most appropriate method for our singularly-perturbed Vlasov problem (2.4), to secondorder schemes and test more thoroughly its AP-properties.
In Figure 10 we first plot the condition-number cond(A) := ||A −1 || 2 ||A|| 2 associated to the two schemes. As for the first toy-model, one remarks the ǫ-independent conditionnumber of the Lagrange-multiplier-scheme, whereas, as expected, the IMP scheme has an 1/ǫ-dependent condition-number. The two curves correspond to the IMP and Lagrange-multiplier schemes.
Then, in Figure 12, we show the numerical solution f ǫ at the final time T and computed for several values of ǫ, with both IMP and Lagrange-multiplier schemes. For ǫ = 1 and ǫ = 0.1, we do not distinguish any difference. However for smaller ǫ values, the solution obtained with the Lagrange-multiplier scheme seems to be ǫ-independent, contrary to the IMP scheme, which diffuses more and more as ǫ → 0. Indeed, the IMP solution is completely damped as ǫ → 0 and leads towards the zero-solution, whereas the Lagrangemultiplier scheme keeps the form of the Gaussian, with a usual ǫ-independent (∆x, ∆y)diffusion. This permits to conclude that the Lagrange-multiplier scheme is an APscheme contrary to the IMP scheme.
In order to distinguish much better this AP-property, we plot on Figure 11 a cut of the previous curves at the point x = 0. We observe clearly the diffusion in the IMP scheme which depends of ǫ contrary to the Lagrange-multiplier scheme. 7.5. Numerical analysis. The aim of this section is to explain the plots presented before. In particular we will investigate why the IMP scheme does not work for small ǫ-values, whereas the Lagrange-multiplier scheme preserves the asymptotics. First of all, we compute the local truncation error of both schemes. We shall consider only the case x ≥ 0 and y ≥ 0, the remaining cases changing nothing in the following reasoning. 7.5.1. IMP scheme. We begin by recalling the expression of the full discretized expression of this scheme : i,j ∆y = 0 , ∀(n, i, j) ∈ Q h . Theorem 7.1. The IMP scheme (7.37) is consistent with the second Vlasov problem (7.31), first order accurate in time and in space. Moreover, the local truncation error writes where Proof: This proof is very similar to the proof of Theorem 6.1.
Remark 7.2. Contrary to the first toy-model, where the diffusion was 1/ǫ-dependent only in the anisotropy-direction, which was aligned with the coordinate system, in the present case, the diffusion-matrix is scaled by a 1/ǫ factor, meaning that this time we have a very strong 1/ǫ-dependent diffusion in all directions. This large diffusion leads rapidly to a damping of the solution towards zero, as ǫ becomes smaller, and leads thus to completely erroneous results. 7.5.2. Lagrange-multiplier scheme. We use the same reasoning for the Lagrange-multiplier scheme Supposing y 0 and x 0, we have the full discretization of (La) ǫ (7.39) Theorem 7.3. The Lagrange-multiplier scheme (7.39) is consistent with the second Vlasov model (7.31) and first order accurate in time and in space. Furthermore, the local truncation error writes Proof: We begin by the computation of the T L1 term. Supposing sufficient regularity for the functions f ǫ and q ǫ , we use Taylor series expansion to get T L1 (t n , x i , y j , ∆t, ∆x, ∆y) = ∆t 2 Since the first equation (7.38) is verified, we can write ∂ tt f ǫ = −y j ∂ xt q ǫ + x i ∂ yt q ǫ , And we differentiate in time the second equation of (7.38) to obtain T L1 (t n , x i , y j , ∆t, ∆x, ∆y) = ∆t 2ǫ We have the following relations ∂ tx f ǫ = x ∂ yx q ǫ − y ∂ xx q ǫ + ∂ y q ǫ , ∂ ty f = x ∂ yy q ǫ − ∂ x q ǫ − y ∂ xy q ǫ .
Remark 7.4. In contrast to the first Vlasov toy-model (3.11), the IMP and Lagrangemultiplier schemes do not have the same behavior with respect to the local truncation error. More particularly, the dependence on ǫ is very different. The IMP-scheme is diffusing in all directions, diffusion proportional to 1/ǫ. The only 1/ǫ-dependent diffusion in the Lagrange-multiplier scheme arises in relation with the auxiliary unknown q ǫ , i.e. in the term ∇ (D L 1 ∇q ǫ ). And one can immediately verify that the 1/ǫ-dependence arises only aligned with the anisotropy field lines, and not perpendicular to them. Indeed, one gets immediately for the diffusion along resp. perpendicular to the field lines: (y , −x) D L 1 (y , −x) T = y 3 ∆x 2 + x 3 ∆y 2 + ∆t 2 (x 2 + y 2 ) 2 (x , y) D L 1 (x , y) T = x y 2 [x ∆x + y ∆y] .

Concluding remarks
To conclude, let us summarize here the knowledge we acquired about the resolution of anisotropic Vlasov equations of the type (2.4) arising in fusion plasma modelisation. Two types of techniques can be adopted from the beginning. One can decide to pass directly to polar coordinates in velocity and get hence a field-aligned formulation as for ex. (2.7). In this case, a simple IMEX-scheme is the most appropriate scheme to be used, being simple enough and giving rise to accurate results up a sufficiently small ǫ-value. However, the disadvantage is that one has to change coordinate system, which can be rather cumbersome if the magnetic field is variable, in time and space. The second technique is rather simple, as it avoids to pass to field-aligned coordinates and remains in a nice Cartesian framework. The drawback is that in this case it is no more sufficient to implicit the stiff term and take the other terms explicitly. Indeed, for small ǫ-values (already ǫ = 10 −4 ), meaning strong magnetic fields as in tokamak plasmas, an IMEX scheme would lead to erroneous results. An Asymptotic-Preserving reformulation like our "Lagrange-multiplier-method" is more adequate, leading in the limit ǫ → 0 towards the right Limit-problem. This Lagrange-multiplier-method is indeed usable for all ǫ ≥ 0 and gives accurate and stable results independently on ǫ. However there is a disadvantage, namely the fact that it is more time-consuming, as it involves an additional unknown q ǫ . Solving an anisotropic Vlasov equation of the type (2.4) needs hence an a priori decision, which of these two techniques to follow. The first technique is at the moment the basis of several codes. The second technique has not be tested up to now, and its rigorous validation and comparison with the first one will be the aim of a forthcoming paper, in a more physical context. E-mail address: baptiste.fedele@math.univ-toulouse.fr, claudia.negulescu@math.univ-toulouse.fr