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 ANOMALOUS DIFFUSION LIMIT. PART II: DEGENERATE COLLISION FREQUENCY

Keywords: BGK equation, Anomalous diffusion limit, Fractional diffusion equation, Asymptotic preserving schemes, Degenerate collision frequency AMS subject classifications. 35B25, 41A60, 65L04, 65M22

In this work, which is the continuation of [9], we propose numerical schemes for linear kinetic equation which are able to deal with the fractional diffusion limit. When the collision frequency degenerates for small velocities it is known that for an appropriate time scale, the small mean free path limit leads to an anomalous diffusion equation. From a numerical point of view, this degeneracy gives rise to an additional stiffness that must be treated in a suitable way to avoid a prohibitive computational cost. Our aim is therefore to construct a class of numerical schemes which are able to undertake these stiffness. This means that the numerical schemes are able to capture the effect of small velocities in the small mean free path limit with a fixed set of numerical parameters. Various numerical tests are performed to illustrate the efficiency of our methods in this context.

1. Introduction. The theoretical and numerical study of particle system equations has major applications in plasma physics, galactic dynamics, chemotaxis or biological models. Depending on the physical properties of the system, different scales can be considered to describe it. When the mean free path of the particles is large compared to a typical length, a kinetic description of the system is required. Kinetic equations describe the time evolution of a distribution of particles depending on the time, the space and the velocity. This distribution of particles represents the probability for a particle to be at a given point of the phase space (x, v) at the time t. When the mean free path of the particles is small, a macroscopic description can be used. In this regime, the former distribution function is close to an equilibrium and can be averaged in the velocity variable. The relevant quantity is then the density of the particles, which only depends on the time and the space.

The passage from the microscopic to the macroscopic scale is mathematically performed by asymptotic analysis. From a numerical point of view, the small mean free path limit corresponds to a singular perturbation problem which is known to be difficult to simulate (although a lot of works has already been achieved in this context). Indeed, this limit induces the presence of multiple scales in the system, so that the numerical parameters have to solve the smallest ones for stability reasons. As a consequence, the construction of efficient schemes in this context is of great interest.

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 ∈ V ⊂ R d with d = 1, 2, 3. If f stands for a distribution of particles, they undergo the effect of collisions which are modelized here by a linear operator L acting on f through

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

ρ ν (t, x) = ν(v)f (t, x, v) V ν(v)M (v) V =: V ν(v)f (t, x, v)dv V ν(v)M (v)dv , (1.1)
is not the usual density ρ(t, x) = f (t, x, v) V but has been modified to ensure the mass conservation for L: L(f ) V = 0. In the sequel, we will always denote by brackets the integration over v. The equilibrium M (v) will be taken constant on V = {v, |v| ≤ δ} (we choose δ = 1 in the sequel)

M (v) = M 0 , |v| ≤ 1, M 0 > 0, (1.2)
such that M (v) V = 1 and we will consider the case of a degenerate collision frequency

ν(v) = ν 0 |v| d+2+β , |v| ≤ 1, (1.3)
for a given β > 0. Hence, the second order moment of M is not finite

D = |v|≤1 v ⊗ v ν(v) M (v)dv = +∞. (1.4)
In order to capture a nontrivial asymptotic model, a suitable scaling has to be considered (see for instance [START_REF] Ben Abdallah | Anomalous diffusion limit for kinetic equations with degenerate collision frequency[END_REF]). We consider the kinetic equation

ε α ∂ t f + εv • ∇ x f = L(f ), with α = 2d + 2 + β d + 1 + β ∈ (1, 2), (1.5)
where ε > 0 is the Knudsen number, which is proportional to the mean free path of the particles. Equation (1.5) has to be supplemented with an initial condition f (0, x, v) = f 0 (x, v) and spatial periodic boundary conditions are considered. Considering α = 2, an equilibrium function M and a collision frequency ν such that the integral (1.4) is finite, (1.5) degenerates into a diffusion equation for ρ ∂ t ρ(t, x) -∇ x • (D∇ x )ρ(t, x) = 0, when ε goes to zero. There is a huge literature on this subject, we refer for example to [START_REF] Larsen | Asymptotic solution of neutron transport problems for small mean free paths[END_REF][START_REF] Wigner | Nuclear reactor theory[END_REF][START_REF] Bensoussan | Boundary layers and homogenization of transport processes[END_REF][START_REF] Degond | Diffusion limit for non homogeneous and non-microreversible processes[END_REF] for details. When the diffusion coefficient (1.4) is no longer finite, this asymptotic breaks down which can be interpreted by saying that the time scale t ∼ ε -2 is too long. Then, it is necessary to consider an appropriate scaling to get a nontrivial limit for (1.5). For instance, it is known (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]) that when the collision frequency is not degenerated for small velocities and when the equilibrium M is an heavy-tailed function (equivalent to a power of v for large velocities) (1.5) degenerates into a fractional diffusion equation when ε goes to zero. As the scaling ε α is not an integer power of ε this fractional diffusion limit is called an anomalous diffusion limit. This scaling arises 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]), plasmas (see [START_REF] Del Castillo-Negrete | Non-diffusive transport in plasma turbulence: A fractional diffusion approach[END_REF]) or even economy ( [START_REF] Kleiber | Statistical size distributions in economics and actuarial sciences[END_REF]). The case with degenerate collision frequency we are considering has been studied in [START_REF] Ben Abdallah | Anomalous diffusion limit for kinetic equations with degenerate collision frequency[END_REF]. It does not correspond to a particles dynamics physical situation but rather to the modelization of some chains of oscillators, see [START_REF] Jara | Limit theorems for additive functionals of a Markov chain[END_REF][START_REF] Mellet | Anomalous energy transport in FPU-β chain[END_REF].

When ε goes to 0, the solution of (1.5) converges to ρ(t, x)M (v) where ρ is the solution of the fractional diffusion equation that can be written in Fourier variable ∂ t ρ(t, k) = -κ|k| α ρ(t, k), (1.6) where ρ stands for the space Fourier transform of ρ, k is the Fourier variable and κ is a constant which depends on M, ν and α and is defined by

κ = M 0 ν 1-α 0 d + 1 + β 1 |v| d+α (v • e) 2 1 + (v • e) 2 R d , (1.7)
for any e ∈ R d such that |e| = 1 (note that κ does not depend on e). The fractional diffusion equation can also be written in the space variable

∂ t ρ(t, x) = - 1 c d,α M 0 ν 1-α 0 Γ(α + 1) d + 1 + β (-∆ x ) α 2 ρ(t, x), (1.8)
where Γ is the usual Gamma 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 )
. The fractional Laplacian (-∆ x ) α 2 can be defined with its Fourier variable

(-∆ 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. is the principal value distribution. From a numerical point of view, the presence of ε makes the numerical simulations challenging. Indeed, if one uses a naive approach, a relation linking the space and time numerical parameters to the stiffness ε (typically ∆t ∼ ε α , where ∆t is the time step) is required for stability reasons, which implies a redhibitory cost when ε is small. The Asymptotic Preserving (AP) schemes are designed to overcome this strong restriction on the numerical parameters. Indeed, if we consider the continuous problem P ε degenerating into P 0 when ε goes to zero, we ask the scheme S h ε to be consistent with the problem P ε when the discretization parameter h tends to 0 and to degenerate into a scheme S h 0 solving the problem P 0 when ε → 0. An AP scheme can also enjoy the stronger property of being Uniformly Accurate (UA). It means that the accuracy in h of the scheme does not depend of ε. Note that it is very desirable for an AP scheme to enjoyy the stronger property of being Uniformly Accurate but this is not automatic and is a difficult issue in general. For kinetic equations, AP scheme has already been developed (see [START_REF] Jin | Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations[END_REF]) and there is an important literature for the diffusion asymptotic, for example [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][START_REF] Mieussens | On the asymptotic preserving property of the unified gas kinetic scheme for the diffusion limit of linear kinetic models[END_REF][START_REF] Lemou | Micro-macro schemes for kinetic equations including boundary layers[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] Klar | An asymptotic-induced scheme for nonstationary transport equations in the diffusive limit[END_REF][START_REF] Carrillo | Numerical schemes of diffusion asymptotics and moment closures for kinetic equations[END_REF][START_REF] Bardos | Diffusion approximation and computation of the critical size[END_REF]. In the case of a fractional diffusion scaling, the previous schemes do not enjoy the AP property. Indeed, as the velocities must be discretized with a finite number of points, the diffusion coefficient (1.4), which is an infinite integral in the case of the anomalous diffusion, appears to be always finite at the numerical level. Hence, the effects of the low velocities are missed and the approaches for the classical diffusion show only a numerical diffusion which does not match with the expected fractional diffusion.

In the case of the fractional diffusion induced by a heavy-tailed equilibrium, we set in a previous work ([9]) three numerical schemes enjoying the AP property. As the fractional diffusion limit did not naturally appear in the formulation of the schemes, we had to suitably modify them to make the fractional diffusion arise numerically. Let us remark that when α ∈ (1, 2), an alternative approach has been proposed in [START_REF] Wang | An asymptotic-preserving scheme for linear kinetic equation with fractional diffusion limit[END_REF] for the case of heavy tail equilibria. In this paper, we present an extension of the strategy presented in [START_REF] Crouseilles | Numerical schemes for kinetic equations in the anomalous diffusion limit. Part I: the case of heavy-tailed equilibrium[END_REF] when the fractional diffusion comes from the degeneracy of the collision frequency. In addition to the technical difficulties which are inherent to the velocity dependency of ν there are many important differences with the heavy-tail case that we have to face. First, to numerically capture the important effects of the low velocities (instead of large velocities in the heavy-tail case). Second, in the present case, the renormalized density ρ ν defined in (1.1) is different from the usual density ρ while these two quantities coincide in the heavy-tail case. This property induces an additional difficulty, as we will see, and a specific treatment is required in order to cover an asymptotic preserving property for the scheme. Lastly, the convergence rate to the fractional diffusion limit is very slow, so that the intermediate regime involves a large range of ε; as a consequence, the numerical validation turns out to be more difficult.

We will present three different Asymptotic Preserving schemes. First, we consider a fully implicit scheme for (1.5) in which Fourier space variable is considered. The delicate issue of the velocity discretization is highlighted to explain how the fractional diffusion arises. Then, a second scheme is proposed which allows to avoid, when desired, both the use of the space Fourier transform and that of a fully implicit scheme. A micro-macro formulation of (1.5) is then studied to write an semi-implicit scheme enjoying the AP property. Once again, the terms giving theoretically the anomalous diffusion limit will have to be treated with care to ensure the AP character of the scheme. Finally, a scheme based on a Duhamel formulation of (1.5) is proposed, which enjoys the stronger property of being uniformly accurate (UA) with respect to ε. These results were announced in [START_REF] Crouseilles | Multiscale numerical schemes for kinetic equations in the anomalous diffusion limit[END_REF]. Moreover, a proof of the convergence of the solution of (1.5) to the solution of (1.8) is given in [START_REF] Ben Abdallah | Anomalous diffusion limit for kinetic equations with degenerate collision frequency[END_REF]. We propose a formal derivation of the convergence rate in ε for this limit. This convergence rate is illustrated in the numerical tests that we present in the last part.

The paper is organized as follows. In the next section, we will start by establishing a formal derivation of the fractional diffusion limit equation (1.8), and the rate of this convergence is formally derived (the related computation is postponed in the appendix). Then, in Section 3, we present the three asymptotic preserving schemes for (1.5). Eventually, we will present in the last section some numerical tests to highlight the properties of our schemes.

2. Formal derivation of the fractional diffusion limit. In this section, we formally derive the fractional diffusion (1.6) from the kinetic equation (1.5), where the computations are based on spatial Fourier transform. This will be the basis of the numerical methods proposed in Section 3. The proof has been done in [START_REF] Ben Abdallah | Anomalous diffusion limit for kinetic equations with degenerate collision frequency[END_REF].

We remind that f (t, k, v) (resp. ρ(t, k)) denotes the spatial Fourier transform of f (t, x, v) (resp. ρ(t, x))

f (t, k, v) = R d e -ik•x f (t, x, v)dx.
Proposition 2.1. Let f be the solution of (1.5) in the Fourier variable for a given initial condition f0 . Then, when ε → 0, f converges to ρM with ρ solution of the fractional diffusion equation

∂ t ρ(t, k) + κ|k| α ρ(t, k) = 0,
with κ defined by (1.7) and with initial condition ρ(0, k) = f0 (k, v)

V .
Proof. Considering the Fourier transform (in space) of (1.5), we get

∂ t f + ε -α (ν + iεk • v) f = ε -α ν ρν M.
(2.1)

Integrating with respect to the velocity variable leads to

∂ t ρ + iε 1-α k • v f V = 0. (2.2)
From (2.1), we can write

f = ν ρν M ν + iεk • v -ε α ∂ t f ν + iεk • v , (2.3)
so that f -→ ε→0 ρM and (2.2) becomes

∂ t ρ + iε 1-α νk • vM ν + iεk • v V ρν = O(ε γ ), γ > 0. (2.4)
The fact that the right hand side is negligible compared to the left hand side will be detailed in the Appendix. Then, the fractional diffusion operator is obtained by considering the limit ε → 0 in the second term of the left hand side. Let us focus on this term in the sequel. First, since M and ν are even, we write

iε 1-α νk • vM ν + iεk • v V = ε 2-α ν(k • v) 2 M ν 2 + ε 2 (k • v) 2 V .
We can remark that the limit ε → 0 is not defined since the term into brackets goes to infinity when ε → 0 because of the degeneracy of ν for small velocities. Then, to capture the effect of small velocities, we perform the change of variable w = ε|k|v/ν to get

lim ε→0 iε 1-α νk • vM ν + iεk • v V = M 0 ν 1-α 0 d + 1 + β 1 |w| d+α (w • e) 2 1 + (w • e) 2 R d |k| α = κ|k| α
where κ has been defined in (1.7). Moreover, since when ε → 0, we have ρ → ρν , we then deduce the fractional diffusion asymptotic model from (2.4).

Furthermore, we have the following result for the convergence rate (at the formal level). The proof is postponed to the Appendix.

Proposition 2.2. Let f be the solution of (1.5) in the Fourier variable for a given initial condition f0 . We get the following rate of convergence

f -ρM = O ε d d+1+β + O ε β d+1+β .
3. Numerical schemes. In this section, we present three different numerical schemes, based respectively on a fully implicit scheme, a micro-macro formulation and an integral formulation of (1.5), designed to approximate the solution of (1.5). We will denote t n = n∆t, 0 ≤ n ≤ N the time discretization, such that N ∆t = T and we will set f n ≃ f (t n ). The space domain Ω is bounded and we consider periodic conditions, allowing the use of the Fourier variable. At the discrete level, we will use the discrete Fourier transform. The computation of the numerical solution requires the use of a velocity discretization. In the tests, we will consider a first order quadrature method based on a uniform fixed set of N v points v j , symmetrically distributed in the domain |v| ≤ 1. In the sequel, we will denote by • Nv,D the discrete integration in v on a domain D :

f Nv,D = ∆v j∈[[1,Nv ]] v j ∈D f (v j ). (3.1)
We want these schemes to be Asymptotic Preserving (AP), indeed that for any initial condition f 0 (i) when ε is fixed, the solution given by the numerical scheme is consistent with the solution of (1.5), (ii) when the numerical parameters are fixed, the scheme degenerates into a scheme solving the anomalous diffusion equation (1.8) when ε tends to 0.

Since the fractional Laplacian of the asymptotic equation (1.6) comes from the small velocities of the degenerate collision frequency, it is necessary to take it into account in the schemes to ensure the AP property. One idea would be to consider an adaptive grid which is refined around v = 0 but this will impose to link the velocity discretization to ε, which breaks down the AP property. We propose a method based on an analytical modification of the terms degenerating into the fractional Laplacian in the scheme. This modification consists into applying a change of variable in some integrals in v of the scheme before discretizing them. Once it is done, we are able to show that the three schemes we propose preserve the asymptotic of fractional diffusion for a fixed set of numerical parameters.

3.1. Implicit scheme. The first idea to write a scheme for (1.5) which has the AP property is to use a fully implicit formulation of (1.5) in Fourier variable. It is well-known that in the case of the classical diffusion limit, such a scheme preserves the asymptotics, but here in the case of the fractional diffusion limit, it is not true, and a suitable modification of the scheme has to be performed to ensure the AP property.

We start with (1.5) written in Fourier variable and we consider a fully implicit time discretization

f n+1 -f n ∆t + ν ε α 1 + i εk • v ν f n+1 = ν ε α ρn+1 ν M, that is f n+1 = 1 -λ(v) 1 + iλ(v) εk•v ν f n + λ(v) 1 + iλ(v) εk•v ν ρn+1 ν M, (3.2)
with λ(v) = ∆tν(v)/(ε α + ∆tν(v)). Remark that 0 ≤ λ(v) ≤ 1, λ goes to 1 when ε → 0 and λ goes to 0 as ∆t → 0. So when ε → 0 with fixed ∆t, we get

f n+1 = ρn+1 ν M, (3.3)
which is the expected limit for f . Now, it remains to find an expression for ρn+1 ν . To do so, we multiply (3.2) by ν and integrate in v to get

ρn+1 ν = 1 ε α +∆tν ν f n 1 + iλ(v) εk•v ν V 1 ε α +∆tν νM 1 + iλ(v) εk•v ν V + I , with I = 1 ε α νM λ(v) 2 εk•v ν 2 1 + λ(v) 2 εk•v ν 2 V . (3.4)
At the continuous level, the term I gives the anomalous diffusion limit when ε goes to zero. Indeed, the change of variable w = ε|k|v/ν(v) in I enables to capture the small velocity effects so that lim ε→0 I = κ|k| α , where κ has been defined in (1.7). However, the situation is completely different at the numerical level, namely when a discretization in velocity with N v points has to be used. Indeed, a naive quadrature of I given by (3.4) would lead to

I ≈ ε 2-α Nv j=1 ν(v j )M (v j )∆t 2 (k • v j ) 2 (ε α + ∆tν(v j )) 2 + ε 2 ∆t 2 (k • v j ) 2 ∼ ε→0 ε 2-α Nv j=1 (k • v j ) 2 ν(v j ) M (v j ).
This implies that I goes to zero when ε goes to zero since the velocity discretization is fixed (remember that α < 2). Therefore, the asymptotic numerical scheme becomes ρn+1 ν = ρn ν when ε → 0, which is not a consistent numerical scheme for the asymptotic model. To recover the correct asymptotic behavior, the main idea is to transform the expression of I at the continuous level before applying the discretization in order to make the anomalous diffusion operator appear in it. To do so, we perform a change of variables w = ε|k|v/ν(v) in I to get

I = |k| α ν 1-α 0 M 0 d + 1 + β |w|≥ ε|k| ν 0 1 |w| d+α ∆tϕ ε α +∆tϕ 2 (w • e) 2 1 + ∆tϕ ε α +∆tϕ 2 (w • e) 2 dw,
where e denotes any unitary vector and ϕ = ν 0 (ε|k|/(ν 0 |w|)) d+2+β d+1+β . Note that I does not depend on the choice of e. Once the change of variables has been done, the velocity discretization can be performed so that we get

I ≈ |k| α ν 1-α 0 M 0 d + 1 + β 1 |w| d+α ∆tϕ ε α +∆tϕ 2 (w • e) 2 1 + ∆tϕ ε α +∆tϕ 2 (w • e) 2 Nv,|wmax|≥|w|≥ ε|k| ν 0 (3.5)
where • Nv ,D has been defined in (3.1). Let us remark that when ε goes to zero, this approximation of I degenerates towards

|k| α ν 1-α 0 M 0 d + 1 + β 1 |w| d+α (w • e) 2 1 + (w • e) 2 Nv,|w|≤|wmax| =: κ, (3.6)
which is a consistent approximation of κ given by (1.7). Let us remark that the numerical integration is done up to a maximum discrete velocity |w max |. Eventually, we have the following proposition Proposition 3.1. We consider the following scheme defined for all k and for all time indices 0 ≤ n ≤ N, N ∆t = T by

ρn+1 ν = 1 ε α +∆tν ν f n 1 + iλ(v) εk•v ν Nv ,V 1
ε α +∆tν νM 1 + iλ(v) εk•v ν Nv,V + I , f n+1 = (1 -λ(v)) f n + λ(v)ρ n+1 ν M 1 + iλ(v) εk • v ν , (3.7)
with I given by

(3.5), ρn+1 = f n+1 Nv,V , λ(v) = ∆tν(v)/(ε α + ∆tν(v)
) and 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. (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), with κ given by (3.6).

(3.8)

Remark 1. We can go further, modifying the approximation (3.5) of I by

I ≈ ∆t 2 (ε α-d+2+β d+1+β + ∆t) 2 κ,
where κ is defined in (3.6). The so-obtained scheme enjoys the same properties stated in the above proposition.

Proof. The first point of (i) is a direct consequence of the use of a fully implicit Euler scheme for (1.5). For the conservation of the total mass, let us write (3.7) for k = 0 and integrate with respect to 3.2. Micro-macro decomposition based scheme. In the previous part, we presented a fully implicit scheme designed to solve (1.5); however the use of Fourier transform may be restrictive in some cases and one may have to use the original space variable. Therefore, making the transport part implicit in time may induce high computational cost especially in the multidimensional case. Here we write a scheme based on a micro-macro decomposition of the distribution function f in which the transport part is explicit. As for the previous scheme, suitable modifications have to be done to ensure the AP character.

v νM ε α + ∆tν Nv,V ρn+1 ν (0) = ν ε α + ∆tν f n (0) Nv,V , (3.9)
ρn+1 (0) = ρn (0) -∆t ν ε α + ∆tν f n (0) Nv ,V + ∆tν ε α + ∆tν M Nv ,V ρn+1
Denoting

f (t, x, v) = ρ(t, x)M (v) + g(t, x, v
) such that g V = 0, we insert this decomposition into (1.5) before integrating in v to get

∂ t ρ + ε 1-α v • ∇ x g V = 0. (3.11)
Then we multiply this last equality by M and subtract it from (1.5) to write an expression for g

∂ t g + ε 1-α v • ∇ x ρM + ε 1-α (v • ∇ x g -v • ∇ x g V M ) = - 1 ε α ν g - νg V M νM V .
(3.12)

From (3.11) and (3.12) we can write the following semi-discrete micro-macro scheme (see [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] Lemou | Micro-macro schemes for kinetic equations including boundary layers[END_REF])

                   ρ n+1 -ρ n ∆t + ε 1-α v • ∇ x g n+1 V = 0 g n+1 -g n ∆t + ε 1-α v • ∇ x ρ n M + ε 1-α (v • ∇ x g n -v • ∇ x g n V M ) = - 1 ε α ν g n+1 - νg n+1 V νM V M . (3.13)
As previously, if the velocity integrations are done directly using the approximation

• V ≈ • Nv,V
where the discrete bracket • Nv,V is given by (3.1), the so-obtained scheme does not enjoy the AP property. Indeed, as g n+1 goes to zero when ε goes to zero, we have

εv • ∇ x ρ n M = -ν g n+1 - νg n+1 Nv,V νM Nv,V M + o(ε), that is, for nonzero v g n+1 = νg n+1 Nv,V νM Nv,V M - ε ν v • ∇ x ρ n M + o(ε).
Then, we insert this expression in the first line of (3.13), reminding that M is even, to get

ρ n+1 -ρ n ∆t -ε 2-α M ν v • ∇ x (v • ∇ x ρ n ) Nv,V = o(ε).
As in the previous part, the term into brackets is finite although it is not finite at the continuous level. Hence, when ε → 0, (3.13) degenerates into ρ n+1 = ρ n which of course does not solve the fractional diffusion equation. To obtain the right asymptotic scheme, we must modify (3.13) to make the anomalous diffusion limit appear. To do so, we first express the term v

• ∇ x g n+1 V = v • ∇ x f n+1
V appearing in the first line of (3.13) with a semi-discrete implicit formulation of (1.5)

f n+1 = λ(v) I + ε λ(v) ν v • ∇ x -1 ρ n+1 ν M + (1 -λ(v)) I + ε λ(v) ν v • ∇ x -1 f n , with λ(v) = ∆tν(v)/(ε α + ∆tν(v)). As ελ(v)/ν(v) = O(∆t)
we will rather inject the following approximation for f n+1 in (3.13) (recalling that v ∈ V )

f n+1 = λ(v) I + ε λ(v) ν v • ∇ x -1 ρ n+1 ν M + (1 -λ(v)) I -ε λ(v) ν v • ∇ x f n .
Indeed, since the scheme is first order accurate, it is possible to remove terms of order ∆t 2 with no incidence on its global accuracy. Hence, we decide to keep only the ones which play a role in the AP character of the scheme. With additional simplifications, the first line of (3.13) becomes

ρ n+1 -ρ n ∆t + 1 ε α ελ(v)v • ∇ x I + ελ(v) ν v • ∇ x -1 ρ n+1 ν M V + ε 1-α (1 -λ(v)) v • ∇ x g n V = 0, (3.14) because (1 -λ(v)) v • ∇ x I - ελ(v) ν v • ∇ x f n V = (1 -λ(v)) v • ∇ x g n V + O(∆t).
As the anomalous diffusion limit comes from the second term of (3.14), we have to rewrite it in order to make the fractional Laplacian appear. First, we have to express ρ ν in terms of ρ; to do that, we inject the decomposition f = ρM + g in the definition of ρ ν

ρ ν = νf V νM V = ν(ρM + g) V νM V = ρ + νg V νM V .
Then, using space Fourier transform, we rewrite (3.14) as

ρ n+1 -ρ n ∆t + F -1 I ρn+1 + F -1 I νĝ n+1 V νM V + ε 1-α (1 -λ(v)) v • ∇ x g n V = 0, (3.15
) where F -1 denotes the inverse Fourier transform and

I = ε 1-α iλ(v)k • vM (v) 1 + i ελ(v) ν k • v V = ε -α λ(v) 2 (εk•v) 2 ν M 1 + λ(v) 2 εk•v ν 2 V
.

Let us remark that for k = 0, we have I = 0. For k = 0, the change of variables w = ε|k|v/ν(v) is performed in I before applying the discretization • Nv ,V (using the definitions (1.2) and (1.3) of M and ν). Then, we have

I ≈ |k| α ν 1-α 0 M 0 d + 1 + β 1 |w| d+α (w • e) 2 (w • e) 2 + (1 + ψ) 2
Nv ,|wmax|≥|w|≥ ε|k| ν 0 (3.16) where ψ = ε q ∆t ν r 0 (|w|/|k|) p , p = d+2+β d+1+β and q = d d+1+β , r = 1 d+1+β . To compute (ρ n+1 , g n+1 ) recursively with the second line of (3.13) and (3.15), we need to determine νg n+1 Nv,V . To do that, we start with the with the second line of (3.13) which gives an expression for g n+1

g n+1 = g n -∆t ε α T n 1 + ∆t ε α ν + ∆tνM νM Nv (ε α + ∆tν) νg n+1 Nv,V ,
with

T n = εv • ∇ x ρ n M + ε v • ∇ x g n -v • ∇ x g n Nv,V M . (3.17)
Multiplying the expression for g n+1 by ν and integrating it using • Nv,V , we find an expression for νg n+1

Nv,V

νg n+1 Nv,V = νg n 1+ ∆t ε α ν(v) Nv,V -∆t ε α νT n 1+ ∆t ε α ν Nv,V 1 - ∆t ε α νM Nv ,V ν 2 M 1+ ∆t ε α ν Nv,V
.

Finally, we have the following proposition.

Proposition 3.2. We consider the following scheme (discretized in time and velocity) defined for all x ∈ R d , v ∈ V and all time indices 0 ≤ n ≤ N with N ∆t = T (T > 0), by

                               νg n+1 Nv,V = νM Nv,V νg n + T n ε α + ∆tν Nv ,V ν 2 M ε α + ∆tν Nv ,V g n+1 -g n ∆t + T n = - 1 ε α ν g n+1 - νg n+1 Nv ,V νM Nv ,V M ρ n+1 -ρ n ∆t + F -1 I ρn+1/2 + F -1 I νĝ n+1 Nv ,V νM Nv ,V +ε 1-α (1 -λ(v)) v • ∇ x g n Nv,V = 0, (3.18)
with T n defined in (3.17), λ(v) = ∆tν(v)/(ε α + ∆tν(v)), where F -1 denotes the inverse Fourier transform in space and with the initial conditions

ρ 0 (x) = ρ(0, x) = f 0 (x, v) Nv,V and g 0 (x, v) = f 0 (x, v) -ρ(0, x)M (v).
The quantity ρn+1/2 can be chosen equal to ρn or to ρn+1 depending on the desired asymptotic scheme (explicit or implicit in time) and I is given by (3.16). This scheme has the following properties: (i) The scheme is of order 1 for any fixed ε and preserves the total mass. (ii) The scheme is AP: for a fixed ∆t, it solves the anomalous diffusion equation when ε goes to zero

ρn+1 (k) -ρn (k) ∆t = -κ|k| α ρn+1/2 (k),
where κ is defined by (3.6). Remark 2. Assuming that the space and velocity discretizations have been fixed, (ρ n , g n ) are computed recursively with (3.18). The first line of the scheme enables to compute νg n+1 Nv,V which is reported in the second and third lines to compute g n+1 and ρ n+1 .

As the term I defined in (3.16) has been precomputed to make the anomalous diffusion operator clearly appear, the velocity discretization is not a difficult issue anymore. Once it is defined, the integrals in v are computed with any quadrature method, such as rectangle formula. Finally, a standard upwind scheme is used to approximate the spatial derivative appearing in the term T n defined in (3.17).

Remark 3. To guarantee the structure of the micro-macro model, it is important that the scheme satisfies g n+1 Nv,V = 0. To do that, g n+1 is modified consistently as g n+1 ← g n+1 -g n+1 Nv,V M . Remark 4. As in the fully implicit scheme, we can go further by approximating consistently I given by (3.16) as

I = ∆t 2 ∆t 2 + ε q (∆t + ε q ) |k| α κ,
where κ is given by (3.6) and q = d/(d + 1 + β). The so-obtained micro-macro based scheme still enjoys the same properties stated in the above proposition.

Proof. The proof of (i) is immediate from the derivation of the scheme, let us prove (ii). At first, we remark that T n goes to 0 as ε → 0, then the first line of (3.18) gives that ν(v)g n+1 Nv ,V goes to 0 as ε → 0 because g n Nv = 0. After that, the second line of (3.18) implies that g n+1 = o(1). When reported in the fourth line of (3.18), as ε

1-α (1 -λ(v)) = O(ε), we get ρ n+1 -ρ n ∆t + F -1 I ρn+1/2 = o(1).
Finally, I is equivalent to κ|k| α for small ε, with κ given by (3.6).

3.3. Scheme based on an integral formulation of the equation (1.5). In the previous parts, we wrote two AP schemes solving (1.5). In this section, in the spirit of [START_REF] Crouseilles | Numerical schemes for kinetic equations in the anomalous diffusion limit. Part I: the case of heavy-tailed equilibrium[END_REF], we present a scheme based on an integral formulation of (1.5). However, the presence of the collision frequency ν requires specific developments to ensure the AP property. In particular, for general initial data, the strategy proposed in [START_REF] Crouseilles | Numerical schemes for kinetic equations in the anomalous diffusion limit. Part I: the case of heavy-tailed equilibrium[END_REF] leads to an AP scheme for ρ ν instead of ρ (which are different if f is not at equilibrium). Then, the present section is decomposed into two parts: first, the case of equilibrium initial data is tackled, for which the scheme is proved to be uniformly accurate, and second the more delicate case of general initial data is studied, for which the scheme is proved to be AP.

Equilibrium initial data.

We start from the Duhamel form for f

f (t, k, v) = e -t ε α (ν+iεk•v) f0 (k, v) + νM t ε α 0 e -s(ν+iεk•v) ρν (t -ε α s, k)ds, (3.19)
which we multiply by ν, integrate in v and evaluate at t = t n+1 to get

K ρν (t n+1 , k) = Â0 (t n+1 , k) + t n+1 ε α 0 ν 2 M e -s(ν+iεk•v) V ρν (t n+1 -ε α s, k)ds, (3.20)
where

K = νM V , and Â0 (t, k) = νe -t ε α (ν+iεk•v) f0 (k, v) V . (3.21)
Then, we get

K ρν (t n+1 , k) = n j=0 t j+1 ε α t j ε α ν 2 M e -s(ν+iεk•v) V ρν (t n+1 -ε α s, k)ds + Â0 (t n+1 , k), (3.22)
where Â0 and K are defined in (3.21). We perform a quadrature of order 2 in the integrals in s. Assuming that the time derivatives of ρν are uniformly bounded in ε,

we have ∀j ∈ [[0, N -1]], ∀s ∈ tj ε α , tj+1 ε α , ρν (t n+1 -ε α s, k) = 1 - ε α s -t j ∆t ρν (t n+1-j , k) + ε α s -t j ∆t ρν (t n-j , k) + O ∆t 2 .
(3.23)

Inserting this expression in the integral terms of (3.22) leads to

t j+1 ε α t j ε α ν 2 M e -s(ν+iεk•v) V ρν (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.24)
where the remainder is bounded independently of ε and where we denoted

b j (ε, k) = t j+1 ε α t j ε α ε α s -t j ∆t ν 2 M e -s(ν+iεk•v) V ds, (3.25)
c j (ε, k) = t j+1 ε α t j ε α 1 - ε α s -t j ∆t ν 2 M e -s(ν+iεk•v) V ds. (3.26)
Then, we use the quadrature (3.24) in (3.22) to write (using the obvious notations b j := b j (ε, k)) and c j := c j (ε, k))

ρν (t n+1 , k) = 1 K n j=0 c j ρν (t n+1-j , k) + b j ρν (t n-j , k) + O(∆t 2 ) + 1 K Â0 (t n+1 , k).
As N ∆t = T , then n∆t 2 ≤ T ∆t and the last expression permits to write a first order scheme for ρν 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) K -c 0 . (3.27)
To ensure the AP property for this scheme, the discretization of the terms b j and c j is crucial: as in the previous cases, their limit as ε goes to zero makes the anomalous diffusion operator appear, but when computed numerically it does not arise because of the effects of the low velocities missed with the discretization. Indeed, we consider the velocity integration appearing in (3.25)-(3.26)

ν 2 M e -s(ν+iεk•v) V = ν 2 M e -s(ν+iεk•v) -e -ν(v)s V + ν 2 M e -ν(v)s V . (3.28)
At the continuous level, the first term in the right hand side gives the anomalous diffusion when ε goes to zero. Indeed, with M and ν defined by (1.2) and (1.3), the change of variables w = ε|k|v/ν in the integral in v reads

ν 2 M e -s(ν+iεk•v) -e -ν(v)s V = ε α |k| α ν 1-α 0 M 0 d + 1 + β ϕ 1 |w| d+α e -ϕs(1+iw•e) -e -sϕ
j (ε, k) = ε α ∆t M e -t j ε α ν 1 -e -∆t ε α ν Nv,V -νM e -t j+1 ε α ν Nv,V + ε α |k| α ν 1-α 0 M 0 d + 1 + β 1 |w| d+α e -t j+1 ε α ϕ - e -t j+1 ε α ϕ(1+iw•e) 1 + iw • e (3.30) + ε α ∆t   e -t j ε α ϕ(1+iw•e) 1 -e -∆t ε α ϕ(1+iw•e) ϕ(1 + iw • e) 2 - e -t j ε α ϕ 1 -e -∆t ε α ϕ ϕ     Nv,W c j (ε, k) = νM e -t j ε α ν Nv ,V - ε α ∆t M e -t j ε α ν 1 -e -∆t ε α ν Nv,V + ε α |k| α ν 1-α 0 M 0 d + 1 + β 1 |w| d+α e -t j ε α ϕ(1+iw•e) 1 + iw • e -e -t j ε α ϕ (3.31) - ε α ∆t   e -t j ε α ϕ(1+iw•e) 1 -e -∆t ε α ϕ(1+iw•e) ϕ(1 + iw • e) 2 - e -t j ε α ϕ 1 -e -∆t ε α ϕ ϕ     Nv,W , with W = w ∈ R d , |w max | ≥ |w| ≥ ε|k|/ν 0 .
The so-obtained scheme (3.27) for ρν where b j , c j are given by (3.30)-(3.31) enables us to find a way to recover the distribution function f in Fourier variable. Indeed, with a Duhamel formulation of (1.5) integrated between t n and t n+1 it is possible to write an expression for f (t n+1 , k, v)

f (t n+1 , k, v) = e -∆t ε α (ν+iεk•v) f (t n , k, v) + ∆t ε α 0 e -s(ν+iεk•v) ρν (t n+1 -ε α s, k)ds νM.
(3.32)

Then, the quadrature (3.23) leads to

f (t n+1 , k, v) = e -∆t ε α (ν+iεk•v) f (t n , k, v) + νM (γ ρν (t n+1 , k) + β ρν (t n , k)) + O(∆t 2 ),
where

γ = ∆t ε α 0 1 - ε α s ∆t e -s(ν+iεk•v) ds, β = ∆t ε α 0 ε α s ∆t e -s(ν+iεk•v) ds. (3.33)
This provides a first order scheme for f using the values ρn ν determined with (3.27)

f n+1 = e -∆t ε α (ν+iεk•v) f n + νM γ ρn+1 ν + β ρn ν ,
with the values of β and γ computed exactly to ensure the AP property of the scheme

γ = 1 ν + iεk • v - ε α ∆t 1 -e -∆t ε α (ν+iεk•v) (ν + iεk • v) 2 , β = ε α ∆t 1 -e -∆t ε α (ν+iεk•v) (ν + iεk • v) 2 - e -∆t ε α (ν+iεk•v) ν + iεk • v .
(3.34)

A simple integration of f n+1 with • Nv,V gives an expression for ρn+1 . Eventually, we have the following proposition.

Proposition 3.3. We consider the following scheme defined for all k and for all time index 0 ≤ n ≤ N, N ∆t = T by For an initial data f0 (k, v) at equilibrium, this scheme enjoys the following properties: (i) The scheme is first order in time and preserves the total mass. (ii) The scheme is AP: for a fixed ∆t, it solves the anomalous diffusion equation when ε goes to zero

           ρ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) K -c 0 f n+1 = e -∆t ε α (ν+iεk•v) f n + νM γ ρn+1 ν + β ρn ν , (3.35
ρn+1 (k) -ρn (k) ∆t = -κ|k| α ρn+1 (k)
,

where ρn = f n N v,V
and κ is defined by (3.6).

(iii) Moreover, the semi-discrete-in-time scheme enjoys the UA property: it is first order uniformly in ε

∃C > 0, sup ε∈( 0,1 ] f N (•, •) -f (T, •, •) L ∞ k,v ≤ C∆t.
Remark 5. The numerical tests (see Fig. 4.3) suggests that this scheme is of order 2 for a fixed value of ε and of order 1 uniformly in ε.

Proof. For the mass conservation (i), we integrate (at the discrete level) with respect to v the second equation of (3.35) and evaluate it at k = 0. For an initial data at equilibrium, we have ρ0

ν (k = 0) = ρ0 (k = 0) so that ρ1 = e -∆t/ε α ν M Nv,V ρ0 + νM γ Nv ,V ρ1 ν + νM β Nv,V ρ0 ν , for k = 0.
From the first equation of (3.35), we easily get ρ1 ν (k = 0) = ρ0 (k = 0) so that the previous equation leads to ρ1 (k = 0) = ρ0 (k = 0). The mass conservation is deduced by induction, assuming that ρℓ ν (k = 0) = ρℓ (k = 0) for ℓ = 0, . . . , n. For the AP character (ii), let us remark that with (3.30) and (3.31), the formula for ρν in (3.35)

degenerates into ρn+1 ν (k) -ρn ν (k) ∆t = -κ|k| α ρn+1 ν (k).
Moreover, the second line of (3.35) gives f n+1 = ρn+1 ν M (v) for n ≥ 0 when ε → 0, that implies ρn+1 = ρn+1 ν for small ε. As the initial condition is at equilibrium, we also have ρ0 = ρ0 ν . Eventually, when ε → 0, (3.35) degenerates into

ρn+1 (k) -ρn (k) ∆t = -κ|k| α ρn+1 (k),
for all n ≤ N , which proves the AP property. Now, let us prove that the semi-discrete (in time) scheme is first order accurate uniformly in ε (iii). Using the notations (3.30)-(3.31), (3.22) writes

ρν (t n+1 , k) = 1 K Â0 (t n+1 , k) + 1 K n j=0 c j ρν (t n+1-j , k) + b j ρν (t n-j , k) + ∆t 2 K F, (3.36)
where the remainder F is given by

F = - 1 2 n j=0 t j+1 ε α t j ε α 1 - ε α s -t j ∆t ε α s -t j ∆t 2 ∂ 2 t ρν (t n+1 -ε α ξ 1 (s), k) (3.37) + ε α s -t j ∆t t j+1 -ε α s ∆t 2 ∂ 2 t ρν (t n+1 -ε α ξ 2 (s), k) ν 2 M e -s(ν+iεk•v) V ds, with ξ 1 (s), ξ 2 (s) ∈ (0, ∆t/ε α ). As 1 ≥ 1 - ε α s-tj ∆t ≥ 0 and 1 ≥ ε α s-tj ∆t ≥ 0 for s ∈ tj ε α , tj+1 ε α , we have |F | ≤ ∂ 2 t ρν ∞ 2 n j=0 t j+1 ε α t j ε α 1 - ε α s -t j ∆t ε α s -t j ∆t 2 + ε α s -t j ∆t t j+1 -ε α s ∆t 2 ν 2 M e -sν V ds ≤ ∂ 2 t ρν ∞ 2 n j=0 t j+1 ε α t j ε α ε α s -t j ∆t ν 2 M e -sν V ds. (3.38)
Now, let us denote E n = ρν (t n , k) -ρn ν (k) and suppose that E 0 = 0. From (3.27) and (3.36) we have

1 - c 0 K E n+1 = b 0 K E n + 1 K n j=1 (c j E n+1-j + b j E n-j ) + ∆t 2 K F. (3.39)
In the sequel, to get an estimate for E n+1 , we will need two different bounds for (1 -c0 K ) -1 = K/(K -c 0 ). From the definition (3.26) of c 0 , we directly get

|c 0 | ≤ ∆t ε α 0 1 - ε α s ∆t ν 2 M e -sν V ds
from which we deduce the first estimate for

|K/(K -c 0 )| K K -c 0 ≤ K νM V - ∆t ε α 0 1 - ε α s ∆t ν 2 M e -sν V ds . (3.40) Now, using 0 ≤ 1 -ε α s ∆t ≤ 1 for s ∈ 0, ∆t ε α , we get |c 0 | ≤ ν 2 M ∆t ε α 0 e -sν ds V ≤ K -νM e -∆t ε α ν V so that we deduce the second estimate for |K/(K -c 0 )| K K -c 0 ≤ K νM e -∆t ε α ν -1 V . (3.41)
In (3.39) let us focus on the term containing F . From (3.38), we have

1 K F 1 -c0 K = F K -c 0 ≤ ∂ 2 t ρν ∞ 2 ∆t ε α 0 ε α s ∆t ν 2 M e -sν V ds |K -c 0 | + ∂ 2 t ρν ∞ 2 n j=1 t j+1 ε α t j ε α ε α s -t j ∆t ν 2 M e -sν V ds |K -c 0 | , (3.42) =: ∂ 2 t ρν ∞ 2 (T 1 + T 2 ).
To find an estimate for T 1 , we use (3.40) to write

T 1 ≤ ∆t ε α 0 ν 2 M e -sν(v) V ds - ∆t ε α 0 1 - ε α s ∆t ν 2 M e -sν V ds K - ∆t ε α 0 1 - ε α s ∆t ν 2 M e -sν V ds , (3.43)
that implies that the first line of (3.42) is bounded by ∂ 2 t ρν ∞ /2. Now, using (3.41), T 2 can be estimated as

T 2 ≤ t n+1 ε α ∆t ε α ν 2 M e -sν V ds |K -c 0 | -1 ≤ 1. ( 3
1 K b 0 E n + n j=1 (c j E n+1-j + b j E n-j ) ≤ 1 K ∆t ε α 0 ε α s ∆t ν 2 M e -s(ν+iεk•v) V E n ds + n j=1 t j+1 ε α t j ε α 1 - ε α s -t j ∆t ν 2 M e -s(ν+iεk•v) V E n+1-j ds + t j+1 ε α t j ε α ε α s -t j ∆t ν 2 M e -s(ν+iεk•v) V E n-j ds ≤ 1 K K - ∆t ε α 0 1 - ε α s ∆t ν 2 M e -νs V ds E n .
Hence, using (3.40), we get

|b 0 E n + n j=1 (c j E n+1-j + b j E n-j )| |K -c 0 | -1 ≤ E n . (3.46)
Eventually (3.39) gives, using (3.45) and (3.46),

E n+1 ≤ E n + ∂ 2 t ρν ∞ ∆t 2 . Hence E n+1 ≤ E n + ∂ 2 t ρν ∞ ∆t 2 .
As we supposed E 0 = 0, it writes

E n+1 ≤ ∂ 2 t ρν ∞ T ∆t, (3.47)
where T = N ∆t is the final time. This proves that the scheme in ρν is first order accurate uniformly in ε. Now, to prove that the scheme in f is first order uniformly in ε we consider the Duhamel formulation for f

f (t n+1 , k, v) = e -∆t ε α (ν+iεk•v) f (t n , k, v) + νM ∆t ε α 0 e -s(ν+iεk•v) ρν (t n+1 -ε α s, k)ds,
and with the same Taylor expansion as previously, we have

f (t n+1 , k, v) = νM (γ ρν (t n+1 , k) + β ρν (t n , k)) + νM ∆t 2 G + e -∆t ε α (ν+iεk•v) f (t n , k, v),
with β and γ defined in (3.33) and

G = - 1 2 ∆t ε α 0 1 - ε α s ∆t ε α s ∆t 2 ∂ 2 t ρν (t n+1 -ε α ξ 1 (s), k) + ε α s ∆t 1 - ε α s ∆t 2 ∂ 2 t ρν (t n+1 -ε α ξ 2 (s), k) e -s(ν+iεk•v) ds, (3.48)
where ξ 1 (s), ξ 2 (s) ∈ (0, ∆t/ε α ) . As previously, we can estimate |νM G|

|ν M G| ≤ νM ∂ 2 t ρν ∞ 2 ∆t ε α 0 1 - ε α s -t j ∆t ε α s -t j ∆t 2 + ε α s -t j ∆t ε α s -t j+1 ∆t 2 e -s(ν+iεk•v) ds ≤ νM ∂ 2 t ρν ∞ 2 ∆t ε α 0 e -sν ds ≤ C 1 -e -∆t ε α ν .
One can now look at the truncation error

E f n = f (t n , k, v) -f n (k, v) which satisfies E f n+1 = e -∆t ε α (ν+iεk•v) E f n + νM (γE n+1 + βE n ) + νM ∆t 2 G. (3.49)
The term νM (γE n+1 + βE n ) can be estimated by (using definition (3.33))

νM (γE n+1 + βE n ) ≤ νM ∆t ε α 0 e -sν ds E n+1 ≤ C 1 -e -∆t ε α ν ∆t,
where we used (3.47) for the last inequality. Hence, denoting

Λ n = max v =0,k,ε∈ (0,1 ] |E f n | we get from (3.49): Λ n+1 ≤ e -∆t ε α ν Λ n + C 1 -e -∆t ε α ν ∆t 2 + ∆t . With Λ 0 = 0, it finally comes that Λ n+1 ≤ C 1 -e -∆t ε α ν ∆t + ∆t 2 n j=0 e -j ∆t ε α ν ≤ C∆t.

General initial data.

In the previous subsection, we observed that the asymptotic scheme of (3.35) as ε → 0 writes

ρn+1 ν (k) -ρn ν (k) ∆t = -κ|k| α ρn+1 ν (k), ∀n ≥ 0,
and moreover, we have ρn+1 = ρn+1 ν for all n ≥ 0. But, for general initial data, since ρ0

ν is different from ρ0 , the asymptotic scheme is not the correct one because of this initial mismatch. Then when the initial data is not at equilibrium, the previous scheme is not AP and a suitable modification has to be done to recover the AP property. To do so, we modify (3.35) to ensure that the first step of the asymptotic scheme writes (when ε → 0)

ρ1 ν (k) -ρ0 (k) ∆t = -κ|k| α ρ1 ν .
Hence, we have to find an expression for the term ρ0 ν appearing in the first equation of (3.35) for n = 0. Starting from (3.32) for n = 0, we get, up to first order

ρ1 (k) = e -∆t ε α (ν+iεk•v) f 0 Nv,V + dρ 0 ν (k),
where d is given by

d = 1 -e -∆t ε α (ν+iεk•v) ν + iεk • v νM Nv,V . (3.50)
Up to an error of order ∆t, it is possible to replace ρ1 by ρ0 so that ρ0 ν becomes

ρ0 ν = 1 d ρ0 -e -∆t ε α (ν+iεk•v) f 0 Nv,V . (3.51)
This expression is injected in (3.35) for n = 0 to compute ρ1 ν . Then, for any ε, (3.51) is a consistent approximation and when ε goes to zero, ρ0 ν tends to ρ0 which ensures the AP property for any initial data.

Another defect in the previous scheme (3.35) is that, for a general initial data, it does not preserve the total mass. To overcome this, we add the following consistent step for f n+1

f n+1 (k = 0) ← f n+1 (k = 0) -f n+1 (k = 0) Nv,V M + f 0 (k = 0) Nv,V M. (3.52)
Eventually, for a general initial condition, we have the following proposition Proposition 3.4. We consider the scheme defined for all k and all time indices 1 ≤ n ≤ N, n∆t = T by (3.35)-(3.52) and for the first step by (3.50). This scheme has the following properties: (i) The scheme is first order in time and preserves the total mass. (ii) The scheme is AP: for a fixed ∆t, the scheme solves the anomalous diffusion equation when ε goes to 0

                 ρ1 ν (k) = 1 K -c 0 Â0 (t 1 , k) + b 0 d ρ0 -e -∆t ε α (ν+iεk•v) f 0 Nv,V f 1 (k, v) = e -∆t ε α (ν+iεk•v) f 0 (k, v) + νM γ ρ1 ν (k) + β ρ0 ν (k) f 1 (0, v) := f 1 (0, v) -f 1 (0, v) Nv,V M + f 0 (0, v) Nv,V M, (3.53
ρn+1 (k) -ρn (k) ∆t = -κ|k| α ρn+1 (k),
where κ is defined in (3.6). Remark 6. for a general initial data, the time derivatives at t = 0 are not bounded uniformly in ε. Hence, the uniform accuracy of the scheme is not a relevant property to be considered. However, the numerical tests tends to suggest that the consistency error is of order ∆t uniformly in ε.

Proof. To prove that the scheme is first order in time, we have to show that when ε is fixed

ρν (t 1 , •) -ρ1 ν (•) ∞ ≤ C ε ∆t,
where C ε is a constant when ε is fixed. We will then be able to use the proof of Prop. 3.3. At the continuous level, (3.32) gives

ρ(t 1 , k) = e -∆t ε α (ν+iεk•v) f0 (k, v) V + ∆t ε α 0 νM e -s(ν+iεk•v) V ρν (t 1 -ε α s, k)ds.
Denoting by

G(s) = t1-ε α s 0 ∂ t ρν (u)du and H = - ∆t 0 ∂ t ρ(u)du the remainders of the Taylor expansions of ρν (t 1 -ε α s, k) and ρ(t 1 , k), we have ρ(t 0 , k) = e -∆t ε α (ν+iεk•v) f0 (k, v) V + dρ ν (t 0 , k) + H + ∆t ε α 0 νM e -s(ν+iεk•v) V G(s)ds,
with d defined in (3.50). Hence, the error made by replacing ρν (t 0 , k) by its approximation (3.51) in the first iteration of (3.35) is of magnitude

1 d H + ∆t ε α 0 νM e -s(ν+iεk•v) V G(s)ds .
Eventually, the error done in the first iteration of (3.53) writes Moreover, since

ρν (t 1 , •) -ρ1 ν (•) ∞ ≤ 1 d b 0 K -c 0 H + ∆t ε α 0 νM e -s(ν+iεk•v) V G(s)ds (3.54) + F K -c 0 ∆t 2 ,
|b 0 | ≤ ∆t ε α 0 ν 2 M e -s(ν+iεk•v) V ds ,
an exact computation of this integral shows that when ε is fixed, the term b 0 /d is bounded. Assuming that the time derivatives of ρ and ρν are bounded, there exists a constant C ε depending a priori of ε such that ρν (t

1 , •) -ρ1 ν (•) ∞ ≤ C ε ∆t.
The proof of Prop. 3.3 can then be used to show that the scheme is of order 1.

To show that (3.53) enjoys the AP property, we freeze ∆t and let ε go to 0 in the expression of the scheme. The first iteration degenerates into

ρ1 ν (k) -ρ0 (k) ∆t = -κ|k| α ρ1 ν (k), ρ1 (k, v) = ρ1 ν (k)M,
then, when ε goes to zero ρ1 (k) = ρ1 ν (k) which ensures that the first iteration of the scheme enjoys the AP property. As the following iterations are the same as in Prop.

3.3, this proves (ii).

4. Numerical results. In this section, we present some numerical tests to validate the three schemes we proposed in Section 3. For simplicity, we denote by IS the the implicit scheme, MMS the micro-macro scheme and DS the scheme based on the Duhamel formulation of the kinetic equation. Eventually, the implicit Euler scheme for the anomalous diffusion equation (3.8) is denoted adiff. We use the equilibrium (1.2), the collision frequency (1.3) and the initial data

f 0 (x, v) = (1 + sin(πx))M 0 , (4.1)
for x ∈ [-1, 1] and v ∈ V = [-1, 1] and w max = 1 for the computation of κ in (3.6).

The normalization constant M 0 is chosen such that M Nv,V = 1 and ν 0 is fixed ν 0 = 1. We consider periodic conditions for the space, then the discretization used for the space variable uses N x + 1 points and is given by

x i = -1 + i∆x, 0 ≤ i ≤ N x , (4.2)
where ∆x = 2 Nx . This space discretization is linked to the Fourier modes we use for the discrete Fourier transform in the numerical computations, that are the integers k such that -N x /2 ≤ k ≤ N x /2. The uniform discretization of the velocities is done with N v = 2N ′ v + 1 points symmetrically distributed to ensure that the property vM Nv ,V = 0 is conserved at the discrete level. Denoting ∆v = 1/N ′ v , we set

v j = -1 + j∆v, 0 ≤ j ≤ 2N ′ v . (4.3)
Since the schemes do not require a refined velocity grid to enjoy the AP property, we fix N ′ v = 10 in all the following tests. All the tests are performed using a final time T = 1 and β = 1, the time step ∆t being specified in each case.

4.1. The implicit scheme (IS). In this part, we test the properties of the IS scheme. To highlight its AP character, we present in the left hand side of Fig. 4.1 the densities ρ(T, x) obtained by IS for ∆t = 10 -1 , N x = 32 and a range of ε. For this figure, we compare our results to the anomalous diffusion limit given by the adiff scheme with the same time step. It appears that the densities converge slowly to the anomalous diffusion limit when ε goes to 0. In the right part of Fig. 4.1, we plot the L ∞ norm (in space) of the error between ρ(T, x) given by IS and ρ adiff (T, x) given by adiff as a function of ε (in log scale). Then, we can compare our numerical results to the rate of convergence computed in Prop. 2.2. We observe that the numerical rate 0.334 is in a very good agreement with the theoretical rate d/(d + 1 + β) ≈ 0.333. Let us remark that we showed in Prop. 2.2 that theoretically, the convergence happens rather with rate ε min( d d+1+β , β d+1+β ) , but the terms (5.1)-(5.2) giving the O(ε β d+1+β ) do not arise numerically. Indeed, it is given by a divergent integral multiplied by a power of ε. As we use a bounded domain for v and we did not apply a specific treatment for these terms, this divergent integral is numerically finite and then it vanishes when ε becomes small. Left: For ∆t = 10 -1 , the densities given by IS for a range of ε and by adiff. Right: For ∆t = 10 -1 , the error between the densities given by IS and by adiff as a function of ε (log scale).

Micro-macro based scheme (MMS).

In this part, we test the properties of MMS. In the left hand side of Fig. 4.2 we present the densities ρ(T, x) obtained with the MMS scheme for a range of ε and we compare it to the anomalous diffusion limit given by adiff. We used N x = 32 points to discretize the space. We use Fourier variable for the equation on ρ (macro part), but the space derivatives of the micro part are discretized with a classical first order upwind scheme. Hence, a CFL condition has to respected so that ∆t = 10 -3 . When ε goes to zero, the densities obtained by MMS converges towards the density obtained by adiff. On the right hand side of Fig. 4.2, we plot the L ∞ norm (in space) of the error between ρ(T, x) obtained by MMS and ρ adiff (T, x) obtained by adiff (using ∆t = 10 -3 ), as a function of ε. We can observe that the convergence happens with a rate 0.337 which is very close to the theoretical value ε d d+1+β = 0.333. , the densities given by the MMS scheme for a range of ε and the anomalous diffusion limit given by the adiff scheme. Right: For ∆t = 10 -3 , the error between the MMS scheme and the DSA scheme as a function of ε (log scale).

The integral formulation scheme (DS).

In this section, we test the properties of the DS scheme. First of all, we remark that for ε = 1 the scheme appears to be of second order, as suggested by Fig. 4.3. For this test, we consider N x = 8 and we consider the L ∞ norm (in space) of the error between the densities given by DS for a range of ∆t ∈ [2 × 10 -3 , 0.5] and the density given by DS for ∆t = 10 -3 (which is the reference). Then, we study the AP character of the scheme, fixing ∆t = 10 -1 , considering a range of ε. As previously, on the left hand side of Fig. 4.4, the densities ρ(T, x) obtained by DS and by adiff (using ∆t = 10 -1 ) are plotted as a function of space with N x = 32. On the right hand side of Fig. 4.4, the L ∞ norm (in space) of the error between ρ(T, x) obtained by DS and ρ adiff (T, x) obtained by adiff (using ∆t = 10 -1 ) is plotted as a function of ε, for N x = 8. As in the previous cases, we observe that the numerical convergence rate is close to the theoretical one ε d d+1+β . Eventually, we illustrate that DS is first order in time uniformly in ε. Indeed, the left hand side of Fig. 4.5 displays the L ∞ error (in space)

E(ε, ∆t) = ρ ε ∆t (T, •) -ρ ε ∆t=10 -3 (T, •) ∞ , (4.4)
where ρ ε ∆t=10 -3 (T, x) is obtained by DS with ∆t = 10 -3 and ρ ε ∆t (T, •) is obtained by DS using larger ∆t. The error E(ε, ∆t) is computed for a range of ∆t ∈ [5 × 10 -3 , 0.5] and of ε ∈ [10 -7 , 1]. For this test, we fixed N x = 8. These errors are stratified with respect to ∆t, confirming the uniform accuracy of DS with respect to ε.

We now consider the following non equilibrium initial data and the same numerical parameters as in the last test. As previously, we also look at (4.4) (see the right hand side of Fig. 4.5). It turns out that even if the initial condition is out of equilibrium, the Duhamel based scheme is uniformly accurate. As a consequence, considering the maximum with respect to ε ∈ [10 -7 , 1] of E(ε, ∆t) enables to highlight that DS is of order 1 uniformly in ε (obviously, the same occurs for the well-prepared initial data).

f 0 (x, v) = (1 + sin(πx))e -|v| , (4.5)

Appendix.

In this section, we deal with the proof of Prop. 2.1.

The proof is a direct consequence of the following lemma: Lemma 5.1. Let f be the solution of (1.5) in the Fourier variable for a given initial condition f0 and let t > 0 be fixed. Then its modified density ρν satisfies formally

∂ t ρν (t, k) + |k| α κρ ν (t, k) = O ε d d+1+β + O ε β d+1+β
, with κ given by (1.7) and with the initial condition ρν (0, k) = ν f0 (k, v) and consider each line separately. First of all, let us remark that in (5.4), the integral in s can be directly computed and that the initial layer term Â0 in (5.1) is exponentially small. Let us denote I 1 the integral term appearing in (5.1) and show that is equivalent to -ε α ∂ t ρν (t, k) for small ε up to a rest of order O ε α+ αd d+2+β

. We rewrite I 1 as

I 1 = - t ε α 0 ν 2 M e -νs V t t-ε α s ∂ t ρν (ξ, k)dξ ds,
and we permute the integrals before computing exactly the integral in s to get

I 1 = - t 0 ∂ t ρν (t -ξ, k) νM e -ν ξ ε α V dξ + O(ε ∞ ).
Hence, up to an exponentially small error in ε, I 1 reads

I 1 = -ε α ∂ t ρν (t, k) -ε α t ε α
0 νM e -ν(v)ξ V (∂ t ρν (t -ε α ξ, k) -∂ t ρν (t, k)) dξ.

It remains to show that the second term appearing in I 1 is smaller that ε α . Indeed in the brackets we do the change of variables w = (ν 0 ξ) 1 d+2+β v to make appear a fractional time derivative. Then with the change of variables ξ → ε α ξ in the integral in ξ, I 1 writes (5.8)

I 1 = -ε α ∂ t ρν (t, k) + O (ε ∞ ) -ε α( 2d+2+β d+2+β ) M 0 ν
The exact computation of the integral in ξ leads to

I 2 2 = O ε d d+1+β

,

(5.9)

and we now focus on I 1 2 . In the brackets we do the change of variables v = ν0 |k|ε (ξν 0 ) -d+1+β d+2+β such that ξϕ = .

(5.12)

  in (3.10) leads to the result. For the AP character of the scheme (ii), letting ε to 0 with fixed ∆t in the first equation of (3.7) leads to ρn+1 ν -ρn ∆t = -κ|k| α ρn+1 ν and with (3.3) we have ρn = ρn ν , ∀n ≥ 1 when ε → 0. It implies that (3.7) degenerates into (3.8).

  e denotes any unitary vector. Here, the fractional diffusion limit clearly appears. We then inject (3.29) computed with discrete brackets • Nv,|wmax|≥|w|≥ε|k|/ν0 in (3.28) and eventually compute explicitly the integrals in s of (3.25)-(3.26) to get the following expression of b j and c j b

  ) with K, Â0 , b j , c j , γ and β defined in (3.21)-(3.30)-(3.31) and (3.34).

  the other terms in (3.39), we introduceE n = max j=0•••n ε∈ (0,1 ] |E j | .From the definition (3.30) and (3.31), we have for the first two terms of the right hand side of(3.39)

  ) with K, Â0 , b 0 , c 0 , γ, β and d defined in (3.21)-(3.30)-(3.31)-(3.34) and

  with K, b 0 , c 0 defined in (3.21)-(3.25)-(3.26) and F in (3.37) for n = 0. The proof of Prop. 3.3 gives estimates for 1/(K -c 0 ) and F/(K -c 0 ) in (3.41) and (3.45).

  Fig.[START_REF] Bensoussan | Boundary layers and homogenization of transport processes[END_REF].1. Left: For ∆t = 10 -1 , the densities given by IS for a range of ε and by adiff. Right: For ∆t = 10 -1 , the error between the densities given by IS and by adiff as a function of ε (log scale).

Fig. 4 . 2 .

 42 Fig.4.2. left: For ∆t = 10 -3 , the densities given by the MMS scheme for a range of ε and the anomalous diffusion limit given by the adiff scheme. Right: For ∆t = 10 -3 , the error between the MMS scheme and the DSA scheme as a function of ε (log scale).

Fig. 4 . 3 .

 43 Fig. 4.3. For ε = 1, the error between DS computed for ∆t = 10 -3 and DS computed for a range of ∆t ∈ [2 × 10 -3 , 0.5] (log scale).

Fig. 4 . 4 .

 44 Fig. 4.4. Left: For ∆t = 10 -1 , the densities given by DS for a range of ε and by adiff. Right: For ∆t = 10 -1 , the error between DS and the adiff as a function of ε (log scale).

Fig. 4 . 5 .

 45 Fig. 4.5. The error between DS computed for ∆t = 10 -3 and DS computed for a range of ∆t, as functions of ε (log scale). Left: equilibrium initial data (4.1). Right: non equilibrium initial data (4.5).

Vε α 0 ν 2 ν 2 ν 2 ν 2 M

 02222 / νM V . Furthermore, the density and modified density ρ and ρν satisfy formally ρ(t, k) = ρν (t, k) + O (ε α-γ ) , for all 0 < γ < α.Proof. (Lemma 5.1) The proof relies on the integral formulation (3.19) which gives the expression (3.20) after an integration in v ∈ V . Using the notations K and Â0 defined in (3.21), let us rewrite (3.20) at time t as followsK ρν (t, k) = Â0 (t, k) + t M e -νs V (ρ ν (t -ε α s, k) -ρν (t, k)) ds (M e -s(ν+iεk•v) -e -νs V (ρ ν (t -ε α s, k) -ρν (t, k)) ds (M e -s(ν+iεk•v) -e -νsV dsρ ν (t, k) e -νs V dsρ ν (t, k),(5.4)

|w| 1 d+2+β. ( 5 . 5 )ε α 0 ϕ 2 I 2 = 5

 1550225 d+2+β e -|w| d+2+β |w|≤ δ ∂ t ρν (t -ξ, k) -∂ t ρν (t, and • |w|≤ δ denotes the integration in w on the domain |w| ≤ δ. The term into brackets is integrable on R d and 2d+2+β d+2+β= 1 + d d+2+β ∈ (1, 2). Supposing that ρν is regular, the integral in ξ is well defined and is bounded by a constant for any fixed t. HenceI 1 = -ε α ∂ t ρν (t, k) + O ε α+ αd d+2+βNow, let us denote I 2 the term (5.2) and remark that it vanishes for k = 0. We want to show that it is of magnitude O ε α+ d d+1+β when ε goes to 0. For nonzero k we apply the change of variable w = ε|k|v ν(v) , such that v = to make the limit appear in the brackets. Then, with a Taylor expansion of ρν , I 2 writesI 2 = -ε α C t |w| d+α e -sϕ(1+iw•e) -e -sϕ |w|≥∆ t t-ε α s ∂ t ρν (ξ, k)dξ ds, ε|k| ν0 and where e denotes any unitary vector. As for I 1 , we permute the integrals before computing the integral in s to get the following expression of I ε 2α CI 1 2 + ε α C∂ t ρν (t, k)I 2 2 + O (ε ∞ ) , e -ξϕ(1+iw•e) 1 + iw • e -e -ξϕ |w|≥∆ (∂ t ρν (t -ε α ξ, k) -∂ t ρν (t, k)) dξ (

0 3 = t ε α 0 ν 2 MI 3 = - ε α |k| α ν 1 3 I 3 =

 3023133 |∂ t ρν (t -ε α ξ, k) -∂ t ρν (t, k)| |ξ| α d+1+β d+2+β dξ,we then apply the change of variables ξ → ε α ξ to getε α I 1 2 = O ε αd d+2+β, which combined with (5.9) in (5.6) leads toI 2 = O ε α+ d d+1+β + O ε α+ αd d+2+β = O ε α+ d d+1+β .(5.10) It remains to consider (5.3), let us denote I e -s(ν+iεk•v) -e -νs V ds, that writes after the change of variables w = ε|k|v ν(v) and an exact integration in s w • e) 2 |w|≥∆ + O(ε ∞ ), where we see the anomalous diffusion coefficient appear. Eventually, we get the following equivalent for I -ε α κ + O ε α+ β d+1+β , (5.11) with κ given by (1.7). Finally replacing the lines (5.1)-(5.2)-(5.3)-(5.4) with their equivalents leads to the result ∂ t ρν + |k| α κρ ν = O ε d d+1+β + O ε β d+1+β


Thomas Bonald
email: thomas.bonald@telecom-paristech.fr

Salah Eddine Elayoubi
email: salaheddine.elayoubi@orange.com

Yu-Ting Lin
email: yuting.lin@orange.com

A Flow-Level Performance Model for Mobile Networks Carrying Adaptive Streaming Traffic

This paper proposes a performance model for mobile networks carrying adaptive streaming traffic. The proposed model takes into account the flow dynamics in addition to the main parameters influencing the performance of adaptive streaming, such as the playout buffer and the video bit rates. We show how to compute several performance metrics like the average video bit rate, the deficit rate, defined as the probability of having an instantaneous throughput lower than the chosen video bit rate, and the average buffer surplus, related to the amount of data accumulated in the buffer. Considering the coexistence of multiple services and heterogeneous radio conditions that make the exact solution intractable, we propose a simple yet accurate approximation that is easy to integrate in the operator's dimensioning tools. Our numerical results investigate the performance trade-offs between the different parameters of the adaptive streaming service.

I. INTRODUCTION

According to Cisco forecast 2013 [START_REF]Cisco Visual Networking Index: Global Mobile Data Traffic Forecast Update[END_REF], video services will account for increasing proportion of data transmitted in future mobile networks and adaptive streaming will represent most of the video streaming services. Several implementations of adaptive streaming exist, such as Microsoft Smooth Streaming, Adobe HTTP Dynamic Streaming, Apple HTTP Live Streaming (HLS) and the MPEG Forum specified a standard Dynamic Adaptive Streaming over HTTP (MPEG-DASH). In this paper, we do not look on a specific protocol implementation but rather develop a general performance model that aims at helping operators assessing the quality of service perceived by their users and properly dimensioning their networks.

In [START_REF] Oyman | Quality of experience for HTTP adaptive streaming services[END_REF], authors show typical implementations for adaptive streaming, where the video files are divided into several chunks (segments) [START_REF]Information Technology -Dynamic Adaptive Streaming Over HTTP (DASH) -part 1: Media Presentation Description and Segment Formats[END_REF] [START_REF] Sodagar | The MPEG-DASH standard for multimedia streaming over the internet[END_REF] and the server has several encoded versions of each chunk. At the client side, a playout buffer stores the downloaded video segments that are not yet played. Even if this buffer is of finite size, we suppose in this paper that it is able to store the whole video files and that the video chunks are downloaded by the user at the maximum achievable throughput so that the amount of buffered chunks is maximized. We also consider small chunk sizes, ensuring an instantaneous adaptation of the video bit rate to the available throughput, when the number of active flows changes in the system. A traffic flow stands for a download task such as data transfer or HTTP video streaming.

Users are satisfied if they watch the video at a high coded rate and if the video playout is smooth, i.e. the playout buffer never gets empty. The Key Performance Indicators (KPIs) that directly influence the users' Quality of Experience (QoE) are the average video bit rate observed during the video session and the starvation probability, i.e. the probability that the buffer gets empty. While our flow-level model is able to compute the average video bit rate, the starvation probability computation needs a packet-level analysis as it depends on the behavior of the player in terms of prefetching policy and the detailed buffer state. As our objective is to provide simple models that can be used for mobile network dimensioning purposes, we provide two KPIs that are related to the starvation probability and can be computed using the flow-level model: the deficit rate defined as the probability that the instantaneous user throughput is lower than the chosen video bit rate, and the buffer surplus representing the average buffer variation during the video download.

Related literature

An important set of works considered flow-level modeling in mobile networks. Note that a flow represents for instance a file download or a video streaming session and that it may be subject to different radio conditions due to the position of the user in the cell and to throughput variations due to the dynamics of arrivals and departures of other users. In [START_REF] Bonald | Wireless downlink data channels: user performance and cell dimensioning[END_REF], a flow-level model has been proposed for elastic traffic in mobile networks assuming static users. This model is extended to mobile users in [START_REF] Bonald | Flowlevel performance and capacity of wireless networks with user mobility[END_REF] and to advanced radio features in [START_REF] Khlass | Flow-level performance of intrasite coordination in cellular networks[END_REF]. Another set of works considers the performance of real-time streaming services at the flow level [START_REF] Blaszczyszyn | Quality of Real-Time Streaming in Wireless Cellular Networks -Stochastic Modeling and Analysis[END_REF] and combining flow and packet levels [START_REF] Lin | Capacity dimensioning for realtime video services in wireless mobile networks[END_REF]. The integration of elastic and realtime streaming services is considered in [START_REF] Bonald | On performance bounds for the integration of elastic and adaptive streaming flows[END_REF] that provides performance bounds for both services. This integration has been also considered in [START_REF] Borst | Integration of streaming and elastic traffic in wireless networks[END_REF] where the streaming performance is evaluated through the blocking rate.

The above mentioned flow-level models focus on classical elastic and real-time streaming services and do not consider HTTP video streaming with buffering. A flow-level model is proposed in [START_REF] Xu | Impact of Flowlevel Dynamics on QoE of Video Streaming in Wireless Networks[END_REF] for HTTP video streaming with infinite buffers and KPIs like starvation probability are computed using a detailed buffer analysis. However, this work does not consider the video bit rate adaptability and the integration of elastic traffic.

Main contribution

We develop a flow-level model for the performance of adaptive streaming taking into account the main characteristics of this service like the presence of a buffer and the video bit rate limitations. We extend this model to consider heterogeneous radio conditions and a mix between streaming and elastic traffic. We then propose and validate an approximation model to facilitate the calculation of performance metrics. Finally, we demonstrate the performance tradeoffs of different maximum video bit rates.

Paper organization

The remainder of the paper is organized as follows. Section II describes the system and develops the flow level performance model for adaptive streaming traffic. In section III, we extends the model with the consideration of heterogeneous radio conditions and the integration of elastic and streaming services. The approximation model is introduced in section IV. In section V, numerical analysis shows the validation of approximation model and the impacts of different video bit rates. Section VI concludes the paper.

II. MODEL FOR ADAPTIVE STREAMING

A. System description

In this section, we consider a cellular network and focus on the performance of a typical cell as described in Fig. 1. We will begin, for the ease of understanding of the models, by a homogeneous scenario where all users have the same physical rate R corresponding to the average radio condition over the cell. Section III will show how to extend this model to multiple radio conditions and how to include elastic traffic.

Moreover, to facilitate the flow-level modeling, for the streaming system configuration, we assume that each flow has an infinite buffer size so that each leaves the system as soon as its corresponding video is fully downloaded. During its stay in the system, each flow continuously downloads the video chunks according to the allocated resource. The streaming chunk size is assumed to be significantly small, enabling users to adapt their video bit rates at any time.

B. Markov model

We make the classical assumption that traffic flows arrive according to a Poisson process with intensity λ. We only consider adaptive streaming services in this section, the mix with other services is considered later. The video durations are assumed to be independent and exponentially distributed with mean T .

Let X(t) be the number of flows at time t. When X(t) = x, which means that there are x flows served in the system, the flow departure rate, µ(x), can be expressed as

µ(x) = ϕ(x) σ(x) , (1)
where ϕ(x) stands for the physical throughput allocated to all UEs of the cell at state x and σ(x) stands for the remaining flow size at state x. As we consider an infinite buffer size at the user side, users can fully utilize the throughput allocated to them, occupying thus the whole throughput, ϕ(x) = R.

As for the flow size, as we consider a small chunk size and an instantaneous adaptation to the radio conditions, the video bit rate depends on state x only and not on the history of the system. Furthermore, as the video duration is assumed exponential, the memoryless property applies and the flow size at state x is also exponentially distributed with mean

σ(x) = v(x)T, ( 2
)
where v(x) is the video bit rate at state x. X(t) is thus a Markov process whose transition rates depend on the selection of video bit rate. We will show in the following sections how the selected video bit rate is computed and how to derive the performance model in different scenarios of video bit rate limitations.

C. No video bit rate limitation

We first assume that number of available video bit rates is not limited so that the bit rate assigned to each flow is exactly equal to the instantaneous throughput, γ(x), that the user gets. This corresponds to an extremely adaptive streaming service where very low and very high video bit rates are allowed. Then, we have

v(x) = γ(x) = R x . ( 3
)
With equation ( 3), the flow departure rate then can be expressed as

µ(x) = ϕ(x) v(x)T = R T × R/x = x T . ( 4
)
The system can be easily shown to have a product form and the stationary distribution is

π(x) = π(0) ρ x x! , ( 5
)
where π(0) = e -ρ and ρ = λT . The number of users in the system evolves here like realtime communications as flows are always admitted in the system and the video duration is independent of the system state. Indeed, the buffer never fills as the user always watches the video with a rate equal to its available throughput, and starvation never happens as the video bit rate can go to zero without a playout interruption.

D. Video bit rate limitation

We now consider a more realistic case where minimal v min and maximal v max video bit rates are setup. However, we still assume a continuous set of video bit rates between these limits, although our model is sufficiently general to consider a discrete set. Therefore, the bit rate can be expressed as

v(x) = max ( min ( γ(x), v max ) , v min ) . ( 6
)
The flow departure rate becomes

µ(x) = ϕ(x) v(x)T = R max ( min ( γ(x), v max ) , v min ) T . ( 7
)
When the available throughput γ(x) = R x is larger than v max , the buffer starts filling as the download capacity is larger than the playout rate. The flow behaves thus as an elastic one as there is no limitation on the buffer size. On the other hand, when this available throughput is smaller than v min , the behavior is also elastic but with possible starvation if the buffer empties due to the constant playout rate. Between v max and v min , the flow behaves as a real time one and the buffer occupation remains constant. With the reversibility property [START_REF] Bonald | Network Performance Analysis[END_REF], the stationary distribution π(x) is

π(n) = λ n ∏ n x=1 µ(x) π(0), (8)
where π(0) can be obtained by calculating

∑ ∞ x=0 π(x) = 1.

E. Stability condition

If v min = 0, the system is always stable as the average flow duration is always equal to T . When v min > 0, the system may become unstable for a large number of arrivals. The flow arrival rate should be smaller than the max flow departure rate, leading to the following stability condition:

λ < R v min T .

F. KPIs definition

To evaluate the performance of adaptive streaming service, we propose three key performance indicators, mean video bit rate, deficit rate and buffer surplus. All of them are defined based on the stationary distribution π(x).

1) Mean video bit rate: The mean video bit rate stands for the average bit rate that a user experiences while watching the video. When the mean video bit rate is high, user has a better video experience. Due to the observers' paradox [START_REF] Bonald | Network Performance Analysis[END_REF], users see a biased distribution of the state, proportional to xπ(x) in state x. Thus, we define the overall mean video bit rate as

v = ∑ x:x>0 xπ(x) x v(x), (9)
with x = ∑ x>0 xπ(x).
2) Deficit rate: A popular QoE indicator used to evaluate streaming performance is the starvation probability [START_REF] Xu | Impact of Flowlevel Dynamics on QoE of Video Streaming in Wireless Networks[END_REF]. Even if starvation happens only when the video bit rate is larger than the instantaneous throughput, the latter condition is not a sufficient condition for starvation as the buffer may counteract the impact of short periods of low throughput. The computation of the starvation probability has to take into account the memory of the system by introducing the buffer size in the Markovian analysis, as in [START_REF] Xu | Impact of Flowlevel Dynamics on QoE of Video Streaming in Wireless Networks[END_REF].

Here, we introduce and examine the performance of a metric called deficit rate. The deficit rate is equal to the probability an ongoing flow sees its instantaneous throughput lower than its chosen bit rate. As the rate adaptation is assumed to occur instantaneously in reaction to the variations of the observed throughput, the deficit rate is defined by the probability that the instantaneous throughput, γ(x) = ϕ(x)

x , is smaller than v min . The overall deficit rate is defined by weighting the stationary distribution at different state x with the number of flows:

D = P{γ(x) < v min } = ∑ x:x>0 xπ(x)
x {γ(x)<vmin} , [START_REF] Bonald | On performance bounds for the integration of elastic and adaptive streaming flows[END_REF] where stands for the indicator function.

3) Buffer surplus: We also introduce another performance metric called buffer surplus, which represents the average relative buffer variation of each flow. It is expressed as

B = ∑ x:x>0 xπ(x) x ( γ(x) -v(x) v(x) ) . ( 11
)
When γ(x) > v(x), the video contents accumulates in the buffer. When γ(x) < v(x), then user starts to consume the video packets stored in the buffer. When B > 0, the video buffer of each flow increases on average. We shall see that B = 0 is a good objective for dimensioning the maximum traffic load.

III. MODEL EXTENSION

Section II proposed performance metrics for adaptive streaming and showed how to compute them when only adaptive streaming flows share the capacity. Here, we generalize the model to heterogeneous radio conditions and discuss the performance of these flows in the presence of elastic traffic. We also show how to integrate opportunistic scheduling in the model.

A. Heterogeneneous radio conditions and mixed service

Based on the 3GPP LTE-A standards [START_REF] Dahlman | 4G: LTE/LTE-Advanced for Mobile Broadband[END_REF], users with various positions have different discrete CQI (Channel Quality Indicator), for example from CQI-1 to CQI-15. Therefore, traffic flows can be separated into several classes, i ∈ I = {1, • • • }. Each class has a radio condition R i and accounts for a proportion of flow arrivals, p i , with ∑ i∈I p i = 1. We also consider an extended setting where streaming flows share the cell capacity with elastic flows. Therefore, the system can be represented as a network of two groups of coupled processor-sharing queues. Flows in queue e, i correspond to the elastic traffic with radio condition, R i , and flows in queue s, i correspond to the streaming ones with R i . The flow arrival rates at each queue are calculated as λ e,i = q e p i λ, λ s,i = q s p i λ, where q e and q s stand for the proportions of elastic and streaming arrival rates and q e + q s = 1.

Applying the same concept to formulate the departure rate as equation ( 4), the flow departure rate of elastic and adaptive streaming services can be expressed as

µ e,i (x) = ϕ e,i (x) σ , µ s,i (x) = ϕ s,i (x) v i (x)T , ( 12
)
where x = (x e,1 , • • • , x e,i , x s,1 , • • • , x s,i ) represents the number of flows for each service and σ is the mean flow size of elastic data. ϕ e,i (x) and ϕ s,i (x) stand for the allocated wireless resources for each class. Assuming the Round-Robin scheduling, these capacity shares are

ϕ e,i (x) = x e,i R i |x| , ϕ s,i (x) = x s,i R i |x| , ( 13
)
where |x| = ∑ i (x e,i +x s,i ). Here, the video bit rate, v i (x), is chosen depending on the instantaneous throughput as follows,

v i (x) = max ( min ( γ s,i (x), v max ) , v min ) . ( 14
)
where γ s,i (x) = ϕs,i(x) xs,i = Ri |x| . With the formulation above, note that the balanced property in [START_REF] Bonald | On performance bounds for the integration of elastic and adaptive streaming flows[END_REF] is not valid as µ s,i (xe e,j )µ e,j (x) ̸ = µ e,j (x -e s,i )µ i (x), where i, j ∈ I and e i is a vector with only 1 user belonging to class i. The Markov chain is not reversible and we have to solve the balance equations (by truncating the state space and inversing the corresponding matrix) for the stationary distribution π(x): ∀x,

( ∑ j=e,s ∑ i λ j,i + µ j,i (x) ) π(x) (15) = ∑ j=e,s ∑ i ( λ j,i π(x -e j,i ) + µ j,i (x + e j,i )π(x + e j,i )
) ,

B. Stability condition

We begin by assessing the stability region of the system. When the system approaches the stability limit, the video bit rate of streaming users decreases to v min . The system behaves like a system of elastic traffic. Stability holds only when the sum of offered loads for both services is less than 1, that is:

∑ i ( λq e p i σ R i + λq s p i v min T R i ) ≤ 1.
We deduce the maximum flow arrival rate:

λ max = ( ∑ i R i q e p i σ + q s p i v min T ) -1 . ( 16
)

C. KPIs definition

We extend the adaptive streaming related KPIs to this general case and define a KPI for the performance of elastic traffic.

1) Mean video bit rate:

We define the overall mean video bit rate of general model, v, by summing up all the classes and weighting by the flow number of each class at state x and the mean video bit rate of class i, vi as

v = ∑ x:|xs|>0 π(x) xs ∑ i x s,i v i (x), (17) vi
= ∑ x:xs,i>0 π(x) xs,i x s,i v i (x), (18)
where

xs = ∑ x π(x)|x s |, xs,i = ∑ x x s,i π(x) and |x s | = ∑ i x s,i .
2) Deficit rate: By using the same concept, we define the overall deficit rate of multiple class model, D and the deficit rate of class i, D i as

D = ∑ x:|xs|>0 π(x) xs ∑ i x s,i {γs,i(x)<vmin} , ( 19
)
D i = ∑ x:xs,i>0 π(x) xs,i x s,i {γs,i(x)<vmin} , ( 20
)
where is the indicator function equal to 1 when the condition is satisfied, otherwise the indicator function will become 0.

3) Buffer surplus: Same as the previous concept, we define the overall buffer surplus of multiple class model, B and the buffer surplus of class i, B i as

B = ∑ x:|xs|>0 π(x) xs ∑ i x s,i ( γ s,i (x) -v i (x) v i (x) ) , ( 21
)
B i = ∑ x:xs,i>0 π(x) xs,i x s,i ( γ s,i (x) -v i (x) v i (x) ) . ( 22
)
4) Average elastic throughput: The average elastic throughput is chosen as the performance metric for elastic flows, with γe and γe,i standing for the overall metric and mean throughput for each class i,

γe = ∑ x:|xe|>0 π(x) xe ∑ i x e,i ϕ e,i (x) x e,i , (23)
γe,i =

∑ x:xe,i>0 π(x) xe x e,i ϕ e,i (x) x e,i , (24)
where

xe = ∑ x π(x)|x e |, xe,i = ∑ x x e,i π(x) and |x e | = ∑ i x e,i .

D. Impact of opportunistic scheduling

The capacity shares in equation ( 13) are computed assuming a round robin scheduling. Now consider a channel aware scheduling, like the proportional fair scheduler, that operates at the fast fading time scale [START_REF] Tse | Multiuser diversity in wireless networks[END_REF]. In this case, when there are |x| flows that are active in the LTE cell, the throughput of a user of radio condition i that gets a proportion Φ of the cell resources is equal to ΦR i G(|x|), where G(|x|) is the opportunistic scheduling gain that depends on many parameters such as the channel model, the receiver and the Multiple Input Multiple Output (MIMO) scheme [START_REF] Combes | Cross-layer analysis of scheduling gains: Application to lmmse receivers in frequency-selective rayleigh-fading channels[END_REF]. Note that, contrary to real time streaming that does not benefit from the opportunistic scheduling gain due to its stringent delay constraints, http streaming with buffering takes advantage of this type of scheduling like elastic traffic. The performance model proposed above can thus be extended to the opportunistic scheduling case by introducing the fast fading gain, G(|x|) into the formulation of allocated resource in Eq. ( 13) and video bit rate in Eq.( 14). Note that, G(|x|) = 1 when round-robin scheduling is applied.

IV. APPROXIMATION

Section III introduces the general model with heterogeneous radio conditions and elastic flows. However, the complexity of solving numerically the balance equations mentioned in Eq. ( 15) increases exponentially with the number of classes, making the resolution of the stationary distribution very difficult. In this section, we propose an approximation model simplifying multiple queues into two coupled processor-sharing queues for elastic and streaming flows respectively.

A. Markov model

In our approximation, we make use of the intuition indicating that the contribution of a radio class to the load of the system is proportional to the ratio of the data volume it generates to the throughput it sees:

• For elastic traffic, class-i users generate a traffic volume of λ e p i and observe a throughput of R i , their contribution to the load of the elastic traffic queue, λepi Ri , then can be used to calculate the equivalent throughput R e .

• For the adaptive streaming traffic, class-i users generate a volume proportional to the video bit rate, i.e. that depends on the state of the system, |x| = x s + x e . The load it generates in the streaming queue is proportional to λspivi(x)

Ri

, where we use to calculate the equivalent throughput, R s (x). As we mentioned, system is approximated as the one carrying only two equivalent classes of adaptive streaming and elastic traffic and their flow departures are calculated as

μe (x) = φe (x) σ , μs (x) = φs (x) v min T , (25)
where x = (x e , x s ) denotes the number of elastic and streaming flows and the allocated resources for elastic and streaming are expressed as φe (x) = Re xe+xs , φs (x) = Rs(x) xe+xs with equivalent physical throughput

R e = ( ∑ i p i R i ) -1 , R s (x) = ( ∑ i p i α i (x) R i ) -1 . ( 26
)
Moreover, α i (x) is calculated by Eq. ( 14) as

α i (x) = v i (x) v min = max ( min ( γ s,i (x), v max ) , v min ) v min (27)
With the flow departure rate mentioned above and flow arrival rate shown as λ e = p e λ and λ s = p s λ, the approximated stationary distribution of π(x) can be obtained using the same concept of balance equations shown in Eq. [START_REF] Tse | Multiuser diversity in wireless networks[END_REF]. However, the complexity is much lower. Note that the stability condition of this approximation is the same as that gotten in Eq.( 16). Approximation also holds when applying opportunistic scheduling; it is thus sufficient to multiply the throughput by the state dependent scheduling gain G(|x|).

B. KPIs definition

As for the approximated performance metrics, the mean video bit rate, vs , the deficit rate, Ds and the buffer surplus, Bs for adaptive streaming traffic and the mean throughput for elastic traffic, γe , can be computed with the newly calculated stationary distribution π(x) as

vs = ∑ x:xs>0 x s π(x) xs ∑ i β i (x) β(x) v i (x), (28)

Ds = ∑

x:xs>0

x s π(x) xs

∑ i β i (x) β(x) {γs,i<vmin} , (29)

Bs = ∑

x:xs>0

x s π(x) xs

∑ i β i (x) β(x) ( γ s,i (x) -v i (x) v i (x) ) , (30) γe =

∑

x:xe>0

x e π(x) xe

ϕ e (x) x e = ∑ x π(x)ϕ e (x) xe , ( 31
)
where

β i (x) = p i α i (x) R i , β(x) = ∑ i β i (x), γ s,i (x) = R i |x| , (32)
and βi(x) β(x) represents the fraction of load volume of classi streaming users when there are x users in the system. In addition, xs = ∑ x:xs>0 x s π(x) and xe = ∑ x:xe>0 x e π(x) stand for the average number of streaming and elastic calls in the cell. It is also worth mentioning that the approximation model can also predict all the metrics for each class i.

V. NUMERICAL ANALYSIS

In our numerical analysis, we first validate the proposed approximation model, showing that the streaming and elastic performance with heterogeneous radio conditions can be accurately modeled by two classes. We then illustrate how to use our model for network dimensioning and study the impact of the video bit rates on the performance.

A. Validation of the approximation

To make the exact model tractable, we consider two radio conditions only representing the cell center and the cell edge, respectively. Therefore, there are four classes for the general model and two for the approximate model. Then we examine the performance of the approximate model by setting the parameters as (R C , R E ) = (10, 4)Mbps, (p C , p E ) = ( 12 , 1 2 ) for cell center and cell edge, (v min , v max ) = (0.5, 2), T = 10s, (p e , p s ) = ( 12 , 1 2 ), σ = 5Mbits. The simulation results in Fig. 2 show the four defined metrics for both models when the offered traffic increases. These figures demonstrate that the predictions of the approximate model are very close to the ones obtained by the exact model for four different metrics, also for the metrics of each class.

B. System dimensioning

Based on the results, operators can obtain the proper traffic intensity by setting some QoS constraints, which might be any combination of the performance metrics we defined. Take the buffer surplus, B > 0 as an example of QoS constraint in Fig. 2c, the traffic load should be lower than 0.655λ max . Otherwise, the average buffer surplus is smaller than zero, meaning the starvation events happen very often. Likewise, more QoS constraints can be considered together.

C. Performance trade-offs of v max configuration

After having validated the approximation using a simple radio model, we use it in the following scenario with realistic radio conditions based on measurement data from a 4G network in a large European city, with an average cell radius of 350 meters. The concerned frequency band is LTE 1800 MHZ. Figure 3 shows the measured probability distribution function of the Channel Quality Indicator (CQI) obtained from base station measurements collected using an O&M tool. Each CQI is associate to an MCS, determining its spectral efficiency. Using the CQI-MCS association figures of [START_REF] Bouguen | LTE et les réseaux 4G[END_REF] and considering a bandwidth of 10 MHZ, the corresponding harmonic capacity of an LTE cell is computed as equal to R e = 16.82 Mbps. Other traffic-related parameters are configured as T = 20s, σ = 5Mbits and v min = 0.5Mbps. Based on the stability condition, we have λ max = 2.189 flows/s. As of the opportunistic scheduling gain, we make use of the scheduling gain calculated in [START_REF] Wang | System level analysis of lte-advanced: with emphasis on multi-component carrier management[END_REF] for a MIMO 2 × 2 LTE system and an Additive White Gaussian Noise (AWGN) channel and that converges to G(∞) = 1.7 starting from a number of active uses in the LTE cell equal to 15.

In Fig. 4, four performance metrics are shown with different v max configurations. It can be observed that when v max increases, vs and D s increases, but B s and γe decreases, meaning that in the case of v min = 0.5Mbps, decreasing v max can benefit vs , D s and B s regardless of the trade-offs of γe reduction. We can also observe that the deficit rate is not highly influenced when v max is large. Therefore, we believe that buffer surplus is a better metric than deficit rate for dimensioning purposes.

VI. CONCLUSIONS

In this paper, we develop a flow-level model to describe the performance of adaptive streaming in wireless mobile networks carrying elastic traffic. We propose an approximate model which reduces the complexity of solving system metrics with heterogeneous conditions. Several performance metrics like mean video bit rate, deficit rate, buffer surplus and elastic mean throughput are taken as the main KPIs. It is shown that good buffer surplus and mean elastic throughput can be obtained by properly controlling the maximum video bit rate. In the future, the impacts of other system parameters such as the chunk duration and the buffer size need to be examined.

Fig. 1 :

 1 Fig. 1: System Model for Mobile Wireless Networks

  (a) Exact versus approximate model: mean video bit rate. (b) Exact versus approximate model: deficit rate. (c) Exact versus approximate model: buffer surplus. (d) Exact versus approximate model: mean elastic throughput.

Fig. 2 :

 2 Fig. 2: Validation of approximation model

Fig. 3 :

 3 Fig. 3: Measured CQI probability distribution function on a live LTE network.

  (a) Mean streaming video bit rate (b) Streaming deficit rate (c) Streaming buffer surplus (d) Mean elastic throughput.

Fig. 4 :

 4 Fig. 4: System performance with different v max settings.

ACKNOWLEDGEMENT

This work has been carried out in the framework of IDEFIX project, funded by the ANR under the contract number ANR-13-INFR-0006.