NUMERICAL SCHEMES FOR KINETIC EQUATION WITH DIFFUSION LIMIT AND ANOMALOUS TIME SCALE

. In this work, we propose numerical schemes for linear kinetic equation, which are able to deal with a diﬀusion limit and an anomalous time scale of the form ε 2 (1 + | ln( ε ) | ). When the equilibrium distribution function is a heavy-tailed function, it is known that for an appropriate time scale, the mean-free-path limit leads either to diﬀusion or fractional diﬀusion equation, depending on the tail of the equilibrium. This work deals with a critical expo- nent between these two cases, for which an anomalous time scale must be used to ﬁnd a standard diﬀusion limit. Our aim is to develop numerical schemes which work for the diﬀerent regimes, with no restriction on the numerical parameters. Indeed, the degeneracy ε → 0 makes the kinetic equation stiﬀ. From a numerical point of view, it is necessary to construct schemes able to under-take this stiﬀness to avoid the increase of computational cost. In this case, it is crucial to capture numerically the eﬀects of the large velocities of the heavy-tailed equilibrium. Moreover, we prove that the convergence towards the diﬀusion limit happens with two scales, the second being very slow. The schemes we propose are designed to respect this asymptotic behavior. Various numerical tests are performed to illustrate the eﬃciency of our methods in this context.

(Communicated by Axel Klar) Abstract. In this work, we propose numerical schemes for linear kinetic equation, which are able to deal with a diffusion limit and an anomalous time scale of the form ε 2 (1 + |ln(ε)|). When the equilibrium distribution function is a heavy-tailed function, it is known that for an appropriate time scale, the mean-free-path limit leads either to diffusion or fractional diffusion equation, depending on the tail of the equilibrium. This work deals with a critical exponent between these two cases, for which an anomalous time scale must be used to find a standard diffusion limit. Our aim is to develop numerical schemes which work for the different regimes, with no restriction on the numerical parameters. Indeed, the degeneracy ε → 0 makes the kinetic equation stiff. From a numerical point of view, it is necessary to construct schemes able to undertake this stiffness to avoid the increase of computational cost. In this case, it is crucial to capture numerically the effects of the large velocities of the heavy-tailed equilibrium. Moreover, we prove that the convergence towards the diffusion limit happens with two scales, the second being very slow. The schemes we propose are designed to respect this asymptotic behavior. Various numerical tests are performed to illustrate the efficiency of our methods in this context.

1.
Introduction. The modeling of a large amount of particles, from a theoretical or numerical point of view, is a very active field of research. The direct application of Newton's laws leads to a large system of coupled equations, one for each particle of the system. Since such a huge number of unknowns is beyond the reach of numerical computations, the so-called microscopic scale is not adapted to the numerical analysis of these systems. Instead, an approach based on statistical physics is preferred, where the distribution function of particles f ε (t, x, v) depending on the time t ≥ 0, the position x ∈ R d , and the velocity v ∈ R d is considered. The parameter ε ∈ (0, 1] denotes here the Knudsen number, which is proportional to the mean free path of the particle. Provided an initial condition f (0, x, v) = f in (x, v), we consider the following kinetic equation for f ε Θ(ε)∂ t f ε + εv · ∇ x f ε = L(f ε ). (1) The quantity Θ(ε) = ε 2 (1 + ln (ε)), is a suitable scaling parameter to be chosen according to the nature of L, in order to capture a non-trivial dynamic when ε goes to 0. The linear operator L describes the collisions of the particles. We will consider the particular case of the BGK operator with Here, the equilibrium function M is even, positive and normalized to 1 In what follows, the integral in v of f on a domain D will always be denoted f D . The integral is on the whole space R d when no domain is specified. The asymptotic analysis of (1) when ε → 0 can be found in [29]. Our aim in this paper is to construct numerical schemes for (1) that fit with this asymptotic analysis. The properties of the equilibrium function M , and the scaling Θ(ε) of (1) make various asymptotic behaviors arise. For instance, considering the scaling Θ(ε) = ε with a local Maxwellian equilibrium in (1) and the Boltzmann operator instead of the linear collision operator L, leads when ε goes to 0 to a fluid limit for (1) (see [2,1]). With the BGK collision operator, the scaling Θ(ε) = ε 2 , and a global Maxwellian equilibrium M (v) = (2π) −d/2 e −v 2 /2 , the solution f ε of (1) degenerates when ε → 0 to a distribution at equilibrium f = ρM , where ρ solves a diffusion equation ∂ t ρ − ∇ x · (D∇ x ρ) = 0, with initial condition ρ in = f in . A rigorous derivation of diffusion-type equations in this last case was first investigated in [34,3,5,15]. When M is a Maxwellian and Θ(ε) = ε 2 , the diffusion coefficient D is finite because of the exponential decay of M for large v, and is given by However, for non Maxwellian equilibrium, it may happen in many situations that the coefficient D becomes infinite. An important example is the case of a so-called heavy-tailed function M (v) ∼ 1/|v| d+α for large v, α ∈ (0, 2). In this case, the scaling Θ(ε) = ε 2 is not the suitable choice and does not lead to a non trivial dynamics when ε goes to 0. It is well-known that, in this case, the time scale Θ(ε) should be modified according to α in order to capture a non-trivial dynamics, which turns out to be a fractional diffusion model [29,4]. Heavy-tailed equilibria arise in the study of granular media (see [17,7,6]), astrophysical plasmas (see [30,31]), tokamaks (see [16]) or in economy (see [25]). Usually, the fractional diffusion equation describes the motion of the particles driven by a Levy process (see [14,13]), while the classical diffusion is governed by a Brownian process.
From a numerical point of view, the case α ∈ (0, 2) has been treated in [11,10,33]. In this paper, we consider the so-called critical case α = 2, of the above heavy-tailed equilibrium, which induces a different asymptotic behavior when ε goes to 0. For a sake of simplicity we will consider here the following particular case where m is such that M = 1. With this assumption, it has been shown in [29] that if Θ(ε) ∼ ε 2 | ln(ε)| in (1) the solution f ε of (1) converges weakly when ε → 0 to ρ 0 M , where ρ 0 solves a classical diffusion equation Here c d is a coefficient that only depends on the dimension d (c 1 = 2, c 2 = π, and c 3 = 8π/3). To have a non-vanishing of Θ(ε) for ε = 1, we shall consider In this work, we prove that when the initial data f in is at equilibrium, i. e. when f in = ρ in M , the convergence when ε → 0 of the solution f ε of (1) to the solution of (5) is very slow. However, we will show that f ε approaches much faster an intermediate equilibriumρ ε M , whereρ ε is the solution of the following equation (in Fourier variable) and whereρ ε is the space Fourier transform ofρ ε . The main goal of this work is to construct numerical schemes for (1) which do not need to be refined in order to capture the right asymptotic behavior, when ε goes to 0. There are, in fact, some stiffness in the problem (1). First, the smallness of ε imposes severe conditions on the numerical parameters, if a naive approach is used. We will follow an Asymptotic Preserving (AP) strategy [20,23,24] to overcome this difficulty. Namely, if we consider a problem P ε which degenerates into a problem P 0 when ε tends to 0, then an AP scheme S h ε should enjoy the following properties: 1. at fixed ε, S h ε must be consistent with P ε when the discretization parameter h tends to 0. 2. S h ε must degenerate into a scheme S h 0 consistent with P 0 when ε goes to 0. An AP scheme can even enjoy the stronger property of being Uniformly Accurate (UA), meaning that its accuracy does not depend on ε.
The construction of AP schemes for the case of classical diffusion with classical time scale Θ(ε) = ε 2 has been widely investigated (see [20,21,21,22,23,8,9,28,27]), but the heavy-tail of the equilibrium brings an additional stiffness, which is not usual in the classical cases. Indeed, the anomalous time scale (6) in (1) is designed to capture the effect of the high velocities in the asymptotic analysis. From a numerical point of view, taking into account these high velocities is then crucial to ensure the AP property of the scheme. In previous works (see [10,11,12]), we investigated AP schemes for kinetic equations in the fractional diffusion limit. In the critical case, the same methodology cannot rendered word for word to construct AP schemes for this problems. Indeed, the very slow rate of the convergence of (1) towards the diffusion limit (5) makes the asymptotic regime hard to reach numerically. Instead of requiring only the AP property, we then design our schemes to make them respect the asymptotic behavior of the solution, i. e. to converge at first into schemes solving (7) and then into schemes solving (5).
Three numerical schemes will be investigated in this paper. The first one is based on a fully implicit in time scheme for (1) using a Fourier variable in space. A suitable modification is introduced in this case, and its AP behavior is proved. Then, a scheme based on a micro-macro decomposition of the distribution function f ε is presented. It allows to avoid, if needed, the use of the Fourier transform and of the time implicit character of the previous scheme. Moreover, it can be adapted to deal with more general collision operators in (1). Once again, the convergence towards (7) has to be treated with care to ensure the good asymptotic behavior of the scheme. In the end, we propose an approach based on an integral formulation of (1) in Fourier variable, which enjoys the stronger UA property.
To summarize, we will show that (1) approaches (7) with a rate ε 1 + |ln(ε)|, while (7) approaches the limit model (5) with the slower rate 1/ (1 + |ln(ε)|), and we will construct numerical schemes which respect this asymptotic behavior. This can be illustrated by the following diagram The schemes S h ε andS h ε are respectively consistent with the problems P ε andP ε when ε > 0 is fixed, and S h 0 is consistent with P 0 . The scheme S h ε approachesS h ε with rate ε 1 + |ln(ε)| whileS h ε approaches S 0 with rate 1/ (1 + |ln(ε)|) when the discretization parameter h is fixed and when ε goes to 0. The asymptotic behavior P ε →P ε → P 0 will be proved in suitable functional spaces and illustrated by several numerical tests.
The paper is organized as follows. In the next section, we will start by establishing the convergence rates of (1) towards (5), and of (7) to (5). Then, in Section 3, we present an asymptotic preserving scheme for (7), and the three asymptotic preserving schemes for (1), which are tested in Section 4.
2. Degeneracy to the diffusion limit. In this section, we show that the convergence of (1) towards the diffusion limit (5) can be quantified by two steps as described by the diagram (9). This will be the basis of the numerical methods proposed in Section 3. In what follows, the space Fourier transform will be often used. We will denotef (t, k, v) (resp.ρ(t, k)) the space Fourier transform of We will use the following weighted L 2 norm where M is given by (4). We have the following propositions: H N d R d denotes the usual Sobolev space, and let T > 0. Then, there exists a constant C such that whereρ ε ∈ L ∞ 0, T ; L 2 solves (7) with initial condition ρ in , and with M defined in (4).
These two propositions establish that the convergence ε → 0 of (1) makes appear two dynamics with very different rates in ε. Indeed, in a first step, the solution f ε approaches (relatively) quickly a distribution function at equilibriumρ ε (t, x)M (v), whereρ ε is solution of (7). Then, in a second step, this densityρ ε goes much more slowly to the solution ρ 0 of the limit diffusion equation (5). Hence, the approximation of f ε by its diffusion limit (5) is valid only for very small ε. At the numerical level, it implies that a large range of ε may stay out of reach if no appropriate strategy is used in the numerical schemes. The next parts present numerical methods for (1) degenerating when ε tends to 0 into numerical methods solving (7). In particular, the schemes enjoy the AP property that is • when ε is fixed, the scheme is consistent with (1), • the difference between the numerical scheme S h ε for (1) and the numerical schemeS h ε for (7) converges to 0 when ε goes to 0. Afterwards, letting ε go to zero in the schemeS h ε for (7) makes it degenerate in a scheme S h 0 solving the limit diffusion equation (5). This approach enables to capture numerically the two scales of convergence of the solution of (1) to its diffusion limit, as described by (9).
The proofs of Prop. 1 and Prop. 2 follow the ideas developed in [4] in the case of a fractional diffusion limit for kinetic equation, coming from a heavy-tailed equilibrium function. We start by proving the following lemma, which gives the limit of a ε when ε goes to 0.
2. If d = 2, 3, with c 2 = π and c 3 = 8π 3 . Proof of Lemma 2.1. In the 1-dimensional case, with the change of variables w = εv, the coefficient a ε reads a ε (k) = m 1 + |ln(ε)| where the last term in (15) can be computed for nonzero k with the change of variables u = 1/w and the first term is bounded by Ck 2 /(1 + |ln(ε)|). Finally, the inequality ln(x)/2 ≤ 1 + x for x > 0, gives (13). When d = 2, 3, with the change of variables w = ε|k|v for nonzero k, (8) reads where e denotes any unitary vector. The first term is lower than C|k| 2 , it remains to consider the second one. Depending on the dimension, a change of variables in polar or spherical coordinates can be applied in it. Since no additional difficulty arise when d = 3, we treat here the case of polar coordinates, when d = 2. Choosing e such that w · e = |w| cos(θ), it can be rewritten as At this stage, the last two integrals can be bounded and the first one can be computed explicitly In the end, the inequality ln(|k|) ≤ |k| yields (14). Remark 1. The proof of Lemma 2.1 provides, for nonzero k, a rewriting of a ε in which the diffusion limit clearly appears. In dimension 1, it is and in dimension 2 it reads Finally, a change of variables in spherical coordinates provides the case of dimension 3 Proof of Prop. 2. We denote byρ 0 andρ ε the solutions of (5) and (7) in space Fourier variable, with initial conditionρ in , and a 0 (k) = c d mk 2 . They satisfŷ and gives (12).
The convergence of the solution of (7) to the solution of (5) is very slow in ε, but Prop. 1 states that the degeneracy of (1) to (7) happens much faster. The proof of Prop. 1 is based on the proof of the similar result stated in [4] for the case of a kinetic equation with heavy-tailed equilibrium degenerating to a fractional diffusion equation. It uses the following properties of the collision operator, established in [15] and [29] in a more general case: Proposition 3. The collision operator L given by (2) has the following properties Proof of Prop. 1. As in [4], the proof relies on the following Hilbert expansion of the solution f ε of (1) where ρ, g 1 , g 2 and r solve Considering (23) in Fourier variable leads tô Furthermore, the left-hand side of (24) is such that g 2 M − g 2 = 0. Still in Fourier variable, (24) then reads − g 1 + Θ(ε)∂ tρε = 0, that is using (26) and the symmetry of M ∂ tρε (t, k) + a ε (k)ρ ε (t, k) = 0, with a ε defined in (8). Note that solving this equation leads tỗ Asρ ε solves (7), the proof of the proposition relies on establishing estimates for . Since a ε ≥ 0, (27) and the expression ofĝ 1 in (26) give with Θ(ε) defined in (6). The estimation of a ε given by Lemma 2.1 yields with N 1 = 1, and N 2 = N 3 = 2.
The solution g 2 of (24) is unique if g 2 = 0. It can be rewritten as and (26) and (7) imply g 2 = 0. The estimation of r can be done with (25) integrated on the characteristic curves, multiplied by r/M and integrated in velocity. Since g 2 = 0 and the term in r M −r gives a negative integral with the second point of Prop. 3, it reads After an integration in time and space, and using a Cauchy-Schwarz inequality, it writes , so that an upper bound for the right-hand side is needed to complete the proof of the proposition. A differentiation in time of (26), together with (27), gives 3. Numerical schemes. This section presents three numerical schemes designed to approximate the solution of (1), and to preserve the asymptotic behavior of Prop. 1 and Prop. 2 when ε becomes small. These schemes are based respectively on an implicit formulation of (1) in Fourier variable, a micro-macro decomposition, and an integral formulation of (1). These three strategies have been investigated in [11,12] in the case of the fractional diffusion limit of the kinetic equation, α ∈ (0, 2). Here we extend them to the case α = 2. In this critical case, we required that not only the limit equation (5), but the asymptotic behavior (7) of the numerical solution is prescribed.
In what follows, we will denote by t n = n∆t, 0 ≤ n ≤ N , such that N ∆t = T the time discretization, and we will set f n ≈ f (t n ). All the numerical tests will be performed at time T = 0.1, the time step ∆t being specified in each case. The schemes require a velocity discretization, we will denote the quadrature formula for the integration in velocity on the domain D. When no domain is specified, the integration is for all the index 1 ≤ j ≤ N v . In the tests, we will consider a uniform discretization of the domain |v| ≤ 10, with N v points symmetrically distributed on both sides of the origin to ensure that vM Nv = 0. Namely, we will consider N v = 2N v with N v = 100 and use the grid with ∆v = 10/N v . The normalization constant m of the equilibrium M is chosen such that M Nv = 1. The space domain is bounded, and we consider periodic conditions, allowing the use of the Fourier variable. At the discrete level, the discrete Fourier variable will be used. In the tests, we will consider x ∈ [−1, 1], discretized with N x points such that where ∆x = 2/N x . This space discretization is linked to the Fourier modes we use for the computation of the discrete Fourier transform, that are the integers k such that −N x /2 ≤ k ≤ N x /2. At the discrete level, the usual L 2 norm in space, and the norm L 2 M −1 defined in (10) will be computed using the space and velocity discretization Finally, in the following numerical tests, we will take N x = 32 and consider the initial condition The schemes are required to enjoy a stronger property than the AP one. Indeed, they respect the two steps of the convergence of (1) towards (5) following the diagram (9).
The peculiar asymptotic behavior of (1) comes from the effect of the high velocities in the equilibrium function M . Indeed, its heavy-tailed character makes the usual diffusion coefficient (3) be infinite. The anomalous scaling Θ(ε) is chosen to counterbalance the weight of M for large velocities. It is then necessary to take them into account in the numerical computations. One idea would be to consider an adaptive grid for large velocities, to use as much large velocities as needed. However, linking the stiffness ε and the discretization parameters is not in the spirit of AP strategies. We propose a method based on analytical modifications of the terms degenerating into (8) in the schemes. These modifications are done consistently in the schemes, and (8) is carefully discretized to ensure the degeneracy to the limit equation (5). Once it is done, we are able to prove that the schemes respect the asymptotic behavior of (1).
3.1. Discretization of a ε . Since the degeneracy of (1) to the diffusion limit (5) goes through the equation (7), it is necessary to handle the degeneracy of (7) to (5) numerically to design AP schemes for (1). In Fourier variable, both (5) and (7) are differential equations, which can be solved analytically, and that can also be solved with classical schemes. For instance, a semidiscrete-in-time implicit scheme for (7) readsρ and it degenerates into a scheme solving (5) if a ε (k) → a 0 (k) = mc d |k| 2 when ε goes to 0. If the velocity integrations are done continuously, this is ensured thanks to (13). Unfortunately, the discrete velocity integration breaks this property. Indeed, with (29) the discrete version of a ε (k) writes Since the quadrature uses a finite number of points, the sum above is finite whereas its continuous version is infinite. Then, the naive numerical version of a ε (k) tends to 0 when ε → 0, which is not satisfactory. Consequently, it is necessary to modify analytically a ε before applying the discretization (29), to ensure the correct asymptotic behavior of a ε . In Remark 1, the expressions a ε (k) rewritten in (20)-(21)-(22) make the limit behavior of a ε clearly appear. Hence, at the discrete level, we rather use a discrete version of these terms, where the integrals are computed with an usual quadrature (29). It ensures the discrete version a Nv ε of a ε to degenerate to the expected diffusion coefficient when ε tends to 0. To do so, in dimension 1, we use the following discrete version of a ε , from (20) Nv,|v|≤1 which respects the estimates (13) at the discrete level. With a Nv ε , the Euler scheme for (7) degenerates when ε goes to 0, to a Euler scheme for (5).
The numerical tests highlight this behavior. Indeed, Fig. 1 and 2 display the solutions given by the scheme (34), with a ε computed with the discretizations (35) and (36). The solutions are computed with ∆t = 10 −2 , for a range of ε. We also plot the numerical solution of the limit model approximated bŷ When a ε is computed directly with (35), we observe that the solution of (34) does not tend to the limit diffusion when ε goes to 0, but to the initial condition ρ in = f in Nv , with f in defined in (33). Conversely, the solution of (34) with a ε computed with (36), comes close to the solution of (37) when ε goes to 0. Moreover, Fig. 3 suggests that this convergence in ε happens with the speed 1/(1 + |ln(ε)|), proved in Prop. 2. To obtain this figure, we computed the solution ρ ε,∆t of (34) with (36), and the solution ρ ∆t of (37), with ∆t = 10 −2 . Then, we plotted the relative error with · L 2 Nx defined in (32), as a function of ε ∈ (0, 1).   3.2. Implicit scheme (IS). The first idea to write a scheme for (1) which enjoys the AP property is to use a fully implicit scheme for (1) written in Fourier variable. Indeed, it is well-known that in the case of the classical diffusion limit, such a scheme degenerates when ε → 0 into a scheme solving the limit equation. In our case, such an approach does not work, due to the effects of large velocities. Then, it is necessary to suitably write the scheme to make this degeneracy appear. We start with (1) written in space Fourier variable, and we consider a fully implicit time discretization where the dependence in ε of f n+1 , f n and ρ n+1 has been omitted to simplify the notations. This scheme degenerates when ε → 0 intô Therefore, provided that the densityρ n+1 is known, the scheme (38) enjoys the AP property, andf n+1 readsf An expression forρ n+1 is given by a fully implicit time discretization of (1) integrated in velocityρ which writes, using the expression off n+1 in (39) (40) Using the symmetry of M , the last term of the previous expression can be rewritten as Let us remark that the expression (42) converges to the coefficient a ε when ε decreases. Moreover, the second term of (40) vanishes for small ε. As a consequence, if the integrals in velocity are computed exactly, the scheme (40) converges into a scheme solving (7) when ε goes to 0. To ensure the scheme to enjoy the AP property at the fully discrete level, we slightly modify the expression (42) to make the coefficient a ε appear. This modification is done consistently, since the expression (41) indicates that is of order ∆t 2 for fixed ε. Thus, it can be removed or modified into a term of the same order with no incidence on the accuracy of the scheme. We then replace (42) consistently by ∆t Θ(ε) + ∆t whereρ n+1/2 can be taken equal toρ n orρ n+1 , depending on the desired limit scheme (implicit or explicit). Then the coefficient a ε is numerically approximated as in Section 3.1 to ensure the AP property. Eventually, we have the following proposition Proposition 4. We consider the following scheme defined for all k and all time indices with Θ(ε) and a Nv ε defined in (6) and (36), and with initial conditionf Nv . The quantityρ n+1/2 can be chosen equal toρ n or ρ n+1 depending on the desired asymptotic scheme (implicit or explicit in time). This scheme has the following properties: 1. The scheme is of order 1 for fixed ε > 0 and preserves the total mass. 2. The scheme is AP: for a fixed ∆t, the scheme solves the diffusion equation with c 1 = 2, c 2 = π, c 3 = 8π 3 and where m has been defined in (4). Remark 2. The degeneracy of (44) to a scheme solving (5) respects the two dynamics discussed in the diagram (9). Indeed, the numerical tests show that (44) degenerates with speed ε 1 + |ln(ε)| tô ρ n+1 −ρ n ∆t + a Nv ε (k)ρ n+1/2 = 0, which solves (7). Then letting ε going to 0 in the previous scheme makes it degenerate to (45) with speed 1/(1 + |ln(ε)|).
Proof. The first point of Prop. 4 is a direct consequence of the fact that we used an Euler scheme for (1), and that the modifications we applied to it were done consistently. Taking k = 0 in the second line of (44) yields the conservation of the total mass. Finally, with the modification of the term (42), the AP character is straightforward. Indeed, letting ε becomes small with fixed ∆t in (44) leads to   f n+1 =ρ n+1 M ρ n+1 −ρ n ∆t + a Nv ε (k)ρ n+1/2 = 0, which degenerates into a scheme solving (5) when ε → 0.

Micro-macro scheme (MMS).
In the previous part, we proposed an implicit scheme for (1), which has the same behavior as the solution of (1) when ε goes to 0. However, the use of the Fourier variable may be restrictive in some cases, and its extension to more general collision operators seems difficult. Moreover, the implicit treatment of the transport operator induces high computational cost, especially in multi-dimensional cases. We propose here a scheme based on a micromacro decomposition of the distribution function f ε , which respects the behavior of f ε when ε → 0. It treats the transport terms explicitly, and can be extended to the case of general collision operators of [29] with a strategy similar to [26].
Denoting ρ ε = f ε , we consider the decomposition f ε (t, x, v) = ρ ε (t, x)M (v) + g ε (t, x, v), such that g ε = 0. Inserting it in (1) and integrating in velocity gives Then, we multiply (47) by M and we subtract it from (1) to get an equation for g ε A semi-discrete-in-time numerical scheme can be designed following [28], in which the stiffest terms are treated implicitly (49) Here, the dependence in ε of ρ n and g n is omitted to simplify the notations. For fixed ε, it is a consistent scheme, it remains to check wether it preserves the correct asymptotic. The first line of (49) rewrites which implies that g n+1 = O(ε) when ∆t is fixed. After that, we remark that in the right hand side of the equality, all the terms in g n contains a power of ε. Hence, g n+1 = −εv · ∇ x ρ n M + o(ε) for fixed ∆t. When reported in the second line of (49), this gives the following limit scheme for ρ which does not corresponds to the expected diffusion limit. Indeed, the term between brackets is infinite when the integration in velocity is done continuously, but it is finite when the discretization (29) is applied. Hence, in this latter case, when ε → 0, the limit scheme is ρ n+1 = ρ n , which is not consistent with (5). It is then necessary to modify (49) to ensure the correct asymptotic behavior. To do so, we modify the macro equation by making appear the inverse of the transport operator. We first remark that v · ∇ x g n+1 = v · ∇ x f n+1 . Indeed, one can show by induction that g n = 0 at each step of the computations, using (50). The micro-macro decomposition of the density f n is then ρ n M + g n for all n. We then use a semi discrete implicit formulation of f as in the previous part (see (38) for instance) This expression is injected in the flux of the macro equation of (49) to get Since the scheme is first order in time, the bracket in (52) can be simplified consistently in v · ∇ x f n to avoid the inversion of a differential operator. For the same reason, we could remove (53) with no incidence on the consistency of the scheme; however, it degenerates when ε → 0 to the diffusion of (5). Therefore, this term is kept to ensure the AP property of the scheme, but it must be treated with care to take correctly the effects of the high velocities into account. In Fourier variable, and using the evenness of M , (53) reads This term is of order ∆t 2 when ε is fixed, and it degenerates to a ε (k)ρ n+1 (with a ε defined in (8)) when ε becomes small with ∆t fixed. Thus, (53) can be modified consistently in ∆t 2 /(Θ(ε) + ∆t) 2 a ε (k)ρ n+1/2 , which ensures both the AP character of the scheme, and the fact that the numerical solution solves (7) when ε belongs to the intermediate regime in which the solution f ε of (1) solves (7). Here,ρ n+1/2 can be chosen equal toρ n orρ n+1 . The passage from the semi-discrete-in-time scheme to the fully discretized scheme can then be done with the quadrature (29), and the suitable discretization of a ε defined in (8), as in section 3.1. The derivatives in space are treated with a classical first order upwind scheme; for a sake of simplicity in the notations, we still write here the spatial derivatives as continuous. Eventually, we have the following proposition.
Proposition 5. We consider the following scheme (discretized in time and velocity) defined for all x ∈ R d , v j , 1 ≤ j ≤ N v and all time indices 0 ≤ n ≤ N, N ∆t = T > 0 by with Θ(ε) and a Nv ε defined in (6) and in (36), where F −1 denotes the inverse of the Fourier transform, and with initial condition ρ 0 = f in Nv , and g 0 = f in − ρ 0 M . The quantityρ n+1/2 can be chosen equal toρ n orρ n+1 depending on the desired asymptotic scheme (implicit or explicit in time). This scheme has the following properties: 1. The scheme is of first order in time for any fixed ε > 0 and preserves the total mass. 2. The scheme is AP: for a fixed ∆t, the scheme solves the diffusion equation (5) when ε goes to 0ρ with c 1 = 2, c 2 = π, c 3 = 8π 3 , and where m has been defined in (4). (49), for which the integrals with respect to velocities were not discretized, the definition of g n+1 in (55) ensures that g n+1 Nv = 0 for all 0 ≤ n ≤ N − 1.

Remark 3. As in
Remark 4. As in the case of the implicit scheme, this scheme enjoys a property stronger than being AP, since it respects the diagram (9). Hence, it respects the dynamics of the convergence to the solution of the diffusion equation (5) when ε → 0. Indeed, the numerical tests show that (55) degenerates with speed ε 1 + |ln(ε)| to the scheme (46), which solves (7). The degeneracy to the limit scheme (56) happens much more slowly, with speed 1/(1 + |ln(ε)|).
Remark 5. If we want to avoid the use of Fourier variable, we can replace F −1 a Nv ερ n+1/2 by −mc d ∆ x ρ n+1/2 . This modification is consistent with (1), and still enjoys the AP property, but the convergence to the limit scheme does not make a scheme for (7) arise in the intermediate regime in ε.
Remark 6. This micro-macro method can even be adapted to deal with more general linear collision operators. Indeed, the micro-macro decomposition of the kinetic equation (47)-(48) holds, with a slight modification on the right-hand side of (48) which becomes L(g)/Θ(ε), where L denotes the collision operator. The scheme for the micro equation can then be written using a relaxed method, as presented in [26], and the scheme for the macro equation comes very similarly as for the BGK case, see [18]. Indeed, the key point is again to inject an implicit formulation of (1) in (47), and to deal with the velocity integrations to make sure the high velocities are correctly taken into account. Note that this gives an AP scheme for the problem, and that the two scales convergence stated in Section (2) is a priori not valid for more general collision operator, so that it is not relevant to require the scheme to respect it. In the similar case of a fractional diffusion limit for a kinetic equation, micro-macro schemes for kinetic equation with general collision operator have been proposed in [19,32].
Proof. Since the first point comes immediately from the derivation of the scheme, let us prove the second point. The first line of (55) gives that g n+1 = o(1) when ε → 0. Once injected in the second line, it yields which is a consistent approximation of (7). Finally, a Nv ε tends to mc d |k| 2 when ε → 0, which concludes the proof.

Integral formulation based scheme (DS).
In the previous parts, we proposed two numerical schemes solving (1), which in addition respect the asymptotic behavior of the continuous solution when ε goes to 0. They enjoy the AP property, but their precision is not uniform in ε, as it is suggested in Fig. 6 and 10 of Section 4. Indeed, these figures present the densities given by the two previous schemes, for different values of ε, and with fixed discretization parameters. In particular, one can remark that the curves on the top left picture are different. This phenomenon comes from the non-uniform accuracy of the previous schemes, and from the fact that the densities have been computed for ε = 10 −2 , which belongs to the intermediate regime in which AP schemes are known to suffer from accuracy reduction. No refinement on the discretization parameter has been done, since the test that is performed in these figures consists in highlighting the convergence towards the intermediate equation (7) for fixed discretization parameters (AP character). The AP/not-UA character of these schemes can be explained by the treatment applied on the integrals in velocity arising in the coefficient a ε defined in (8). Indeed, it is numerically computed with (36) after manipulations, such as a change of variable depending on ε in its expression, intended to ensure the AP property of the schemes. It preserves the AP property, since there is no need of refining the grids to make the limit scheme appear, but it breaks the UA property. Another reason why the two schemes do not enjoy the UA property is the way the expressions (42) and (54) have been treated to construct the schemes. In both cases, they were modified consistently to make the coefficient a ε (k) defined in (8) appear explicitly in the scheme, and thus to make easier the construction of an AP scheme respecting the diagram (9). However, this modification breaks the UA property. Considering for instance the case of the implicit scheme, the modified expression (43) has the same asymptotic behavior as a ε when ε goes to 0 with fixed ∆t, it corresponds to the AP property. But, if it is considered with ∆t = Θ(ε), it does not have the same asymptotic behavior as a ε when ε goes to 0, which corresponds to a lack of UA property.
In this section, we propose a scheme based on an integral formulation of (1) in Fourier variable. This approach was investigated in [11]- [12] in the case of the fractional diffusion limit. This scheme was shown to be uniformly accurate (UA): its precision does not depend on the stiffness parameter ε. A similar strategy is performed here, to deal with the diffusion case with an anomalous scaling.
We start from the Duhamel form forf ε , solution of (1) in space Fourier variablê Evaluating at time t n+1 and integrating in velocity, we get (58) To write a scheme onρ ε , we perform a quadrature of order 2 in the integral in s; Assuming that the time derivatives ofρ ε are uniformly bounded in ε, the remainder is uniformly bounded by a term of magnitude O ∆t 2 independently of ε. Inserting (59) in (58) yields the following scheme forρ ε where the dependence in ε ofρ ε has been omitted to simplify the notations. The coefficients b j , c j , 0 ≤ j ≤ n are given by The consistency of this scheme comes directly from its derivation, but the discretization of the integrals in velocity appearing in b j and c j is crucial to ensure the AP property. Indeed, as in the previous parts, their continuous limit when ε goes to 0 makes the diffusion (5) However, it is still not sufficient to ensure the AP character of the scheme when the integrals are discretized since, when j = 0, the line (63) of the coefficient c j (62) reads which is equal to Θ(ε)a ε , with a ε defined in (8) and rewritten in (15) (note that the equality (15) is valid for any dimension). In Section 3.1, we established that the discretization of this term is crucial to ensure the AP property of our schemes. Following the strategy presented previously, we then replace the line (63) by an appropriate discretization Θ(ε)a Nv ε defined in (36) when j = 0. Finally, c 0 rewrites The expressions b j (j ≥ 0), c j (j ≥ 1) and c 0 , defined in (61)-(62)-(64), computed with the discretization (29) make the scheme (60) AP. Moreover, let us remark that the scheme only uses ρ and not f . It implies that the computation of the whole distribution f is not necessary in the scheme. Moreover, as they do not depend on the distribution f , but only on the discretization parameters and on ε, the coefficients b j and c j can be precomputed. If several computations are done using the same velocity grid, time step, Fourier modes and ε, the coefficients b j and c j can be computed only once. In this case, as the problem for f is of higher dimension than the problem on ρ, this scheme written for ρ only represents a gain of computational cost, compared to a scheme for which the computation of f is needed. Let us note that f can even so be computed in Fourier variable with a scheme similar to (60). Indeed, Duhamel formula of (1) between the times t n and t n+1 leads tô As in the scheme forρ, we apply the quadrature (59) in the integral in s, and the scheme forf readsf where the dependence in ε off has been omitted, and with Eventually, we have the following proposition.
, we consider the following scheme, defined for all k and all time indices 0 ≤ n ≤ N, N ∆t = T , by with A 0 , β, γ defined in (57)-(66)-(67), and b j (j ≥ 0), c j (j ≥ 1) and c 0 defined in (61)-(62)-(64) and computed with (29). This scheme has the following properties: 1. The scheme is first order in time and preserves the total mass. 2. The scheme is AP: for a fixed ∆t, it solves the diffusion equation (5) when ε goes to zeroρ n+1 −ρ n ∆t + mc d |k| 2ρn+1 = 0, with c 1 = 2, c 2 = π, c 3 = 3π 8 , and where m has been defined in (4). 3. Moreover, the semi-discrete-in-time scheme enjoys the UA property: it is first order uniformly in ε ∃C > 0, sup Remark 7. The UA property of the scheme ensures that it respects the dynamics of the convergence towards the diffusion equation (5).
Proof. The conservation of the mass can be proved by induction on the expression ofρ n+1 for k = 0. Moreover, since the proof of the third point of Prop. 6 is very similar to the one presented in [11]- [12] in the case of the fractional diffusion limit, we just prove the AP character of the scheme (68). Let us just remind that the UA property come from the fact that the quadrature (59) is a convex combination, and from the fact that all the integrals in velocity arising in the scheme are computed with the same quadrature, including the part modified to ensure the AP property. Hence the scheme forρ degenerates when ε goes to 0 intô The asymptotic behavior of a Nv ε has already been studied, it remains to consider the limit of the expression off n+1 when ε goes to 0. From (66) and (67), we have and in the end, the scheme forf n+1 degenerates for small ε intof n+1 =ρ n+1 M.
4. Numerical tests. In this section, we propose numerical tests in dimension 1 to validate the schemes of the previous parts. For simplicity, we will denote IS the implicit scheme of Prop. 4, MMS the micro-macro scheme of Prop. 5, and DS the integral formulation based scheme of Prop. 6. The implicit Euler scheme (46) for the equation (7), with a ε computed with (36), will be denoted quasi-diff. At last, the scheme (45) for the limit diffusion equation (5) will be denoted diff. To highlight the consistency of the schemes, we will compare their solutions to the solution of a reference scheme which is an explicit Euler scheme in Fourier variable for (1) As detailed in the beginning of Section 3, all the tests are performed at time T = 0.1, the time step being specified in each case. We will consider N x points in space on the interval [−1, 1] with periodic boundary conditions, and the velocities are discretized with N v = 199 points symmetrically distributed on the interval [−10, 10]. The initial condition is (33). In the DS scheme, the velocity integrations with the variable w in (61)-(62)-(64) are computed with the same grid as the integrals with the v variable.
4.1. Implicit scheme (IS). In this section, we test the properties of the IS scheme. First of all, its consistency is tested. Fig. 4 shows that, for ε = 1, its solution coincides with the solution of the explicit scheme (69). For this figure, the two solutions are computed with the time step ∆t = 10 −2 . Then, the solution ρ ∆t of the IS scheme for a range of ∆t is compared to the solution ρ ∆t ref =10 −6 of the IS scheme for ∆t ref = 10 −6 at time T = 0.1. The relative error between these densities with · L 2 Nx defined in (32), is displayed in Fig. 5 in function of ∆t. As it is a line with slope 0.98, the numerical order of the method is 1.
Then, the AP character of the IS scheme is tested. The dynamics of the convergence towards the diffusion limit (5) is highlighted in Fig. 6, where the densities obtained with the IS scheme are presented for a range of ε, and compared to the densities given by quasi-diff and diff schemes. These solutions are computed with ∆t = 10 −2 . For intermediate ε, the solution of the IS scheme sticks to the solution of quasi-diff, and the two densities goes together to the solution of diff when ε tends to 0. The convergence towards the solution of quasi-diff happens with speed ε 1 + |ln(ε)|, as suggested by the error study presented in Fig. 7. Denoting by f  the distribution function obtained with the IS scheme, and ρ quasi−dif f the density obtained with the quasi-diff scheme for ∆t = 10 −2 , the quantity

4.2.
Micro-macro scheme (MMS). In this section, we focus on the MMS scheme and we highlight its properties with the same tests as in the previous section. As the transport operator is treated with a classical upwind scheme in the MMS scheme, the CFL condition imposes to take time steps smaller than in the previous section.    Fig.  9, where the consistency error (70) for ε = 1 is displayed in function of ∆t. For this figure, the density ρ ref erence is the density obtained with the MMS scheme for ∆t = 10 −6 and ρ ∆t are densities obtained with the MMS scheme for a range of ∆t. As expected, the numerical order of the method is 1.
Concerning the AP character of the MMS scheme, Fig. 10 shows, for a range of ε, the densities obtained with the MMS, quasi-diff and diff schemes with ∆t = 10 −4 . Once again, the solution of the MMS scheme first joins the solution of the quasi-diff scheme when ε becomes small and the two reach the solution of diff scheme together when ε tends to 0. The speed of the convergence towards the solution of quasi-diff scheme is ε 1 + |ln(ε)| stated in Prop. 1. Indeed, Fig. 11 displays the error (71),  where ρ quasi−dif f , and the solution f given by the MMS scheme are computed with ∆t = 10 −4 .

4.3.
Integral formulation based scheme (DS). In this section, we test the integral formulation scheme DS. For ε = 1, we check Fig. 12, that the solutions of the explicit and DS schemes are close. This figure was obtained with ∆t = 10 −2 .
The order of accuracy of the DS scheme is studied in Fig. 13. For ε = 1, it displays the error (70) between ρ ref erence defined as the solution of the DS scheme for ∆t = 5 · 10 −5 , and the solution of the DS scheme ρ ∆t for a range of ∆t. In a logarithmic scale, the slope of the line obtained is close to 1, as expected. Then, we can study the AP character of the DS scheme. The dynamics of the convergence towards the diffusion limit (5) is highlighted in Fig. 14. It presents the densities obtained with the DS, quasi-diff and diff schemes for ∆t = 10 −2 . The solution of the DS scheme has the right behavior, since it is very close to the solution  To highlight the fact that the DS scheme is of order 1 uniformly in ε, we compare the results given by the DS scheme for ∆t ref = 5 · 10 −5 to the results given by the same DS scheme for different values of ∆t and ε. With these densities, the error (70) is displayed in Fig. 16 as a function of ε. We observe that the error curves are stratified with respect to ε, showing the uniformity of the scheme with respect to ε.    property of the schemes. To deal with this stiffness, we proposed a method based on analytical transformation of the terms leading to the asymptotic model in the schemes. A similar approach was proposed in [11,12] in the case of a fractional diffusion limit for the kinetic equation.
The first scheme we presented is based on a fully implicit discretization of the kinetic equation in Fourier variable, and the second one uses a micro-macro decomposition of the distribution function and treats the transport explicitely. Both of them enjoy the AP property and respect the dynamics of the convergence towards the limit model. The last scheme we proposed, which is based on an integral formulation of the kinetic equation in space Fourier variable, even enjoys the UA property. In a near future, we aim at extending this work to more general context, considering an integral collision operator. In this context, the study can not be simplified by the use of the space Fourier variable. However, schemes based on a micro-macro decomposition of the distribution can be written. Moreover, the case of fractional diffusion limit of kinetic equations with an integral collision operator can also be treated with such an approach.