# Hierarchical dynamic domain decomposition for the multiscale Boltzmann equation

Domenico Caparello, Lorenzo Pareschi, Thomas Rey

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

## To cite this version:

Domenico Caparello, Lorenzo Pareschi, Thomas Rey. Hierarchical dynamic domain decomposition for the multiscale Boltzmann equation. 2025. hal-05057640v2

HAL Id: hal-05057640

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

Preprint submitted on 25 Aug 2025

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# HIERARCHICAL DYNAMIC DOMAIN DECOMPOSITION FOR THE MULTISCALE BOLTZMANN EQUATION

DOMENICO CAPARELLO, LORENZO PARESCHI, AND THOMAS REY

**ABSTRACT.** In this work, we present a hierarchical domain decomposition method for the multi-scale Boltzmann equation based on moment realizability matrices, originally introduced by Levermore, Morokoff, and Nadiga [42]. This criterion dynamically partitions the two-dimensional spatial domain into three regimes: the Euler regime, an intermediate kinetic regime governed by the ES-BGK model, and the full Boltzmann regime. The key advantage of this approach is the selective use of solvers: the Euler equations efficiently handle near-equilibrium regions, the ES-BGK model captures moderately non-equilibrium states, and the full Boltzmann solver is applied only where strong kinetic effects dominate, such as near shocks and boundary layers. This results in both high accuracy and substantial computational savings. To guarantee robust coupling across regimes, we employ asymptotic-preserving (AP) schemes and fast spectral solvers for the Boltzmann collision operator. The main novelties of the present work are: (i) a full 2D-in-space and 3D-in-velocity hierarchical dynamic domain decomposition, (ii) the integration of three distinct physical regimes within a single solver framework, and (iii) a parallelized implementation exploiting CPU multithreading. Taken together, these features make our approach one of the first to combine the accuracy of deterministic full Boltzmann solvers with the efficiency of adaptive domain decomposition, enabling scalable simulations of multiscale kinetic flows in complex geometries.

## CONTENTS

<table><tr><td>1. Introduction</td><td>2</td></tr><tr><td>2. The Boltzmann equation</td><td>3</td></tr><tr><td>2.1. The Boltzmann operator</td><td>4</td></tr><tr><td>2.2. The ES-BGK operator</td><td>5</td></tr><tr><td>3. Regime indicators</td><td>6</td></tr><tr><td>3.1. The Chapmann-Enskog expansion</td><td>6</td></tr><tr><td>3.2. The moment realizability criterion</td><td>9</td></tr><tr><td>3.3. Three layers hierarchical model</td><td>10</td></tr><tr><td>4. Numerical schemes</td><td>12</td></tr><tr><td>4.1. Spatial, velocity and time settings</td><td>12</td></tr><tr><td>4.2. Euler equations - Layer 0</td><td>13</td></tr><tr><td>4.3. ES-BGK model - Layer 1</td><td>13</td></tr><tr><td>4.4. Boltzmann equation - Layer 2</td><td>15</td></tr><tr><td>4.5. Evolving the solution and coupling the models</td><td>20</td></tr><tr><td>5. Numerical simulations</td><td>23</td></tr><tr><td>5.1. Test 1: Sod shock tube</td><td>24</td></tr><tr><td>5.2. Test 2: Fluid flux towards a fixed rectangular obstacle</td><td>25</td></tr><tr><td>5.3. Test 3: Moving object in a fluid at rest</td><td>26</td></tr><tr><td>5.4. Test 4: Fluid flux towards a fixed round obstacle</td><td>29</td></tr><tr><td>5.5. Test 5: Non uniform Knudsen number</td><td>29</td></tr></table>

<footnote n="2010">2010 Mathematics Subject Classification. Primary: 76P05, 82C40, Secondary: 65N08, 65N35.</footnote>

<footnote>Key words and phrases. multiscale flows, Boltzmann equation, compressible Euler system, ES-BGK operator, Burnett transport coefficients, asymptotic-preserving methods, domain decomposition method, spectral methods.</footnote>
5.6. Discussion on speed-up 34

6. Conclusions 35

Appendix A. Moment realizability matrix 35

Acknowledgments 37

References 38

# 1. INTRODUCTION

Many problems in engineering and physics involve fluid flows in transitional regimes, such as atmospheric reentry of spacecraft, micro-electro-mechanical systems (MEMS), and industrial heat transfer systems. In these settings, classical macroscopic models like the Euler equations are often inadequate due to the presence of shocks, kinetic boundary layers, or other non-equilibrium effects. Accurately capturing such dynamics requires kinetic models that resolve the microscopic behavior of the gas.

Kinetic theory is a cornerstone of modern applied mathematics and physics. Its roots trace back to foundational thermodynamic works by Bernoulli, Avogadro, and Carnot in the 18th and 19th centuries. One of the earliest mathematical contributions was by James Clerk Maxwell, who, in 1859, introduced a statistical model for particle interactions in his study of Saturn's rings [45]. This was followed by Ludwig Boltzmann's celebrated formulation of the Boltzmann equation [5], which remains central to the modeling of rarefied gas dynamics.

Kinetic equations describe the evolution of a particle distribution function in a high-dimensional phase space, typically six-dimensional for a three-dimensional configuration domain. These models combine linear transport with stiff, nonlinear collision operators, making them significantly more computationally expensive than their fluid-dynamical counterparts. This “curse of dimensionality” motivates the use of macroscopic equations such as Euler or Navier-Stokes models in regions where the flow is close to equilibrium, while reserving kinetic solvers for regions where non-equilibrium effects are dominant.

This leads to the concept of hybrid multiscale models and hybrid numerical schemes, which have been widely developed in the literature. These methods blend kinetic and fluid models, often using domain decomposition and switching criteria based on physical indicators. Relevant examples include domain decomposition methods [13, 15, 54, 56], asymptotic-preserving (AP) schemes [17, 19, 20, 33, 34], and Boltzmann solvers with particle/DSMC coupling [7, 9, 12, 14–16, 18, 57]. Our approach is inspired by these developments and by the foundational works of Levermore, Morokoff, and Nadiga [42] and Filbet and Rey [22], but differs by introducing a three-regime decomposition strategy implemented in a fully adaptive, non-intrusive fashion.

The proposed hierarchical method offers two main advantages:

* **High accuracy in multiscale regimes compared to classical fluid solvers:** The method dynamically applies the appropriate physical model in each region—Euler equations in equilibrium zones, ES-BGK in transitional regimes, and the full Boltzmann equation in highly non-equilibrium areas—ensuring accurate resolution of shocks, boundary layers, and other kinetic effects where fluid models fail.

* **Increased computational efficiency compared to full Boltzmann solvers:** By confining the expensive kinetic computations to a minimal subset of the domain and relying on efficient fluid solvers elsewhere, the method achieves significant speedups over fully kinetic simulations without compromising physical fidelity.

The ES-BGK model plays a crucial role in this hierarchy: it allows the method to handle moderately non-equilibrium regions where the fluid description is insufficient, yet applying the full Boltzmann model would be unnecessarily costly. This three-layer decomposition offers a computational advantage over classical two-regime methods (e.g., Euler-Boltzmann [22]) while preserving greater physical accuracy.
To maximize computational efficiency, we minimize the size of the full Boltzmann region during simulations. The algorithm is designed to be non-intrusive with respect to the solvers, meaning that it can work with different implementations of both kinetic and fluid models.

We emphasize that this hierarchical structure makes feasible simulations with the full Boltzmann operator within reasonable computational times, since it is applied only in regions where it is strictly necessary. This selective use ensures both significant speedups and high accuracy, thus providing a practical and efficient tool for simulating non-equilibrium multiscale flows.

It is well known that solving the Boltzmann equation in stiff regimes is especially challenging. In these regimes, frequent particle collisions rapidly drive the system toward local thermodynamic equilibrium, allowing for a reduced description via macroscopic or diffusive equations [18]. However, in transitional regimes, collisions are sufficiently strong to introduce stiffness but not strong enough to restore equilibrium, leading to severe time step restrictions for explicit schemes.

Asymptotic-preserving (AP) schemes offer a solution: these numerical methods remain stable and accurate across a range of Knudsen numbers, including the stiff limit, without the need for mesh or time step refinement [8, 11, 27, 32–36]. Our method incorporates AP strategies to ensure consistent behaviour across regimes and uses fast spectral solvers for the Boltzmann operator [21, 46] to improve performance.

To our knowledge, these spectral AP techniques represent the most efficient deterministic approach currently available for solving the full Boltzmann collision operator with high accuracy. The novelty of the present work is that these state-of-the-art techniques are further enhanced by a dynamic domain decomposition strategy. This hierarchical coupling allows computational resources to be concentrated only where a kinetic description is required, thereby achieving substantial savings. Other strategies for multiscale fluid problems include recent developments with neural-network-based solvers [25, 40, 60], or Discrete Unified Gas Kinetic Schemes (DUGKS) [29–31, 59]. While these methods offer valuable alternatives, we stress that our approach explicitly incorporates the full Boltzmann hard-sphere operator, thereby achieving higher physical accuracy than simplified collision models such as BGK-type operators, while still optimizing computational resources through domain decomposition.

In this work, we develop and implement this hierarchical strategy in a fully two-dimensional spatial domain with three-dimensional velocity dependence. We describe in detail the moment realizability-based switching criterion, the structure of the hybrid solver, and its coupling across regimes. Through benchmark problems with strong shocks and boundary layers, we demonstrate that the approach achieves the accuracy of full Boltzmann solvers while delivering substantial computational savings relative to existing spectral-AP schemes. This makes the proposed framework one of the most efficient deterministic strategy currently available for multiscale kinetic flows in complex geometries.

The rest of the paper is organized as follows: in Section 2, we introduce the kinetic and fluid models; in Section 3, we describe the hybrid decomposition indicators whereas in Section 4 we describe the details of the numerical schemes adopted; in Section 5, we present numerical experiments; and in Section 6, we conclude with a discussion of future directions.

## 2. THE BOLTZMANN EQUATION

Let us now present the general kinetic model we are interested in. For a given nonnegative initial condition $f_0$, we will study a particle distribution function $f^\epsilon = f^\epsilon(t, x, v)$, for $t \ge 0$, $x \in \Omega \subset \mathbb{R}^{d_x}$ and $v \in \mathbb{R}^3$, solution to the initial-boundary value problem

$$ (2.1) \qquad \begin{cases} \frac{\partial f^\epsilon}{\partial t} + v \cdot \nabla_x f^\epsilon = \frac{1}{\epsilon} Q(f^\epsilon), \\ f^\epsilon(0, x, v) = f_0(x, v), \end{cases} $$

where the collision operator $Q$ is a Boltzmann-like operator, and $\epsilon > 0$ is the Knudsen number (it is equal to the ratio between the mean free path of particles and the physical length scale of the system).
This quantity gives an indication of the rarefaction of the gas: if $\epsilon \ll 1$ it means that the particle system is subject to many collisions and then the gas is dense and in fluid regime; vice-versa, if $\epsilon \sim 1$ it means that the number of collisions is small and then the gas can be considered rarefied [22].

The open set $\Omega$ is a bounded Lipschitz-continuous domain of $\mathbb{R}^{d_x}$, which means that the model (2.1) has to be supplemented with boundary conditions described later.

We assume that the collision operator fulfils the three following assumptions

(H_1) It preserves mass, momentum and kinetic energy. This properties can be, respectively, expressed as

$$ \int_{\mathbb{R}^3} Q(f)(v) dv = 0, \quad \int_{\mathbb{R}^3} Q(f)(v) v dv = 0, \quad \int_{\mathbb{R}^3} Q(f)(v) |v|^2 dv = 0; $$

(H_2) It dissipates the Boltzmann entropy (H-theorem)

$$ \int_{\mathbb{R}^3} Q(f)(v) \log(f)(v) dv \le 0; $$

(H_3) Its equilibria are given by Maxwellian distributions

$$ Q(f) = 0 \Leftrightarrow f = \mathcal{M}_{\rho, \mathbf{u}, T} := \frac{\rho}{(2\pi T)^{3/2}} \exp\left(-\frac{|v - \mathbf{u}|^2}{2T}\right), $$

where the *density*, *velocity* and *temperature* of the gas $\rho$, $\mathbf{u}$ and $T$ are computed from the distribution function $f$ as

$$ \rho = \int_{v \in \mathbb{R}^3} f(v) dv, \quad \mathbf{u} = \frac{1}{\rho} \int_{v \in \mathbb{R}^3} v f(v) dv, \quad T = \frac{1}{3\rho} \int_{v \in \mathbb{R}^3} |\mathbf{u} - v|^2 f(v) dv. $$

According to assumptions (H_2)-(H_3), it can be proved that, when $\epsilon \to 0$, the distribution $f^\epsilon$ converges (at least formally) to a Maxwellian distribution, whose moments are solution to the compressible Euler system

$$ (2.2) \qquad \left\{ \begin{array}{l} \partial_t \rho + \operatorname{div}_x (\rho \mathbf{u}) = 0, \\ \\ \partial_t (\rho \mathbf{u}) + \operatorname{div}_x (\rho \mathbf{u} \otimes \mathbf{u} + \rho T \mathbf{I}) = \mathbf{0}_{\mathbb{R}^3}, \\ \\ \partial_t E + \operatorname{div}_x (\mathbf{u} (E + \rho T)) = 0. \end{array} \right. $$

This limit is typically referred to as an *hydrodynamic limit*.

In the following we are going to present two operators that satisfy these hypothesis, and that will be implemented in our hybrid domain decomposition scheme.

**2.1. The Boltzmann operator.** If we consider only binary elastic collisions between particles, then a dilute gas can be described by the kinetic equation with the Boltzmann operator, and in a velocity space of dimension $d$ the collision operator is

$$ (2.3) \qquad Q(f, f)(v) = \int_{\mathbb{R}^d} \int_{S^{d-1}} B(|v - v_*|, \cos \theta) [f'_* f' - f_* f] d\sigma dv_*, $$

where we used the notation $f = f(v)$, $f_* = f(v_*)$, $f' = f(v')$, $f'_* = f(v'_*)$. The pair $(v, v_*)$ are the pre-collisional velocities, while the quantities $(v', v'_*)$ are the post-collisional velocities. These two pairs are related by the following relations

$$ v' = \frac{v + v_*}{2} + \frac{|v - v_*|}{2} \sigma, \quad v'_* = \frac{v + v_*^*}{2} - \frac{|v - v_*|}{2} \sigma. $$
The collision kernel $B$ is a non-negative function which by physical arguments of invariance only depends on $|v - v_*|$ and $\cos \theta = \hat{u} \cdot \sigma$, where $\hat{u} = (v - v_*)/|v - v_*|$.

It can be proved that the Boltzmann collision operator has the fundamental properties $(\mathbf{H}_1)$, $(\mathbf{H}_2)$ and $(\mathbf{H}_3)$.

In the remainder of this paper we shall use the assumption of short range interactions, and in particular assume that $B$ is locally integrable. This assumption is satisfied by the hard-spheres model, which in dimension $d = 3$ can be written as

$$ (2.4) \qquad B(|v - v_*|, \cos(\theta)) = |v - v_*|. $$

As an important benchmark model for the numerical simulation, and for the development of fast spectral methods for the computation of the Boltzmann operator, we therefore introduce the so-called variable hard-spheres model (VHS), which writes

$$ (2.5) \qquad B(|v - v_*|, \cos(\theta)) = C_{\gamma}|v - v_*|^{\gamma}, $$

for some $\gamma \in [0, 1]$ and a constant $C_{\gamma} > 0$.

It is important to remark that, for this last class of models, it is possible to split the collision operator as follows

$$ (2.6) \qquad Q(f, f) = Q^{+}(f, f) - L(f)f, $$

where

$$ (2.7) \qquad Q^{+}(f, f) = \int_{R^d} \int_{S^{d-1}} B(|v - v_*|, \cos(\theta)) f' f_*' d\sigma dv_*, $$

and

$$ (2.8) \qquad L(f) = \int_{R^d} \int_{S^{d-1}} B(|v - v_*|, \cos(\theta)) f_* d\sigma dv_* $$

The computation of Boltzmann-like collision operator, performed using deterministic numerical methods, requires to work on a bounded velocity space. In most cases, the classical technique is to periodize the particle distribution function and the collision operator (which correspond to add some non physical binary collisions), and then using Fourier series to compute the truncate operator. Unfortunately this causes the loss of some local invariants, but the mass can be preserved by implementing an accurate periodization, which is the basis of spectral methods [50], that we discuss in further details in Section 4.4.1.

## 2.2. The ES-BGK operator.
The BGK operator [4], and its Ellipsoidal Statistical (ES-BGK) extension [3], are two other well known collision operators which have the properties $(\mathbf{H}_1)-(\mathbf{H}_2)-(\mathbf{H}_3)$. The bilinear Boltzmann collision operator $Q$ is replaced by a nonlinear relaxation operator, which matches the same hydrodynamic limit than the Boltzmann operator.

First of all we introduce some macroscopic quantities of the particle distribution function $f$ such as the stress tensor

$$ \rho_f \Theta_f(t, x) = \int_{\mathbb{R}^3} (v - \mathbf{u}_f) \otimes (v - \mathbf{u}_f) f(t, x, v) dv. $$

Therefore the translational temperature is related to the stress tensor as $T_f = \text{Trace}(\Theta_f)/3$. We finally introduce the corrected tensor

$$ \mathcal{T}_f(t, x) = [(1 - \beta) T_f \mathbf{I} + \beta \Theta_f] (t, x), $$

which can be viewed as a linear combination of the initial stress tensor $\Theta_f$ and of the isotropic stress tensor $T_f \mathbf{I}$ developed by a Maxwellian distribution. The parameter $-\infty < \beta < 1$ is used to modify the value of the Prandtl number through the formula

$$ 0 \le \text{Pr} = \frac{1}{1 - \beta} \le +\infty \quad \text{for } \beta \in (-\infty, 1). $$
The correct Prandtl number for a monatomic gas of hard-spheres is equal to 2/3 [53], namely obtained here for $\beta = -1/2$ whereas the classical BGK operator, obtained for $\beta = 0$, has a Prandtl number equal to 1.

To define the ES-BGK operator, we introduce a corrected Gaussian $G[f]$ defined by

$$G[f] = \frac{\rho_f}{\sqrt{\det(2\pi T_f)}} \exp\left(-\frac{(v - \mathbf{u}_f) T_f^{-1} (v - \mathbf{u}_f)}{2}\right),$$

and the corresponding collision operator is now

$$(2.9) \qquad Q_{ES-BGK}(f) = \nu(\rho_f, T_f)(G[f] - f),$$

where $\nu$ is the collision frequency from the Boltzmann operator. It can be shown [53] that it depends only on the kinetic density $\rho_f$ and temperature $T_f$.

### 3. REGIME INDICATORS

Most of the papers in literature about decomposition methods rely on the domain decomposition introduced by Boyd, Chen and Chandler in [6], in which a macroscopic criterion is used to switch between the hydrodynamic description (easy to compute numerically, but inaccurate near shocks or boundary layers) and the kinetic one (computationally expansive but accurate in most of the situations).

This criterion is based on the value of the local Knudsen number: if it is below a (problem-dependent) threshold, then the kinetic description is used. In fact, it is important to remark that the Knudsen number is defined as the ratio between the mean free path of the particles and the characteristic physical length scale of the system: so when the Knudsen number is very small, then it means that the gas is rarefied and so the usage of kinetic models is necessary.

Several authors have implemented methods of this kind, considering also the spatial derivatives in the numerical solution itself, among which we recall Kolobov, et al. in [37], and Degond and Dimarco in [12]. In literature it is possible to find other works where different, but similar, indicators are used. For example Tiwari in [54] introduced a criterion based on the viscous and heat fluxes of the Navier-Stokes equation, through a Grad's 13-moments expansion. Similar strategies have been used, among the others, by Tiwari, Klar and Hardit in [57], and by Alaia and Puppo in [1, 2].

In this work we will use two criteria: one to detect when the hydrodynamic description is not valid any more (and then it is necessary to switch to the kinetic one), and the other one to understand when to switch from the kinetic regime to the hydrodynamic one.

The idea of using hybrid decomposition solver is not limited only to the Boltzmann equation, but other recent works exploiting the lower dimensionality of the asymptotic fluid model for different kinetic equations, like the Vlasov-BGK or the Vlasov-Poisson-BGK equations in the diffusive scaling are [13, 38, 39]. This shows the growing interest and need for the development of hybrid domain decomposition schemes to speed up the calculation of the solution for kinetic equations.

In the following, we present the criteria, and we start introducing a key mathematical tool: the Chapmann-Enskog expansion.

3.1. **The Chapmann-Enskog expansion.** The criterion used in this article was introduced by Levermore, Morokoff and Nadiga in [42], and it has already been used by other authors, like Tiwari in [55], Filbet and Rey in [22], Filbet and Xiong in [23], Li, Song and Wang in [44] and Xiong and Qiu in [61]. In particular, in this last paper, Xiong and Qiu propose a hierarchical numerical scheme based on the following three levels (Euler equations, Navier-Stokes equation and kinetic equation). This approach has the main interest to depend on the closure made for obtaining the hydrodynamic model. For the hydrodynamic regions it is possible to consider either the compressible Euler equations, or the compressible Navier-Stokes, or other models, such as Burnett or Super-Burnett [53]. In our code we implemented the Euler equations.
Assuming that the distribution $f^\epsilon$ is close to equilibrium thanks to the relaxation property ($\mathbf{H}_2$), we
can do formally the *Chapman-Enskog* expansion

$$
(3.1) \qquad f^{\epsilon} = \mathcal{M}_{\rho^{\epsilon}, \mathbf{u}^{\epsilon}, T^{\epsilon}} \left[ 1 + \epsilon g^{(1)} + \epsilon^2 g^{(2)} + \dots \right],
$$

where the fluctuations $g^{(i)}$ for $i \ge 0$ designate a function that depends smoothly on

$$
(\rho^{\epsilon}, \rho^{\epsilon} \mathbf{u}^{\epsilon}, E^{\epsilon}) = \int_{\mathbb{R}^3} f^{\epsilon}(v) \left(1, v, \frac{|v|^2}{2}\right) dv, \quad T^{\epsilon} = \frac{1}{3} \rho^{\epsilon} \int_{\mathbb{R}^3} f^{\epsilon}(v) |v - \mathbf{u}^{\epsilon}|^2 dv,
$$

and any finite number of its derivatives with respect to the $x$-variable at the same point $(t, x)$, and on
the $v$-variable. According to $(\mathbf{H}_1)$, it verifies

$$
\int_{\mathbb{R}^3} g^{(i)}(v) \left(1, v, \frac{|v|^2}{2}\right) dv = \mathbf{0}_{\mathbb{R}^5}^T.
$$

Without any closure, according to the conservative properties ($\mathbf{H}_1$) of the collision operator $\mathcal{Q}$, we have

$$
(3.2) \qquad \left\{
\begin{aligned}
& \partial_t \rho^\epsilon + \operatorname{div}_x (\rho^\epsilon \mathbf{u}^\epsilon) = 0, \\
& \partial_t (\rho^\epsilon \mathbf{u}^\epsilon) + \operatorname{div}_x (\rho^\epsilon \mathbf{u}^\epsilon \otimes \mathbf{u}^\epsilon + \rho^\epsilon T^\epsilon (\mathbf{I} + \bar{\mathbf{A}}^\epsilon)) = \mathbf{0}_{\mathbb{R}^3}, \\
& \partial_t E^\epsilon + \operatorname{div}_x \left( \frac{1}{2} \rho^\epsilon |\mathbf{u}^\epsilon|^2 \mathbf{u}^\epsilon + \rho^\epsilon T^\epsilon \left( \frac{5}{2} \mathbf{I} + \bar{\mathbf{A}}^\epsilon \right) \mathbf{u}^\epsilon + \rho^\epsilon (T^\epsilon)^{3/2} \bar{\mathbf{B}}^\epsilon \right) = 0,
\end{aligned}
\right.
$$

where the traceless matrix $\bar{\mathbf{A}}^\epsilon \in M_3$ and the vector $\bar{\mathbf{B}}^\epsilon \in \mathbb{R}^3$ are given by [22]

$$
(3.3) \qquad \left\{
\begin{array}{ll}
\bar{\mathbf{A}}^{\epsilon} & := \displaystyle \frac{1}{\rho^{\epsilon}} \int_{\mathbb{R}^3} \mathbf{A}(\mathbf{V}) f^{\epsilon}(v) \, dv, \quad \mathbf{A}(\mathbf{V}) = \mathbf{V} \otimes \mathbf{V} - \frac{|\mathbf{V}|^2}{3} \mathbf{I}, \\
\\
\bar{\mathbf{B}}^{\epsilon} & := \displaystyle \frac{1}{\rho^{\epsilon}} \int_{\mathbb{R}^3} \mathbf{B}(\mathbf{V}) f^{\epsilon}(v) \, dv, \quad \mathbf{B}(\mathbf{V}) = \frac{1}{2} [|\mathbf{V}|^2 - 5] \mathbf{V},
\end{array}
\right.
$$

and where we used the shorthand

$$
\mathbf{V}(v) = \frac{v - \mathbf{u}}{\sqrt{T}}.
$$

Therefore, depending on the order in $\epsilon$ of the truncation of the series (3.1), we will obtain different
hydrodynamic description of the fluid.

3.1.1. *Zeroth order: compressible Euler system.* At zeroth order with respect to $\epsilon$, we have $f^\epsilon = \mathcal{M}_{\rho,\mathbf{u},T}$.
This distribution is isotropic in $v-\mathbf{u}$ and its odd moments with respect to $(v-\mathbf{u})$ are all equal to zero.
Since the matrix $\bar{\mathbf{A}}^\epsilon$ is traceless, we then have that

$$
\bar{\mathbf{A}}_{Euler} = \mathbf{0}_{M_3}.
$$

Moreover, since $\bar{\mathbf{B}}^\epsilon$ involves odd, centered moments of $f^\epsilon$, we also obtain

$$
\bar{\mathbf{B}}_{Euler} = \mathbf{0}_{\mathbb{R}^3}.
$$

Hence, the moments $(\rho, \mathbf{u}, T)$ are solution to the compressible Euler system (2.2). We notice in particular
that the Maxwellian distribution in (3.1) is independent of $\epsilon$.
3.1.2. First order: compressible Navier-Stokes system. Going to the next order in $\epsilon$, we plug-in the expansion (3.1) in the Boltzmann equation (2.1). Since the Maxwellian distribution is an equilibrium of the collision operator (according to $(\mathbf{H}_3)$), the fluctuation $g^{(1)}$ is given by

$$ (3.4) \qquad \partial_t \mathcal{M}_{\rho, \mathbf{u}, T} + v \cdot \nabla_x \mathcal{M}_{\rho, \mathbf{u}, T} = \mathcal{L}_{\mathcal{M}_{\rho, \mathbf{u}, T}} g^{(1)} + \mathcal{O}(\epsilon), $$

where $(\rho, \mathbf{u}, T)$ are the unknowns to be determined and $\mathcal{L}_{\mathcal{M}}$ is the linearized<fnref n="1" /> collision operator around the Maxwellian distribution.

Besides, using a logarithmic derivative, one has<fnref n="2" />

$$ \partial_t \mathcal{M}_{\rho, \mathbf{u}, T} + v \cdot \nabla_x \mathcal{M}_{\rho, \mathbf{u}, T} = \mathcal{M}_{\rho, \mathbf{u}, T} \left[ \frac{1}{\rho} \partial_t \rho + \frac{1}{\rho} v \cdot \nabla_x \rho + \frac{1}{\sqrt{T}} (\mathbf{V} \cdot \partial_t \mathbf{u} + \mathbf{V} \otimes v : \nabla_x \mathbf{u}) + \frac{1}{2T} (|\mathbf{V}|^2 - 3) (\partial_t T + v \cdot \nabla_x T) \right]. $$

Then, using the conservation laws (3.2), we replace the time derivatives by spatial ones, and drop the terms of order $\epsilon$ in (3.4), hence we find after some computations that [41, 42]

$$ (3.5) \qquad \mathcal{L}_{\mathcal{M}_{\rho, \mathbf{u}, T}} g^{(1)} = \mathcal{M}_{\rho, \mathbf{u}, T} \left[ \frac{1}{2} \mathbf{A}(\mathbf{V}) : \mathbf{D}(\mathbf{u}) + 2 \mathbf{B}(\mathbf{V}) \cdot \nabla_x \sqrt{T} \right], $$

where $\mathbf{A}$, $\mathbf{B}$ and $\mathbf{V}$ are defined in (3.3) and the traceless deformation tensor $\mathbf{D}$ of $\mathbf{u}$ is given by [42]

$$ \mathbf{D}(\mathbf{u}) = \nabla_x \mathbf{u} + (\nabla_x \mathbf{u})^\top - \frac{2}{3} (\mathrm{div}_x \mathbf{u}) \mathbf{I}. $$

Moreover, using the hypothesis $(\mathbf{H}_1)$ on the conservation laws of the collision operator, it is possible to show that linear combinations of collisional invariants form exactly the kernel of the linear operator $\mathcal{L}_{\mathcal{M}_{\rho, \mathbf{u}, T}}$. In particular, we have the orthonormal family

$$ \ker \mathcal{L}_{\mathcal{M}_{\rho, \mathbf{u}, T}} = \mathrm{Span} \left\{ \frac{1}{\rho}, \frac{\mathbf{V}}{\rho}, \frac{1}{2\rho} (|\mathbf{V}|^2 - 3) \right\}. $$

Using the orthogonality properties of the moments of a Maxwellian distribution, we have that on $L^2(\mathcal{M}_{\rho, \mathbf{u}, T})$

$$ \mathbf{A}(\mathbf{V}), \mathbf{B}(\mathbf{V}) \perp \ker \mathcal{L}_{\mathcal{M}_{\rho, \mathbf{u}, T}}. $$

Since the operator $\mathcal{L}_{\mathcal{M}_{\rho, \mathbf{u}, T}}$ is invertible on the orthogonal of its kernel and finally using (3.5), it yields

$$ (3.6) \qquad g^{(1)} = \frac{1}{2} \mathcal{L}_{\mathcal{M}_{\rho, \mathbf{u}, T}}^{-1} (\mathcal{M}_{\rho, \mathbf{u}, T} \mathbf{A}) : \mathbf{D}(\mathbf{u}) + 2 \mathcal{L}_{\mathcal{M}_{\rho, \mathbf{u}, T}}^{-1} (\mathcal{M}_{\rho, \mathbf{u}, T} \mathbf{B}) \cdot \nabla_x \sqrt{T}. $$

We can then plug this expression into the definition (3.3) to obtain using some classical symmetry properties of the collision operator [42] that

$$ (3.7) \qquad \left\{ \begin{array}{l} \bar{A}_{NS}^{\epsilon} := \frac{1}{\rho} \int_{\mathbb{R}^3} \mathbf{A}(\mathbf{V}) \mathcal{M}_{\rho, \mathbf{u}, T}(v) [1 + \epsilon g^{(1)}(v)] dv = -\epsilon \frac{\mu}{\rho T} \mathbf{D}(\mathbf{u}), \\ \\ \bar{B}_{NS}^{\epsilon} := \frac{1}{\rho} \int_{\mathbb{R}^3} \mathbf{B}(\mathbf{V}) \mathcal{M}_{\rho, \mathbf{u}, T}(v) [1 + \epsilon g^{(1)}(v)] dv = -\epsilon \frac{\kappa}{\rho T^{3/2}} \nabla_x T. \end{array} \right. $$

<footnote n="1">¹Namely the Frechet derivative of the collision operator.</footnote>

<footnote n="2">²The Double-Dot product for matrices is the Frobenius inner product, and it is calculated computing the Hadamard product and adding up all the elements of the resulting matrix.</footnote>
The scalar quantities $\mu$ and $\kappa$ in (3.7), respectively the viscosity and the thermal conductivity, are given by [41, 42]

$$ \mu := \frac{T}{10} \int_{\mathbb{R}^3} \mathbf{A}(\mathbf{V}) : \mathcal{M}_{\rho, \mathbf{u}, T}(v) \mathcal{L}_{\mathcal{M}_{\rho, \mathbf{u}, T}}^{-1} (\mathcal{M}_{\rho, \mathbf{u}, T} \mathbf{A}) (v) \, dv, $$

$$ \kappa := \frac{T}{3} \int_{\mathbb{R}^3} \mathbf{B}(\mathbf{V}) \cdot \mathcal{M}_{\rho, \mathbf{u}, T}(v) \mathcal{L}_{\mathcal{M}_{\rho, \mathbf{u}, T}}^{-1} (\mathcal{M}_{\rho, \mathbf{u}, T} \mathbf{B}) (v) \, dv. $$

They depend on the collision kernel of the model. For example, for the Boltzmann operator in the 3D hard-sphere case, it can be shown [26] that there exists some positive constants $\mu_0$, $\kappa_0$ such that

$$ \mu = \mu_0 \sqrt{T} \quad \text{and} \quad \kappa = \kappa_0 \sqrt{T}. $$

Finally, the evolution of the macroscopic quantities at first order with respect to $\varepsilon$ is given by the compressible Navier-Stokes equations

$$ (3.8) \qquad \left\{ \begin{aligned} & \partial_t \rho + \mathrm{div}_x (\rho \mathbf{u}) = 0, \\ & \partial_t (\rho \mathbf{u}) + \mathrm{div}_x (\rho \mathbf{u} \otimes \mathbf{u} + \rho T \mathbf{I}) = \varepsilon \mathrm{div}_x (\mu \mathbf{D}(\mathbf{u})), \\ & \partial_t E + \mathrm{div}_x (\mathbf{u} (E + \rho T)) = \varepsilon \mathrm{div}_x (\mu \mathbf{D}(\mathbf{u}) \cdot \mathbf{u} + \kappa \nabla_x T). \end{aligned} \right. $$

3.1.3. Second order: Burnett equations. Pushing the expansion (3.1) at second order in $\varepsilon$, we can use the same type of argument that for the compressible Navier-Stokes system to obtain another correction of the compressible Euler equations: the Burnet system. Although this system is ill-posed [26], the computation of its coefficients is still possible. We have for the BGK case [22, 53]

$$ (3.9) \qquad \begin{aligned} \bar{\mathbf{A}}_{\text{Burnett}}^{\varepsilon} := & \frac{1}{\rho} \int_{\mathbb{R}^3} \mathbf{A}(\mathbf{V}) \mathcal{M}_{\rho, \mathbf{u}, T}(v) [1 + \varepsilon g^{(1)}(v) + \varepsilon^2 g^{(2)}(v)] \, dv \\ = & -\varepsilon \frac{\mu}{\rho T} \mathbf{D}(\mathbf{u}) - 2\varepsilon^2 \frac{\mu^2}{\rho^2 T^2} \left\{ \begin{aligned} & -\frac{T}{\rho} \mathrm{Hess}_x(\rho) + \frac{T}{\rho^2} \nabla_x \rho \otimes \nabla_x \rho - \frac{1}{\rho} \nabla_x T \otimes \nabla_x \rho \\ & + (\nabla_x \mathbf{u}) (\nabla_x \mathbf{u})^{\top} - \frac{1}{3} \mathbf{D}(\mathbf{u}) \mathrm{div}_x(\mathbf{u}) + \frac{1}{T} \nabla_x T \otimes \nabla_x T \end{aligned} \right\}; \end{aligned} $$

$$ (3.10) \qquad \begin{aligned} \bar{\mathbf{B}}_{\text{Burnett}}^{\varepsilon} := & \frac{1}{\rho} \int_{\mathbb{R}^3} \mathbf{B}(\mathbf{V}) \mathcal{M}_{\rho, \mathbf{u}, T}(v) [1 + \varepsilon g^{(1)}(v) + \varepsilon^2 g^{(2)}(v)] \, dv \\ = & -\varepsilon \frac{\kappa}{\rho T^{3/2}} \nabla_x T - \varepsilon^2 \frac{\mu^2}{\rho^2 T^{5/2}} \left\{ \begin{aligned} & + \frac{25}{6} (\mathrm{div}_x \mathbf{u}) \nabla_x T \\ & - \frac{5}{3} [T \mathrm{div}_x (\nabla_x \mathbf{u}) + (\mathrm{div}_x \mathbf{u}) \nabla_x T + 6 (\nabla_x \mathbf{u}) \nabla_x T] \\ & + \frac{2}{\rho} \mathbf{D}(\mathbf{u}) \nabla_x (\rho T) + 2T \mathrm{div}_x (\mathbf{D}(\mathbf{u})) + 16 \mathbf{D}(\mathbf{u}) \nabla_x T \end{aligned} \right\}. \end{aligned} $$

3.2. The moment realizability criterion. The matrix $\bar{\mathbf{A}}^{\varepsilon}$ and the vector $\bar{\mathbf{B}}^{\varepsilon}$ will allow us to define our hydrodynamic break down criterion. Let us set the vector of the reduced collisional invariants for $\mathbf{V} = (v - \mathbf{u}^{\varepsilon}) / \sqrt{T^{\varepsilon}}$,

$$ \mathbf{m} := \left( 1, \mathbf{V}, \left(\frac{2}{5}\right)^{1/2} \left(\frac{|\mathbf{V}|^2}{2} - \frac{5}{2}\right) \right). $$

Similarly to [42, 55], now we define the so-called moment realizability matrix $\mathbf{M}$ by setting

$$ (3.11) \qquad \mathbf{M} := \frac{1}{\rho^{\varepsilon}} \int_{\mathbb{R}^3} (\mathbf{m} \otimes \mathbf{m}) f^{\varepsilon}(v) \, dv. $$
Indeed, if the distribution function is nonnegative, then the quantity $\int \phi^2(v) f \, dv$ is also nonnegative [55], for any function $\phi(v)$. The same is true if we consider $\mathbf{m} = \mathbf{m}(v)$ a columns vector of $c$ polynomials, and we set $\phi = \mathbf{a}^\top \mathbf{m}$ for a generic $\mathbf{a} \in \mathbb{R}^c$, and we have

$$ (3.12) \qquad \mathbf{a}^\top \left( \int_{\mathbb{R}^3} \mathbf{m} \otimes \mathbf{m} f \, dv \right) \mathbf{a} = \int_{\mathbb{R}^3} (\mathbf{a}^\top \mathbf{m})^2 f \, dv \geq 0. $$

Then if $f$ is nonnegative and not identically zero, then for every $\mathbf{m}$ the $c \times c$ matrix

$$ (3.13) \qquad \mathcal{M} = \frac{1}{\rho} \int_{\mathbb{R}^3} \mathbf{m} \otimes \mathbf{m} f \, dv, $$

is positive semi-definite.

In the Appendix A we prove that under the assumption of being sufficiently close enough to the equilibrium distribution function (i.e. the Maxwellian) then the quantity $\mathcal{M}$ is positive semi-definite if and only if the following quantity is also positive semi-definite:

$$ (3.14) \qquad \mathcal{V}_{\varepsilon} = \mathbf{I} + \bar{\mathbf{A}}^{\varepsilon} - \frac{2}{3} \bar{\mathbf{B}}^{\varepsilon} \otimes \bar{\mathbf{B}}^{\varepsilon}. $$

Then the idea is to monitor the evolution of the matrix $\mathcal{V}_{\varepsilon}$ [42, 54]: indeed the moment realizability criterion states that the fluid dynamic description has broken down when the perturbation is too large that $\mathcal{V}_{\varepsilon}$, and thus $\mathcal{M}$, is no longer positive semi-definite. More details on how to build a hierarchical solver are discussed in the next section.

### 3.3. Three layers hierarchical model.
Now we present a hybrid decomposition numerical method able to automatically switch between three different layers: Euler equations (Layer 0), kinetic model with ES-BGK operator (Layer 1) and kinetic model with Boltzmann operator (Layer 2). We decided to not implement the Navier-Stokes equation because the CFL condition would have required an extremely small timestep.

<table><thead><tr><th>LAYER</th><th>MODEL</th><th>Computational Cost & Accuracy</th></tr></thead><tbody><tr><td>2</td><td>Kinetic model with Boltzmann operator</td><td>Most accurate & highest computational cost</td></tr><tr><td>1</td><td>Kinetic model with ES-BGK operator</td><td>Most accurate & highest computational cost</td></tr><tr><td>0</td><td>Euler equations</td><td>Least accurate & lowest computational cost</td></tr></tbody></table>

Now from the remarks of the previous section, let us define a criterion to determine the appropriate model – fluid or kinetic – to be used. The key idea is the following: we evaluate the quantity $\mathcal{V}_{\varepsilon}$ defined in Equation 3.14 for three orders in $\varepsilon$ of the Chapman-Enskog expansion, and, since the greater the order of expansion, then the greater the accuracy of the model, we associate to each expansion order a different model. In particular, to 0-th order in the Chapman-Enskog expansion we associate the Euler equations (lowest accuracy), to 1-st order we associate the ES-BGK equation, and to 2-nd order we associate the full Boltzmann equation (highest accuracy). Then, we will have three different matrices accordingly to the expansion order in $\varepsilon$, that are: $\mathcal{V}_{\varepsilon^0}$ (associated to 0-th expansion order, and to Euler equations), $\mathcal{V}_{\varepsilon^1}$ (associated to 1-st expansion order, ES-BGK equation) and $\mathcal{V}_{\varepsilon^2}$ (associated to 2-nd expansion order, full Boltzmann equation). So, in each cell of our spatial domain we can evaluate these three matrices, and we compare them in order to decide if we have to switch from one level to another one. If in a given cell the matrix $\mathcal{V}_{\varepsilon^i}$ is too different from matrix $\mathcal{V}_{\varepsilon^{i-1}}$, then it means that the layer $i-1$ is not accurate, and then we have to use the layer $i$. In order to establish if the previous matrices are too different, we will compare their eigenvalues. This is the main idea behind our hybrid domain decomposition numerical scheme.

Let's consider the zeroth order Chapman-Enskog expansion with respect to $\varepsilon$: we get that $\bar{\mathbf{A}}_0 = \mathbf{0}_{M_3}$ and $\bar{\mathbf{B}}_0 = \mathbf{0}_{\mathbb{R}^3}$ and then

$$ \mathcal{V}_{\varepsilon^0} = \mathbf{I}. $$
On the other hand, consider the first order model. By cutting the Chapman-Enskog expansion (3.1) at the first order with respect to $\epsilon$, we can compute explicitly the matrix $V_{\epsilon^1}$ obtaining

$$ (3.15) \qquad V_{\epsilon^1} = \mathbf{I} - \epsilon \frac{\mu}{\rho T} \mathbf{D}(\mathbf{u}) - \epsilon^2 \frac{2}{3} \frac{\kappa^2}{\rho^2 T^3} \nabla_x T \otimes \nabla_x T, $$

where $(\rho, \mathbf{u}, T)$ are solution to the Navier-Stokes equations (3.8).

Hence, we claim that the layer 0 (Euler equations) is correct when the matrix $V_{\epsilon^0}$ behaves like the matrix $V_{\epsilon^0} = \mathbf{I}$, that is, it is positive definite and if its eigenvalues are close to 1 or not: the Layer 0 description of the fluid will be considered incorrect if the difference $\lambda_A$ between the eigenvalues $\lambda_{\epsilon^1}$ of $V_{\epsilon^1}$ and the eigenvalues $\lambda_{\epsilon^0} = 1$ of $V_{\epsilon^0}$ is too big. Then our criterion to switch from Euler description (Layer 0) to ES-BGK (Layer 1) is the following inequality

$$ (3.16) \qquad |\lambda_A| = |\lambda_{\epsilon^1} - 1| > \eta_0, \quad \forall \lambda_{\epsilon^1} \in \text{Sp}(V_{\epsilon^1}), $$

where $\eta_0$ is a small parameter (and it is a problem dependent quantity); in fact if this conditions it is satisfied then it means that the eigenvalues associated to the first order of the expansions are too far from those associated to the zeroth order, and so the Layer 0 is not accurate. Then we decide to use this criterion in our code to automatically switch between Euler equations (Layer 0) and ES-BGK equation (Layer 1).

For the switching criterion between the Layer 1 (ES-BGK equation) and Layer 2 (Boltzmann equation), we follow the same strategy as before. So we evaluate $V_{\epsilon^1}$ and $V_{\epsilon^2}$, we compute the respective eigenvalues, and we compare them. If they are too different, then we switch from Layer 1 to Layer 2. In the paper [23], the authors evaluated the dominant differences between the matrices associated to $V_{\epsilon^1}$ and $V_{\epsilon^2}$. In their work, the $V_{\epsilon^1}$ was associated to Navier-Stokes equations, and $V_{\epsilon^2}$ to the ES-BGK equation, while in our case, they are associated respectively to ES-BGK and Boltzmann equations. We report now the main differences between the matrices of the two levels in the Chapman-Enskog expansion, as reported in [23]

$$ (3.17) \qquad \lambda_B(t, \mathbf{x}) := \epsilon^2 \left( \frac{|\nabla_x T|^2}{T} + |\nabla_x \mathbf{u}|^2 + \sqrt{\left( |\Delta_x \mathbf{u}|^2 + \left| \frac{\Delta_x \rho}{\rho} \right|^2 \right) (1 + T^2)} \right). $$

We will use this quantity as indicator to switch from the ES-BGK equation to the full Boltzmann equation. In particular, we will say that the ES-BGK equation is not appropriate at some point $(t, \mathbf{x})$ if:

$$ (3.18) \qquad |\lambda_B(t, \mathbf{x})| > \eta_1, $$

where $\eta_1$ is a small parameter (and it is a problem dependent quantity). If this inequality is satisfied then we switch to the Boltzmann equation.

**Remark 1.** Note that this last two criteria must be also verified when going from Layer 2 to Layer 1, and from Layer 1 to Layer 0, respectively. In particular, for moving from Layer 2 the following condition must be satisfied

$$ (3.19) \qquad |\lambda_B(t, \mathbf{x})| \le \eta_1. $$

In our code this last condition is sufficient to switch from Layer 2 to Layer 1. Analogously, for moving from Layer 1 to Layer 0 the following condition must be satisfied

$$ (3.20) \qquad |\lambda_A| = |\lambda_{\epsilon^1}(t, \mathbf{x}) - 1| \le \eta_0. $$

This last condition is necessary but not sufficient to switch from Layer 1 to Layer 0, as we explain in the following.
In fact, knowing the full kinetic description of a gas (Layer 1 or Layer 2), there exists several methods [51] to decide how far this gas is from the thermal equilibrium, i.e. the fluid regime (Layer 0). In our code we compare $f^\epsilon$, solution to the collisional equation (2.1) with the truncated Chapman-Enskog distribution $f_k^\epsilon$ given by (3.1), whose moments match the one of $f^\epsilon$, and whose order $k$ corresponds to the order of the macroscopic model considered.

Our criterion is then the following: the kinetic description given by the ES-BGK equation (Layer 1) at point $(t, x)$ corresponds to an hydrodynamic closure of order 0 (Layer 0) if

$$ (3.21) \qquad \|f^\epsilon(t, x, \cdot) - M_{\rho, u, T}^\epsilon(t, x, \cdot)\|_{L^1} \le \delta_0, $$

where $\delta_0$ is a small parameter (which is a problem dependent quantity) and if the eigenvalues of the moment realizability matrix $V_{\epsilon 1}$ computed using the moments of $f^\epsilon$ does not verify the criterion (3.16).

## 4. NUMERICAL SCHEMES

### 4.1. Spatial, velocity and time settings.
We consider a 2D spatial cartesian grid, of intervals $[0, L_x]$ and $[0, L_y]$, and discretized with $N_x$ and $N_y$ points in the $x$ and $y$ direction respectively.

The spatial cartesian grid is defined by nodes $(x_i, y_j) = (i dx, j dy)$ and cells

$$ (4.1) \qquad I_{ij} = \left(x_{i-\frac{1}{2}}, x_{i+\frac{1}{2}}\right) \times \left(y_{j-\frac{1}{2}}, y_{j+\frac{1}{2}}\right). $$

The discrete velocity grid which represents the set of admissible velocities in the discrete model. Let $\mathcal{K} \subseteq \mathbb{Z}^3$ be a set of $N_v$ integer vectors, and let

$$ (4.2) \qquad \nu = \{v_k \in \mathbb{R}^3, k \in \mathcal{K}\}, $$

be a discrete-velocity grid of $N_v$ points indexed by $k = (k_1, k_2, k_3) \in \mathcal{K}$, and defined by

$$ (4.3) \qquad v_k = (v_{k_1}, v_{k_2}, v_{k_3}) = (k_1 \Delta v_1, k_2 \Delta v_2, k_3 \Delta v_3), $$

where $\Delta v_1, \Delta v_2, \Delta v_3$ are three positive numbers characterising the size of the velocity mesh. To simplify notations we denote by $\Delta v_{\mathcal{K}} = \Delta v_1 \Delta v_2 \Delta v_3$. The velocity distribution $f$ (Layer 1 or Layer 2) is then replaced by a vector $\mathbf{f} = (f_k(t, x, y))_{k \in \mathcal{K}} \in \mathbb{R}^{N_x}$ where each component $f_k(t, x, y)$ is assumed to be an approximation of $f(t, x, y, v_k)$.

Each cell $(x_i, y_j)$ is in one of three following regimes: Layer 0 is hydrodynamic (described by Euler equations), Layer 1 is the ES-BGK equation, or Layer 2 is the full Boltzmann equation. We identify with the symbol $r_{ij}^n = \{0, 1, 2\}$ the regime of the cell $(x_i, y_j)$ at the timestep $n$. If $r_{ij}^n = 2$ then the fluid in the cell $(x_i, y_j)$, at time $t^n$, is described by a distribution function $f_k^2(t^n, x_i, y_j)$ which is updated using the spectral method for the Boltzmann operator; if $r_{ij}^n = 1$ then the fluid in the cell $(x_i, y_j)$ is described by a distribution function $f_k^1(t^n, x_i, y_j)$ which is updated using the ES-BGK operator; if $r_{ij}^n = 0$ then the fluid in the cell $(x_i, y_j)$ is described by the moments $U_{t^n, x_i, y_j} = (\rho, \vec{u}, T)_{ij}$ and they are updated using the Euler equations. To each cell we can associate a Maxwellian $M_{\rho, \vec{u}, T}^n(t^n, x_i, y_j)$, computed using the moments $\vec{m}(t^n, x_i, y_j) = [\rho, \vec{u}, T] (t^n, x_i, y_j)$ (density, velocity and temperature) obtained by integrating the distribution function in the cell itself.

Our implementation is characterised by one single time step $\Delta t = 0.1 dx$, used to update both the hydrodynamic solution, and the kinetic ones (ES-BGK and Boltzmann operator), where $dx$ is the spatial discretization (in our simulations we use the same discretization in both spatial directions, thus $dx = dy$). The time integration, discussed in further detail in the next sections, is performed using a TVD Runge Kutta 2 [28, 52] for the hydrodynamic regime and the so-called IMEX PR(2,2,2) [48] numerical scheme for the kinetic levels.
4.2. **Euler equations - Layer 0.** In this subsection, we shall focus on the space and time discretization of the system of N conservation laws

$$ (4.4) \qquad \left\{ \begin{array}{l} \displaystyle \frac{\partial U}{\partial t} + \operatorname{div}_x F(U) = 0, \quad \forall (t, x) \in \mathbb{R}_+ \times \Omega, \\[1.5ex] U(0, x) = U_0(x), \end{array} \right. $$

for a smooth function $F: \mathbb{R}^N \to M_{N \times d_x}(\mathbb{R})$ and a Lipschitz-continuous domain $\Omega \subset \mathbb{R}^{d_x}$.

Here we apply a finite volume scheme, where the reconstruction of the variables at the interfaces is performed using a CWENO scheme [43] and the flux is reconstructed using the Lax-Friedrich flux [49]. In particular, we use the 1D reconstruction in both the x and y directions<fnref n="3" />.

The Euler equations in a 2D space can be expressed using the following conservation laws

$$ (4.5) \qquad \frac{\partial U}{\partial t} + \frac{\partial \mathbf{h}}{\partial x} + \frac{\partial \mathbf{g}}{\partial y} = 0, $$

where

$$ (4.6) \qquad U = \begin{bmatrix} \rho \\ \rho u_x \\ \rho u_y \\ E \end{bmatrix}, \quad \mathbf{h} = \begin{bmatrix} \rho u_x \\ \rho u_x^2 + P \\ \rho u_x u_y \\ u_x(E + P) \end{bmatrix}, \quad \mathbf{g} = \begin{bmatrix} \rho u_y \\ \rho u_y u_x \\ \rho u_y^2 + P \\ u_y(E + P) \end{bmatrix}, $$

and $P = (\Xi - 1)(E - \frac{1}{2}\rho(u_x^2 + u_y^2))$ is the pressure, $\rho, u_x, u_y, E$ and $\Xi$ are, respectively, the density, velocity along x-direction, velocity along y-direction, energy and specific heat ratio.

The time integration is performed using the second order TVD Runge Kutta scheme presented in [28, 52], which is the following one

$$ (4.7) \qquad U^{(1)} = U^n + \Delta t L(U^n), $$

$$ (4.8) \qquad U^{n+1} = \frac{1}{2} U^n + \frac{1}{2} U^{(1)} + \frac{1}{2} \Delta t L(U^{(1)}), $$

where, in our case, $L(U) = -\frac{\partial \mathbf{h}}{\partial x} - \frac{\partial \mathbf{g}}{\partial y}$.

4.3. **ES-BGK model - Layer 1.** We now focus briefly on the evolution of the ES-BGK equation

$$ (4.9) \qquad \left\{ \begin{array}{l} \displaystyle \frac{\partial f^{\varepsilon}}{\partial t} + v \cdot \nabla_x f^{\varepsilon} = \frac{\nu}{\varepsilon} (G[f] - f), \\[1.5ex] f^{\varepsilon}(0, x, v) = f_0(x, v), \end{array} \right. $$

Here we will consider $f = f^1$, i.e. the distribution function associated to the Layer 1, which is updated using the ES-BGK operator. In the following we denote by $f_{ij,k}^n = f_k(t^n, x_i, y_j), \forall k \in \mathcal{K}$ and by $M_{ij,k}^n[f] = M_k[f](t^n, x_i, y_j)$ the discrete Maxwellian equilibrium. We further denote with $\nu_{ij}^n = \nu(\rho_{ij}^n, T_{ij}^n)$ the collision frequency.

We numerically integrate the ES-BGK equation using the same IMEX strategy presented in [20]. The advantage of such a time discretization is that it is uniformly stable with respect to small Knudsen number, thus removing the stiffness of the relaxation term, yet the implicit relaxation term can be solved explicitly, thanks to the special structure of the relaxation term. In addition, the authors have proved the asymptotic-preserving property of this scheme, which allows it to capture the fluid dynamic

<footnote n="3">³Such an approach is also used to compute the linear free transport part of the kinetic layers.</footnote>
behaviour without resolving numerically the small Knudsen number. Here we now present the IMEX strategy as illustrated in [20], but in our code we integrate in time using the PR(2,2,2) scheme [48], that we will discuss at the end of this section. Applying a first order implicit-explicit IMEX scheme to the ES-BGK equation one obtains [20]

$$ (4.10) \qquad \left\{ \begin{array}{l} \displaystyle \frac{f^{n+1} - f^n}{\Delta t} + v \cdot \nabla_x f^n = \frac{\nu^{n+1}}{\epsilon} (G[f^{n+1}] - f^{n+1}), \\ f^0(x, v) = f_0(x, v). \end{array} \right. $$

After some manipulations, one can obtain the following scheme [20]

$$ (4.11) \qquad \left\{ \begin{aligned} U^{n+1} &= \int_{\mathbb{R}^{d_v}} \phi(v) (f^n - \Delta t v \cdot \nabla_x f^n) dv, \\ \Sigma^{n+1} &= \frac{\epsilon}{\epsilon + (1-\beta)\nu^{n+1}\Delta t} \left( \Sigma^n - \Delta t \int_{\mathbb{R}^{d_v}} v \otimes v v \cdot \nabla_x f^n dv \right) + \\ & \qquad \frac{(1-\beta)\nu^{n+1}\Delta t}{\epsilon + (1-\beta)\nu^{n+1}\Delta t} \rho^{n+1} (T^{n+1}I + u^{n+1} \otimes u^{n+1}), \\ f^{n+1} &= \frac{\epsilon}{\epsilon + \nu^{n+1}\Delta t} [f^n - \Delta t v \cdot \nabla_x f^n] + \frac{\nu^{n+1}\Delta t}{\epsilon + \nu^{n+1}\Delta t} G[f^{n+1}], \end{aligned} \right. $$

where $\phi(v) := (1, v, |v|^2)$ (collisional invariants),

$$ (4.12) \qquad \Sigma^{n+1} = \int_{\mathbb{R}^{d_v}} v \otimes v f^{n+1} dv = \rho^{n+1}(\Theta^{n+1} + u^{n+1} \otimes u^{n+1}), $$

and

$$ (4.13) \qquad \Theta(t, x) = \frac{1}{\rho} \int_{\mathbb{R}^3} (v - u) \otimes (v - u) f(t, x, v) dv. $$

Thus, even if (4.10) is nonlinearly implicit, it can be solved explicitly.

The system (4.10) is a first order IMEX scheme, however we used the same strategy to implement the PR(2,2,2) second order IMEX [48], which is adoptable for any stiff systems of differential equation in the form

$$ (4.14) \qquad y' = q(y) + \frac{1}{\epsilon} g(y), $$

where $y = y(t) \in \mathbb{R}^N$, $q, g: \mathbb{R}^N \to \mathbb{R}^N$. In our case $y$ is the distribution function, when considering only the temporal dependence, $q(y)$ is the transport part of the Boltzmann equation, and $g(y)$ is the ES-BGK operator.

An Implicit-Explicit IMEX Runge Kutta scheme for system (4.14) is of the form

$$ (4.15) \qquad \left\{ \begin{aligned} Y_i &= y_n + \Delta t \sum_{j=1}^{i-1} \tilde{a}_{ij} q(t_n + \tilde{c}_j \Delta t, Y_j) + \Delta t \sum_{j=1}^{\Lambda} a_{ij} \frac{1}{\epsilon} g(t_n + c_j \Delta t, Y_j), \\ y_{n+1} &= y_n + \Delta t \sum_{i=1}^{\Lambda} \tilde{\omega}_i q(t_n + \tilde{c}_i \Delta t, Y_i) + \Delta t \sum_{i=1}^{\Lambda} \omega_i \frac{1}{\epsilon} g(t_n + c_i \Delta t, Y_i). \end{aligned} \right. $$

The matrices $\tilde{A} = (\tilde{a}_{ij})$, $\tilde{a}_{ij} = 0$ for $j \ge i$ and $A = (a_{ij})$ are $\Lambda$ matrices such that the resulting scheme is explicit in $q$, and implicit in $g$. An IMEX Runge-Kutta scheme is characterized by these two matrices and the coefficient vectors $\tilde{c} = (\tilde{c}_1, ..., \tilde{c}_\Lambda)^\mathrm{T}$, $\tilde{\omega} = (\tilde{\omega}_1, ..., \tilde{\omega}_\Lambda)^\mathrm{T}$, $c = (c_1, ..., c_\Lambda)^\mathrm{T}$, $\omega = (\omega_1, ..., \omega_\Lambda)^\mathrm{T}$. They can be represented by a double tableau in the usual Butcher notation,

$$ \begin{array}{c|c} \tilde{c} & \tilde{A} \\ \hline \tilde{\omega}^\mathrm{T} & c \end{array}, \qquad \begin{array}{c|c} A & \tilde{\omega}^\mathrm{T} \\ \hline \omega^\mathrm{T} & c \end{array}. $$
The PR(2,2,2) [48] is characterized by the following Butcher tableaus

<table><thead><tr><th>0</th><th>0</th><th>0</th><th>1-C</th><th>1-C</th><th>0</th></tr></thead><tbody><tr><td>1</td><td>1</td><td>0</td><td>C</td><td>C-δ</td><td>δ</td></tr><tr><td colspan="2">1/2</td><td>1/2</td><td colspan="2">1/2</td><td>1/2</td></tr></tbody></table>

where $\delta = 1 - 1/(2C)$. In our scheme we use $C = 1/\sqrt{2}$.

Now it remains only to discuss the spatial integration of the transport part in the Boltzmann equation, which is performed using the same CWENO scheme [43,49] used for the hydrodynamic setting.

**4.4. Boltzmann equation - Layer 2.** We now focus on the evolution of the kinetic equation where we use the Boltzmann operator

$$ (4.16) \qquad \begin{cases} \frac{\partial f^{\epsilon}}{\partial t} + v \cdot \nabla_x f^{\epsilon} = \frac{1}{\epsilon} Q(f^{\epsilon}), \\ f^{\epsilon}(0, x, v) = f_0(x, v). \end{cases} $$

Here we consider $f = f^2$, i.e. the distribution function associated to the Layer 2. We shall adopt the strategy proposed in [19]. When the collision operator $Q$ is the BGK operator, it is well known that even an implicit collision term can be solved explicitly, using the property that $Q$ preserves mass, momentum and energy. The main idea in [19] is to use this property penalizing the Boltzmann collision operator $Q$ by the BGK operator as follows

$$ (4.17) \qquad Q = [Q - \lambda(M - f)] + \lambda[M - f], $$

where $M$ is the Maxwellian associated to the distribution function $f$ and $\lambda$ is the spectral radius of the linearized collision operator of $Q$ around $M$.

In this last equation, the first term on the right-hand side is less stiff (or not stiff at all) compared to the second one. So the first term can be discretized explicitly, thus avoiding to invert the nonlinear Boltzmann operator $Q$, and the second term can be treated explicitly.

In [19] the authors propose the following split procedure in the equation 4.16

$$ (4.18) \qquad \frac{Q(f)}{\epsilon} = \frac{Q(f) - P(f)}{\epsilon} + \frac{P(f)}{\epsilon}, $$

where the first term on the RHS is less stiff than $\frac{P(f)}{\epsilon}$. The quantity $P(f)$ is a linear operator asymptotically close to the source term $Q(f)$, and that preserves the steady state $P(M) = 0$.

It is possible to perform a Taylor expansion leading to

$$ (4.19) \qquad Q(f) = Q(M) + \nabla Q(M)(f - M) + O(||f - M||_H^2), $$

and thus we can choose

$$ (4.20) \qquad P(f) := \nabla Q(M)(f - M). $$

However, since it is not always possible to compute $\nabla Q(M)$ we can introduce an upper bound (or some approximation) $\beta$ of $||\nabla Q(M)||$.

At this point we have split the Boltzmann collision operator $Q$ in a stiff and in a non-stiff term. Thus, we can adopt an IMEX approach like we have done in the previous section for the ES-BGK operator, where in this case, the explicit part is constituted by

$$ (4.21) \qquad -\nabla_x f + \frac{Q(f) - P(f)}{\epsilon}, $$
where $P(f) = \beta(f - \mathcal{M})$; while the implicit term is

$$ (4.22) \qquad \frac{P(f)}{\epsilon}. $$

A first order IMEX scheme as presented in [19] takes the following form

$$ (4.23) \qquad \begin{cases} \frac{f^{n+1} - f^n}{\Delta t} + v \cdot \nabla_x f^n = \frac{\mathcal{Q}(f^n) - P(f^n)}{\epsilon} + \frac{P(f^{n+1})}{\epsilon} \\ f^0(x, v) = f_0(x, v). \end{cases} $$

This system can be rewritten as

$$ (4.24) \qquad f^{n+1} = \frac{\epsilon}{\epsilon + \beta^{n+1}\Delta t} [f^n - \Delta t v \nabla_x f^n] + \Delta t \frac{\mathcal{Q}(f^n) - P(f^n)}{\epsilon + \beta^{n+1}\Delta t} + \frac{\beta^{n+1}\Delta t}{\epsilon + \beta^{n+1}\Delta t} \mathcal{M}^{n+1}, $$

where $\beta^{n+1} = \beta(\rho^{n+1}, T^{n+1})$ and $\mathcal{M}^{n+1}$ is the Maxwellian associated to $f^{n+1}$. In our simulations we considered $\beta = 2\pi\rho$ [19].

The equation (4.24) appears to be nonlinearly implicit, but it can be computed explicitly [19]. Indeed, multiplying it by $\phi(v) := (1, v, |v|^2)$ (collisional invariants) and introducing the macroscopic quantity $U := (\rho, \rho u, T)$ one gets

$$ (4.25) \qquad U^{n+1} = \int \phi(v)(f^n - \Delta t v \cdot \nabla_x f^n) dv. $$

From this last relation, the quantity $U^{n+1}$ can be obtained explicitly from $f^n$, and consequently it is possible to compute $\mathcal{M}^{n+1}$.

Summing up, it is possible to compute explicitly $f^{n+1}$ using the first order IMEX scheme (4.23).

At this point it is possible to generalize this approach to higher order IMEX scheme, and in particular, in our code we have adopted it for the PR(2,2,2) scheme, already discussed in the previous section. The spatial part of the Boltzmann equation has been computed using the already discussed CWENO scheme [43, 49].

### 4.4.1. Fourier-Galerkin spectral methods for the Boltzmann equation.
The aim is to build a numerical method to compute the Boltzmann operator (2.3), which depends only on the velocity variable $v$. We shall then neglects the time and space variables in this section.

The starting point is the truncation of the integration domain of the Boltzmann integral (2.3), then we suppose that the distribution function $f$ has compact support on the ball $B_0(R)$, of radius $R$ centred in the origin. We have the classical result from [47]

**Proposition 4.1.** Let the distribution function $f$ be compactly supported on the ball $B_0(R)$ of radius $R$ centered in the origin, then

$$ (4.26) \qquad \mathrm{Supp}_v(\mathcal{Q}(f, f)) \subset B_0(\sqrt{2}R). $$

It is possible to introduce a collision operator with cut-off in the following way

$$ (4.27) \qquad \mathcal{Q}^R(f,g)(v) = \iint_{B_0(R) \times \mathbb{S}^{d-1}} [f(v')g(v'_*) - f(v)g(v_*)] B(v - v_*, \sigma) dv_* d\sigma, \quad \forall v \in \mathbb{R}^d $$

It is important to have a spectral approximation without superposition of periods. In order to reach this goal, it is sufficient to have a distribution function $f(v)$ restricted on the cube $[-T, T]^d$, with $T \ge (2 + \sqrt{2})R$. Then, we assume $f(v) = 0$ on $[-T, T]^d \setminus B_0(R)$ and extend $f(v)$ to a periodic function on the set $[-T, T]^d$. We want to remark that the lower bound for $T$ can be improved: in fact, the choice $T = (3 + \sqrt{2})R/2$ guarantees the absence of intersection between periods where the distribution function $f$ is different from 0. Unfortunately, the support of $f$ will increase in time, so we can just minimise the errors due to aliasing [50].
In the following we will assume $T = \pi$ and hence $R = 2\pi/(3 + \sqrt{2})$ in order to simplify the notation. Projecting the collision operator onto the space of trigonometric polynomials we obtain

$$ (4.28) \qquad \hat{Q}_k = \int_{[-\pi, \pi]^d} Q^R(f_N, f_N)e^{-ik \cdot v} dv, \quad k = -N, \dots, N $$

where $f_N$ is the truncated Fourier series of the distribution function $f$,

$$ (4.29) \qquad f_N(v) = \sum_{k=-N}^{N} \hat{f}_k e^{ik \cdot v}, $$

$$ \hat{f}_k = \frac{1}{(2\pi)^d} \int_{D_t} f(v)e^{-ik \cdot v} dv, $$

where we have used just one index to denote the $d$-dimension sums with respect to the vector $k = (k_1, \dots, k_d) \in \mathbb{Z}^d$, and so we have

$$ (4.30) \qquad \sum_{k=-N}^{N} := \sum_{k_1, \dots, k_d=-N}^{N}. $$

Now we can substitute the truncated Fourier series $f_N$ inside the expression (4.28), and after some computations we obtain [47]

$$ (4.31) \qquad \hat{Q}_k = \sum_{\substack{l,m=-N \\ l+m=k}}^{N} \hat{f}_l \hat{f}_m \beta(l, m), \quad k = -N, \dots, N, $$

where $\beta(l, m) = \hat{B}(l, m) - \hat{B}(m, m)$, and

$$ (4.32) \qquad \hat{B}(l, m) = \int_{B_0(2\lambda\pi)} \int_{\mathbb{S}^{d-1}} |q|\sigma(|q|, \cos(\theta))e^{-i(l \cdot q^+ + m \cdot q^-)} d\omega dq, $$

with

$$ (4.33) \qquad \lambda = 2/(3 + \sqrt{2}), \quad q^+ = \frac{1}{2}(q + |q|\sigma), \quad q^- = \frac{1}{2}(q - |q|\sigma), \quad q = v - v_* $$

We would like to remark that the naive evaluation of (4.31) requires $O(n^2)$ operations, where $n = N^d$, which means that this spectral method is computationally expensive. The computation can be speed up by reducing the number of operations needed to evaluate the collision integral, introducing a FFT-based convolution approach in our spectral method.

It is possible to find a new representation, called *Carleman-like representation*, of the collision operator which conserves more symmetries when approximated on a bounded domain.

**Lemma 4.1.** *Introducing the change of variables $x = r\sigma/2$ and $y = v_* - v - x$ the collision operator can be rewritten as [46]*

$$ (4.34) \qquad Q(f, f)(v) = \int_{x \in \mathbb{R}^d} \int_{y \in \mathbb{R}^d} \tilde{B}(x, y) \delta(x \cdot y) [f(v+y)f(v+x) - f(v+x+y)f(v)] dx dy, $$

where

$$ (4.35) \qquad \tilde{B}(x, y) = \tilde{B}(|x|, |y|) = 2^{d-1} B \left( \sqrt{|x|^2 + |y|^2}, \frac{|x|}{\sqrt{|x|^2 + |y|^2}} \right) (|x|^2 + |y|^2)^{-\frac{d-2}{2}}, $$

where the last expression of $\tilde{B}(x, y)$ is valid on the manifold defined by $x \cdot y = 0$.
Let us notice that, if we assume the validity of the VHS model (2.4), then we obtain constant $\tilde{B}(x, y)$ if $\gamma = 0$ and $d = 2$ or if $\gamma = 1$ and $d = 3$.

Now let's discuss two possible truncation strategies [50]. We start considering a bounded domain $D_T = [-T, T]^d$ ($0 < T < +\infty$). A first approach is removing the collisions related with some points out of the box. This is the natural first stage for deriving conservative schemes based on the discretization of the velocity. A truncation on the modulus of $x$ and $y$ is not necessary in this case, since we are imposing them to be in the box. So we obtain

$$ (4.36) \quad Q^{tr}(f, f)(v) = \iint_{\{x, y \in \mathbb{R}^d \mid v+x, v+y, v+x+y \in D_T\}} \tilde{B}(x, y)\delta(x \cdot y) [f(v+y)f(v+x) - f(v+x+y)f(v)] dx dy, $$

where $v \in D_T$.

This last operator satisfies the following weak form

$$ (4.37) \quad \int Q^{tr}(f, f)(v)\zeta(v)dv = \frac{1}{4} \iint \int_{\{v, x, y \in \mathbb{R}^d \mid v, v+x, v+y, v+x+y \in D_T\}} \tilde{B}(x, y)\delta(x \cdot y) f(v+x+y)f(v) [\zeta(v+y) + \zeta(v+x) - \zeta(v+x+y) - \zeta(v)] dv dx dy, $$

from which we can deduce the conservation of mass, momentum, energy and the validity of the H-theorem for the entropy. Doing this derivation we have changed the collision kernel by adding artificial dependencies on $v, v_*, v', v'_*$, and so the convolution-like properties are lost. This is the problem of this kind of truncation.

A second possible approach is truncating the integration with respect to variables $x$ and $y$ by setting them to be in $B_0(R)$. In order to obtain all possible collisions we consider $R = 2S$ for a compactly supported function $f$ whose support is $B_0(S)$.

Our goal is to implement the FFT algorithm to evaluate the quadrature approximation, so we will make use of periodic distribution functions, and then we have to properly take into account the aliasing effect due to periods superposition in the Fourier space. A geometrical argument shows that using the periodicity of the function it is enough to take $T \ge (3 + \sqrt{2})S/2$ to prevent intersections of the region where $f$ is different from zero. The operator now can be written as

$$ (4.38) \quad Q^R(f, f)(v) = \int_{x \in B_0(R)} \int_{y \in B_0(R)} \tilde{B}(x, y) \delta(x \cdot y) [f(v+y)f(v+x) - f(v+x+y)f(v)] dx dy, $$

for $v \in D_T$. The main interest of this representation is preserving the real collision kernel and its properties.

The main idea to develop a more efficient algorithm is to use the previous representation, where taking $\zeta(v) = e^{ik \cdot v}$ we obtain the following spectral quadrature formula

$$ (4.39) \quad \hat{Q}_k = \sum_{\substack{l,m=-N \\ l+m=k}}^{N} \hat{f}_l \hat{f}_m \hat{\beta}_F(l,m), \quad k = -N, \dots, N, $$

where $\hat{\beta}_F(l, m) = \hat{B}_F(l, m) - \hat{B}_F(m, m)$ are now given by

$$ (4.40) \quad \hat{B}_F(l, m) = \int_{B_0(R)} \int_{B_0(R)} \tilde{B}(x, y) \delta(x \cdot y) e^{i(l \cdot x + m \cdot y)} dx dy. $$
Finding a convolution structure in the previous equation reduces the number of operations needed to evaluate it. In particular, we want to approximate each $\hat{\beta}_F(l, m)$ as follows

$$ (4.41) \qquad \hat{\beta}_F(l, m) \approx \sum_{p=1}^{A} \alpha_p(l) \alpha'_p(m), $$

where A represents the number of finite possible directions of collisions. In the end, we will obtain a sum of A discrete convolutions, which can be computed in $O(AN \log_2 N)$ operations using the standard FFT technique.

Now we make the decoupling assumption

$$ (4.42) \qquad \tilde{B}(x, y) = a(|x|)b(|y|), $$

in order to obtain the convolution form. We would like to remark that this assumption is satisfied when $\tilde{B}$ is constant.

This is the case of the Maxwellian molecules in dimension two, and hard-spheres in dimension three, the two most used cases. Indeed, putting the variable hard-sphere kernel inside the equation (4.35) we obtain

$$ (4.43) \qquad \tilde{B}(x, y) = 2^{d-1} C_{\gamma} (|x|^2 + |y|^2)^{-\frac{d-\gamma-2}{2}}, $$

and so we obtain that $\tilde{B}$ is constant if $d=2$, $\gamma=0$ (Maxwellian molecules in dimension two) or if $d=3$, $\gamma=1$ (hard-sphere molecules in dimension three).

The 2d case. We start with Maxwellian molecules in dimension two, assuming $\tilde{B} = 1$. We move to spherical coordinate $x = \rho e$ and $x' = \rho' e'$ inside (4.41), and we get

$$ (4.44) \qquad \tilde{B}_F(l, m) = \frac{1}{4} \int_{S^1} \int_{S^1} \delta(e \cdot e') \left[ \int_{-R}^{R} e^{i\rho(l \cdot e)} d\rho \right] \left[ \int_{-R}^{R} e^{i\rho'(m \cdot e')} d\rho' \right] de \, de'. $$

Now we introduce $\phi_R^2(s) = \int_{-R}^R e^{i\rho s} d\rho$, for $s \in \mathbb{R}$, and we obtain the explicit formula

$$ (4.45) \qquad \phi_R^2(s) = 2R \operatorname{sinc}(Rs), $$

where $\operatorname{sinc}(x) = (\sin(x)/x)$. Using the parity properties of the sinc function inside the expression of $\hat{B}_F(l, m)$ we get

$$ (4.46) \qquad \hat{B}_F(l, m) = \int_0^\pi \phi_R^2(l \cdot e_\theta) \phi_R^2(m \cdot e_{\theta+\pi/2}) \, d\theta. $$

Finally a regular discretization of A equally spaced points $\theta_p = \pi p/A$ of the unit sphere, which is spectrally accurate [10] because of the periodicity of the function, gives

$$ (4.47) \qquad \hat{B}_F(l, m) = \frac{\pi}{A} \sum_{p=1}^{A} \alpha_p(l) \alpha'_p(m), $$

with

$$ (4.48) \qquad \alpha_p(l) = \phi_R^2(l \cdot e_{\theta_p}), \qquad \alpha'_p(m) = \phi_R^2(m \cdot e_{\theta_p + \pi/2}). $$

The 3d case. Now we deal with the hard-sphere collision kernel in the three dimensional case. First of all we consider a spherical parametrization $(\theta, \phi)$ of the vector $e \in S_+^2$ and a uniform grid of respective size $A_1$ and $A_2$ for the variables $\theta$ and $\phi$ (spectrally accurate thanks to the periodicity of the function) which leads to the following quadrature for $\tilde{B}_F(l, m)$ [50]

$$ (4.49) \qquad \tilde{B}_F(l, m) = \frac{\pi^2}{A_1 A_2} \sum_{p,q=0}^{A_1, A_2} \alpha_{p,q}(l) \alpha'_{p,q}(m), $$
where

$$ (4.50) \qquad \alpha_{p,q}(l) = \phi_{R,a}^3 (l \cdot e_{(\theta_p, \phi_q)}), \quad \alpha'_{p,q}(m) = \psi_{R,b}^3 \left( \Pi_{e_{\tilde{\theta}_p, \phi_q}} (m) \right), $$

$$ (4.51) \qquad \phi_{R,a}^3(s) = \int_{-R}^{R} \rho a(\rho) e^{i\rho s} d\rho, \quad \psi_{R,b}^3(s) = \int_{0}^{\pi} \sin(\theta) \phi_{R,b}^3(s \cos(\theta)) d\theta, $$

and for all $p$ and $q$

$$ (4.52) \qquad (\theta_p, \phi_q) = \left( \frac{p\pi}{A_1}, \frac{q\pi}{A_2} \right). $$

### 4.5. Evolving the solution and coupling the models.
We are now interested in evolving in time the domain decomposition scheme. At a given time $t^n$, we denote by $K_{ij}$ a control volume, the space domain $\Omega = \Omega_0^n \sqcup \Omega_1^n \sqcup \Omega_2^n$ is decomposed in

* Fluid cells $K_{ij} \subset \Omega_0^n$, described by the hydrodynamic fields

$$ U_{ij}^n := (\rho_{ij}^n, u_{ij}^n, T_{ij}^n) \approx (\rho(t^n, x_i, y_j), u(t^n, x_i, y_j), T(t^n, x_i, y_j)); $$

* Kinetic cells evolved with ES-BGK equation $K_{ij} \subset \Omega_1^n$, described by the particle distribution function $f^1$ evolved using the ES-BGK equation

$$ f_{ij}^{1,n}(v) \approx f^1(t^n, x_i, y_j, v), \quad \forall v \in \mathbb{R}^3. $$

* Kinetic cells evolved with Boltzmann equation $K_{ij} \subset \Omega_2^n$, described by the particle distribution function $f^2$ evolved using the Boltzmann operator by mean of the spectral method

$$ f_{ij}^{2,n}(v) \approx f^2(t^n, x_i, y_j, v), \quad \forall v \in \mathbb{R}^3. $$

Now let's introduce some important quantities. Since we are in a two dimensional space, the matrix $V_{\varepsilon 1}$ is equal to

$$ (4.53) \qquad V_{\varepsilon 1} = \begin{pmatrix} 1 - \xi_1 a_1 - \xi_2 a_2 & -\xi_1 b_1 - \xi_2 b_2 & 0, \\ -\xi_1 b_1 - \xi_2 b_2 & 1 - \xi_1 c_1 - \xi_2 c_2 & 0, \\ 0 & 0 & 1 - \xi_1 h_1, \end{pmatrix} $$

where $a_1 = \frac{4}{3} \frac{\partial u_x}{\partial x} - \frac{2}{3} \frac{\partial u_y}{\partial y}$, $a_2 = (\frac{\partial T}{\partial x})^2$, $b_1 = \frac{\partial u_x}{\partial y} + \frac{\partial u_y}{\partial x}$, $b_2 = (\frac{\partial T}{\partial x} \frac{\partial T}{\partial y})$, $c_1 = \frac{4}{3} \frac{\partial u_y}{\partial y} - \frac{2}{3} \frac{\partial u_x}{\partial x}$, $c_2 = (\frac{\partial T}{\partial y})^2$, $h_1 = -\frac{2}{3} (\frac{\partial u_x}{\partial x} + \frac{\partial u_y}{\partial y})$, $\xi_1 = \frac{\epsilon\mu}{\rho T}$ and $\xi_2 = \epsilon^2 \frac{2k^2}{3\rho^2 T^3}$.

The three eigenvalues are then:

$$ (4.54) \qquad \left\{ \begin{array}{l} \lambda_{\varepsilon 1}^a = \frac{1}{2} \left( \sqrt{[\xi_1 (a_1 - c_1) + \xi_2 (a_2 - c_2)]^2 + 4d^2} - a_1 \xi_1 - a_2 \xi_2 - c_1 \xi_1 - c_2 \xi_2 + 2 \right), \\ \lambda_{\varepsilon 1}^b = \frac{1}{2} \left( -\sqrt{[\xi_1 (a_1 - c_1) + \xi_2 (a_2 - c_2)]^2 + 4d^2} - a_1 \xi_1 - a_2 \xi_2 - c_1 \xi_1 - c_2 \xi_2 + 2 \right), \\ \lambda_{\varepsilon 1}^c = 1 - \xi_1 h_1 \end{array} \right. $$

where $d = \xi_1 b_1 + \xi_2 b_2$.

Now we briefly discuss the computation of the first and second spatial derivative that appear inside the expressions of $\lambda_{\varepsilon 1}^a$, $\lambda_{\varepsilon 1}^b$ and $\lambda_B$ (Equation (3.17)), taking into account that in general there are cells close to each other that are of different type (Layer 0, Layer 1 and Layer 2). For computing these quantities it is necessary to consider a stencil of one cell in positive and negative direction for both $x$ and $y$. A schematic example representation is illustrated in Figure 1. We are interested in
computing the first and second spatial derivative in each direction (x and y) of the following quantities: {$\rho, u_x, u_y, T$} := $\xi$; and $\nabla_x T$, $\Delta_x u$ and $\Delta_x \rho$. Then these quantities have been computed as follows

$$ (4.55) \qquad \left\{ \begin{aligned} \frac{\partial \xi}{\partial x} &= \frac{\xi_{i+1,j} - \xi_{i-1,j}}{2\Delta x}, \\ \frac{\partial \xi}{\partial y} &= \frac{\xi_{i,j+1} - \xi_{i,j-1}}{2\Delta y}, \\ \Delta_x u &= \Delta_x u_x + \Delta_x u_y, \\ \Delta_x \rho &= \frac{\rho_{i-1,j} + \rho_{i+1,j} + \rho_{i,j-1} + \rho_{i,j+1} - 4\rho_{i,j}}{\Delta x \Delta y}, \\ \nabla_x T &= \hat{\mathbf{x}} \frac{\partial T}{\partial x} + \hat{\mathbf{y}} \frac{\partial T}{\partial y}. \end{aligned} \right. $$

where

$$ (4.56) \qquad \left\{ \begin{aligned} \Delta_x u_x &= \frac{u_{x,i-1,j} + u_{x,i+1,j} + u_{x,i,j-1} + u_{x,i,j+1} - 4u_{x,i,j}}{\Delta x \Delta y}, \\ \Delta_x u_y &= \frac{u_{y,i-1,j} + u_{y,i+1,j} + u_{y,i,j-1} + u_{y,i,j+1} - 4u_{y,i,j}}{\Delta x \Delta y}. \end{aligned} \right. $$

It is evident that for applying these differential operators to one of the $\xi$ in the cell $(x_i, y_j)$ it is necessary to know the value of $\xi$ in the cells $(x_i, y_{j-1})$, $(x_i, y_{j+1})$, $(x_{i-1}, y_j)$ and $(x_{i+1}, y_j)$. The code properly estimate the value of $\xi$ in each of these cells, according to its Layer. For example, in Figure 1 the cell $(x_{i+1}, y_j)$ is of type Layer 1, so the code estimate $\xi_{i+1,j}$ as a moment associated to the distribution function $f^1(i+1, j)$; while, in the same Figure, the cell $(x_i, y_{j+1})$ is of type Layer 0, so the quantity $\xi_{i,j+1}$ is obtained directly from the hydrodynamic solution $U(i, j+1)$.

The Algorithm 1 summarises the switching procedure between different regimes, and the following quantities are used: $\lambda_{\varepsilon_1}^k(i, j) = \lambda_{\varepsilon_1}^k(x_i, y_j)$ for $k = \{a, b, c\}$, indicate the three eigenvalues associated to the matrix $V_{\varepsilon_1}$; and $\lambda_B(i, j) = \lambda_B(x_i, y_j)$.

Now we explain what happens between two cells of different types when performing the spatial integration. Let us consider the situation at time $t^n$ depicted in Figure 1. Since the cell $(x_i, y_j)$ is hydrodynamic (Layer 0), the code uses the CWENO TVD RK2 scheme to integrate it. This method requires a stencil in both x and y direction. As we can see from the same Figure, not all the cells of the stencil are of the correct type (hydrodynamic) Layer 0, so in order to use the numerical method we estimate the value of the hydrodynamic solution $U$ in the other cells using the distribution function $f^2$ for the cells of Layer 2 and the distribution function $f^1$ for the cells of Layer 1. For example, in the cells $(i-1, j)$, $(i+1, j)$ and $(i, j-1)$, the hydrodynamic solution is computed using, respectively the moments of the distribution function $f^2(i-1, j)$, $f^1(i+1, j)$ and $f^2(i, j-1)$. We perform the same operation for all cells constituting the stencil.

If now we consider a small variation of the Figure 1, using the Layer 2 in the cell $(i, j)$, then a consideration similar to previous one is still valid. In particular the method that the code uses to evolve in time the distribution functions $f^1$ and $f^2$, uses as well a stencil in both x and y directions. As we can see from the same Figure, not all the cells of the stencil are of the correct type Boltzmann operator (Layer 2), so in order to use the numerical method we estimate the value of the distribution function $f^2$ in the cells $(i, j+1)$, $(i-2, j)$, using the Chapman-Enskog expansion up to the zeroth order computed with the hydrodynamic solutions $U(i, j+1)$ and $U(i-2, j)$; while in the cells $(i+1, j)$, $(i+2, j)$, $(i, j-2)$, we replace the distribution function as follows $f^2(i+1, j) = f^1(i+1, j)$, $f^2(i+2, j) = f^1(i+2, j)$ and
<figure><img src="image_3.png" /><figcaption>FIGURE 1. Example of schematic view of the stencil for the computation of the first and second spatial derivatives, and for the application of a finite volume method.</figcaption></figure>

$f^2(i, j - 2) = f^1(i, j - 2)$. We perform the same operation for all cells constituting the stencil. So after these substitutions we execute the numerical scheme.

It now remains to discuss what happens to the solution in a specific cell of coordinate $(x_i, y_j)$ when its regime changes. We have four possible situations:

* the first situation is a cell that switches from Layer 0 (Euler equations) to Layer 1 (ES-BGK equation). It means that in the previous timestep the cell has been updated using Euler equations, but in the actual timestep it must be updated using the kinetic solver, so we have to associate a distribution function to the cell. In this situations we decide that the distribution function associated to the cell is the Maxwellian (corresponding to 0-th order in the Chapman-Enskog expansion) computed using the moments of the solutions of the Euler equations obtained in the previous timestep using the CWENO TVD RK2 scheme. Better approximation of the distribution function could be obtained using higher order in the Chapman-Enskog expansion (see [22]);

* cell switching from Layer 1 (ES-BGK equation) to Layer 0 (Euler equations). In this situation the solution computed in the previous timestep was the distribution function updated with the ES-BGK equation, but now we have to update it using the Euler equations. So we compute the moments of the distribution function and we use them as input for the CWENO TVD RK2 scheme;

* cell switching from Layer 1 (ES-BGK equation) to Layer 2 (Boltzmann equation). In this situation we decide to use for the Layer 2 the same distribution function that we obtained using the ES-BGK equation in the previous timestep, but of course, in the actual timestep the distribution function is updated using the Boltzmann operator;

* the last situation is a cell switching from Layer 2 (Boltzmann equation) to Layer 1 (ES-BGK equation). Similarly to the previous case, in this configuration we decide to use for the Layer 1
<caption>**Algorithm 1 Regime update**</caption>

1: Actual time: $t^n$
2: **if** $r_{ij}^n == 2$ **then**
3: compute $\lambda_B(i, j)$
4: **if** $|\lambda_B(i, j)| > \eta_1$ **then**
5: $r_{ij}^{n+1} = 2$
6: **else**
7: $r_{ij}^{n+1} = 1$
8: **end if**
9: **else if** $r_{ij}^n == 1$ **then**
10: compute $\lambda_B(i, j), \lambda_{\varepsilon_1}^a(i, j), \lambda_{\varepsilon_1}^b(i, j)$ and $\lambda_{\varepsilon_1}^c(i, j)$
11: **if** $|\lambda_B(i, j)| > \eta_1$ **then**
12: $r_{ij}^{n+1} = 2$
13: **else**
14: **if** $|\lambda_{\varepsilon_1}^a - 1| \le \eta_0$ and $|\lambda_{\varepsilon_1}^b - 1| \le \eta_0$ and $|\lambda_{\varepsilon_1}^c - 1| \le \eta_0$ **then**
15: estimate $\tilde{U}_{ij} = (\rho, \mathbf{u}, T)$ as moments of $f^1(i, j)$
16: compute $\xi = \|f^1(i, j) - \mathcal{M}_{\tilde{U}_{ij}}^n\|_{L^1}$
17: **if** $\xi \le \delta_0$ **then**
18: $r_{ij}^{n+1} = 0$
19: **else**
20: $r_{ij}^{n+1} = 1$
21: **end if**
22: **else**
23: $r_{ij}^{n+1} = 1$
24: **end if**
25: **end if**
26: **else**
27: compute $\lambda_{\varepsilon_1}^a(i, j)$ and $\lambda_{\varepsilon_1}^b(i, j)$ and $\lambda_{\varepsilon_1}^c(i, j)$
28: **if** $|\lambda_{\varepsilon_1}^a - 1| > \eta_0$ or $|\lambda_{\varepsilon_1}^b - 1| > \eta_0$ or $|\lambda_{\varepsilon_1}^c - 1| > \eta_0$ **then**
29: $r_{ij}^{n+1} = 1$
30: **else**
31: $r_{ij}^{n+1} = 0$
32: **end if**
33: **end if**

the same distribution function that we updated using the Boltzmann operator in the previous timestep, but of course, in the actual timestep the distribution function is updated using the ES-BGK equation.

## 5. NUMERICAL SIMULATIONS

The following simulations show the stability of the method implemented and the speed-up offered by the hybrid domain decomposition method compared to using only the Boltzmann operator. The solution computed with the hybrid scheme is compared with the two solutions obtained using, respectively, only the Euler equations and only the Boltzmann operator. The hybrid code (and the full Euler one) is written using C++20 ISO, it is executed only on CPU, and it takes advantage of the multithreads capabilities: in particular in order to speed up the execution, the space domain is divided in subdomains, each of which is assigned to a specific thread to apply the numerical integration method, and all threads are
executed in parallel. All the simulations have been executed on a non-dedicated computer equipped with a Intel (R) Core (TM) Ultra 7 155H CPU, 16.0 GB of RAM and a NVIDIA(R) RTX 1000 Ada Generation Laptop GPU. This shows the great scalability, and numerical optimization of the code. The CPU FFT and IFFT numerical operations required by the Spectral Methods for the Boltzmann operator have been performed using the FFTW library. Since some operations of this library are not multithreads safe, the semaphores [24] have been implemented in order to control the access to shared resources, while the respective GPU operations have been implemented using the cuFFT library.

The solutions computed using only the Boltzmann operator have been obtained using CUDA Programming Language, by massively parallelizing the software. This allows us to highly speed up the execution taking advantage of the intrinsic multithreads capability of the GPUs. Of course, the computation of the kinetic solution using GPUs is faster than the computation of the hybrid solution using CPUs. So in order to show the speed-up offered by the hybrid-CPU solution compared to the full kinetic-CPU solution, we have measured the execution time of a few timesteps of the full kinetic-CPU solution, and then we have estimated the total duration of this integration by multiplying this quantity by the ratio between the total number of time steps and the already computed time steps.

In all the simulations we used the simplest boundary condition implementing the “ghost cells”: on the four rows and four columns closest to the border, we copy the value of the solution and the domain typology of the closest physical cell. In all simulations the velocity space of the distribution functions has 3 dimensions.

In all simulations, except for the last one (Test 5), the regime is initialized everywhere as Euler equations (Layer 0), and then we evolve the system according to our scheme.

In all our simulations we used a collision frequency for the ES-BGK operator equals to $\nu_{ij}^n = \nu(\rho_{ij}^n, T_{ij}^n) = \frac{\pi}{2} \rho_{ij}^n$.

**5.1. Test 1: Sod shock tube.** A classical benchmark test for the 1D Euler equations is the so-called Sod shock tube test problem [58]: its solution consists of a left rarefaction, a contact and a right shock. In the 1D case, on a spatial domain $[0, 1]$ the initial condition is constituted by a left state $U_L = (\rho_L, u_{xL}, P_L) = (1.0, 0.0, 1.0)$ for $x \le 0.5$ and by a right state $U_R = (\rho_R, u_{xR}, P_R) = (0.125, 0.0, 0.1)$ for $x > 0.5$. In this test we simulate these initial conditions, considering also the $y$-direction, along which the solution is always uniform in space. Also the velocity along the $y$-direction is null. In Figure 2 we compare the density and the temperature obtained using the hybrid scheme, with the Euler solution, with the Hard-Sphere solution and the analytic one (associated to the Euler equations). We discretize the spatial domain with 100 points in the $x$-direction, and 16 in the $y$-direction, imposing $\Delta x = \Delta y$. The velocity domain is discretized using 32 points in the three directions, and we considered a cut-off domain $\Omega_v = [-8, 8]^3$. The value of $R$ has been taken approximately equal to 3.3, and both numbers of points $A_1$ and $A_2$ discretizing the spherical integration of the Boltzmann operator (see Section 2.1) have been taken equal to 4.

The distribution function is initialized everywhere as a Maxwellian, whose moments, correspond to the ones of the Euler equations (i.e. $U_L$ and $U_R$), and null velocity in the other directions. We consider a Knudsen number $\epsilon = 10^{-6}$ and a time step $\Delta t = 0.1dx$. It is remarkable that the usage of an IMEX approach for the integration of the ES-BGK equation and of the full Boltzmann one allowed us to use the timestep associated to the transport part of the Boltzmann equation, instead of the one constrained by the collision operator, which in that case, would be much smaller than the one implemented thanks to the IMEX approach. This property significantly accelerates the code execution.

The thresholds used are: $\eta_0 = 1 \times 10^{-5}$, $\eta_1 = 3.5 \times 10^{-10}$, $\delta_0 = 10^{-3}$.

This test allows us to check that the coupling and the transitions between different regimes are executed correctly and that all the numerical schemes are working properly and converging towards the
<figure><img src="image_4.png" /><figcaption>FIGURE 2. Time evolution of density for **Test 1**. In the first row the regime adaptation and the density at three different times are displayed. Red, green and black dots, identify, respectively, cells updated using Euler equations, Boltzmann equations with ES-BGK operator and Boltzmann equations with Boltzmann (Hard-Sphere) operator, using the new hybrid scheme. The cyan line is the analytical solution, the magenta one is the solution obtained using only the Euler equations, and the blue line is computed using only the Boltzmann equation. The same plots, but for the temperature, are showed in the second row. The time evolution of the density for the hybrid scheme in the 2D space (solution constant along the *y*-direction) is displayed in the last row.</figcaption></figure>

right solution. Indeed, since this test is performed in a very stiff regime (i.e. $\epsilon \ll 1$) then the solution of the kinetic levels must converge towards the hydrodynamic limit (Euler equations). It is clear that the numerical solution is strongly consistent with the analytical one. For this specific test case, the hybrid solution is approximately 1.8 times faster than the full kinetic one.

**5.2. Test 2: Fluid flux towards a fixed rectangular obstacle.** One of the aims of this test is to show the capability of the scheme to handle more complex geometries and real world phenomena. In
particular we use our code to simulate the flow around a fixed object with rectangular shape, which is embedded in a computational domain of size $L_x = L_y = 1.0$ discretized with 100 points in each direction. We have discretized the velocity domain using 16 points in the three directions, and we considered a cut-off domain $\Omega_v = [-8, 8]^3$. The value of $R$ has been taken approximately equal to 3.3, and both numbers of points $A_1$ and $A_2$ discretizing the spherical integration of the Boltzmann operator (see Section 2.1) have been taken equal to 4.

To simulate the presence of a rectangular object, we have defined a rectangular domain (white cells in Figure 3 and Figure 4) in which the integrators are never executed. These cells, are part of the stencils of the closest cells of the integration domain, thus we have assigned to them a constant and uniform value of temperature, pressure, density and velocity, given by: $P_* = 1.5, \rho_* = 1.0, T_* = P_*/\rho_*$ and $\vec{u}_* = \vec{0}$. This rectangular shape domain, is surrounded by a small region (2 cells thick) in which the integrator is always the one associated to the full Boltzmann equation with Boltzmann operator. The fluid in the integration domain is initialized with the following parameters: $P = 1.0, \rho = 1.0, \vec{u} = (3.0, 0.0, 0.0)$. The distribution functions are initialized everywhere as a Maxwellian whose moments corresponds to the initial values of the Euler equations.

In Figures 3 and 4, are showed, respectively, the density and temperature profiles, with $\epsilon = 10^{-4}$, at time $t = 0.075, 0.15, 0.225$ for the hybrid scheme, the full Euler integrator and for the Boltzmann operator (computed using CUDA GPU parallelization), as well as the domain indicators for the hybrid scheme (red, green and black dots indicate cells computed using, respectively, Euler equations, ES-BGK operator and Boltzmann operator). In this simulation we considered $\Xi = 5/3$ as specific heat ratio. We estimated, as explained in Section 5, the total duration of the numerical integration performed using just the full Boltzmann operator on CPU over the entire domain, and the ratio between the execution time for the full Boltzmann and the hybrid scheme is equal to 3. The thresholds used are: $\eta_0 = 1 \times 10^{-3}$, $\eta_1 = 2.8 \times 10^{-5}$, $\delta_0 = 10^{-4}$. The time step is equal to $0.1dx$.

It is important to remark that the kinetic regions tend to be confined in time, and in particular the one computed with the Boltzmann operator, which is the most computationally expensive. Indeed the Boltzmann operator is triggered only really close to the object and along the frontal shocks. This feature show the robustness and effectiveness of the implemented criterion.

### 5.3. Test 3: Moving object in a fluid at rest.
The aim of this test is to show the capability of the scheme to handle complex moving geometries, real world phenomena and moving boundaries inside the integration domain. In particular we use our code to simulate the flow generated by a moving object in a fluid at rest, which is embedded in a computational domain of size $L_x = L_y = 1.0$ discretized with 100 points in each direction. We have discretized the velocity domain using 16 points in the three directions, and we considered a cut-off domain $\Omega_v = [-8, 8]^3$. The value of $R$ has been taken approximately equal to 3.3, and both numbers of points $A_1$ and $A_2$ discretizing the spherical integration of the Boltzmann operator (see Section 2.1) have been taken equal to 4.

To simulate the presence of a rectangular object, we have defined a rectangular domain (white cells in Figures 5 and 6) in which the integrators are never executed. These cells, are part of the stencils of the closest cells of the integration domain, thus we have assigned to them a constant and uniform value of pressure, density, temperature and velocity, given by: $P_* = 1.3, \rho_* = 0.3, T_* = P_*/\rho_*$ and $\vec{u}_* = \vec{0}$. This rectangular shape domain, is surrounded by a small region (2 cells thick) in which the integrator is always the one associated to the full Boltzmann equation with Boltzmann operator. The objects moves in the negative direction of the y-axis, by one cell, at each timestep. The fluid in the integration domain is initialized with the following values of pressure, density, temperature and velocity: $P = 0.9, \rho = 0.3, T = P/\rho$ and $\vec{u} = \vec{0}$. The distribution functions are initialized everywhere as a Maxwellian, whose moments correspond to the initial values of the Euler equations.

In Figure 5 we show the density profile with $\epsilon = 10^{-4}$, at time $t = 0.032, 0.054, 0.080$ for the hybrid scheme, for the full Euler equations and for the Boltzmann (Hard-Sphere) operator (computed using
<figure><img src="image_5.png" /><figcaption>FIGURE 3. Time evolution of density for **Test 2**. The first row displays the regime adaptation at three different times. Red, green and black dots, identify, respectively, cells updated using Euler equations, ES-BGK equation and Boltzmann equation. In the second, third and fourth row are displayed, respectively, the time evolution of density computed using hybrid scheme, Euler equations and Boltzmann equation.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>FIGURE 4. Time evolution of temperature for **Test 2**. In the first, second and third row are displayed, respectively, the time evolution of temperature computed using hybrid scheme, Euler equations and Boltzmann equation.</figcaption></figure>

CUDA GPU parallelization), as well as the domain indicators for the hybrid scheme (red, green and black dots indicate cells computed using, respectively, Euler equations, ES-BGK equation and Boltzmann equation). It is interesting to note that some cells that were initially hydrodynamic (red) become kinetic, and then they return hydrodynamic again at the end of the simulations, since the distribution function converges to equilibrium (Maxwellian). In Figure 6 is reported the temperature profile for the same instants in time, for the hybrid scheme, for the Euler and Boltzmann equations. In this simulation we considered $\Xi = 5/3$ as specific heat ratio. We estimated, as explained in Section 5, the total duration of the numerical integration performed using just the full Boltzmann operator on CPU over the entire domain, and the ratio between the execution time for the full Boltzmann and the hybrid scheme (Euler equations, ES-BGK operator and Boltzmann operator) is equal to 2.4.

The thresholds used are: $\eta_0 = 1 \times 10^{-3}$, $\eta_1 = 2.8 \times 10^{-5}$, $\delta_0 = 10^{-4}$. The timestep $t_H$ associated to the hydrodynamic domain of the hybrid scheme has been updated in such a way to guarantee $CFL = 0.95$,
and the timestep is 0.1dx. It is important to remark that the kinetic regions tend to be confined in time, and in particular the one computed with the Boltzmann operator, which is the most computationally expensive. Indeed the Boltzmann operator is triggered only really close to the object and along the frontal shocks. This feature show the robustness and effectiveness of the implemented criterion.

5.4. **Test 4: Fluid flux towards a fixed round obstacle.** One of the aim of this test is to show the capability of the scheme to handle round geometries. In particular we use our code to simulate the flow around a fixed round object, which is embedded in a computational domain of size $L_x = L_y = 1.0$ discretized with 100 points in each direction. We have discretized the velocity domain using 16 points in the three directions, and we considered a cut-off domain $\Omega_v = [-8, 8]^3$. The value of $R$ has been taken approximately equal to 3.3, and both numbers of points $A_1$ and $A_2$ discretizing the spherical integration of the Boltzmann operator (see Section 2.1) have been taken equal to 4.

To simulate the presence of the object, we have defined a circular domain (white cells in Figure 7 and Figure 8) in which the integrators are never executed. These cells, are part of the stencils of the closest cells of the integration domain, thus we have assigned to them a constant and uniform value of pressure, density, temperature and velocity, given by: $P_* = 1.5, \rho_* = 1.0, T_* = P_*/\rho_*$ and $\vec{u}_* = \vec{0}$. This circular shape domain, is surrounded by a small region (2 cells thick) in which the integrator is always the one associated to the full Boltzmann equation with Boltzmann operator. We identify by $(x_c, y_c) = (3N_x/10, N_y/2)$ the coordinates of the centre of the round shape, and we identify all the cells satisfying the following condition $(x_i - x_c)^2 + (y_i - y_c)^2 \le R_c^2$ as part of the object (white cells), where $R_c = N_x/10$ represents the radius of the round obstacle.

The fluid in the integration domain is initialized with the following values of pressure, density, temperature and velocity: $P = 1.0, \rho = 1.0, T = P/\rho$ and $\vec{u}_* = (1.0, 0.0, 0.0)$. The distribution functions are initialized everywhere as a Maxwellian whose moments correspond to the initial values of the Euler equations.

In Figure 7 we show the density profile with $\epsilon = 10^{-4}$, at time $t = 0.050, 0.100, 0.150$ for the hybrid scheme, for the Euler equations, and for the Boltzmann operator (this last one computed using CUDA GPU parallelization), as well as the domain indicators for the hybrid scheme (red, green and black dots indicate cells computed using, respectively, Euler, ES-BGK and Boltzmann equations). It is important to note that some kinetic cells (green or black), after a certain amount of time, become hydrodynamic (red) again as the fluid moves toward equilibrium (Maxwellian distribution). In Figure 8 is reported the temperature profile for the same instants in time, for the hybrid scheme, the Euler and Boltzmann equations. We estimated, as explained in Section 5, the total duration of the numerical integration performed using just the full Boltzmann operator on CPU over the entire domain, and the ratio between the execution time for the full Boltzmann and the hybrid scheme (Euler equations, ES-BGK operator and Boltzmann operator) is equal to 2.6.

The thresholds used are: $\eta_0 = 1 \times 10^{-3}$, $\eta_1 = 2.8 \times 10^{-5}$, $\delta_0 = 10^{-4}$. The time step is 0.1dx.

5.5. **Test 5: Non uniform Knudsen number.** The aim of this test is to check if the scheme is able to deal with a spatially non-uniform Knudsen number. This simulation is inspired by one of those presented in [22]. Here we consider a 2D spatial domain discretized with $N_x = 50$ and $N_y = 16$ points, with a physical space $[-L_x, L_x]$ and $[-L_y, L_y]$ of dimensions $L_x = 1/2$ and $L_y$ chosen such that $dy = dx$. The macroscopic initial conditions are the following:

$$ (5.1) \qquad \begin{cases} \rho(x, y) = 1 + 0.1 \sin(\pi x), \\ u_x(x, y) = 0.125, \\ u_y(x, y) = 0, \\ T(x, y) = 1 + 0.2 \cos(\pi x), \end{cases} $$
<figure><img src="image_7.png" /><figcaption>FIGURE 5. Time evolution of density for **Test 3**. The first row displays the regime adaptation at three different times. Red, green and black dots, identify, respectively, cells updated using Euler equations, ES-BGK equation and Boltzmann equation. In the second, third and fourth row are displayed, respectively, the time evolution of density computed using hybrid scheme, Euler equations and Boltzmann equation.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>FIGURE 6. Time evolution of temperature for **Test 3**. In the first, second and third row are displayed, respectively, the time evolution of temperature computed using hybrid scheme, Euler equations and Boltzmann equation.</figcaption></figure>

and the distribution function is initialized as

$$f(0, x, \cdot) = \frac{1}{2} (M_{\rho(x),u(x),T(x)} + M_{\rho(x),-u(x),T(x)}),$$

so it is everywhere far from equilibrium.

The Knudsen number is a function of the physical space accordingly to the following smooth relation $\epsilon(x) = 10^{-6} + 0.5(\arctan(1+30x) + \arctan(1-30x))$, which means that the more external region of the domain is where we expect the Euler regime to be triggered by the hybrid scheme due to the small Knudsen number. The regime is initially set as ES-BGK everywhere, and then we let it evolves accordingly to our scheme. The velocity domain $[-L, L]^3$, where $L = 8$ is discretized with 16 points in each direction. The value of $R$ has been taken approximately equal to 3.3, and both numbers of points $A_1$ and $A_2$ discretizing the spherical integration of the Boltzmann operator (see Section 2.1) have been taken equal to 4.
<figure><img src="image_9.png" /><figcaption>FIGURE 7. Time evolution of density for **Test 4**. In the first row is displayed the regime adaptation at three different times. Red, green and black dots, identify, respectively, cells updated using Euler equations, ES-BGK equation and Boltzmann equation. In the second, third and fourth row are displayed, respectively, the time evolution of density computed using hybrid scheme, Euler equations and Boltzmann equation.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>FIGURE 8. Time evolution of temperature for **Test 4**. In the first, second and third row are displayed, respectively, the time evolution of temperature computed using hybrid scheme, Euler equations and Boltzmann equation.</figcaption></figure>

In Figure 9 it is reported the density and temperature evolution computed with the hybrid scheme (red, green and black crosses are, respectively, the cells computed using Euler equations, ES-BGK model and Hard-Sphere Boltzmann model), using only Euler equations (blue dashed line) and using only Full Boltzmann Hard-Sphere model (magenta dashed line). It is evident that the Euler solution is far from the reference one (Boltzmann Hard-Sphere model) since the Knudsen number is not everywhere sufficiently small. The hybrid scheme correctly triggers the Euler regime only where the Knudsen number is small (close to the extrema of the physical space), and the region associated to the Boltzmann model (black crosses) is well confined in space in the central region. In addition the solution computed with the hybrid scheme is in perfect agreement with the reference solution (magenta line) which is a clear indicator of the correctness of the implemented methods and of the switching procedure between different regimes.

The thresholds used are: $\eta_0 = 4 \times 10^{-2}$, $\eta_1 = 8 \times 10^{-2}$, $\delta_0 = 10^{-3}$. The time step is $0.1dx$. The ratio between the execution time for the full Boltzmann and the hybrid scheme (Euler equations, ES-BGK
<figure><img src="image_11.png" /><figcaption>FIGURE 9. Time evolution of density for **Test 5**. In the first row the regime adaptation and the density at three different times are displayed. Red, green and black dots, identify, respectively, cells updated using Euler equations, Boltzmann equations with ES-BGK operator and Boltzmann equations with Boltzmann (Hard-Sphere) operator, using the new hybrid scheme. The magenta line is the solution obtained using only the Euler equations, and the blue line is computed using only the Boltzmann equation. The same plots, but for the temperature, are showed in the second row. It is clearly visible that the hybrid solution matches almost perfectly the one obtained with the hard-sphere solver.</figcaption></figure>

operator and Boltzmann operator) is equal to 1.3, which means that the full Boltzmann solver is 30% slower than the hybrid scheme.

**5.6. Discussion on speed-up.** As already observed, the speed-up offered by our numerical domain adaptation compared to the simulations computed using only the full hard-sphere Boltzmann solver is almost equal to 3. This is a quantity which depended on several factors, among which, the principles are, the choice of the thresholds which directly determines the sizes of the domains and the actual numerical implementation of the code.

In Table 5.6 we report, for each simulation, the total percentage of cells used by each domain computed over the total number of cells of the numerical integration (i.e. 100% corresponds to the total number of grid cell in a single timestep multiplied by the total number of timesteps). In the last column is indicated the ratio of time execution of a given solver by the time execution of the Euler solver. These execution times are computed for the **Test 2**, and using each solver over the entire spatial domain for the entire duration of the simulation. Of course, the percentage values depends on the thresholds implemented, and the ratio between the execution time strongly depend on the numerical implementation of the code and the machine performances.
<table><caption>TABLE 1. It is indicated for each simulation the total percentage of cells used by each domain computed over the total number of cells of the numerical integration (i.e. 100% corresponds to the total number of grid cell in a single timestep multiplied by the total number of timesteps). In the last column is indicated the ratio of time execution of a given solver by the time execution of the Euler solver. These execution time are computed for the **Test 2**, and using each solver over the entire spatial domain for the entire duration of the simulation.</caption><thead><tr><th rowspan="2"></th><th colspan="5">Total percentage</th><th rowspan="2">Ratio Execution Time</th></tr><tr><th>TEST 1</th><th>TEST 2</th><th>TEST 3</th><th>TEST 4</th><th>TEST 5</th></tr></thead><tbody><tr><th>Euler</th><td>82.5%</td><td>94.1%</td><td>92.6%</td><td>93.0%</td><td>25.1%</td><td>1</td></tr><tr><th>ES-BGK</th><td>6.7%</td><td>3.4%</td><td>6.6%</td><td>4.8%</td><td>61.8%</td><td>4</td></tr><tr><th>Hard-Sphere</th><td>10.8%</td><td>2.5%</td><td>0.8%</td><td>2.2%</td><td>13.1%</td><td>14</td></tr></tbody></table>

## 6. CONCLUSIONS

In this work, we have proposed a hierarchical domain decomposition strategy for solving the Boltzmann equation, and, more generally, kinetic models in multiscale regimes. The approach relies on two complementary switching criteria. The first, inspired by the works of Levermore, Morokoff, Nadiga [42], Filbet and Rey [22], Filbet and Xiong [23], Tiwari [55], Li, Son and Wang [44], and Xiong and Qiu [61], uses Chapman-Enskog expansions to determine when to move from a fluid description to either the ES-BGK approximation or the full Boltzmann model. This criterion relies solely on macroscopic quantities derived from the closure of the kinetic model and does not require evaluating the distribution function itself. The second governs the reverse transition from kinetic to hydrodynamic regimes, by comparing truncated Chapman-Enskog expansions with the equilibrium state.

The key novelty lies in the implementation of a 2D-in-space and 3D-in-velocity hybrid solver that integrates three distinct physical regimes within a single framework. By combining high-order asymptotic-preserving (AP) schemes with fast spectral solvers for the Boltzmann operator, the method ensures stability and accuracy across a wide range of Knudsen numbers. Numerical experiments confirm that kinetic regions are automatically and dynamically confined in both space and time, allowing the solver to retain the accuracy of a full Boltzmann simulation while achieving nearly a threefold speedup. Further gains can be expected through parallelization and memory management.

Beyond the immediate performance improvements, the proposed methodology demonstrates that deterministic solvers for the full Boltzmann equation can be made practical for multiscale simulations in nontrivial geometries. This opens new perspectives for applying kinetic theory to complex physical and engineering systems. Future work will address multi-species mixtures and collisional plasma models, as well as the integration of advanced software engineering high-performance computing strategies to further extend the scope and efficiency of the approach.

## APPENDIX A. MOMENT REALIZABILITY MATRIX

In this Appendix we give more details about the Moment Realizability Matrix and in particular on the derivation of the Equation (3.14). The matrix $\tilde{A}^\epsilon$ and the vector $\tilde{B}^\epsilon$ will allow us to define our hydrodynamic break down criterion. Let us set the vector of the reduced collisional invariants for $V = (v - u^\epsilon)/\sqrt{T^\epsilon}$,

$$m := \left( 1, V, \left( \frac{2}{5} \right)^{1/2} \left( \frac{|V|^2}{2} - \frac{5}{2} \right) \right).$$
Similarly to [42,55], now we define the so-called moment realizability matrix $M$ by setting

$$ (A.1) \qquad M := \frac{1}{\rho^{\varepsilon}} \int_{\mathbb{R}^3} (m \otimes m) f^{\varepsilon}(v) \, dv. $$

Indeed, if the distribution function is nonnegative, then the quantity $\int \phi^2(v) f \, dv$ is also nonnegative [55], for any function $\phi(v)$. The same is true if we consider $m = m(v)$ a columns vector of $c$ polynomials, and we set $\phi = a^T m$ for a generic $a \in \mathbb{R}^c$, and we have

$$ (A.2) \qquad a^T \left( \int_{\mathbb{R}^3} m \otimes m f \, dv \right) a = \int_{\mathbb{R}^3} (a^T m)^2 f \, dv \geq 0. $$

Then if $f$ is nonnegative and not identically zero, then for every $m$ the $c \times c$ matrix

$$ (A.3) \qquad M = \frac{1}{\rho} \int_{\mathbb{R}^3} m \otimes m f \, dv, $$

is positive semi-definite.

Then, the quantity given by Equation A.1 is also positive semi-definite. After some computations we obtain

$$ M = \begin{pmatrix} 1 & \mathbf{0}_{\mathbb{R}^3}^\top & -\left(\frac{2}{5}\right)^{\frac{1}{2}} \\ \mathbf{0}_{\mathbb{R}^3} & I + \bar{A}^\varepsilon & \left(\frac{2}{5}\right)^{1/2} \bar{B}^\varepsilon \\ -\left(\frac{2}{5}\right)^{\frac{1}{2}} & \left(\frac{2}{5}\right)^{1/2} (\bar{B}^\varepsilon)^\top & \bar{C}^\varepsilon \end{pmatrix}, $$

where $\bar{C}^\varepsilon$ is the dimensionless fourth order moment of $f^\varepsilon$

$$ \bar{C}^{\varepsilon} := \frac{2}{5\rho} \int_{\mathbb{R}^3} \left[ \frac{|V|^2}{2} - \frac{5}{2} \right]^2 f^{\varepsilon}(v) \, dv. $$

The matrix $M$ is of the form

$$ (A.4) \qquad M = \begin{pmatrix} M_A & M_B \\ M_B^\top & M_C \end{pmatrix}, $$

where

$$ (A.5) \qquad M_A = \begin{pmatrix} 1 & \mathbf{0}_{\mathbb{R}^3}^\top \\ \mathbf{0}_{\mathbb{R}^3}^\top & I + \bar{A}^\varepsilon \end{pmatrix}, \quad M_B = \begin{pmatrix} -\left(\frac{2}{5}\right)^{\frac{1}{2}} \\ \left(\frac{2}{5}\right)^{1/2} \bar{B}^\varepsilon \end{pmatrix}, \quad M_C = \bar{C}^\varepsilon. $$

Assuming that the scalar quantity $M_C = \bar{C}^\varepsilon > 0$ (we will verify later that this assumption is verified in our circumstances) then the matrix $M$ is positive semi-definite if and only if its Schur complement $X$ is positive semi-definite

$$ (A.6) \qquad X = M_A - M_B M_C^{-1} M_B^\top = \begin{pmatrix} 1 - \frac{2}{5\bar{C}^\varepsilon} & \frac{2}{5} \frac{(\bar{B}^\varepsilon)^\top}{\bar{C}^\varepsilon} \\ \frac{2}{5} \frac{\bar{B}^\varepsilon}{\bar{C}^\varepsilon} & I + \bar{A}^\varepsilon - \frac{2}{5} \frac{\bar{B}^\varepsilon (\bar{B}^\varepsilon)^\top}{\bar{C}^\varepsilon} \end{pmatrix}. $$

For the sake of simplicity, let's introduce the following matrix $Q$ for the change of basis

$$ (A.7) \qquad Q = \begin{pmatrix} \alpha & -\frac{\chi_1^\top}{\chi_1^2} \\ \mathbf{0}_{\mathbb{R}^3} & I \end{pmatrix}, $$

where $\chi_1 = \frac{2}{5} \frac{\bar{B}^\varepsilon}{\bar{C}^\varepsilon}$, $\chi_2 = 1 - \frac{2}{5\bar{C}^\varepsilon}$ and $\alpha^2 = \frac{\bar{C}^\varepsilon}{1 - \frac{2}{5\bar{C}^\varepsilon}}$. Then the operation $Q^\top X Q := W$ corresponds to a change of coordinate (since $X$ can be viewed as quadratic form) and this operation does not change the
semipositiveness of the matrix $\mathbf{X}$, i.e. the signs of the eigenvalues of the matrices $\mathbf{X}$ and $\mathbf{W}$ are the same. So we obtain that

$$ (A.8) \qquad \mathbf{W} := \mathbf{Q}^\top \mathbf{X} \mathbf{Q} = \begin{pmatrix} \bar{C}^\epsilon & \mathbf{0}_{\mathbb{R}^3}^\top \\ \mathbf{0}_{\mathbb{R}^3} & \mathbf{I} + \bar{\mathbf{A}}^\epsilon - \frac{2}{5} \frac{\bar{\mathbf{B}}^\epsilon (\bar{\mathbf{B}}^\epsilon)^\top}{\bar{C}^\epsilon} - \frac{\chi_1 \chi_1^\top}{\chi_2} \end{pmatrix}. $$

It is important to note that doing that we have assumed the validity of the following relation $\alpha^2 = \frac{\bar{C}^\epsilon}{1 - \frac{2}{5\bar{C}^\epsilon}}$, which requires that $\bar{C}^\epsilon > \frac{2}{5}$. Now we want to prove that, in our framework, $\bar{C}^\epsilon = 1 > \frac{2}{5}$. If the distribution function is at equilibrium, i.e. $f = \mathcal{M}_{\rho,u,T}$ then a direct computation gives that $\bar{C}^\epsilon = 1$. More in general, if we assume that the distribution function $f$ has a small deviation from the local Maxwellian $\mathcal{M}_{\rho,u,T}$ then we can approximate it by [55]

$$ (A.9) \qquad f = \mathcal{M}_{\rho,u,T}(1 + \phi), $$

where

$$ (A.10) \qquad \phi = \frac{\langle \mathbf{v} - \mathbf{u}, \mathbf{q} \rangle}{\rho T^2} \left[ \frac{|\mathbf{v} - \mathbf{u}|^2}{5T} - 1 \right] + \frac{1}{2\rho T^2} (\mathbf{v} - \mathbf{u})^\top \mathbf{\tau}(\mathbf{v} - \mathbf{u}), $$

where $\mathbf{\tau}$ and $\mathbf{q}$ are, respectively, the stress tensor and the heat flux

$$ (A.11) \qquad \tau_{ij} = \int_{\mathbb{R}^3} (v_i - u_i)(v_j - u_j) f(t, x, v) dv - p \delta_{ij}, $$

$$ (A.12) \qquad q_j = \frac{1}{2} \int_{\mathbb{R}^3} (v_j - u_j) |v - u|^2 f(t, x, v) dv. $$

Using this last expression for $\phi$, it can be proved [55] that if $f = \mathcal{M}_{\rho,u,T}(1 + \phi)$ then $\bar{C}^\epsilon = 1$. This confirms that in our framework $\bar{C}^\epsilon = 1 > \frac{2}{5}$.

Since $\bar{C}^\epsilon = 1$, then from Equation (A.8) the matrix $\mathbf{W}$ is positive semi-definite, if and only if, the following quantity is semipositive definite

$$ (A.13) \qquad \mathcal{V}_\epsilon := \mathbf{I} + \bar{\mathbf{A}}^\epsilon - \frac{2}{5} \frac{\bar{\mathbf{B}}^\epsilon (\bar{\mathbf{B}}^\epsilon)^\top}{\bar{C}^\epsilon} - \frac{\chi_1 \chi_1^\top}{\chi_2} = \mathbf{I} + \bar{\mathbf{A}}^\epsilon - \frac{\frac{2}{5\bar{C}^\epsilon}}{1 - \frac{2}{5\bar{C}^\epsilon}} \bar{\mathbf{B}}^\epsilon \otimes \bar{\mathbf{B}}^\epsilon. $$

Substituting $\bar{C}^\epsilon = 1$ in this last equation we get

$$ (A.14) \qquad \mathcal{V}_\epsilon = \mathbf{I} + \bar{\mathbf{A}}^\epsilon - \frac{2}{3} \bar{\mathbf{B}}^\epsilon \otimes \bar{\mathbf{B}}^\epsilon. $$

From this derivation, the initial matrix $\mathbf{M}$ defined in Equation (A.3), is positive semi-definite if and only if the matrix $\mathcal{V}_\epsilon$ is positive semi-definite. Then the idea is to monitor the evolution of the matrix $\mathcal{V}_\epsilon$ [42,54]: indeed the moment realizability criterion states that the fluid dynamic description has broken down when the perturbation is too large that $\mathcal{V}_\epsilon$, and thus $\mathbf{M}$, is no longer positive semi-definite.

## ACKNOWLEDGMENTS

The authors received funding from the European Union's Horizon Europe research and innovation program under the Marie Skłodowska-Curie Doctoral Network DataHyking (Grant No. 101072546). DC and TR are also supported by the French government, through the UniCA<sub>JEDI</sub> Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR-15-IDEX-01. DC would like to thank Tommaso Tenna for the useful discussions about the manuscript. The research of LP has been supported by the Royal Society under the Wolfson Fellowship "Uncertainty quantification, data-driven simulations and learning of multiscale complex systems governed by PDEs". This work has been written within the activities of GNCS group of INdAM (Italian National Institute of High Mathematics). LP also acknowledges the partial support by European Union - NextGenerationEU
through the Italian Ministry of University and Research as part of the PNRR – Mission 4 Component 2, Investment 1.3 (MUR Directorial Decree no. 341 of 03/15/2022), FAIR “Future” Partnership Artificial Intelligence Research”, Proposal Code PE00000013 - CUP DJ33C22002830006) and by MIUR-PRIN Project 2022, No. 2022KKJP4X “Advanced numerical methods for time dependent parametric partial differential equations with applications”.

## REFERENCES

[1] ALAIA, A., AND PUPPO, G. A hybrid method for hydrodynamic-kinetic flow Part I: A particle-grid method for reducing stochastic noise in kinetic regimes. *J. Comput. Phys.* **230**, 14 (2011), 5660–5683.

[2] ALAIA, A., AND PUPPO, G. A hybrid method for hydrodynamic-kinetic flow – Part II – Coupling of hydrodynamic and kinetic models. *J. Comput. Phys.* **231**, 16 (2012), 5217–5242.

[3] ANDRIES, P., LE TALLEC, P., PERLAT, J.-P., AND PERTHAME, B. The Gaussian-BGK model of Boltzmann equation with small Prandtl number. *Eur. J. Mech. B Fluids* **19**, 6 (Nov. 2000), 813–830.

[4] 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. *Phys. Rev.* **94** (May 1954), 511–525.

[5] BOLTZMANN, L. Weitere studien über das wärmegleichgewicht unter gasmoleklen. *Wiener Berichte* **66** (1872), 275–370.

[6] BOYD, I. D., CHEN, G., AND CANDLER, G. V. Predicting failure of the continuum fluid equations in transitional hypersonic flows. *Phys. Fluids* **7**, 1 (01 1995), 210–219.

[7] BURT, J. M., AND BOYD, I. D. A hybrid particle approach for continuum and rarefied flow simulation. *J. Comput. Phys.* **228**, 2 (2009), 460–475.

[8] CAFLISCH, R., JIN, S., AND RUSSO, G. Uniformly Accurate Schemes For Hyperbolic Systems With Relaxation. *SIAM J. Numer. Anal.* **34** (07 1995).

[9] CAFLISCH, R., WANG, C., DIMARCO, G., COHEN, B., AND DIMITS, A. A hybrid method for accelerated simulation of Coulomb collisions in a plasma. *Multiscale Model. Simul.* **7**, 2 (2008), 865–887.

[10] CANUTO, C., HUSSAINI, M., QUARTERONI, A., AND ZANG, T. *Spectral Methods in Fluid Dynamics*. Springer Series in Computational Physics. Springer-Verlag, New York, 1988.

[11] CORON, F., AND PERTHAME, B. Numerical passage from kinetic to fluid equations. *SIAM J. Numer. Anal.* **28** (1991), 26–42.

[12] DEGOND, P., AND DIMARCO, G. Fluid simulations with localized boltzmann upscaling by direct simulation Monte-Carlo. *J. Comput. Phys.* **231**, 6 (2012), 2414–2437.

[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. *J. Comput. Phys.* **274** (2014), 122–139.

[14] DIMARCO, G., AND PARESCHI, L. Hybrid multiscale methods. II. Kinetic equations. *Multiscale Model. Simul.* **6**, 4 (2007/08), 1169–1197.

[15] DIMARCO, G., AND PARESCHI, L. Domain decomposition techniques and hybrid multiscale methods for kinetic equations. In *Hyperbolic Problems: Theory, Numerics, Applications*. Springer, Berlin, 2008, pp. 457–464.

[16] DIMARCO, G., AND PARESCHI, L. Fluid solver independent hybrid methods for multiscale kinetic equations. *SIAM J. Sci. Comput.* **32**, 2 (2010), 603–634.

[17] DIMARCO, G., AND PARESCHI, L. Asymptotic preserving implicit-explicit Runge-Kutta methods for nonlinear kinetic equations. *SIAM J. Numer. Anal.* **51**, 2 (2013), 1064–1087.

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

[19] FILBET, F., AND JIN, S. A class of asymptotic-preserving schemes for kinetic equations and related problems with stiff sources. *J. Comput. Phys.* **229**, 20 (2010), 7625–7648.

[20] FILBET, F., AND JIN, S. An Asymptotic Preserving Scheme for the ES-BGK Model of the Boltzmann Equation. *J. Sci. Comput.* **46** (02 2011), 204–224.

[21] FILBET, F., MOUHOT, C., AND PARESCHI, L. Solving the Boltzmann Equation in N log2 N. *SIAM J. Sci. Comput.* **28**, 3 (2007), 1029–1053.

[22] FILBET, F., AND REY, T. A Hierarchy of Hybrid Numerical Method for Multi-Scale Kinetic Equations. *SIAM J. Sci. Comput.* **37**, 3 (2015), A1218–A1247.

[23] FILBET, F., AND XIONG, T. A hybrid discontinuous Galerkin scheme for multi-scale kinetic equations. *J. Comput. Phys.* **372** (2018), 841–863.

[24] FRIGO, M., AND JOHNSON, S. G. The Design and Implementation of FFTW3. *Proc. IEEE* **93**, 2 (2005), 216–231. Special issue on “Program Generation, Optimization, and Platform Adaptation”.

[25] GARG, G., MANKODI, T., ESMAEILIFAR, E., AND MYONG, R. S. Neural network-based finite volume method and direct simulation Monte Carlo solutions of non-equilibrium shock flow guided by nonlinear coupled constitutive relations. *Phys. Fluids* **36** (10 2024), 106113.
[26] GOLSE, F. The Boltzmann equation and its hydrodynamic limits. In *Handbook of Differential Equations: Evolutionary Equations Vol. 2*, C. Dafermos and E. Feireisl, Eds. North-Holland, 2005, pp. 159–303.

[27] GOSSE, L., AND TOSCANI, G. An asymptotic-preserving well-balanced scheme for the hyperbolic heat equations. *Cr. Math. 334*, 4 (2002), 337–342.

[28] GOTTLIEB, S., AND SHU, C.-W. Total Variation Diminishing Runge-Kutta Schemes. *Math. Comput.* **67** (08 1996).

[29] GUO, Z., WANG, R., AND XU, K. Discrete unified gas kinetic scheme for all Knudsen number flows: II. Compressible case. *Phys. Rev. E* **91** (06 2014).

[30] GUO, Z., AND XU, K. Progress of discrete unified gas-kinetic scheme for multiscale flows. *Advances in Aerodynamics* **3** (12 2021).

[31] GUO, Z., XU, K., AND WANG, R. Discrete unified gas kinetic scheme for all Knudsen number flows: Low-speed isothermal case. *Phys. Rev. E* **88** (09 2013), 033305.

[32] JIN, S. Runge-Kutta Methods for Hyperbolic Conservation Laws with Stiff Relaxation Terms. *J. Comput. Phys.* **122**, 1 (1995), 51–67.

[33] JIN, S. Efficient Asymptotic-Preserving (AP) Schemes For Some Multiscale Kinetic Equations. *SIAM J. Sci. Comput.* **21**, 2 (Jan. 1999), 441–454.

[34] JIN, S., AND PARESCHI, L. Asymptotic-Preserving (Ap) Schemes for Multiscale Kinetic Equations: a Unified Approach. In *Hyperbolic Problems: Theory, Numerics, Applications*, H. Freistühler and G. Warnecke, Eds. Birkhäuser Basel, Basel, 2001, pp. 573–582.

[35] JIN, S., PARESCHI, L., AND TOSCANI, G. Diffusive Relaxation Schemes For Multiscale Discrete-Velocity Kinetic Equations. *SIAM J. Numer. Anal.* **35** (12 1998), 2405–2439.

[36] KLAR, A. An Asymptotic-Induced Scheme for Nonstationary Transport Equations in the Diffusive Limit. *SIAM J. Numer. Anal.* **35**, 3 (1998), 1073–1094.

[37] KOLOBOV, V., ARSLANBEKOV, R., ARISTOV, V., FROLOVA, A., AND ZABELOK, S. Unified solver for rarefied and continuum flows with adaptive mesh and algorithm refinement. *J. Comput. Phys.* **223**, 2 (2007), 589–608.

[38] LAIDIN, T. Hybrid kinetic/fluid numerical method for the Vlasov-BGK equation in the diffusive scaling. *Kinet. Relat. Mod.* **16** (2023).

[39] LAIDIN, T., AND REY, T. Hybrid Kinetic/Fluid Numerical Method for the Vlasov-Poisson-BGK Equation in the Diffusive Scaling. In *Finite Volumes for Complex Applications X—Volume 2, Hyperbolic and Related Problems* (Cham, 2023), E. Franck, J. Fuhrmann, V. Michel-Dansac, and L. Navoret, Eds., Springer Nature Switzerland, pp. 229–237.

[40] LEE, J. Y., SCHOTTHÖFER, S., XIAO, T., KRUMSCHEID, S., AND FRANK, M. Structure-Preserving Operator Learning: Modeling the Collision Operator of Kinetic Equations. arXiv preprint arXiv:2402.16613, 2024.

[41] LEVERMORE, C. Moment closure hierarchies for kinetic theories. *J. Stat. Phys.* **83** (06 1996), 1021–1065.

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

[43] LEVY, D., PUPPO, G., AND RUSSO, G. Compact Central WENO Schemes for Multidimensional Conservation Laws. *SIAM J. Sci. Comput.* **22** (09 2000).

[44] LI, W., SONG, P., AND WANG, Y. A Hybrid Moment Method for Multi-scale Kinetic Equations Based on Maximum Entropy Principle. *J. Sci. Comput.* **89**, 2 (Nov. 2021).

[45] MAXWELL, J. C. On the Stability of the Motion of Saturn’s Rings; an Essay which obtained the Adams’ Prize for the Year 1856, in the University of Cambridge. *Mon. Not. R. Astron Soc.* **19**, 8 (1859).

[46] MOUHOT, C., AND PARESCHI, L. Fast algorithms for computing the Boltzmann collision operator. *Math. Comput.* **75**, 256 (2006), 1833–1852 (electronic).

[47] PARESCHI, L., AND RUSSO, G. Numerical Solution of the Boltzmann Equation I : Spectrally Accurate Approximation of the Collision Operator. *SIAM J. Numer. Anal.* **37**, 4 (2000), 1217–1245.

[48] PARESCHI, L., AND RUSSO, G. *Implicit-Explicit Runge-Kutta schemes for stiff systems of differential equations*, vol. 3. Nova Sci. Publ., 01 2001, pp. 269–288.

[49] PUPPO, G., SEMPLICE, M., AND VISCONTI, G. Quinpi: Integrating Conservation Laws with CWENO Implicit Methods. *Commun. Appl. Math. Comput.* **5** (02 2022).

[50] REY, T. *Contributions to the Mathematical and Numerical Analysis of Multiscale Kinetic Equations*. Habilitation à diriger des recherches, Université de Lille, Feb. 2023. 128 pages.

[51] SAINT-RAYMOND, L. *Hydrodynamic Limits of the Boltzmann Equation*. Springer, 2009.

[52] SHU, C.-W., AND OSHER, S. Efficient implementation of essentially non-oscillatory shock-capturing schemes, II. *J. Comput. Phys.* **83**, 1 (1989), 32–78.

[53] STRUCHTRUP, H. *Macroscopic Transport Equations for Rarefied Gas Flows*. Springer Berlin, Heidelberg, 2005.

[54] TIWARI, S. Coupling of the Boltzmann and Euler equations with automatic domain decomposition. *J. Comput. Phys.* **144**, 2 (Aug. 1998), 710–726.
[55] TIWARI, S. Application of Moment Realizability Criteria for the Coupling of the Boltzmann and Euler Equations. Transport Theor. Stat. 29, 7 (Nov. 2000), 759–783.

[56] TIWARI, S., AND KLAR, A. An adaptive domain decomposition procedure for Boltzmann and Euler equations. J. Comput. Appl. Math. 90, 2 (1998), 223–237.

[57] TIWARI, S., KLAR, A., AND HARDT, S. A particle–particle hybrid method for kinetic and continuum equations. J. Comput. Phys. 228, 18 (2009), 7109–7124.

[58] TORO, E. *Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction*. Springer Science & Business Media, 01 2009.

[59] WU, C., SHI, B., CHAI, Z., AND WANG, P. Discrete unified gas kinetic scheme with a force term for incompressible fluid flows. *Computers & Mathematics with Applications* 71, 12 (2016), 2608–2629.

[60] XIAO, T., AND FRANK, M. Using neural networks to accelerate the solution of the Boltzmann equation. *J. Comput. Phys.* 443 (10 2021), 110521.

[61] XIONG, T., AND QIU, J.-M. A hierarchical uniformly high order DG-IMEX scheme for the 1D BGK equation. *J. Comput. Phys.* 336 (2017), 164–191.

DOMENICO CAPARELLO, UNIVERSITÉ CÔTE D'AZUR, CNRS, LJAD, PARC VALROSE, F-06108 NICE, FRANCE, DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE UNIVERSITY OF FERRARA

Email address: Domenico.CAPARELLO@univ-cotedazur.fr

Email address: domenico.caparello@unife.it

LORENZO PARESCHI, MAXWELL INSTITUTE FOR MATHEMATICAL SCIENCES AND DEPARTMENT OF MATHEMATICS HERIOT-WATT UNIVERSITY, EDINBURGH, DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE UNIVERSITY OF FERRARA,

Email address: L.Pareschi@hw.ac.uk

THOMAS REY, UNIVERSITÉ CÔTE D'AZUR, CNRS, LJAD, PARC VALROSE, F-06108 NICE, FRANCE

Email address: Thomas.REY@univ-cotedazur.fr