# Hybrid Kinetic/Fluid numerical method for the Vlasov-BGK equation in the diffusive scaling

Tino Laidin

► To cite this version:

Tino Laidin. Hybrid Kinetic/Fluid numerical method for the Vlasov-BGK equation in the diffusive scaling. Kinetic and Related Models, 2023, 16 (6), pp.913-947. 10.3934/krm.2023013. hal-03560784v2

HAL Id: hal-03560784

https://hal.science/hal-03560784v2

Submitted on 30 Sep 2022

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
# HYBRID KINETIC/FLUID NUMERICAL METHOD FOR THE VLASOV-BGK EQUATION IN THE DIFFUSIVE SCALING

TINO LAIDIN

Univ. Lille, CNRS, Inria, UMR 8524 - Laboratoire Paul Painlevé, F-59000 Lille,
France

**ABSTRACT.** This work presents a hybrid numerical method for linear collisional kinetic equations with diffusive scaling. The aim of the method is to reduce the computational cost of kinetic equations by taking advantage of the lower dimensionality of the asymptotic fluid model while reducing the error induced by the latter approach. It relies on two criteria motivated by a perturbative approach to obtain a dynamic domain adaptation. The first criterion quantifies distance between a local equilibrium in velocity and the distribution function of particles. The second one depends only on the macroscopic quantities that are available on the whole computing domain. A key idea is the use of a micro-macro decomposition to deal with interface conditions. The method is significantly more efficient than a standard full kinetic approach. Some properties, such as the conservation of mass, are also investigated and illustrated through various examples.

## 1. INTRODUCTION

Systems of particles in interactions arise in many fields of science such as gas theory, plasmas physic or even semiconductors. The mathematical models for such systems can be classified in three major scales: particles (microscopic), kinetic (mesoscopic), and fluid (macroscopic). The first type of model is about describing the system as point particles interacting with each other via collisions or electromagnetic forces. Although the resulting system is the most realistic, it is in practice extremely large. Its study both theoretically and numerically becomes unattainable. In this work we consider a kinetic description of the system modelled by the Vlasov-BGK equation. The unknown is the probability distribution $f(t, x, v)$ solution to:

$$ \begin{cases} \frac{\partial}{\partial t} f + v \cdot \nabla_x f + E \cdot \nabla_v f = Q(f), \\ f(0, x, v) = f_0(x, v). \end{cases} $$

The short-range interactions between particles are taken into account through the collision operator $Q(f)$, and the function $E(x)$ is a given exterior electrical field. While attainable, the simulation of this equation remains expensive in computational resources. Moreover, using the kinetic description of the system may not

MATHEMATICS SUBJECTS CLASSIFICATION: Primary 65M08, 82M12; Secondary 35B40, 65M55
KEYWORDS: Kinetic equations; Diffusion scaling; Asymptotic preserving scheme; Micro-macro decomposition; Hybrid solver
be necessary on the whole computing domain. The fluid one, that is less precise but much less costly, can be used where it is accurate enough. The aim of the paper is therefore to design a hybrid kinetic/fluid scheme with an automatic domain adaptation method. It relies on a robust numerical scheme for the kinetic equation, on relevant criteria to carefully determine fluid and kinetic regions and on a smart implementation.

**Diffusive scaling.** In some applications, it is relevant to consider a scaled version of the Vlasov-BGK equation. Let us introduce the scaling parameter $\epsilon > 0$. It is related to the Knudsen number: the ratio between the mean free path of the particles and the length scale of observation. This work will focus on the diffusive scaling. Let $d_x \ge 1$ and $d_v \in \mathbb{N}^*$ be integers. We denote by $\Omega_x$ a subset of $\mathbb{R}^{d_x}$. Let $t \ge 0$, $x \in \Omega_x$ and $v \in \mathbb{R}^{d_v}$. We look for a particle distribution function $f^\epsilon$ solution to the following scaled equation:

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial}{\partial t} f^{\epsilon} + \frac{v}{\epsilon} \cdot \nabla_x f^{\epsilon} + \frac{E}{\epsilon} \cdot \nabla_v f^{\epsilon} = \frac{1}{\epsilon^2} Q(f^{\epsilon}), \\ f^{\epsilon}(0, x, v) = f_0(x, v). \end{array} \right. \qquad (P^{\epsilon}) $$

We assume that the initial condition $f_0$ is nonnegative and does not depend on $\epsilon$. In practice, the Knudsen number can be of order 1 down to 0 depending on the physics being modelled. On the one hand, when $\epsilon \sim 1$, the system is said to be in the kinetic regime. It models a system with few collisions between particles. On the other hand, when $\epsilon \ll 1$, the system reaches the fluid regime. This scaling and its asymptotic limit have first been studied in [5]. The asymptotic expansion of the distribution function $f^\epsilon$ in $\epsilon$ is justified in [3] for the neutron transport and in [36] for the linear Boltzmann equation. In [11] a large class of linear collision operators is dealt with and in [21], the authors justified an approximation of the kinetic equation by diffusion using homogenization. In our setting, the limit case $\epsilon = 0$ is described by a drift-diffusion equation on the density $\rho(t, x)$:

$$ \left\{ \begin{array}{l} \displaystyle \partial_t \rho - \mathrm{div}_x (\nabla_x \rho - E\rho) = 0, \\ \rho(0, x) = \rho_0(x). \end{array} \right. \qquad (P) $$

**Asymptotic preserving scheme.** In order to design a numerical method that is efficient, one first needs a numerical scheme that performs well for any value of the parameter $\epsilon$. Indeed, as $\epsilon$ tends to 0, the transport velocity in $(P^\epsilon)$ formally goes to infinity. Numerically, it translates to smaller and smaller time steps to guarantee the stability of a naive scheme and, consequently, a reasonable computation time cannot be ensured. A solution to this issue is to use schemes that remain stable in the diffusive limit $\epsilon \to 0$. These schemes fall into the framework of Asymptotic Preserving (AP) schemes, a notion introduced in [29] and [26]. We also refer to the recent review articles [27, 15]. This AP property can be summarized by the diagram in Figure 1. In the diagram, $\rho$ corresponds to a solution to the problem (P) and $\rho_h$ is an approximation of $\rho$. On the other hand, $f^\epsilon$ is a solution to the problem $(P^\epsilon)$ and $f_h^\epsilon$ is an approximation of $f^\epsilon$. The idea behind AP scheme is threefold. Firstly, the scheme for $(P^\epsilon)$ has to be a consistent discretization of the limit model as $\epsilon \to 0$. Secondly, a scheme is considered truly AP only if the stability criterion on the time step is independent on the parameter $\epsilon$. Thirdly, one can explicitly take $\epsilon = 0$ in the scheme. The need of an AP scheme in the kinetic domain of the hybrid scheme is crucial. One the one hand, for computation time considerations. One the other hand, the limit scheme is used in the fluid regions of the domain
adaptation to ensure good transitions between kinetic and fluid states. While AP schemes are designed to resolve both the mesoscopic and the macroscopic scales automatically, it often implies more expensive computations even in a fluid regime because of the resolution of the kinetic scale. By using a hybrid method, one can effectively take advantage of the properties of an AP scheme while limiting its use and therefore reduce the computation time.

**Hybrid method.** A strategy to take advantage of both the kinetic and fluid scales of description and reduce the computational time is to use a hybrid method. The notion of multiscale coupling has already been studied and a wide range of techniques has been developed. We will rely on the technique that consists in adapting the domain in the position variable. Kinetic and fluid regions are created and move throughout the simulations and the appropriate solver can then be used in each subdomains. The key idea is the definition of the interfaces between subdomains. Various criteria have been investigated and we refer to [37, 39, 40, 10, 31] and the references therein. These criteria vary from being based purely on macroscopic quantities (macroscopic criteria) to considering the deviation of the distribution function from a local equilibria in velocity (kinetic criteria). In particular, the idea to use the asymptotic limit of the kinetic model to achieve a domain adaptation can be found in [28]. In this work, we shall adapt to the diffusive scaling the criterion introduced in [34] and later used in [38] and [19]. In particular, we will use a combination of macroscopic and kinetic criteria to achieve the dynamic domain adaptation. We also refer to [31] where in addition to the coupling between the kinetic and fluid scales, mesh refinement is also considered. Another technique, first introduced for the hydrodynamic scaling in [12], also relies on domain adaptation but consists in adding a buffer zone between kinetic and fluid subdomains using a transition function in the continuous model. It was extended to a dynamic setting with moving interfaces in [10] and applied to the diffusive scaling in [13]. Another approach that do not rely on the adaptation of the domain in position [14, 8] consists in splitting the distribution function into a macroscopic part solved via the finite volume method and a microscopic one solved using a Monte-Carlo method. In [24, 25] the splitting of the distribution function is also considered in a stationary setting. While the microscopic part is again solved via Monte-Carlo methods, the macroscopic one is solved using a moment method. Another technique introduced in [43, 42] relies on meshfree methods where the microscopic part is solved using particles and a lagrangian scheme is used for the macroscopic part. Finally, an alternative approach to reduce the computational cost relies on dynamical low-rank methods [18] where the idea is to reduce the matrix of the linear system associated to the discrete problem.

**Main contribution.** In this work, we develop a hybrid kinetic/fluid numerical method with an automatic domain adaptation. In particular, no buffer zones are introduced and we consider a full finite volume approach. The criteria for the dynamic interfaces are adapted from [19]. The domain adaptation is accomplished cell-wise and the mesh is fixed throughout the simulation. Domain adaptation methods are heavily dependent on the treatment of interfaces. One of the main novelty of this work is the way to deal with interface conditions between kinetic and fluid. This is achieved by using a micro-macro reformulation of the kinetic equation. Moreover, the conservation of mass of the hybrid method is thoroughly investigated for a relevant toy model. The accuracy and long-time behaviour of the
<figure><img src="image_2.png" /><figcaption>FIGURE 1. The AP diagram (h denotes the size of the discretization)</figcaption></figure>

method are observed numerically. Finally, the speedup provided by the method is significant and illustrated through several test cases.

**Plan of the paper.** The outline is as follows. Section 2 is dedicated to the derivation of a hierarchy of macroscopic models based on the Chapman-Enskog expansion of the distribution function. In Section 3 we recall the micro-macro reformulation of the Vlasov-BGK equation. This reformulation is then used to develop an Asymp totic Preserving scheme with a finite volume approach. Section 4 is dedicated to the hybrid method. The coupling indicators based on the hierarchy introduced in Section 2 are presented and the implementation of the hybrid scheme is discussed. Finally, numerical experiments are performed in Section 5.

## 2. CHAPMAN-ENSKOG EXPANSION

The aim of this section is to derive a hierarchy of macroscopic models from which we will deduce a macroscopic criterion.

### 2.1. Notations and functional setting.
From now on, we consider periodic boundary conditions in position and let $v \in \mathbb{R}^{d_v}$, $d_v \in \mathbb{N}^*$. We also assume that the electrical field $E$ is periodic on $[0, x_*)$. From now on, the collision operator will be the linear BGK operator [7]:

$$Q(f) = \rho M - f, \quad \text{where } \rho = \langle f \rangle \text{ and } \langle f \rangle = \int_{\mathbb{R}^{d_v}} f \, dv. \qquad (2.1)$$

The notation $\langle \cdot \rangle$ will be used either for scalars or component-wise for higher order tensors. Here $M(v)$ denotes the so-called Maxwellian. A standard function that we consider in this work is the multidimensional centered Gaussian:

$$M(v) = \frac{e^{-|v|^2/2}}{(2\pi)^{d/2}}.$$

Let us recall some properties of $M$:

$$\left\{ \begin{array}{l} M(v) > 0, \quad \forall v \in \mathbb{R}^{d_v}, \\ \langle M \rangle = 1. \end{array} \right. \qquad (2.2)$$

Moreover, since $M$ can be expressed as the product of $d_v$ 1-dimensional Gaussians, it is isotropic and even in each direction. As a consequence of the symmetric domain in velocity and the symmetry of Gaussian, its odd moments vanish. The
1-dimensional Gaussian also admits finite zeroth, second, and fourth moments in velocity. We denote by $m_k$ its k-th moment:

$$m_k = \frac{1}{\sqrt{2\pi}} \int_{\mathbb{R}} v_x^k \exp\left(-\frac{v_x^2}{2}\right) dv_x.$$

Let us recall the integro-differential problem we are interested in:

$$\left\{ \begin{array}{l} \partial_t f^\varepsilon + \frac{1}{\varepsilon} \mathcal{T}(f^\varepsilon) = \frac{1}{\varepsilon^2}(\rho^\varepsilon \mathcal{M} - f^\varepsilon), \\ f^\varepsilon(0, x, v) = f_0(x, v), \end{array} \right. \qquad (\text{VBGK})$$

where

$$\mathcal{T}(f) = v \cdot \nabla_x f + E \cdot \nabla_v f.$$

Let us set $\gamma(v) = \frac{1}{\mathcal{M}(v)}$ and introduce the measure

$$d\gamma = \gamma(v) dv = \frac{dv}{\mathcal{M}(v)}.$$

Let $L^2(dx d\gamma)$ the space of square integrable functions against the measure $dx d\gamma$ equipped with the scalar product

$$(f_1, f_2)_{L^2(dx d\gamma)} = \int_{\Omega_x \times \mathbb{R}^{d_v}} f_1 f_2 \, dx \, d\gamma.$$

With an initial data in $L^2(dx d\gamma)$, there is a unique solution to (VBGK) (see, e.g., [2]) which conserves mass and nonnegativity.

One can define the null space of the linear BGK operator (2.1) as

$$\mathcal{N} = \{f = \rho \mathcal{M} \text{ where } f \in L^2(dx d\gamma), \rho = \langle f \rangle\}.$$

The space $\mathcal{N}$ is sometimes referred to as the equilibrium manifold. In particular, one has that

$$\mathcal{N}^{\perp} = \{f \in L^2(dx d\gamma) \text{ such that } \langle f \rangle = 0\}.$$

With these notations, one can decompose $f$ as its equilibrium part in $\mathcal{N}$ plus a perturbative part in $\mathcal{N}^{\perp}$. Note that the perturbative part is not necessarily small. Let us now introduce the so-called Chapman-Enskog expansion of the distribution function $f^{\varepsilon}$:

$$f^{\varepsilon}(t, x, v) = \rho^{\varepsilon}(t, x)\mathcal{M}(v) + \sum_{k=1}^{\infty} \varepsilon^k h^{(k)}(t, x, v). \qquad (2.3)$$

This expansion comes with the following assumptions. First, the functions $h^{(k)}$ do not depend on the parameter $\varepsilon$. Secondly, we assume that $h^{(k)} \in \mathcal{N}^{\perp}$ for all $k$ and is therefore such that:

$$\langle h^{(k)} \rangle = 0, \quad \forall k \ge 1.$$

In particular, we will show that these functions can be expressed using the density $\rho^{\varepsilon}$, the electric field $E$, the velocity variable $v$ and the Maxwellian $\mathcal{M}$.
2.2. **Hierarchical truncations.** To derive a hierarchy of models, let us consider truncations of order $K \in \mathbb{N}$ of the Chapman-Enskog expansion:

$$f^{\epsilon}(t, x, v) \approx \rho^{\epsilon}(t, x)\mathcal{M}(v) + \sum_{k=1}^{K} \epsilon^{k}h^{(k)}(t, x, v). \qquad (2.4)$$

Plugging this expansion in (VBGK) leads to

$$\partial_t (\rho^{\epsilon} \mathcal{M}) + \partial_t \sum_{k=1}^{K} \epsilon^k h^{(k)} = -\frac{1}{\epsilon} \mathcal{T}(\rho^{\epsilon} \mathcal{M}) - \sum_{k=1}^{K} \epsilon^{k-1} \mathcal{T}(h^{(k)}) - \frac{1}{\epsilon} \sum_{k=1}^{K} \epsilon^{k-1} h^{(k)}.$$

Multiplying by $\epsilon$ and rearranging the terms, one obtains

$$\sum_{k=0}^{K-1} \epsilon^k h^{(k+1)} = -\mathcal{T}(\rho^{\epsilon} \mathcal{M}) - \sum_{k=1}^{K} \epsilon^k \mathcal{T}(h^{(k)}) - \partial_t \sum_{k=2}^{K+1} \epsilon^k h^{(k-1)} - \epsilon \partial_t (\rho^{\epsilon} \mathcal{M}).$$

We now identify powers of $\epsilon$:

$$k=0: \quad h^{(1)} = -\mathcal{T}(\rho^{\epsilon} \mathcal{M}), \qquad (2.5a)$$

$$k=1: \quad h^{(2)} = -\partial_t (\rho^{\epsilon} \mathcal{M}) - \mathcal{T}(h^{(1)}), \qquad (2.5b)$$

$$2 \le k \le K-1: \quad h^{(k+1)} = -\partial_t h^{(k-1)} - \mathcal{T}(h^{(k)}). \qquad (2.5c)$$

2.2.1. **Macroscopic model.** To derive the fluid model, let us truncate the Chapman-Enskog expansion at first order $K=1$:

$$f^{\epsilon} \approx \rho^{\epsilon} \mathcal{M} + \epsilon h^{(1)}. \qquad (2.6)$$

We start by integrating (VBGK) in velocity to obtain:

$$\partial_t \rho^{\epsilon} + \frac{1}{\epsilon} \mathrm{div}_x \langle v f^{\epsilon} \rangle = 0.$$

Then $f^{\epsilon}$ is replaced by its expression (2.6):

$$\partial_t \rho^{\epsilon} + \frac{1}{\epsilon} \mathrm{div}_x \left( \rho^{\epsilon} \langle v \mathcal{M} \rangle + \epsilon \langle v h^{(1)} \rangle \right) = 0.$$

The function $h^{(1)}$ is given by the identification (2.5a) and can be simplified using the identity $\nabla_v \mathcal{M} = -v \mathcal{M}$:

$$h^{(1)} = -\mathcal{T}(\rho^{\epsilon} \mathcal{M}) = -v \mathcal{M} \cdot J^{\epsilon} \quad \text{where} \quad J^{\epsilon} = \nabla_x \rho^{\epsilon} - E \rho^{\epsilon}.$$

Using the fact that $\rho^{\epsilon}$ does not depend on the velocity and that odd moments of the Maxwellian are zero, we obtain by plugging in the expression of $h^{(1)}$:

$$\partial_t \rho^{\epsilon} - \langle v \otimes v \mathcal{M} \rangle : \nabla_x J^{\epsilon} = 0, \qquad (2.7)$$

where "$\otimes$" denotes the tensor product and ":" is the tensor contraction of order 2. The moment tensor is then computed:

$$\langle v \otimes v \mathcal{M} \rangle = m_2 I$$

with $I$ the identity matrix. Finally, assuming that $\rho^{\epsilon} \to \rho$ as $\epsilon \to 0$ we formally obtain the drift-diffusion model:

$$\partial_t \rho - m_2 \mathrm{div}_x J = 0, \quad \text{where} \quad J = \nabla_x \rho - E \rho. \qquad (DD)$$
2.2.2. Higher order macroscopic model. To derive the third order fluid model, let us place ourselves in the 1D-1D setting to lighten the computations and focus on the various steps involved.

Let us start by truncating the Chapman-Enskog expansion at third order K = 3:

$$f^{\epsilon} \approx \rho^{\epsilon} \mathcal{M} + \epsilon h^{(1)} + \epsilon^{2} h^{(2)} + \epsilon^{3} h^{(3)}. \quad (2.8)$$

We will see in the derivation process that the second order yields no additional information.

Again, we start by integrating (VBGK) in velocity and we replace $f^{\epsilon}$ by its expansion (2.8):

$$\partial_t \rho^{\epsilon} + \frac{1}{\epsilon} \partial_x \langle v \rho^{\epsilon} \mathcal{M} \rangle + \partial_x \langle v h^{(1)} + \epsilon v h^{(2)} + \epsilon^2 v h^{(3)} \rangle = 0. \quad (2.9)$$

At this point, we use the identification (2.5) to compute the perturbations $h^{(1)}$, $h^{(2)}$ and $h^{(3)}$. One obtains

$$h^{(1)} = -v \mathcal{M} J^{\epsilon},$$

$$h^{(2)} = -\mathcal{M} \partial_t \rho^{\epsilon} + v^2 \mathcal{M} \partial_x J^{\epsilon} + (1 - v^2) \mathcal{M} E J^{\epsilon}.$$

We replace $h^{(1)}$ and $h^{(2)}$ by their expressions in $h^{(3)} = -\partial_t h^{(1)} - \mathcal{T}(h^{(2)})$ to obtain:

$$h^{(3)} = 2v\mathcal{M} \partial_t J^{\epsilon} - v^3\mathcal{M} \partial_{xx} J^{\epsilon} - (v\mathcal{M} - v^3\mathcal{M}) \partial_x (EJ^{\epsilon}) \\ - (2v\mathcal{M} - v^3\mathcal{M})E \partial_x J^{\epsilon} + (3v\mathcal{M} - v^3\mathcal{M})E^2 J^{\epsilon}. \quad (2.10)$$

We simplify (2.9) by using the fact that $v\mathcal{M}$ and $v h^{(2)}$ are odd in $v$, obtaining

$$\partial_t \rho^{\epsilon} + \partial_x \langle v h^{(1)} + \epsilon^2 v h^{(3)} \rangle = 0. \quad (2.11)$$

Equation (2.11) shows that the truncation K = 2 gives no further information. We have already shown in (2.7) that

$$\partial_x \langle v h^{(1)} \rangle = -m_2 \partial_x J^{\epsilon}.$$

Therefore, (2.11) gives $\partial_t \rho^{\epsilon} = m_2 \partial_x J^{\epsilon} + \mathcal{O}(\epsilon^2)$. It follows that

$$\partial_t J^{\epsilon} = m_2(\partial_{xx} J^{\epsilon} - E \partial_x J^{\epsilon}) + \mathcal{O}(\epsilon^2).$$

We use this relation to replace the time derivative of $J^{\epsilon}$ in (2.10) which gives:

$$h^{(3)} = 2m_2 v \mathcal{M} \partial_{xx} J^{\epsilon} - 2m_2 v \mathcal{M} E \partial_x J^{\epsilon} - v^3 \mathcal{M} \partial_{xx} J^{\epsilon} - (v \mathcal{M} - v^3 \mathcal{M}) \partial_x (E J^{\epsilon}) \\ - (2v \mathcal{M} - v^3 \mathcal{M}) E \partial_x J^{\epsilon} + (3v \mathcal{M} - v^3 \mathcal{M}) E^2 J^{\epsilon} + \mathcal{O}(\epsilon^2).$$

The choice to replace the time derivative is motivated by the discrete setting. Indeed, we want to avoid the discretization of mixed derivatives to lighten the cost of the macroscopic criterion. Next, the remaining integral is computed:

$$\partial_x \langle v h^{(3)} \rangle = \partial_x \left[ 2m_2^2 \partial_{xx} J^{\epsilon} - 2m_2^2 E \partial_x J^{\epsilon} - m_4 \partial_{xx} J^{\epsilon} - (m_2 - m_4) \partial_x (E J^{\epsilon}) \right. \\ \left. - (2m_2 - m_4) E \partial_x J^{\epsilon} + (3m_2 - m_4) E^2 J^{\epsilon} + \mathcal{O}(\epsilon^2) \right].$$

With our choice of $\mathcal{M}(v)$, we can explicitly compute $m_2 = 1$ and $m_4 = 3$. Therefore one has:

$$\partial_x \langle v h^{(3)} \rangle = \partial_x \left[ 2 \partial_x (E J^{\epsilon}) - E \partial_x J^{\epsilon} - \partial_{xx} J^{\epsilon} + \mathcal{O}(\epsilon^2) \right].$$

Finally, we obtain a higher order model in the drift-diffusion limit.
**Proposition 1.** (formal) Let $f^\epsilon$ be a solution of (VBGK). Assuming that $f^\epsilon$ admits a Chapman-Enskog expansion of order $K=3$, the truncated model up to order 2 in $\epsilon$ is given by a higher order drift-diffusion equation. The macroscopic density $\rho^\epsilon = \langle f^\epsilon \rangle$ is a solution to:

$$ \partial_t \rho^\epsilon - \partial_x J^\epsilon + \epsilon^2 \partial_x (2\partial_x(EJ^\epsilon) - E\partial_x J^\epsilon - \partial_{xx}J^\epsilon) = \mathcal{O}(\epsilon^4), \quad (\widetilde{DD}) $$

where $J^\epsilon = \partial_x \rho^\epsilon - E\rho^\epsilon$.

**Remark 1.** $(\widetilde{DD})$ is, as expected, a second order correction of (DD).

### 3. MICRO-MACRO MODEL

In this part, we recall the derivation of the micro-macro model for (VBGK). Then, we introduce a micro-macro finite volume scheme that enjoys the property of being Asymptotic Preserving which is a crucial point of the hybrid method we want to construct. The micro-macro approach was first used to derive AP schemes for the radiative heat transfer in [30]. It was then applied to the Boltzmann equation in [4, 33] and to the Vlasov-Poisson-BGK equation in [9].

#### 3.1. Continuous setting.
Let us decompose the distribution $f^\epsilon$ as follows:

$$ f^\epsilon = \rho^\epsilon \mathcal{M} + g^\epsilon. \quad (3.1) $$

We introduce the orthogonal projector $\Pi$ in $L^2(dx d\gamma)$ on $\mathcal{N}$ defined for all $f \in L^2(dx d\gamma)$ by:

$$ \Pi f = \langle f \rangle \mathcal{M}. $$

To help us in the derivation of the micro-macro model, let us first recall the following lemma

**Lemma 1.** Let $f^\epsilon = \rho^\epsilon \mathcal{M} + g^\epsilon$ be a solution of (VBGK). One has:

$$ \Pi(g^\epsilon) = \Pi(\partial_t g^\epsilon) = \Pi(\mathcal{T}(\rho^\epsilon \mathcal{M})) = (I - \Pi)(\partial_t (\rho^\epsilon \mathcal{M})) = 0. $$

Moreover, one has $\Pi(\mathcal{T}g^\epsilon) = \text{div}_x \langle v g^\epsilon \rangle \mathcal{M}$.

To derive the micro-macro model, we start by injecting (3.1) in (VBGK):

$$ \partial_t (\rho^\epsilon \mathcal{M}) + \partial_t g^\epsilon + \frac{1}{\epsilon}(\mathcal{T}(\rho^\epsilon \mathcal{M}) + \mathcal{T}g^\epsilon) = \frac{-1}{\epsilon^2} g^\epsilon. \quad (3.2) $$

We then apply $(I - \Pi)$ to (3.2), and simplify using Lemma 1 to obtain the micro part of the model:

$$ \partial_t g^\epsilon + \frac{1}{\epsilon}(\mathcal{T}g^\epsilon - \Pi(\mathcal{T}g^\epsilon)) + \frac{1}{\epsilon}v\mathcal{M} \cdot J^\epsilon = \frac{-1}{\epsilon^2}g^\epsilon. \quad (3.3) $$

The macro part is obtained by applying $\Pi$ to (3.2) and using Lemma 1. It leads to:

$$ \partial_t \rho^\epsilon \mathcal{M} + \frac{1}{\epsilon} \Pi(\mathcal{T}g^\epsilon) = 0. \quad (3.4) $$

Finally, the micro-macro model is given by:

$$ \partial_t g^\epsilon + \frac{1}{\epsilon}(\mathcal{T}g^\epsilon - \text{div}_x \langle v g^\epsilon \rangle \mathcal{M} + v\mathcal{M} \cdot J^\epsilon) = \frac{-1}{\epsilon^2} g^\epsilon, \quad (\text{Micro}) $$

$$ \partial_t \rho^\epsilon + \frac{1}{\epsilon}\text{div}_x \langle v g^\epsilon \rangle = 0. \quad (\text{Macro}) $$

The following proposition states the equivalence between the (Micro)-(Macro) model and the original equation (VBGK) [9].
<figure><img src="image_3.png" /><figcaption>FIGURE 2. Discretization of the velocity domain.</figcaption></figure>

**Proposition 2. (formal)**

(1) If $f^\epsilon$ is a solution to (VBGK) with an initial data in $L^2(\text{d}x \, \text{d}\gamma)$, then $(\rho^\epsilon, g^\epsilon) = (\langle f^\epsilon \rangle, f^\epsilon - \langle f^\epsilon \rangle \mathcal{M})$ is a solution to (Micro)-(Macro) with the associated initial data

$$ \rho_0 = \langle f_0 \rangle, \quad g_0 = f_0 - \rho_0 \mathcal{M}. $$

(2) Conversely, if $(\rho^\epsilon, g^\epsilon)$ is a solution to (Micro)-(Macro) with initial data $\rho^\epsilon(t=0) = \rho_0$ and $g^\epsilon(t=0) = g_0$ with $\langle g_0 \rangle = 0$ then $\langle g^\epsilon(t) \rangle = 0$, for all $t > 0$ and $f^\epsilon = \rho^\epsilon \mathcal{M} + g^\epsilon$ is a solution to (VBGK) with initial data $f_0 = \rho_0 \mathcal{M} + g_0$.

3.2. **Discrete setting.** Let us now tackle the discretization of the (Micro)-(Macro) model. We shall adopt a finite volume approach to discretize the phasespace. From now on, we restrict ourselves to the 1D-3D problem, namely one dimension in position and three in velocity. In particular, let $\Omega_x = [0, x_\star]$ with periodic boundary conditions. In this setting, the (Micro)-(Macro) equations and the transport operator $\mathcal{T}$ reduce to

$$ \left\{ \begin{array}{l} \partial_t g^\epsilon + \frac{1}{\epsilon^2} g^\epsilon + \frac{1}{\epsilon} (\mathcal{T} g^\epsilon - \partial_x \langle v_x g^\epsilon \rangle \mathcal{M} + v_x \mathcal{M} J_x^\epsilon) = 0 \\ \partial_t \rho^\epsilon + \frac{1}{\epsilon} \partial_x \langle v_x g^\epsilon \rangle = 0, \end{array} \right. $$

and

$$ \mathcal{T}(g) = v_x \partial_x g + E_x \partial_{v_x} g, $$

where $v = \begin{pmatrix} v_x \\ v_y \\ v_z \end{pmatrix}$, $E = \begin{pmatrix} E_x \\ 0 \\ 0 \end{pmatrix}$ and $J_x^\epsilon = \partial_x \rho^\epsilon - E_x \rho^\epsilon$. In the following, we shall omit the subscript $x$ for $E_x$ and $J_x^\epsilon$.

**The mesh.** The velocity domain is restricted to a bounded symmetric cube $[-v_\star, v_\star]^3$ as it is impractical to implement a numerical scheme on an unbounded domain. We consider a Cartesian mesh of the cube composed of $N_v = 2L$ velocity cells in each direction arranged symmetrically around $v = 0$. Let $\mathcal{J} = \{-L+1, \dots, L\}$ and let us denote by $j = (j_x, j_y, j_z) \in \mathcal{J}^3$ a multi-index. The cells of the velocity mesh are given by

$$ \mathcal{V}_j = (v_{j_x-\frac{1}{2}}, v_{j_x+\frac{1}{2}}) \times (v_{j_y-\frac{1}{2}}, v_{j_y+\frac{1}{2}}) \times (v_{j_z-\frac{1}{2}}, v_{j_z+\frac{1}{2}}), \quad j \in \mathcal{J}^3. $$

Each cell $\mathcal{V}_j$ has a constant volume $\Delta v^3$ and midpoint $v_j$. The velocity mesh is illustrated in one dimension in Figure 2.
In position, because of the periodic boundary conditions, we consider a discretization of the 1-D torus $\mathbb{T}$ into $N_x$ primal cells

$$ \mathcal{X}_i = (x_{i-\frac{1}{2}}, x_{i+\frac{1}{2}}), \quad i \in \mathcal{I} = \mathbb{Z}/N_x\mathbb{Z}, $$

of constant length $\Delta x$ and centers $x_i$. We also define the dual cells

$$ \mathcal{X}_{i+\frac{1}{2}} = (x_i, x_{i+1}), \quad i \in \mathcal{I}, $$

of constant length $\Delta x$ and centers $x_{i+\frac{1}{2}}$. The primal control volumes in the phase space are defined by

$$ K_{ij} = \mathcal{X}_i \times \mathcal{V}_j, \quad \forall (i, j) \in \mathcal{I} \times \mathcal{J}^3, $$

while the dual control volumes are given by

$$ K_{i+\frac{1}{2},j} = \mathcal{X}_{i+\frac{1}{2}} \times \mathcal{V}_j, \quad \forall (i,j) \in \mathcal{I} \times \mathcal{J}^3. $$

Finally, we set a time step $\Delta t > 0$ and we define $t^n = n\Delta t$, $n \in \mathbb{N}$.

**The discrete Maxwellian.** We assume that we are given cell values of the uni-dimensional Maxwellian $(M_l)_{l \in \mathcal{J}}$ satisfying:

$$ \left\{ \begin{array}{ll} M_l > 0, & M_l = M_{-l+1}, \quad \forall l \in \mathcal{J}, \\ M_L = M_{L+1}, & M_{-L} = M_{-L+1}, \\ \sum_{l \in \mathcal{J}} M_l \Delta v = 1. & \end{array} \right. \qquad (3.5) $$

These properties are a discrete version of the continuous ones, namely, the positivity, the parity, and the unit mass. The third assumption is a zero boundary flux on the Maxwellian. For a sufficiently large domain in velocity, it is relevant due to the fast decay of the Gaussian. We now define the discrete multidimensional Maxwellian as:

$$ \mathcal{M}_j = M_{j_x} M_{j_y} M_{j_z}, \qquad (3.6) $$

where $(M_{j_x})_{j_x \in \mathcal{J}}$, $(M_{j_y})_{j_y \in \mathcal{J}}$ and $(M_{j_z})_{j_z \in \mathcal{J}}$ satisfy (3.5) and $j$ is the multi-index $(j_x, j_y, j_z)$.

Let us introduce the discrete integration operator in velocity: for $f = (f_j)_{j \in \mathcal{J}^3}$

$$ \langle f \rangle_{\Delta} = \sum_{j \in \mathcal{J}^3} f_j \Delta v^3. \qquad (3.7) $$

We also introduce the discrete moments of the discrete 1D Maxwellian (3.5):

$$ m_k^{\Delta v} = \sum_{l \in \mathcal{J}} v_l^k M_l \Delta v. \qquad (3.8) $$

**Semi-discretization in the phase-space.** We start by considering a semi-discretization in the phase-space of the (Micro)-(Macro) model. Let $(i, j) \in \mathcal{I} \times \mathcal{J}^3$. We choose to approximate the perturbation $g^\epsilon$ on the dual cells while the density $\rho^\epsilon$ is approximated on the primal mesh:

$$ g_{i+\frac{1}{2},j}^{\epsilon}(t) \approx \frac{1}{\Delta x \Delta v^3} \int_{K_{i+\frac{1}{2},j}} g^{\epsilon}(t, x, v) \, dx \, dv \quad \text{and} \quad \rho_{i}^{\epsilon}(t) \approx \frac{1}{\Delta x} \int_{\mathcal{X}_{i}} \rho^{\epsilon}(t, x) \, dx. \quad (3.9) $$

This choice of staggered meshes will result in a more compact stencil for the asymptotic scheme and is quite standard in the literature [33, 35, 9]. We start by
integrating (Macro) on $\mathcal{X}_i$:

$$ \int_{\mathcal{X}_i} \left[ \partial_t \rho^\varepsilon + \frac{1}{\varepsilon} \partial_x \langle v_x g^\varepsilon \rangle \right] dx = 0. $$

After integrating the space derivative, one then obtains a continuous in time finite volume scheme for (Macro):

$$ \frac{d}{dt} \rho_i^\varepsilon = -\frac{1}{\Delta x} \left( \frac{1}{\varepsilon} \langle \xi g_{i+\frac{1}{2}}^\varepsilon \rangle_\Delta - \frac{1}{\varepsilon} \langle \xi g_{i-\frac{1}{2}}^\varepsilon \rangle_\Delta \right), \quad \text{where } \xi_j = \xi_{(j_x, j_y, j_z)} = v_{j_x} \quad \forall j \in \mathcal{J}^3. \tag{3.10} $$

Next, we deal with the (Micro) equation which is integrated on $K_{i+\frac{1}{2},j}$:

$$ \int_{K_{i+\frac{1}{2},j}} \left( \partial_t g^\varepsilon + \frac{1}{\varepsilon^2} g^\varepsilon \right) dx \, dv + \frac{1}{\varepsilon} \int_{K_{i+\frac{1}{2},j}} [\mathcal{T}g^\varepsilon - \partial_x \langle v_x g^\varepsilon \rangle \mathcal{M} + v_x \mathcal{M} J^\varepsilon] \, dx \, dv = 0. \tag{3.11} $$

One then obtains

$$ \begin{aligned} & -\varepsilon \Delta x \Delta v^3 \left( \frac{d}{dt} g_{i+\frac{1}{2},j}^\varepsilon + \frac{1}{\varepsilon^2} g_{i+\frac{1}{2},j}^\varepsilon \right) \\ & = \underbrace{\int_{K_{i+\frac{1}{2},j}} \mathcal{T} g^\varepsilon \, dx \, dv}_{\text{(A)}} - \underbrace{\int_{K_{i+\frac{1}{2},j}} \partial_x \langle v_x g^\varepsilon \rangle \mathcal{M} \, dx \, dv}_{\text{(B)}} + \underbrace{\int_{K_{i+\frac{1}{2},j}} v_x \mathcal{M} J^\varepsilon \, dx \, dv}_{\text{(C)}}. \end{aligned} \tag{3.12} $$

Using the definition of the transport operator $\mathcal{T}$, one has:

$$ \begin{aligned} \text{(A)} &= \int_{V_j} v_x (g^\varepsilon(t, x_{i+1}, v) - g^\varepsilon(t, x_i, v)) \, dv \\ &\quad + \int_{\mathcal{X}_{i+\frac{1}{2}} \times V_{j_y} \times V_{j_z}} E(g^\varepsilon(t, x, v_{j_x+\frac{1}{2}}, v_y, v_z) - g^\varepsilon(t, x, v_{j_x-\frac{1}{2}}, v_y, v_z)) \, dx \, dv_y \, dv_z, \\ \text{(B)} &= \int_{V_j} \mathcal{M} (\langle v_x g^\varepsilon(t, x_{i+1}, v) \rangle - \langle v_x g^\varepsilon(t, x_i, v) \rangle) \, dv, \\ \text{(C)} &= \int_{V_j} v_x \mathcal{M} \, dv \int_{\mathcal{X}_{i+\frac{1}{2}}} J^\varepsilon \, dx. \end{aligned} $$

We now denote by $(\mathcal{F}_{i,j}^\varepsilon)_{ij}$ an approximation of the microscopic flux in position at interfaces $(x_i)_i$, namely

$$ \mathcal{F}_{i,j}^\varepsilon \approx \int_{V_j} [v_x g^\varepsilon(t, x_i, v) - \mathcal{M} \langle v_x g^\varepsilon(t, x_i, v) \rangle] \, dv. \tag{3.13} $$

We also denote by $(\mathcal{G}_{i+\frac{1}{2},j_x+\frac{1}{2},j_y,j_z}^\varepsilon)_{ij}$ an approximation of the microscopic flux in velocity, namely

$$ \mathcal{G}_{i+\frac{1}{2}, j_x+\frac{1}{2}, j_y, j_z}^\varepsilon \approx \int_{\mathcal{X}_{i+\frac{1}{2}} \times V_{j_y} \times V_{j_z}} E g^\varepsilon(t, x, v_{x,j+\frac{1}{2}}, v_y, v_z) \, dx \, dv_x \, dv_y. \tag{3.14} $$
In the following, to lighten the notation, we will denote by $G_{i+\frac{1}{2},j+\frac{1}{2}}^{\epsilon}$ this flux. Let us now present our choice of numerical fluxes. In position, we choose a first-order upwind approximation for both terms of (3.13):

$$ \mathcal{F}_{i,j}^{\epsilon} = \left( \xi_j^+ g_{i-\frac{1}{2},j}^{\epsilon,n} + \xi_j^- g_{i+\frac{1}{2},j}^{\epsilon} \right) \Delta v^3 - \mathcal{M}_j \left\langle \xi^+ g_{i-\frac{1}{2}}^{\epsilon} + \xi^- g_{i+\frac{1}{2}}^{\epsilon} \right\rangle_{\Delta} \Delta v^3, \quad (3.15) $$

where the notation $r^{\pm} = \frac{r^{\pm}|r|}{2}$ is used. At the boundaries in position, the periodic setting implies

$$ \mathcal{F}_{0,j}^{\epsilon} = \mathcal{F}_{N_x,j}^{\epsilon}. \quad (3.16) $$

In velocity, a first-order upwind approximation is used and since $E$ is given, it is explicitly discretized on the dual mesh, $E_{i+\frac{1}{2}} = E(x_{i+\frac{1}{2}})$. The numerical flux in velocity then reads:

$$ \mathcal{G}_{i+\frac{1}{2},j+\frac{1}{2}}^{\epsilon} = \left( E_{i+\frac{1}{2}}^{+} g_{i+\frac{1}{2},j_x,j_y,j_z}^{\epsilon,n} + E_{i+\frac{1}{2}}^{-} g_{i+\frac{1}{2},j_x+1,j_y,j_z}^{\epsilon} \right) \Delta x \Delta v^2. \quad (3.17) $$

Zero flux boundary conditions are applied in velocity and therefore we set

$$ \mathcal{G}_{i+\frac{1}{2},-L+\frac{1}{2}}^{\epsilon} = \mathcal{G}_{i+\frac{1}{2},L+\frac{1}{2}}^{\epsilon} = 0. \quad (3.18) $$

Finally, (C) is treated as a source term and approximated using first-order centered finite differences:

$$ (C) \approx \xi_j \mathcal{M}_j J_{i+\frac{1}{2}}^{\epsilon} \Delta x \Delta v^3, $$

with $J_{i+\frac{1}{2}}^{\epsilon} = \frac{\rho_{i+1}^{\epsilon} - \rho_i^{\epsilon}}{\Delta x} - E_{i+\frac{1}{2}} \rho_{i+\frac{1}{2}}^{\epsilon}$, and $\rho_{i+\frac{1}{2}}^{\epsilon} = \frac{1}{2}(\rho_i^{\epsilon} + \rho_{i+1}^{\epsilon})$. A continuous in time finite volume scheme for equation (Micro) finally reads:

$$ \frac{d}{dt} g_{i+\frac{1}{2},j}^{\epsilon} + \frac{1}{\epsilon^2} g_{i+\frac{1}{2},j}^{\epsilon} = -\frac{1}{\epsilon} \left( \frac{T_{i+\frac{1}{2},j}^{\epsilon}}{\Delta x \Delta v^3} + \xi_j \mathcal{M}_j J_{i+\frac{1}{2}}^{\epsilon} \right), \quad (3.19) $$

where $T_{i+\frac{1}{2},j}^{\epsilon} = \mathcal{F}_{i+1,j}^{\epsilon} - \mathcal{F}_{i,j}^{\epsilon} + \mathcal{G}_{i+\frac{1}{2},j+\frac{1}{2}}^{\epsilon} - \mathcal{G}_{i+\frac{1}{2},j-\frac{1}{2}}^{\epsilon}$.

**Full discretization.** In order to obtain an AP scheme, one must carefully choose the discretization in time. Following [32], we adapt the so-called relaxed micro-macro scheme to our finite volume setting. This method falls into the framework of exponential time integrators [23]. Let $n \in \mathbb{N}$ and $(i, j) \in \mathcal{I} \times \mathcal{J}^3$. Let $(g_{i+\frac{1}{2},j}^{\epsilon,n})_{ij}$ be an approximation of $(g_{i+\frac{1}{2},j}^{\epsilon}(t^n))_{ij}$ and $(\rho_i^{\epsilon,n})_i$ an approximation of $(\rho_i^{\epsilon}(t^n))_i$. The first step is to multiply (3.19) by $e^{t/\epsilon^2}$ which gives:

$$ \frac{d}{dt} \left( g_{i+\frac{1}{2},j}^{\epsilon}(t)e^{t/\epsilon^2} \right) = -\frac{e^{t/\epsilon^2}}{\epsilon} \left( \frac{T_{i+\frac{1}{2},j}^{\epsilon}}{\Delta x \Delta v^3} + \xi_j \mathcal{M}_j J_{i+\frac{1}{2}}^{\epsilon} \right). \quad (3.20) $$

Let us then integrate between $t^n$ and $t^{n+1}$ and divide by $e^{t^{n+1}/\epsilon^2}$:

$$ g_{i+\frac{1}{2},j}^{\epsilon}(t^{n+1}) = g_{i+\frac{1}{2},j}^{\epsilon}(t^n)e^{-\Delta t/\epsilon^2} + \int_{t^n}^{t^{n+1}} -\frac{e^{(t-t^{n+1})/\epsilon^2}}{\epsilon} \left( \frac{T_{i+\frac{1}{2},j}^{\epsilon}}{\Delta x \Delta v^3} + \xi_j \mathcal{M}_j J_{i+\frac{1}{2}}^{\epsilon} \right) dt. \quad (3.21) $$
Then, the transport and source terms are approximated at time $t^n$ and the integral can be computed explicitly:

$$ 
\begin{aligned}
& \int_{t^n}^{t^{n+1}} -\frac{e^{(t-t^{n+1})/\epsilon^2}}{\epsilon} \left( \frac{T_{i+\frac{1}{2},j}^{\epsilon}}{\Delta x \Delta v^3} + \xi_j \mathcal{M}_j J_{i+\frac{1}{2}}^{\epsilon} \right) dt \\
& = -\epsilon(1 - e^{-\Delta t/\epsilon^2}) \left( \frac{T_{i+\frac{1}{2},j}^{\epsilon,n}}{\Delta x \Delta v^3} + \xi_j \mathcal{M}_j J_{i+\frac{1}{2}}^{\epsilon,n} \right).
\end{aligned} 
\quad (3.22) $$

Finally, the fully discretized microscopic equation reads:

$$ g_{i+\frac{1}{2},j}^{\epsilon,n+1} = g_{i+\frac{1}{2},j}^{\epsilon,n} e^{-\Delta t/\epsilon^2} - \epsilon(1 - e^{-\Delta t/\epsilon^2}) \left( \frac{T_{i+\frac{1}{2},j}^{\epsilon,n}}{\Delta x \Delta v^3} + \xi_j \mathcal{M}_j J_{i+\frac{1}{2}}^{\epsilon,n} \right), \quad (3.23) $$

where

$$ T_{i+\frac{1}{2},j}^{\epsilon,n} = \mathcal{F}_{i+1,j}^{\epsilon,n} - \mathcal{F}_{i,j}^{\epsilon,n} + \mathcal{G}_{i+\frac{1}{2},j+\frac{1}{2}}^{\epsilon,n} - \mathcal{G}_{i+\frac{1}{2},j-\frac{1}{2}}^{\epsilon,n}. \quad (3.24) $$

The discretization in time of (3.10) is quite standard, with an implicit discretization of the stiff term:

$$ \rho_i^{\epsilon,n+1} = \rho_i^{\epsilon,n} - \frac{\Delta t}{\epsilon \Delta x} \left( \langle \xi g_{i+\frac{1}{2}}^{\epsilon,n+1} \rangle_{\Delta} - \langle \xi g_{i-\frac{1}{2}}^{\epsilon,n+1} \rangle_{\Delta} \right) \quad (3.25) $$

Note that (3.23) defines an explicit scheme. Moreover, (3.25) does not require the inversion of a system. Indeed, (3.23) is explicitly computed at time $t^{n+1}$ and is then used to update the density in (3.25). In practice, the method is therefore fully explicit.

Before stating the next proposition, let us introduce the following assumption:

**Assumption 1.** Let $(\rho_i^{\epsilon,n})_{i \in \mathcal{I}}$ be given by (3.25). Then,

$$ \rho_i^{\epsilon,n} \rightarrow \rho_i^n, \quad \forall i \in \mathcal{I}. \quad (3.26) $$

This assumption corresponds to the convergence of $\rho^{\epsilon}$ to $\rho$ as $\epsilon \to 0$. Such property is not trivial to obtain in the discrete setting. A rigorous proof of this result requires, among other things, uniform estimates in $\epsilon$ of the discrete $L^2$-norm of $\rho^{\epsilon}$, $g^{\epsilon}$, and moments of $g^{\epsilon}$. It is outside the scope of this article and may be thoroughly investigated in upcoming work.

The following proposition states the AP property of our discretization of the (Micro)-(Macro) model.

**Proposition 3.** Let $n \in \mathbb{N}$. Let $(g_{i+\frac{1}{2},j}^{\epsilon,n})_{ij}$ and $(\rho_i^{\epsilon,n})_i$ be given by the following micro-macro finite volume scheme:

$$ g_{i+\frac{1}{2},j}^{\epsilon,n+1} = g_{i+\frac{1}{2},j}^{\epsilon,n} e^{-\Delta t/\epsilon^2} - \epsilon(1 - e^{-\Delta t/\epsilon^2}) \left( \frac{T_{i+\frac{1}{2},j}^{\epsilon,n}}{\Delta x \Delta v^3} + \xi_j \mathcal{M}_j J_{i+\frac{1}{2}}^{\epsilon,n} \right), \quad (S_{\text{Micro}}) $$

$$ \rho_i^{\epsilon,n+1} = \rho_i^{\epsilon,n} - \frac{\Delta t}{\epsilon \Delta x} \left( \langle \xi g_{i+\frac{1}{2}}^{\epsilon,n+1} \rangle_{\Delta} - \langle \xi g_{i-\frac{1}{2}}^{\epsilon,n+1} \rangle_{\Delta} \right), \quad (S_{\text{Macro}}) $$

where $T_{i+\frac{1}{2},j}^{\epsilon,n}$ is given by (3.24).

Assuming that Assumption 1 holds and for a fixed mesh size $\Delta x, \Delta v > 0$, the scheme enjoys the AP property in the diffusion limit. This property does not depend on the initial data, and the associated limit scheme reads

$$ \rho_i^{n+1} = \rho_i^n + m_2^{\Delta v} \frac{\Delta t}{\Delta x} (J_{i+\frac{1}{2}}^n - J_{i-\frac{1}{2}}^n), \quad (S_{\text{Lim}}) $$
with the limit flux

$$J_{i+\frac{1}{2}}^{n} = \frac{\rho_{i+1}^{n} - \rho_{i}^{n}}{\Delta x} - E_{i+\frac{1}{2}} \rho_{i+\frac{1}{2}}^{n}, \quad (3.27)$$

where $m_{\frac{1}{2}}^{\Delta v}$ is given by (3.8).

*Proof.* The mesh size $\Delta x$, $\Delta v > 0$ being set, let us emphasize that we consider only the pointwise convergence of the scheme as $\varepsilon$ tends to 0. The first step is to study the asymptotic behaviour of the perturbation $(g_{i+\frac{1}{2},j}^{\varepsilon,n+1})_{i,j}$. By induction on $n$, let us show that $g_{i+\frac{1}{2},j}^{\varepsilon,n+1} \xrightarrow{\varepsilon \to 0} 0$ for any initial data $(\rho^0, g^0)$ and for all $(i, j) \in \mathcal{I} \times \mathcal{J}^3$. At $n=0$, one has

$$g_{i+\frac{1}{2},j}^{\varepsilon,1} = g_{i+\frac{1}{2},j}^{0} e^{-\Delta t/\varepsilon^2} - \varepsilon(1 - e^{-\Delta t/\varepsilon^2}) \left( \frac{T_{i+\frac{1}{2},j}^{\varepsilon,0}}{\Delta x \Delta v^3} + \xi_j \mathcal{M}_j J_{i+\frac{1}{2}}^{\varepsilon,0} \right). \quad (3.28)$$

As $e^{-\Delta t/\varepsilon^2} \xrightarrow{\varepsilon \to 0} 0$ and since $(\frac{T_{i+\frac{1}{2},j}^{\varepsilon,0}}{\Delta x \Delta v^3} + \xi_j \mathcal{M}_j J_{i+\frac{1}{2}}^{\varepsilon,0})_{ij}$ depends only on the initial data which itself is independent of $\varepsilon$,

$$g_{i+\frac{1}{2},j}^{\varepsilon,1} \xrightarrow{\varepsilon \to 0} 0 \quad \forall (i,j) \in \mathcal{I} \times \mathcal{J}^3.$$

Let us now assume that

$$g_{i+\frac{1}{2},j}^{\varepsilon,n} \xrightarrow{\varepsilon \to 0} 0 \quad \forall (i,j) \in \mathcal{I} \times \mathcal{J}^3. \quad (3.29)$$

Under the hypothesis (3.29), one obtains that the transport term, that depends only on the perturbation, vanishes in the limit:

$$T_{i+\frac{1}{2},j}^{\varepsilon,n} \xrightarrow{\varepsilon \to 0} 0 \quad \forall (i,j) \in \mathcal{I} \times \mathcal{J}^3.$$

From Assumption 1 one obtains the convergence of the source term:

$$\xi_j \mathcal{M}_j J_{i+\frac{1}{2}}^{\varepsilon,n} \xrightarrow{\varepsilon \to 0} \xi_j \mathcal{M}_j J_{i+\frac{1}{2}}^{n}. \quad (3.30)$$

Now, the asymptotic limit of $(S_{Micro})$ can be computed. Since $\varepsilon(1 - e^{-\Delta t/\varepsilon^2}) \xrightarrow{\varepsilon \to 0} 0$ and $e^{-\Delta t/\varepsilon^2} \xrightarrow{\varepsilon \to 0} 0$ we can use (3.30) and our induction hypothesis (3.29) to obtain:

$$g_{i+\frac{1}{2},j}^{\varepsilon,n+1} \xrightarrow{\varepsilon \to 0} 0, \quad \forall n, \quad \forall (i,j) \in \mathcal{I} \times \mathcal{J}. \quad (3.31)$$

As a consequence, one also has that for all $n$,

$$\frac{T_{i+\frac{1}{2},j}^{\varepsilon,n}}{\Delta x \Delta v^3} + \xi_j \mathcal{M}_j J_{i+\frac{1}{2}}^{\varepsilon,n} \xrightarrow{\varepsilon \to 0} \xi_j \mathcal{M}_j J_{i+\frac{1}{2}}^{n}, \quad \forall (i,j) \in \mathcal{I} \times \mathcal{J}. \quad (3.32)$$

The next step is to plug $(S_{Micro})$ into $(S_{Macro})$:

$$\rho_i^{\varepsilon,n+1} = \rho_i^{\varepsilon,n} - \frac{\Delta t}{\varepsilon \Delta x} \left\langle \xi \left[ e^{-\Delta t/\varepsilon^2} (g_{i+\frac{1}{2}}^{\varepsilon,n} - g_{i-\frac{1}{2}}^{\varepsilon,n}) - \varepsilon(1 - e^{-\Delta t/\varepsilon^2}) \left( \frac{T_{i+\frac{1}{2}}^{\varepsilon,n} - T_{i-\frac{1}{2}}^{\varepsilon,n}}{\Delta x \Delta v^3} + (\xi \mathcal{M}_j J_{i+\frac{1}{2}}^{\varepsilon,n} - \xi \mathcal{M}_j J_{i-\frac{1}{2}}^{\varepsilon,n}) \right) \right] \right\rangle_{\Delta} \quad (3.33)$$
We can then take the limit $\epsilon \to 0$ in (3.33) using the previous asymptotic limits (3.31) and (3.32):

$$ \rho_i^{n+1} = \rho_i^n + \langle \xi^2 \mathcal{M} \rangle_{\Delta} \frac{\Delta t}{\Delta x} \left( J_{i+\frac{1}{2}}^n - J_{i-\frac{1}{2}}^n \right). \quad (3.34) $$

Then, using the definitions of the discrete Maxwellians (3.5), (3.6) one can obtain:

$$ \begin{aligned} \langle \xi^2 \mathcal{M} \rangle_{\Delta} &= \sum_{j \in \mathcal{J}^3} \xi_j^2 \mathcal{M}_j \Delta v^3 \\ &= \sum_{j_x \in \mathcal{J}} \sum_{j_y \in \mathcal{J}} \sum_{j_z \in \mathcal{J}} v_{j_x}^2 M_{j_x} M_{j_y} M_{j_z} \Delta v^3 \\ &= \sum_{j_x \in \mathcal{J}} v_{j_x}^2 M_{j_x} \Delta v \\ &= m_2^{\Delta v}. \end{aligned} \quad (3.35) $$

Finally, we obtain the asymptotic scheme (S$_{Lim}$):

$$ \rho_i^{n+1} = \rho_i^n + m_2^{\Delta v} \frac{\Delta t}{\Delta x} \left( J_{i+\frac{1}{2}}^n - J_{i-\frac{1}{2}}^n \right), $$

which concludes the proof. □

In order to show that the scheme (S$_{Micro}$)-(S$_{Macro}$) is truly AP, one also needs the stability condition to be independent (or at least does not degenerate) on $\epsilon$. While we do not prove the stability of the scheme, in practice, we can indeed use the same time-step for both large and small values of $\epsilon$.

## 4. HYBRID METHOD

The aim of this section is to introduce a hybrid method between kinetic and fluid schemes. The goal is to obtain a coupled solver that is faster than a full kinetic one to solve ($P^\epsilon$) while still being accurate. These methods come naturally when designing accurate numerical codes while guarantying reasonable computation times.

Following [19] we first construct a hybrid kinetic/fluid solver with a dynamic domain adaptation method and present its implementation. In the second part, we are interested in understanding the conservative aspect of the method. More precisely, we give a result on the mass variation induced by the coupling.

### 4.1. Coupling criteria.
The idea of the dynamic domain adaptation method is twofold. First, the subdomains must accurately describe the state of the solution. In particular, the fluid model is only valid where the solution is near the local equilibrium in velocity. Secondly, we want the method to be dynamic in the sense that the subdomains are adapted at each time step. For this purpose, let us introduce $\Omega_K^n$ the kinetic domain and $\Omega_F^n$ the fluid one at time $t^n$. To determine in which domain each cell lies, we introduce criteria based on the higher order fluid model introduced in Section 2.2.2 and the norm of the perturbation $g^\epsilon = f^\epsilon - \rho^\epsilon \mathcal{M}$. Indeed, when $g^\epsilon$ is close to 0, it means that the solution is close to the local equilibrium in velocity.

Let us consider a fluid subdomain. In this subdomain, one only has access to the macroscopic quantity $\rho$ and the given electrical field $E$. Therefore, one cannot consider the perturbation $g^\epsilon$. The solution we propose is to use the higher order
<table><caption>TABLE 1. Central finite differences coefficients.</caption><thead><tr><th>Derivative</th><th>-3</th><th>-2</th><th>-1</th><th>0</th><th>1</th><th>2</th><th>3</th></tr></thead><tbody><tr><th>1</th><td>-1/60</td><td>3/20</td><td>-3/4</td><td>0</td><td>3/4</td><td>-3/20</td><td>1/60</td></tr><tr><th>2</th><td>1/90</td><td>-3/20</td><td>3/2</td><td>-49/18</td><td>3/2</td><td>-3/20</td><td>1/90</td></tr><tr><th>3</th><td>1/8</td><td>-1</td><td>13/8</td><td>0</td><td>-13/8</td><td>1</td><td>-1/8</td></tr><tr><th>4</th><td>-1/6</td><td>2</td><td>-13/2</td><td>28/3</td><td>-13/2</td><td>2</td><td>-1/6</td></tr></tbody></table>

model (DD̃) and derive a macroscopic criterion. We have formally shown that (DD̃) can be written in the form:

$$ \partial_t \rho^\epsilon - m_2 \operatorname{div}_x J^\epsilon = \mathcal{R}^\epsilon, $$

where $\mathcal{R}^\epsilon$ is a remainder that depends only on the density $\rho^\epsilon$ and the electrical field $E$. During the coupling procedure it will be computed using both the kinetic density $\rho^\epsilon$ in kinetic cells and using the fluid density $\rho$ in fluid cells. In the 1D-3D setting, it is given by

$$ \mathcal{R}^{\epsilon} = -\epsilon^2 \partial_x (2 \partial_x (E J^{\epsilon}) - E \partial_x J^{\epsilon} - \partial_{xx} J^{\epsilon}), \quad \text{where } J^{\epsilon} = \partial_x \rho^{\epsilon} - E \rho^{\epsilon}. \quad (4.1) $$

Expanding $\mathcal{R}^\epsilon$ shows that one needs derivatives of the density up to fourth order and of $E$ up to third order:

$$ \begin{aligned} \mathcal{R}^{\epsilon} = & -\epsilon^2 \left( -\partial_{xxxx} \rho^{\epsilon} \right. \\ & \qquad \left. + E(2 \partial_{xxx} \rho^{\epsilon} - E \partial_{xx} \rho^{\epsilon}) \right. \\ & \qquad \left. + \partial_x E(-3\rho^{\epsilon} \partial_x E - 5E \partial_x \rho^{\epsilon} + 6 \partial_{xx} \rho^{\epsilon}) \right) \\ & + \partial_{xx} E(-3\rho^{\epsilon}E + 5 \partial_x \rho^{\epsilon}) \\ & \qquad \left. + \rho^{\epsilon} \partial_{xxx} E \right). \end{aligned} \quad (4.2) $$

Let us denote by $\mathcal{R}_i^{\epsilon,n}$ a discretization of the remainder $\mathcal{R}^\epsilon$ at time $t^n$ in cell $X_i$. High order finite difference schemes are used (See Table 1).

Let $\eta_0, \delta_0 > 0$ be the coupling thresholds. In a fluid domain, when $\mathcal{R}^{\epsilon,n}$ is large, the model (DD̃) is far from the limit model (DD) and one must use the kinetic one instead. More specifically, consider a fluid cell $X_i \subset \Omega_F^n$.

* If $|\mathcal{R}_i^{\epsilon,n}| \le \eta_0$, then the cell stays fluid at $t^{n+1}$.

* If $|\mathcal{R}_i^{\epsilon,n}| > \eta_0$, then the cell becomes kinetic at $t^{n+1}$:

$$ X_i \not\subset \Omega_F^{n+1} \quad \text{and} \quad X_i \subset \Omega_K^{n+1}. $$

In a kinetic subdomain, unlike the previous case, one has access to the perturbation $g^\epsilon$. When this perturbation is small, it means that the solution is near a local equilibrium in velocity. As a consequence, the behaviour of the system is close to the fluid one and one can use the limit model instead. Moreover, we also use the criterion that the remainder $\mathcal{R}^\epsilon$ must be small. Consider now a kinetic cell $X_i \subset \Omega_K^n$:

* If $||g_{i-1/2}^{\epsilon,n}||_2 > \delta_0$ and $||g_{i+1/2}^{\epsilon,n}||_2 > \delta_0$ then the cell stays kinetic at $t^{n+1}$.

* If $||g_{i-1/2}^{\epsilon,n}||_2 \le \delta_0$, $||g_{i+1/2}^{\epsilon,n}||_2 \le \delta_0$ and $|\mathcal{R}_i^{\epsilon,n}| > \eta_0$, then the cell stays kinetic at $t^{n+1}$.
* If $||g_{i-\frac{1}{2}}^{\epsilon,n}||_2 \le \delta_0$, $||g_{i+\frac{1}{2}}^{\epsilon,n}||_2 \le \delta_0$ and $|\mathcal{R}_i^{\epsilon,n}| \le \eta_0$, then the cell becomes fluid at $t^{n+1}$:

$$ \mathcal{X}_i \not\subset \Omega_K^{n+1} \quad \text{and} \quad \mathcal{X}_i \subset \Omega_F^{n+1}. $$

The discrete norm $||g_{i+\frac{1}{2}}^{\epsilon,n}||_2$ is the classic $l^2$-norm. It is defined for $(g_{i+\frac{1}{2},j}^{\epsilon,n})_{j \in \mathcal{J}^3}$ by:

$$ ||g_{i+\frac{1}{2}}^{\epsilon,n}||_2 = \sum_{j \in \mathcal{J}^3} (g_{i+\frac{1}{2},j}^{\epsilon,n})^2 \Delta v^3. \quad (4.3) $$

**Remark 2.** Note that in a kinetic cell, the criterion on the norm of $g^\epsilon$ is mandatory. Indeed, the remainder $\mathcal{R}_i^{\epsilon,n}$ could be small because of small gradients, but the perturbation large. In this situation, one does not want to change from kinetic to fluid. As an example, one could take a distribution function at equilibrium in position and far from the Maxwellian in velocity.

### 4.2. Implementation.
We now present in more details the implementation of the hybrid method. An important part of this approach is the management of boundary conditions. When solving on the whole space domain, periodic boundary conditions are applied. However, when solving in the subdomains $\Omega_K^n$ and $\Omega_F^n$, we need to adapt our solvers. Our strategy is to use ghost cell values that are chosen appropriately. The difficulty lies in the fact that the limit scheme only computes the density $\rho$ and not the pair $(\rho^\epsilon, g^\epsilon)$. Since the hybrid method is dynamic, one does not know in advance the state of the cells. As a consequence, one must be able to access all unknowns on the whole domain at any time. Our solution is to take advantage of the structure of the micro-macro scheme. Indeed, aside from visualisation and diagnostics, an explicit discretization of the distribution function isn't necessary. We are working only with $\rho^{\epsilon,n}$ and $g^{\epsilon,n}$. Therefore, we have access to the macro unknown on the whole domain and there is no information missing in the arrays. The distribution $f^\epsilon$ is reconstructed using $f_{i,j}^{\epsilon,n} = \rho_i^{\epsilon,n} M_j + \frac{1}{2} (g_{i-\frac{1}{2},j}^{\epsilon,n} + g_{i+\frac{1}{2},j}^{\epsilon,n})$ only for posttreatment. However, the kinetic solver may still need values of $g^\epsilon$ on the whole space domain. In theory, the array storing $g^\epsilon$ must be filled with zeros in the fluid domain. However, to improve the performance, $g^\epsilon$ is in practice not updated in fluid regions and it is set to 0 only when needed. In particular, it only occurs when a fluid cell becomes kinetic and when saving datas.

Another important remark is that since $g^\epsilon$ is approximated on the dual mesh, one must be careful at the interfaces between kinetic and fluid subdomains. To avoid any ambiguity on the state of an interface when updating the perturbation $g^\epsilon$, we impose that a fluid subdomain is at least two cells wide. Under this condition the state of the ghost interface is well determined. See Figure 3 for an illustration of such a situation.

The algorithm can be summarized as follows:
<figure><img src="image_4.png" /><figcaption>FIGURE 3. Transition between kinetic and fluid cell for the micro-macro scheme.</figcaption></figure>

**Algorithm 1** Hybrid scheme

(1) Set $\varepsilon, \delta_0, \eta_0$ and a final time $T$.

(2) Initialize micro-macro unknowns using the relations $\rho^0 = \langle f^0 \rangle$ and $g^0 = f^0 - \rho^0 M$.

(3) Initialize $\Omega_K^0$ as the whole space ($\Omega_K^0 = \emptyset$).

(4) Compute $g^{\varepsilon,n+1}$ and $\rho^{\varepsilon,n+1}$ in $\Omega_K^n$ using the kinetic scheme ($S_{Micro}$)-($S_{Macro}$).

(5) Compute $\rho^{n+1}$ in $\Omega_F^n$ using the limit scheme ($S_{Lim}$).

(6) Set $g^{n+1} = 0$ in $\Omega_F^n$.

(7) Update $\Omega_K^n$ and $\Omega_F^n$ to $\Omega_K^{n+1}$ and $\Omega_F^{n+1}$ using the criteria presented above.

(8) Increment time and repeat until $t^{n+1} = T$.

In particular, Algorithm 1 explicitly defines a numerical scheme on the hybrid density $\tilde{\rho}$:

$$ \tilde{\rho}_i^{n+1} = \tilde{\rho}_i^n + \frac{\Delta t}{\Delta x} J_i^{H,n}, \quad (4.4) $$

where

$$ J_i^{H,n} = \begin{cases} -\frac{1}{\varepsilon} \left( \langle \xi g_{i+\frac{1}{2}}^{\varepsilon,n+1} \rangle_{\Delta} - \langle \xi g_{i-\frac{1}{2}}^{\varepsilon,n+1} \rangle_{\Delta} \right) & \text{if } \chi_i \in \Omega_K^n, \\ m_2^{\Delta v} \left( J_{i+\frac{1}{2}}^n - J_{i-\frac{1}{2}}^n \right) & \text{if } \chi_i \in \Omega_F^n. \end{cases} \quad (4.5) $$

Note that we want to start the resolution with the approach containing the full information on the system. Hence, it makes sense to initialize our domain as fully kinetic. Moreover, let us emphasize again that the kinetic fluxes are in practice explicitly computed. Therefore the hybrid setting remains an explicit method.

### 4.3. Mass conservation.
This section is dedicated to investigate the mass conservation of the hybrid method. This property being satisfied in the continuous case, one expects conservation in the discrete setting. Each of the standard schemes is conservative on its own by construction. However, the question arises when considering the hybrid scheme.

In order to understand the variation of mass, we consider a toy model that isn't relevant in practice but will highlight the key elements to constrain the mass variation. Let us set the state of cells for every time step into two domains:

$$ \Omega_K = \bigcup_{i=1}^{s-1} \chi_i \quad \text{and} \quad \Omega_F = \bigcup_{i=s}^{Nx} \chi_i. \quad (4.6) $$
<figure><img src="image_5.png" /><figcaption>FIGURE 4. Zoom on the interface of a steady domain adaptation.</figcaption></figure>

Note that in the next result, we neglect what happens at the boundary. Our primary focus is to understand what happens at the interface $x_{s-\frac{1}{2}}$ between the two domains. Moreover, in that context and with periodic boundary conditions in position, the same analysis can be done at the interface $x_{\frac{1}{2}}$. Figure 4 illustrates this framework. We define the mass of the hybrid system as:

$$m^n = \sum_{i \in \mathcal{I}} \tilde{\rho}_i^n \Delta x. \qquad (4.7)$$

Note that this definition suggests that the quantity $\langle g_{i+\frac{1}{2}}^n \rangle_\Delta$ is zero. This property holds in the continuous setting and will be numerically investigated in the next section. The following lemma quantifies the mass variation between two time steps.

**Lemma 2.** Let $(\tilde{\rho}_i^n)_i$ and $(g_{i+\frac{1}{2},j}^{\varepsilon,n})_{ij}$ be computed using the hybrid scheme $(S_{Micro})$- (4.4). Let the mass variation between $t^n$ and $t^{n+1}$ be defined as:

$$\Delta m^{n+\frac{1}{2}} = \sum_{i \in \mathcal{I}} \Delta x \frac{(\tilde{\rho}_i^{n+1} - \tilde{\rho}_i^n)}{\Delta t}. \qquad (4.8)$$

In the context of the steady domain adaptation (4.6) and neglecting the boundaries, one has:

$$\Delta m^{n+\frac{1}{2}} = -\left\langle \xi g_{s-\frac{1}{2}}^{\varepsilon,n} \right\rangle_{\Delta} \frac{e^{-\Delta t/\varepsilon^2}}{\varepsilon} + \frac{1 - e^{-\Delta t/\varepsilon^2}}{\Delta x \Delta v^3} \left\langle \xi T_{s-\frac{1}{2}}^{\varepsilon,n} \right\rangle_{\Delta} \\ - m_2^{\Delta v} e^{-\Delta t/\varepsilon^2} J_{s-\frac{1}{2}}^{\varepsilon,n} + m_2^{\Delta v} (J_{s-\frac{1}{2}}^{\varepsilon,n} - J_{s-\frac{1}{2}}^n). \qquad (4.9)$$

An important consequence of this lemma is that thanks to (3.31) and Assumption 1, the mass variation converges to 0 as $\varepsilon$ tends to 0.
*Proof.* Let us consider the hybrid scheme (4.4)-(4.5) on the density. Using the fixed domain adaptation (4.6) and neglecting the boundary, the mass variation writes:

$$ 
\begin{aligned}
\Delta m^{n+\frac{1}{2}} &= \sum_{i \in \mathcal{I}} \frac{\Delta x}{\Delta t} (\tilde{\rho}_i^{n+1} - \tilde{\rho}_i^n) \\
&= \sum_{i \in \mathcal{I}} J_i^{H,n} \\
&= -\frac{1}{\varepsilon} \sum_{i=1}^{s-1} \left( \langle \xi g_{i+\frac{1}{2}}^{\varepsilon,n+1} \rangle_{\Delta} - \langle \xi g_{i-\frac{1}{2}}^{\varepsilon,n+1} \rangle_{\Delta} \right) + m_2^{\Delta v} \sum_{i=s}^{N_x} \left( J_{i+\frac{1}{2}}^n - J_{i-\frac{1}{2}}^n \right) \\
&= -\frac{1}{\varepsilon} \left\langle \xi g_{s-\frac{1}{2}}^{\varepsilon,n+1} \right\rangle_{\Delta} - m_2^{\Delta v} J_{s-\frac{1}{2}}^n.
\end{aligned} 
\quad (4.10) $$

Similarly as in the proof of Proposition 3, $g_{s-\frac{1}{2},j}^{\varepsilon,n+1}$ is replaced by its expression $(S_{Micro})$. The quantity $\frac{1}{\varepsilon} \langle \xi g_{s-\frac{1}{2}}^{\varepsilon,n+1} \rangle_{\Delta}$ then reads:

$$ \frac{1}{\varepsilon} \langle \xi g_{s-\frac{1}{2}}^{\varepsilon,n+1} \rangle_{\Delta} = \langle \xi g_{s-\frac{1}{2}}^{\varepsilon,n} \rangle_{\Delta} \frac{e^{-\Delta t / \varepsilon^2}}{\varepsilon} - \langle \xi T_{s-\frac{1}{2}}^{\varepsilon,n} \rangle_{\Delta} (1 - e^{-\Delta t / \varepsilon^2}) \frac{1}{\Delta x \Delta v^3} \\ - \langle \xi^2 \mathcal{M}_j \rangle_{\Delta} J_{s-\frac{1}{2}}^{\varepsilon,n} (1 - e^{-\Delta t / \varepsilon^2}). \quad (4.11) $$

Finally, plugging (4.11) into (4.10) and using the computation (3.35) for the term $\langle \xi^2 \mathcal{M}_j \rangle_{\Delta}$, one obtains:

$$ \Delta m^{n+\frac{1}{2}} = - \langle \xi g_{s-\frac{1}{2}}^{\varepsilon,n} \rangle_{\Delta} \frac{e^{-\Delta t / \varepsilon^2}}{\varepsilon} + \frac{1 - e^{-\Delta t / \varepsilon^2}}{\Delta x \Delta v^3} \langle \xi T_{s-\frac{1}{2}}^{\varepsilon,n} \rangle_{\Delta} \\ - m_2^{\Delta v} e^{-\Delta t / \varepsilon^2} J_{s-\frac{1}{2}}^{\varepsilon,n} + m_2^{\Delta v} (J_{s-\frac{1}{2}}^{\varepsilon,n} - J_{s-\frac{1}{2}}^n). \quad \square $$

**Remark 3.** The proof only holds in the context of the toy problem (4.6). However, it can be extended to a more general setting seeing that the mass variation occurs at all interfaces between kinetic and fluid subdomains. Namely,

$$ \Delta t \Delta m^{n+\frac{1}{2}} = \sum_{\alpha \in S} \beta \left( - \langle \xi g_{\alpha}^{\varepsilon,n} \rangle_{\Delta} \frac{e^{-\Delta t / \varepsilon^2}}{\varepsilon} + \frac{1 - e^{-\Delta t / \varepsilon^2}}{\Delta x \Delta v^3} \langle \xi T_{\alpha}^{\varepsilon,n} \rangle_{\Delta} \right. \\ \left. - m_2^{\Delta v} e^{-\Delta t / \varepsilon^2} J_{\alpha}^{\varepsilon,n} + m_2^{\Delta v} (J_{\alpha}^{\varepsilon,n} - J_{\alpha}^n) \right), \quad (4.12) $$

where $S$ is the set of interfaces between kinetic and fluid subdomains and $\beta = \pm 1$ depends on the orientation of the subdomains.

## 5. NUMERICAL SIMULATIONS

In the following, unless specified otherwise, the phase-space is discretized as follows:

$$ N_v = 16, \quad N_x = 50, \quad v_* = 10, \quad x_* = \pi, \quad \Delta t = 10^{-4}. $$

The same time step is used for all schemes. Note that since the limit scheme is explicit, its stability is therefore guaranteed under a parabolic condition: $\Delta t \le C \Delta x^2$. Let us assume that the electrical field is the gradient of a potential $V: E =$
$-\partial_x V$. To satisfy the periodicity of the domain, we choose $V(x) = -\frac{\sin(2x)}{4} \times 10^{-2}$ so $E(x) = \frac{1}{2} \cos(2x) \times 10^{-2}$. We also set

$$f_0^1 = \frac{1}{(2\pi)^{3/2}} e^{-|v|^2/2} (1 + \cos(2x)), \qquad (5.1)$$

an initial data at local equilibrium in velocity and

$$f_0^2 = \frac{1}{(2\pi)^{3/2}} |v|^4 e^{-|v|^2/2} (1 + \cos(2x)), \qquad (5.2)$$

an initial data far from the local equilibrium in velocity. Finally, we consider four configurations:

* Case 1: $E = 0$, with initial data (5.1);

* Case 2: $E \neq 0$, with initial data (5.1);

* Case 3: $E = 0$, with initial data (5.2);

* Case 4: $E \neq 0$, with initial data (5.2).

The properties of the fully kinetic implementation are shown in the 1D1D setting. The performance and properties of the hybrid method are presented in the full 1D3D one.

### 5.1. The full kinetic scheme.

**Convergence towards the drift-diffusion equation.** Let us first numerically investigate the AP property of the (Micro)-(Macro) scheme. We consider this analysis for the Cases 1 and 2. The results can be found in Figure 5 for Case 1 and in Figure 6 for Case 2. We can observe a convergence of the kinetic scheme to the limit one as $\varepsilon \to 0$. In particular, the curves for $\varepsilon = 0.05$ and $10^{-4}$ overlap and are close to the limit case. This validates the asymptotic consistency of the (Micro)-(Macro) scheme. The stability is numerically verified as the same $\Delta t$ is used for every $\varepsilon$.

**Long time behaviour.** The long time behaviour of solutions to $(P^\varepsilon)$ have been extensively studied in the past decades. In a more general setting, the electric field $E$ is the gradient of a potential $V \in C^2(\Omega_x)$, $E = -\nabla V$, and $(P^\varepsilon)$ admits a global equilibrium given by

$$F(x, v) = \frac{M_0}{\mu_0} e^{-(V(x) + \frac{|v|^2}{2})}, \quad (x, v) \in \Omega_x \times \mathbb{R}^{d_v} \qquad (5.3)$$

where $\mu_0 = (2\pi)^{-d/2} \int_{\Omega_x \times \mathbb{R}^{d_v}} e^{-(V(x) + \frac{|v|^2}{2})} dx \, dv$ and $M_0 = \int_{\Omega_x \times \mathbb{R}^{d_v}} f_0(x, v) dx \, dv$ is the mass of the initial condition. In particular, $F$ can be written under a separate variable form:

$$F(x, v) = M_0 \phi(x) \mathcal{M}(v), \quad \text{where} \quad \phi = \frac{e^{-V(x)}}{\int_{\Omega_x} e^{-V(x)} dx}. \qquad (5.4)$$

The functions $\mathcal{M}$ and $\phi$ are called local equilibria in velocity and position respectively. When one considers models such as $(P^\varepsilon)$, there are various ways to show that there exists $\kappa(\varepsilon) > 0$ and $C(\varepsilon) > 0$, such that if $f^\varepsilon$ is solution to $(P^\varepsilon)$,

$$\|f^\varepsilon(t) - F\|_V \le C(\varepsilon) \|f_0 - F\|_V e^{-\kappa(\varepsilon)t}, \qquad (5.5)$$

where $V$ is an appropriate functional space. A proof of (5.5) was done in [22] in a setting without electric field. In recent years, the literature on the subject expanded a lot. Robust and systematic methods were developed to show the convergence to an equilibrium. Those are called hypocoercivity methods. A general
<figure><img src="image_6.png" /><figcaption>FIGURE 5. Case 1. Comparison of the solution of the limit scheme ($S_{Lim}$) with the solution obtained with the (Micro)-(Macro) scheme with different $\epsilon$, $t = 0.0, 0.2, 1.0$ and $5.0$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIGURE 6. Case 2. Comparison of the solution of the limit scheme ($S_{Lim}$) with the solution obtained with the (Micro)-(Macro) scheme with different $\epsilon$, $t = 0.0, 0.2, 1.0$ and $5.0$.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>FIGURE 7. Case 3. Snapshots of the distribution function f computed with the scheme $(S_{Micro})-(S_{Macro})$, $\epsilon = 1.0$, $N_x = 100$, $N_v = 128$.</figcaption></figure>

abstract framework for the $H^1$ norm has been given in [41]. These methods often have the drawback of requiring regularity on the initial data. However the techniques have been adapted to consider only weighted $L^2$ initial data. More recently, an $L^2$-hypocoercivity method has been developed in [16] for linear kinetic equations. We also refer to [1] where $L^2$-hypocoercivity is shown for a more general kinetic equation. Both a self-consistent potential given by the Poisson equation and an exterior potential are considered and such a model is closer to the physics of semiconductors. From a numerical point of view, recovering such long-time behaviour at the discrete level is a significant property to obtain. In recent papers, hypocoercivity methods were adapted to the discrete setting using finite differences [17], finite elements [20] and finite volumes [6].

Following these ideas, we want to observe the convergence of the $(S_{Micro})-(S_{Macro})$ scheme to equilibrium in a large time scale. Figure 7 shows the evolution of the marginal distribution $f^\epsilon = \int_{\mathbb{R}^2} f^\epsilon dv_y dv_z$ as time increases (Case 3, $\epsilon = 1.0$). In particular, the numerical solution indeed seems to converge to equilibrium. Let us introduce the following discrete norm for $f = (f_{ij})_{ij}$:

$$ \|f\|_{\Delta}^{2} = \sum_{(i,j) \in \mathcal{I} \times \mathcal{J}^{3}} f_{ij}^{2} F_{ij}^{-1} \Delta x \Delta v^{3}, \qquad (5.6) $$

where $(F_{ij})_{ij}$ is a discretization of the global equilibrium $F$, $F_{ij} = F(x_i, v_j)$. For $(\rho_i)_{i \in \mathcal{I}}$, we also denote by $\|\rho\|_2 = \sum_{i \in \mathcal{I}} \rho_i^2 \Delta x$ the discrete $L^2$-norm in position. We now investigate the rate of convergence of the following discrete norms:

$$ \|f - F\|_{\Delta}, \quad \|g\|_{2}, \quad \|\rho^{\epsilon} - \langle F \rangle_{\Delta}\|_{2} \quad \text{and} \quad \|\rho - \langle F \rangle_{\Delta}\|_{2}, \qquad (5.7) $$

where $\rho^\epsilon$ is the density obtained with the kinetic scheme and $\rho$ is obtained with the limit scheme. We consider Case 1. On Figure 8 we choose $\epsilon = 0.5$ and $0.1$,
<figure><img src="image_9.png" /><figcaption>FIGURE 8. Case 1. Time evolution of the norms (5.7) computed with the fully kinetic scheme and limit scheme, $\varepsilon = 0.5$ (Left), $\varepsilon = 0.1$ (Right).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIGURE 9. Case 4. Time evolution of the norm $||f^{\varepsilon} - F||_{\Delta}$ computed with the fully kinetic scheme for $N_x = 20, 50, 100$ and $200$, $\varepsilon = 0.1$.</figcaption></figure>

and show the norms (5.7) as functions of time in semilog scale. The exponential convergence of the various norms is clear. Moreover, the rates $\kappa(\varepsilon)$ observed are $\kappa(0.5) = 2.43$ and $\kappa(0.1) = 3.72$. The rate $\kappa(\varepsilon)$ increases as the Knudsen number gets smaller. In particular we observe the same rate of convergence between the fully kinetic scheme and the limit one for small values of $\varepsilon$.

Let us point out that in the case of a nonzero electric field, we do not recover the same convergence to equilibrium. Indeed, our numerical scheme is not well-balanced, i.e. designed to preserve steady states. As a consequence, the numerical solution only converges to an equilibrium that is an approximation of the steady state. Figure 9 shows the convergence to the equilibrium as the number of cells in position increases fore Case 4.

## 5.2. Properties of the hybrid scheme.

**Choice of the coupling parameters.** Before investigating the properties of the hybrid scheme, a natural question is the choice of the coupling parameters $\eta_0$ and $\delta_0$. Indeed, as we have seen earlier, that choice has an impact on the conservation of mass. The smaller the parameters, the latter the coupling occurs and the more
<figure><img src="image_11.png" /><figcaption>FIGURE 10. Case 3. Snapshots of the difference between the densities computed with the kinetic and the limit schemes (Top), macroscopic indicator (Middle), $L^2(d\gamma)$ norm of the perturbation $g$ (Bottom), $\epsilon = 0.5$.</figcaption></figure>

one can control this variation. However, the bigger the parameters are, the faster is the resulting hybrid scheme as one allows more fluid cells to appear. Therefore, one must find a good balance between accuracy and computation time. To illustrate how the macroscopic indicator behaves, we compute it without updating the state of the cells. Figure 10 shows the indicator compared to the difference between the kinetic and fluid densities. One can observe that this indicator behaves as expected. When the kinetic and limit densities are close, the indicator is also small. Regarding the norm of $g^\epsilon$, its behaviour is also expected. Indeed, we chose an initial data far from the local equilibrium in velocity and therefore, the norm can be high even if the densities are close (See first column, third row in Figure 10). Lastly, both the macroscopic indicator and the norm of $g^\epsilon$ tend to 0 as time increases. As a consequence, the closer to the equilibrium the solution is, the more fluid cells will appear.

**Qualitative comparison.** Let us now compare the kinetic and the hybrid schemes. Figures 11 and 12 show the densities computed by the kinetic, hybrid and limit schemes for Cases 3 and 4 with $\epsilon = 0.1$. We can see a good agreement between the three schemes. The domain adaptation works for $E \neq 0$ which was not investigated in previous works on the method. As time increases, the solution relaxes to an equilibrium and the domain becomes fully fluid. One can also observe that when both type of cells co-exists the hybrid density slightly deviates from the full kinetic model.

**Conservation of mass.** Let us now numerically investigate the conservation of mass. Indeed, we have shown in Lemma 2 that the hybrid method is not exactly
<figure><img src="image_12.png" /><figcaption>FIGURE 11. Case 3. Snapshots of the densities computed using the full kinetic, hybrid and limit schemes, $\epsilon = 0.1$, $\eta_0 = \delta_0 = 10^{-4}$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIGURE 12. Case 4. Snapshots of the densities computed using the full kinetic, hybrid and limit schemes, $\epsilon = 0.1$, $\eta_0 = \delta_0 = 10^{-4}$.</figcaption></figure>

conservative. However, it becomes conservative asymptotically. In addition, the hybrid method was constructed so that the cells become fluid when the solution is close to a local equilibrium in velocity. As a consequence, the perturbation is small when the coupling occurs and so is the mass variation. In practice, we can observe a mass variation of the order of the machine accuracy. In addition, the zero-mass property of the perturbation

$$ \mu_g = \sum_{(i,j) \in \mathcal{I} \times \mathcal{J}^3} g_{i+\frac{1}{2},j} \Delta x \Delta v^3 $$
<figure><img src="image_14.png" /><figcaption>FIGURE 13. Case 3. Time evolution of the state of the cells (Left), mass variation and mass of $g^\epsilon$ (Right), $\epsilon = 0.5$, $\eta_0 = \delta_0 = 10^{-4}$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>FIGURE 14. Case 3. Time evolution of the state of the cells (Left), mass variation and mass of $g^\epsilon$ (Right), $\epsilon = 0.1$, $\eta_0 = \delta_0 = 10^{-4}$.</figcaption></figure>

is preserved. We illustrate the state of the cells, the corresponding mass variation as well as $\mu_g$ on Figures 13 and 14. Case 3 is considered for $\epsilon = 0.5$ and $0.1$.

**Error analysis.** We are now interested in the error induced by the hybrid method. In particular, we investigate the error between the full kinetic scheme and the hybrid method. The goal is to be more efficient than the full kinetic solver. However, the gain in computation time comes with a slight loss in accuracy. Let $f_{\text{Kinetic}}$ be the distribution computed using the full kinetic scheme and $f_{\text{Hybrid}}$ be the distribution obtained from the hybrid scheme. On Figure 15 we compute the error between the two distributions in $l^2$-norm: $||f_{\text{Kinetic}} - f_{\text{Hybrid}}||_2$ at several time steps with $\epsilon = 0.1$. The corresponding state of the cells can be found in Figure 14 for $\eta_0 = \delta_0 = 10^{-4}$. Quite expectedly, there indeed is a slight loss in accuracy as soon as the coupling occurs. However, it quickly diminishes as the coupled solution relaxes to equilibrium. Moreover, one can lessen this error by tuning down the coupling parameters but at the expense of the computational gain.

**Long time behaviour.** Similarly as for the full kinetic scheme, we are interested in the long time behaviour of the hybrid scheme. We shall focus on the case $E = 0$. Figure 16 shows the convergence of the norms (5.7) in the hybrid setting. As the
<figure><img src="image_16.png" /><figcaption>FIGURE 15. Case 3. Time evolution of the $l^2$-error between the kinetic distribution $f_{\text{Kinetic}}$ and the hybrid one $f_{\text{Hybrid}}$, $\varepsilon = 0.1$, $\eta_0 = \delta_0 = 10^{-4}$ (Left), $\eta_0 = \delta_0 = 10^{-3}$ (Right).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>FIGURE 16. Case 1. Time evolution of the norms (5.7) computed with the hybrid and limit schemes, $\varepsilon = 0.5$ (Left), $\varepsilon = 0.1$ (Right), $\eta_0 = \delta_0 = 10^{-4}$.</figcaption></figure>

perturbation is set to 0 in fluid subdomains, the norms of $g$ is directly projected to 0 when all cells are fluid. In addition, one can observe the convergence of the density towards the global mass. The rates of convergence are not exactly recovered compared to the full kinetic scheme but remains close: $\kappa_{\text{Hybrid}}(0.5) = 2.65$ and $\kappa_{\text{Hybrid}}(0.1) = 3.36$.

**Computation time.** Let us now consider the efficiency of the hybrid method. We set the final time $T = 5.0$ to compare the computation time. Tables 2-3-4-5 show the computation time of the full kinetic, hybrid and limit scheme for different test cases with two sets of coupling parameters: $\eta_0 = \delta_0 = 10^{-4}$ and $\eta_0 = \delta_0 = 10^{-3}$. We recall that the same time step, $\Delta t = 10^{-4}$, is used for the three schemes. We can make several observations. First, the fluid solver is as expected, much faster than the full kinetic one. Moreover it is also always faster than the hybrid method. This can easily be explained by the additional cost of computing the indicators and the added cost of dealing with interfaces between kinetic and fluid. Table 6 shows the computational gain for the previous tests. In particular, the hybrid method offers no significant gain in very low collision regimes. Because of the slow convergence rate towards the equilibrium, the coupling occurs very late or not at all
<table><caption>TABLE 2. Test 1. Comparison of the computation time (sec), ε = 1.0, T = 5.0, Nₓ = 50, Nᵧ = 16.</caption><thead><tr><th rowspan="2">Case</th><th colspan="4">η₀ = δ₀ = 10⁻⁴</th><th colspan="4">η₀ = δ₀ = 10⁻³</th></tr><tr><th>1</th><th>2</th><th>3</th><th>4</th><th>1</th><th>2</th><th>3</th><th>4</th></tr></thead><tbody><tr><th>MM</th><td>298.2</td><td>285.8</td><td>287.7</td><td>281.0</td><td>298.2</td><td>285.8</td><td>287.7</td><td>281.0</td></tr><tr><th>Hybrid</th><td>291.5</td><td>282.2</td><td>283.5</td><td>303.1</td><td>278.0</td><td>281.2</td><td>276.0</td><td>280.3</td></tr><tr><th>Limit</th><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td></tr></tbody></table>

<table><caption>TABLE 3. Test 2. Comparison of the computation time (sec), ε = 0.1, T = 5.0, Nₓ = 50, Nᵧ = 16.</caption><thead><tr><th rowspan="2">Case</th><th colspan="4">η₀ = δ₀ = 10⁻⁴</th><th colspan="4">η₀ = δ₀ = 10⁻³</th></tr><tr><th>1</th><th>2</th><th>3</th><th>4</th><th>1</th><th>2</th><th>3</th><th>4</th></tr></thead><tbody><tr><th>MM</th><td>290.1</td><td>295.7</td><td>282.0</td><td>295.8</td><td>290.1</td><td>295.7</td><td>282.0</td><td>295.8</td></tr><tr><th>Hybrid</th><td>181.4</td><td>210.7</td><td>152.3</td><td>210.1</td><td>116.8</td><td>131.0</td><td>110.4</td><td>131.0</td></tr><tr><th>Limit</th><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td></tr></tbody></table>

<table><caption>TABLE 4. Test 2. Comparison of the computation time (sec), ε = 0.05, T = 5.0, Nₓ = 50, Nᵧ = 16.</caption><thead><tr><th rowspan="2">Case</th><th colspan="4">η₀ = δ₀ = 10⁻⁴</th><th colspan="4">η₀ = δ₀ = 10⁻³</th></tr><tr><th>1</th><th>2</th><th>3</th><th>4</th><th>1</th><th>2</th><th>3</th><th>4</th></tr></thead><tbody><tr><th>MM</th><td>289.0</td><td>281.7</td><td>299.4</td><td>293.6</td><td>289.0</td><td>281.7</td><td>299.4</td><td>293.6</td></tr><tr><th>Hybrid</th><td>110.6</td><td>116.6</td><td>111.9</td><td>130.1</td><td>0.09</td><td>0.09</td><td>79.4</td><td>73.0</td></tr><tr><th>Limit</th><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td></tr></tbody></table>

<table><caption>TABLE 5. Test 4. Comparison of the computation time (sec), ε = 10⁻⁴, T = 5.0, Nₓ = 50, Nᵧ = 16.</caption><thead><tr><th rowspan="2">Case</th><th colspan="4">η₀ = δ₀ = 10⁻⁴</th><th colspan="4">η₀ = δ₀ = 10⁻³</th></tr><tr><th>1</th><th>2</th><th>3</th><th>4</th><th>1</th><th>2</th><th>3</th><th>4</th></tr></thead><tbody><tr><th>MM</th><td>278.6</td><td>296.5</td><td>286.7</td><td>290.1</td><td>278.6</td><td>296.5</td><td>286.7</td><td>290.1</td></tr><tr><th>Hybrid</th><td>0.11</td><td>0.10</td><td>0.10</td><td>0.10</td><td>0.09</td><td>0.10</td><td>0.11</td><td>0.10</td></tr><tr><th>Limit</th><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td></tr></tbody></table>

depending on the chosen final time. However, we can point out that the coupling algorithm doesn't add much cost relatively to the full kinetic scheme. This shows that the implementation is quasi-optimal. In addition, since the cost of the fluid model is negligible compared to the kinetic one, if fluid cells are used half the time, the computation time is essentially also reduced by half.

Another observation is that the speedup is testcase dependent as it can be seen in Table 6, ε = 0.05 and δ₀ = η₀ = 10⁻³. When the parameter is small the gain becomes extremely significant and the hybrid method becomes competitive with the fluid solver. A final observation is that the choice of larger coupling parameters indeed speeds up the method. If we consider Test 3, with ε = 0.1 and δ₀ = η₀ = 10⁻⁴ the speedup is 1.94 while a choice of δ₀ = η₀ = 10⁻³ offers a
<table><caption>TABLE 6. Speedup of the hybrid method compared to the full kinetic scheme, $N_x = 50$, $N_v = 16$.</caption><thead><tr><th rowspan="2">ε</th><th colspan="4">η₀ = δ₀ = 10⁻⁴</th><th colspan="4">η₀ = δ₀ = 10⁻³</th></tr><tr><th>1</th><th>2</th><th>3</th><th>4</th><th>1</th><th>2</th><th>3</th><th>4</th></tr></thead><tbody><tr><th>1.0</th><td>1.02</td><td>1.01</td><td>1.01</td><td>0.92</td><td>1.07</td><td>1.01</td><td>1.04</td><td>1.00</td></tr><tr><th>0.1</th><td>1.59</td><td>1.40</td><td>1.94</td><td>1.41</td><td>2.48</td><td>2.26</td><td>2.68</td><td>2.26</td></tr><tr><th>0.05</th><td>2.6</td><td>2.42</td><td>2.68</td><td>2.26</td><td>3211</td><td>3130</td><td>3.77</td><td>4.02</td></tr><tr><th>10⁻⁴</th><td>2533</td><td>2965</td><td>2867</td><td>2901</td><td>3096</td><td>2965</td><td>2606</td><td>2901</td></tr></tbody></table>

speedup of 2.68. Figure 15 compares the error made as time increases for these two sets of parameters. One can observe that the error reaches a maximum at two different orders of magnitude: $3.10^{-8}$ with $δ_0 = η_0 = 10^{-4}$ and $3.10^{-6}$ with $δ_0 = η_0 = 10^{-3}$. This last point raises the question of optimal parameters relatively to the error between hybrid and full kinetic. This shall be addressed in the future but seems, again, problem-dependant.

**Non homogeneous Knudsen number.** In this last numerical experiment, we place ourselves in the 1D-1D setting and consider a non homogeneous Knudsen number in the physical domain. Let us define the function

$$e(x) = \frac{1}{2}\left(\arctan\left(5 + 10\left(x - \frac{\pi}{2}\right)\right) + \arctan\left(5 - 10\left(x - \frac{\pi}{2}\right)\right)\right). \quad (5.8)$$

In particular, we choose $ε = ε(x)$ as

$$\varepsilon(x) = \frac{e(x)}{\max(e(x))}. \qquad (5.9)$$

Such a function admits a maximum of 1 in the center of the domain and decays to 0 near the boundaries. Physically, it corresponds to few collisions in the center of the domain and to a fluid behaviour elsewhere. In the following simulations, $Δt = 5.10^{-5}$ and the coupling parameters are $δ_0 = η_0 = 10^{-4}$. Note that depending on the choice of $ε(x)$ one may need to decrease the time step to ensure stability. From an implementation point of view, the constant $ε$ is simply replaced by $ε_{i+\frac{1}{2}} = ε(x_{i+\frac{1}{2}})$ without any change to the indicators.

Figure 17 shows that the hybrid scheme captures well the behaviour of the distribution. Indeed, we observe a fast relaxation where $ε(x)$ is small and a much slower one in the center of the domain where $ε(x)$ is around 1. Regarding the state of the cells, one can see on Figures 18 and 19 that the fluid solver is quickly used where $ε(x)$ is small. Moreover, the last cells to become fluid are the ones where the gradient of $ε(x)$ is large. It is explained by the nature of the macroscopic indicator which uses derivatives up to order 4. One can also observe that the hybrid density starts deviate from the full kinetic one as more fluid cells appear. However, this deviation occurs at a small scale: between $10^{-8}$ and $10^{10}$. In this setting, the variation of mass was of order $10^{-11}$. Finally, we looked at the convergence in time to the global equilibrium. On Figure 20, one can again observe the exponential convergence to equilibrium and the rate is slightly higher than the one obtained with an homogeneous value of $ε = 1$. Up to stability considerations, this experiment shows the robustness of the hybrid method. Performance-wise, considering Case 3
<figure><img src="image_18.png" /><figcaption>FIGURE 17. Case 3. Snapshots of distribution $f$ obtained with the hybrid scheme for a non homogeneous Knudsen number.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>FIGURE 18. Case 3. Snapshots of the densities computed using the full kinetic and hybrid schemes for a non homogeneous Knudsen number.</figcaption></figure>

with $N_x = 200$ and $N_v = 256$, the full kinetic scheme takes 183.7 seconds to run while the hybrid one takes 147.6 seconds offering a speedup of 1.24.

## 6. CONCLUSION

In this work, a new hybrid numerical method for linear kinetic equations in the diffusive scaling was presented. The method relies on two criteria motivated by a perturbative approach. The first one quantifies how far from a local equilibrium
<figure><img src="image_20.png" /><figcaption>FIGURE 19. Case 3. Time evolution of the state of the cells (Top) and mass variation (Bottom) for a non homogeneous Knudsen number.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>FIGURE 20. Case 3. Convergence to a global equilibrium, full kinetic(Left), hybrid (Right) for a non homogeneous Knudsen number.</figcaption></figure>

the distribution function is. The second criterion depends on the macroscopic quantities that are available on the whole computing domain. We have managed to quantify the mass variation induced by the method and we have shown that it is in practice very small. The method has proven to be efficient through various numerical experiments: the computational gain compared to a full kinetic scheme is significant. Moreover, the method performs well with a non-homogeneous Knudsen number which is encouraging to tackle more physically motivated problems.

In future works, a more general and physically relevant context will be considered. In particular, the multi-dimensional setting requires a smart implementation. Moreover, the load distribution over several cores of computation of the hybrid method is not straightforward. From a modelization point of view, a more general collision operator and a coupling with the Poisson equation shall be investigated. We are confident that the computational gain will be even more worthwhile in a full 3D-3D setting in the case of the Boltzmann operator which is known to be costly numerically.
ACKNOWLEDGEMENTS

T.L. was partially funded by Labex CEMPI (ANR-11-LABX-0007-01) and ANR
Project MoHyCon (ANR-17-CE40-0027-01). He also would like to thank Thomas
Rey and Marianne Bessemoulin-Chatard for the fruitful discussions and suggestions
that improved the paper.

REFERENCES

[1] ADDALA, L., DOLBEAULT, J., LI, X., AND TAYEB, M. L. L2-Hypocoercivity and large time asymptotics of the linearized Vlasov-Poisson-Fokker-Planck system. *Journal of Statistical Physics* (2019), 4.

[2] ALLAIRE, G., BLANC, X., DESPRES, B., AND GOLSE, F. Transport et diffusion. Lecture notes, downloaded from http://paestel.fr/content/transport-et-diffusion-g-allaire-x-blanc-b-despres-f-golse in August 2021.

[3] BARDOS, C., SANTOS, R., AND SENTIS, R. Diffusion Approximation and Computation of the Critical Size. *Transactions of the American Mathematical Society* (1984).

[4] BENNOUNE, M., LEMOU, M., AND MIEUSSENS, L. Uniformly stable numerical schemes for the Boltzmann equation preserving the compressible Navier-Stokes asymptotics. *Journal of Computational Physics* 227 (2008).

[5] BENSOUSSAN, A., LIONS, J.-L., AND PAPANICOLAU, G. C. Boundary Layers and Homogenization of Transport Processes. *Publications of The Research Institute for Mathematical Sciences* (1979).

[6] BESSEMOULIN-CHATARD, M., HERDA, M., AND REY, T. Hypocoercivity and diffusion limit of a finite volume scheme for linear kinetic equations. *Mathematics of Computation* 89, 323 (nov 2019), 1093–1133.

[7] BHATNAGAR, P. L., GROSS, E. P., AND KROOK, M. A Model for Collision Processes in Gases. I. Small Amplitude Processes in Charged and Neutral One-Component Systems. *Physical Review* 94 (May 1954), 511–525.

[8] CRESTETTO, A., CROUSEILLES, N., DIMARCO, G., AND LEMOU, M. Asymptotically complexity diminishing schemes (ACDS) for kinetic equations in the diffusive scaling. *Journal of Computational Physics* 394 (OCT 1 2019), 243–262.

[9] CROUSEILLES, N., AND LEMOU, M. An asymptotic preserving scheme based on a micro-macro decomposition for collisional Vlasov equations: diffusion and high-field scaling limits. *Kinetic and Related Models* 4, 2 (2011), 441–477.

[10] DEGOND, P., DIMARCO, G., AND MIEUSSENS, L. A moving interface method for dynamic kinetic-fluid coupling. *Journal of Computational Physics* 227 (2007), 1176–1208.

[11] DEGOND, P., GOUDON, T., AND POUPAUD, F. Diffusion limit for non homogeneous and non-micro-reversible processes. *Indiana University Mathematics Journal* (2000).

[12] DEGOND, P., JIN, S., AND MIEUSSENS, L. A smooth transition model between kinetic and hydrodynamic equations. *Journal of Computational Physics* 209 (2005), 665–694.

[13] DIMARCO, G., MIEUSSENS, L., AND RISPOLI, V. An asymptotic preserving automatic domain decomposition method for the Vlasov-Poisson-BGK system with applications to plasmas. *Journal of Computational Physics* (2014).

[14] DIMARCO, G., AND PARESCHI, L. Hybrid Multiscale Methods II. Kinetic Equations. *Multiscale Modeling and Simulation* 6 (2008), 1169–1197.

[15] DIMARCO, G., AND PARESCHI, L. Numerical methods for kinetic equations. *Acta Numerica* 23 (2014), 369–520.

[16] DOLBEAULT, J., MOUHOT, C., AND SCHMEISER, C. Hypocoercivity for linear kinetic equations conserving mass. *Transactions of the American Mathematical Society* (2015), 3807–3828. 21 pages.

[17] DUJARDIN, G., HÉRAU, F., AND LAFITTE-GODILLON, P. Coercivity, hypocoercivity, exponential time decay and simulations for discrete Fokker-Planck equations. *Numerische Mathematik* 144 (2018).

[18] EINKEMMER, L., HU, J., AND WANG, Y. An asymptotic-preserving dynamical low-rank method for the multi-scale multi-dimensional linear transport equation. *Journal of Computational Physics* (2021).
[19] FILBET, F., AND REY, T. A hierarchy of hybrid numerical methods for multi-scale kinetic equations. *SIAM Journal on Scientific Computing* 37, 3 (May 2015), A1218–A1247.

[20] GEORGOULIS, E. H. Hypocoercivity-compatible Finite Element Methods for the Long-time Computation of Kolmogorov's Equation. *SIAM Journal on Numerical Analysis* 59 (2018), 173–194.

[21] GOUDON, T., AND POUPAUD, F. Approximation by Homogenization and Diffusion of Kinetic Equations. *Communications in Partial Differential Equations* (2001).

[22] HÉRAU, F., AND NIER, F. Isotropic hypoellipticity and trend to the equilibrium for the Fokker-Planck equation with high degree potential. *Archive for Rational Mechanics and Analysis* 171, 2 (2004), 151–218.

[23] HOCHBRUCK, M., AND OSTERMANN, A. Exponential integrators. *Acta Numerica* 19 (may 2010), 209–286.

[24] HORSTEN, N., SAMAEY, G., AND BAELMANS, M. Hybrid fluid-kinetic model for neutral particles in the plasma edge. *Nuclear Materials and Energy* (2019).

[25] HORSTEN, N., SAMAEY, G., AND BAELMANS, M. A hybrid fluid-kinetic model for hydrogenic atoms in the plasma edge of tokamaks based on a micro-macro decomposition of the kinetic equation. *J. Comput. Phys.* 409 (2020), 109308.

[26] JIN, S. Efficient Asymptotic-Preserving (AP) Schemes For Some Multiscale Kinetic Equations. *SIAM Journal on Scientific Computing* 21 (1999), 441–454.

[27] JIN, S. Asymptotic-Preserving Schemes for Multiscale Physical Problems. *Acta Numerica* (2022). In Press.

[28] KLAR, A. Asymptotic-Induced Domain Decomposition Methods for Kinetic and Drift Diffusion Semiconductor Equations. *SIAM Journal on Scientific Computing* 19 (1998), 2032–2050.

[29] KLAR, A. A Numerical Method for Kinetic Semiconductor Equations in the Drift-Diffusion Limit. *SIAM Journal on Scientific Computing* 20 (1999), 1696–1712.

[30] KLAR, A., AND SCHMEISER, C. Numerical Passage from Radiative Heat Transfer to Nonlinear Diffusion Models. *Mathematical Models and Methods in Applied Sciences* 11 (2001), 749–767.

[31] KOLOBOV, V. I., ARSLANBEKOV, R. R., ARISTOV, V. V., FROLOVA, A. A., AND ZABELOK, S. A. Unified solver for rarefied and continuum flows with adaptive mesh and algorithm refinement. *Journal of Computational Physics* 223 (2007), 589–608.

[32] LEMOU, M. Relaxed micro-macro schemes for kinetic equations. *Comptes Rendus Mathématiques* 348, 7-8 (apr 2010), 455–460.

[33] LEMOU, M., AND MIEUSSENS, L. A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit. *SIAM Journal on Scientific Computing* 31, 1 (2008), 334–368.

[34] LEVERMORE, C. D., MOROKOFF, W. J., AND NADIGA, B. T. Moment realizability and the validity of the Navier-Stokes equations for rarefied gas dynamics. *Physics of Fluids* 10 (1998), 3214–3226.

[35] LIU, J.-G., AND MIEUSSENS, L. Analysis of an asymptotic preserving scheme for linear kinetic equations in the diffusion limit. *SIAM Journal on Numerical Analysis* 48, 4 (2010), 1474–1491.

[36] POUPAUD, F. Diffusion approximation of the linear semiconductor Boltzmann equation: analysis of boundary layers. *Asymptotic Analysis* (1991).

[37] TIWARI, S. Coupling of the Boltzmann and Euler Equations with Automatic Domain Decomposition. *Journal of Computational Physics* 144 (1998), 710–726.

[38] TIWARI, S. Application of Moment Realizability Criteria for the Coupling of the Boltzmann and Euler Equations. *Transport Theory and Statistical Physics* 29 (2000), 759–783.

[39] TIWARI, S., AND KLAR, A. An adaptive domain decomposition procedure for Boltzmann and Euler equations. *Journal of Computational and Applied Mathematics* 90 (1998), 223–237.

[40] TIWARI, S., KLAR, A., AND HARDT, S. A particle-particle hybrid method for kinetic and continuum equations. *Journal of Computational Physics* 228 (2009), 7109–7124.

[41] VILLANI, C. Hypocoercivity. *Memoirs of the American Mathematical Society* (Sept. 2009).

[42] WEINAN, E. *Principles of Multiscale Modeling*. Cambridge University Press, Cambridge, 2011.

[43] WEINAN, E., AND ENGQUIST, B. The Heterogeneous Multiscale Methods. *Communications in Mathematical Sciences* 1, 1 (2003), 87–132.