Nicolas Crouseilles
email: nicolas.crouseilles@inria.fr

É L Ène Hivert ‡ ¶

Mohammed Lemou
email: mohammed.lemou@univ-rennes1.fr

NUMERICAL SCHEMES FOR KINETIC EQUATIONS IN THE DIFFUSION AND ANOMALOUS DIFFUSION LIMITS. PART I: THE CASE OF HEAVY-TAILED EQUILIBRIUM

Keywords: BGK equation, Diffusion limit, Anomalous diffusion equation, Asymptotic preserving scheme AMS subject classifications. 35B25, 41A60, 65L04, 65M22

In this work, we propose some numerical schemes for linear kinetic equations in the diffusion and anomalous diffusion limit. When the equilibrium distribution function is a Maxwellian distribution, it is well known that for an appropriate time scale, the small mean free path limit gives rise to a diffusion type equation. However, when a heavy-tailed distribution is considered, another time scale is required and the small mean free path limit leads to a fractional anomalous diffusion equation. Our aim is to develop numerical schemes for the original kinetic model which works for the different regimes, without being restricted by stability conditions of standard explicit time integrators. First, we propose some numerical schemes for the diffusion asymptotics; then, their extension to the anomalous diffusion limit is studied. In this case, it is crucial to capture the effect of the large velocities of the heavy-tailed equilibrium, so that some important transformations of the schemes derived for the diffusion asymptotics are needed. As a result, we obtain numerical schemes which enjoy the Asymptotic Preserving property in the anomalous diffusion limit, that is: they do not suffer from the restriction on the time step and they degenerate towards the fractional diffusion limit when the mean free path goes to zero. We also numerically investigate the uniform accuracy and construct a class of numerical schemes satisfying this property. Finally, the efficiency of the different numerical schemes is shown through numerical experiments.

1. Introduction. The modeling and numerical simulation of particles systems is a very active field of research. Indeed, they provide the basis for applications in neutron transport, thermal radiation, medical imaging or rarefied gas dynamics. According to the physical context, particles systems can be described at different scales. When the mean free path of the particles (i.e. the crossed distance between two collisions) is large compared to typical macroscopic length, the system is described at a microscopic level by kinetic equations; kinetic equations consider the time evolution of a distribution function which gives the probability of a particle to be at a given state in the six dimensional phase space at a given time. Conversely, when the mean free path is small, a macroscopic description (such as diffusion or fluid equations) can be used. It makes evolve macroscopic quantities which depends only on time and on the three dimensional spatial variable. In some situations, this description can be sufficient and leads to faster numerical simulations.

Mathematically, the passage from kinetic to macroscopic models is performed by asymptotic analysis. From a numerical point of view, considering a small mean free path, the kinetic equation then contains stiff terms which make the numerical simulations very expensive for stability reasons. In fact, a typical example is the presence of multiple spatial and temporal scales which intervene in different positions and at different times. These behaviors make the construction of efficient numerical methods a real challenge.

In this work, we are interested in the time evolution of the distribution function f which depends on the time t ≥ 0, the space variable x ∈ Ω ⊂ R d and the velocity v ∈ R d , with d = 1, 2, 3. Particles undergo the effect of collisions which are modelized here by a linear operator L acting on f through

L(f )(t, x, v) = ρ(t, x)M β (v) -f (t, x, v),
where ρ(t, x) = f (t, x, v) =:

R d f (t, x, v)dv,
and where the equilibrium M β is defined by

M β (v) =      1 (2π) d/2 exp -|v| 2 /2 if β = 2 + d, ∼ |v|→∞ m |v| β if β ∈ (d, d + 2), (1.1)
where m is a normalization factor. In the sequel, we will always denote by brackets the integration over v. In order to capture a nontrivial asymptotic model, a suitable scaling has to be considered. The good scaling of the kinetic equation (see [START_REF] Mellet | Fractional diffusion limit for collisional kinetic equations[END_REF][START_REF] Mellet | Fractional diffusion limit for collisional kinetic equations: A moments method[END_REF][START_REF] Ben Abdallah | Fractional diffusion limit for collisional kinetic equations : A Hilbert expansion approach[END_REF]) satisfied by the distribution function f is given by

ε α ∂ t f + εv • ∇ x f = L(f ), with α = β -d ∈ (0, 2], (1.2)
where ε > 0 is the Knudsen number (ratio between the mean free path and a typical macroscopic length). The case α = 2 refers to the classical diffusion limit (the equilibrium corresponds to the first equation in (1.1)) whereas the case α ∈ (0, 2) refers to the anomalous diffusion limit (the equilibrium corresponds to the second equation in (1.1)). Although our analysis could be applied to general operators L, we will consider for a sake of simplicity a BGK type operator. We suppose that M β given by (1.1) is an even positive function such that M β (v) = 1 and vM β (v) = 0. Equation (1.2) has to be supplemented with an initial condition f (0, x, v) = f 0 (x, v) and spatial periodic boundary conditions are considered.

The main goal of this work is to construct numerical schemes for (1.2) which are uniformly stable along the transition from kinetic to macroscopic regime (i.e. for all ε > 0). More precisely, two asymptotics are considered here: the diffusion limit (α = 2) and the anomalous diffusion limit (α ∈ (0, 2)). For the diffusion limit, several numerical schemes have already been proposed but in the anomalous diffusion limit, no numerical strategy has been proposed in the literature for (1.2), to the best of author's knowledge.

Mathematically, the derivation of diffusion type equation from kinetic equations such as (1.2) when α = 2 was first investigated in [START_REF] Wigner | Nuclear reactor theory[END_REF], [START_REF] Bensoussan | Boundary layers and homogenization of transport processes[END_REF], [START_REF] Larsen | Asymptotic solution of neutron transport problems for small mean free paths[END_REF] or [START_REF] Degond | Diffusion limit for non homogeneous and non-microreversible processes[END_REF]. When M β decreases quickly enough for large |v|, the solution f of (1.2) converges, when ε goes to zero, to an equilibrium function f (t, x, v) = ρ(t, x)M β (v) where ρ(t, x) is solution of the diffusion equation ∂ t ρ(t, x) -∇ x • (D∇ x ρ(t, x)) = 0, (1.3) and where the diffusion matrix D is given, in the simple case of the BGK operator, by the formula

D = R d v ⊗ vM β (v)dv. (1.4)
In particular, a crucial assumption on M β for D to be finite is that

R d |v| 2 M β (v)dv < +∞;
this is the case when M β is Maxwellian (first case of (1.1)). Hence, if the equilibrium has no finite second order moment, the asymptotic of rescaled kinetic equation with α = 2 is not able to capture a nontrivial dynamics and the value of α should be tuned with M β in order to recover this macroscopic equation when ε → 0. Typically, let M β be a heavy-tailed distribution function, corresponding to the second case of (1.1). An example of such a distribution function we will often consider is

M β (v) = m 1+|v| β with m chosen such that M β (v) = 1.
In the case of astrophysical plasmas, it often occurs that the equilibrium M β is a heavy-tailed distribution of particles which typically generates an anomalous diffusion behaviour (see [START_REF] Mendis | Cosmic dusty plasma[END_REF][START_REF] Summers | The modified plasma dispersion function[END_REF]). In particular, the diffusion matrix D is no longer well-defined and this requires to adapt the value of α in terms of the tail decreasing rate β, in order to capture a nontrivial dynamics. It has been shown in [START_REF] Mellet | Fractional diffusion limit for collisional kinetic equations[END_REF][START_REF] Mellet | Fractional diffusion limit for collisional kinetic equations: A moments method[END_REF] and [START_REF] Ben Abdallah | Fractional diffusion limit for collisional kinetic equations : A Hilbert expansion approach[END_REF] that the appropriate scaling is α ∈ (0, 2) depends on the size β of the tail of the equilibrium through the relation β = α + d.

When ε goes to 0, the solution of this equation with BGK operator converges to ρ(t, x)M β (v) where ρ is the solution of the anomalous diffusion equation which can be written in Fourier variable

∂ t ρ(t, k) = -κ|k| α ρ(t, k), (1.5)
where ρ stands for the space Fourier variable of ρ, k is the Fourier variable and κ is a constant which depends only on the size of the tail of the equilibrium M β (v) and is defined by

κ = R d (w • e) 2 1 + (w • e) 2 m |w| β dw, (1.6)
for any e ∈ R d such that |e| = 1 (note that κ does not depend on e). The anomalous diffusion equation can also be written in the space variable

∂ t ρ(t, x) = - mΓ(α + 1) c d,α (-∆ x ) α 2 ρ(t, x), (1.7)
where m is the normalization factor of M β in (second case in 1.1), Γ is the usual function defined by

Γ(x) = +∞ 0 t x-1 e -t dt,
and c d,α is a normalization constant given by

c d,α = αΓ d+α 2 2π d 2 +α Γ 1 -α 2 .
(1.8)

The fractional operator of the anomalous diffusion is easily defined by its Fourier transform

(-∆ x ) α 2 ρ (k) = |k| α ρ(k),
but has also an integral definition

(-∆ x ) α 2 ρ(x) = c d,α P.V. R d ρ(x + y) -ρ(x) |y| d+α dy,
where P.V. denotes the principal value of the integral. The anomalous diffusion occurs not only in astrophysical plasmas but also in the study of granular media (see [START_REF] Ernst | Scaling solutions of inelastic Boltzmann equations with overpopulated high energy tails[END_REF][START_REF] Bobylev | Boltzmann equations for mixtures of maxwell gases: Exact solutions and power like tails[END_REF][START_REF] Bobylev | On some properties of kinetic and hydrodynamic equations for inelastic interactions[END_REF]), tokamaks (see [START_REF] Del Castillo-Negrete | Non-diffusive transport in plasma turbulence: A fractional diffusion approach[END_REF]) or even in economy or social science: the Pareto distribution is a power tail distribution that satisfies the second case of (1.1). This kind of distribution is a common object used to modelise the repartition of sizes in a given set (asterods, cities, income, opinions. . . ), see [START_REF] Kleiber | Statistical size distributions in economics and actuarial sciences[END_REF] for detailed explanations. At a microscopic level, these velocity distributions arise when the motion of the particles is no longer governed by a brownian process but by a Levy process (see [START_REF] Moor | Fractional diffusion limit for a stochastic kinetic equation[END_REF], [START_REF] Moor | Limites diffusives pour des équations cinétiques stochastiques[END_REF]).

In this work, we want to construct a numerical scheme over the kinetic equation using a heavy-tailed distribution equilibrium, which is robust when the stiffness parameter ε goes to 0. Asymptotic Preserving (AP) schemes already exist in the case of the diffusion limit where the equilibrium is a Maxwellian (see [START_REF] Jin | Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations[END_REF][START_REF] Buet | Asymptotic preserving scheme and numerical methods for radiative hydrodynamic models[END_REF][START_REF] Carrillo | Numerical schemes of diffusion asymptotics and moment closures for kinetic equations[END_REF][START_REF] Jin | Discretization of the multiscale semiconductor Boltzmann equation by diffusive relaxation schemes[END_REF][START_REF] Klar | An asymptotic-induced scheme for nonstationary transport equations in the diffusive limit[END_REF][START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF][START_REF] Jin | Uniformly accurate diffusive relaxation schemes for multiscale transport equations[END_REF][START_REF] Naldi | Numerical schemes for kinetic equations in diffusive regimes[END_REF]) but the strategy used to obtain them cannot be rendered word for word to the anomalous diffusion case. As in the diffusion case, the difficulties when numerically solving the kinetic equation in the case of anomalous diffusion come from the stiff terms which require a severe condition on the numerical parameters. But in the anomalous diffusion case, there is an additional difficulty since one needs to take into account the large velocities arising in the equilibrium M β . Indeed, the role of these high velocities is crucial and they have to be captured numerically in order to produce the anomalous diffusion operator when ε → 0.

To derive a numerical scheme which enjoys the AP property, different strategies are investigated in this work. The first one is based on a fully implicit scheme in time (both the transport and the collision term are considered implicit). Using this brute force strategy, the so-obtained numerical scheme is obviously AP in the diffusion limit but not in the anomalous diffusion one. Therefore, a suitable transformation of this implicit scheme is proposed in this work and leads to a scheme enjoying the AP property in the anomalous diffusion limit.

From a computational point of view, a fully implicit scheme may be very expensive especially if one deals with high-dimensional problems. To avoid this implicitation, we propose another strategy which is based on a micro-macro decomposition. The distribution function f is written as a sum of an equilibrium part ρM β plus a remainder g. A model equivalent to (1.2) is then derived, composed of a kinetic equation satisfied by g and a macroscopic equation satisfied by ρ. The macroscopic flux is then reformulated to obtain a numerical scheme which enjoys the AP property and which does not require the implicitation of the transport term.

The last strategy proposed in this work is based on a Duhamel formulation of the equation from which a AP numerical scheme is obtained; this approach bears similarities with [START_REF] Hochbruck | Exponential integrators[END_REF]. This numerical scheme is efficient for the diffusion case and has to be adapted also to get an AP scheme in anomalous diffusion limit. Moreover, there are two main advantages in using such strategy. First we derive a closed equation on ρ and the computation of the full distribution function is not needed at any time. Second, this strategy ensures that the numerical scheme has uniform accuracy in ε. Eventually, it would be easier to construct high order (in time) schemes on this formulation, this will be done in [START_REF] Crouseilles | Uniformly accurate numerical scheme for kinetic equation in the diffusion and anomalous diffusion limit[END_REF].

The paper is organized as follows. In the next section we shall start with formal computations recalling how the diffusion and of anomalous diffusion limits can be obtained from the kinetic equation. Section 3 is devoted to the derivation of the different numerical schemes: a fully implicit scheme, a micro-macro scheme and a scheme based on a Duhamed formulation of (1.2). After the theoretical study of their properties, several numerical tests are presented in Section 4 in a one-dimensional case with periodic boundary conditions in space to illustrate and compare the efficiency of the various numerical schemes.

2.

A closed equation for the density and its diffusion asymptotics. This section is devoted to the derivation from (1.2) of a closed equation satisfied by ρ. This equation will be the basic equation from which we formally derive the diffusion and the anomalous diffusion equations. Of course, the formal derivation of diffusion and anomalous diffusion equations is classical but, for a sake of clarity, we present one way (among many others) to perform such derivation. In particular, the formal computation below will drive the construction of a class of our numerical schemes. We recall here the hypothesis on the equilibrium M β (v) introduced in (1.1): it is an even function such that M β (v) = 1, vM β (v) = 0. Then we will consider the two following cases mentioned in (1.1): α = 2 and 0

< α < 2 Case 1. v ⊗ vM β (v) < +∞ and α = 2. Case 2. v ⊗ vM β (v) = ∞ and α ∈ (0, 2).
For example, M β has a heavy tail of the form

M β (v) ∼ |v|→∞ m |v| β , β ∈ (d, d + 2), (2.1)
with a positive normalization factor m ensuring M β = 1.

2.1. Space Fourier transform based computations. Starting from (1.2), we propose here a method based on spatial Fourier transform to formally derive the asymptotic equation in both diffusion and anomalous diffusion limits. It will also be the basis of the numerical methods we set in Section 3.3.

Eventually, we remind that we denote by f (t, k, v) (resp. ρ(t, k))

f (t, k, v) = R d e -ik•v f (t, x, v)dx,
the space Fourier transform of f (t, x, v) (resp. ρ(t, x)). We have Proposition 2.1.

(i) Equation (1.2) formally implies the following exact equation on ρ

ρ(t, k) = e -t ε α (1+iεk•v) f (0, k, v) + t ε α 0 e -s e -iεsk•v M β (v) ρ(t -ε α s, k)ds. (2.2) (ii)
In case 1 (diffusion scaling), when ε → 0, ρ solves the diffusion equation

∂ t ρ(t, k) = -(k • v) 2 M β (v) ρ(t, k).
(iii) In case 2 (anomalous diffusion scaling), when ε → 0, ρ solves the anomalous diffusion equation

∂ t ρ(t, k) = -κ |k| α ρ(t, k),
with κ defined by (1.6) and β -d = α.

Proof.

(i) Formal derivation of the expression (2.2) from (1.2) We start by taking the Fourier transform of (1.2)

∂ t f + iε 1-α v • k f = 1 ε α ρM β -f .
We solve in f (assuming ρ is given)

f (t, k, v) = e -t ε α (1+iεk•v) f (0, k, v) + 1 ε α t 0 e -t-s ε α (1+iεk•v) ρ(s, k)M β (v)ds,
which can be written after a change of variables s

→ (t -s)/ε α f (t, k, v) = e -t ε α (1+iεk•v) f0 + t ε α 0 e -s(1+iεk•v) ρ(t -ε α s, k)M β (v)ds.
Still denoting by brackets the integration over the velocity space, we can integrate the previous expression with respect to v to get (2.2). Then, we can write it to make the limit equation appear

ρ = Â0 (t, k) + t ε α 0 e -s e -iεsk•v M β dsρ + ε α t ε α 0 se -s ρ(t -ε α s, k) -ρ(t, k) ε α s e -iεsk•v M β (v) ds, (2.3)
where

Â0 (t, k) = e -t ε α (1+iεk•v) f0 , (2.4)
denotes the initial layer term. As we will see in the next parts of the proof, this formulation gives us the two asymptotic equations of the proposition. (ii) The diffusion limit Considering the diffusion scaling summarized in case 1, we will show that the limit equation is the diffusion equation. In the second integral of (2.3), we expand the expression of ρ in a power series of ε up to the second order. Note that the first term of the expression is exponentially small. Indeed for all t > 0, Â0 = o(ε ∞ ) when ε goes to zero. Then we can write, using a Taylor expansion at ε = 0 for the exponential

ρ = t ε 2 0 e -s M β dsρ -iε t ε 2 0 se -s k • vM β dsρ - ε 2 2 t ε 2 0 s 2 e -s (k • v) 2 M β dsρ + ε 2 t ε 2 0 se -s ρ(t -ε 2 s, k) -ρ(t, k) ε 2 s M β ds + o(ε 2 ).
Reminding that M β = 1 and that M β is even, we can simplify this expression to obtain

ρ = ρ -ε 2 (k • v) 2 M β ρ -ε 2 ∂ t ρ + o(ε 2 ),
where we assume that ρ is sufficiently smooth in time to write

ρ(t -ε 2 s, k) -ρ(t, k) ε 2 s = ε→0 -∂ t ρ(t, k) + o(1).
Finally, when ε goes to 0 we get the expected diffusion equation

∂ t ρ = -(k • v) 2 M β ρ.
(iii) The anomalous diffusion limit Now we will show that when M β fulfills the conditions of case 2 (anomalous diffusion scaling), the solution of (1.2) converges to the solution of the anomalous diffusion equation when ε goes to zero. As in the case of the diffusion limit, the first term has exponential decay when ε goes to zero, Â0 = o(ε ∞ ), for all t > 0. With the same assumption of smoothness in time for ρ, we get the following behaviour for the third term of (2.3)

t ε α 0 se -s ρ(t -ε α s, k) -ρ(t, k) ε α s e -iεsk•v M β ds ∼ ε→0 -∂ t ρ(t, k) + o(1)
.

We now have to make appear the fractional diffusion in the second term of (2.3).

Since the moment of order 2 of M β is not finite, we cannot expand the exponential term into a power series in ε anymore. We decompose this term as

t ε α 0 e -s e -iεsk•v M β ds = t ε α 0 e -s ds + t ε α 0 e -s e -iεsk•v -1 M β ds,
and we focus on the second term on the right hand side of this last equality.

In order to simplify the computations, we will consider the following specific form for the equilibrium

M β : M β (v) = m 1+|v| β , where m is chosen to ensure that M β (v) = 1.
The general case, dealing with M β satisfying (2.1) can be done similarly, see supplementary materials for details. We consider the integral e -iεsk•v -1 M β and, for nonzero k, we perform the change of variables w = ε|k|v

e -iεsk•v -1 m 1 + |v| β = (ε|k|) β-d R d e -is k |k| •w -1 m (ε|k|) β + |w| β dw.
Since the last integral has rotational symmetry, for any unitary vector e of R d we have

e -iεsk•v -1 m 1 + |v| β = (ε|k|) β-d R d e -ise•w -1 m (ε|k|) β + |w| β dw (2.5) = (ε|k|) β-d C(s) + (ε|k|s) β-d R(ε, s, k), where
C(s) = R d e -ise•w -1 m |w| β dw, (2.6)
and

R(ε, s, k) = -m (εs|k|) β R d e -ie•w -1 1 |w| β (εs|k|) β + |w| β dw.
In particular, R(ε, s, k) tends to 0 when ε go to zero and is bounded. If we set R(ε, s, 0) = 0, equality (2.5) is also true for k = 0. From (2.3) and (2.5), we have to set α = β -d to get the fractional diffusion equation when ε → 0 and then we get

ρ = ε→0 ρ + ε α |k| α +∞ 0 w∈R d e -isw•e -1 m |w| β dw e -s dsρ -ε α ∂ t ρ + o(ε α ), ∀t > 0. From the equality +∞ 0 w∈R d e -isw•e -1 m |w| β dw e -s ds = - w∈R d (w • e) 2 1 + (w • e) 2 m |w| d+α dw := -κ,
we derive, when ε goes to 0, the fractional diffusion equation (1.5) given in [START_REF] Mellet | Fractional diffusion limit for collisional kinetic equations[END_REF][START_REF] Ben Abdallah | Fractional diffusion limit for collisional kinetic equations : A Hilbert expansion approach[END_REF].

2.2. Computations in the space variable. The aim of this section is to generalize the previous computations to the original space variable. We have the following proposition Proposition 2.2.

(i) Equation (1.2) formally implies the following exact equation on ρ

ρ(t, x) = e -t ε α f 0 (x -ε 1-α tv, v) + t ε α 0 e -s ρ(t -ε α s, x -εsv)M β (v) ds. (2.7) (ii)
In case 1 (diffusion scaling), when ε → 0, ρ solves the diffusion equation

∂ t ρ(t, x) = ∇ x • (D∇ x ρ(t, x)) , with D = v ⊗ vM β (v) .
(iii) In case 2 (anomalous diffusion scaling), when ε → 0, ρ solves the anomalous diffusion equation

∂ t ρ(t, x) = - mΓ(α + 1) c d,α (-∆ x ) α 2 ρ(t, x),
where Γ is the usual Euler Gamma function, m the normalization constant appearing in the expression of the equilibrium M β , c d,α is defined by (1.8) and β -d = α.

Proof.

(i) Formal derivation of the expression (2.7) from (1.2) As in the previous section, we can integrate (1.2) to get

ρ = A 0 + t ε α 0 e -s ρ(t -ε α s, x -εsv)M β ds, with A 0 (t, x) = e -t ε α f 0 (x -ε 1-α tv, v) . (2.8)
(ii) The diffusion limit Since ρ(t, x -εsv) = ρ(t, x + εsv) , we can rewrite (2.7) in order to make the limit equation appear

ρ(t, x) = A 0 (t, x) + ε 2 t ε 2 0 se -s ρ(t -ε 2 s, x -εsv) -ρ(t, x -εsv) ε 2 s M β (v) ds + ε 2 2 t ε 2 0 s 2 e -s ρ(t, x -εvs) + ρ(t, x + εvs) -2ρ(t, x) ε 2 s 2 M β (v) ds + t ε 2 0 e -s dsρ(t, x).
(2.9)

First, we again notice that the initial layer term A 0 defined in (2.8) has exponential decay when ε goes to zero. Second, Taylor expansions in the two integrals in time enable to write for ε sufficiently small

ρ = ε 2 t ε α 0 se -s (-∂ t ρ) M β ds + ε 2 2 t ε α 0 s 2 e -s t vH x ρvM β ds + 1 -e -t ε 2 ρ + o(ε 2 ),
where H x ρ denotes the Hessian of ρ. This last expression can be simplified into

∂ t ρ = t vH x ρvM β (v) , when ε goes to 0. With D defined in (1.4), we have t vH x ρvM β (v) = ∇ x • (D∇ x ρ)
so this last equation is the diffusion equation as expected.

(iii) The anomalous diffusion limit We rewrite (2.7) in the case of the anomalous diffusion scaling as follows

ρ(t, x) = A 0 (t, x) + ε α t ε α 0 se -s ρ (t -ε α s, x -εvs) -ρ (t, x -εvs) ε α s M β (v) ds + t ε α 0 e -s (ρ(t, x -εvs) -ρ(t, x)) M β (v) ds + t ε α 0 e -s ds M β (v) ρ(t, x).
(2.10)

Under smoothness assumptions on ρ, the first integral of (2.10) becomes, for small ε

t ε α 0 se -s ρ (t -ε α s, x -εvs) -ρ (t, x -εvs) ε α s M β (v) ds = -∂ t ρ(t, x) + o(1).
Since the last term of (2.10) can be explicitly computed, let us focus on the second integral of (2.10). We perform the change of variables y = x -εvs in the integral over v and exchanging the integrals in t and y to get

t ε α 0 e -s (ρ(t, x -εvs) -ρ(t, x)) M β (v) ds = R d ρ(t, y) -ρ(t, x) |x -y| β a(ε, x -y)dy, (2.11)
where

a(ε, z) = t ε α 0 |z| β e -s (εs) d M β z εs ds.
(2.12)

Now, we want to find an equivalent of a(ε, z) for small values of ε. Here we consider the equilibrium M β (v) = m 1+|v| β , the general case with M β satisfying (2.1) can be done similarly (see supplementary materials for more details). The integral a becomes

a(ε, z) = m t ε α 0 |z| β (εs) β (εs) β + |z| β e -s (εs) d ds,
and we consider, for a nonzero z, the following quantity

a(ε, z) -mε β-d Γ(β -d + 1) = -m +∞ 0 (εs) β-d e -s (εs) β (εs) β + |z| β ds -m +∞ t ε α |z| β (εs) β (εs) β + |z| β e -s (εs) d ds. Then, for almost every z ∈ R d , 1 ε β-d a(ε, z) -mε β-d Γ(β -d + 1
) converges to 0 when ε goes to zero. We also have the following estimation

a(ε, z) -mε β-d Γ(β -d + 1) ≤ K ε β-d + ε β-d e -t 2ε α .
(2.13)

Let us remark that for a reasonable regularity of ρ, the integral appearing in (2.11) is well defined as a principal value because β ∈ (d, d + 2)

P.V. R d ρ(t, y) -ρ(t, x) |y -x| β dy = R d ρ(t, y) -ρ(t, x) -(y -x) • ∇ x ρ(t, x) |y -x| β dy.
Hence, using (2.13) we obtain that for almost all y ∈ R d

1 ε β-d ρ(t, y) -ρ(t, x) -(y -x) • ∇ x ρ(t, x) |y -x| β a(ε, z) -mε β-d Γ(β -d + 1) ,
is dominated by an integrable function. We also obtained that for almost all y ∈ R d it tends to zero when ε goes to zero, by using the dominated convergence theorem for this function.

We have, from (2.11)

t ε α 0 e -s (ρ(t, x -εvs) -ρ(t, x)) M β (v) ds = R d ρ(t, y) -ρ(t, x) -(y -x) • ∇ x ρ(t, x) |x -y| β a(ε, x -y)dy,
where we used that a is even. So we have

t ε α 0 e -s (ρ(t, x -εvs) -ρ(t, x)) M β (v) ds ∼ ε→0 mε β-d Γ(β -d + 1) R d ρ(t, y) -ρ(t, x) -(y -x)∇ x ρ(t, x) |y -x| β dy ∼ ε→0 mε β-d Γ(β -d + 1)P.V. R d ρ(t, y) -ρ(t, x) |y -x| β dy.
We remark that we have to set β -d = α to recover the limit equation, so we get

ρ(t, x) = ε→0 -ε α ∂ t ρ(t, x) + mε α Γ(α + 1)P.V. R d ρ(t, y) -ρ(t, x) |y -x| α+d dy + ρ(t, x) + o (ε α ) ,
that is when ε goes to 0

∂ t ρ(t, x) = -mΓ(α + 1)P.V. R d ρ(t, x) -ρ(t, y) |y -x| α+d dy,
which is the expected anomalous diffusion equation.

3. Numerical schemes. This section is devoted to the presentation of appropriate numerical schemes to capture the solution of (1.2). We want these schemes to be Asymptotic Preserving (AP), that is, for arbitrary initial condition f 0 1. to be consistent with the kinetic equation (1.2) when the time step ∆t tends to 0, with a fixed ε, 2. to degenerate into a scheme solving the asymptotic equation (diffusion or anomalous diffusion) when ε goes to 0 with a fixed time step ∆t. In the sequel, three different numerical schemes are presented: a fully implicit scheme, a micro-macro decomposition based scheme and a Duhamel based scheme. We will consider a time discretization t n = n∆t, n = 0, . . . , N such that N ∆t = T where T is the final time and we will set

f n ≃ f (t n ).
As we used Fourier transform in the previous formal computations, we will consider a bounded domain Ω for the spatial domain with periodic conditions. Hence, we will be able to use the discrete Fourier transform in the algorithm.

3.1. Implicit scheme. The first idea to design a scheme for the kinetic equation (1.2) is to set an implicit scheme over the Fourier formulation of the kinetic equation. We will see that in the case of the diffusion limit it preserves easily the asymptotic equation. But in the case of the anomalous diffusion limit, the effect of large velocities must be taken into account to obtain the good asymptotic and this needs a suitable modification of this fully implicit scheme.

We start with (1.2) written in the spatial Fourier variable and consider a fully implicit time discretization

f n+1 -f n ∆t + 1 ε α (1 + iεk • v) f n+1 = 1 ε α ρn+1 M β .
We have

f n+1 = 1 -λ 1 + iλεk • v f n + λ 1 + iλεk • v ρn+1 M β , (3.1)
with

λ = ∆t ε α + ∆t . (3.2) Note that 0 < λ < 1, λ -→ ∆t→0 0 and that λ → ε→0 1.
To compute f n+1 from (3.1), ρn+1 has to be determined first. To do that we integrate (3.1) with respect to v to get

ρn+1 = λ M β 1 + iλεk • v ρn+1 + (1 -λ) f n 1 + iλεk • v , (3.3)
that is

ρn+1 = f n 1 + iλεk • v 1 M β 1+iλεk•v + 1 1-λ iλεk•vM β 1+iλεk•v . (3.4)
We have the following proposition Proposition 3.1. In the case 1 (diffusion scaling), we consider the scheme defined for all k and for all time index 0 ≤ n ≤ N, N ∆t = T by

       ρn+1 (k) = f n 1 + iλεk • v M β 1 + iλεk • v + 1 1 -λ iλεk • vM β 1 + iλεk • v f n+1 = 1 1 + iλεk • v (1 -λ) f n + λρ n+1 M β , (3.5)
with λ = ∆t/(ε 2 + ∆t) and with the initial condition f 0 (k, v) = f0 (k, v). This scheme has the following properties: (i) The scheme is of order 1 for any fixed ε and preserves the total mass ∀n ∈ [[1, N ]], ρn (0) = ρ0 (0).

(ii) The scheme is AP: for a fixed ∆t, the scheme solves the diffusion equation when ε goes to zero

ρn+1 (k) -ρn (k) ∆t = -(k • v) 2 M β (v) ρn+1 (k). (3.6)
Proof. As (i) is straightforward, let us prove (ii). In the case of the diffusion asymptotic, we must check that the scheme is AP. We rewrite (3.4) to get

(1 -λ) ρn+1 = -λ iλεk • vM β 1 + iλεk • v ρn+1 + (1 -λ) f n 1 + iλεk • v ,
that is, with the equality

1 ∆t λ 1-λ = 1 ε 2 , ρn+1 -ρn ∆t = - 1 ε 2 iελk • v 1 + iλεk • v M β ρn+1 - 1 ∆t iλεk • v 1 + iλεk • v f n . (3.7)
Hence we can write

ρn+1 -ρn ∆t = - (k • v) 2 λ 2 1 + λ 2 ε 2 (k • v) 2 M β ρn+1 - 1 ∆t iλεk • v 1 + iλεk • v f n ,
and so, when ε → 0 with a fixed ∆t, the scheme degenerates into

ρn+1 -ρn ∆t = ε→0 -(k • v) 2 M β (v) ρn+1 ,
that is an implicit scheme for (1.3) in Fourier variable. Now we consider the case of an anomalous diffusion scaling 0 < α < 2. In this case, the equilibrium M β is heavy-tailed and D = v ⊗ vM β (v) = +∞. However, we have

lim ε→0 1 ε α ε 2 (k • v) 2 λ 2 1 + λ 2 ε 2 (k • v) 2 M β (v) ρn+1 = |k| α κ,
with κ defined in (1.6). Unfortunately a discretization in velocity of the integral involved in this limit does not see the heavy tail of the equilibrium, which means that the effect of large velocities is completely missed. In fact, the limit when ε goes to zero of the implicit kinetic scheme we wrote in Prop. 3.1 degenerates into the following scheme

ρn+1 = ρn , f n+1 = ρn+1 M β ,
which is not the correct asymptotics. Hence we have to transform its expression to make the right limit clearly appear in the scheme. The scheme is presented in the following proposition. Proposition 3.2. In the case 2 (anomalous diffusion scaling), we consider the following scheme defined for all k and for all time index 0 ≤ n ≤ N, N ∆t = T by

                   ρn+1 (k) = f n 1 + iλεk • v M β 1 + iλεk • v + (ελ|k|) α 1 -λ m (ελ|k|) α+d + |v| α+d (v • e) 2 1 + (v • e) 2 f n+1 = 1 1 + iλεk • v (1 -λ) f n + λρ n+1 M β , with λ = ∆t/(ε α + ∆t)
, where e is any unitary vector and with the initial condition f 0 (k, v) = f0 (k, v). This scheme has the following properties: (i) The scheme is of order 1 for any fixed ε and preserves the total mass

∀n ∈ [[1, N ]], ρn (0) = ρ0 (0).
(ii) The scheme is AP: for a fixed ∆t, the scheme solves the anomalous diffusion equation when ε goes to zero

ρn+1 (k) -ρn (k) ∆t = -κ|k| α ρn+1 (k), (3.8)
where κ has been defined by (1.6). Proof. As (i) is straightforward, let us prove (ii). We start by considering the expression of ρn+1 in the scheme of Prop. 3.1, and more precisely the integral iελk

• v M β 1 + iελk • v = R d M β (v) ε 2 λ 2 (k • v) 2 1 + ε 2 λ 2 (k • v) 2 dv.
As in the previous section, we perform the change of variables w = ελ|k|v to obtain

R d M β (v) ε 2 λ 2 (k • v) 2 1 + ε 2 λ 2 (k • v) 2 dv = (ελ|k|) -d R d M β w ελ|k| w • k |k| 2 1 + w • k |k| 2 dw, (3.9) = (ελ|k|) -d R d M β w ελ|k| (w • e) 2 1 + (w • e) 2 dw
, with e = k/|k|; note that the last integral does not depend on this unitary vector e thanks to its rotational invariance. We consider explicitly the case M β (v) = m 1+|v| α+d with m such that M β = 1 and with α ∈]0, 2[. We have

M β ε 2 λ 2 (k • v) 2 1 + ε 2 λ 2 (k • v) 2 = (ελ|k|) α R d m (ελ|k|) α+d + |w| α+d (w • e) 2 1 + (w • e) 2 dw
, where e is any unitary vector. Inserting this formula in (3.5), we get the expression for ρn+1 . The equation satisfied by f n+1 in Prop. 3.1 is unchanged, then the numerical scheme of the proposition is derived. Now, it remains to show that this scheme enjoys the AP property. From the following form of the equation on ρn+1

ρn+1 -ρn ∆t = -λ α |k| α m (ελ|k|) α+d + |v| α+d (v • e) 2 1 + (v • e) 2 ρn+1 - iλεk • v 1 + iλεk • v f n ∆t , (3.10
) we easily observe that it degenerates when ε goes to zero to an implicit scheme for the anomalous diffusion equation. This concludes the proof.

Scheme based on a micro-macro decomposition.

In the previous part, we constructed a fully implicit scheme enjoying the AP property. The implicit character of the transport may induce a high computational cost. Therefore, we propose another scheme, which is based on a micro-macro decomposition of the kinetic equation and in which the transport part is explicit in time. In the diffusion case, such a scheme has been set in [START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF]. We first recall here the way to derive it in the diffusion case. Then, we consider the case of the anomalous diffusion and show how to take into account the effect of the large velocities induced by the heavy-tailed structure of the equilibrium M β (v).

The distribution f is decomposed as

f = f M β + g = ρM β + g,
where g = 0. Injecting this decomposition into (1.2), we have

∂ t ρ M β + ∂ t g + ε ε α v • ∇ x ρ M β + ε ε α v • ∇ x g = - 1 ε α g. (3.11)
To derive an equation on ρ, we integrate with respect to v to obtain

∂ t ρ + ε ε α v • ∇ x g = 0.
(3.12)

Now we replace ∂ t ρ in (3.11) to get an equation satisfied by g

∂ t g + ε ε α v • ∇ x ρ M β + ε ε α (v • ∇ x g -v • ∇ x g M β ) = - 1 ε α g. (3.13)
From this formulation the semi-implicit scheme proposed in [START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF] writes

     ρ n+1 -ρ n ∆t + ε ε α v • ∇ x g n+1 = 0 g n+1 -g n ∆t + ε ε α v • ∇ x ρ n M β + ε ε α (v • ∇ x g n -v • ∇ x g n M β ) = - 1 ε α g n+1 , (3.14)
and we have the following proposition Proposition 3.3. In case 1 (diffusion scaling), we consider the following scheme defined for all x ∈ Ω, v ∈ V and all time index 0 ≤ n ≤ N by

     ρ n+1 -ρ n ∆t + 1 ε v • ∇ x g n+1 = 0 g n+1 -g n ∆t + 1 ε v • ∇ x ρ n M β + 1 ε (v • ∇ x g n -v • ∇ x g n M β ) = - 1 ε 2 g n+1 ,
with initial conditions ρ 0 (x) = ρ(0, x) and g 0 (x, v) = f (0, x, v) -ρ(0, x)M β (v). This scheme has the following properties:

(i) The scheme is of order 1 for any fixed ε and preserves the total mass ∀n ∈ [[1, N ]], ρn (0) = ρ0 (0).

(ii) The scheme is AP: for a fixed ∆t, the scheme solves the diffusion equation when ε goes to zero

ρn+1 (k) -ρn (k) ∆t = -(k • v) 2 M β (v) ρn (k).
For the proof, we refer to [START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF]. We now consider the anomalous diffusion scaling. Formulation (3.14) does not work because the fractional laplacian does not appear when ε goes to 0, so we have to modify the micro-macro scheme. As M β is even we have

v • ∇ x g n+1 = v • ∇ x ρ n+1 M β + v • ∇ x g n+1 = v • ∇ x f n+1 ,
and so the first equation of (3.14) can be rewritten as

ρ n+1 -ρ n ∆t + ε ε α v • ∇ x f n+1 = 0.
Now we use an implicit formulation of the kinetic equation to express f n+1

ε α f n+1 -f n ∆t + εv • ∇ x f n+1 = ρ n+1 M β -f n+1 .
As we did in the previous part we introduce the variable λ defined by (3.2) and write this last equation as

(I + ελv • ∇ x ) f n+1 = λρ n+1 M β + (1 -λ)f n ,
that is

f n+1 = λ (I + ελv • ∇ x ) -1 ρ n+1 M β + (1 -λ) (I + ελv • ∇ x ) -1 f n .
As λ = ∆t→0 O(∆t), we will use the following approximated expression for f n+1

f n+1 = λ (I + ελv • ∇ x ) -1 ρ n+1 M β + (1 -λ) (I -ελv • ∇ x ) f n + O(∆t),
to avoid a costly inversion of the transport operator. In this expression, the terms we removed also vanish when ε goes to zero. The scheme presented in Prop. 3.3 can be rewritten in the anomalous diffusion limit as

           ρ n+1 -ρ n ∆t + ε ε α λ v • ∇ x (I + ελv • ∇ x ) -1 ρ n+1 M β + ε ε α (1 -λ) v • ∇ x (I -ελv • ∇ x ) f n = 0, g n+1 -g n ∆t + ε ε α v • ∇ x ρ n M β + ε ε α (v • ∇ x g n -v • ∇ x g n M β ) = -1 ε α g n+1 .
(3.15) We can make additional simplifications. First, we remark that,

v • ∇ x (I -ελv • ∇ x ) f n = v • ∇ x f n + O(∆t) = v • ∇ x g n + O(∆t),
where we used vM β = 0 and λ = O(∆t) at fixed ε. Then, the last term in the first line of the scheme (3.15) can be replaced by ε 1-α (1 -λ) v • ∇ x g n . Second, we will reformulate the second term of the first line of the scheme (3.15) (which produces the fractional diffusion when ε goes to zero) in order to get the right limit when ε → 0. To simplify the presentation, we present the reformulation using the Fourier variable (because of the non-local character of the fractional laplacian in the space coordinates, the case in these coordinates is more delicate). The second term of (3.15) writes in the Fourier variable

ε ε α λ ik • v 1 + iελk • v M β ρn+1 (k) = 1 ε α ε 2 λ 2 (k • v) 2 1 + ε 2 λ 2 (k • v) 2 M β ρn+1 (k).
Then, as in (3.9) we make the change of variables w = ελ|k|v in the velocity integration to the following semi-discrete scheme (with

M β (v) = m 1+|v| β )            ρ n+1 -ρ n ∆t +λ α F -1 |k| α (v • e) 2 1 + (v • e) 2 m (ελ|k|) β + |v| β ρn+1 (k) + ε ε α + ∆t v • ∇ x g n = 0, g n+1 -g n ∆t + ε ε α v • ∇ x ρ n M β + ε ε α (v • ∇ x g n -v • ∇ x g n M β ) = - 1 ε α g n+1 , (3.16
) where e is any unitary vector and where we denoted by F -1 the inverse of the Fourier transform. The following proposition summarizes the main properties of this micromacro scheme: Proposition 3.4. In case 2 (anomalous diffusion scaling), we introduce the following micro-macro the scheme defined for all x ∈ Ω, v ∈ R d and all time index

0 ≤ n ≤ N, N ∆t = T by            ρ n+1 -ρ n ∆t +λ α F -1 |k| α (v • e) 2 1 + (v • e) 2 m (ελ|k|) β + |v| β ρn+1 (k) + ε ε α + ∆t v • ∇ x g n = 0 g n+1 -g n ∆t + ε ε α v • ∇ x ρ n M β + ε ε α (v • ∇ x g n -v • ∇ x g n M β ) = - 1 ε α g n+1 ,
where e is any unitary vector, and where initial conditions are ρ 0 (x) = ρ(0, x) and g 0 (x, v) = f (0, x, v) -ρ(0, x)M β . This scheme has the following properties: (i) The scheme is of order 1 for any fixed ε and preserves the total mass

∀n ∈ [[1, N ]], ρn (0) = ρ0 (0).
(ii) The scheme is AP: for a fixed ∆t, the scheme solves the diffusion equation when ε goes to zero

ρn+1 (k) -ρn (k) ∆t = -κ|k| α ρn+1 (k),
where κ is defined by (1.6).

Proof. The proof of (i) is immediate, let us prove (ii). We have when ε goes to zero with a fixed ∆t, the scheme solves the anomalous diffusion equation. The equation on g gives, when ε goes to zero.

g n+1 = O(min(ε, ε α )).
The equation on ρ gives, in the Fourier variable

ρn+1 -ρn ∆t + λ α |k| α (v • e) 2 1 + (v • e) 2 m (ελ|k|) β + |v| β ρn+1 = O(min(ε 2 , ε 1+α )).
Thanks to the following relations,

λ α = ∆t ε α + ∆t α = 1 + O(ε α ), and m (ελ|k|) β + |v| β = m |v| β + O(ε α+d ),
we get the implicit discretization of the anomalous diffusion equation (3.8) when ε → Hence the scheme solves the anomalous diffusion equation when ε goes to zero with a fixed ∆t.

Scheme based on an integral formulation of the equation.

In the previous parts, we wrote two AP schemes solving (1.2) in the cases of the diffusion and anomalous diffusion limits, both of them were of order 1 in time. Here we present a scheme based on a Duhamel formulation of (1.2) which has uniform accuracy in ε. This approach bears similarities with the UGKS scheme (see [START_REF] Xu | A gas-kinetic BGK scheme for the Navier-Stokes equations and its connection with artificial dissipation and Godunov method[END_REF][START_REF] Mieussens | On the asymptotic preserving property of the unified gas kinetic scheme for the diffusion limit of linear kinetic models[END_REF]). As in the previous parts, the case of the classical diffusion gives easily an AP scheme but the large velocities effects require a specific treatment for the anomalous diffusion case.

For both diffusion and anomalous diffusion, we start from (2.2)

ρ(t, k) = e -t ε α (1+iεk•v) ρ(0, k, v) + t ε α 0 e -s e -iεsk•v M β (v) ρ(t -ε α s, k)ds, hence evaluating at time t = t n+1 leads to ρ(t n+1 , k) = Â0 (t n+1 , k) + n j=0 t j+1 ε α t j ε α e -s e -iεsk•v M β ρ(t n+1 -ε α s, k)ds, (3.17)
where Â0 is defined by (2.4). We perform a quadrature of order 2 in the integrals.

Assuming that the time derivatives of ρ are uniformly bounded in ε, we have: (3.18) uniformly in ε, with a(ε, s) = 1 -ε α s-tj ∆t . Inserting (3.18) in the integral term of (3.17) leads to

∀j ∈ [[1, N -1]], ∀s ∈ tj ε α , tj+1 ε α , ρ(t n+1 -ε α s, k) = a(ε, s)ρ(t n+1 -t j , k)+(1-a(ε, s))ρ(t n+1 -t j+1 , k)+O(∆t 2 ),
t j+1 ε α t j ε α e -s e -iεsk•v M β ρ(t n+1 -ε α s, k)ds = c j (k)ρ(t n+1 -t j , k) + b j (k)ρ(t n+1 -t j+1 , k) + O ∆t 2 , (3.19)
uniformly in ε, where we used the following notations ∀j ∈

[[0, N ]] b j (k) = t j+1 ε α t j ε α ε α s -t j ∆t e -s e -iεsk•v M β ds, (3.20) c j (k) = t j+1 ε α t j ε α 1 - ε α s -t j ∆t e -s e -iεsk•v M β ds. (3.21)
3.3.1. AP scheme for diffusion. We consider the case of the diffusion limit and we use the previous quadrature to write an AP scheme for the kinetic equation (1.2). We use the quadrature (3.19) in (3.17) to write

ρ(t n+1 , k) = Â0 (t n+1 , k) + n j=0 c j ρ(t n+1-j , k) + b j ρ(t n-j , k) + O ∆t 2 , (3.22)
and as n∆t ≤ T then n j=0 ∆t 2 = O(∆t) and we get a first order scheme that writes

ρn+1 (k) = Â0 (t n+1 , k) + n j=1 c j ρn+1-j (k) + b j ρn-j (k) + b 0 ρn (k) 1 -c 0 .
To ensure the AP property, the time integration in the last two terms are computed exactly to get

b j (k) =   ε α ∆t e -t j ε α (1+εik•v) 1 -e -∆t ε α (1+εik•v) (1 + εik • v) 2 - e -t j+1 ε α (1+εik•v) 1 + εik • v   M β (v) , (3.23)
c j (k) =   e -t j ε α (1+εik•v) 1 + εik • v - ε α ∆t e -t j ε α (1+εik•v) 1 -e -∆t ε α (1+εik•v) (1 + εik • v) 2   M β (v) . (3.24)
We have the following proposition: Proposition 3.5. In case 1 (diffusion scaling), we consider the following the scheme defined for all k and for all time index 0 ≤ n ≤ N such that N ∆t = T by

ρn+1 (k) = Â0 (t n+1 , k) + n j=1 c j ρn+1-j (k) + b j ρn-j (k) + b 0 ρn (k) 1 -c 0 ,
with the initial condition ρ0 (k) = ρ(0, k) and Â0 , b j and c j defined in (2.4)-(3.23)- (3.24). This scheme has the following properties: (i) The scheme is of order 1 and the total mass is conserved

∀n ∈ [[1, N ]], ρn (0) = ρ0 (0).
(ii) The scheme is AP: for a fixed ∆t, the scheme solves the diffusion equation when ε goes to zero

ρn+1 (k) -ρn (k) ∆t = -(k • v) 2 M β (v) ρn+1 (k).
Remark 1. From the numerical test, the scheme appears to be of order one uniformly in epsilon, as highlighted in Fig. 4.7. This uniform property can be proved; however, since the proof is rather technical, this is postponed to a work in preparation (see [START_REF] Crouseilles | Uniformly accurate numerical scheme for kinetic equation in the diffusion and anomalous diffusion limit[END_REF]).

Proof. From the construction of the scheme developed above, this scheme is clearly consistent and of order 1 in time, uniformly in ε. By induction we also remark that it conserves the numerical mass ρn (0). Let us prove that, for a fixed ∆t the scheme solves the diffusion equation when ε goes to zero. Indeed, for small ε the terms b j (k), c j (k) for j ≥ 1 and Â0 are exponentially small, hence the scheme becomes

(1 -c 0 (k))ρ n+1 (k) = b 0 (k)ρ n (k) + o(ε ∞ ), when ε → 0.
(k) = ε 2 ∆t + o(ε 2 ), c 0 (k) = 1 - ε 2 ∆t -ε 2 (k • v) 2 M β (v) + o(ε 2 ).
Hence when ε goes to zero with a fixed ∆t, the scheme of Prop. 3.5 degenerates into

ρn+1 (k) -ρn (k) ∆t = -(k • v) 2 M β (v) ρn+1 (k),
that is an implicit scheme solving the diffusion equation (1.3) in Fourier variable.

3.3.2. AP scheme for anomalous diffusion. Now, we want to set an AP scheme for the kinetic equation (1.2) in the anomalous diffusion scaling. Still using the notations defined in (2.4)-(3.20)-(3.21), the scheme in Prop. 3.5 writes as

(1 -c 0 )ρ n+1 (k) = Â0 (t n+1 , k) + n j=1 c j ρn+1-j (k) + b j ρn-j (k) + b 0 ρn (k).
We remind that in the case of classical diffusion, the term c 0 ρn+1 produces the right limit for small ε. If we consider the same algorithm, as we are forced to consider a bounded velocity space, the numerical second order moment of M β will be large but finite. Therefore, we have to suitably transform the integrals c j and b j before discretizing in velocity. We have

b j = t j+1 ε α t j ε α ε α s -t j ∆t e -s e -iεsk•v -1 M β ds + t j+1 ε α t j ε α ε α s -t j ∆t e -s ds, c j = t j+1 ε α t j ε α 1 - ε α s -t j ∆t e -s e -iεsk•v -1 M β ds + t j+1 ε α t j ε α 1 - ε α s -t j ∆t e -s ds,
and in the velocity integrations we set w = ε|k|v as in (2.5), so when

M β (v) = m/ 1 + |v| β we have b j = ε α |k| α m (ε|k|) β + |w| β t j+1 ε α t j ε α ε α s -t j ∆t e -s e -isw•e -1 ds + t j+1 ε α t j ε α ε α s -t j ∆t e -s ds, c j = ε α |k| α m (ε|k|) β + |w| β t j+1 ε α t j ε α 1 - ε α s -t j ∆t e -s e -isw•e -1 ds + t j+1 ε α t j ε α 1 - ε α s -t j ∆t e -s ds,
with e any unitary vector. To ensure the AP property, the time integrations are computed exactly to get

b j = e -t j ε α 1 -e -∆t ε α ∆t ε α -e -∆t ε α 1 -ε α |k| α m (ε|k|) β + |v| β (3.25) +ε α |k| α - e -t j+1 ε α (1+iv•e) 1 + iv • e + ε α ∆t e -t j ε α (1+iv•e) 1 -e -∆t ε α (1+iv•e) (1 + iv • e) 2 m (ε|k|) β + |v| β , c j = e -t j ε α 1 - 1 -e -∆t ε α ∆t ε α 1 -ε α |k| α m (ε|k|) β + |v| β (3.26) +ε α |k| α e -t j ε α (1+iv•e) 1 + iv • e - ε α ∆t e -t j ε α (1+iv•e) 1 -e -∆t ε α (1+iv•e) (1 + iv • e) 2 m (ε|k|) β + |v| β .
We have the following proposition Proposition 3.6. In case 2 (anomalous diffusion scaling), with the notations Â0 , b j and c j defined in (2.4)-(3.25)-(3.26), we consider the following scheme defined for all k and for all time index 0 ≤ n ≤ N such that N ∆t = T by

ρn+1 (k) = Â0 (t n+1 , k) + n j=1 c j ρn+1-j (k) + b j ρn-j (k) + b 0 ρn (k) 1 -c 0 ,
with the initial condition ρ0 (k) = ρ(0, k). This scheme has the following properties: (i) The scheme is of order 1 and preserves the total mass

∀n ∈ [[1, N ]], ρn (0) = ρ0 (0).
(ii) The scheme is AP: for a fixed ∆t, the scheme solves the anomalous diffusion equation when ε goes to zero

ρn+1 (k) -ρn (k) ∆t = -κ|k| α ρn+1 (k),
where κ has been defined by (1.6).

Remark 2. As in the case of diffusion, the numerical tests (see Fig. 4.10) suggests that this scheme is order one uniformly in epsilon. This uniform property can also be proved and this will be done in a future work (see [START_REF] Crouseilles | Uniformly accurate numerical scheme for kinetic equation in the diffusion and anomalous diffusion limit[END_REF]).

Proof. As for the case of classical diffusion, the conservation of the mass is obtained by induction and the fact that it is of order 1 comes from the Taylor expansion we performed in the integrals. Let us prove that for a fixed ∆t the scheme solves the anomalous diffusion equation when ε goes to zero. From the exponential decay of A 0 and b j , c j , j ≥ 1 we have when ε goes to zero

(1 -c 0 ) ρn+1 (k) = b 0 (k)ρ n (k) + o(ε ∞ ).
On the one side we have

b 0 (k) = ε α ∆t + o(ε α ),
and on the other side we have

c 0 = ε α |k| α 1 1 + iv • e -1 m (ε|k|) β + |v| β + 1 - ε α ∆t + o(ε α ) = -ε α |k| α m |v| β (v • e) 2 1 + (v • e) 2 + 1 - ε α ∆t + o(ε α ) = -ε α |k| α κ + 1 - ε α ∆t + o(ε α ).
We deduce that when ε goes to zero, the scheme degenerates into

ρn+1 (k) -ρn (k) ∆t = -κ|k| α ρn+1 (k),
where κ is defined by (1.6). This is an implicit scheme for the anomalous diffusion equation (1.5). Remark 3. The previous schemes of Prop. 3.5-3.6 can be modified into a numerical scheme which degenerates into a second order time approximation of the asymptotic model. Indeed, we can modify (3.22) to get the following scheme

ρn+1 (k) = A 0 (t n+1 , k) + b 0 ρn (k) + n j=1 c j ρn+1-j (k) + b j ρn-j (k) + (c 0 + b 0 -1) ρn+1 (k) + ρn (k) 2 + (1 -b 0 )ρ(t n+1 , k).
By construction, this scheme degenerates when ε goes to zero, into a Crank-Nicolson numerical scheme for the diffusion or anomalous diffusion equations.

on Duhamel Formulation in the Prop. 3.5 and 3.6. Finally, DS and ADS denote the numerical schemes (3.6) and (3.8) of the asymptotic models (diffusion and anomalous diffusion equations).

In this part, we will consider t > 0, x ∈ [-1, 1], v ∈ R and the following initial data

f 0 (x, v) = (1 + sin(πx))M β (v),
and periodic boundary conditions are imposed in x. In the diffusion case, we will consider the equilibrium

M β (v) = 1 √ 2π e -v 2 /2 ,
whereas in the anomalous diffusion case, we will take

M β (v) = m 1 + |v| β ,
with m chosen such that R M β (v)dv = 1. In the sequel, unless other discretizations are mentioned, the following numerical parameters will hold:

(i) time: the final time will be set T = 0.1 and ∆t = 10 -3 .

(ii) space: we will consider a uniform mesh in x considering N x = 64 points.

(iii) velocity: in both diffusion and anomalous diffusion case, a uniform velocity grid of the truncated domain [-v max , v max ] is considered, with N v points. In the diffusion case, we will take v max = 10 and N v = 20. In the anomalous diffusion case, we will consider v max = 50 and N v = 200. The discretization ensures at the discrete level that R vM β (v)dv = 0. In the tests for the anomalous diffusion, we will consider α = 1.5.

In the sequel, we will be interested in the relative error. The relative error is defined by the difference between a reference solution (obtained by a scheme using a refined mesh) and the solution obtained by the tested scheme through the following formula

Error(T ) = ρ reference (T ) -ρ scheme (T ) 2 ρ reference (T ) 2 , (4.1)
where • 2 denotes the discrete L 2 norm and the reference scheme will be precised in each case.

4.1. The implicit scheme in the case of the anomalous diffusion limit. In this section, we test the ISA scheme. We check that it is of order 1 for a fixed value of ε and the AP property when ε goes to 0 for a fixed value of ∆t.

Firstly, we remark that the ISD scheme used for the diffusion case does not preserve the asymptotic of anomalous diffusion: the left hand side of Fig. 4.1 shows that for ε = 10 -6 the implicit scheme returns a result very different from the result given by (3.8). Indeed, when ε < < 1, the ISD scheme produces a stationary solution. So it appears that the change of variables in the velocity integrals is necessary to capture the anomalous diffusion limit, as shown in the right hand side of Fig. 4.1.

The left hand side of Fig. 4.2 shows the error between the reference scheme (defined as the ADS scheme) and the ISA scheme as a function of ε, for α = 0.8, 1, 1.5. We observe that the convergence of the kinetic equation to the anomalous diffusion when ε goes to zero arises with a speed α. Theorically, the convergence to the anomalous diffusion equation arises with a speed ε min(α,2-α) , but as we consider a finite domain for v in the tests, the convergence rate always appears to be ε α . This result will be proved in [START_REF] Crouseilles | Uniformly accurate numerical scheme for kinetic equation in the diffusion and anomalous diffusion limit[END_REF]. Left: For ∆t = 10 -3 and α = 0.8, 1, 1.5, the error as a function of ε between the reference ADS scheme and the ISA scheme (log scale). Right: The error between the ISA scheme for a range of ∆t = ε α and the ADS scheme as a function of ∆t, which is fixed equal to ε α (log scale).

However the convergence in time of the ISA scheme to the DSA scheme is not uniform with respect to ε. Setting ∆t = ε α in the scheme and letting ∆t go to 0 shows that the densities does not converge to the density given by the anomalous diffusion equation. On the right hand side of Fig. 4.2, we plot the error in time between the ADS scheme as reference and the ISA scheme computed for a different time steps satisfying ∆t = ε α . We observe that the error does not converge to 0 when ∆t goes to 0, illustrating the lack of uniformity.

4.2. The micro-macro scheme. As we saw in Section 3.2, the MMSD scheme is known to be efficient in the diffusion case (see [START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF]). In this part, we focus on the MMSA scheme in the case of the kinetic equation with anomalous diffusion limit. As it does not use only the Fourier variable, we must consider a grid for the space variable. Here we will take x ∈ [-1, 1] discretized with N x = 32 points. As we solve the transport of the micro part (equation on g) in an explicit way (using finite differences), a CFL condition has to be imposed; then the time step is chosen small enough.

We start by testing the consistency of the scheme. We fix α = 1.5 and for a decreasing range of time steps, we compute the solution given by the MMSA scheme. Then, for ε = 1, 0.5, we compare it to the reference solution given by MMSD with ∆t = 10 -6 . For ε = 10 -7 , 10 -8 , we compare to the solution given by the ADS scheme. 4.3.1. Case of diffusion limit. In the case of diffusion limit, we compare the results obtained by the DSD scheme to the results obtained with the ISD scheme. The first thing we check is the convergence of the algorithm when ∆t goes to zero. We choose a range of ε and we compare the densities obtained with the two algorithms as a function of ∆t. Fig. 4.5 displays the errors associated to these comparisons. It appears that for ε ∼ 1, the scheme is of order 2 and of order 1 for small ε. In order to show the AP character of the scheme, we compute the densities given by the scheme DSD for a decreasing range of ε (∆t being fixed) and we check that they converge to the solution of the diffusion equation as ε goes to 0. The densities obtained are presented in Fig. 4.6. To highlight the fact that the scheme is first order in time uniformly in ε, we compare the results given by the ISD scheme with ∆t = 10 -8 to the results given by the DSD scheme for a range of ∆t. These errors are displayed in the left hand side of Fig. 4.7 as a function of ε where we observe that the error curves are stratified with respect to ∆t, showing the first order uniform accuracy of the scheme. This property appears also in the right hand side of Fig. 4.7 where we plotted the error between the DSD scheme computed for ∆t = ε 2 and the limit DS scheme computed for ∆t = 10 -7 . This error tends to zero when ∆t and ε go to zero, thanks to the uniformity of the DSD scheme. 4.3.2. Case of anomalous diffusion limit. In this section we test the DSA scheme. In the case of anomalous diffusion limit, we compare the results given by the

Dt=ε2

Error Fig. 4.7. Left: The difference between the ISD scheme computed for ∆t = 10 -8 and the DSD scheme computed for a range of ∆t as a function of ε (log scale). Right: The difference between the DSD scheme computed for ∆t = ε α and the DS scheme computed for ∆t = 10 -7 as a function of ∆t (log scale).

DSA scheme to the results given by the DSA scheme computed with a smaller time step, to make the same properties appear. We start by highlighting the consistency of the scheme. Considering different values of ε, we compute a reference solution with the scheme DSA with ∆t = 10 -5 and we compare it to the results given by the DSA scheme for some larger time steps. Fig. 4.8 show error study of the convergence of the densities obtained. We observe that for ε = 1, 0.5 the scheme seems to be of order 2 whereas it is of order 1 for smaller values of ε. Then, we can study the AP character of the DSA scheme. The time step being fixed to ∆t = 10 -3 , we check that the results given by the DSA scheme converge to the result given by the ADS scheme when ε goes to 0. Fig. 4.9 presents the densities ρ(t = 0.1, x) obtained by these two schemes for different values of ε and then, the error as a function of ε is plotted for α = 0.8, 1, 1.5. As previously, the expected numerical speed of the convergence to the anomalous diffusion is recovered.

To highlight the fact that the scheme is first order in time uniformly in ε, we compare the results given by the DSA scheme with ∆t = 10 -5 to the results given by the same DSA scheme for a range of ∆t. These errors are displayed in the left hand side os Fig. 4.10 as a function of ε where we observe that the error curves are stratified with respect to ∆t, showing the uniformity of the scheme with respect to ε. As for the case of diffusion limit, the right hand side of Fig. 4.10 presents the error between the DSA scheme computed for ∆t = ε α and the limit ADS scheme computed for ∆t = 10 -7 . Since the DSA scheme is of order 1 uniformly in ε, this error tends to Fig. 4.9. Left: For ∆t = 10 -3 the densities ρ(t = 0.1, x) given by the DSA scheme for different values of ε and the ADS scheme. Right: For ∆t = 10 -3 and α = 0.8, 1, 1.5, the error in ε between the reference ADS scheme and the DSA scheme (log scale).

zero when ∆t and ε go to zero. The difference between the DSA scheme computed for ∆t = 10 -5 and the DSA scheme computed for a range of ∆t as a function of ε (log scale). Right: The difference between the DSA scheme computed for ∆t = ε α and the DSA scheme computed for ∆t = 10 -7 as a function of ∆t (log scale).

Conclusion.

In this paper, we have first presented a formal derivation of diffusion and anomalous diffusion asymptotics from a BGK kinetic equation. This analysis enables us to understand the role of the large velocities induced by the heavytailed equilibrium in the anomalous diffusion asymptotics. Moreover, this formal derivation paves the way to construct three different numerical schemes for the kinetic equation, which all enjoy the asymptotic preserving property in both the diffusion and anomalous diffusion regimes.

The asymptotic preserving character in the case of anomalous diffusion comes from the fact that we managed to take into account the effect of large velocities of the equilibrium. As we saw in the simplest case of a direct fully implicit scheme, the asymptotic is not preserved when these large velocities are truncated. The schemes using the micro-macro formulation as well as the Duhamel formulation of the equations require the use of the same trick on the velocities. Moreover, this last scheme seems to enjoy a uniform accuracy property which will be proved in [START_REF] Crouseilles | Uniformly accurate numerical scheme for kinetic equation in the diffusion and anomalous diffusion limit[END_REF].

In the near future, we aim at extending this work to more general context, considering higher dimensions, non periodic boundary conditions. The case of singular collision frequency also may generate anomalous diffusion (see [START_REF] Ben Abdallah | Anomalous diffusion limit for kinetic equations with degenerate collision frequency[END_REF]) and this also deserves a numerical study that we plan to do in a forthcoming work [START_REF] Crouseilles | Numerical schemes for kinetic equations in the anomalous diffusion limit. Part II: the case of degenerate collision frequency[END_REF].

From

  the explicit expressions of b 0 (k) and c 0 (k) given in (3.23)-(3.24) we get b 0

Fig. 4 . 1 .

 41 Fig. 4.1.For ∆t = 10 -3 , the densities ρ(t = 0.1, x). Left: ISD scheme (with ε = 10 -6 ), by the ADS scheme and the initial data. Right: ISA scheme for different values of ε and the ADS scheme.

  Fig.[START_REF] Bobylev | On some properties of kinetic and hydrodynamic equations for inelastic interactions[END_REF].2. Left: For ∆t = 10 -3 and α = 0.8, 1, 1.5, the error as a function of ε between the reference ADS scheme and the ISA scheme (log scale). Right: The error between the ISA scheme for a range of ∆t = ε α and the ADS scheme as a function of ∆t, which is fixed equal to ε α (log scale).

Fig. 4 .

 4 Fig. 4.3 shows the error curves obtained for different values of ε. It appears that the MMSA scheme is of order 1 in time.

Fig. 4 . 5 .

 45 Fig.4.5. For ε = 1, 0.1, 10 -5 , 10 -6 , the error in time between the DSD scheme and the reference ISD scheme computed for ∆t = 10 -6 (log scale).

Fig. 4 . 6 .

 46 Fig. 4.6. For ∆t = 10 -3 the densities ρ(t = 0.1, x) given by the DSD scheme for different values of ε and the DS scheme.

Fig. 4 . 8 .

 48 Fig.4.8. For ε = 1, 0.1, 10 -4 , 10 -5 , the error as a function of ∆t between the DSA scheme and the DSA scheme computed for ∆t = 10 -5 (log scale).

Fig. 4 . 10 .

 410 Fig.4.10. Left: The difference between the DSA scheme computed for ∆t = 10 -5 and the DSA scheme computed for a range of ∆t as a function of ε (log scale). Right: The difference between the DSA scheme computed for ∆t = ε α and the DSA scheme computed for ∆t = 10 -7 as a function of ∆t (log scale).

For ε = 1, 0.5, 10 -7 , 10 -8 , the error in time between the MMSA scheme computed for a range of ∆t and the reference scheme computed with ∆t = 10 -6 (log scale).

In order to show that the scheme preserves the asymptotic of anomalous diffusion we compute the densities ρ(t = 0.1, x) obtained by the MMSA scheme for a range of ε; we compare them to the density obtained by the ADS scheme computed with the same ∆t. In order to respect the CFL condition, we took ∆t = 10 -4 for ε = 1, 0.5, 0.1, 0.01 and ∆t = 10 -3 for the smallest ε.The left hand side of Fig. 4.4 shows these densities in the case α = 1.5, we observe that they are very close to the anomalous diffusion limit for small ε. On the right hand side of Fig. 4.4, the errors associated to this latter study are displayed in the cases α = 0.8, α = 1 and α = 1.5, the ADS scheme being considered as reference. We observe that the convergence happens with speed α, as expected. 4.3. The integral formulation based scheme. In this section, we focus on the Duhamel formulation based schemes DSD and DSA in both diffusion and anomalous diffusion regimes. We put their properties of consistency, AP character and uniformity with respect to ε, written in Prop. 3.5 and Prop. 3.6, in evidence in these two cases.