{'fraction_non_alphanumeric': 0.05832063092318908, 'fraction_numerical': 0.022035920206773148, 'mean_word_length': 3.7986961361106695, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Uniformly accurate methods for Vlasov equations with non-homogeneous strong magnetic field

In this paper, we consider the numerical solution of highly-oscillatory Vlasov and Vlasov-Poisson equations with non-homogeneous magnetic field. Designed in the spirit of recent uniformly accurate methods, our schemes remain insensitive to the stiffness of the problem, in terms of both accuracy and computational cost. The specific difficulty (and the resulting novelty of our approach) stems from the presence of a non-periodic oscillation, which necessitates a careful ad-hoc reformulation of the equations. Our results are illustrated numerically on several examples.


Introduction
In this article, we are concerned with the numerical solution of the 2-dimensional Vlasov equation with non-homogeneous magnetic field [2,3,15,16,24]. More specifically, if b : x ∈ R 2 → b(x) ∈ R and f 0 : (x, v) ∈ R 2 ×R 2 → f 0 (x, v) ∈ R are given functions, and if 0 < ε ≤ 1 denotes a dimensionless parameter and ]0, T ] a non-empty interval of time, we shall consider the Cauchy problem for the distribution function f ε : (t, x, v) ∈ [0, T ] × R 2 × R 2 → f ε (t, x, v) ∈ R given by (ii) the electric-field function E ε : (t, x) ∈ R + × R 2 → E ε (t, x) ∈ R 2 is either external and explicitly given or self-consistent. In the latter case, E ε solves the Poisson equation where n i denotes the ion density of the background.
Solving equation (1.1) with standard methods is notoriously difficult for vanishing values of the parameter ε, as the Lorentz term then creates high-oscillations in the solution: this indeed imposes to use tiny time-steps (usually of the order of ε) and leads to formidable computational costs. Hence, it is now admitted that specific techniques are required that can cope with this particular regime of small values of ε and which, as logic dictates, preserve the asymptotics of f ε in the limit where ε goes to zero. Numerical methods obeying to this paradigm (i.e. consistent with the limit equation for f 0 ) and that are consistent with (1.1) when ε = O(1) have been called asymptotic preserving methods and may be found in various publications [17,18,19,21,23].
Nevertheless, if the value of ε is not known prior to the simulation, it is often observed that the error behaviour of asymptotic-preserving methods is largely deteriorated for certain (not so small) values of ε. As a consequence, it appears highly desirable to design numerical methods for (1.1) which are uniformly accurate (UA) with respect to the parameter ε ∈]0, 1]. That is to say, p th -order methods which, when used with time step ∆t, deliver approximate solutions f ε ∆t such that f ε ∆t − f ε ≤ C(∆t) p in an appropriate function-norm, where the constant C as well as the computational cost are independent of ε.
It is precisely the aim of this paper to introduce UA schemes for equation (1.1), which, as we shall illustrate numerically, are indeed able to capture the various scales occurring in the system while keeping numerical parameters (in particular the time step) independent of the degree of stiffness ε. Although alternative options are possible [4,5,9,10,11], the strategy we develop to reach this goal is very much inspired by the recent papers [7,14]: its main underlying idea consists in separating explicitly the two time scales naturally present in (1.1), namely the slow time t and the fast time t/ε. This is done at the level of the characteristic equations (resulting from the use of the Particle-In-Cell method, see e.g. [1,26,22]), which are, for each macro-particles, stiff ordinary differential equations of the form ẋ = v, where the term b(x) ε Jv is the source of high-oscillations and at the origin of numerical difficulties. As compared to our previous works [7,9], the main obstacle we are confronted with (and accordingly the main novelty of the proposed solution) is the fact that the aforementioned oscillations are not per se time-periodic (see also [6]). However, we will show that the trajectory in the physical space x(t) remains confined within an ε-neighbourhood of the initial condition x 0 , allowing to regard a system in a form which is now amenable to the embedding ofx(s) andỹ(s) into the functions X(s, τ ) and Y (s, τ ) periodic with respect to τ and such that X(s, s/ε) =x(s) and Y (s, s/ε) =ỹ(s). The resulting transport equations then need to be complemented with initial conditions X(0, τ ) and Y (0, τ ). Their choice is the fundamental ingredient of the two-scale strategy proposed in [7,12] and it requires to be handled here with additional care owing to the presence of the 1/ε-term in the right-hand side of (1.6). Under this form, the problem shares similarities with the model analyzed in [7]. However, (1.6) contains two main additional difficulties due to the presence of the term (b(X)/b(x 0 ) − 1)JY : first, this nonlinear term prevents from a direct application of Gronwall lemma; second, this term is not smooth with respect to the unknown. Finally, we will see that the numerical solution also enjoys at a discrete confinement property, i.e. it remains confined within an ε-neighbourhood of the initial data. The organisation of the paper follows closely the steps exposed above. The two-scale formulation of the characteristics is introduced in Section 2, which includes in particular Subsection 2.1 devoted to the scaling and filtering operations and Subsection 2.3 devoted to the detailed derivation of the initial conditions of (1.6). The section concludes with rigorous estimates of the derivatives of X and Y (see Subsection 2.4) and a numerical confirmation of the expected smoothness of the solution is brought in Subsection 2.5. Section 3 is concerned with the effective derivation of numerical schemes of first and second orders for solving equations (1.6) and the proof of their convergence (Theorem 3.1) which constitutes the main result of this paper. Within Subsection 3.3 the adaptation of our numerical strategy to the situation of a coupling of Vlasov equation with Poisson equation is considered: although no rigorous statement is established at this stage, we provide empirical evidence of the efficiency of our method in this case. Section 4 describes several examples and the corresponding numerical experiments, confirming the interest of the technique.

Two-scale formulation of the characteristics equations
Using the Particle-In-Cell (PIC) discretisation, we get the characteristic equation for 1 ≤ k ≤ N p , We see (2.2) is a solution dependent highly oscillatory problem. As a basic requirement throughout the paper, we consider the magnetic field function b(x) is uniformly above zero, i.e. for some constant Scaling of time and filtering. For each k, introduce the scaled time and definex k (s) := x k (t),ṽ k (s) := v k (t), s ≥ 0, where we omit the subscript k in s k for simplicity of notations. Note that under the assumption (2.3), s k = s k (t) is a monotone increasing function, which is interpreted as a time for particle k. Through this transformation, we make each particle living in its own time. Then we can rewrite the characteristics equation (2.2) aṡ and then filter it out by introducing We shall analyse and solve (2.7) up to any fixed time For technical reasons, hereafter we shall assume that the given electric field E ε (t, x) and the magnetic field b(x) are globally Lipschitz functions, i.e.
for two constants C E , C b > 0 independent of ε. Here and after, the norm | · | of a vector always refers to the standard euclidian norm in R 2 , whereas it refers to the absolute value when it is applied to a scalar quantity.
, for the solution of (2.7), we have that for each 1 ≤ k ≤ N p , for some C 3 > 0. The constants C 1 , C 2 , C 3 depend on k but are independent of ε.
Proof. Based on the assumption for E ε (t, x) and b(x), we have the global well-posedness of (2.7) andx k (s),ỹ k (s) ∈ C(R + ). Indeed, taking the inner product on both sides of (2.7a) and (2.7b) with x k (s) andỹ k (s) and applying Cauchy-Schwarz inequality respectively gives Hence we find Then, adding the two last inequalities and by Gronwall's inequality, one can get an a priori estimate for boundedness of the solution, i.e.
By applying the Duhamel's principle to (2.7a) and then integrating by parts, we havẽ Hence from (2.7b), we havẽ Then we can see for all 0 ≤ s ≤ S k , By Gronwall's inequality, we get estimate for a constant C 1 > 0 independent of ε. The last assertion (2.10) then follows from Thanks to Lemma 2.1, we observe that the right-hand-side of (2.7) is bounded as ε → 0. If we consider E ε (t, x) is a given external field and contains no fast frequency in the t-variable, the formulation (2.7) fits the requirement of the two-scale strategy.
2.2. Two-scale formulation. Let us first of all consider the case that E ε (t, x) is a given external field without any fast frequencies in the t-variable. We shall address the Vlasov-Poisson case later. Now we perform the two-scale formulation on (2.7). Denote the fast variable τ = s/ε and separate it out in (2.7), then we have , we recover the original unknown by considering the two-scale unknown on the diagonal (2.12) In the next section, numerical schemes will be proposed for the two-scale system (2.11), and our aim will be to prove that these schemes enjoy uniform accuracy with respect to ε. This property requires a preliminary analysis. Indeed, one can observe that no initial condition for (2.11) is evident since only the condition X k (0, 0) = x k (0), Y k (0, 0) = v k (0) is required. This degree of freedom will be used to derive initial conditions X k (0, τ ), Y k (0, τ ) such that the two-scale unknown X k (s, τ ), Y k (s, τ ) and its time derivative are uniformly bounded. This will be the objective of the rest of this section.
First, we start with the following elementary lemma Lemma 2.2. Let H be a Banach algebra space of real-valued functions and let H 2 denote the cartesian product H × H. Consider the following system of ordinary differential equations in H 2 (i.e. X ε and Y ε are considered as functions from [0, with α ε (s), β ε (s) are 2 × 2 matrices with coefficients in H and γ ε : [0, S k ] → H 2 . Assume that there exists a constant C > 0 independent of ε such that α ε (s) H ≤ C, β ε (s) H ≤ C and γ ε (s) H ≤ C, ∀s ∈ [0, S k ]. Then there exists a constant M > 0 independent of ε such that Proof. We integrate the equation on X ε and perform an integration by parts to get where we used the equation on Y ε . We always use C in proofs to denote a positive constant independent of ε, and its value may change from one line to the next. Considering the norm · 2 H leads to which reported in the former inequality gives The two last inequalities clearly lead to A standard Gronwall lemma enables to prove the result of the lemma.
For convenience, we shall denote L ∞ s := L ∞ s ([0, S k ]) and L ∞ τ := L ∞ τ (T) the functional spaces in s and τ variables. Now, for a smooth periodic function u(τ ) on T, we introduce For a smooth vector field U(τ ) on T, we define its W 1,∞ τ -norm as . For the two-scale problem (2.11), if initially for some constant C 0 > 0 independent of ε, then 14) for some C 1 > 0 independent of ε. Hence,
Proof. Consider τ as a parameter in T and define The two scale problem (2.11) then readṡ Using the same strategy as in the proof of Lemma 2.1, we have Taking now the derivative with respect to τ of (2.19) and denotinĝ we geṫ by assumption, which concludes the proof.
2.3. Suitable initial data for the two-scale formulation. In this subsection, we look for an initial data X k (0, τ ), Y k (0, τ ) of the two-scale formulation (2.11), which will ensure that the time derivatives of the solutions X k (s, τ ), Y k (s, τ ) are uniformly bounded. This will be done using Chapman-Enskog expansion of the solution. First order preparation. We perform the Chapman-Enskog expansion to get the full initial data X k (0, τ ), Y k (0, τ ) for (2.11). This will be done by formal arguments and a rigorous statement will be proved in the next subsection. Denote with the average operator Π defined for some periodic function u(τ ) on T as Πu = 1 and Taking the inverse of L, which for a zero average function u(τ ) is computed as τ 0 u(θ)dθ, on the above equations for h k and r k and denoting A : Thanks to Lemma 2.3, we have b(X k )/b k − 1 = O(ε), and consequently from (2.25b) we get We now take the time derivative of (2.25) We then obtain the first order asymptotic expansions from (2.25) To determine h k and r k at s = 0, we then need to compute X k (0) and Y k (0). These quantities will be determined from the initial conditions x k,0 and v k,0 . Indeed, we recall that so that we get the following first order expansion for X k We then derive a first order expansion for Y k We combine this identity with (2.27) to get the first order approximation of Y k Second order preparation. We continue the preparation of initial data to the second order in ε. (2.27). We then get the following second order expansion for X k where, using (2.27) we define where h 2nd k is given by (2.34) and where X 1st k , Y 1st k are given by (2.29) and (2.32). Let us deal with the second order expansion of Y k . From (2.25b), we get Therefore, it remains to find an expansion of ∂ s r k (given by (2.26b)) up to order 1 in ε. To that purpose, we use (2.26b), (2.26a) and the first equation of (2.24). We find We now insert this expression in (2.37) to get We then get the following third order expansion for Y k

Estimates of the time derivatives.
In this subsection, we prove that the time derivatives of X k (s, τ )/ε, Y k (s, τ ) are uniformly bounded when the initial data is chosen following the Chapman-Enskog procedure presented previously. Note that due to the 1/ε factor for X k , the expansion for X k has to be performed one order further compared to the expansion of Y k . This is stated in the following proposition.
for some constant C 0 > 0 independent of ε.
Proof. The time derivatives of the unknown, denoted in this proof asX k (s, With X k (0, τ ) = X 1st k (τ ) (given by (2.28)-(2.29)) and Y k (0, τ ) = v k,0 and using equation (2.11), we find the following initial data for the previous system Now we fix τ ∈ T as a parameter and together with (2.16), we definẽ We can apply Lemma 2.2 with H = L ∞ τ ([0, 2π]) since all the assumptions of this lemma are fullfilled. Indeed,X k,τ (0) andỸ k,τ (0) are uniformly bounded in ε, the functions X k,τ and Y k,τ enjoy some boundedness properties (thanks to Lemma 2.3) and the functions b and E are smooth. This enables to derive the following estimate We now deal with W 1,∞ τ estimates. This is done by differentiating the above system with respect to τ . The so-obtained system is still of the form of the Lemma 2.2 and we have to check that the initial data remains bounded. Clearly the τ -derivative ofX k (0, τ ) is equal to zero ; concerning the τ -derivative ofỸ k (0, τ ), we have From the definition (2.28)-(2.29) of X 1st ) ends the proof of (i). Let us now prove (ii). We denoteX k := ∂ 2 s X k ,Y k := ∂ 2 s Y k which are solutions of the following system In order to apply Lemma 2.2, we first check that the initial data for this system is uniformly bounded (in ε). To do so, we use the system (2.42)-(2.43) at s = 0 to obtaiň Firstly, we find which leads toX k (0, τ ) = O(ε) thanks to (2.44), so that we deduce 1 ε ∂ τXk (0, τ ) = O(1). Concerning the second term in (2.46), we look atỸ k (0, τ ) This enables to prove thatX k (0, τ ) given by (2.46) is uniformly bounded.
Again, considering new unknown 1 εX k (s, s ε + τ ) and e τ J b kY k (s, s ε + τ ) enables to recast the previous system so that, using the previous estimates, we can use Lemma 2.2 with H = L ∞ τ ([0, 2π]) provided that the initial data 1 εX k (0, τ ) = O(1). Then, we compute 1 εX We focus on the first term Then, we compute the second term 1 Gathering the two term leads to . Similar computations for ∂ τXk and ∂ τYk leads to the required estimate.  s X, ∂ 2 s Y (second and third rows) under first order initial data (X 1st , v 0 ) and respectively second order initial data (X 2nd , Y 1st ).

Numerical illustrations.
To end this section, we illustrate the effect of the preparation of the initial data on the behaviour of the time derivative of X k and Y k . To do so, we consider an example of a single particle (we then omit subscript k) characteristics (2.2) with initial condition In Figures 1, we plot the time history of ∂ s X, ∂ s Y, ∂ 2 s X and ∂ 2 s Y (by accurate numerical solver) under norm for different initial data (first order initial data (X 1st , v 0 ) with X 1st given by (2.29) and second order initial data (X 2nd , Y 1st ) with X 2nd given by (2.36) and Y 1st given by (2.32)). The numerical results confirm the results of Proposition 2.4.

Numerical method
This section is devoted to the construction of numerical schemes for the two-scaled system (2.11). We will perform the analysis of a first order numerical scheme: we will prove that this numerical scheme enjoy the uniform accuracy property with respect to ε. In addition, we will prove that the scheme is able to reproduce the confinement property (2.14) at the discrete level. Then, we will propose a strategy to reach the second order accuracy and to handle the coupling with Poisson equation.
Let ∆t > 0 be the time step and denote t n = n∆t for n ≥ 0 as the descretisation of the t-variable. For each particle 1 ≤ k ≤ N p , the discretisation of the scaled time s-variable is consequently as ∆s k = b k ∆t, s n k = n∆s k , n = 0, 1, . . . .
We certainly omit this subscript k for brevity, i.e. ∆s = ∆s k , s n = s n k . Denote the numerical solution as 3.1. First order numerical scheme. A first order implicit-explicit (IMEX1) finite difference scheme reads for n ≥ 0, In the Fourier space in τ , the above scheme is easily diagonalized. By discretizing the τ -direction as τ j = j∆τ, j = 0, . . . , N τ with ∆τ = 2π/N τ , N τ being some positive even integer, one can use the Fourier transform in τ to get a fully discretized scheme. By doing so, let us remark that the IMEX scheme (3.1) is explicit from a computational point of view and the error in τ is uniformly (with respect to ε) spectrally uniform. By assuming that E ε (t, x) and b(x) are smooth given functions, we analyse the first order IMEX scheme (3.1) for which we have the following uniform convergence results.

2)
and The error estimate shows that the scheme with well-prepared initial data offers super-convergence in x k . As a consequence, this super-convergence is also true for space dependent macroscopic quantities such as ρ ε (t, x). The estimate (3.3) indicates the confinement property at the discrete level.
We are going to prove this theorem by first introducing two lemmas concerning local truncation error and error propagation. To simplify the notations, we will always use C to denote a positive constant independent of ∆s, n or ε and it could change from line to line. We shall omit the subscript k from now on.
Plugging them into (3.4b), we have ∆t∂ t E ε (t n + ∆tθ, X(s n , τ )) + ∆s∇ x E ε (t n+1 , X(s n + ∆sθ, τ ))∂ s X(s n + ∆sθ, τ ) dθ, + ∆s bε ∆t∂ t E ε (t n + ∆tθ, X(s n , τ )) + ∆s∇ x E ε (t n+1 , X(s n + ∆sθ, τ ))∂ s X(s n + ∆sθ, τ ) dθ, Hence again by Proposition 2.4 and Lemma 2.3, we get ξ n We now focus on the propagation error. First, we denote the error function as It is obvious by the choice of the initial data that Taking the difference between the local error (3.4) and scheme (3.1), we get the error equations In the following, for any τ -function ϕ(τ ), its Fourier coefficients are defined by For the IMEX1 scheme (3.1), we have the following formula for the error function on the l ∈ Z Fourier coefficients in τ and for n ≥ 1, Proof. Taking the Fourier transform (3.8) of (3.6), we have Let n = 0 in the above relations, and then inserting (3.10b) into (3.10a), noting that e 0 Therefore, the formula (3.9) with n = 1 is true. Let us assume (3.9) is true for 1 ≤ m ≤ n and we check the case m = n + 1. Plugging First, we use the relation which enables to recover the second term of (3.9a) for n + 1 . Second, we remark than p l α 0,l = 0 so that the first term in the previous expression of (e n+1 X ) l can be reformulated as Taking the square of (3.14) and using Cauchy-Schwarz inequality lead to Then, to get H 1 which completes the induction proof by considering ∆s ≤ ∆s 0 = 1/C.

3.2.
Second order numerical scheme. A second order scheme (IMEX2) could be written down as follows. For n ≥ 0, The above IMEX2 scheme is also explicit. From practical results (as can be seen in the next section), we observe that the IMEX2 with the 3rd order prepared initial data (derived in Appendix A) gives second order uniform accuracy. However, the rigorous error estimates would be more involved than the first order IMEX1 scheme and it is still under-going. We will address it in a future work.

3.3.
A strategy for the Vlasov-Poisson case. When the Vlasov equation (1.1a) is coupled to Poisson equation, the electric field E ε is a self-consistent field and the problem becomes nonlinear. Under PIC discretisation, E ε in (2.2) is given by Hence by under the scale of time s j = b j t, the electric field evaluated at one particular particle x j (t) (1 ≤ j ≤ N p ) solves (as we see in (2.5)) Each particle carries its own frequency and now all the particle are coupled to each other through Poisson equation, which as a result mixes all the frequencies. Thus, (2.5) is a multiple-frequency system with a large number N p of degrees [8] and the proposed two-scale formulation is not rigorously working. Here, we give a practical strategy that works well based on our numerical experiments. Note the above strategy relies on the key confinement property (see Lemma 2.1) as well as the two-scale version in Lemma 2.2. In order to have a better control of C(s) in the oscillatory case, we consider the scale of time (2.4) dynamically. Discretise time with ∆t > 0 and denote t n = n∆t. For n ≥ 0 and t n ≤ t ≤ t n+1 , define and we solveẋ for one step with ∆s = b n k ∆t. Again we isolate the leading order oscillation term, filter out this main oscillation (2.6) and then consider the two-scale formulation but leave the high-frequency character of the electric field part alone [13]. We then obtain (3.19b) For the initial data, we formally choose X k (0, τ ) = X 1st k (τ ), Y k (0, τ ) = Y 1st k (τ ) given by (2.29) and (2.32) by replacing x k,0 and v k,0 with respectively x k (t n ) and v k (t n ). This initial data enables to offer second order uniform accuracy when an exponential integrator scheme (as in [14]) is used for (3.19).
We shall briefly derive the scheme. For the simplicity of notations, we put (3.19a)-(3.19b) into the following compact form: where we denote (3.21b) Applying Fourier transform in τ on (3.20) and then by Duhamel's principle from 0 to ∆s (n ≥ 0), A first order uniformly accurate scheme, shorted as EI1 in the following, is obtained as, ∆s, l = 0. A second order scheme, shorted as EI2, is given as, Suppose now we have computed numerically X k (∆s, τ ) as the two-scale solution for system (3.18), we update the electric field for the next time level t n+1 as x ∈ R 2 .

Numerical results
This section is devoted to numerical illustrations of the numerical schemes introduced above. We consider (1.1) with the following initial data to test convergence order. The spatial domain is x = (x 1 , x 2 ) ∈ Ω = [0, 2π/k] × [0, 2π] for some k, η > 0. We choose η = 0.05, k = 0.5 and discretise Ω with 64 points in x 1 -direction and 32 points in x 2 -direction. As a diagnostic, we consider the following two quantities: We then compute the relative errors of the different numerical schemes with respect to ρ ε and ρ ε v at the final time t f = 1 in maximum space norm. We numerically solve (1.1) with two configurations. For the first one, we consider an external electric field given by which will be addressed as 'given E' in the numerical results. For the second case, we consider the nonlinear Vlasov-Poisson equation (1.1)-(1.3). The reference solution is obtained by using a fourth order Runge-Kutta method on the original problem (2.2) with step size ∆t = 10 −6 .
For the PIC method, we choose N p = 204800 particles and the projection of the particles on the uniform spatial grid is done by quintic splines. The time step ∆s is determined by fixing ∆t (recall the relation s = b k t and ∆s = b k ∆t) so that after N time steps (such that t f = N ∆t, every particles stops at the same time t f . Finally, we denote by N τ the number of points in the τ -direction: ∆τ = 2π/N τ . In the sequel, second order initial data will refer to (X 2nd k , Y 1st k ) with X 2nd k given by (2.36) and Y 1st k given by (2.32), whereas third order initial data will refer to (X 3rd k , Y 2nd k ) with X 3rd k given by (A.1) and Y 2nd k given by (2.41). External electric field. We first study the 'given E' case. In Figure 2, the errors (in L ∞ t norm) in time of the IMEX1 scheme in ρ ε with second order initial data are displayed for different values of ε. As shown in the numerical analysis (see Proposition 2.4 and Theorem 3.2), the scheme IMEX1 has uniform first order accuracy in time. Moreover, we can observe that the error decreases as ε goes to zero, in agreement with theoretical results.
In Figure 3, the errors (in L ∞ t norm) in time of the IMEX1 scheme with second order initial data is plotted regarding the quantities ρ ε /ε and ρ ε v . We can observe the curves are almost superimposed confirming the theoretical error estimates derived previously.
The influence of the discretization in the τ direction of the IMEX1 (using ∆t = 10 −5 ) is presented in Figure 4. We computed the difference between the numerical solution obtained with several N τ and the one using N τ = 64, for the two quantities ρ ε and ρ ε v . The discretization error in τ has a spectral behavior with respect to N τ , and when ε becomes small, we observe that the method reaches machine accuracy with very few number of grids (typically N τ = 8 is sufficient when ε = 10 −3 ).
Next, we study the convergence of the second order IMEX2 scheme with third order initial data. In Figures 5-7, we can see second order uniform accuracy of the scheme with respect to ε, regarding both ρ ε /ε and ρ ε v . We also study the convergence rate of the Vlasov equation (1.1) to the asymptotic model (derived in Appendix B) on the characteristics level when ε → 0. To do so, we measure the difference between the solution of (1.1) for several ε and the one obtained with ε = 10 −4 . In Figure 11, we show the convergence of the model (1.1) in the limit regime (with ∆t = 5 −3 ), for which the rate is equal to one. Finally, we study the dynamics of the solution for a fixed ε = 0.1, aiming to see the effect from the non-constant magnetic field. The quantity ρ ε (t, x) is plotted as a function of x at different times in Figure 8 (with ∆t = 5 −3 ).
Vlasov-Poisson case. For the Vlasov-Poisson case, we apply the strategy presented in subsection 3.3, namely the dynamical scaling EI2 scheme with initial data (X 1st k , Y 1st k ). In Figure 9, we show the convergence results in   Figure 2. Temporal error of IMEX1 with second order initial data (X 2nd k , Y 1st k ) for given E case: relative maximum error in ρ ε . time regarding both ρ ε and ρ ε v . For a long-time diagnostic test of the scheme, we consider the energy of the Vlasov-Poisson equation which is conserved as We compute the numerical energy H(t) by the EI2 scheme with ∆t = 0.05, N τ = 16. In Figure 10 3) to the asymptotic model when ε → 0. We proceed as in the linear case to plot in Figure 11 the convergence of the model (1.1) towards the limit regime. Here again, the rate is close to one.

Conclusion
We proposed a multi-scale numerical scheme for the Vlasov equation with a strong non-homogeneous magnetic field by using a Particle-in-Cell strategy. The solution of the problem is highly oscillatory in time, space and velocity with non-periodic oscillation. Making use of the fact that the positions of the particles in this regime are confined around the initial position, we transformed the characteristics into a suitable form which enabled us to perform the separation of scales techniques. A uniformly accurate first order scheme was then proposed and rigorously analyzed for the Vlasov equation with external electric field; for this scheme, we also proved that it enjoys the confinement property at the discrete level. Practical extensions are performed to achieve the second order accuracy, and also to deal with the case of the Vlasov-Poisson equation. In the later case, it turns out that the characteristics equations are a huge highly oscillatory system with multiple frequencies.
Numerical results are then presented to confirm the theoretical results and illustrate the efficiency of the proposed schemes.

Appendix A. Third order preparation
In this appendix, we derive the third order initial data, which will ensure (following the same strategy used in Proposition 2.4 that the quantities 1 ε ∂ 3 s X k and ∂ 3 s Y k are uniformly bounded. We first derive (2.26a) with respect to s and assuming ∂ 3 s h k = O(ε), we find   Figure 3. Temporal error of IMEX1 with second order initial data (X 2nd k , Y 1st k ) for given E case: relative maximum error in ρ ε /ε (first row) and in ρ ε v (second row).  . Error in τ of IMEX1 with second order initial data (X 2nd k , Y 1st k ) for given E case: relative maximum error in ρ ε (left) and in ρ ε v (right) with respect to the number of grids N τ .
where we used the following notations    Figure 6. Temporal error of IMEX2 with third order initial data (X 3rd k , Y 2nd k ) for given E case: relative maximum error in ρ ε /ε.
Proof. The proof is a recursive process of the known results and it is very similar to that of Proposition 2.4. We omit the details here for brevity.  Figure 11. Convergence the Vlasov equation (1.1) in the given E case (left) and in the Vlasov-Poisson case (right) to in the limit regime: relative maximum error in ρ ε (first row) and ρ ε v (second row) at t f = 1.