# An admissible asymptotic-preserving numerical scheme for the electronic M1 model in the diffusive limit

Sébastien Guisset, Stéphane Brull, Bruno D. Dubroca, Rodolphe Turpault

► To cite this version:

Sébastien Guisset, Stéphane Brull, Bruno D. Dubroca, Rodolphe Turpault. An admissible asymptotic-preserving numerical scheme for the electronic M1 model in the diffusive limit. Communications in Computational Physics, 2018, 24 (5), pp.1326-1354. 10.4208/CICP.OA-2017-0188. hal-03993091

HAL Id: hal-03993091

https://hal.science/hal-03993091

Submitted on 14 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# An admissible asymptotic-preserving numerical scheme for the electronic $M_1$ model in the diffusive limit.

S. Guisset<fnref n="1" />,*, S. Brull<fnref n="2" />, B. Dubroca<fnref n="3" />, R. Turpault<fnref n="2" />

**Abstract.** This work is devoted to the derivation of an admissible asymptotic-preserving scheme for the electronic $M_1$ model in the diffusive regime. A numerical scheme is proposed in order to deal with the mixed derivatives which arise in the diffusive limit leading to an anisotropic diffusion. The derived numerical scheme preserves the re-alisibility domain and enjoys asymptotic-preserving properties correctly handling the diffusive limit recovering the relevant limit equation. In addition, the cases of non constants electric field and collisional parameter are naturally taken into account with the present approach. Numerical test cases validate the considered scheme in the non-collisional and diffusive limits.

## 1 Introduction

In order to initiate nuclear fusion reactions, it was proposed to use laser pulses in order to ignite a deuterium-tritium target. During this process the energy is transported from the critical surface to denser parts through the electron transport. This transport plays a key role in the understanding of plasma physical phenomena such as, parametric [27, 48] and hydrodynamic [20, 55, 62] instabilities, laser-plasma absorption [36, 53], wave damping [18, 41], energy redistribution and hot spot formation [10, 46]. Spitzer and Härm were the first to propose an electron transport theory in a fully ionised plasma without magnetic field. They derived the electron plasma transport coefficients by solving the electron kinetic equation by using the expansion of the electron mean free path over the temperature scale length (denoted $\epsilon$ in this paper). The results of Spitzer and Härm have been reproduced in other works [3, 9, 56] using the early works of Chapman [15, 16] and Enskog [25] for neutral gases. However in the case of non-local regimes [51], the Spitzer-Härm theory is no more valid. Indeed the electron transport plasma coefficients were derived in the case where the isotropic part of the electron distribution function

<footnote n="1">¹ Univ. St-Quentin-en-Yvelines, LMV, UMR 8100, F-78035 Versailles, France.</footnote>

<footnote n="2">² Bordeaux INP, UMR 5251, F-33405 Talence, France.</footnote>

<footnote n="3">³ Univ. Bordeaux, CELIA, UMR 5107, F-33400 Talence, France.</footnote>

<footnote>*Corresponding author. Email addresses: sebastien.guisset@uvsq.fr</footnote>
remains close to the Maxwellian. For example, in the context of inertial confinement fu-
sion, the plasma particles may have an energy distribution far from the thermodynamic
equilibrium so that the fluid description is not adapted. Moreover kinetic effects like the
non local transport [10, 46], wave damping or the development of instabilities [20] can
be important over time scales shorter than the collisional time so that fluid simulations
are insufficient and kinetic codes have to be considered to capture the physical processes.
Therefore, a kinetic description seems unavoidable for the study of inertial confinement
fusion processes. However such a kinetic description is accurate but also computation-
ally expensive for describing most of real physical applications. Kinetic codes are often
limited to time and length much shorter than those studied with fluid simulations. It
is therefore an essential issue to describe kinetic effects by using reduced kinetic codes
operating on fluid time scales.

Angular moments models can be seen as a compromise between kinetic and fluid
models. On the one hand, they have the advantage to be less computationally expensive
than kinetic description since less variables are involved in the models and on the other
hand they provide results with a higher accuracy than fluid models [26, 59]. Grad [31],
initially proposed a moment closure hierarchy which leads to a hyperbolic set of equa-
tions for close equilibrium flows. The hierarchy proposed is based on a polynomial series
expansion of a distribution function close to the Maxwellian equilibrium. However, the
truncation of this expansion leads to a loss of the positivity of the distribution function
and to unrealisable moments, ie moments which can not be derived from a positive dis-
tribution function. In [1, 43, 49, 50, 58], closures based on entropy minimisation principles
are investigated. It has been shown that this closure choice enables to recover funda-
mental properties such as the positivity of the underlying distribution function, the hy-
perbolicity of the model and an entropy dissipation property [32, 43, 47]. In this work,
the moment model is based on an angular moments extraction. The kinetic equation is
only integrated with respect to the velocity direction while the velocity modulus is kept
as a variable. The closure used is based on an entropy minimisation principle and gives
the angular M₁ model. This model is used in numerous applications such as radiative
transfer [5, 17, 22, 52, 60] or electron transport [21, 34, 44]. It satisfies fundamental proper-
ties and recovers the asymptotic diffusion equation in the long time and small mean free
path regimes [23]. In order to perform numerical simulations, the HLL scheme [37] is
often used for the M₁ electronic model because it ensures the positivity of the first angu-
lar moment and the flux limitation property. However, this scheme does not degenerate
accordingly in the diffusive limit as most of the schemes. Even if the scheme is consistent
and one could use a very refined mesh with space step smaller than the mean free path,
such a solution would be far too computationally expensive to be used in practice. Over-
coming this major drawback a class of numerical schemes has emerged over the years
called asymptotic-preserving schemes (AP). Asymptotic-preserving schemes in the sense
of Jin-Levermore [38, 39] are designed to handle multi-scale situations and correctly cap-
ture the asymptotic limit considered. In this context many works have been performed
following different approaches in a one dimensional framework [2, 8, 19, 28, 40, 42, 54]. In
particular, one of the most productive approach from the work of Gosse-Toscani [30] and which has been largely extended [5,7,11,12,14], is based on the modification of approximate Riemann solvers. Some works also deal with the two dimensional case [6, 13]. In the present paper, the difference and the main difficulty comes from the mixed derivatives arising in the diffusive limit. In the present paper, we consider the $M_1$ model for the electronic transport [21, 35, 36, 44, 45]. Ions are supposed fixed and electron-electron collisions are not considered. The angular moment model studied reads

$$ egin{cases} 
\partial_t f_0(t,x,\zeta) + \zeta \partial_x f_1(t,x,\zeta) + E(x) \partial_\zeta f_1(t,x,\zeta) = 0, \\
\partial_t f_1(t,x,\zeta) + \zeta \partial_x f_2(t,x,\zeta) + E(x) \partial_\zeta f_2(t,x,\zeta) \\
\qquad - \frac{E(x)}{\zeta} (f_0(t,x,\zeta) - f_2(t,x,\zeta)) = - \frac{2\alpha^{ei}(x) f_1(t,x,\zeta)}{\zeta^3},
\end{cases} 
\quad (1.1) $$

where $f_0, f_1$ and $f_2$ are the first three angular moments of the electron distribution function $f$. Omitting the $x$ and $t$ dependency, they are given by

$$ egin{aligned} f_0(\zeta) &= \zeta^2 \int_{-1}^{1} f(\mu, \zeta) d\mu, \\ f_1(\zeta) &= \zeta^2 \int_{-1}^{1} f(\mu, \zeta) \mu d\mu, \\ f_2(\zeta) &= \zeta^2 \int_{-1}^{1} f(\mu, \zeta) \mu^2 d\mu. \end{aligned} 
\quad (1.2) $$

The coefficient $\alpha^{ei}$ is a positive physical function which may depend of $x$, $E$ represents the electrostatic field as function of $x$ and $\zeta$ the velocity modulus. The term $\alpha^{ei}/\zeta^3$ comes from the angular integration of the Landau electron-ion collision operator [44]. This parameter is homogeneous to a frequency. The fundamental point of the moments models is the definition of the closure which gives an expression of the highest moment as a function of the lower ones. This closure relation corresponds to an approximation of the underlying distribution function, which the moments system is constructed from. For the $M_1$ model the closure relation originates from an entropy minimisation principle [43,49]. The moment $f_2$ can be computed [21,22] as a function of $f_0$ and $f_1$

$$ f_2(t,x,\zeta) = \chi\left(\frac{f_1(t,x,\zeta)}{f_0(t,x,\zeta)}\right) f_0(t,x,\zeta), \quad \text{with} \quad \chi(\alpha) = (1+\alpha^2+\alpha^4)/3. \quad (1.3) $$

The set of admissible states [21] is defined by

$$ A = \left( (f_0, f_1) \in \mathbb{R}^2, f_0 \ge 0, |f_1| \le f_0 \right). \quad (1.4) $$

In [33], a numerical scheme was proposed for the electronic $M_1$ model without electric field and in the homogeneous case. The scheme derived using the consistency with the
integral form of the approximate Riemann solver ensures the admissibility conditions (1.4) and correctly captures the limit diffusion equation. However, the general model considering the $x$ and $\zeta$ dependencies has not been considered. In such a general case, mixed derivatives arise in the diffusion limit leading to a complex diffusion equation. In addition, the source term $-E(x)(f_0 - f_2)/\zeta$ also contributes in the limit equation. The expression of the limit equation is detailed in section 2.

In this paper, the general electronic $M_1$ model (1.1) is considered. The aim is to propose a numerical scheme, extending the ideas of [33], in order to take into account the mixed derivatives in the diffusive limit. Such a scheme must ensure the admissibility conditions (1.4) and include the contribution of the source term in the diffusion $-E(x)(f_0 - f_2)/\zeta$ limit.

We first introduce the electronic $M_1$ model with its diffusion limit in Section 2. In Section 3, extending the ideas of [33], a numerical scheme is proposed. The scheme is modified to ensure the admissibility conditions (1.4) and to capture the non isotropic diffusion then the asymptotic-preserving property is exhibited. The contribution of the term $-E(x)(f_0 - f_2)/\zeta$ is finally included in the scheme. In Section 4, numerical examples are presented to testify of the efficiency of the method. Finally, Section 5 presents our conclusions.

## 2 Model and diffusion limit

In this section, the diffusion limit of the electronic $M_1$ model (1.1) is introduced. After considering a diffusive scaling, a formal Hilbert expansion is used to derive the limit model. Even if such a procedure is not strictly rigorous from a mathematical point of view, it is known that this approach gives an easy way to derive the limit model. In addition, mathematical rigorous methods used for continuous descriptions can not be easily adapted in this numerical context.

The scaled variables are defined by

$$ \tilde{t} = t/t^*, \quad \tilde{x} = x/x^*, \quad \tilde{\zeta} = \zeta/v_{th}, \quad \tilde{E} = Ex^*/v_{th}^2, $$

where $t^*$ is the typical time of the problem, $x^*$ the typical length and $v_{th}$ the electron thermal velocity. Also, we assume that there exists $\sigma$ such that $\alpha^{ei}(x) = \bar{\alpha}^{ei}\sigma(x)$ where $\bar{\alpha}^{ei}/v_{th}^3$ is the typical electron-ion collisional frequency.

A diffusion scaling is now considered. More precisely the parameters $t^*$ and $x^*$ are chosen such that $\tau^{ei}/t^* = \epsilon^2$, $\lambda^{ei}/x^* = \epsilon$, where the electron-ion collisional period is given by $\tau^{ei} = v_{th}^3/\bar{\alpha}^{ei}$ and the mean free path by $\lambda^{ei} = v_{th}\tau^{ei}$. In this case system (1.1) rewrites

$$ \left\{ \begin{array}{l} \epsilon\partial_t f_0(t,x,\zeta) + \zeta\partial_x f_1(t,x,\zeta) + E(x)\partial_\zeta f_1(t,x,\zeta) = 0, \\ \epsilon\partial_t f_1(t,x,\zeta) + \zeta\partial_x f_2(t,x,\zeta) + E(x)\partial_\zeta f_2(t,x,\zeta) \\ \qquad -\frac{E(x)}{\zeta}(f_0(t,x,\zeta) - f_2(t,x,\zeta)) = -\frac{2\sigma(x)}{\zeta^3}\frac{f_1(t,x,\zeta)}{\epsilon}. \end{array} \right. \qquad (2.1) $$
Introducing the following Hilbert expansion of $f_0^\epsilon$ and $f_1^\epsilon$

$$ \begin{cases} f_0^\epsilon = f_0^0 + \epsilon f_0^1 + O(\epsilon^2), \\ f_1^\epsilon = f_1^0 + \epsilon f_1^1 + O(\epsilon^2), \end{cases} \qquad (2.2) $$

into the second equation of (2.1) taken at order $\epsilon^{-1}$ leads to

$$ f_1^0 = 0. \qquad (2.3) $$

Using the definition of $f_2$ in (1.3), it follows that

$$ f_2^0 = f_0^0 / 3. $$

Inserting the Hilbert expansion (2.2) into the second equation of (1.1) gives at order $\epsilon^0$

$$ f_1^1 = -\frac{\zeta^4}{6\sigma} \partial_x f_0^0 - \frac{E\zeta^3}{6\sigma} \partial_\zeta f_0^0 + \frac{E\zeta^2}{3\sigma} f_0^0. \qquad (2.4) $$

Finally, using the previous equation in the first equation of (1.1) at order $\epsilon^1$, the following limit equation is obtained

$$ \begin{aligned} & \partial_t f_0^0 + \zeta \partial_x \left( -\frac{\zeta^4}{6\sigma} \partial_x f_0^0 - \frac{E\zeta^3}{6\sigma} \partial_\zeta f_0^0 + \frac{E\zeta^2}{3\sigma} f_0^0 \right) \\ & + E \partial_\zeta \left( -\frac{\zeta^4}{6\sigma} \partial_x f_0^0 - \frac{E\zeta^3}{6\sigma} \partial_\zeta f_0^0 + \frac{E\zeta^2}{3\sigma} f_0^0 \right) = 0. \end{aligned} \qquad (2.5) $$

In the case $E=0$, one recognises a classical diffusion equation involving a second order space derivative with a diffusion coefficient of $-\zeta^5/6\sigma$. However, in the general case this limit equation involves mixed $x$ and $\zeta$ derivatives leading to an anisotropic diffusion. In addition, the source term $E(f_0-f_2)/\zeta$ also contributes in the diffusive limit leading to the term $(E\zeta^2/(3\sigma))f_0^0$ in the right side of (2.4) and in the $x$ and $\zeta$ derivatives of (2.5). Such an asymptotic limit is unusual compared to what has been studied in radiative transfer for example [4, 5]. The difference lies in the fact that here charged particles are considered. Then, the contribution of the electric field must be taken into account leading to these unexpected limit involving mixed derivatives.

# 3 Numerical scheme

## 3.1 Previous results and limit of the approach

In [33] a numerical scheme was proposed for the simplified case without electric field $E$ and for the homogeneous case with non-constant electric field. The scheme is based on an approximate Riemann solver whose the intermediate states are chosen to obtain the asymptotic-preserving property while preserving the admissibility of the solution. However the generalisation of such approach to the general model (1.1) is challenging since mixed-derivatives in space and velocity modulus arise in the asymptotic limit (2.5). In this section we briefly summarise the results obtained in [33] and explain the difficulty encounters when generalising to the general model (1.1).
### 3.1.1 Previous results

We start recalling the results obtained in [33]. For conciseness we only recall the case with no electric field which was studied in the first part of [33]. In this case the model writes

$$ \begin{cases} \partial_t f_0(t,x,\zeta) + \zeta \partial_x f_1(t,x,\zeta) = 0, \\ \partial_t f_1(t,x,\zeta) + \zeta \partial_x f_2(t,x,\zeta) = -\frac{2\alpha^{ei}(x)f_1(t,x,\zeta)}{\zeta^3}, \end{cases} \qquad (3.1) $$

and the corresponding diffusive limit reads

$$ \partial_t f_0^0 + \zeta \partial_x \left(-\frac{\zeta^4}{6\sigma}\partial_x f_0^0\right) = 0. \qquad (3.2) $$

We consider an uniform mesh with a constant space step $\Delta x = x_{i+1/2} - x_{i-1/2}$ and a time step $\Delta t$. For this case the following numerical scheme has been proposed to compute the numerical solution at time $t^{n+1}$

$$ \left\{ \begin{aligned} f_{0i}^{n+1} &= \frac{a_x \Delta t}{\Delta x} f_{0i-1/2}^{R*} + \left(1 - \frac{2a_x \Delta t}{\Delta x}\right) f_{0i}^{n} + \frac{a_x \Delta t}{\Delta x} f_{0i+1/2}^{L*}, \\ f_{1i}^{n+1} &= \frac{a_x \Delta t}{\Delta x} f_{1i-1/2}^{*} + \left(1 - \frac{2a_x \Delta t}{\Delta x}\right) f_{1i}^{n} + \frac{a_x \Delta t}{\Delta x} f_{1i+1/2}^{*}, \end{aligned} \right. \qquad (3.3) $$

where the intermediate states are given by

$$ f_1^* = \frac{2a_x \zeta^3}{2a_x \zeta^3 + 2\alpha^{ei}\Delta x} \left( \frac{f_1^L + f_1^R}{2} - \frac{1}{2a_x} (\zeta f_2^R - \zeta f_2^L) \right), $$

and

$$ \begin{cases} f_0^{L*} = \tilde{f}_0 - \Gamma\theta, \\ f_0^{R*} = \tilde{f}_0 + \Gamma\theta, \end{cases} $$

with

$$ \tilde{f}_0 = \frac{f_0^L + f_0^R}{2} - \frac{1}{2a_x}(\zeta f_1^R - \zeta f_1^L), $$

and

$$ \Gamma = \frac{1}{2} \left( f_0^R - f_0^L - \frac{\zeta}{a_x} (f_1^L - 2f_1^* + f_1^R) \right). $$

The coefficient $\theta \in [0,1]$ is fixed to ensure the admissibility conditions and is chosen as $\theta = \min(\tilde{\theta}, 1)$ where

$$ \tilde{\theta} = \frac{\tilde{f}_0 - |f_1^*|}{|\Gamma|} \geq 0. $$

We associate the standard CFL condition

$$ \Delta t \leq \Delta x / a_x. \qquad (3.4) $$

The properties obtained in [33] can be summarised as the following
**Theorem 3.1.** *The numerical scheme (3.3)*

(i) is consistent with (3.1),

(ii) preserves the admissibility of the numerical solutions under the CFL condition (3.4),

(iii) considered in diffusive regimes, is consistent with the limit diffusion equation (3.2).

The spatially homogeneous case but taking into account the electric field has been investigated in the second part of [33].

**Remark:** The scheme (3.3) is explicit in the sense that, in order to update the numerical solution, there is no need to solve an equation or a system. More precisely, in the same spirit as in [7] (see also [29]), the approximate Riemann solver considered here is explicit and, this is the key point, the intermediate states carefully combine the hyperbolic part and the source term. To go further, an implicit treatment of the source term can been seen "inside" the definition of the intermediate states. Indeed, in the derivation of the intermediate states $f_1^*$ in [32], right after the time integration of the source term (writing the consistency relations in the integral sense) a time implicit approximation was chosen. This implication makes appear the term $2a_x\zeta^3/(2a_x\zeta^3+2\alpha^{ei}\Delta x)$ in the definition of $f_1^*$.

### 3.1.2 Limit of the approach

In this part we show that the generalisation of this approach to the general model (1.1) is not straightforward. Indeed, if one simply sums the schemes proposed in [33], for the simplified case without electric field $E$ and for the homogeneous case with non-constant electric field, a direct discrete Hilbert expansion shows that the wrong limit model is obtained since the mixed derivatives are not recovered in the limit regime. More precisely, only the second order space and velocity modulus derivatives are captured by the scheme and the anisotropic diffusion is missing. This point is proved in this section.

We first consider the case without the source term $\frac{E}{\zeta}(f_0 - f_2)$. The general model is studied in part 3.4. For clarity, we start without considering it. In this case the electronic $M_1$ model reads

$$ \left\{ \begin{aligned} \partial_t f_0(t, x, \zeta) + \zeta \partial_x f_1(t, x, \zeta) + E(x) \partial_\zeta f_1(t, x, \zeta) &= 0, \\ \partial_t f_1(t, x, \zeta) + \zeta \partial_x f_2(t, x, \zeta) + E(x) \partial_\zeta f_2(t, x, \zeta) &= -\frac{2\alpha^{ei}(x) f_1(t, x, \zeta)}{\zeta^3}, \end{aligned} \right. \qquad (3.5) $$

and its diffusive limit equation writes

$$ \partial_t f_0^0 + \zeta \partial_x \left( -\frac{\zeta^4}{6\sigma} \partial_x f_0^0 - \frac{E\zeta^3}{6\sigma} \partial_\zeta f_0^0 \right) + E \partial_\zeta \left( -\frac{\zeta^4}{6\sigma} \partial_x f_0^0 - \frac{E\zeta^3}{6\sigma} \partial_\zeta f_0^0 \right) = 0. \qquad (3.6) $$

Let us consider a constant velocity modulus step $\Delta\zeta = \zeta_{i+1/2} - \zeta_{i-1/2}$. From now the indices $i$ and $j$ refer to the numerical solution considered in $x_i$ and $\zeta_j$. Specular reflection
<figure><img src="image_2.png" /><figcaption>Figure 1: Structure of the approximate Riemann solver considered.</figcaption></figure>

boundary conditions are considered in $\zeta=0$. At the discrete level the following conditions are imposed in the ghost cell denoted by the index 0

$$f_{0i0}^n = f_{0i1}^n, \quad f_{1i0}^n = -f_{1i1}^n \quad \text{for all } i.$$

By summing the schemes proposed in [33], for the simplified case without electric field $E$ and for the homogeneous case with non-constant electric field ones obtains the following scheme for (1.1)

$$\frac{U_{ij}^{n+1} - U_{ij}^{n}}{\Delta t} = \frac{a_x}{\Delta x} U_{i-1/2j}^{R*} + \frac{2a_x}{\Delta x} U_{ij}^{n} + \frac{a_x}{\Delta x} U_{i+1/2j}^{L*} \\ + \frac{a_\zeta}{\Delta \zeta} U_{ij-1/2}^{R*} + \frac{2a_\zeta}{\Delta \zeta} U_{ij}^{n} + \frac{a_\zeta}{\Delta \zeta} U_{ij+1/2}^{L*}, \quad (3.7)$$

where the intermediate states of the approximated Riemann solver (see Figure 1) $U_{i+1/2j}^{L*}$, $U_{i-1/2j}^{R*}$, $U_{ij+1/2}^{L*}$ and $U_{ij-1/2}^{R*}$ are defined by

$$U_{i-1/2j}^{R*} = \begin{pmatrix} f_{0i-1/2j}^{R*} \\ f_{1i-1/2j}^{*} \end{pmatrix}, \quad U_{i+1/2j}^{L*} = \begin{pmatrix} f_{0i+1/2j}^{L*} \\ f_{1i+1/2j}^{*} \end{pmatrix},$$

$$U_{ij-1/2}^{R*} = \begin{pmatrix} f_{0ij-1/2}^{R*} \\ f_{1ij-1/2}^{*} \end{pmatrix}, \quad U_{ij+1/2}^{L*} = \begin{pmatrix} f_{0ij+1/2}^{L*} \\ f_{1ij+1/2}^{*} \end{pmatrix}.$$

Following [32] the second components of the intermediate states at each interface are chosen equal, ie $f_{1i+1/2j}^{L*} = f_{1i+1/2j}^{R*} = f_{1i+1/2j}^{*}$ and $f_{1ij+1/2}^{L*} = f_{1ij+1/2}^{R*} = f_{1ij+1/2}^{*}$. Following [4, 5, 33], the velocity waves $a_x$ and $a_\zeta$ are fixed such that

$$a_x = \zeta_j, \quad a_\zeta = |E_i|. \quad (3.8)$$
For clarity, in the following, we omit the dependency of the speed $a_x$ in velocity modulus and $a_ζ$ in space. However, the results presented hold in the general case. Following the ideas introduced in [33] one can naturally proposes the following intermediate states

$$f_{1i+1/2j}^* = \alpha_{i+1/2j} \left( \frac{f_{1i+1j}^n + f_{1ij}^n}{2} - \frac{(\zeta_j f_{2i+1j}^n - \zeta_j f_{2ij}^n)}{2a_x} \right), \quad (3.9)$$

$$f_{1ij+1/2}^* = \beta_{ij+1/2} \left( \frac{f_{1ij+1}^n + f_{1ij}^n}{2} - \frac{(E_i f_{2ij+1}^n - E_i f_{2ij}^n)}{2a_ζ} \right), \quad (3.10)$$

with

$$\alpha_{i+1/2j} = \frac{2a_x \zeta_j^3}{2a_x \zeta_j^3 + \alpha_{i+1/2}^{ei} \Delta x}, \quad \beta_{ij+1/2} = \frac{2a_ζ \zeta_{j+1/2}^3}{2a_ζ \zeta_{j+1/2}^3 + \alpha_i^{ei} \Delta ζ}.$$

We introduce the following notations

$$\begin{aligned} \tilde{f}_{0i+1/2j} &= \frac{f_{0i+1j}^n + f_{0ij}^n}{2} - \zeta_j \frac{(f_{1i+1j}^n - f_{1ij}^n)}{2a_x}, \\ \tilde{f}_{0ij+1/2} &= \frac{f_{0ij+1}^n + f_{0ij}^n}{2} - E_i \frac{(f_{1ij+1}^n - f_{1ij}^n)}{2a_ζ}. \end{aligned} \quad (3.11)$$

and

$$\begin{aligned} \tilde{f}_{1i+1/2j} &= \frac{f_{1i+1j}^n + f_{1ij}^n}{2} - \zeta_j \frac{(f_{2i+1j}^n - f_{2ij}^n)}{2a_x}, \\ \tilde{f}_{1ij+1/2} &= \frac{f_{1ij+1}^n + f_{1ij}^n}{2} - E_i \frac{(f_{2ij+1}^n - f_{2ij}^n)}{2a_ζ}. \end{aligned} \quad (3.12)$$

In [33], the intermediate states of the considered approximate Riemann solvers were defined using consistency relations and a corrective coefficient to ensure the admissibility conditions. Following these ideas the intermediate states $f_{0i+1/2j}^{R*}$ and $f_{0i+1/2j}^{L*}$ are defined by

$$\begin{cases} f_{0i+1/2j}^{L*} = \tilde{f}_{0i+1/2j} - \Gamma_{i+1/2j} \theta_{1i+1/2j}, \\ f_{0i+1/2j}^{R*} = \tilde{f}_{0i+1/2j} + \Gamma_{i+1/2j} \theta_{1i+1/2j}, \end{cases} \quad (3.13)$$

with

$$\Gamma_{i+1/2j} = \frac{1}{2} \left( f_{0i+1j}^n - f_{0ij}^n - \frac{\zeta_j}{a_x} (f_{1ij}^n - 2f_{1i+1/2j}^* + f_{1i+1j}^n) \right), \quad (3.14)$$

and the coefficient $\theta_{1i+1/2j}$ is fixed in order to ensure the admissibility conditions (1.4). Similarly, the definitions of $f_{0ij+1/2}^{R*}$ and $f_{0ij-1/2}^{L*}$ read

$$\begin{cases} f_{0ij+1/2}^{L*} = \tilde{f}_{0ij+1/2} - \Gamma_{ij+1/2} \theta_{2ij+1/2}, \\ f_{0ij+1/2}^{R*} = \tilde{f}_{0ij+1/2} + \Gamma_{ij+1/2} \theta_{2ij+1/2}, \end{cases} \quad (3.15)$$
with

$$ \Gamma_{ij+1/2} = \frac{1}{2} \left( f_{0ij+1}^n - f_{0ij}^n - \frac{\zeta_j}{a_\zeta} (f_{1ij}^n - 2f_{1ij+1/2}^* + f_{1ij+1}^n) \right). \quad (3.16) $$

Then the corrective coefficients $\theta_{1i+1/2j}$ and $\theta_{2ij+1/2}$ are fixed in the interval $[0,1]$, the larger possible such that the admissibility requirements (1.4) are fulfilled. A simple calculation gives the following conditions

$$ \tilde{\theta}_{1i+1/2j} = \frac{\tilde{f}_{0i+1/2j} - \alpha_{i+1/2j} |\tilde{f}_{1i+1/2j}|}{|\Gamma_{i+1/2j}|}, \quad (3.17) $$

and

$$ \tilde{\theta}_{2ij+1/2} = \frac{\tilde{f}_{0ij+1/2} - \beta_{ij+1/2} |\tilde{f}_{1ij+1/2}|}{|\Gamma_{ij+1/2}|}. \quad (3.18) $$

Finally, $\theta_{1i+1/2j} = \min(\tilde{\theta}_{1i+1/2j}, 1)$ and $\theta_{2ij+1/2} = \min(\tilde{\theta}_{2ij+1/2}, 1)$.

However as shown with the following result the wrong asymptotic behaviour is obtained with this scheme.

**Theorem 3.2. (Wrong consistency in the diffusion regime)**
In the diffusion limit, the numerical scheme (3.7) degenerates into

$$ \frac{f_{0ij}^{n+1,0} - f_{0ij}^{n,0}}{\Delta t} = \frac{\zeta_j}{\Delta x} \left( \frac{\zeta_j^4}{6\sigma_{i+1/2}\Delta x} (f_{0i+1j}^{n,0} - f_{0ij}^{n,0}) - \frac{\zeta_j^4}{6\sigma_{i-1/2}\Delta x} (f_{0i1j}^{n,0} - f_{0i-1j}^{n,0}) \right) + \frac{E_i}{\Delta \zeta} \left( \frac{E_i \zeta_{j+1/2}^3}{6\sigma_i \Delta \zeta} (f_{0ij+1}^{n,0} - f_{0ij}^{n,0}) - \frac{E_i \zeta_{j-1/2}^3}{6\sigma_i \Delta \zeta} (f_{0i1j}^{n,0} - f_{0ij-1}^{n,0}) \right), \quad (3.19) $$

which is not consistent with the limit diffusion model (3.6).

*Proof.* Following the same approach as in [5,7,33], using the diffusion scaling and equation (3.7) leads to

$$ \begin{aligned} \frac{\epsilon U_{ij}^{n+1,\epsilon} - U_{ij}^{n,\epsilon}}{\Delta t} = & \frac{a_x}{\Delta x} U_{i-1/2j}^{R*,\epsilon} - \frac{2a_x}{\Delta x} U_{ij}^{n,\epsilon} + \frac{a_x}{\Delta x} U_{i+1/2j}^{L*,\epsilon} \\ & + \frac{a_\zeta}{\Delta \zeta} U_{ij-1/2}^{R*,\epsilon} - \frac{2a_\zeta}{\Delta \zeta} U_{ij}^{n,\epsilon} + \frac{a_\zeta}{\Delta \zeta} U_{ij+1/2}^{L*,\epsilon} \end{aligned} \quad (3.20) $$

and equations (3.9) and (3.10) give

$$ \begin{aligned} f_{1i+1/2j}^{*,\epsilon} &= \alpha_{i+1/2j}^{\epsilon} \left( \frac{f_{1i+1j}^{n,\epsilon} + f_{1ij}^{n,\epsilon}}{2} - \frac{(\zeta_j f_{2i+1j}^{n,\epsilon} - \zeta_j f_{2ij}^{n,\epsilon})}{2a_x} \right), \quad (3.21) \\ f_{1ij+1/2}^{*,\epsilon} &= \beta_{ij+1/2}^{\epsilon} \left( \frac{f_{1ij+1}^{n,\epsilon} + f_{1ij}^{n,\epsilon}}{2} - \frac{(E_i f_{2ij+1}^{n,\epsilon} - E_i f_{2ij}^{n,\epsilon})}{2a_\zeta} \right), \end{aligned} $$
with

$$ \alpha_{i+1/2j}^{\epsilon} = \frac{2a_x \zeta_j^3}{2a_x \zeta_j^3 + \sigma_{i+1/2} \Delta x / \epsilon}, \qquad \beta_{ij+1/2}^{\epsilon} = \frac{2a_\zeta \zeta_{j+1/2}^3}{2a_\zeta \zeta_{j+1/2}^3 + \sigma_i \Delta \zeta / \epsilon}. \qquad (3.22) $$

Consequently as $\epsilon$ tends to zero, (3.22) shows that $\alpha^\epsilon$ and $\beta^\epsilon$ also tends to zero. Therefore from (3.21) one obtains that

$$ f_{1i+1/2j}^{*,0} = 0 \quad \text{and} \quad f_{1ij+1/2}^{*,0} = 0. \qquad (3.23) $$

Now, since the second component of (3.20) reads

$$ \epsilon \frac{f_{1ij}^{n+1,\epsilon} - f_{1ij}^{n,\epsilon}}{\Delta t} = \frac{a_x}{\Delta x} f_{1i-1/2j}^{*,\epsilon} - \frac{2a_x}{\Delta x} f_{1ij}^{n,\epsilon} + \frac{a_x}{\Delta x} f_{1i+1/2j}^{*,\epsilon} \\ + \frac{a_\zeta}{\Delta \zeta} f_{1ij-1/2}^{*,\epsilon} - \frac{2a_\zeta}{\Delta \zeta} f_{1ij}^{n,\epsilon} + \frac{a_\zeta}{\Delta \zeta} f_{1ij+1/2}^{*,\epsilon} $$

it follows that at order $\epsilon^0$ the previous equation gives

$$ f_{1ij}^{n,0} = 0. \qquad (3.24) $$

Therefore one correctly recovers the condition (2.3) and from the definition (1.3), it follows that

$$ f_{2ij}^{n,0} = f_{0ij}^{n,0} / 3. $$

The consistency with the wrong limit diffusion equation is now proved. Firstly, one remarks that in the limit $\epsilon$ tends to zero, the definitions (3.17) and (3.18) and the result (3.24) give

$$ \theta_{1i+1/2j} = 1, \qquad \theta_{2ij+1/2} = 1. \qquad (3.25) $$

Indeed, when $\epsilon$ tends to zero, the definitions (3.17) and (3.18) lead to

$$ \tilde{\theta}_{1i+1/2j} = \frac{f_{0i+1j}^{n,0} + f_{0ij}^{n,0}}{|f_{0i+1j}^{n,0} - f_{0ij}^{n,0}|} \ge 1, \qquad \tilde{\theta}_{2ij+1/2} = \frac{f_{0ij+1}^{n,0} + f_{0ij}^{n,0}}{|f_{0ij+1}^{n,0} - f_{0ij}^{n,0}|} \ge 1. $$

Then from the definitions (3.13)-(3.14)-(3.15)-(3.16), we obtain

$$ \begin{cases} f_{0i+1/2j}^{L*} = f_{0ij}^{n} + (f_{1i+1/2j}^{L*} - f_{1ij}^{n}), \\ f_{0i-1/2j}^{R*} = f_{0ij}^{n} + (f_{1ij}^{n} - f_{1i-1/2j}^{R*}), \end{cases} \qquad (3.26) $$

and

$$ \begin{cases} f_{0ij+1/2}^{L*} = f_{0ij}^{n} + (f_{1ij+1/2}^{L*} - f_{1ij}^{n}), \\ f_{0ij-1/2}^{R*} = f_{0ij}^{n} + (f_{1ij}^{n} - f_{1ij-1/2}^{R*}), \end{cases} \qquad (3.27) $$
Secondly, the first component of (3.20) reads

$$\epsilon \frac{f_{0ij}^{n+1,\epsilon} - f_{0ij}^{n,\epsilon}}{\Delta t} = \frac{a_x}{\Delta x} f_{0i-1/2j}^{R*,\epsilon} - \frac{2a_x}{\Delta x} f_{0ij}^{n,\epsilon} + \frac{a_x}{\Delta x} f_{0i+1/2j}^{L*,\epsilon} \\ + \frac{a_\zeta}{\Delta \zeta} f_{0ij-1/2}^{R*,\epsilon} - \frac{2a_\zeta}{\Delta \zeta} f_{0ij}^{n,\epsilon} + \frac{a_\zeta}{\Delta \zeta} f_{0ij+1/2}^{L*,\epsilon}$$

and by using (3.26)-(3.27) it follows that

$$\epsilon \frac{f_{0ij}^{n+1,\epsilon} - f_{0ij}^{n,\epsilon}}{\Delta t} = a_x \frac{f_{1i+1/2j}^{L*,\epsilon} - f_{1i-1/2j}^{R*,\epsilon}}{\Delta x} + a_\zeta \frac{f_{1ij+1/2}^{L*,\epsilon} - f_{1ij-1/2}^{R*,\epsilon}}{\Delta \zeta}$$

Finally, the previous equation considered at order $\epsilon^0$ brings no information but at the order $\epsilon^1$ one obtains

$$\frac{f_{0ij}^{n+1,0} - f_{0ij}^{n,0}}{\Delta t} = a_x \frac{f_{1i+1/2j}^{L*,1} - f_{1i-1/2j}^{R*,1}}{\Delta x} + a_\zeta \frac{f_{1ij+1/2}^{L*,1} - f_{1ij-1/2}^{R*,1}}{\Delta \zeta},$$

which gives, by using the definitions (3.21), the numerical scheme (3.19). □

This wrong behaviour was expected since standard asymptotic-preserving corrections do not give anisotropic numerical viscosity. Therefore an innovative approach must be investigated. In the present work, in order to capture the complete anisotropic model while ensuring the admissibility requirement, new intermediate states are considered.

## 3.2 Derivation of the new scheme

In this part the derivation of a new numerical scheme for the model (3.5) is detailed based on the scheme introduced in section 3.1. The objective is to extend the scheme introduced in [33] to the model (1.1) in the sense that one should recover the scheme of [33] in the homogeneous case or in the case without electric field.

We propose to modify the scheme introduced in section 3.1 to capture the correct asymptotic limit. Therefore the following states are proposed

$$f_{1i+1/2j}^* = \alpha_{i+1/2j} \left( \frac{f_{1i+1j} + f_{1ij}}{2} - \frac{(\zeta_j f_{2i+1j} - \zeta_j f_{2ij})}{2a_x} \right) \quad (3.28) \\ - c_{i+1/2j} \left( \frac{\partial f_0}{\partial \zeta} \right)_{i+1/2j} (1 - \alpha_{i+1/2j}),$$

$$f_{1ij+1/2}^* = \beta_{ij+1/2} \left( \frac{f_{1ij+1} + f_{1ij}}{2} - \frac{(E_i f_{2ij+1} - E_i f_{2ij})}{2a_\zeta} \right) \quad (3.29) \\ - \bar{c}_{ij+1/2} \left( \frac{\partial f_0}{\partial x} \right)_{ij+1/2} (1 - \beta_{ij+1/2}),$$
with

$$ \alpha_{i+1/2j} = \frac{2a_x \zeta_j^3}{2a_x \zeta_j^3 + \alpha_{i+1/2}^{ei} \Delta x'}, \qquad \beta_{ij+1/2} = \frac{2a_\zeta \zeta_{j+1/2}^3}{2a_\zeta \zeta_{j+1/2}^3 + \alpha_i^{ei} \Delta \zeta}. \qquad (3.30) $$

The differences with the intermediate states of [33] come from the last terms in (3.28) and (3.29). In the present study it is proved that these intermediate states give the asymptotic-preserving property for the general model (3.6). In the present case the numerical viscosity contributes in the x and $\zeta$ directions. The terms $(\frac{\partial f_0}{\partial \zeta})_{i+1/2j}$, $(\frac{\partial f_0}{\partial x})_{ij+1/2}$ and the coefficients $c$ and $\bar{c}$ are fixed in order to obtain the relevant limit equation (2.5) in the diffusion regime. We set

$$ c_{i+1/2j} = \frac{E_{i+1/2} \Delta x}{3a_x}, \qquad \bar{c}_{ij+1/2} = \frac{\zeta_{j+1/2} \Delta \zeta}{3a_\zeta}. \qquad (3.31) $$

We use an upwind scheme for the discretisation of the terms $(\frac{\partial f_0}{\partial \zeta})_{i+1/2j}$ and $(\frac{\partial f_0}{\partial x})_{ij+1/2}$. Since the coefficient $\bar{c}$ is always positive the term $(\frac{\partial f_0}{\partial x})_{ij+1/2}$ is always upwinded in the same direction

$$ \bar{c}_{ij+1/2} (\frac{\partial f_0}{\partial x})_{ij+1/2} \approx \bar{c}_{ij+1/2} \frac{f_{0i+1j+1/2} - f_{0ij+1/2}}{\Delta x} \\ \approx \bar{c}_{ij+1/2} \frac{f_{0i+1j+1} - f_{0ij+1} + f_{0i+1j} - f_{0ij}}{2\Delta x}. $$

Similarly one obtains

$$ c_{i+1/2j} (\frac{\partial f_0}{\partial \zeta})_{i+1/2j} \approx \begin{cases} c_{i+1/2j} \frac{f_{0i+1j} - f_{0i+1j-1} + f_{0ij} - f_{0ij-1}}{2\Delta\zeta} & \text{if } c_{i+1/2j} < 0, \\ c_{i+1/2j} \frac{f_{0i+1j+1} - f_{0i+1j} + f_{0ij+1} - f_{0ij}}{2\Delta\zeta} & \text{if } c_{i+1/2j} > 0. \end{cases} $$

The previous two conditions rewrite

$$ c_{i+1/2j} (\frac{\partial f_0}{\partial \zeta})_{i+1/2j} = c_{i+1/2j}^{-} \frac{f_{0i+1j} - f_{0i+1j-1} + f_{0ij} - f_{0ij-1}}{2\Delta\zeta} + c_{i+1/2j}^{+} \frac{f_{0i+1j+1} - f_{0i+1j} + f_{0ij+1} - f_{0ij}}{2\Delta\zeta}, $$

with $(c)^+ = \max(c, 0)$ and $(c)^- = \min(c, 0)$.

In order to ensure the admissibility conditions (1.4), the definitions of the intermediate states $f_{1i+1/2j}^*$ and $f_{1ij+1/2}^*$ given in (3.28) and (3.29) are modified such that

$$ f_{1i+1/2j}^* = \alpha_{i+1/2j} (\tilde{f}_{1i+1/2j} - \theta_{1i+1/2j} c_{i+1/2j} (\frac{\partial f_0}{\partial \zeta})_{i+1/2j} (1 - \alpha_{i+1/2j})), \qquad (3.32) $$

$$ f_{1ij+1/2}^* = \beta_{ij+1/2} (\tilde{f}_{1ij+1/2} - \theta_{2ij+1/2} \bar{c}_{ij+1/2} (\frac{\partial f_0}{\partial x})_{ij+1/2} (1 - \beta_{ij+1/2})). \qquad (3.33) $$
**Remark 3.1.** In the case $\theta_{1i+1/2j}=0$ and $\theta_{2ij+1/2}=0$, the admissibility requirements (1.4) are fulfilled.

Then $\theta_{1i+1/2j}$ and $\theta_{2ij+1/2}$ are now fixed in the interval $[0,1]$, the larger possible such that the admissibility requirements (1.4) are fulfilled. A simple calculation gives the following conditions

$$ \tilde{\theta}_{1i+1/2j} = \frac{\tilde{f}_{0i+1/2j} - \alpha_{i+1/2j} |\tilde{f}_{1i+1/2j}|}{|\Gamma_{i+1/2j}| + |\alpha_{i+1/2j} (\frac{\partial f_0}{\partial \zeta})_{i+1/2j} c_{i+1/2j}|}, \quad (3.34) $$

and

$$ \tilde{\theta}_{2ij+1/2} = \frac{\tilde{f}_{0ij+1/2} - \beta_{ij+1/2} |\tilde{f}_{1ij+1/2}|}{|\Gamma_{ij+1/2}| + |\beta_{ij+1/2} (\frac{\partial f_0}{\partial \zeta})_{ij+1/2} \tilde{c}_{ij+1/2}|}. \quad (3.35) $$

Finally, $\theta_{1i+1/2j} = \min(\tilde{\theta}_{1i+1/2j}, 1)$ and $\theta_{2ij+1/2} = \min(\tilde{\theta}_{2ij+1/2}, 1)$.

We point out that the definitions of the intermediates states $f_0^{L,*}$, $f_0^{R,*}$ are left unchanged and are given by (3.13)-(3.15).

### 3.2.1 Properties of the new scheme

In this part, the admissibility, the consistency in the classical regime and the asymptotic-preserving property in the diffusive regime of the scheme are proved.

**Theorem 3.3.** (Admissibility) If for all $(i,j) \in \mathbb{N}^2$, $U_{i,j}^n \in \mathcal{A}$, then for all $(i,j) \in \mathbb{N}^2$, $U_{i,j}^{n+1} \in \mathcal{A}$ as soon as the following CFL condition holds

$$ \Delta t \le \frac{\Delta \zeta \Delta x}{(2a_x \Delta \zeta + 2a_\zeta \Delta x)}. \quad (3.36) $$

*Proof.* The numerical scheme (3.7) also writes as a convex combination of vectors of $\mathcal{A}$

$$ U_{ij}^{n+1} = (1 - \frac{2a_x \Delta t}{\Delta x} - \frac{2a_\zeta \Delta t}{\Delta \zeta}) U_{ij}^n + \frac{a_x \Delta t}{\Delta x} U_{i-1/2j}^{R*} + \frac{a_\zeta \Delta t}{\Delta \zeta} U_{i+1/2j}^{L*} \\ + \frac{a_\zeta \Delta t}{\Delta \zeta} U_{i-1/2j}^{R*} + \frac{a_\zeta \Delta t}{\Delta \zeta} U_{i+1/2j}^{L*}, \quad (3.37) $$

Using the definitions of $\theta_1$ and $\theta_2$ given in (3.34) and (3.35) the intermediate states $U_{i-1/2j}^{R*}$, $U_{i+1/2j}^{L*}$, $U_{ij-1/2}^{R*}$ and $U_{ij+1/2}^{L*}$ belong to $\mathcal{A}$. Since $\mathcal{A}$ is a convex space it follows that the updated states $U_i^{n+1}$ belongs to $\mathcal{A}$. $\square$

**Theorem 3.4.** (Consistency in the classical regime)

The numerical scheme (3.7) is consistent, when $\Delta t$ and $\Delta x$ tend to zero, with the set of equations (3.5).
*Proof.* Using the definitions (3.28) and (3.29), the second component of (3.7) reads

$$ \frac{f_{1ij}^{n+1} - f_{1ij}^{n}}{\Delta t} = \frac{a_x}{\Delta x} \left[ \alpha_{i+1/2j} \left( \tilde{f}_{1i+1/2j} - \theta_{1i+1/2j} c_{i+1/2j} \left( \frac{\partial f_0}{\partial \zeta} \right)_{i+1/2j} (1-\alpha_{i+1/2j}) \right) \right] \\ - \frac{2a_x}{\Delta x} f_{1ij}^n \\ + \frac{a_x}{\Delta x} \left[ \alpha_{i-1/2j} \left( \tilde{f}_{1i-1/2j} - \theta_{1i-1/2j} c_{i-1/2j} \left( \frac{\partial f_0}{\partial \zeta} \right)_{i-1/2j} (1-\alpha_{i-1/2j}) \right) \right] \quad (3.38) \\ + \frac{a_\zeta}{\Delta \zeta} \left[ \beta_{ij+1/2} \left( \tilde{f}_{1ij+1/2} - \theta_{2ij+1/2} \bar{c}_{ij+1/2} \left( \frac{\partial f_0}{\partial x} \right)_{ij+1/2} (1-\beta_{ij+1/2}) \right) \right] \\ - \frac{2a_\zeta}{\Delta \zeta} f_{1ij}^n \\ + \frac{a_\zeta}{\Delta \zeta} \left[ \beta_{ij-1/2} \left( \tilde{f}_{1ij-1/2} - \theta_{2ij-1/2} \bar{c}_{ij-1/2} \left( \frac{\partial f_0}{\partial x} \right)_{ij-1/2} (1-\beta_{ij-1/2}) \right) \right]. $$

Inserting the definitions (3.11) into (3.38) and using the following expressions for $\alpha_{i+1/2j}$ and $\beta_{ij+1/2}$

$$ \alpha_{i+1/2j} = \frac{2a_x \zeta_j^3}{2a_x \zeta_j^3 + \alpha_{i+1/2}^{ei} \Delta x} = 1 - \frac{\sigma_{i+1/2} \Delta x}{2a_x \zeta_j^3 + \alpha_{i+1/2}^{ei} \Delta x}, $$

and

$$ \beta_{ij+1/2} = \frac{2a_\zeta \zeta_{j+1/2}^3}{2a_\zeta \zeta_{j+1/2}^3 + \alpha_i^{ei} \Delta \zeta} = 1 - \frac{\alpha_i^{ei} \Delta \zeta}{2a_\zeta \zeta_{j+1/2}^3 + \alpha_i^{ei} \Delta \zeta}, $$

leads to the following expression

$$ \frac{f_{1ij}^{n+1} - f_{1ij}^{n}}{\Delta t} = - \frac{(f_{2i+1j} - f_{2i-1j}) - a_x (f_{1i+1j} - f_{2ij} + f_{1i-1j})}{2\Delta x} \\ - \frac{(f_{2ij+1} - f_{2ij-1}) - a_\zeta (f_{1ij+1} - f_{2ij} + f_{1ij-1})}{2\Delta \zeta} \\ + \frac{a_x}{2\Delta x} \frac{\alpha_{i+1/2}^{ei} \Delta x}{2a_x \zeta_j^3 + \alpha_{i+1/2}^{ei} \Delta x} \left( (f_{1i+1j} + f_{1ij}) - \frac{f_{2i+1j} - f_{2ij}}{a_x} \right) \\ + \frac{a_x}{2\Delta x} \frac{\alpha_{i-1/2}^{ei} \Delta x}{2a_x \zeta_j^3 + \alpha_{i-1/2}^{ei} \Delta x} \left( (f_{1ij} + f_{1i-1j}) - \frac{f_{2ij} - f_{2i-1j}}{a_x} \right) \\ + \frac{a_\zeta}{2\Delta \zeta} \frac{\alpha_i^{ei} \Delta \zeta}{2a_\zeta \zeta_{j+1/2}^3 + \alpha_i^{ei} \Delta \zeta} \left( (f_{1ij+1} + f_{1ij}) - \frac{f_{2ij+1} - f_{2ij}}{a_\zeta} \right) \\ + \frac{a_\zeta}{2\Delta \zeta} \frac{\alpha_i^{ei} \Delta \zeta}{2a_\zeta \zeta_{j-1/2}^3 + \alpha_i^{ei} \Delta \zeta} \left( (f_{1ij} + f_{1ij-1}) - \frac{f_{2ij} - f_{2ij-1}}{a_\zeta} \right) \\ + \mathcal{O}(\Delta x) + \mathcal{O}(\Delta \zeta), $$
which is correctly consistent with the second equation of (3.5) as $\Delta x$ and $\Delta t$ tend to zero. Similarly, the first equation of (3.7) writes

$$ \frac{f_{0ij}^{n+1} - f_{0ij}^{n}}{\Delta t} = \frac{a_x}{\Delta x} f_{0i+1/2j}^{L*} + \frac{2a_x}{\Delta x} f_{0ij}^{n} + \frac{a_x}{\Delta x} f_{0i-1/2j}^{R*}. $$

Considering the definitions (3.13)-(3.15) this equation rewrites

$$ \frac{f_{0ij}^{n+1} - f_{0ij}^{n}}{\Delta t} = \frac{a_x}{\Delta x} (\tilde{f}_{0i+1/2j} - 2f_{0ij}^{n} + \tilde{f}_{0i-1/2j}) - \frac{a_x}{\Delta x} (\Gamma_{i+1/2j} \theta_{i+1/2j} - \Gamma_{i-1/2j} \theta_{i-1/2j}) \\ + \frac{a_\zeta}{\Delta\zeta} (\tilde{f}_{0ij+1/2} - 2f_{0ij}^{n} + \tilde{f}_{0ij-1/2}) - \frac{a_\zeta}{\Delta\zeta} (\Gamma_{ij+1/2} \theta_{ij+1/2} - \Gamma_{ij-1/2} \theta_{ij-1/2}) $$

and by using (3.11)-(3.14)-(3.16), one obtains the consistency with the first equation of (3.5). □

Now, in order to show the asymptotic-preserving property, in the limit $\tau^{e}/t^{*}$ tends to zero with the limit diffusion equation, one considers the diffusion scaling and a discrete Hilbet expansion is used.

**Theorem 3.5. (Consistency in the diffusion regime)**

In the diffusion limit, the numerical scheme (3.7) degenerates into

$$ \frac{f_{0ij}^{n+1,0} - f_{0ij}^{n,0}}{\Delta t} = \frac{\zeta_j}{\Delta x} \left[ \frac{\zeta_j^4}{6\sigma_{i+1/2}\Delta x} (f_{0i+1j}^{n,0} - f_{0ij}^{n,0}) - \frac{\zeta_j^4}{6\sigma_{i-1/2}\Delta x} (f_{0i1j}^{n,0} - f_{0i-1j}^{n,0}) \right. \\ \left. + \frac{\zeta_j^3 E_{i+1/2}}{6\sigma_{i+1/2}} \left(\frac{\partial f_0^{n,0}}{\partial \zeta}\right)_{i+1/2j} - \frac{\zeta_j^3 E_{i-1/2}}{6\sigma_{i-1/2}} \left(\frac{\partial f_0^{n,0}}{\partial \zeta}\right)_{i-1/2j} \right] \\ + \frac{E_i}{\Delta \zeta} \left[ \frac{E_i \zeta_{j+1/2}^3}{6\sigma_i \Delta \zeta} (f_{0ij+1}^{n,0} - f_{0ij}^{n,0}) - \frac{E_i \zeta_{j-1/2}^3}{6\sigma_i \Delta \zeta} (f_{0i1j}^{n,0} - f_{0ij-1}^{n,0}) \right. \\ \left. + \frac{\zeta_{j+1/2}^4}{6\sigma_i} \left(\frac{\partial f_0^{n,0}}{\partial x}\right)_{ij+1/2} - \frac{\zeta_{j-1/2}^4}{6\sigma_i} \left(\frac{\partial f_0^{n,0}}{\partial x}\right)_{ij-1/2} \right]. \quad (3.39) $$

*Proof.* Following the same approach as in section 3.1.2, equations (3.32) and (3.33) give

$$ \begin{aligned} f_{1i+1/2j}^{*,\epsilon} &= \alpha_{i+1/2j}^{\epsilon} \left[ \tilde{f}_{1i+1/2j}^{\epsilon} - \theta_{1i+1/2j} c_{i+1/2j} \left( \frac{\partial f_0^{\epsilon}}{\partial \zeta} \right)_{i+1/2j} (1 - \alpha_{i+1/2j}^{\epsilon}) \right], \\ f_{1ij+1/2}^{*,\epsilon} &= \beta_{ij+1/2}^{\epsilon} \left[ \tilde{f}_{1ij+1/2}^{\epsilon} - \theta_{2ij+1/2} \bar{c}_{ij+1/2} \left( \frac{\partial f_0^{\epsilon}}{\partial x} \right)_{ij+1/2} (1 - \beta_{ij+1/2}^{\epsilon}) \right], \end{aligned} \quad (3.40) $$

where $\alpha_{i+1/2j}^{\epsilon}$ and $\beta_{ij+1/2}^{\epsilon}$ are given by (3.22). Consequently it follows that

$$ f_{1i+1/2j}^{*,0} = 0 \quad \text{and} \quad f_{1ij+1/2}^{*,0} = 0. \quad (3.41) $$
Now, since the second component of (3.20) reads

$$\epsilon \frac{f_{1ij}^{n+1,\epsilon} - f_{1ij}^{n,\epsilon}}{\Delta t} = \frac{a_x}{\Delta x} f_{1i-1/2j}^{*,\epsilon} - \frac{2a_x}{\Delta x} f_{1ij}^{n,\epsilon} + \frac{a_x}{\Delta x} f_{1i+1/2j}^{*,\epsilon} + \frac{a_\zeta}{\Delta \zeta} f_{1ij-1/2}^{*,\epsilon} - \frac{2a_\zeta}{\Delta \zeta} f_{1ij}^{n,\epsilon} + \frac{a_\zeta}{\Delta \zeta} f_{1ij+1/2}^{*,\epsilon}$$

it follows that at order $\epsilon^0$ the previous equation gives

$$f_{1ij}^{n,0} = 0. \qquad (3.42)$$

Therefore one correctly recovers the condition (2.3) and from the definition (1.3), it follows that

$$f_{2ij}^{n,0} = f_{0ij}^{n,0} / 3.$$

The consistency with the limit diffusion equation is now proved. Firstly, one remarks that in the limit $\epsilon$ tends to zero, the results (3.41) and (3.42) give

$$\theta_{1i+1/2j} = 1, \qquad \theta_{2ij+1/2} = 1. \qquad (3.43)$$

Indeed, when $\epsilon$ tends to zero, the definitions (3.34) and (3.35) lead to

$$\tilde{\theta}_{1i+1/2j} = \frac{f_{0i+1j}^{n,0} + f_{0ij}^{n,0}}{|f_{0i+1j}^{n,0} - f_{0ij}^{n,0}|} \ge 1, \qquad \tilde{\theta}_{2ij+1/2} = \frac{f_{0ij+1}^{n,0} + f_{0ij}^{n,0}}{|f_{0ij+1}^{n,0} - f_{0ij}^{n,0}|} \ge 1.$$

Then from the definitions (3.13)-(3.15)-(3.11), we obtain

$$\left\{ \begin{aligned} f_{0i+1/2j}^{L*} &= f_{0ij}^{n} + (f_{1i+1/2j}^{L*} - f_{1ij}^{n}), \\ f_{0i-1/2j}^{R*} &= f_{0ij}^{n} + (f_{1ij}^{n} - f_{1i-1/2j}^{R*}), \end{aligned} \right. \qquad (3.44)$$

and

$$\left\{ \begin{aligned} f_{0ij+1/2}^{L*} &= f_{0ij}^{n} + (f_{1ij+1/2}^{L*} - f_{1ij}^{n}), \\ f_{0ij-1/2}^{R*} &= f_{0ij}^{n} + (f_{1ij}^{n} - f_{1ij-1/2}^{R*}), \end{aligned} \right. \qquad (3.45)$$

Secondly, the first component of (3.20) reads

$$\epsilon \frac{f_{0ij}^{n+1,\epsilon} - f_{0ij}^{n,\epsilon}}{\Delta t} = \frac{a_x}{\Delta x} f_{0i-1/2j}^{R*,\epsilon} - \frac{2a_x}{\Delta x} f_{0ij}^{n,\epsilon} + \frac{a_x}{\Delta x} f_{0i+1/2j}^{L*,\epsilon} + \frac{a_\zeta}{\Delta \zeta} f_{0ij-1/2}^{R*,\epsilon} - \frac{2a_\zeta}{\Delta \zeta} f_{0ij}^{n,\epsilon} + \frac{a_\zeta}{\Delta \zeta} f_{0ij+1/2}^{L*,\epsilon}$$

and by using (3.44)-(3.45) it follows that

$$\epsilon \frac{f_{0ij}^{n+1,\epsilon} - f_{0ij}^{n,\epsilon}}{\Delta t} = a_x \frac{f_{1i+1/2j}^{L*,\epsilon} - f_{1i-1/2j}^{R*,\epsilon}}{\Delta x} + a_\zeta \frac{f_{1ij+1/2}^{L*,\epsilon} - f_{1ij-1/2}^{R*,\epsilon}}{\Delta \zeta}.$$
Finally, the previous equation considered at order $\epsilon^0$ brings no information but at the order $\epsilon^1$ ones obtains

$$ \frac{f_{0ij}^{n+1,0} - f_{0ij}^{n,0}}{\Delta t} = a_x \frac{f_{1i+1/2j}^{L*,1} - f_{1i-1/2j}^{R*,1}}{\Delta x} + a_\zeta \frac{f_{1ij+1/2}^{L*,1} - f_{1ij-1/2}^{R*,1}}{\Delta \zeta}, $$

which gives, by using the definitions (3.32)-(3.33), the numerical scheme (3.39). □

### 3.3 General case with the term $\frac{E}{\zeta}(f_0 - f_2)$

As specified in part 3.1, in order to take into account the contribution of the source term $\frac{E}{\zeta}(f_0 - f_2)$, we simply propose to modify the intermediate states $f_{1i+1/2j}^*$ and $f_{1ij+1/2}^*$ given in (3.32) and (3.33) such that

$$ f_{1i+1/2j}^* = \alpha_{i+1/2j} \left[ \tilde{f}_{1i+1/2j} - \theta_{1i+1/2j} c_{i+1/2j} \left( \left( \frac{\partial f_0}{\partial \zeta} \right)_{i+1/2j} - \frac{\tilde{S}_{i+1/2j}}{2} \right) (1 - \alpha_{i+1/2j}) \right], \quad (3.46) $$

$$ f_{1ij+1/2}^* = \beta_{ij+1/2} \left[ \tilde{f}_{1ij+1/2} + \frac{\Delta \zeta}{2a_\zeta} S_{ij+1/2} - \theta_{2ij+1/2} \bar{c}_{ij+1/2} \left( \frac{\partial f_0}{\partial x} \right)_{ij+1/2} (1 - \beta_{ij+1/2}) \right], $$

with

$$ \tilde{S}_{i+1/2j} = \frac{\zeta_j^2}{3\alpha_i^{ei}} \frac{f_{0i+1j} + f_{0ij}}{2} \quad \text{and} \quad S_{ij+1/2} = \frac{E_i}{2} \left( \frac{f_{0ij+1} - f_{2ij+1}}{\zeta_{j+1}} + \frac{f_{0ij} - f_{2ij}}{\zeta_j} \right). $$

In this case, as in the previous part the coefficients $\theta_1$ and $\theta_2$ are also fixed to ensure the admissibility requirements.

**Theorem 3.6.** In the diffusive limit, the numerical scheme given by (3.7)-(3.13)-(3.15)-(3.46) is consistent with the limit equation (2.5).

*Proof.* The proof is the same than for Theorem 3, considering the intermediate states $f_{1i+1/2j}^*$ and $f_{1ij+1/2}^*$ given in (3.46). A direct calculation using the Hilbert expansions leads to the result. The terms $S_{ij+1/2}$ are consistent with the term $\frac{E}{\zeta}(f_0 - f_2)$ while the terms $\tilde{S}_{i+1/2j}$ enable to correctly recover the contribution of the two terms $\frac{E\zeta^2}{3\sigma} f_0^0$ in the $x$ and $\zeta$ derivatives of the limit equation. □

## 4 Numerical examples

In this section, the asymptotic-preserving scheme (3.7) is compared with the HLL scheme [37] and an explicit discretisation of the diffusion equation (2.5) in collisional regimes. We point out that no exact solution is available for a given value of the collisional parameter $\alpha^{ei}$. However, in very collisional regime, one expects the solution obtained with
the asymptotic-preserving scheme (3.7) to be close to the one obtained with an explicit discretisation of the diffusion equation (2.5). Similarly, in non-collisional regimes the solution obtained with the asymptotic-preserving scheme (3.7) is expected to be close to the one obtained with a HLL scheme which is known to behave correctly in such regimes. The time step used for the asymptotic-preserving scheme is computed with the CFL condition (3.36).

## 4.1 Relaxation of a Gaussian profile in the diffusion regime

In this example, the numerical scheme (3.7)-(3.13)-(3.15)-(3.46) is validated in the diffusive regime considering an inhomogeneous plasma with electric field. In this case, the initial conditions are the following

$$ \begin{cases} f_0(t=0,x,\zeta) = \zeta^2 \exp(-x^2) \exp(2(\zeta-3)^2), \\ f_1(t=0,x,\zeta) = 0. \end{cases} $$

The profile of $f_0$ at initial time as a function of $x$ and $\zeta$ is displayed in Figure 2.

<figure><img src="image_3.png" /><figcaption>Figure 2: Representation of the $f_0$ profile at the initial time as a function of $x$ and $\zeta$.</figcaption></figure>

For this test we have set $E=1$, $\alpha^{ei}=10^4$, the space range chosen is [-10,10] and the modulus velocity range [0,6]. The collisional parameter $\alpha_{ei}$ is large and chosen equals to $10^4$. Therefore we can consider that the limit $t_{ei}/t$ is reached particularly quickly and the particle transport can be described with the limit diffusion equation. In Figure 3, the solution obtained with the numerical scheme (3.7)-(3.13)-(3.15)-(3.46) is compared with the solution obtained with the HLL scheme and with an explicit discretisation of the limit diffusion equation (2.5) at different times. For all schemes we have set $\Delta x = 0.125$ and $\Delta \zeta = 0.2$. At time $t=1$, one remarks that the $f_0$ profile obtained with the HLL scheme is already seriously spread out while the profiles obtained with the AP scheme and the diffusion equation do not have changed. At time $t=50$, the AP scheme and diffusion equation discretisation $f_0$ profiles are spread out while the profile obtained with the HLL scheme has vanished. As observed at time $t=100$, in the long time regime, the AP scheme and
the discretisation of the diffusion equation behave identically. In addition a convergence study is presented in Figure 4. The $L^1$ and $L^2$ errors (in $x$ and $\zeta$) between the AP scheme and the diffusion scheme are displayed in green while the $L^1$ and $L^2$ errors between the HLL scheme and the diffusion scheme are displayed in red as functions of the collisional parameter $\alpha^{ei}$ at time $t=50$. Here this numerical convergence study shows that when the collisional parameter becomes large the $L^1$ and $L^2$ errors obtained with the AP scheme correctly tend to zero as expected while the $L^1$ and $L^2$ errors obtained with the HLL scheme miss the diffusion limit.

<figure><img src="image_4.png" /><figcaption>Figure 3: Representation of the $f_0$ profile as function of $x$ and $\zeta$ at time $t=1$ (top), $t=50$ (middle), $t=100$ (bottom), for the HLL scheme (left), AP scheme (middle) and the diffusion equation.</figcaption></figure>

## 4.2 Relaxation of a temperature profile in the diffusive regime with a self-consistent electric field

In this example, we consider the relaxation of a temperature profile in the diffusive regime considering a self-consistent electric field. The space range is [-40,40] and the velocity modulus range [0,6]. The collisional parameter $\alpha^{ei}$ is chosen equal to $10^4$. For all
<figure><img src="image_5.png" /><figcaption>Figure 4: Representation of the L¹ error (left) and L² error (right) (in x and ζ) between the AP scheme and the diffusion scheme (green) and the HLL scheme and the diffusion scheme (red) as a function of the collisional parameter a^ei at time t=50.</figcaption></figure>

schemes we have set Δx = 0.2 and Δζ = 0.125. The initial conditions are the following

$$ \begin{cases} f_0(t=0, x, \zeta) = \sqrt{\frac{2}{\pi}} \frac{\zeta^2}{T^{\text{ini}}(x)^{3/2}} \exp\left(-\frac{\zeta^2}{2T^{\text{ini}}(x)}\right), \\ f_1=0, \end{cases} \quad (4.1) $$

with $T^{\text{ini}}(x) = 2 - \arctan(x)$.

In this case the electric field is self-consistent meaning that at each time step it is calculated from the plasma profile. In this case we consider a Spitzer type model [9, 57], to evaluate the electric field

$$ E(x) = -\frac{dT(x)}{dx}, \quad (4.2) $$

where

$$ T(x) = \frac{1}{3n_e} \left( \int_{0}^{+\infty} \zeta^2 f_0 d\zeta - u^2 n_e \right), $$

with $n_e = \int_{0}^{+\infty} f_0 d\zeta$ and $u = \frac{1}{n_e} \int_{0}^{+\infty} f_1 \zeta d\zeta$.

In Figure 5, the temperature profile is displayed at the initial time and at time t=80. The temperature profiles obtained with the HLL scheme, the AP scheme and a discretisation of the diffusion equation (2.5) are compared at time t=80. On one hand, one remarks that
the HLL temperature profile is excessively spread out compared to the AP and diffusion profiles while on the other hand the AP and diffusion profiles match exactly at time t=80. This example demonstrates the inability of the HLL scheme in capturing the correct temperature profile while the AP scheme presented handle perfectly the diffusion limit regime.

<figure><img src="image_6.png" /><figcaption>Figure 5: Representation of the temperature profile as a function of x at time t=80.</figcaption></figure>

## 4.3 Two electron beams interaction

In this example the interaction between two electron beams is considered. This collisionless test case enables us to validate the AP scheme (3.7)-(3.13)-(3.15)-(3.46) in a regime where electrostatic effects are predominant compared to the collisional effects, therefore we set $\alpha^{ei} = 0$. For the AP and the HLL schemes we have set $\Delta x = 0.125$ and $\Delta \zeta = 0.2$. Consider two electron beams propagating at velocity $v_0$ and $v_1$. In that case, the dispersion relation [36] is given by

$$1 - \frac{1}{(\omega - k v_0)^2} - \frac{1}{(\omega - k v_1)^2} = 0,$$

where $v_0$ and $v_1$ denote the beams velocities.

This configuration can lead to electrostatic instabilities [18, 36]. Indeed, the solutions of the form $Ae^{i\omega t+ikx}$ are unstable when $\omega_1$ the imaginary part of $\omega$ is strictly positive. In the case $v_0 = -v_1$ we can show that the solution is stable if $kv_0 \ge \sqrt{2}$.

This test is problematic for the $M_1$ model. Indeed, if we consider two electron beams propagating with opposite velocities the distribution function is well defined. Nevertheless, the $M_1$ model considers only the angular moments $f_0$ and $f_1$. For the calculation
of $f_1$ the two populations contributions cancel and we get $f_1 = 0$. The $M_1$ model sees an unrealistic isotropic configuration. To overcome this problem the superposition principle is used since the model is linear [36, 61]. Two particle populations (one per beam) are considered. For each time step the $M_1$ problem is solved for the first population then for the second one. Hence the electrostatic field is calculated using the Maxwell-Ampere solved taking into account the two distribution functions. This approach was validated for the present test case in [34].

In the case of two streams propagating with opposite velocities $v_d$ and $-v_d$, the initial electron distribution function is the following

$$f(t=0,x,v) = 0.5[(1+A\cos(kx))M_{v_d}(v) + (1-A\cos(kx))M_{-v_d}(v)],$$

with

$$M_{\pm v_d}(v) = \exp\left(-\frac{(v \mp v_d)^2}{2}\right).$$

The first corresponding angular moments $f_0^1$ and $f_0^2$ of the first and second population read

$$\begin{cases} f_0^1(t=0, x, \zeta) = 0.5(1+A\cos(kx))\frac{\zeta}{v_d}\left(\exp\left(-\frac{(\zeta-v_d)^2}{2}\right) - \exp\left(-\frac{(\zeta+v_d)^2}{2}\right)\right), \\ f_0^2(t=0, x, \zeta) = 0.5(1-A\cos(kx))\frac{\zeta}{v_d}\left(\exp\left(-\frac{(\zeta-v_d)^2}{2}\right) - \exp\left(-\frac{(\zeta+v_d)^2}{2}\right)\right). \end{cases}$$

The second angular moments $f_1^1$ and $f_1^2$ of the first and second population read

$$\begin{cases} f_1^1(t=0, x, \zeta) = 0.5(1+A\cos(kx))\frac{1-\zeta v_d}{v_d^2}\left(\exp\left(-\frac{(\zeta-v_d)^2}{2}\right) - \exp\left(-\frac{(\zeta+v_d)^2}{2}\right)\right), \\ f_1^2(t=0, x, \zeta) = -0.5(1-A\cos(kx))\frac{1-\zeta v_d}{v_d^2}\left(\exp\left(-\frac{(\zeta-v_d)^2}{2}\right) - \exp\left(-\frac{(\zeta+v_d)^2}{2}\right)\right). \end{cases}$$

At each time step, the electrostatic field is computed using the Maxwell-Ampere equation considering the contribution of the two population of particles

$$\frac{dE}{dt} = \int_{0}^{+\infty} f_{1}^{1}\zeta d\zeta + \int_{0}^{+\infty} f_{1}^{2}\zeta d\zeta.$$

The parameter $A$ is introduced to perturb the initial condition in order to enable the development of the electrostatic instability. The velocity modulus range chosen is [0,12] and the space range is [0,25]. In this example we set $v_d=4$, $A=0.001$ and periodical boundary conditions are used. The results have been compared with a kinetic code [24]. In Figure 6, the evolution of the electrostatic energy is displayed as a function of time using the AP scheme in red and the kinetic code in dashed blue. The AP scheme and the kinetic code give analogous results. This numerical experiment shows the good behaviour of the AP scheme in a regime where electrostatic effects are predominant.
<figure><img src="image_7.png" /><figcaption>Figure 6: Representation of the temporal evolution of the electrostatic energy.</figcaption></figure>

## 4.4 Relaxation of a temperature profile in the diffusive regime with a self-consistent electric field and non-constant collisional parameter

In this example, the initial conditions are the same as for the previous example where the initial temperature profile is given by (4.1) and the electric field is computed using (4.2). In this case the collisional parameter $\alpha^{ei}$ is not constant and follows the linear profile

$$ \alpha^{ei}(x) = ax + b, $$

with $\alpha^{ei}(x_{min} = -40) = 5 \cdot 10^3$ and $\alpha^{ei}(x_{max} = 40) = 10^5$. It follows that the coefficients $a$ and $b$ read

$$ a = \frac{10^5 - 5 \cdot 10^3}{x_{max} - x_{min}}, \qquad b = 5 \cdot 10^3 - ax_{min}. $$

The space range is [-40,40] and the velocity modulus range [0,6]. For this test case we have set $\Delta x = 0.4$ and $\Delta \zeta = 0.25$ for all the schemes. In Figure 7, the temperature profile is displayed at the initial time and at time t=5000 for the AP scheme and an explicit discretisation of the diffusion equation (2.5). After a long time (t=5000) and despite the strong spatial variation of the function $\alpha^{ei}$ the AP and diffusion profiles give very close result. One remark on the space interval [-40,0] the AP curve in red is slightly different to the diffusion curve in dashed blue while on the interval [0,40] the results match perfectly. This could be explained as the collisional parameter $\alpha^{ei}$ becomes larger for important $x$, therefore, the limit diffusive regime is fully reached for large $x$ where the comparison with the diffusion equation is valid.

## 4.5 Case a non-constant self-consistent collisional parameter

When considering physical relevant configurations occurring in plasma physics, the collisional parameter depends of the state of the plasma. The knowledge of the ionic and electronic distribution function is required to compute the collisional parameter. Therefore in this test case, we choose to consider a nonlinear collisional parameter which de-
<figure><img src="image_8.png" /><figcaption>Figure 7: Representation of the temperature profile at time t=5000.</figcaption></figure>

pends of the solution itself

$$ \alpha^{ei}(t,x,\zeta) = \exp(f_0(t,x,\zeta) + f_1(t,x,\zeta)). $$

In this case, $E=1$, the space range chosen is [-10,10] and the velocity modulus range [0,6]. For this test case we have set $\Delta x = 0.33$ and $\Delta \zeta = 0.125$. The initial condition is given by

$$ \begin{cases} f_0(t=0, x, \zeta) = \zeta^2 \exp(-(\zeta-3)^2) \exp(-x^2/10), \\ f_1 = 0. \end{cases} $$

We consider periodical boundary conditions. In Figure 8, the initial profile of $f_0$ is displayed at the initial time and at time t=3.

## 5 Conclusion

In this work, an asymptotic-preserving scheme has been proposed for the electronic $M_1$ model in the diffusive limit. In order to deal with the mixed derivatives which arise in the diffusive limit an anisotropic numerical viscosity has been considered. The numerical scheme preserves the realisibility domain and captures the correct limit equation. The contribution of the source term $E(f_0 - f_2)/\zeta$ is integrated and the cases of non constant electric field and collisional parameter are naturally included. Numerical examples have been performed in non-collisional and diffusive regimes. It has been observed that the present scheme behaves correctly in both regimes.
<figure><img src="image_9.png" /><figcaption>Figure 8: Representation of the $f_0$ profile as function of $x$ and $\zeta$ at the initial time (left) and $t=3$ (right).</figcaption></figure>

A possible perspective could be to consider an electron-electron collisional operator or the study of the coupling with the Maxwell's equations.

References

[1] G.W. Alldredge, C.D. Hauck, and A.L. Tits. High-order entropy-based closures for linear transport in slab geometry II: A computational study of the optimization problem. SIAM Journal on Scientific Computing Vol. 34-4 (2012), pp. B361-B391.

[2] D. Aregba-Driollet, M. Briani, and R. Natalini. Asymptotic High-Order Schemes for 2x2 Dissipative Hyperbolic Systems. SIAM Journal on Numerical Analysis 46 (2008), 869-894.

[3] R. Balescu. Transport Processes in Plasma, Elsevier, (Amsterdam, 1988), Vol. 1.

[4] C. Berthon, P. Charrier, and B. Dubroca. An asymptotic preserving relaxation scheme for a moment model of radiative transfer. C. R. Acad. Sci. Paris, Ser. I 344 (2007).

[5] C. Berthon, P. Charrier, and B. Dubroca. An HLLC Scheme to Solve The M1 Model of Radiative Transfer in Two Space Dimensions. Journal of Scientific Computing, Vol. 31, No. 3, (2007).

[6] C. Berthon, G. Moebs, and R. Turpault. An Asymptotic-Preserving Scheme for Systems of Conservation Laws with Source Terms on 2D Unstructured Meshes. Finite Volumes for Complex Applications VII-Methods and Theoretical Aspects Volume 77 of the series Springer Proceedings in Mathematics and Statistics pp 107-115.

[7] C. Berthon and R.Turpault. Asymptotic preserving HLL schemes. Numerical Methods for Partial Differential Equations, 27 (6) (2011) 1396-1422.

[8] F. Bouchut, H. Ounaissa, and B. Perthame. Upwinding of the source term at interfaces for Euler equations with high friction. Computers and Mathematics with Applications 53 (2007), 361-375.

[9] S.I. Braginskii. Reviews of Plasma Physics. M.A Leontovich, Ed., Consultants Bureau (New York, 1965), Vol. 1, p.205.

[10] A.V. Brantov, V.Yu. Bychenkov, O.V. Batishchev, and W.Rozmus. Nonlocal heat wave propagation due to skin layer plasma heating by short laser pulses. Computer Physics communications 164 67, (2004).

[11] C. Buet and S. Cordier. Conservative and entropy decaying numerical scheme for the isotropic Fokker-Planck-Landau equation. J. Comput. Phys. 145, No.1, 228-245 (1998).
[12] C. Buet and B. Després. Asymptotic preserving and positive schemes for radiation hydro-dynamics. J. Compt. Phy., 215, 717740 (2006).

[13] C. Buet, B. Després, and E. Franck. Design of asymptotic preserving finite volume schemes for hyperbolic heat equation on unstructured meshes. Numerische Mathematik, Volume 122, Issue 2, pp 227-278 (2012).

[14] C. Chalons, F. Coquel, E. Godlewski, P.-A. Raviart, and N. Seguin. Godunov-type schemes for hyperbolic systems with parameter-dependent source. The case of Euler system with friction. Math. Models Methods Appl. Sci. 20 (2010), no. 11, 21092166. MR 2740716 (2011m:65179).

[15] S. Chapman. Phil. Trans. Roy. Soc. London 216 (1916) 279.

[16] S. Chapman and T. G. Cowling. The Mathematical Theory of Non-Uniform Gases. Cambridge University Press, Cambridge, England, 1995.

[17] P. Charrier, B. Dubroca, G. Duffa, and R. Turpault. Multigroup model for radiating flows during atmospheric hypersonic re-entry. Proceedings of International Workshop on Radiation of High Temperature Gases in Atmospheric Entry, pp. 103110. Lisbonne, Portugal. (2003).

[18] F. Chen. *Introduction to Plasma Physics and Controlled Fusion*, (1984). Plenum Press, New York.

[19] G Dimarco and L Pareschi. Asymptotic Preserving Implicit-Explicit Runge-Kutta Methods for Nonlinear Kinetic Equations. SIAM Journal on Numerical Analysis 51 (2013), 1064-1087.

[20] J.F. Drake, P.K. Kaw, Y.C. Lee, G. Schmidt, C.S. Liu, and M.N. Rosenbluth. Parametric instabilities of electromagnetic waves in plasmas. Phys. Fluids 17, 778, 1974.

[21] B. Dubroca, J.-L. Feugeas, and M. Frank. Angular moment model for the Fokker-Planck equation. European Phys. Journal D, 60, 301, (2010).

[22] B. Dubroca and J.L. Feugeas. Entropic moment closure hierarchy for the radiative transfert equation. C. R. Acad. Sci. Paris Ser. I, 329 915, (1999).

[23] B. Dubroca and J.L. Feugeas. Étude théorique et numérique d'une hiéarchie de modèles aux moments pour le transfert radiatif. C. R. Acad. Sci. Paris, t. 329, SCrie I, p. 915-920, (1999).

[24] R. Duclous, B. Dubroca, F. Filbet, and V. Tikhonchuk. High order resolution of the Maxwell-Fokker-Planck-Landau model intended for ICF applications. J. Comput. Phys. 228(14): 5072-5100 (2009).

[25] D. Enskog. Kinetische Theorie der Vorgnge in Mssig Verdnnten Gasen. Uppsala, 1917.

[26] Reduced entropic model for studies of multidimensional nonlocal transport in high-energy-density plasmas. D. Del Sorbo and J.-L. Feugeas and Ph. Nicolai and M. Olazabal-Loumé and B. Dubroca and S. Guisset and M. Touati and V. Tikhonchuk. Phys. Plasmas 22, 082706 (2015).

[27] E. Epperlein and R. Short. Phys. Fluids B 4, 2211 (1992).

[28] F. Filbet and S. Jin. A class of asymptotic preserving schemes for kinetic equations and related problems with stiff sources. J. Comp. Phys. vol. 229, no 20 (2010).

[29] L. Gosse and G. Toscani. An asymptotic-preserving well-balanced scheme for the hyperbolic heat equations. C. R. Math. Acad. Sci. Paris 334 (4) (2002) 337342.

[30] L. Gosse and G. Toscani. Space localizaion and well-balanced schemes for discrete kinetic models in diffusive regimes. SIAM J. Numer. Anal. 41 (2) (2003) 641658.

[31] H. Grad. On the kinetic theory of rarefied gases. Commun. Pure Appl. Math. 2, 331-407 (1949).

[32] C.P.T. Groth and J.G. McDonald. Towards physically-realizable and hyperbolic moment closures for kinetic theory. Continuum Mech. Thermodyn. 21, 467-493 (2009).

[33] S. Guisset, S. Brull, E. d'Humières, and B. Dubroca. Asymptotic-preserving well-balanced
scheme for the electronic $M_1$ model in the diffusive limit: particular cases. ESAIM: M2AN, 51 5 (2017) 1805-1826.

[34] S. Guisset, S. Brull, B. Dubroca, E. d'Humières, S. Karpov, and I. Potapenko. Asymptotic-preserving scheme for the Fokker-Planck-Landau-Maxwell system in the quasi-neutral regime. Communications in Computational Physics, volume 19, issue 02, pp. 301-328 (2016).

[35] S. Guisset, S. Brull, E. d'Humires, B. Dubroca, and V. Tikhonchuk. Classical transport theory for the collisional electronic M1 model. Physica A: Statistical Mechanics and its Applications, Volume 446, Pages 182-194 (2016).

[36] S. Guisset, J.G. Moreau, R. Nuter, S. Brull, E. d'Humières, B. Dubroca, and V.T. Tikhonchuk. Limits of the M1 and M2 angular moments models for kinetic plasma physics studies. J. Phys. A: Math. Theor. 48, 335501 (2015).

[37] A. Harten, P.D. Lax, and B. Van Leer. On upstream differencing and Godunov-type schemes for hyperbolic conservation laws. SIAM Review 25 (1983), 35-61.

[38] S. Jin and C.D. Levermore. Fully discrete numerical transfer in diffusive regimes. Transport Theory Statist. Phys. 22 (6) 739791. (1993).

[39] S. Jin and C.D. Levermore. The discrete-ordinate method in diffusive regimes. Transport Theory Statist. Phys. 20 (56) 413439. (1991).

[40] P. Lafitte and G. Samaey. Asymptotic-preserving projective integration schemes for kinetic equations in the diffusion limit. SIAM Journal on Scientific Computing, 34(2):A579 A602, 2012.

[41] L. Landau. On the vibration of the electronic plasma. J. Phys. USSR 10 (1946).

[42] M. Lemou and L. Mieussens. A New Asymptotic Preserving Scheme Based on Micro-Macro Formulation for Linear Kinetic Equations in the Diffusion Limit. SIAM J. Sci. Comput., 31(1), 334368, 2008.

[43] C.D. Levermore. Moment closure hierarchies for kinetic theories. J. Stat. Phys. 83, 1021-1065 (1996).

[44] J. Mallet, S. Brull, and B. Dubroca. An entropic scheme for an angular moment model for the classical Fokker-Planck-Landeau equation of electrons. Comm. Comput. Phys., 422, (2013).

[45] J. Mallet, S. Brull, and B. Dubroca. General moment system for plasma physics based on minimum entropy principle. Kinetic and Related Models, vol. 8, No.3, 533-558, (2015).

[46] A. Marocchino, M. Tzoufras, S. Atzeni, A. Schiavi, Ph. D. Nicola, J. Mallet, V. Tikhonchuk, and J.-L. Feugeas. Nonlocal heat wave propagation due to skin layer plasma heating by short laser pulses. Phys. Plasmas 20, 022702, (2013).

[47] J.G. McDonald and C.P.T. Groth. Towards realizable hyperbolic moment closures for viscous heat-conducting gas flows based on a maximum-entropy distribution. Continuum Mech. Thermodyn. 25, 573-603 (2012).

[48] N. Meezan, L. Divol, M. Marinak, G. Kerbel, L. Suter, R. Stevenson, G. Slark, and K. Oades. Phys. Plasmas 11, 5573 2004.

[49] G.N. Minerbo. Maximum entropy Eddington Factors. J. Quant. Spectrosc. Radiat. Transfer, 20, 541, (1978).

[50] I. Muller and T. Ruggeri. Rational Extended Thermodynamics. Springer, New York (1998).

[51] Ph. D. Nicola, J.-L. A. Feugeas, and G. P. Schurtz. A practical nonlocal model for heat transport in magnetized laser plasmas. Phys. Plasmas 13, 032701, (2006).

[52] J.-F. Ripoll. An averaged formulation of the M1 radiation model with presumed probability density function for turbulent flows. J. Quant. Spectrosc. Radiat. Trans. 83 (34), 493517. (2004).

[53] W. Rozmus, V. T. Tikhonchuk, and R. Cauble. A model of ultrashort laser pulse absorption
in solid targets. Phys. Plasmas 3, 360 (1996).

[54] S.Boscarino, P.G. LeFloch, and G. Russo. High-order asymptotic-preserving methods for fully nonlinear relaxation problems. SIAM J. Sci. Comput. Vol. 36, No.2, pp.A377-A395.

[55] K. Shigemori, H. Azechi, M. Nakai, M. Honda, K. Meguro, N. Miyanaga, H. Takabe, and K. Mima. Phys. Rev. Lett. 78, 250 (1997).

[56] I.P. Shkarofsky, T.W. Johnston, and The Particle Kinetics of Plasmas M.P. Bachynski. Addison-Wesley (Reading, Massachusetts, 1966).

[57] L. Spitzer and R. Härm. Phys. Rev. 89 (1953) 977.

[58] H. Struchtrup. Macroscopic Transport Equations for Rarefied Gas Flows. Springer, Berlin (2005).

[59] M. Touati, J.L. Feugeas, P. Nicolai, J.J. Santos, L. Gremillet, and V.T. Tikhonchuk. New Journ. of Physics 16 (2014).

[60] R. Turpault. A consistent multigroup model for radiative transfer and its underlying mean opacity. J. Quant. Spectrosc. Radiat. Transfer 94, 357-371 (2005).

[61] R. Turpault, M. Frank, B. Dubroca, and A. Klar. Multigroup half space moment approximations to the radiative heat transfer equations. J. Comput. Phys. 198 363 (2004).

[62] A. Velikovich, J. Dahlburg, J. Gardner, and R. Taylor. Phys. Plasmas 5, 1491 (1998).