A particle micro-macro decomposition based numerical scheme for collisional kinetic equations in the diffusion scaling

TL;DR

This paper introduces a particle-based numerical scheme for linear kinetic equations in the diffusion limit, reformulating the micro-macro system to reduce noise and computational cost, while maintaining asymptotic-preserving properties.

Contribution

The paper develops a non-stiff reformulation of the micro-macro system enabling standard particle approximations, improving efficiency and noise reduction in kinetic equation simulations.

Findings

Scheme maintains asymptotic-preserving property.

Significantly reduces noise compared to traditional methods.

Computational cost decreases as the system approaches the diffusion limit.

Abstract

In this work, we derive particle schemes, based on micro-macro decomposition, for linear kinetic equations in the diffusion limit. Due to the particle approximation of the micro part, a splitting between the transport and the collision part has to be performed, and the stiffness of both these two parts prevent from uniform stability. To overcome this difficulty, the micro-macro system is reformulated into a continuous PDE whose coefficients are no longer stiff, and depend on the time step $\Delta t$ in a consistent way. This non-stiff reformulation of the micro-macro system allows the use of standard particle approximations for the transport part, and extends the work in [5] where a particle approximation has been applied using a micro-macro decomposition on kinetic equations in the fluid scaling. Beyond the so-called asymptotic-preserving property which is satisfied by our schemes, they significantly reduce the inherent noise of traditional particle methods, and they have a computational cost which decreases as the system approaches the diffusion limit.

Equations

$${\partial }_{t}f+\frac{1}{\epsilon }v{\partial }_{x}f=\frac{1}{{\epsilon }^2}(\rho M-f),f(t=0,x,v)=f_0(x,v),$$

$${\partial }_t\bar{\rho }-\frac{1}{3}{\partial }_{xx}\bar{\rho }=0,\bar{\rho }(t=0,x)=\frac{1}{2}{\int }_V{f}_0(x,v)\text{d}v.$$

$${\partial }_{t}f+\frac{1}{\epsilon }v{\partial }_{x}f+\frac{1}{\epsilon }E{\partial }_{v}f=\frac{1}{{\epsilon }^2}(\rho M-f),$$

$${\partial }_{x}E=\rho -1,$$

$${\partial }_t\bar{\rho }-{\partial }_x({\partial }_x\bar{\rho }-\bar{E}\bar{\rho })=0,\bar{\rho }(t=0,x)={\int }_{\mathbb{R}}{f}_0(x,v)\text{d}v,$$

$$\{\begin{array}{rl}{\partial }_t\rho & +\frac{1}{\epsilon }{\partial }_x⟨vg⟩=0,\\ {\partial }_{t}g& +\frac{1}{\epsilon }(I-\Pi )\left[v{\partial }_x(\rho M+g)+E{\partial }_v(\rho M+g)\right]=-\frac{1}{{\epsilon }^2}g.\end{array}$$

$${\partial }_{t}g+\frac{1}{\epsilon }\left[vM{\partial }_x\rho -vME\rho +(I-\Pi )(v{\partial }_{x}g+E{\partial }_{v}g)\right]=-\frac{1}{{\epsilon }^2}g.$$

$${\partial }_t\rho -{\partial }_x\left[⟨v^{2}M⟩{\partial }_x\rho -⟨v^{2}M⟩E\rho \right]=\mathcal{O}(\epsilon ),$$

$${\partial }_t\bar{\rho }-{\partial }_{xx}\bar{\rho }+{\partial }_x(\bar{E}\bar{\rho })=0,{\partial }_x\bar{E}=\bar{\rho }-1.$$

$$\{\begin{array}{rl}{\partial }_t\rho & +\frac{1}{\epsilon }{\partial }_x⟨vg⟩=0,\\ {\partial }_{t}g& +\frac{1}{\epsilon }(I-\Pi )\left[v{\partial }_x(\rho M+g)\right]=-\frac{1}{{\epsilon }^2}g,\end{array}$$

$${\partial }_t\bar{\rho }-\frac{1}{3}{\partial }_{xx}\bar{\rho }=0,$$

$${\partial }_t(e^{t/{\epsilon }^2}g)=-\frac{e^{t/{\epsilon }^2}}{\epsilon }\mathcal{F}\left(\rho ,g\right),$$

$$\mathcal{F}\left(\rho ,g\right)=v{\partial }_x\rho +v{\partial }_{x}g-{\partial }_x⟨vg⟩.$$

$$g(t^{n+1})=e^{-\Delta t/{\epsilon }^2}g(t^n)-\epsilon (1-e^{-\Delta t/{\epsilon }^2})\mathcal{F}\left(\rho (t^n),g(t^n)\right)+\mathcal{O}(\Delta t^2).$$

$$\frac{g(t^{n+1})-g(t^n)}{\Delta t}=\frac{e^{-\Delta t/{\epsilon }^2}-1}{\Delta t}g(t^n)-\epsilon \frac{1-e^{-\Delta t/{\epsilon }^2}}{\Delta t}\mathcal{F}\left(\rho (t^n),g(t^n)\right)+\mathcal{O}(\Delta t^2),$$

$${\partial }_{t}g(t^n)=\frac{e^{-\Delta t/{\epsilon }^2}-1}{\Delta t}g(t^n)-\epsilon \frac{1-e^{-\Delta t/{\epsilon }^2}}{\Delta t}\mathcal{F}\left(\rho (t^n),g(t^n)\right),\forall n.$$

$${\partial }_t\rho$$

$${\partial }_{t}g$$

$$g\left(t,x,v\right)=\sum _{k=1}^{N_p}{\omega }_k(t)\delta (x-x_k(t))\delta (v-v_k(t)),$$

$${\omega }_k(t)=g(t,x_k(t),v_k(t))\frac{L_x{L}_v}{N_p}.$$

$${\partial }_{t}g=\epsilon \frac{1-e^{-\Delta t/{\epsilon }^2}}{\Delta t}v{\partial }_{x}g,$$

$${\dot{x}}_k=\epsilon \frac{1-e^{-\Delta t/{\epsilon }^2}}{\Delta t}{v}_k,$$

$${\partial }_{t}g=\frac{e^{-\Delta t/{\epsilon }^2}-1}{\Delta t}g-\epsilon \frac{1-e^{-\Delta t/{\epsilon }^2}}{\Delta t}\left[v{\partial }_x\rho -{\partial }_x⟨vg⟩\right],$$

$${\dot{\omega }}_k=\frac{e^{-\Delta t/{\epsilon }^2}-1}{\Delta t}{\omega }_k-\epsilon \frac{1-e^{-\Delta t/{\epsilon }^2}}{\Delta t}\left[v_k({\partial }_x\rho (x_k)-{\partial }_x⟨vg⟩(x_k)\right]\frac{L_x{L}_v}{N_p}.$$

$$x_k^{n+1}=x_k^n+\epsilon (1-e^{-\Delta t/{\epsilon }^2})v_k.$$

$$⟨vg⟩(\text{mbox}{x}_i)\approx \sum _{k=1}^{N_p}{\omega }_k^n{B}_{\ell }(\text{mbox}{x}_i-x_k^{n+1})v_k,$$

$$B_{\ell}(x)=(B_{0}*B_{\ell-1})(x),\;\;\mbox{ with }\;\;B_{0}(x)=\left\{\begin{array}[]{llcc}\frac{1}{\Delta x}&\mbox{ if }|x|<\Delta x/2,\\ 0&\mbox{ else}.\end{array}\right.$$

$${\omega }_k^{n+1}=e^{-\Delta t/{\epsilon }^2}{\omega }_k^n-\epsilon (1-e^{-\Delta t/{\epsilon }^2})\left[{\alpha }_k^n+{\beta }_k^n\right],$$

$${\alpha }_k^n=v_k{\partial }_x{\rho }^n(x_k^{n+1})\frac{L_x{L}_v}{N_p}\text{and}{\beta }_k^n=-{\partial }_x⟨vg⟩(x_k^{n+1})\frac{L_x{L}_v}{N_p}.$$

$${\partial }_x{\rho }^n(x_k^{n+1})\approx \sum _{i=1}^{N_x}\frac{{\rho }_{i+1}^n-{\rho }_{i-1}^n}{2\Delta x}{B}_{\ell }(\text{mbox}{x}_i-x_k^{n+1}).$$

$${\rho }_i^{n+1}={\rho }_i^n-\frac{\Delta t}{\epsilon }\frac{⟨v{g}^{n+1}{⟩}_{i+1}-⟨v{g}^{n+1}{⟩}_{i-1}}{2\Delta x},$$

Full text

1 Introduction

Particle systems appearing in plasma physics or radiative transfer can be described at different scales. When the system is far from its thermodynamical equilibrium, a kinetic description is necessary. Particles are then represented by a distribution function $f$ which depends on time $t\geq 0$, position $x\in\mathbb{R}^{d}$ and velocity $v\in V\subset\mathbb{R}^{d}$, $d\geq 1$. The distribution $f\left(t,x,v\right)$ satisfies a collisional kinetic equation. Particle methods are often used for simulating kinetic problems, especially in realistic $3$-dimensional situations, $d=3$. However, they are affected by numerical noise due to their probabilistic character. A simple way to reduce this noise is to increase the number of particles, but then the numerical cost increases as well. Other standard kinetic descriptions, as phase space grid methods, may require too much memory in the two or three dimensional framework. Otherwise, macroscopic descriptions depending only on $t$ and $x$ can be sufficient if the system stays near its thermodynamical equilibrium, and are less expensive since their unknown does not depend on the velocity variable anymore. Beside the noisy character of standard particle methods, there is an additional difficulty in kinetic descriptions which is linked to the presence of various scales in the system. Multi-scale phenomena may indeed appear in plasma devices or radiative transfer applications, depending on some physical parameters as for example the mean free path of particles or the Knudsen number denoted here by $\varepsilon$. This multi-scale character is often represented by stiff terms in the kinetic equation, and the general challenge is to construct efficient numerical methods for these multiscale kinetic equations: this means that, without numerically resolving the stiffness, the numerical method must solve accurately the kinetic regime, must have the right asymptotics in the high-stiffness limit (the so-called asymptotic preserving property) and its computational cost should decrease as the system approaches the equilibrium (a time diminishing property). Note that direct numerical methods whose parameters resolve the smallest scale of size $\varepsilon$ are impossible to use, since they automatically involve an extremely high computational cost.

Several strategies have been proposed to overcome this strong constraint. Domain decomposition methods can be applied when we have different regions with different values of the scaling parameter, see [9, 13, 22, 25]. When the different scales are less clearly delimited, we have to develop kinetic schemes that naturally reduce to good approximations of the macroscopic problem when the system goes near its equilibrium, and overcome the stiffness. Such schemes are often called Asymptotic-Preserving (AP), see [14, 17, 16, 12, 8, 18, 19, 20, 21, 15, 24, 4]. Mainly, the numerical cost remains comparable to the one of the non-stiff kinetic problem, even when $\varepsilon\ll 1$.

Our goal is to design an efficient AP scheme, using particles, for the following kinetic radiative transport equation (RTE) in the diffusion scaling

[TABLE]

where $x\in\Omega\subset\mathbb{R}$, $\rho(t,x)=\frac{1}{2}\int_{V}f(t,x,v)\textnormal{d}v$, $V=[-1,1]$, $M=1$ and $f_{0}(x,v)$ is a given initial condition. Periodic boundary conditions are considered. It is well-known (see [18, 11]) that when $\varepsilon$ goes to zero, the distribution function $f(t,x,v)$ converges towards $\bar{\rho}(t,x)M(v)$, where $\bar{\rho}$ satisfies the following diffusion equation

[TABLE]

An extension to the Vlasov-Poisson-BGK case is presented in Section 5. The kinetic equation is coupled to a Poisson equation for the electric field denoted by $E(t,x)$. More precisely, we consider

[TABLE]

where $x\in\Omega\subset\mathbb{R}$, $\rho(t,x)=\int_{V}f(t,x,v)\textnormal{d}v$, $V=\mathbb{R}$, $M\left(v\right)=\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{v^{2}}{2}\right)$ is the absolute Maxwellian and we consider periodic boundary conditions. Note that an additional condition $\int_{\Omega}E\textnormal{d}x=0$ is imposed to obtain a well-posed problem. When $\varepsilon$ goes to zero, the asymptotic model is a drift-diffusion equation satisfied by $\bar{\rho}(t,x)$ (see [1])

[TABLE]

where $\bar{E}$ is linked to $\bar{\rho}$ by the Poisson equation $\partial_{x}\bar{E}=\bar{\rho}-1$.

The strategy will be the use of the micro-macro decomposition (see [23, 21, 2, 8]). It consists in writing the distribution function as the sum of the equilibrium part and a rest. One can then derive a system of two equations: a kinetic one for the rest $g(t,x,v)$ and a macroscopic one for the equilibrium $\rho(t,x)M(v)$. AP micro-macro schemes for (1.1) have been proposed in [21, 2, 8]. These schemes consist in a semi-implicit phase space grid method for the kinetic part, coupled to a classical spatial grid method for the macro part. Our strategy in this work follows the strategy of [5] in the case of a fluid scaling: we use particles to sample the kinetic part whereas an Eulerian solver is used to discretize the macro unknown. The main motivation of this strategy lies in the fact that the micro part $g$ converges to zero when $\varepsilon$ goes to zero, so that a very few number of particles can sample it. As a consequence, in this regime, the cost of the global micro-macro solver is almost the same as the cost of an asymptotic solver for (1.2).

In this work, we focus on a diffusion type scaling (as in (1.1) or (1.3)) so that an additional scale is involved compared to the fluid scaling considered in [5]. In [21, 8], a diffusion scaling was studied, but using a fully grid based solver. Hence, the stiffest term (of order $1/\varepsilon^{2}$) is considered implicit in time in the micro equation, which enables to stabilize the transport term (of order $1/\varepsilon$) and then to derive an AP scheme for (1.1) and (1.3). The use of particles for the micro part prevents from a similar strategy since a splitting between the transport term (of order $1/\varepsilon$) and the source term needs to be done. Then, a uniform stable scheme is hard to obtain in this context. To overcome this difficulty, a suitable formulation of the original model (1.1) is performed so that the stiff transport term $(1/\varepsilon)\,v\partial_{x}g$ becomes $\varepsilon/\Delta t(1-e^{-\Delta t/\varepsilon^{2}})v\partial_{x}g$. This reformulation is correct up to $\Delta t^{2}$ (for fixed $\varepsilon>0$) and has the good behavior when $\varepsilon$ goes to zero (for fixed $\Delta t>0$). This formulation is the starting point of the design of micro-macro-particle based numerical schemes which enjoy the AP property and for which the numerical cost diminishes as $\varepsilon$ goes to zero. This approach is extended to the second-order (in time) and to the Vlasov-Poisson-BGK case (1.3)-(1.4).

The sequel of the paper is organized as follows. In Section 2, we recall the formal derivation of the asymptotic model of (1.3) and (1.1). The first-order (in time) reformulation of (1.1) is presented in Subsection 3.1 and its Lagrangian discretization in Subsection 3.2. Its extension to a second-order in time model is detailed in Section 4: the continuous model is presented in Subsection 4.1 and its discretization is developed in Subsection 4.2. Section 5 proposes an extension of our strategy to the Vlasov-Poisson-BGK system. Finally, Section 6 is devoted to numerical simulations.

2 Diffusion asymptotics

In this section, we recall the main steps of the derivation of the model obtained from (1.3) when $\varepsilon$ goes to zero. To do so, we consider the micro-macro decomposition (see [21, 2, 23]) of $f$: $f(t,x,v)=\rho(t,x)M(v)+g(t,x,v)$, with $\rho(t,x)=\langle f\rangle$, $M(v)$ is the Maxwellian equilibrium and the rest $g$ satisfies $\langle g\rangle=0$. Here $\langle f\rangle=\int_{V}f(v)\textnormal{d}v$, with $V=\mathbb{R}$ and we also use the notation $\Pi f=\langle f\rangle M$. The following micro-macro model is equivalent to the original model (1.3)

[TABLE]

Since $\langle vM\rangle=0$ and $\partial_{v}M=-vM$, the micro equation can be rewritten as

[TABLE]

When $\varepsilon$ goes to zero, one gets from (2.2): $g=-\varepsilon(vM\partial_{x}\rho-vME\rho)+{\cal O}(\varepsilon^{2})$. Then the macro equation becomes

[TABLE]

which gives, using $\langle v^{2}M\rangle=1$ the following drift-diffusion equation satisfied by the limit $\bar{\rho}$

[TABLE]

The same calculations enables to derive the micro-macro model equivalent to (1.1)

[TABLE]

from which we derive the corresponding asymptotic model

[TABLE]

with $M=1$, $V=\left[-1,1\right]$ and $\langle f\rangle=\frac{1}{2}\int_{V}f\textnormal{d}v$.

3 First-order in time reformulation and its discretization

In this part, a first-order reformulation of the micro part is proposed, which enables to avoid the stiff transport term in space. The strategy is presented in the case of the equation (1.1) and its corresponding micro-macro model (2.4).

3.1 First-order in time reformulation

We start with (1.1) (with periodic boundary condition in space) and consider the micro-macro decomposition of $f=\rho+g$ (here $M(v)=1$ for all $v\in[-1,1]$) and the micro-macro model (2.4).

First, we rewrite the micro part of (2.4) as

[TABLE]

where $\mathcal{F}\left(\rho,g\right)$ is given by

[TABLE]

We denote $\Delta t>0$ the time step, $t^{n}=n\Delta t$ with $n\in\mathbb{N}$. Then, a second step consists in integrating (3.1) on $[t^{n},t^{n+1}]$ to get

[TABLE]

To derive a continuous (in time) equation, we make appear a discrete time derivative on the left-hand side

[TABLE]

which can be rewritten, up to terms of order ${\cal O}(\Delta t^{2})$, as

[TABLE]

We finally obtain the first-order reformulation of (1.1)

[TABLE]

with $\mathcal{F}\left(\rho,g\right)$ given by (3.2). We remark that the micro equation does not contain any stiff term and then has a suitable form for a numerical discretization using particles. Moreover, this first-order reformulation satisfies the following properties.

• •

Consistency: For all fixed $\varepsilon>0$, equation (3.5) is consistent with the initial micro equation (2.2) as $\Delta t$ goes to zero.

• •

Asymptotic behaviour: For all fixed $\Delta t>0$, as $\varepsilon$ goes to zero, we get from (3.5) $g=-\varepsilon v\partial_{x}\rho+O(\varepsilon^{2})$, which injected in the macro equation (3.4) provides the limit model (1.2).

3.2 Lagrangian discretization

This subsection is devoted to the derivation of an AP-particle based numerical scheme for (3.4)-(3.5).

Explicit in time discretization

We propose now a Lagrangian discretization of (3.5). More precisely, we adopt a particle method, see [3], and consider a set of $N_{p}\in\mathbb{N}$ macro particles. The position of particle $k$, $1\leq k\leq N_{p}$, is denoted by $x_{k}(t)\in\Omega=\left[0,L_{x}\right]$, with $L_{x}>0$, its velocity by $v_{k}(t)\in V=\left[-1,1\right]$ and its weight by $\omega_{k}(t)\in\mathbb{R}$. Let $L_{v}=|V|=2$. The function $g$ is then assumed to be of the form

[TABLE]

where $\delta$ denotes the Dirac mass function. Weights $\omega_{k}(t)$ are related to the distribution function $g$ through

[TABLE]

Initially, particles are randomly distributed in the phase-space domain $\left[0,L_{x}\right]\times V$ and their weights are computed following (3.7).

The density $\rho$ is computed on a uniform spatial grid defined by $\mbox{x}_{i}=i\Delta x$, $i=0,\dots,N_{x}$, $N_{x}\in\mathbb{N}^{\star}$ and $\Delta x=L_{x}/N_{x}$. We denote by $\rho_{i}^{n}$ the approximation at time $t^{n}=n\Delta t$ and position $\mbox{x}_{i}$ of $\rho(t^{n},\mbox{x}_{i})$, with $\Delta t>0$ the time step. Moreover, $g^{n}(x,v)\approx g(t^{n},x,v)$, $x_{k}^{n}\approx x_{k}(t^{n})$, $v_{k}^{n}\approx v_{k}(t^{n})$ and $w_{k}^{n}\approx w_{k}(t^{n})$. Let us remark that $\dot{v}_{k}=0$, so that the velocities $v_{k}(t)$ are constant in time and we will note $v_{k}^{n}=v_{k}^{0}=:v_{k}$ for all $n$.

Our goal is then to extend the particle discretization of [5] to diffusion scaling. To that purpose, we exploit the reformulation (3.5). As already said in [5], we have to use a splitting procedure between the transport part and the source part. Then, the (first-order) splitting writes

• •

start with an initial repartition of the $N_{p}$ particles $(x_{k}^{0},v_{k}^{0})$, with $\omega_{k}^{0}=g(t=0,x_{k}^{0},v_{k}^{0})L_{x}L_{v}/N_{p}$,

• •

solve the transport part

[TABLE]

with the (non stiff) characteristics

[TABLE]

• •

solve the source part

[TABLE]

using the equation satisfied by the weights

[TABLE]

Now, we detail the time discretization of the two steps. First, (3.8) is approximated by a simple forward Euler scheme

[TABLE]

Second, we compute the last term in (3.9). The term $\langle vg\rangle$ is approximated on the spatial grid $\mbox{x}_{i}$ using

[TABLE]

where $B_{\ell}\geq 0$ is a B-spline function of order $\ell$:

[TABLE]

We then approximate the equation on the weights (3.9) using a first-order explicit integrator

[TABLE]

with

[TABLE]

To compute $\alpha_{k}^{n}$ (resp. $\beta_{k}^{n}$), since $\rho^{n}$ (resp. $\langle vg\rangle$) is known on the spatial grid, we approximate $\partial_{x}\rho^{n}$ (resp. $\partial_{x}\langle vg\rangle$) by centered finite differences and evaluate at $x_{k}^{n+1}$ using an interpolation with B-spline functions, for example

[TABLE]

Finally, the macro equation (3.4) is advanced through

[TABLE]

where $\langle vg^{n+1}\rangle_{i}$ is computed using (3.11). Then, we have the following proposition.

Proposition 3.1**.**

The scheme given by (3.10)-(3.13)-(3.15) enjoys the AP property, i.e. it satisfies the following properties

• •

for fixed $\varepsilon>0$, the scheme is a first-order (in time) approximation of the original model (1.1),

• •

for fixed $\Delta t>0$, the scheme degenerates into an explicit first-order (in time) scheme of (1.2).

Proof.

The consistency follows directly from standard approximation. For the asymptotic behavior, when $\varepsilon$ goes to zero, we get $\omega_{k}^{n+1}=-\varepsilon\alpha_{k}^{n}+{\cal O}(\varepsilon^{2})$ (since $\omega_{k}^{n}={\cal O}(\varepsilon)$ $\forall n\geq 1$). Computing the momentum of $g^{n+1}$ means that we use (3.11) with $g^{n+1}$, or in the limit regime

[TABLE]

Injecting in the macro equation (3.4) then leads to a consistent discretization of (1.2). ∎

Remark 3.1** (Preservation of the micro-macro structure.).**

As detailed in [5], we have to correct the particle’ weights in order to preserve at the numerical level the micro-macro structure. Indeed, the micro-macro decomposition technique uses the zero-mean property $\langle g\rangle=0$. But nothing guarantees that this property is satisfied at the discrete level (on the weights $\omega_{k}$). That is why we have to correct the weights, by applying a discrete approximation of the operator $\left(I-\Pi\right)$ to each weight $\omega_{k}$, which is consistent with the continuous model.

We do not give the details of this correction (called projection step in following algorithms) and refer the reader to [5].

Implicit time discretization

In this part, we want to make the previous scheme (3.10)-(3.13)-(3.15) degenerate into an implicit time discretization of (1.2) (as done in [19, 8]). To do so, we decompose (3.13) into two parts so that the macro flux becomes

[TABLE]

with $B_{\ell}$ given by (3.12) and

[TABLE]

Since $\alpha_{k}^{n}$ is the weight of $v\partial_{x}\rho^{n}$ (see (3.14)) and $\langle v^{2}\rangle=1/3$, we can write

[TABLE]

so that the macro scheme can be

[TABLE]

We can go further by considering now the diffusion term implicit in time to get

[TABLE]

We can write the following proposition.

Proposition 3.2**.**

The scheme given by (3.10)-(3.13)-(3.17) enjoys the AP property, i.e. it satisfies the following properties

• •

for fixed $\varepsilon>0$, the scheme is a first-order (in time) approximation of the original model (1.1),

• •

for fixed $\Delta t>0$, the scheme degenerates into an implicit first-order (in time) scheme of (1.2).

Remark 3.2**.**

Instead of an implicit scheme, it is possible to use a Crank-Nicolson method.

The scheme is finally summarized in the following algorithm.

Algorithm 3.1**.**

• •

Initialize $(x_{k}^{0},v_{k}^{0})$, $\omega_{k}^{0}$ and $\rho_{i}^{0}$.

At each time step:

• •

1) Advance micro part:

• –

advance the characteristics with (3.10),

• –

compute $\langle vg\rangle$ with (3.11),

• –

advance the equation on the weights with (3.13).

• •

2) Projection step: compute $(I-\Pi)g^{n+1}$ using **[5]**.

• •

3) Advance macro part:

• –

compute $\langle vg^{n+1}\rangle$ with (3.11),

• –

compute $\rho^{n+1}$ with (3.17).

4 Second-order in time reformulation and its discretization

This section is devoted to the derivation of a second-order scheme for the micro-macro system (2.4). As for the first-order scheme, we will first reformulate the microscopic equation (2.2) in order to suppress stiff terms (see Subsection 3.1 for the first-order case), and then discretize the obtained micro-macro model to get an AP efficient numerical scheme (see Subsection 3.2 for the first-order case).

4.1 Second-order reformulation

Let start from the following (equivalent) reformulation of the micro part of (2.4)

[TABLE]

where $\mathcal{F}\left(\rho,g\right)$ is defined by (3.2). We now integrate with respect to $t\in[t^{n},t^{n+1}]$ and use a second-order mid-point quadrature

[TABLE]

To derive a continuous (in time) equation, we make appear a discrete time derivative on the left-hand side

[TABLE]

We now look for a continuous (in time) equation for which the previous relation is a second numerical scheme. To do so, we perform Taylor expansions of the different terms at $t^{n+1/2}$

[TABLE]

Finally, the microscopic equation of (2.4) is reformulated up to the second-order by

[TABLE]

and we can now consider the second-order reformulated micro-macro system

[TABLE]

4.2 Time discretization

We are now interested in the construction of an AP scheme for system (4.2)-(4.3), based on a second-order Runge-Kutta (RK2) method for the time discretization and a Lagrangian method for the phase space discretization of the micro part.

The RK2 is based on a prediction step on $\Delta t/2$ (first-order) and a correction step on $\Delta t$. Then, a second-order (in time) scheme for (4.2)-(4.3) would read

[TABLE]

We still have to fix $\widetilde{g}$ in (4.4) in order to get a second-order scheme and to ensure the convergence of $g^{n+1}$ to zero as $\varepsilon$ goes to zero. It turns out the choice $\widetilde{g}=\frac{g^{n}+g^{n+1}}{2}$ ensures the two conditions. Indeed, we get for the correction step of the micro part

[TABLE]

Up to now, the micro part converges exponentially fast to zero (when $\varepsilon$ goes to zero), so that the asymptotic behavior of the scheme is $\rho^{n+1}=\rho^{n}$. Hence, the last but not the least step consists in modifying the macro part to capture the correct asymptotic limit.

As done in Section 3, we will modify the discretization in order to make appear the diffusion term directly in the macro part (4.7). However, the modification has to be of order $\Delta t^{3}$ for a fixed $\varepsilon>0$ to not spoil the second-order accuracy of the scheme. The correction we propose consists in computing $\rho_{i}^{n+1}$ using a Crank-Nicolson scheme for the diffusion weighted which is of order ${\cal O}(\Delta t^{3})$ for fixed $\varepsilon>0$ and which degenerates to $\Delta t$ when $\varepsilon$ goes to zero

[TABLE]

4.3 Lagrangian discretization

We consider the same notations as in Subsection 3.2 and detail here the Lagrangian discretization of the micro-macro system (4.2)-(4.3). In the prediction step (4.4), we compute $x_{k}^{n+1/2}$ with a forward Euler integrator

[TABLE]

and advance the weights with

[TABLE]

We end this prediction step by computing the flux $\langle vg^{n+1/2}\rangle$ with (3.11) to get the density

[TABLE]

Now in the correction step, we compute the position at $t^{n+1}$ with

[TABLE]

then the weights are given by

[TABLE]

Now, using (4.8) in the last step, we have

[TABLE]

where $\langle vg^{n+1/2}\rangle_{i}$ is computed using (3.11).

We finally have the following result.

Proposition 4.1**.**

The scheme given by (4.9)-(4.10)-(4.11)-(4.12)-(4.13)-(4.14) enjoys the AP property, i.e. it satisfies the following properties

• •

for fixed $\varepsilon>0$, the scheme is a second-order (in time) approximation of the original model (1.1),

• •

for fixed $\Delta t>0$, the scheme degenerates into an implicit second-order (in time) scheme of (1.2).

Proof.

When $\varepsilon\to 0$, we get from (4.10) $w_{k}^{n+1/2}\to 0$ exponentially fast and then $\langle vg^{n+1/2}\rangle_{i}\to 0$. By injecting it in the macro equation (4.14), we have at the limit $\rho_{i}^{n+1}=\rho_{i}^{n}+\frac{\Delta t}{3}\partial_{xx}\left(\frac{\rho_{i}^{n+1}+\rho_{i}^{n}}{2}\right)$, which is a Crank-Nicolson discretization of the diffusion equation (1.2).∎

Remark 4.1**.**

Let us emphasize that the moments $\langle\cdot\rangle$ have to be computed with B-spline functions of order $\ell\geq 1$ in order to obtain a second-order in time scheme. Taking $\ell=0$ would lead to space discontinuities preventing the time scheme to be of second-order.

Remark 4.2**.**

The projection step being by construction a first-order step, we do it in the prediction step.

The scheme is finally summarized in the following algorithm.

Algorithm 4.1**.**

• •

Initialization of $(x_{k}^{0},v_{k}^{0})$, $\omega_{k}^{0}$ and $\rho_{i}^{0}$.

At each time step:

Prediction step: from $t^{n}$ to $t^{n+1/2}$.**

• •

1) Advance micro part:

• –

advance the characteristics with (4.9),

• –

compute $\langle vg\rangle$ with (3.11) and B-spline functions of order $\ell\geq 1$,

• –

advance the equation on the weights with (4.10).

• •

2) Projection step: compute $(I-\Pi)g^{n+1/2}$ using **[5]**.

• •

3) Advance macro part:

• –

compute $\langle vg^{n+1/2}\rangle$ with (3.11) and B-spline functions of order $\ell\geq 1$,

• –

compute the density with (4.11).

Correction step: from $t^{n}$ to $t^{n+1}$.**

• •

4) Advance micro part:

• –

advance the characteristics with (4.12),

• –

compute $\langle vg\rangle$ with (3.11) and B-spline functions of order $\ell\geq 1$,

• –

advance the equation on the weights with (4.13).

• •

5) Advance macro part with (4.14).

5 Extension to the Vlasov-Poisson-BGK case

This section is devoted to the extension of our method to kinetic equation making appear an electric field in the velocity direction. We consider the Vlasov-Poisson-BGK system in the diffusion scaling

[TABLE]

where $x\in\Omega=\left[0,L_{x}\right]\subset\mathbb{R}$, $\rho(t,x)=\int_{\mathbb{R}}f(t,x,v)\textnormal{d}v$ and $M\left(v\right)=\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{v^{2}}{2}\right)$ is the absolute Maxwellian. Let $f_{0}\left(x,v\right)=f\left(t=0,x,v\right)$ the initial distribution function and let consider periodic boundary conditions in $x$: $f\left(t,0,v\right)=f\left(t,L_{x},v\right)$, $\forall~{}v\in V$, $E\left(t,0\right)=E\left(t,L_{x}\right)$ $\forall t\geq 0$.

We can extend our schemes to this problem by adapting the computations of Subsections 4.1 and 4.2. We do not give all the details of the computations but insist on difficulties coming from the electric field term and write the resulting schemes.

The second-order reformulated micro-macro system corresponding to (1.3) is

[TABLE]

The limit model is here the drift-diffusion equation coupled to Poisson equation (2.3).

Equipped with this system (5.4)-(5.5), we construct the following Lagrangian, second-order in time scheme for (5.1)-(5.2)-(5.3). In the prediction step, characteristics are solved through

[TABLE]

the weights evolve with

[TABLE]

and the macro equation is advanced with

[TABLE]

In the correction step, characteristics are solved through

[TABLE]

the weights evolve with

[TABLE]

and the macro equation is advanced with

[TABLE]

The limit has been directly written in the macroscopic equation and the diffusion term is managed by a Crank-Nicolson method, in the prediction as well as in the correction step.

We have the following proposition.

Proposition 5.1**.**

The scheme given by (5.6)-(5.7)-(5.8)-(5.9)-(5.10)-(5.11) enjoys the AP property, i.e. it satisfies the following properties

• •

for fixed $\varepsilon>0$, the scheme is a second-order (in time) approximation of the original model (1.3),

• •

for fixed $\Delta t>0$, the scheme degenerates into an implicit second-order (in time) scheme of (1.5).

The scheme is finally summarized in the following algorithm.

Algorithm 5.1**.**

• •

Initialize $(x_{k}^{0},v_{k}^{0})$, $\omega_{k}^{0}$, and $\rho_{i}^{0}$.

• •

Compute $E_{i}^{0}$ thanks to FFT or finite differences.

At each time step:

Prediction step: from $t^{n}$ to $t^{n+1/2}$.**

• •

1) Advance micro part:

• –

advance the characteristics with (5.6),

• –

compute $\langle vg\rangle$ with (3.11) and B-spline functions of order $\ell\geq 1$,

• –

advance the equation on the weights with (5.7).

• •

2) Projection step: compute $(I-\Pi)g^{n+1/2}$ using **[5]**.

• •

3) Advance macro part:

• –

compute $\langle vg^{n+1/2}\rangle$ with (3.11) and B-spline functions of order $\ell\geq 1$,

• –

compute $\rho_{i}^{n+1/2}$ with (5.8),

• –

compute $E_{i}^{n+1/2}$ thanks to FFT or finite differences.

Correction step: from $t^{n}$ to $t^{n+1}$.**

• •

4) Advance micro part:

• –

advance the characteristics with (5.9),

• –

compute $\langle vg\rangle$ with (3.11) and B-spline functions of order $\ell\geq 1$,

• –

advance the equation on the weights with (5.10).

• •

5) Advance macro part:

• –

*compute $\rho_{i}^{n+1}$ with (5.11), *

• –

compute $E_{i}^{n+1}$ thanks to FFT or finite differences.

Remark 5.1**.**

We propose to use an upwind discretization of the derivative $\partial_{x}\left(E_{i}^{n+1/2}\rho_{i}^{n+1/2}\right)$:

[TABLE]

where the standard notations $u^{+}=\mathrm{max}(u,0)$ and $u^{-}=\mathrm{min}(u,0)$ are used. The same discretization is done for $\partial_{x}\left(E_{i}^{n}\rho_{i}^{n}\right)$ in the prediction step.

6 Numerical results

This section is devoted to some numerical experiments comparing the here designed micro-macro model with particles (denoted by MiMa-Part-1 for the first-order scheme and by MiMa-Part-2 for the second-order scheme) to the micro-macro Eulerian model (denoted by MiMa-Grid), the moment guided method (denoted by Moment G.) and (i) the Full PIC method in kinetic regimes ($\varepsilon$ of order 1) or (ii) the limit scheme in diffusion regime ($\varepsilon\ll 1$). MiMa-Part-1 corresponds to Proposition 3.2 and MiMa-Part-2 corresponds to Proposition 4.1. The micro-macro Eulerian scheme is presented in Appendix A. The moment guided method was first presented in [10] and is adapted to our context in Appendix B. The Full PIC method [3] consists in applying the particle representation (3.6) to the whole function $f$ (and not only to the perturbation $g$) and to solve the characteristics and the equations on weights coming from equation (1.1) or (1.3).

In the sequel, we consider three families of test cases: radiative transport equation (RTE) test cases with periodic boundary conditions (see equation (1.1)) in Subsection 6.1, Vlasov-Poisson-BGK test cases (see equations (1.3)-(1.4)) of Landau damping type in Subsection 6.2 and two-stream instability (TSI) test cases in Subsection 6.3.

6.1 RTE with periodic boundary conditions

We consider the RTE test case given by the initial condition

[TABLE]

with $M\left(v\right)=1,~{}\forall v\in\left[-1,1\right]$, and periodic boundary conditions in $x$.

We propose here to verify numerically the convergence of MiMa-Part-2 (presented in Subsection 4.3), Lagrangian in space, of second-order in time and with an implicit treatment of the diffusion term (see Algorithm 4.1). In Figure 1, we plot the error in $L^{\infty}$ norm of the density $\rho$ at time $T=0.1$ as a function of $\Delta t$ (from $10^{-4}$ to $0.1$) for the following parameters: $N_{x}=16$, $N_{p}=100$. For $\varepsilon=1$, $0.5$ and $0.1$, the reference solution is computed with MiMa-Part-2 using the same parameters but with $\Delta t=10^{-7}$. Whereas for $\varepsilon=10^{-6}$, the reference is a numerical solution of the diffusion equation (computed on a space grid, with $N_{x}=16$ and $\Delta t=10^{-7}$). In Figure 2, the error in $L^{\infty}$ norm is now represented as a function of $\varepsilon$ for different values of $\Delta t$: $10^{-1}$, $10^{-2}$, $10^{-3}$ and $10^{-4}$.

These plots confirm the fact that MiMa-Part-2 is second-order accurate in time for any fixed $\varepsilon>0$, and also when $\varepsilon\to 0$. However, for intermediate regimes (for instance $\varepsilon=0.1$), order reduction is observed. This is a classical observation for AP schemes. Note that similar behaviour is obtained with $L^{2}$ norm.

In Figure 3, we verify the AP property of the MiMa-Part-2 scheme and plot the density $\rho(T=0.1,x)$ as a function of $x$ for different values of $\varepsilon$: $1$, $0.25$, $10^{-2}$ and $10^{-6}$. We take fixed parameters: $N_{x}=64$, $\Delta t=10^{-3}$ and $N_{p}=10^{4}$. We compare the solutions obtained by MiMa-Part-2 to a numerical solution of the diffusion equation (computed on a space grid, with $N_{x}=512$ and $\Delta t=\Delta x^{2}$) and see that the $\varepsilon$-dependent solutions come closer to the diffusion one when $\varepsilon$ decreases.

Moreover, we illustrate in Figure 4 the fact that the cost of our method is very small at the limit. For that, we plot the density $\rho(T=0.1,x)$ as a function of $x$ for $\varepsilon=10^{-6}$ ($N_{x}=64$ and $\Delta t=10^{-2}$) and see that $N_{p}=100$ is sufficient to represent in a good way (without noise) the density. The numerical cost is then very close to the one of the asymptotic model.

6.2 Landau damping

In this subsection, we present Landau damping test cases in both regimes (kinetic when $\varepsilon=\mathcal{O}(1)$ and diffusive when $\varepsilon\to 0$). The initial distribution function is given by

[TABLE]

with the wave number $k=0.5$ and $\alpha=0.05$. For the micro-macro model (1.3), the initial condition is $\rho\left(t=0,x\right)=1+\alpha\cos\left(kx\right)$ and $g\left(t=0,x,v\right)=0$. For the limit drift-diffusion equation (1.5), we have $\rho\left(t=0,x\right)=1+\alpha\cos\left(kx\right)$.

We first verify the order in time of the MiMa-Part-2 scheme detailed in Section 5 and plot in Figure 5 the error in $L^{\infty}$ norm of the density $\rho$ at time $T=0.1$ as a function of $\Delta t$ (from $10^{-4}$ to $0.1$) for the following parameters: $N_{x}=16$, $N_{p}=100$. For $\varepsilon=1$, $0.5$ and $0.1$, the reference solution is computed with MiMa-Part-2 using the same parameters but with $\Delta t=10^{-7}$. Whereas for $\varepsilon=10^{-6}$ the reference is a numerical solution of the drift-diffusion equation (computed on a space grid, with $N_{x}=16$ and $\Delta t=10^{-7}$).

Results are similar to the RTE case: the second-order in time is preserved for big and small values of $\varepsilon$ but not for intermediate regimes.

We are now interested in more qualitative tests by considering the time history of the electric energy $\mathcal{E}\left(t\right)=\sqrt{\int_{0}^{L_{x}}E\left(t,x\right)^{2}\textnormal{d}x}$ in semi-logarithmic scale for different values of $\varepsilon$. We compare the results obtained by MiMa-Part-2 (detailed in Algorithm 5.1) to other schemes: MiMa-part-1, Moment G. and Full PIC (for $\varepsilon$ of order 1) or the scheme for the drift-diffusion model (for small values of $\varepsilon$).

We expect that the number of particles that is necessary to represent in a good way the perturbation $g$ in MiMa-Part methods decreases when $\varepsilon$ diminishes and consider $N_{p}=10^{5}$ if $\varepsilon\geq 0.5$, $N_{p}=10^{4}$ if $\varepsilon=0.1$ and $N_{p}=10^{2}$ if $\varepsilon=10^{-4}$. For comparison, we take the same $N_{p}$ for moment guided and Full PIC methods.

Results for $\varepsilon=10$ are given in Figure 6. For the four particle methods, we take $\Delta t=0.1$, and $N_{x}=128$. With the same parameters, results for $\varepsilon=1$ are presented in Figure 7. For $\varepsilon=0.5$, we consider $\Delta t=0.01$ and $N_{x}=256$ for the four particle methods. Results are given in Figure 8. For these three values of $\varepsilon$, the reference is given by MiMa-Grid with $N_{x}=N_{v}=512$ and $\Delta t=\Delta x^{2}\approx 6\times 10^{-4}$. First, we note that the behaviour of the electric energy is well described during time by micro-macro schemes (MiMa-Part and MiMa-Grid). As observed in [5], the Full PIC method suffers from numerical noise. This is due to the probabilistic character of particle methods (for instance the random initialization of particles). To reduce this noise, we should consider more particles, which would increase the numerical cost. As expected, the moment guided method gives better results than the Full PIC one, but suffers however also from this noise. In MiMa-Part schemes, only the perturbation $g$ is represented by particles (not the whole distribution function $f$), that is why for the same $N_{p}$, the noise is lower, which enables to capture the reference solution for large times. In addition, MiMa-Part-2 is closer to the reference MiMa-Grid than the first-order version MiMa-Part-1.

For smaller values of $\varepsilon$, we compare the four AP schemes (MiMa-Part-1, MiMa-Part-2, MiMa-Grid and Moment G.) to the limit scheme. Results for $\varepsilon=0.1$ are given in Figure 9. Parameters are the following: $\Delta t=10^{-3}$ and $N_{x}=128$ for particle methods, $\Delta t=0.1\Delta x^{2}\approx 3.5\times 10^{-3}$ and $N_{x}=N_{v}=64$ for MiMa-Grid. We observe that MiMa-Part-2 is the best method since it almost coincides with the reference MiMa-Grid method. Finally, results for $\varepsilon=10^{-4}$ are given in Figure 10, where we have $\Delta t=10^{-2}$ and $N_{x}=128$ for particle methods, $\Delta t=0.1\Delta x^{2}\approx 3.5\times 10^{-3}$ and $N_{x}=N_{v}=64$ for MiMa-Grid. The asymptotic regime is well recovered by all these AP methods.

As in [5], we remark that few particles are sufficient in the particle-micro-macro schemes to describe in a good way the solution when $\varepsilon$ is small. The cost is then reduced at the limit.

6.3 Two stream instability

We propose now a study in which the perturbation $g$ is not zero initially and consider the Two-Stream Instability (TSI) test case in both regimes (kinetic and diffusive). The initial distribution function is given by

[TABLE]

with the wave number $k=0.5$ and $\alpha=0.05$. The initial condition for the micro-macro model (1.3) is $\rho\left(t=0,x\right)=1+\alpha\cos\left(kx\right)$ and $g\left(t=0,x,v\right)=\frac{1}{\sqrt{2\pi}}\left(v^{2}-1\right)\textrm{exp}\left(-\frac{v^{2}}{2}\right)\left(1+\alpha\cos\left(kx\right)\right)$. For the limit drift-diffusion equation (1.5), we have as in the Landau damping case $\rho\left(t=0,x\right)=1+\alpha\cos\left(kx\right)$.

We first verify the order in time of the MiMa-Part-2 scheme detailed in Section 5 and plot in Figure 11 the error in $L^{\infty}$ norm of the density $\rho$ at time $T=0.1$ as a function of $\Delta t$ (from $10^{-4}$ to $0.1$) for the following parameters: $N_{x}=16$, $N_{p}=100$. For $\varepsilon=1$, $0.5$ and $0.1$, the reference solution is computed with MiMa-Part-2 using the same parameters but with $\Delta t=10^{-7}$. Whereas for $\varepsilon=10^{-6}$, the reference is a numerical solution of the drift-diffusion equation (computed on a space grid, with $N_{x}=16$ and $\Delta t=10^{-7}$).

As for the RTE and the Landau damping cases, the second-order in time is preserved for big and small values of $\varepsilon$ but not for intermediate regimes.

We are now interested in the time evolution of the electric energy $\mathcal{E}\left(t\right)=\sqrt{\int_{0}^{L_{x}}E\left(t,x\right)^{2}\textnormal{d}x}$ in all regimes.

Results for $\varepsilon=10$ are given in Figure 12. For the four particle methods, we take $N_{p}=10^{6}$, $\Delta t=0.1$, and $N_{x}=128$. By taking $N_{p}=10^{5}$, $\Delta t=0.1$, and $N_{x}=128$ for particle methods, we obtain results for $\varepsilon=1$ presented in Figure 13. For $\varepsilon=0.5$, we consider $N_{p}=10^{5}$, $\Delta t=0.01$ and $N_{x}=256$ for the four particle methods. Results are given in Figure 14. For these three values of $\varepsilon$, the reference is given by MiMa-Grid with $N_{x}=N_{v}=512$ and $\Delta t=\Delta x^{2}\approx 6\times 10^{-4}$. The behaviour of the electric energy is well described during time by micro-macro schemes (MiMa-Part and MiMa-Grid). As previously, the Full PIC method, as well as moment guided method suffer from numerical noise.

To illustrate the efficiency of the method, we plot $f(T=5,x,v)$ obtained by the reference MiMa-Grid, by MiMa-Part-2 and by Full PIC for $\varepsilon=10$ on Figure 15 and for $\varepsilon=0.5$ on Figure 16. For MiMa-Grid and MiMa-Part-2, $f$ is reconstruted from $g$, $\rho$ and $M$, whereas the approximation of $f$ is directly given by the Full PIC scheme. The numerical parameters are the same as previously (see comments on Figures 12 and 14). On Figure 15, we observe that the result obtained by MiMa-Part-2 and Full PIC are in good agreement with MiMa-Grid; however, some numerical noise can be distinguished on $f$ obtained by the Full PIC method. On Figure 16 ($\varepsilon=0.5$), we can see clearly that the level of the noise is higher for Full PIC, which prevents it from giving good results. On the contrary, MiMa-Part-2 produces good results compared to MiMa-Grid, since the noise only affects the micro part $g$, which is small in this regime.

For smaller values of $\varepsilon$, we compare the four AP schemes (MiMa-Part-1, MiMa-Part-2, MiMa-Grid and Moment G.) to the limit scheme. Results for $\varepsilon=10^{-1}$ are given in Figure 17. Parameters are the following: $\Delta t=10^{-3}$ , $N_{x}=128$ and $N_{p}=10^{4}$ for particle methods, $\Delta t=0.1\Delta x^{2}\approx 3.5\times 10^{-3}$ and $N_{x}=N_{v}=64$ for MiMa-Grid. As in the Landau damping case, MiMa-Part-2 gives the best result comparing to the reference MiMa-Grid. Finally, results for $\varepsilon=10^{-4}$ are given in Figure 18, where we have $\Delta t=10^{-2}$, $N_{x}=128$ and $N_{p}=100$ for particle methods, $\Delta t=0.1\Delta x^{2}\approx 3.5\times 10^{-3}$ and $N_{x}=N_{v}=64$ for MiMa-Grid. The asymptotic regime is well recovered by all these AP methods.

7 Conclusion

In this paper, we have presented new micro-macro models for the kinetic radiative transport equation (RTE), as well as for the Vlasov-Poisson-BGK system, in the diffusion scaling with periodic boundary conditions. First-order in time and second-order in time models are derived, and their Lagrangian discretizations are detailed. The obtained schemes are proved to degenerate into implicit discretizations of the limit model (the diffusion equation in the RTE case and the drift-diffusion equation in the Vlasov-Poisson-BGK case) when $\varepsilon\to 0$. This asymptotic property is shown in the numerical results too.

Moreover, thanks to the use of particle methods for the microscopic equation, the numerical cost is reduced when $\varepsilon$ diminishes. Finally, compared to a standard PIC method (where $f$ is represented by particles, and not $g$), the numerical noise is reduced.

In future works, we would like to extend the Monte-Carlo approach proposed for the hydrodynamic limit of Vlasov-BGK in [7] to the diffusion and the drift-diffusion limits.

Appendix A Time discretization for Eulerian schemes

We present the time discretization of (2.4) having in spirit a Eulerian discretization of the phase space. Obviously, the numerical scheme proposed in [2, 8, 21] works well. Now, (3.5) also provides a numerical scheme that we will exploit in this appendix.

Let us consider staggered grids in the phase-space domain and adopt the following notations: $\mbox{x}_{i}=i\Delta x$ and $\mbox{x}_{i+1/2}=i\Delta x+\Delta x/2$, $i\in\mathbb{N}$, define two grids in space and $\mbox{v}_{j}=j\Delta v$, $j\in\mathbb{N}$, defines a grid in velocity, where $\Delta x$ (resp. $\Delta v$) is the step in space (resp. in velocity). Time is also discretized with a time step $\Delta t$ and we note $t^{n}=n\Delta t$, $n\in\mathbb{N}$. The density $\rho$ is discretized on the first space grid: $\rho_{i}^{n}$ approximates $\rho(t^{n},\mbox{x}_{i})$, whereas the perturbation $g$ is discretized on the second one: $g_{i+1/2,j}^{n}$ approximates $g(t^{n},\mbox{x}_{i+1/2},\mbox{v}_{j})$.

Let an approximation $D$ of the spatial derivative, the numerical scheme we propose consists in computing $g_{i+1/2,j}^{n+1}$ with

[TABLE]

where $(\Pi h)_{i+1/2,j}=(\sum_{j}h_{i+1/2,j}\Delta v)$, and then to compute $\rho_{i}^{n+1}$ with

[TABLE]

Proposition A.1**.**

The scheme given by (A)-(A.2) enjoys the AP property, i.e. it satisfies the following properties

• •

for fixed $\varepsilon>0$, the scheme is a first-order (in time) approximation of the original model (1.1),

• •

for fixed $\Delta t>0$, the scheme degenerates into an explicit first-order (in time) scheme of (1.2).

Proof.

We observe easily that when $\varepsilon$ goes to zero, (A) gives

[TABLE]

which, injected in the time discretization (A.2) for $\rho$, gives up to terms of order ${\cal O}(\varepsilon^{2})$

[TABLE]

Since $\sum_{j}v^{2}_{j}\Delta v$ is an approximation of $\int_{-1}^{1}v^{2}\textnormal{d}v=1/3$, we obtain a consistent discretization of the diffusion equation. ∎

Proposition A.1 is of big interest for impliciting the diffusion term $\partial_{xx}\rho$. Indeed, let us rewrite (A) as follows

[TABLE]

with $h_{i+1/2,j}=e^{-\Delta t/\varepsilon^{2}}g^{n+1}_{i+1/2,j}-\varepsilon(1-e^{-\Delta t/\varepsilon^{2}})\left[(I-\Pi)\left(v^{+}_{j}(D^{-}_{x}g^{n})_{i+1/2,j}+v^{-}_{j}(D^{+}_{x}g^{n})_{i+1/2,j}\right)\right].$ Injecting this relation into the macro part, we get

[TABLE]

Since $h_{i+1/2,j}={\cal O}(\varepsilon^{2})$ as $\varepsilon$ goes to zero after two iterations, the asymptotic preserving property is ensured. Moreover, the diffusion term can now be chosen as implicit, so that the macro equation becomes

[TABLE]

and the scheme is now free from the usual diffusion condition on the time step.

The algorithm finally writes

Algorithm A.1**.**

• •

Initialize $g_{i+1/2,j}^{0}$ and $\rho_{i}^{0}$.

At each time step:

• •

Advance micro part with (A).

• •

Advance macro part with (LABEL:macrod_cn).

And we have the following result.

Proposition A.2**.**

The scheme given by (A)-(LABEL:macrod_cn) enjoys the AP property, i.e. it satisfies the following properties

• •

for fixed $\varepsilon>0$, the scheme is a first-order (in time) approximation of the original model (1.1),

• •

for fixed $\Delta t>0$, the scheme degenerates into an implicit first-order (in time) scheme of (1.2).

We do not present here the extension to the Vlasov-Poisson-BGK case, but it is straightforward.

Appendix B Moment guided

In this section, we present the adaptation of the moment guided particle method proposed in [10] to our context. For the sake of simplicity, we present it in the RTE case but note that these computations can also be extended to the Vlasov-Poisson-BGK case, without difficulty.

The kinetic equation on $f$ has to be reformulated to avoid the singularity linked to the transport term. To do that, we proceed as previously, but from (1.1). Indeed, we rewrite equation (1.1) as

[TABLE]

and we integrate between $t^{n}$ and $t^{n+1}$ to get

[TABLE]

We make the following approximation

[TABLE]

where $f^{n}\approx f(t^{n})$ and $\rho^{n}\approx\rho(t^{n})$, $\forall n$.

Making appear the discrete time derivative enables to write

[TABLE]

which we approximate by

[TABLE]

Following the spirit of the moment guided method (see [10]), this equation is coupled with the macro one, that is

[TABLE]

To derive an AP scheme for this latter system satisfied by $(\rho,f)$, we adapt the strategy presented in [10] to our diffusion framework. To do so, we first remark that $\langle vf\rangle=\langle vg\rangle$ and using the expression of $g$ obtained by (3.5), we get the following approximation for the macro flux (considered implicit in time)

[TABLE]

Then, we get the following scheme for $\rho$

[TABLE]

A Lagrangian method can be used to approximate the equation on $f$. As for the micro-macro scheme, we use a splitting procedure

• •

solve $\partial_{t}f+\varepsilon\frac{1-e^{-\Delta t/\varepsilon^{2}}}{\Delta t}v\partial_{x}f=0$

• •

solve $\partial_{t}f=\frac{e^{-\Delta t/\varepsilon^{2}}-1}{\Delta t}f+\frac{1-e^{-\Delta t/\varepsilon^{2}}}{\Delta t}\rho$.

To do that, the transport part is solved with the (non stiff) characteristics

[TABLE]

The source part is solved using the equation satisfied by the weights

[TABLE]

The last step consists in matching the moment of $f^{n+1}$ obtained by the particle method with $\rho^{n+1}$ obtained with (B.2). This can be done using the techniques proposed in [5]. Indeed, considering the function $g=f-\rho$, its weight can be written as

[TABLE]

Then, we apply the discrete version of $(I-\Pi)$ to the weights $\gamma_{k}$ as in [5]

[TABLE]

Acknowledgements

N. Crouseilles and M. Lemou are supported by the French ANR project MOONRISE ANR-14-CE23-0007-01 and by the Enabling Research EUROFusion project CfP-WP14-ER-01/IPP-03. A. Crestetto is supported by the French ANR project ACHYLLES ANR-14-CE25-0001.