Megala Anandan
email: megalaa@iisc.ac.in

Benjamin Boutin
email: benjamin.boutin@univ-rennes.fr

N Crouseilles
email: nicolas.crouseilles@inria.fr

HIGH ORDER ASYMPTOTIC PRESERVING SCHEME FOR LINEAR KINETIC EQUATIONS WITH DIFFUSIVE SCALING

Keywords: collisional kinetic equation, diffusive scaling, high order Runge-Kutta schemes, asymptotic preserving property . 82C40, 85A25, 65M06, 65L04, 65L06

High order asymptotic preserving scheme for linear kinetic equations with diffusive scaling

1. Introduction. In this work, we are concerned with the numerical approximation of linear kinetic transport equations in a diffusive scaling. Such models are widely used in applications such as rarefied gas dynamics, neutron transport, and radiative transfer. Due to the presence of a small parameter ϵ (which is the normalized mean free path of the particles), standard schemes suffer from a severe restriction on the numerical parameters, making the simulations very costly. In the last decades, the so-called Asymptotic-Preserving (AP) schemes have been proposed to make possible the numerical passage between the micro and macro scale [START_REF] Jin | Efficient asymptotic-preserving (ap) schemes for some multiscale kinetic equations[END_REF][START_REF] Jin | Asymptotic preserving (ap) schemes for multiscale kinetic and hyperbolic equations: a review[END_REF]. Indeed, these AP schemes are uniformly stable and degenerate when ϵ → 0 to a scheme which is consistent with the asymptotic diffusion model. This makes them very attractive to deal with multi-scale phenomena as an alternative to domain decomposition approaches.

The goal of this work is to design high order in time AP schemes for collisional kinetic equations in the diffusive scaling. Several works can be found in the literature on this topic [START_REF] Jin | Uniformly accurate diffusive relaxation schemes for multiscale transport equations[END_REF][START_REF] Jin | Efficient asymptotic-preserving (ap) schemes for some multiscale kinetic equations[END_REF][START_REF] Jin | Asymptotic preserving (ap) schemes for multiscale kinetic and hyperbolic equations: a review[END_REF][START_REF] Klar | Asymptotic-induced domain decomposition methods for kinetic and drift diffusion semiconductor equations[END_REF][START_REF] Klar | Asymptotic-induced scheme for nonstationary transport equations in the diffusive limit[END_REF][START_REF] Klar | Numerical passage from radiative heat transfer to nonlinear diffusion models[END_REF][START_REF] Klar | Numerical passage from radiative heat transfer to nonlinear diffusion models[END_REF][START_REF] Naldi | Numerical schemes for kinetic equations in diffusive regimes[END_REF][START_REF] Dimarco | Asymptotic Preserving Monte Carlo methods for transport equations in the diffusive limit[END_REF][START_REF] Lafitte | Asymptotic-preserving projective integration schemes for kinetic equations in the diffusion limit[END_REF][START_REF] Larsen | Asymptotic solution of neutron transport problems for small mean free paths[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] Lemou | Relaxed micro-macro schemes for kinetic equations[END_REF][START_REF] Peng | Stability-enhanced AP IMEX-LDG schemes for linear kinetic transport equations under a diffusive scaling[END_REF]. Our work is based on a micro-macro decomposition as introduced in [START_REF] Lemou | A New Asymptotic Preserving Scheme Based on Micro-Macro Formulation for Linear Kinetic Equations in the Diffusion Limit[END_REF] where the unknown f of the stiff kinetic equation is split into an equilibrium part ρ plus a remainder g. A micro-macro model (equivalent to the original kinetic one) satisfied by ρ ad g can be derived. This micro-macro strategy turns out to be the starting point of several numerical approximations in phase space (using particles method, Discontinuous Galerkin method or low rank approximation [START_REF] Crouseilles | An asymptotic preserving scheme based on a micro-macro decomposition for collisional Vlasov equations: diffusion and high-field scaling limits[END_REF][START_REF] Crestetto | Asymptotically complexity diminishing schemes (ACDS) for kinetic equations in the diffusive scaling[END_REF][START_REF] Jang | High order asymptotic preserving DG-IMEX schemes for discrete-velocity kinetic equations in a diffusive scaling[END_REF][START_REF] Ding | Dynamical low-rank integrator for the linear boltzmann equation: error analysis in the diffusion limit[END_REF][START_REF] Einkemmer | An asymptotic-preserving dynamical low-rank method for the multi-scale multi-dimensional linear transport equation[END_REF]). In addition, a suitable first order semi-implicit time discretization of the micro-macro model is used as in [START_REF] Lemou | A New Asymptotic Preserving Scheme Based on Micro-Macro Formulation for Linear Kinetic Equations in the Diffusion Limit[END_REF] for which however the asymptotic diffusion equation is solved explicitly. This drawback is overcome following [START_REF] Lemou | Relaxed micro-macro schemes for kinetic equations[END_REF][START_REF] Crouseilles | An asymptotic preserving scheme based on a micro-macro decomposition for collisional Vlasov equations: diffusion and high-field scaling limits[END_REF][START_REF] Crestetto | Asymptotically complexity diminishing schemes (ACDS) for kinetic equations in the diffusive scaling[END_REF] in which the AP scheme degenerates into an implicit treatment of the diffusion equation. This improvement enables to get a numerical scheme which is asymptotically free from the usual parabolic condition.

The derivation of high order in time AP schemes for stiff kinetic problem has been performed by several authors [START_REF] Dimarco | Exponential Runge-Kutta methods for stiff kinetic equations[END_REF][START_REF] Dimarco | Asymptotic Preserving Implicit-Explicit Runge-Kutta Methods for Nonlinear Kinetic Equations[END_REF][START_REF] Dimarco | Implicit explicit linear multistep methods for stiff kinetic equations[END_REF][START_REF] Albi | Implicit-explicit multistep methods for hyperbolic systems with multiscale relaxation[END_REF][START_REF] Jang | High order asymptotic preserving DG-IMEX schemes for discrete-velocity kinetic equations in a diffusive scaling[END_REF] using the so-called high order IMEX methods [START_REF] Carpenter | Additive runge-kutta schemes for convection-diffusionreaction equations[END_REF][START_REF] Ascher | Implicit-explicit runge-kutta methods for time dependent partial differential equations[END_REF][START_REF] Pareschi | Implicit-explicit runge-kutta methods and applications to hyperbolic systems with relaxation[END_REF]. In this work, a family of high order IMEX schemes is proposed for linear collisional kinetic equations in the diffusive scaling, which degenerates when ϵ → 0 to a high order IMEX scheme for the diffusion equation. From the first order semi-implicit AP numerical scheme [START_REF] Crouseilles | An asymptotic preserving scheme based on a micro-macro decomposition for collisional Vlasov equations: diffusion and high-field scaling limits[END_REF], the family of high order schemes proposed in this work is obtained using globally stiffly accurate high order IMEX Runge-Kutta methods, namely type A and type CK [START_REF] Dimarco | Asymptotic Preserving Implicit-Explicit Runge-Kutta Methods for Nonlinear Kinetic Equations[END_REF][START_REF] Jang | High order asymptotic preserving DG-IMEX schemes for discrete-velocity kinetic equations in a diffusive scaling[END_REF].

In addition to the standard diffusion scaling, we also consider two other examples that enter in our framework. First, we consider a modification of the collision operator that enables to derive a transport-diffusion asymptotic model [START_REF] Jin | Uniformly accurate diffusive relaxation schemes for multiscale transport equations[END_REF][START_REF] Jang | High order asymptotic preserving DG-IMEX schemes for discrete-velocity kinetic equations in a diffusive scaling[END_REF]. Second, we discuss half moments micro-macro decomposition which naturally incorporates the incoming boundary conditions [START_REF] Lemou | Micro-Macro Schemes for Kinetic Equations Including Boundary Layers[END_REF].

Lastly, we address the space discretization in order to get a fully high order solver of the kinetic equation. High order space approximation based on finite difference methods is considered. Staggered or non-staggered strategies are proposed to achieve high order accuracy in space.

The paper is organized as follows. First in Section 2, the kinetic and asymptotic diffusion models are introduced. Then in Section 3, high order time integrators (using globally stiffly accurate IMEX Runge-Kutta temporal discretization) are proposed, and their AP property in the diffusive limit is addressed in Section 4. Section 5 is devoted to the space approximation. In Section 6, we discuss some extensions to other collision operators and to half moments. In Section 7, numerical results are presented, illustrating high order accuracy and the main properties of the schemes.

Kinetic equation, diffusion limit and micro-macro decomposition.

In this section, we introduce the kinetic model in the diffusive scaling, and recall the asymptotic limit. Then, the micro-macro decomposition is performed to derive the micro-macro model which serves as a basis for the numerical developments.

2.1. Linear kinetic equation with diffusive scaling. Let Ω ⊂ R d be the position space and V ⊆ R d be the velocity space with measure dµ(v). We consider the linear kinetic equation with diffusive scaling, (2.1)

∂ t f + 1 ϵ v • ∇ x f = 1 ϵ 2 Lf, (t, x, v) ∈ R + × Ω × V where f (t, x, v) ∈ R is the distribution function (depending on time t ∈ R + , space x ∈ Ω ⊂ R d and velocity v ∈ V ⊂ R d
) and ϵ > 0 measures the dimensionless mean free path of particles or the inverse of relaxation time. We consider the initial condition,

(2.2) f (0, x, v) = f init (x, v), (x, v) ∈ Ω × V
and boundary conditions are imposed in space. In this work, we will consider periodic boundary conditions or inflow boundary conditions. The linear collision operator L in (2.1) acts only on the velocity dependence of f , and it relaxes the particles to an equilibrium M (v) which is positive and even. We denote for all velocity dependent distribution functions h,

(2.3) ⟨h⟩ V = V h(v) dµ V M (v) dµ .
In particular, we obtain ⟨M ⟩ V = 1 and ⟨vM ⟩ V = 0. Further, the operator L is nonpositive and self-adjoint in L 2 V, M -1 dµ , with the following null space and range:

(2.4)

N (L) = {f : f ∈ Span (M )}, R(L) = (N (L)) ⊥ = {f : ⟨f ⟩ V = 0}.
This manuscript is for review purposes only.

Therefore, L is invertible on R(L) and we denote its pseudo-inverse by L -1 . We also assume that L is invariant under orthogonal transformations of R d .

Diffusion limit.

In the limit ϵ → 0, it is seen from (2.1) that f → f 0 where f 0 belongs to N (L). Thus, f 0 = ρ(t, x)M where f 0 solves Lf 0 = 0 and where the limiting density ρ is the solution of the asymptotic diffusion equation. To derive the diffusion equation, a Chapman-Enskog expansion has to be performed to get

f = f 0 + ϵL -1 (vM ) • ∇ x ρ + O(ϵ 2
). Integrating with respect to the velocity variable enables to get the diffusion limit

(2.5) ∂ t ρ -∇ x • (κ∇ x ρ) = 0 with κ = -v ⊗ L -1 (vM ) V > 0.
2.3. Micro-macro decomposition. In this part, we derive a micro-macro model which is equivalent to (2.1), and this is the model that will be discretized in the next sections. First, we consider the standard micro-macro decomposition of the unknown f [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 | Relaxed micro-macro schemes for kinetic equations[END_REF],

(2.6)

f = ρM + g, with ρ(t, x) = ⟨f ⟩ V and ⟨g⟩ V = 0.
We introduce the orthogonal projector Π in L 2 V, M -1 dµ onto N (L): Πh = ⟨h⟩ V M , which will be useful to derive the micro-macro model. Substituting (2.6) into (2.1)

and applying successively Π and (I -Π) enables to get the micro-macro model satisfied by (ρ, g)

∂ t ρ + 1 ϵ ∇ x • ⟨vg⟩ V = 0, (2.7) ∂ t g + 1 ϵ (I -Π) (v • ∇ x g) + 1 ϵ vM • ∇ x ρ = 1 ϵ 2 Lg. (2.8)
Initial conditions for macro and micro equations become (2.10) whereas the boundary conditions for ρ and g become periodic if f is periodic. From the micro part (2.8), a Chapman-Enskog expansion of g can be performed to get

ρ(0, x) = ρ init (x) = f init (x, •) V , (2.9) g(0, x, v) = g init (x, v) = f init (x, v) -ρ init (x)M (v),
g = -ϵ ϵ 2 ∂ t -L -1 (I -Π) (v • ∇ x g) + vM • ∇ x ρ = ϵL -1 (vM ) • ∇ x ρ + O(ϵ 2 ),
under some suitable smoothness assumptions. Inserting this expression in (2.7) leads to (2.5) in the limit ϵ → 0.

Time integrators.

In this part, we present the family of high order time integrators for the micro-macro model (2.7)-(2.8). We will keep the phase space variables continuous to ease the reading. We first recall the first order temporal scheme which leads to the implicit treatment of the asymptotic diffusion model before introducing the high order version.

3.1. First order accurate time integrator. Given ρ n , g n that approximate ρ, g at time t = n∆t, we obtain the solution ρ n+1 , g n+1 from the following time integration of (2.7) and (2.8) respectively. We use the following first order implicitexplicit (IMEX) strategy to attain the asymptotic preserving property

ρ n+1 = ρ n - ∆t ϵ ∇ x • vg n+1 V , (3.1)
g n+1 = g n - ∆t ϵ (I -Π) (v • ∇ x g n ) - ∆t ϵ vM • ∇ x ρ n+1 + ∆t ϵ 2 Lg n+1 . (3.2)
This manuscript is for review purposes only.

Let us observe that this scheme is different from the IMEX strategies employed in [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] Jang | High order asymptotic preserving DG-IMEX schemes for discrete-velocity kinetic equations in a diffusive scaling[END_REF], due to our implicit treatment of density gradient in micro equation (3.2) and fully implicit treatment of the macro equation. This strategy enables us to get an implicit scheme for diffusion equation in the limit ϵ → 0.

Although the macro equation is treated in a fully implicit manner, ρ n+1 and g n+1 can be updated using (3.1) and (3.2) in an explicit manner. From (3.2), we get

(3.3) g n+1 = ϵ 2 I -∆tL -1 ϵ 2 g n -ϵ∆t (I -Π) (v • ∇ x g n ) -ϵ∆tvM • ∇ x ρ n+1 .
Inserting this in (3.1), we obtain the following implicit scheme for the macro unknown

ρ n+1 = ρ n -∆t∇ x • ⟨v ϵ 2 I -∆tL -1 ϵg n -∆t (I -Π) (v • ∇ x g n )-∆tvM •∇ x ρ n+1 ⟩ V ,
or, expressing ρ n+1 as quantities at iteration n

ρ n+1 = I -∆t 2 ∇ x • (D ϵ,∆t ∇ x ) -1 ρ n -∆t∇ x • v ϵ 2 I -∆tL -1 (ϵg n -∆t (I -Π) (v•∇ x g n )) V with D ϵ,∆t = ⟨v ⊗ ϵ 2 I -∆tL -1 (vM )⟩ V .
Thanks to this time integrator, ρ n+1 can be updated by inverting a diffusion type operator. Following this, g n+1 can be found explicitly from the knowledge of ρ n+1 . This first order scheme introduced in [START_REF] Lemou | Relaxed micro-macro schemes for kinetic equations[END_REF][START_REF] Crouseilles | An asymptotic preserving scheme based on a micro-macro decomposition for collisional Vlasov equations: diffusion and high-field scaling limits[END_REF] is the basis of the high order scheme presented below.

3.2. High order accurate time integrators. Following previous works [START_REF] Dimarco | Asymptotic Preserving Implicit-Explicit Runge-Kutta Methods for Nonlinear Kinetic Equations[END_REF][START_REF] Jang | High order asymptotic preserving DG-IMEX schemes for discrete-velocity kinetic equations in a diffusive scaling[END_REF][START_REF] Boscarino | Implicit-Explicit Runge-Kutta Schemes for Hyperbolic Systems and Kinetic Equations in the Diffusion Limit[END_REF], we will consider globally stiffly accurate (GSA) IMEX Runge-Kutta (RK) schemes to construct high order time integrators for the micro-macro model (2.7) and (2.8).

An IMEX RK scheme is represented using the double Butcher tableau [START_REF] Carpenter | Additive runge-kutta schemes for convection-diffusionreaction equations[END_REF][START_REF] Ascher | Implicit-explicit runge-kutta methods for time dependent partial differential equations[END_REF] (3. An IMEX RK method is type A if the matrix A is invertible, and it is type CK if the first row of matrix A has zero entries and the square sub-matrix formed by excluding the first column and row of A is invertible. In the special case where the first column of A has zero entries, the scheme is said to be of type CK-ARS. The reader is referred to [START_REF] Dimarco | Asymptotic Preserving Implicit-Explicit Runge-Kutta Methods for Nonlinear Kinetic Equations[END_REF] for more details. In this work, we employ both type A and CK-ARS schemes.

The first order GSA IMEX RK scheme employed in (3.1) and (3. We now use the general IMEX RK scheme from (3.4) with GSA property (3.5) for obtaining high order accurate time integration of macro and micro (2.7) and (2.8) respectively. We introduce the following notations in the presentation of our scheme.

T h (k) = (I -Π) v • ∇ x h (k) , (3.7) D (j) ϵ,∆t = v ⊗ ϵ 2 I -a jj ∆tL -1 (vM ) V , (3.8)
I (j) ϵ,∆t = ϵ 2 I -a jj ∆tL -1 . (3.9)
We will construct high order IMEX RK schemes following the first order guidelines (fully implicit treatment of macro equation, implicit treatment of density gradient and relaxation terms and explicit treatment of transport term in micro equation). Given ρ n , g n that approximate ρ, g at time t = n∆t, we obtain the internal RK stage values ρ (j) and g (j) , j = 1, . . . , s as

ρ (j) = ρ n - j k=1 a jk ∆t ϵ ∇ x • vg (k) V , (3.10) g (j) = g n - j-1 k=1 ãjk ∆t ϵ T g (k) - j k=1 a jk ∆t ϵ vM • ∇ x ρ (k) + j k=1 a jk ∆t ϵ 2 Lg (k) , (3.11)
where, as usual, the summation j-1 k=1 in the explicit term is zero for j = 1.

Although the expressions above are implicit, the stage values ρ (1) , g (1) can be found in an explicit manner by using the known quantities ρ n , g n , and the stage values ρ (j) , g (j) , ∀j ∈ {2, 3, . . . , s} can be found explicitly from ρ n , g n and the previous stage values ρ (l) , g (l) , ∀l ∈ {1, 2, . . . , j -1}. Indeed, proceeding similarly as for the first order scheme, we get the following expression of g (j) , j = 1, . . . , s from (3.11), (3.12)

g (j) = I (j) ϵ,∆t ϵ 2 g n -ϵ j-1 k=1 ãjk ∆tT g (k) -ϵ j k=1 a jk ∆tvM • ∇ x ρ (k) + j-1 k=1
a jk ∆tLg (k) .

Further, we write (3.10) by splitting the summation on k as

ρ (j) = ρ n - j-1 k=1 a jk ∆t ϵ ∇ x • vg (k) V -a jj ∆t ϵ ∇ x • vg (j) V ,
and inserting (3.12) in the last term leads to the update of ρ (j) for j = 1, . . . , s

ρ (j) = I -a 2 jj ∆t 2 ∇ x • D (j) ϵ,∆t ∇ x -1 ρ n - j-1 k=1 a jk ∆t ϵ ∇ x • vg (k) V (3.13) -a jj ∆t∇ x • vI (j) ϵ,∆t ϵg n - j-1 k=1 ãjk ∆tT g (k) - j-1 k=1 a jk ∆tvM • ∇ x ρ (k) + 1 ϵ j-1 k=1 a jk ∆tLg (k) V ,
where the definition of T , D

ϵ,∆t and I (j) ϵ,∆t are given by (3.7)-(3.9). After this reformulation, ρ (j) can be computed from (3.13) by inverting a linear elliptic type problem

This manuscript is for review purposes only.

and following this, g (j) can be found from (3.12). The GSA property in (3.5) guarantees that the solution at t n+1 = (n + 1)∆t is same as the last RK stage values, that is, ρ n+1 = ρ (s) and g n+1 = g (s) .

4. Asymptotic preserving property. In this section, we show that the time integrated scheme (3.13)-(3.12) becomes a consistent scheme for the diffusion equation (2.5) in the limit ϵ → 0. We will discuss the asymptotic preserving property for both CK-ARS and type A time integrators as performed in [START_REF] Dimarco | Asymptotic Preserving Implicit-Explicit Runge-Kutta Methods for Nonlinear Kinetic Equations[END_REF] for the fluid limit. First, we recall the definition of well-prepared initial data in our context.

Definition 4.1 (Well-prepared initial data). The initial data ρ(0, x) and g(0, x, v) in (2.9) and (2.10) are said to be well-prepared if g(0, x, v) = O(ϵ).

Lemma 4.2. Assume that ϵ is sufficiently small. Let ãjk and a jk be the coefficients of the RK method (3.4) applied to the scheme (3.10)- (3.11). Then, the following holds:

1. CK-ARS case: If g n = O(ϵ), then g (1) = g n = O(ϵ) and
g (j) = ϵL -1 (vM ) • ∇ x ρ (j) + O ϵ 2 , ∀j ∈ {2, . . . , s}.
2. Type A case:

g (j) = ϵL -1 (vM ) • ∇ x ρ (j) + O ϵ 2 , ∀j ∈ {1, . . . , s}.
Proof. Let j ∈ {1, . . . , s} such that a jj ̸ = 0. Observe that the operator

I (j) ϵ,∆t
defined in (3.9) admits, for small ϵ, the following expansion:

(4.1)

I (j) ϵ,∆t = -(a jj ∆tL) -1 + O(ϵ 2 ).
Consider now an A-type time integrator, so with a jj ̸ = 0 for any j ∈ {1, . . . , s}, and assume g n = O(1). From (3.12) and the previous expansion, we obtain

g (1) = -(a 11 ∆tL) -1 -ϵa 11 ∆tvM • ∇ x ρ (1) + O(ϵ 2 ) = ϵL -1 (vM ) • ∇ x ρ (1) + O(ϵ 2 ).
Now, the proof is performed by induction on j ∈ {2, . . . , s} assuming that for any

k ∈ {1, . . . , j -1}, g (k) = ϵL -1 (vM ) • ∇ x ρ (k) + O(ϵ 2
). In particular g (k) = O(ϵ) and the formula (3.12) has therefore the following expansion:

g (j) = -(a jj ∆tL) -1 O(ϵ 2 ) -ϵ j k=1 a jk ∆tvM • ∇ x ρ (k) + j-1 k=1 a jk ∆tLg (k) + O(ϵ 2 ).
Inserting the induction hypothesis in the last sum, most of the terms in the two sums eliminate so that finally g

(j) = ϵL -1 (vM ) • ∇ x ρ (j) + O(ϵ 2 ).
The case of a CK-ARS time integrator is slightly different. First a 11 = 0 so that g (1) = g n = O(ϵ) by the particular well-prepared assumption. Now a 22 ̸ = 0 and (3.12)

has the following expansion for j = 2:

g (2) = -(a 22 ∆tL) -1 O(ϵ 2 )-ϵa 22 ∆tvM • ∇ x ρ (2) +O(ϵ 2 ) = ϵL -1 (vM )•∇ x ρ (2) + O(ϵ 2 ).
Again, the proof is by induction on j ∈ {3, . . . , s} assuming for any k ∈ {2, . . . , j -1},

g (k) = ϵL -1 (vM ) • ∇ x ρ (k) + O(ϵ 2
). The same computation as above is available since g (1) = O(ϵ). One has (note that a j1 = 0 for any j so that the sums start at k = 2):

g (j) = -(a jj ∆tL) -1 O(ϵ 2 ) -ϵ j k=2 a jk ∆tvM • ∇ x ρ (k) + j-1 k=2 a jk ∆tLg (k) + O(ϵ 2 ) = ϵL -1 (vM ) • ∇ x ρ (j) + O(ϵ 2 ).
This manuscript is for review purposes only.

Due to the GSA property of both time integrators, we have

g n+1 = g (s) = ϵL -1 (vM ) • ∇ x ρ (s) + O ϵ 2 = ϵL -1 (vM ) • ∇ x ρ n+1 + O ϵ 2
for sufficiently small ϵ. Thus, the following are evident from Lemma 4.2:

1. For type CK-ARS, if the initial data is well-prepared (that is,

g 0 = O(ϵ)), then g n = O(ϵ), ∀n > 0.
2. For type A, if the initial data is such that

g 0 = O(1), then g n = O(ϵ), ∀n > 0.
As observed in [START_REF] Dimarco | Asymptotic Preserving Implicit-Explicit Runge-Kutta Methods for Nonlinear Kinetic Equations[END_REF], the initial data does not need to be well-prepared for type A time integrator, unlike type CK-ARS, to ensure AP property. degenerates to the following scheme for the diffusion equation

(4.2) ρ (j) = ρ n + j k=1 a jk ∆t∇ x • κ∇ x ρ (k) , ∀j = 1, . . . , s, κ = -v ⊗ L -1 (vM ) V .
Proof. Corresponding to each case (CK-ARS or type A), we have the following:

Type CK-ARS Assumptions in criterion 1 of Lemma 4.2 are satisfied, and its implications can be utilised. Hence, inserting

g (ℓ) = ϵL -1 (vM ) • ∇ x ρ (ℓ) + O(ϵ 2 ), ∀ℓ ∈ {2, 3, .
., s} into (3.10), we get (recall that a j1 = 0)

ρ (j) = ρ n - ∆t ϵ j k=2 a jk ∇ x • vϵL -1 (vM ) • ∇ x ρ (k) V + O(ϵ), = ρ n -∆t j k=2 a jk ∇ x • v ⊗ L -1 (vM ) V ∇ x ρ (k) + O(ϵ).
Type A Assumptions in criterion 2 of Lemma 4.2 are satisfied, and its implications can be utilised. Hence, inserting

g (ℓ) = ϵL -1 (vM ) • ∇ x ρ (ℓ) + O(ϵ 2 ), ∀ℓ ∈ {1, 2, .
., s} into (3.10), we get the required result by following the same simplification as before. The only difference is that here j k=1 instead of j k=2 .

Remark 4.4. For type CK-ARS, if the initial data is not well-prepared, computing

g (2) from (3.11) involves ϵ ã21 a22 L -1 (I -Π)(v • ∇ x g (1)
) which is not of O(ϵ 2 ). Thus,

g (2) = ϵ ã21 a 22 L -1 (I -Π)(v • ∇ x g (1) ) + ϵL -1 (vM ) • ∇ x ρ (2) + O ϵ 2 ,
and inserting in the macro equation (3.10) for j = 2 leads to (since a 21 = 0)

ρ (2) = ρ n - ã21 a 22 ∆t v ⊗ L -1 (I -Π)v∇ 2 x g (1) V -a 22 ∆t∇ x • v ⊗ L -1 (vM ) V ∇ x ρ (2) + O(ϵ),
which is not consistent with the diffusion equation. Thus, for CK-ARS, asymptotic consistency cannot be attained if the initial data is not well-prepared.

Space and velocity discretization.

In this section, we present the spatial (for both non-staggered and staggered grids) and velocity discretization strategies that we employ in our numerical scheme.

This manuscript is for review purposes only.

5.1. Discrete velocity method. For the velocity discretization, we will follow the discrete velocity method [START_REF] Jin | The discrete-ordinate method in diffusive regimes[END_REF]. Thus, the velocity domain is truncated as v ∈

[-v max , v max ], and a uniform mesh is used

v k = -v max +k∆v, k = 1, . . . , N v (N v ∈ N ⋆ )
and ∆v = 2v max /N v . Further, f (t, x, v) and M (v) are represented as:

f k (t, x) := f (t, x, v k ), M k := M (v k ) for k = 1, . . . , N v .
Then, according to the definitions (2.3) and (2.6), we have for j

= 1, . . . , N v ρ(t, x) ≈ Nv-1 k=0 f k ∆v Nv-1 k=0 M k ∆v and (Πf (t, x, v)) j ≈ Nv-1 k=0 f k ∆v Nv-1 k=0 M k ∆v M j .
For the presentation, we will skip the velocity part to focus on space discretization.

5.2. Space discretization using staggered grid. First, we will consider staggered grid to approximate g (j) and ρ (j) in space following [START_REF] Lemou | A New Asymptotic Preserving Scheme Based on Micro-Macro Formulation for Linear Kinetic Equations in the Diffusion Limit[END_REF]: the two meshes of the space interval [0, 1] are

x i = i∆x and x i+1/2 = (i+1/2)∆x for i = 0, . . . , N x (N x ∈ N ⋆ ),
with ∆x = L/N x . Periodic boundary conditions will be considered in this section.

The expressions for g (j) and ρ (j) in (3.12)-(3.13) are spatially discretised by considering staggered grid: ρ (j) is stored at x i (ρ

(j) i ≈ ρ (j) (x i ))
, and g (j) is stored at

x i+1/2 (g (j) i+1/2 (v) ≈ g (j) (x i+1/2 , v)). The term v • ∇ x g (k) in (3.12) and (3.13) is discre- tised in an upwind fashion as v • ∇ x ≈ v + • G - upw + v -• G + upw where v ± = (v ± |v|)/2, G ± upw denote the N x × N x matrices that approximate ∇ x .
For instance, the first order version is

(5.1) G - upw = 1 ∆x circ([-1, 1]), G + upw = 1 ∆x circ([-1, 1]),
where the notation circ is defined in Appendix A. With these notations, we get v∂ x g (j)

x i+1/2 ≈ v + g (j) i+ 1 2 -g (j) i-1 2 ∆x + v - g (j) i+ 3 2 -g (j) i+ 1 2 ∆x = v + G - upw + v -G + upw g (j) i ,
where in the last term, the i index has to be understood as the i-th component of the vector. Instead of first order upwind discretization, one can also use high order upwind discretizations so that the matrices G ± upw will be different. Further, the term vM • ∇ x ρ (k) in (3.12)-(3.13) and the terms of the form

∇ x • ⟨(•)⟩ V in (3.13) are
discretised using second order central differences as in [START_REF] Lemou | A New Asymptotic Preserving Scheme Based on Micro-Macro Formulation for Linear Kinetic Equations in the Diffusion Limit[END_REF]. In particular, the term

vM • ∇ x ρ (k) is approximated by (5.2) vM ∂ x ρ (k) x i+1/2 ≈ vM ρ (k) i+1 -ρ (k) i ∆x = vM G ceng ρ (k) i , G ceng = 1 ∆x circ([-1, 1]).
Finally, the gradient terms ∇ x • ⟨(•)⟩ V in (3.13) are approximated as follows

(5.3) (∂ x ⟨•⟩ V ) xi = (⟨•⟩ V ) i+1/2 -(⟨•⟩ V ) i-1/2 ∆x = G cenρ ⟨•⟩ V i , G cenρ = 1 ∆x circ([-1, 1]).
Again, high order centered finite differences methods can be used so that it will give different expressions for G cenρ and G ceng . Let us remark that the term

∇ x • ∇ x = ∇ 2 x in (3.13) is approximated by G cenρ G ceng , ie G cenρ G ceng = 1 ∆x 2 circ([1, -2, 1]), which
gives the standard second order approximation of the Laplacian.

This manuscript is for review purposes only.

To ease the reading, we present the fully discrete scheme for first order ARS(1, 1, 1)

but the generalization to high order can be done using the elements of Section 3

g n+1 = ϵ 2 I -∆tL -1 ϵ 2 g n -ϵ∆t (I -Π) v + G - upw +v -G + upw g n -ϵ∆tvM G ceng ρ n+1 ρ n+1 = I -∆t 2 G cenρ v ⊗ ϵ 2 I -∆tL -1 (vM ) V G ceng -1 × ρ n -∆tG cenρ v ϵ 2 I -∆tL -1 ϵg n -∆t (I -Π) v + G - upw +v -G + upw g n V .
5.3. Space discretization using non-staggered grid. We also address the case of non-staggered grids which may be more appropriate when high dimensions are considered in space since only one spatial mesh is used:

x i = i∆x, for i = 0, 1, .., N x ,
with ∆x = L/N x . Let g (j) and ρ (j) in (3.12)-(3.13) ∀j ∈ {1, 2, .., s} be approximated in space by g

(j) i (v) ≈ g (j) (x i , v) and ρ (j) i ≈ ρ (j) (x i ). The term v•∇ x g (k) in (3.12)-(3.13) is discretised in an upwind fashion as v•∇ x = v + G - upw +v -G + upw , where v ± = (v ± |v|)/2.
Here, G ± upw denote the matrices that represent an upwind approximation of ∇ x . For instance, the definition (5.1) can be used, but also its third order version (5.4)

G - upw = 1 6∆x circ([1, -6, 3, 2]), G + upw = 1 6∆x circ([-2, -3, 6, -1]),
where circ represents the matrix notation described in Appendix A can be used. The term vM • ∇ x ρ (k) in (3.12)-(3.13) and the terms of the form

∇ x • ⟨(•)⟩ V in (3.13) are
discretised in central fashion, since these terms act as source in (3.12) and diffusion in (3.13). Here, ∇ x is approximated by central differences as in (5.3) or (5.2) but in the non-staggered case, the same matrix can be used for both terms. As an example, the fourth order central difference produces:

(5.5)

G cen = 1 12∆x circ([1, -8, 0, 8, -1]).
The term

∇ x • ∇ x = ∇ 2 x in (3.13) is discretised as the matrices product G 2 cen =
G cen G cen . Like in the staggered grid case, we present the fully discrete scheme for first order ARS(1, 1, 1) time discretization to ease the reading:

g n+1 = ϵ 2 I -∆tL -1 ϵ 2 g n -ϵ∆t (I -Π) v + G - upw + v -G + upw g n -ϵ∆tvM G cen ρ n+1 ρ n+1 = I -∆t 2 G cen v ⊗ ϵ 2 I -∆tL -1 (vM ) V G cen -1 × ρ n -∆tG cen v ϵ 2 I -∆tL -1 ϵg n -∆t (I -Π) v + G - upw + v -G + upw g n V Remark 5.1.
We know that the term

j k=1 a jk ∆t ϵ ∇ x • vg (k) V in (3.10
) is split into first j -1 and last j contributions, and g (j) is substituted for the last j contribution, as in (3.13). The gradient in

j-1 k=1 a jk ∆t ϵ ∇ x • vg (k) V of (3.13) is discretised using G cenρ .
Further, the substitution of g (j) for the last j hints the combination of

∇ x • ∇ x as ∇ 2
x for the terms of g (j) involving ∇ x g and ∇ x ρ. However, if we choose a spatial discretization for ∇ 2 x as G diff , then these terms will experience G cenρ G ceng for the first j -1 contributions and G diff for the last j contribution of the ρ (j) update equation. This disrupts the ODE structure present in RK time discretization, and hence reduction to first order time accuracy was observed numerically. Therefore, in order to retain high order time accuracy, it is important to carry out the space discretization carefully. Hence, we do not introduce a different discretization for ∇ 2

x , and we retain G cenρ G ceng even for the last j contribution of ρ (j) equation.

This manuscript is for review purposes only.

Remark 5.2. The matrices introduced for spatial discretization do not change the Chapman-Enskog expansion so that the AP property is still true in the fully discrete form. Thus, we have g (k) = ϵL -1 (vM )G ceng ρ (k) + O(ϵ 2 ) for k ∈ {1, . . . , s} by using type A. For CK-ARS with well-prepared data, we have

g (k) = ϵL -1 (vM )G ceng ρ (k) + O(ϵ 2
) for k ∈ {2, . . . , s}. Inserting this in macro equation, we get the corresponding RK scheme for the diffusion

ρ (j) = ρ n -∆t j k=1 a jk G cenρ v ⊗ L -1 (vM ) V G ceng ρ (k) + O(ϵ).
6. Extensions to other collision operator and inflow boundary problems. In this section, we show that our high order AP schemes can be extended to other problems involving advection-diffusion asymptotics and inflow boundaries.

6.1. Advection-diffusion asymptotics. In this part, an advection-diffusion collision operator is considered (see [START_REF] Jin | Uniformly accurate diffusive relaxation schemes for multiscale transport equations[END_REF][START_REF] Jang | High order asymptotic preserving DG-IMEX schemes for discrete-velocity kinetic equations in a diffusive scaling[END_REF]), (6.1)

Lf := Lf + ϵvM • A ⟨f ⟩ V , A ∈ R d , |ϵA| < 1,
where L denotes a collision satisfying the properties listed in Section 2. A famous

simple example is Lf = ⟨f ⟩ V M -f .
Using the notations introduced in Section 2, we can derive the micro-macro model satisfied by ρ = ⟨f ⟩ V and g = f -ρM by applying Π and I -Π to (2.1) with collision L to get the macro and micro equations in this context

∂ t ρ + 1 ϵ ∇ x • ⟨vg⟩ V = 0, (6.2)
∂ t g + 1 ϵ (I -Π) (v • ∇ x g) + 1 ϵ vM • ∇ x ρ = 1 ϵ 2 Lg + 1 ϵ vM • Aρ. (6.3)
A Chapman-Enskog expansion can be performed to get g = ϵL -1 (vM ) • ∇ x ρ -ϵL -1 (vM ) • Aρ + O(ϵ 2 ). Inserting this in the macro equation (6.2) enables to obtain an advection-diffusion equation in the limit ϵ → 0: (6.4)

∂ t ρ + ∇ x • v ⊗ L -1 (vM ) V ∇ x ρ -∇ x • v ⊗ L -1 (vM ) V Aρ = 0.
The goal is to design a uniformly stable high order time integrators for (6.2)-(6.3) so that they degenerate into a high order time integrator for (6.4) as ϵ → 0. The extension of the schemes introduced in Section 3 will lead to an IMEX discretization of the asymptotic model (6.4), where the advection term is treated explicitely while the diffusion term is implicit.

6.1.1. High order time integrator. In this subsection, we present the discretization of macro and micro equations (6.2)-(6.3). As in Section 3, in the micro equation, we treat 1 ϵ 2 Lg implicitly to ensure uniform stability and the additional term 1 ϵ vM • Aρ explicitly since it will be stabilized by the implicit treatment of the stiffest term. Regarding the macro equation and the remaining terms in micro equation, we follow the lines from previous Section 3. We thus obtain the following high order IMEX RK scheme to approximate (6.2)-(6.3)

ρ (j) = ρ n - j k=1 a jk ∆t ϵ ∇ x • vg (k)
V , (6.5)

g (j) = g n - ∆t ϵ j-1 k=1 ãjk Tg (k) + j k=1 a jk vM•∇ x ρ (k) - j k=1 a jk ϵ Lg (k) - j-1 k=1 ãjk vM•Aρ (k) , (6.6)
This manuscript is for review purposes only.

where the coefficients a jk , ãjk are given by the Butcher tableaux. As in Section 3, some calculations are required to make the algorithm explicit. First, we have (6.7)

g (j) = I (j) ϵ,∆t ϵ 2 g n -ϵ∆t j-1 k=1 ãjk T g (k) + j k=1 a jk vM • ∇ x ρ (k) - 1 ϵ j-1 k=1 a jk Lg (k) - j-1 k=1 ãjk vM • Aρ (k) ,
with

Tg (k) = (I -Π) v•∇ x g (k) and I (j)
ϵ,∆t = ϵ 2 I -a jj ∆tL -1 . Then, ρ (j) is obtained by inserting g (j) given by (6.7) in the macro equation (6.5) to get

ρ (j) = I -a 2 jj ∆t 2 ∇ x • D (j) ϵ,∆t ∇ x -1 ρ n - j-1 k=1 a jk ∆t ϵ ∇ x • vg (k) V (6.8) -a jj ∆t∇ x • vI (j) ϵ,∆t ϵg n - j-1 k=1 ãjk ∆tT g (k) - j-1 k=1 a jk ∆tvM • ∇ x ρ (k) + 1 ϵ j-1 k=1 a jk ∆tLg (k) + j-1 k=1 ãjk ∆tvM • Aρ (k) V ,
where D

ϵ,∆t = ⟨v ⊗ ϵ 2 I -a jj ∆tL -1 (vM )⟩ V . Thus, ρ (j) can be updated by using (6.8) and g (j) can be found explicitly by using (6.7).

6.1.2. Asymptotic preserving property. This part is dedicated to the asymptotic preserving property of the scheme (6.8)-(6.7). We first show the AP property of type A time integrator, and we later remark how this property is true for the CK-ARS time integrator with well-prepared initial data. First we have Lemma 6.1. If g n = O(1) and g (k) = O(ϵ), ∀k ∈ {1, 2, . . . , j -1}, then g (j) = O(ϵ), ∀j ∈ {2, 3, .., s} for small ϵ. In particular, we have ∀j ∈ {2, 3, .., s} (6.9)

g (j) = ϵ j k=1 a jk a jj L -1 (vM ) • ∇ x ρ (k) - j-1 k=1 a jk a jj g (k) -ϵ j-1 k=1 ãjk a jj L -1 (vM ) • Aρ (k) + O(ϵ 2 ).
Proof. Plugging in (6.7) the expansion (4.1) of I (j) ϵ,∆t given by (3.9), along with the assumptions stated in the Lemma, we obtain (6.9) from which we deduce g (j) = O(ϵ) for all j ∈ {2, 3, .., s}.

Remark 6.2. For type A time integrator, if g n = O(1), we have from (6.7):

g (1) = ϵ a 11 a 11 vM • ∇ x ρ (1) + O(ϵ 2 ) = O(ϵ).
This satisfies the induction hypothesis in Lemma 6.1. Further, (6.9) holds by omitting j-1 k=1 terms for j = 1. Thus, (6.9) is true for j ∈ {1, 2, .., s}.

Lemma 6.1 enables to get an expansion of g (j) that can be inserted in (6.8) to identify the time discretization of the asymptotic limit. However, this leads to quite involved calculations which requires to introduce some notations.

This manuscript is for review purposes only.

Definition 6.3. For j ∈ {1, 2, .., s} and k 1 , m ∈ {1, 2, .., j} we define

(6.10) Π m j,k1 = v a jk1 a k1k1 S k0 S k1 S k2 . . . S km-1 R km V , with S k0 = 1, S k l = k l -1 k l+1 =1 a k l k l+1 a k l+1 k l+1 for l ∈ {1, 2, .., m -1}, m ≥ 2, R km = km km+1=1 a kmkm+1 L -1 (vM ) • ∇ x ρ (km+1) - km-1 km+1=1 ãkmkm+1 L -1 (vM ) • Aρ (km+1) .
As usual, we will use the convention

q j=1 ≡ 0 if q ∈ Z\N.
The term Π m j,k1 will be useful in the following study and deserves some remarks:

the index m denotes the depth of the embedded sums, j corresponds to the current stage and k 1 corresponds to the indexing over previous stages. We continue with the following lemma which gives an induction relation on Π m j,k1 .

Lemma 6.4. For j ≥ 2, we have

Π m j,j = j-1 k1=1
Π m-1 j,k1 for m ∈ {2, 3, .., j}, and Π j j,k1 = 0 for k 1 ∈ {1, 2, .., j -1}.

Proof. For the first relation, considering k 1 = j (with j ≥ 2) in (6.10) leads to

Π m j,j = v S k0 S j S k2 . . . S km-1 R km V ,
since a jj ̸ = 0. Further, since S k1=j = j-1 k2=1 a jk2 a k2k2 , we get

Π m j,j = v j-1 k2=1 a jk2 a k2k2 S k0 S k2 ...S km-1 R km V
By employing the change of variables as k ℓ → k ℓ-1 for ℓ ∈ {2, 3, .., m} in the right hand side of above expression, we get

Π m j,j = v j-1 k1=1 a jk1 a k1k1 S k0 S k1 . . . S km-2 R km-1 V = j-1 k1=1 v a jk1 a k1k1 S k0 S k1 . . . S km-2 R km-1 V = j-1 k1=1 Π m-1 j,k1 ,
which proves the first identity.

For the second relation, considering m = j in (6.10) leads to

Π j j,k1 = v a jk1 a k1k1 S k0 S k1 S k2 ...S kj-1 R kj V
We first prove the relation for j = 2. It is clear from Definition 6.3 that the summation in S k1 goes from

k 2 = 1 to k 2 = k 1 -1. For k 1 = 1,

the summation goes to

This manuscript is for review purposes only.

HIGH ORDER AP SCHEME FOR KINETIC EQUATIONS IN DIFFUSIVE SCALING 13

k 2 = k 1 -1 = 0. Thus, since S k1 involves 0 1 for k 1 = 1, it is zero according to the convention. Hence Π j j,k1 = 0 for k 1 = 1.

We now prove the relation for j > 2. From Definition 6.3, it can be seen that the summations in S k1 and S k2 go from

k 2 = 1 to k 2 = k 1 -1 and k 3 = 1 to k 3 = k 2 -1
respectively. Thus, the summation in S k2 can go to atmost

k 3 = k 2 -1 = (k 1 -1)-1 = k 1 -2.
Proceeding in this manner, we see that the summation in S kj-1 can go to atmost k j = k 1 -(j -1).

For k 1 ∈ {1, 2, .., j -1}, k j = k 1 -(j -1) ∈ Z\N so that S kj-1 = 0 and hence Π j j,k1 = 0 for k 1 ∈ {1, 2, .., j -1} which ends the proof. Now, we can use the previous Lemma to identify the asymptotic numerical scheme.

Lemma 6.5. When ϵ → 0, the numerical scheme (6.5)-(6.6) degenerates into (6.11)

ρ (j) = ρ n + ∆t j k1=1 ∇ x • j ℓ=1 (-1) ℓ Π ℓ j,k1
for j ∈ {1, 2, .., s}, where Π ℓ j,k1 is given by Definition 6.3.

Proof. We start with the macro equation in (6.5)

ρ (j) = ρ n - j k1=1 a jk1 ∆t ϵ ∇ x • ⟨vg (k1) ⟩ V ,
in which we insert g (k1) given by (6.9) to get

ρ (j) = ρ n -∆t j k1=1 ∇ x • v a jk1 a k1k1 k1 k2=1 a k1k2 L -1 (vM ) • ∇ x ρ (k2) - k1-1 k2=1 ãk1k2 L -1 (vM ) • Aρ (k2) V + ∆t ϵ j k1=1 ∇ x • v a jk1 a k1k1 k1-1 k2=1 a k1k2 g (k2) V + O(ϵ) = ρ n -∆t j k1=1 ∇ x • v a jk1 a k1k1 S k0 R k1 V + ∆t ϵ j k1=1 ∇ x • v a jk1 a k1k1 k1-1 k2=1 a k1k2 g (k2) V + O(ϵ) = ρ n -∆t j k1=1 ∇ x • Π 1 j,k1 + ∆t ϵ j k1=1 ∇ x • v a jk1 a k1k1 k1-1 k2=1 a k1k2 g (k2) V + O(ϵ).
Inserting g (k2) from (6.9) in the above equation and simplifying as before, we get,

ρ (j) = ρ n -∆t j k1=1 ∇ x • Π 1 j,k1 -Π 2 j,k1 - ∆t ϵ j k1=1 ∇ x • v a jk1 a k1k1 k1-1 k2=1 a k1k2 a k2k2 k2-1 k3=1 a k2k3 g (k3) V +O(ϵ).
This manuscript is for review purposes only.

This procedure can be continued (j -1) times to finally get,

ρ (j) = ρ n + ∆t j k1=1 ∇ x • j-1 ℓ=1 (-1) ℓ Π ℓ j,k1 -(-1) j-1 ∆t ϵ j k1=1 ∇ x • v a jk1 a k1k1   k1-1 k2=1 a k1k2 a k2k2 • • • kj-2-1 kj-1=1 a kj-2kj-1 a kj-1kj-1 kj-1-1 kj =1 a kj-1kj g (kj )   V + O(ϵ) = ρ n + ∆t j k1=1 ∇ x • j-1 ℓ=1 (-1) ℓ Π ℓ j,k1 -(-1) j-1 ∆t ϵ j k1=1 ∇ x • v a jk1 a k1k1   S k0 S k1 . . . S kj-2 kj-1-1 kj =1 a kj-1kj g (kj )   V + O(ϵ).
We know from Definition 6.3 that the summations in S k1 and S k2 go from k 2 = 1 to

k 2 = k 1 -1 and k 3 = 1 to k 3 = k 2 -1
respectively. Thus, the summation in S k2 can go

to atmost k 3 = k 2 -1 = (k 1 -1) -1 = k 1 -2.
Proceeding in this manner, we see that the summations in S kj-2 and

kj-1-1 kj =1 a kj-1kj g (kj ) go to atmost k j-1 = k 1 -(j -2)
and

k j = k 1 -(j -1) respectively.
Since the summation in k 1 goes to atmost j in the above equation, k j in the term

kj-1-1 kj =1
a kj-1kj g (kj ) goes to atmost

k j = k 1 -(j -1) = j -(j -1) = 1, and k j-1 in S kj-2 goes to atmost k j-1 = k 1 -(j -2) = j -(j -2)
= 2 and so on. Thus, only k j = 1 remains in the last summation so that

kj-1-1 kj =1
a kj-1kj g (kj ) = a 21 g (1) =

ϵa 21 L -1 (vM )•∇ x ρ (1) +O(ϵ 2 ) = a21 a11 ϵa 11 L -1 (vM )•∇ x ρ (1) +O(ϵ 2 ) = ϵS kj-1 R kj +O(ϵ 2 ).
Thus, we have

ρ (j) = ρ n + ∆t j k1=1 ∇ x • j-1 ℓ=1 (-1) ℓ Π ℓ j,k1 -(-1) j-1 ∆t j k1=1 ∇ x • v a jk1 a k1k1 S k0 S k1 . . . S kj-1 R kj V + O(ϵ) = ρ n + ∆t j k1=1 ∇ x • j-1 ℓ=1 (-1) ℓ Π ℓ j,k1 + ∇ x • (-1) j Π j j,k1
+ O(ϵ).

We can now prove the asymptotic property of the scheme (6.5)-(6.6).

Theorem 6.6. When ϵ → 0, the scheme (6.5)-(6.6) degenerates into (6.12)

ρ (j) = ρ n -∆t j k=1 a jk ∇ x • v ⊗ L -1 (vM ) V ∇ x ρ (k) + ∆t j-1 k=1 ãjk ∇ x • v ⊗ L -1 (vM ) V Aρ (k)
, for j ∈ {1, 2, . . . , s}.

Proof. From Lemma 6.5, the asymptotic limit ϵ → 0 of the macro equation in This manuscript is for review purposes only.

(6.5) is (for j ∈ {1, 2, .., s})

ρ (j) = ρ n + ∆t j k1=1 ∇ x • j ℓ=1 (-1) ℓ Π ℓ j,k1 = ρ n + ∆t∇ x • j ℓ=1 (-1) ℓ Π ℓ j,j + j-1 k1=1 Π ℓ j,k1 = ρ n + ∆t∇ x • -Π 1 j,j + j ℓ=2 (-1) ℓ Π ℓ j,j + j ℓ=1 (-1) ℓ j-1 k1=1 Π ℓ j,k1
.

Using the recurrence relation given by Lemma 6.4 and a change of indices lead to

ρ (j) = ρ n + ∆t∇ x • -Π 1 j,j + j ℓ=2 (-1) ℓ j-1 k1=1 Π ℓ-1 j,k1 + j ℓ=1 (-1) ℓ j-1 k1=1 Π ℓ j,k1 = ρ n + ∆t∇ x • -Π 1 j,j - j-1 ℓ=1 (-1) ℓ j-1 k1=1 Π ℓ j,k1 + j ℓ=1 (-1) ℓ j-1 k1=1 Π ℓ j,k1 = ρ n + ∆t∇ x • -Π 1 j,j + (-1) j j-1 k1=1 Π j j,k1
.

From Lemma 6.4, we have j-1 k1=1 Π j j,k1 = 0, so that from Definition 6.3 we get

ρ (j) = ρ n + ∆t∇ x • -Π 1 j,j = ρ n -∆t∇ x • v a jj a jj S k0 R k1=j V = ρ n -∆t∇ x •   v k1 k2=1 a k1k2 L -1 (vM ) • ∇ x ρ (k2) - k1-1 k2=1 ãk1k2 L -1 (vM ) • Aρ (k2) V   k1=j = ρ n -∆t j k2=1 a jk2 ∇ x • v ⊗ L -1 (vM ) V ∇ x ρ (k2) + ∆t j-1 k2=1 ãjk2 ∇ x • v ⊗ L -1 (vM ) V Aρ (k2) ,
which ends the proof.

Remark 6.7. For CK-ARS schemes with well-prepared initial data, we obtain g (1) = g n = O(ϵ) and ρ (1) = ρ n . The presentation in this section will apply for CK-ARS from the second RK stage onwards. For instance, Definition 6.3 applies for CK-ARS with the following change in indexes: j ∈ {2, 3, .., s}, k 1 , m ∈ {2, 3, .., j} and all the summations involved start from 2 instead of 1 since a 11 = 0. The lemmas and theorems that follow also undergo the corresponding change in indexes, and the AP property for CK-ARS can be observed for j ∈ {2, 3, .., s}.

Remark 6.8. Upon incorporating the spatial matrices corresponding to staggered grid in place of the continuous gradient operator, we obtain in the limit ϵ → 0, (6.13)

ρ (j) = I + a jj ∆tG cenρ v ⊗ L -1 (vM ) V G ceng -1 × ρ n - j-1 k=1 a jk ∆tG cenρ v ⊗ L -1 (vM ) V G ceng ρ (k) + j-1 k=1 ãjk ∆tG cenρ v ⊗ L -1 (vM ) V G avgg Aρ (k) .
The matrices G cenρ , G ceng are given in subsection 5.2 and G avg g = 1 2 circ([1, 1]). Thus,

A ρ (k) x i+1/2 = 1 2 A(ρ (k) i+1 + ρ (k) i ) = (G avg g Aρ (k) ) i .
This results in a central discretization of the advection term in the macro equation. Thus, we obtain a consistent internal RK stage approximation of the advection-diffusion equation in the limit ϵ → 0.

6.2. Inflow Boundaries. So far, periodic boundary conditions were considered.

In this part, we consider inflow boundary conditions for f solution to (2.1)

(6.14) f (t, x, v) = f b (t, x, v), (x, v) ∈ ∂Ω × V such that v • n(x) < 0, ∀t,
where f b is a given function and n(x) denotes the unitary outgoing normal vector to ∂Ω. As mentioned in [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], such boundary conditions cannot be adapted naturally to the standard micro-macro unknown ρ(t, x) and g(t, x, v) solution to (2.6). To overcome this drawback, another micro-macro decomposition is introduced in [START_REF] Lemou | Micro-Macro Schemes for Kinetic Equations Including Boundary Layers[END_REF] (6.15)

f = ρM + g, ρ(t, x) = ⟨f (t, x, •)⟩ V-, ⟨g(t, x, •)⟩ V-= 0, ⟨f ⟩ V-= V-f dµ V-M dµ ,
where the velocity domain V -is defined by

(6.16) V -(x) = {v ∈ V, ω(x, v) < 0}, V + (x) = V \V -(x).
The function ω(x, v) extends v • n(x) in the interior of domain. Some examples of ω(x, v) for different geometries are provided in [START_REF] Lemou | Micro-Macro Schemes for Kinetic Equations Including Boundary Layers[END_REF]. It can be seen that the boundary conditions for ρ(t, x) and g(t, x, v) can be evaluated from the inflow boundary condition in (6.14). Indeed, for (x, v) ∈ ∂Ω × V such that v • n(x) < 0, ∀t, we define (6.17)

ρ b (t, x) = ⟨f b (t, x, •)⟩ V-, g b (t, x, v) = f b (t, x, v) -ρ b (t, x)M (v).
The derivation of the micro-macro model needs to be adapted to this decomposition. The projector Π -is defined as Π -h = ⟨h⟩ V-M . Then, substituting (6.15) into (2.1) and applying Π -and I -Π -enable to get the macro and micro equations: [START_REF] Klar | Asymptotic-induced scheme for nonstationary transport equations in the diffusive limit[END_REF] where L = (I -Π -) L. Moreover, it can be seen that L = (I -

∂ t ρ + 1 ϵ ⟨vM ⟩ V-• ∇ x ρ + 1 ϵ ∇ x • ⟨vg⟩ V-= 1 ϵ 2 ⟨Lg⟩ V-, (6.18)
∂ t g + 1 ϵ I -Π -(v • ∇ x g) + 1 ϵ I -Π -vM • ∇ x ρ = 1 ϵ 2 Lg, (6.
Π -) L (I -Π -) = (I -Π -) L (I -Π) since Π -h, Πh ∈ N (L), ∀h.
The macro equation (6.18) turns out to be more complicated than the one obtained for standard micro-macro decomposition. It can be made simpler by using ρ = ρ + ⟨g⟩ V , f = ρM -⟨g⟩ V M + g, obtained from the decompositions (2.6) and (6.15). Applying Π to (2.1) instead of Π -, we obtain the simpler macro equation, (6.20)

∂ t ρ + 1 ϵ ∇ x • ⟨vg⟩ V = 0,
and the micro-macro system that we will consider in the sequel is (6.19)- (6.20).

This manuscript is for review purposes only.

6.2.1. Numerical scheme. In this part, we present the fully discretized scheme to approximate (6. [START_REF] Klar | Asymptotic-induced scheme for nonstationary transport equations in the diffusive limit[END_REF])- (6.20). The boundary conditions on ρ b and g b in (6.17) will be utilised along with the relation ρ = ρ + ⟨g⟩ V that allows to link ρ and ρ in the interior of the domain. We will use a staggered grid in space following [START_REF] Lemou | Micro-Macro Schemes for Kinetic Equations Including Boundary Layers[END_REF] and a high order scheme in time, following the strategy developed previously. To ease the reading, only the first order version will be presented.

First, we present the space approximation based on a staggered grid. Let us consider the space interval [0, L] with two grids: x i = i∆x and x i+1/2 = (i + 1/2)∆x, ∆x = L/(N x -1). The 'interior' variables such as ρ, ρ are stored at grid points x i with i = 1, . . . , N x -2) and g is stored at i + 1/2 = 1/2, • • • , N x -3/2. We also use the variable g cl = ḡ ∪ ḡb ∈ R Nx+1 . The whole domain including boundary will be considered for the micro unknown ḡ so that the components of g cl correspond to the

grid indices i + 1/2 = -1/2, • • • , N x -1/2.
The matrices corresponding to spatial operators are given by Nx+1) , (6.21)

B - upw = 1 ∆x circ([-1, 1]) (Nx-1)×(Nx+1) , B + upw = 1 ∆x circ([0, -1, 1]) (Nx-1)×(
B cenρ = 1 ∆x circ([-1, 1]) (Nx-2)×(Nx-1) , B avg = 1 2 circ([1, 1]) (Nx-2)×(Nx-1) , (6.22) B ceng = 1 ∆x circ b ([-1, 1]) (Nx-1)×(Nx-2) . (6.23)
The circ b definition is presented in Appendix A. Further, we also introduce a vector containing the boundary values of ρ as ρ bd = 1 ∆x -ρ bi=0 , 0, 0, ..., 0, ρ b i=Nx -1 T (Nx-1)×1 .

We now present our scheme by using this matrix notation. For simplicity, we assume that ρ bd is time invariant. We also use the following notations:

T h = I -Π -v + B - upw + v -B + upw h, D ϵ,∆t = v ϵ 2 I -∆t L -1 ∆t I -Π -(vM )⟩ V , E ϵ,∆t = ϵ 2 I -∆t L -1 ∆t I -Π -(vM ) V , I ϵ,∆t = ϵ 2 I -∆t L -1 , J = I -Π -(vM ).
The micro equation (6. [START_REF] Klar | Asymptotic-induced scheme for nonstationary transport equations in the diffusive limit[END_REF]) is discretised in time as in the previous (periodic) case (6.24)

g n+1 = I ϵ,∆t ϵ 2 g n -ϵ∆tT g n cl -ϵ∆tJ B ceng ρ n+1 -ϵ∆tJ ρ bd ,
and for the macro equation (6.20), we obtain

ρ n+1 -ρ n ∆t + 1 ϵ vB cenρ g n+1 V = 0
Substituting g n+1 in the above equation, we get (6.25)

ρ n+1 = ρ n -∆tB cenρ vI ϵ,∆t ϵg n -∆tT g n cl -∆tJ B ceng ρ n+1 -∆tJ ρ bd V .
In index notation, we use

ρ n+1 i = ρ n+1 i + 1 2 ⟨g n+1 i-1/2 + g n+1 i+1/2 ⟩ V (since ρ = ρ + ⟨g⟩ V ) to
match the two grids. In matrix notation, this becomes ρ n+1 = ρ n+1 + B avg ⟨g n+1 ⟩ V with B avg given by (6.22). Substituting this into the above equation and inserting the expression for g n+1 into B avg g n+1 V enable to update the interior macro unknown (6.26)

ρ n+1 = I -ϵB avg E ϵ,∆t B ceng -∆tB cenρ D ϵ,∆t B ceng -1 × ρ n -B avg I ϵ,∆t ϵ 2 g n -ϵ∆tT g n cl -ϵ∆tJ ρ bd V -∆tB cenρ vI ϵ,∆t ϵg n -∆tT g n cl -∆tJ ρ bd V .
The right hand side of above expression involves only known quantities so that ρ n+1 can be updated from (6.26) which can then be used to update g n+1 in (6.24). Then, we update g n+1 cl thanks to the boundary conditions (6.17), and finally ρ n+1 can be computed from ρ n+1 = ρ n+1 + B avg g n+1 V . In the limit ϵ → 0, the above equation becomes,

ρ n+1 = I + ∆tB cenρ v ⊗ L-1 J V B ceng -1 ρ n -∆tB cenρ v ⊗ L-1 J V ρ bd
This is a consistent discretization of the diffusion equation in (2.5) since ⟨v⊗ L-1

J ⟩ V = ⟨v ⊗ L -1 (vM )⟩ V = -κ.
Further, the high order scheme in time can be constructed in a similar manner as before.

7. Numerical results. In this section, we present the numerical validation of our high order asymptotic preserving schemes in different configurations.

7.1. Diffusion asymptotics. First, we check time and space accuracy for the micro-macro scheme in the diffusion limit.

ρ(0, x) = 1 + cos(x) Well-prepared data (WP): g(0, x, v) = ϵ 2 (I -Π) v 2 M ρ(0, x) Non-well prepared data (N-WP): g(0, x, v) = (I -Π) v 2 M ρ(0, x), with M (v) = 1 √ 2π e -v 2 /2
. Periodic boundary conditions are used on both ρ and g.

The spatial terms are discretised by using the atmost-third order accurate matrices on non-staggered grid presented in subsection 5. In Figure 1, we plot the time error for the different time integrators in both WP and N-WP cases and for different values of ϵ. Note that the reference solution for each curve is obtained by using the same micro-macro scheme corresponding to that curve with ∆t = 10 -4 . For ϵ = 1, the required orders of accuracy are recovered for type A schemes with both N-WP and WP initial data, as observed in Figures 1a and1b. For ϵ = 10 -4 , due to the asymptotic degeneracy of our scheme into a fully- sections, usage of N-WP initial data for CK-ARS time integrators does not allow the asymptotic accuracy (Figure 1g), as discussed in [START_REF] Dimarco | Asymptotic Preserving Implicit-Explicit Runge-Kutta Methods for Nonlinear Kinetic Equations[END_REF].

Since we proved the asymptotic preserving property, the diffusion solution is used as reference solution in the asymptotic regime (ϵ = 10 -4 ) with ∆t = 10 -4 (in Figure 2)

to check the orders of accuracy of high order integrators. The results are similar to the ones obtained for ϵ = 10 -4 in Figure 1, except that here we observe a plateau for third order scheme and small ∆t. This is due to the O(ϵ 2 ) difference between the schemes based on micro-macro and diffusion models. This error dominates O(∆t 3 ) error, and hence it is observed. 7.1.2. Space order of accuracy. The problem set-up is the same as described in the previous subsection, except for the following changes. Here, we consider the final time to be T = 0.01 and ∆t = 0.001 so that the error in time is small enough to study the spatial accuracy. To do so, we consider the following number of points in space: N x = 20, 24, 30, 40 and 60. The reference solution is obtained with N x = 120.

Since the spatial accuracy plots obtained from different time integrators are quite similar, we present only the plots obtained by using DP1-A(2, 4, 2) and ARS [START_REF] Carpenter | Additive runge-kutta schemes for convection-diffusionreaction equations[END_REF][START_REF] Carpenter | Additive runge-kutta schemes for convection-diffusionreaction equations[END_REF][START_REF] Boscarino | Implicit-Explicit Runge-Kutta Schemes for Hyperbolic Systems and Kinetic Equations in the Diffusion Limit[END_REF] for different values of ϵ (ϵ = 10 -4 , 0.2, 1) in Figures 3a and3b. For the spatial discretization, we only show the results obtained by the third order spatial matrices on non-staggered grid presented in subsection 5.3 so that the scheme is expected to be third order accurate in space. In Figures 3a and3b, the expected order is observed for the two time integrators and for the three considered values of ϵ. 7.1.3. Qualitative results. In this part, we compare the density obtained by the micro-macro equation (MM), the linear kinetic equation with BGK collision operator (BGK) and the asymptotic diffusion equation, for different values of ϵ. The MM scheme described in previous sections is utilised, the BGK is discretized using an IMEX (implicit treatment of collision term and explicit treatment of transport term) scheme whereas for the diffusion model, an implicit scheme is used. For all three models, the Butcher tableau corresponding to DP1-A(2, 4, 2) time integrator is used.

For the spatial discretization, we use third order scheme on non-staggered grid.

This manuscript is for review purposes only. The problem domain L = [0, 2π] is discretised using N x = 20 grid points for all the three models. The final time is T = 0.5, and ∆t = 0.005. We use the same N-WP initial and boundary conditions described in the previous subsection. Further, we also consider the same velocity discretization as before for both MM and BGK models.

In Figure 4a for rarefied regime (ϵ = 1), the MM and BGK models compare very well, while the diffusion model is different as expected. In the intermediate regime (ϵ = 0.2), the BGK and MM models match very well while the diffusion model is slightly different. For ϵ = 10 -4 , we only compare MM and the diffusion in Figure 4c and illustrate the AP property of the time integrators used for MM.

Advection-diffusion asymptotics.

In this subsection, we present the time accuracy of our high order micro-macro scheme for the advection-diffusion case.

As in the diffusion case, the spatial domain L = [0, 2π] is discretised using N x = 20 grid points whereas the velocity domain is [-v max , v max ] with v max = 5 and ∆v = 1.

This manuscript is for review purposes only. The initial condition for the problem is: ρ(0, x) = sin(x) (7.1) Well-prepared data (WP): g(0, x, v) = ϵ 2 (I -Π) v 2 M ρ(0, x) (7.2) Non-well prepared data (N-WP): g(0, x, v) = (I -Π) v 2 M ρ(0, x), (7.3) with M (v) = 1 √ 2π e -v 2 /2 . Periodic boundary conditions are used on both ρ and g.

The spatial terms are discretised by using the atmost-first order accurate matrices on staggered grid presented in subsection 5.2. The final time is T = 0.5, and the following time steps are considered: ∆t = 0.5, 0.1, 0.05, 0.01, 0.005, 0.001. We observe the time order of accuracy for both ϵ = 1 and ϵ = 10 -4 . We choose the highest order time integrator in both type A and CK-ARS schemes for studying the time accuracy.

Hence, we consider DP1-A(2, 4, 2) and ARS(4, 4, 3) with N-WP and WP data respectively.

Asymptotically, our micro-macro scheme degenerates to a consistent scheme for the advection-diffusion equation with advection and diffusion terms being treated explicitly and implicitly respectively. Hence, unlike the case of diffusion asymptotics for which an extra order is observed asymptotically, DP1-A(2, 4, 2) remains second order accurate for ϵ = 10 -4 since both explicit and implicit matrices of the Butcher tableau are involved here (Figure 5a). For ϵ = 1, the required second order accuracy is observed. Further, the required third order accuracy of ARS [START_REF] Carpenter | Additive runge-kutta schemes for convection-diffusionreaction equations[END_REF][START_REF] Carpenter | Additive runge-kutta schemes for convection-diffusionreaction equations[END_REF][START_REF] Boscarino | Implicit-Explicit Runge-Kutta Schemes for Hyperbolic Systems and Kinetic Equations in the Diffusion Limit[END_REF] is observed for both ϵ = 10 -4 , 1 in Figure 5b, since well-prepared initial data is considered.

This manuscript is for review purposes only. interior of the domain. This is consistent with the results observed in the literature [START_REF] Klar | Asymptotic-induced domain decomposition methods for kinetic and drift diffusion semiconductor equations[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] Crouseilles | An asymptotic preserving scheme based on a micro-macro decomposition for collisional Vlasov equations: diffusion and high-field scaling limits[END_REF].

  (ã ij ) and A = (a ij ) are s×s matrices which correspond to the explicit and implicit parts of the scheme (A and à respectively are lower triangular and strictly lower triangular matrices). The coefficients c and c are given by ci = i-1 j=1 ãij , c i = i j=1 a ij , and the vectors b = ( bj ) and b = (b j ) give quadrature weights that combine the stages. For GSA IMEX RK scheme, we have (3.5) c s = cs = 1 and a sj = b j , ãsj = bj , ∀j ∈ {1, 2.., s}.

2 )

 2 follows the type CK-ARS double Butcher tableau (known as ARS(1

Theorem 4 . 3 .

 43 Consider the scheme (3.10)-(3.11) approximating the macro-micro model (2.7)-(2.8), with the RK method (3.4) of type A or of type CK-ARS (with well-prepared initial data g 0 = O(ϵ)). Then in the limit ϵ → 0, the scheme (3.10)-(3.11)

7. 1 . 1 .

 11 Time order of accuracy. The spatial domain L = [0, 2π] of the problem is discretized using N x = 50 grid points. The velocity domain is truncated to [-v max , v max ] with v max = 5 and we take ∆v = 1. The initial condition is:

3 .

 3 The final time is T = 0.5, and the following ∆t are considered to validate the different high order time integrators: ∆t = 0.5, 0.1, 0.05, 0.01, 0.005, 0.001. The type A micro-macro schemes constructed using the Butcher tableau corresponding to DP-A(1, 2, 1), DP2-A(2, 4, 2) and DP1-A(2, 4, 2) are considered. Although DP1-A(2, 4, 2) is second order accurate, the implicit part of it when used separately is third order accurate. Further, we also consider the type CK-ARS micro-macro schemes constructed using Butcher tableau corresponding to ARS(1, 1, 1), ARS(2, 2, 2) and ARS[START_REF] Carpenter | Additive runge-kutta schemes for convection-diffusionreaction equations[END_REF][START_REF] Carpenter | Additive runge-kutta schemes for convection-diffusionreaction equations[END_REF][START_REF] Boscarino | Implicit-Explicit Runge-Kutta Schemes for Hyperbolic Systems and Kinetic Equations in the Diffusion Limit[END_REF]. The Butcher tableau of different time integrators utilised are presented in Appendix B.

  Figures 1e and 1f. However, for ϵ = 10 -4 , orders of accuracy are observed only when WP initial data are used (Figure1h). As shown in the analyses presented in previous

4 Fig. 1 .

 41 Fig. 1. Accuracy in time for different type A and CK-ARS time integrators (both WP and N-WP initial data). The reference solution is obtained from the micro-macro with ∆t = 10 -4 .

Fig. 2 .Fig. 3 .

 23 Fig. 2. Accuracy in time for different type A and CK-ARS time integrators (both WP and N-WP initial data). The reference solution is obtained from the diffusion equation with ∆t = 10 -4 .

4 Fig. 4 .Fig. 5 .

 445 Fig. 4. Qualitative results for diffusion asymptotics

4 Fig. 8 .

 48 Fig. 8. Qualitative results for non-equilibrium inflow at the left boundary.

This manuscript is for review purposes only.

Acknowledgement. Megala Anandan sincerely acknowledges the indispensable support provided by Prof. S. V. Raghurama Rao, whose encouragement was instrumental in the successful completion of this research work.

Fig. 6. Accuracy in time with type A schemes for ϵ = 1 (left) and ϵ = 10 -4 (right). The reference solution is obtained from the micro-macro for inflow boundaries scheme with ∆t = 10 -4 .

Inflow boundary condition.

In this subsection, the high order numerical scheme for micro-macro model that allows inflow boundary conditions is validated numerically. We first present the time accuracy results for high order schemes. Then, some qualitative plots are shown for two tests with zero inflow at the right boundary, and equilibrium and non-equilibrium inflows respectively at the left boundary.

Time order of accuracy.

If the domain of the problem is a half-plane, ω(x, v) = -v, 0, 0, • • • can be chosen ∀x as described in [START_REF] Lemou | Micro-Macro Schemes for Kinetic Equations Including Boundary Layers[END_REF]. Here, for numerical purposes, we consider a domain of L = [0, 2] and assume that the right boundary does not influence the dynamics.

The spatial domain is discretised using N x = 20 grid points and the velocity domain is [-v max , v max ] with v max = 5 with ∆v = 1. The initial conditions at all interior points and right boundary conditions for the variables ρ, ρ and g are considered to be 0. The left boundary conditions (for v k > 0) are:

. The final time is T = 0.1, and the following time steps are considered to check the accuracy in time: ∆t = 0.1, 0.05, 0.01, 0.005, 0.001. Like in the previous problems, we observe the time order of accuracy for both ϵ = 1 and ϵ = 10 -4 .

The time integrators considered are DP-A(1, 2, 1) and DP1-A(2, 4, 2). The reference solution for each curve in Figure 6 is obtained by using the same micro-macro scheme corresponding to that curve with ∆t = 10 -4 . For type A time integrators with ϵ = 1 in Figure 6a, first and second order accuracies of DP-A(1, 2, 1) and DP1-A(2, 4, 2) are observed. In Figure 6b for ϵ = 10 -4 , first and third order accuracies of DP-A(1, 2, 1) and DP1-A(2, 4, 2) respectively are observed. As for the (periodic) diffusion case, DP1-A(2, 4, 2) turns out to be third order accurate since only the implicit part of Butcher tableau is involved asymptotically. For ARS(2, 2, 2) and ARS(4, 4, 3) time integrators (not shown here), order reduction to first order for ϵ = 1 (due to the initial condition). However, for ϵ = 10 -4 , the required second and third orders respectively are observed. we consider the same velocity discretization as before for both MM and BGK models.

The results for MM are obtained by DP1-A(2, 4, 2) time integrator.

In Figure 7a for rarefied regime (ϵ = 1), the MM and BGK results are in good agreement. In the intermediate regime (ϵ = 0.4) in Figure 7b, the MM and BGK results are still close, and still different from the diffusion one. For ϵ = 10 -4 , only MM and the diffusion are plotted and are found to be in very good agreement, thereby illustrating the AP property of the numerical scheme for MM.

7.3.3.

Qualitative results for non-equilibrium inflow. In this part, we consider the same problem as before, but the left boundary condition is chosen as (for

The number of grid points, velocity discretization, final time and time step are the same as in the previous (equilibrium inflow) case. Here, we present a comparison of plots obtained by using schemes based on MM, BGK and diffusion models, for different regimes of ϵ. The scheme described in subsection 6.2.1 is used for the micro-macro model and a standard BGK approximation where only inflow boundary condition is needed serves as a reference. For diffusion, the diffusion term is treated implicitly and the left boundary condition for diffusion model is obtained from [START_REF] Klar | Asymptotic-induced domain decomposition methods for kinetic and drift diffusion semiconductor equations[END_REF] which translates in our context

In Figure 8a for rarefied regime (ϵ = 1), the MM and BGK models compare very well, while the diffusion model is driven by the macro boundary condition. In the intermediate regime (ϵ = 0.4) in Figure 8b, in the MM and BGK results (which are in a good agreement), a boundary layer starts to be created whereas it is not the case for the diffusion model. For ϵ = 10 -4 , one can see that MM model develops a boundary layer at the left boundary before aligning with the diffusion model in the This manuscript is for review purposes only.

The following is the 4-stage third order accurate Butcher tableau ARS(4, 4, 3):

For type A, we use 2-stage first order accurate Butcher tableau DP-A(1, 2, 1) (γ ≥ 1 2 )

This manuscript is for review purposes only.

HIGH ORDER AP SCHEME FOR KINETIC EQUATIONS IN DIFFUSIVE SCALING 25

The following is the 4-stage second order accurate Butcher tableau DP2-A(2, 4, 2):

The following is the 4-stage second order accurate Butcher tableau DP1-A(2, 4, 2)

which achieves third order accuracy on the DIRK part: