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

HAL
open science

# High-Order Asymptotic-Preserving IMEX schemes for an ES-BGK model for Gas Mixtures

Domenico Caparello, Lorenzo Pareschi, Thomas Rey

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

## To cite this version:

Domenico Caparello, Lorenzo Pareschi, Thomas Rey. High-Order Asymptotic-Preserving IMEX schemes for an ES-BGK model for Gas Mixtures. 2025. hal-05385222

HAL Id: hal-05385222

https://hal.science/hal-05385222v1

Preprint submitted on 27 Nov 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_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
High-Order Asymptotic-Preserving IMEX schemes for an
ES-BGK model for Gas Mixtures

Domenico Caparello*¹,², Lorenzo Pareschi†²,³, and Thomas Rey‡¹

¹Université Côte d'Azur, CNRS, LJAD, Parc Valrose, F-06108 Nice, France

²Department of Mathematics and Computer Science University of Ferrara

³Heriot-Watt University, Edinburgh, UK

Abstract

In this work we construct a high-order Asymptotic-Preserving (AP) Implicit-Explicit (IMEX) scheme for the ES-BGK model for gas mixtures introduced in [16]. The time discretization is based on the IMEX strategy proposed in [26] for the single-species BGK model and is here extended to the multi-species ES-BGK setting. The resulting method is fully explicit, uniformly stable with respect to the Knudsen number and, in the fluid regime, it reduces to a consistent and high-order accurate solver for the limiting macroscopic equations of the mixture. The IMEX structure removes the stiffness associated with the relaxation term so that the time step is constrained only by a hyperbolic CFL condition. The full solver couples a high-order space and velocity discretization that includes third-order time integration, a CWENO3 finite-volume reconstruction in space, exact conservation of macroscopic moments in the discrete velocity space, and a multithreaded implementation. The proposed approach can handle an arbitrary number of species. Its accuracy and robustness are demonstrated on a set of multidimensional kinetic tests for gas mixtures, where the AP property and the correct asymptotics are numerically verified across different regimes.

# 1 Introduction

Mixtures of rarefied gas can be found in a wide variety of systems, such as high-altitude gas, nuclear engineering applications, evaporation-condensation processes [37], or in the context of plasma physics, by considering, for example, the evolution of electrons and ions as different species [13]. Here we describe some examples (coming from [37]) in which macroscopic models (like Euler or Navier-Stokes equations) fail to correctly describe the fluid, and thus, kinetic models of the Boltzmann types are needed to describe these systems. In the space shuttle reentry problem, at high altitudes, free streaming effects are significant and thus macroscopic models are not adequate in describing the system evolution. The shock profiles computed by using the Navier-Stokes equations are not accurate for hypersonic flows, while the Boltzmann equation is able to reproduce the correct behaviour. Thus, given the complexity of gas mixtures, a realistic mathematical depiction of such systems must involve multi-species kinetic models.

The analytical formulation of the Boltzmann equation for gas mixtures has been examined in several works, like [14, 15]. In [47] the authors extend the gas mixtures Boltzmann equation framework to granular gases and study its well posedness in the homogeneous case as well as its large time behavior.

Such kind of microscopic kinetic models (like the Boltzmann equation) often contains small parameters, like for example the mean free path of particles, and an asymptotic expansion of the kinetic models on these parameters typically yield macroscopic equations (like the Euler equations), which describe the average behaviour of the fluid. The asymptotic resolution of these scales using kinetic models gives rise to significant numerical challenges due to the high computational cost required, thus in case of validity of the macroscopic equations, it is much less expensive to use them instead of integrating

<footnote>*Domenico.caparello@univ-cotedazur.fr, domenico.caparello@unife.it</footnote>

<footnote>†L.Pareschi@hw.ac.uk</footnote>

<footnote>‡Thomas.rey@univ-cotedazur.fr</footnote>
the microscopic kinetic model. However, we have already seen that in some cases such macroscopic equations can break down in part of the domain, and thus kinetic microscopic models are needed in such cases. A possible approach for dealing with these scenarios is by using multi-physics dynamic domain decomposition strategies, where the scheme employs macroscopic or kinetic equations according to the non-equilibrium behaviour of the fluid, see for example [20, 28, 49]. Another approach for such multiscale problems, which we also adopt in our scheme, is to employ an asymptotic-preserving (AP) strategy. The core idea behind these schemes is to ensure that the macroscopic models are correctly recovered at the discrete level in the asymptotic limit of the underlying microscopic model. Furthermore, for a scheme to be classified as asymptotic-preserving it must fulfil a second condition, namely the implicit collision term must be implemented efficiently, typically allowing the time step to be much larger than the small parameter (typically the Knudsen number, i.e. the ratio between the mean free path of particles and the physical length scale of the system) of the problem. This brief introduction to AP schemes that we have just presented is based on [35]. Several works focused on Asymptotic-Preserving schemes for kinetic equation can be found in literature, as [12, 23, 25, 26, 33, 37].

Gas mixtures models describe the evolution of the system by evolving the distribution of particles $f = f(t, x, v)$ of each species, which is a non-negative function representing the probability to find a particle in a volume $dx$ around $x$, and in an interval $dv$ around $v$, during a small time $dt$. The distribution function is normalized such that $\iint f(t, x, v) dx dv = N$ where $N$ is the total number of particles. Mixture of $P$ species are each described by a particle distribution function $f_p(t, x, v)$ for $p \in \{1, \dots, P\}$ and normalized such that $\iint f_p(t, x, v) dx dv = N_p$, where $N_p$ is the total number of particles of the $p$-th species.

Given $t \ge 0$, $x \in \Omega \subset \mathbb{R}^{Dx}$ and $v \in \mathbb{R}^{Dv}$, each function $f_p$ evolves according the multi-species Boltzmann equation

$$ \left\{ \begin{aligned} \partial_t f_p + v \cdot \nabla_x f_p &= \frac{1}{\varepsilon} Q_p[\mathbf{f}], \\ Q_p[\mathbf{f}] &= \sum_{q=1}^{P} Q_{p,q}[f_p, f_q], \\ f_p(0, x, v) &= f_{p,0}(x, v), \end{aligned} \right. \qquad (1) $$

where $Q_p[\mathbf{f}]$ is the collision operator which describes the changes in momentum due to collisions between particles and $\varepsilon$ is the Knudsen number, which is the ratio between the mean free path of particles and the physical length scale of the system: thus, if $\varepsilon \ll 1$ then the fluid can be described using the so-called hydrodynamic limit (for example, using the compressible Euler equations, Navier-Stokes or Burnett, etc.).

The collision operator is given by [2]

$$ Q_{p,q}[f_p, f_q] = \iint_{\mathbb{R}^{D_V} \times \mathbb{S}^{D_V}-1} B_{p,q}(|v - v_*|, \cos \theta) [f_p(v')f_q(v') - f_p(v)f_q(v_*)] dv_* d\sigma, \quad (2) $$

where the inter-species collisions can be parametrized as

$$ \left\{ \begin{aligned} v' &= \frac{1}{m_p + m_q} (m_p v + m_q v_* + m_q |v - v_*| \sigma), \\ v'_* &= \frac{1}{m_p + m_q} (m_p v + m_q v_* - m_p |v - v_*| \sigma), \end{aligned} \right. \qquad (3) $$

where $m_p$ is the mass of each particle of the $p$-th species.

The quantity $B_{p,q}$ is the inter-species collision kernel: a non-negative function describing the physics of the collisions between particles. More details on the collision kernel for the single species case (and for the Boltzmann equation in general) can be found in [24].

The collision operator must preserve the mass of each species, the total momentum and the total
kinetic energy of the system. In formulas, the following quantities must be conserved

$$ \left\{
\begin{array}{ll}
\displaystyle \int f_p \, dv, & \forall p \in \{1, 2, \dots, P\}, \\[2ex]
\displaystyle \sum_{p=1}^{P} \int \int m_p v f_p(t, x, v) \, dx \, dv, & \\[2ex]
\displaystyle \sum_{p=1}^{P} \int \int \frac{m_p}{2} |v|^2 f_p(t, x, v) \, dx \, dv.
\end{array}
\right.
\qquad (4) $$

Every species has its own global equilibrium distribution, that are all Maxwellians with common average velocity $u$ and temperature $T$

$$ \mathcal{M}_p^{n_p, u, T}(t, x, v) = n_p(t, x) \left(\frac{m_p}{2\pi T}\right)^{\frac{D_V}{2}} \exp\left(-\frac{m_p(u-v)^2}{2T}\right), \qquad (5) $$

where the numerical density $n_p$, the average velocity $v_p$ and the temperature $T_p$ of the p-th species are given by

$$ n_p = \int f_p \, dv, \qquad \rho_p v_p = \int m_p v f_p \, dv, \qquad \frac{D_V}{2} n_p T_p = \int \frac{m_p}{2} (v_p - v)^2 f_p \, dv, \qquad (6) $$

with $\rho_p = m_p n_p$ and $\mathcal{P}_p = n_p T_p$; these last two quantities are, respectively, the density and the pressure of the p-th species, and $D_V$ is the dimensions of the velocity space. The analogous total moments for the mixture are given by $n = \sum n_p$, $\rho = \sum \rho_p$, $\rho u = \sum \rho_p v_p$, $\mathcal{P} = \sum \mathcal{P}_p$ and

$$ \frac{D_V}{2} n T = \sum_{p=1}^{P} \int \frac{m_p}{2} (u - v)^2 f_p \, dv. \qquad (7) $$

It is possible to define also a macroscopic total energy of the mixture as follows

$$ E = \sum_p \int \frac{m_p}{2} |v|^2 f_p dv = \frac{1}{2} \rho |u|^2 + \frac{D_V}{2} n T. \qquad (8) $$

In literature several possible models for the multi-species Boltzmann equation have been proposed. A possible approach consists in a generalization of the single species BGK model [4]: the main idea is that the collision operator drives the distribution function towards the equilibrium Maxwellian through a simple relaxation mechanism. As discussed in [30], the extension of this BGK strategy to gas mixtures is not trivial, since it can give rise to breakdown of positivity for density or temperature fields, and/or of the indifferentiability principle. Another possible approach is to replace the collision operator with a sum of relaxation operators (like [8, 31]) of the BGK kind, each describing the interactions with any other species (itself or different).

In [6] the authors propose a binary mixture model where the collisions between particles of the same species are treated using the full Boltzmann operator, while the BGK relaxation terms are used to describe the interactions between particles of different species. They also validate their model with some interesting numerical simulations. In the recent paper [53], the authors propose an implicit discrete unified gas kinetic scheme (DUGKS) for multiscale steady flows of binary gas mixtures.

Of particular interest is [10], where the authors present a comparison of different models for gas mixtures proposed by [1, 5, 8]. In [37] the authors extends to the multi-species Boltzmann equation the BGK penalization technique of the full Boltzmann operator introduced by [25]: this technique allows to rewrite the Boltzmann equation as a sum of a stiff and a non-stiff term.

In [30] the authors propose an ellipsoidal BGK model for a binary mixture in the kinetic theory framework. Different spectral methods for the full multi-species Boltzmann operator have been proposed, for example, by [34, 51].

In this context, the derivation of the corresponding hydrodynamic equations from kinetic models [39] is of particular interest, for example for inert gas mixtures, reactive mixtures or multi-phase flows [5, 7, 46].
A big challenge in numerical schemes for the multi-species Boltzmann equation (as well as in plasma contexts) arises from the non-unity mass ratios among species. This difficulty becomes more problematic as the mass ratio increases, since larger velocity domains are required, which in turn necessitates more discretization cells and consequently leads to higher computational costs [2]. Our simulations will show that our numerical scheme is able to correctly integrate unequal-mass configurations. In particular we perform some tests with a mass ratio equal to 100 in the spatially homogeneous settings, and up to 20 in the spatially non-homogeneous setting.

In this work we are interested in developing an AP-IMEX scheme for the collision operator presented by Brull in [16]

$$Q_p = \lambda (G_p - f_p), \qquad (9)$$

where $\lambda$ is a free parameter that can be used to fit the transport coefficients, and

$$G_p = \frac{n_p}{\sqrt{\det(2\pi\tilde{\tau}_p)}} \exp\left(-\frac{1}{2}(v-u)^{\top}\tilde{\tau}_p^{-1}(v-u)\right), \qquad (10)$$

and the corrected tensors

$$\left\{ \begin{array}{l} \tilde{\tau}_p = \frac{1}{n_p} \int (v-u) \otimes (v-u) G_p \, dv = \tau \frac{\rho}{n} \frac{1}{m_p}, \\ \tau = \nu\Theta + (1-\nu)T \frac{n}{\rho} I, \end{array} \right. \qquad (11)$$

where $I$ is the identity matrix, the parameter $\nu$ is used to modify the value of the Prandtl number for $\nu \in [-1/2, 1)$ and the following quantity has been introduced

$$\Theta = \frac{1}{\rho} \sum_{p=1}^{P} \int m_p (v-u) \otimes (v-u) f_p \, dv. \qquad (12)$$

This operator [16] is an extension of the Ellipsoidal Statistical Model for the single species Boltzmann equation [21, 32], and its derivation is based on the resolution of a minimization problem, given some relaxation coefficients associated to some moments [17–19]. In particular, in [18] the authors propose a new approach of the ES-BGK model based on an entropy minimization principle [38, 39, 48]: different relaxation rates related to some moments are defined, and they describe the way in which these moments vanish for the trend to equilibrium. This approach has been applied to polyatomic gases in [19].

The AP property of our scheme, which is an extension of the work of [26], is based on an Implicit-Explicit formulation of the problem. This particular structure, allows to solve explicitly both the transport term and the collision term, despite having discretized this last term implicitly. Thanks to the implicit structure, the choice of the timestep is not restricted anymore by the choice of the small parameter of the problem (the Knudsen number in our case), but only by the CFL transport condition. This particular formulation allows for a significant reduction in computational costs (as will be discussed in more detail at the end of the paper), since in stiff regimes the Knudsen number may become very small. In such cases, a fully explicit solver would be subject to a much more restrictive stability condition than the standard transport CFL constraint, as the timestep cannot be larger than the Knudsen number as order of magnitude. Consequently, a fully explicit scheme would result in a prohibitive large computational cost.

To construct a high-order solver for the two-dimensional ES-BGK system for gas mixtures in both physical and velocity space, we combine several advanced numerical ingredients into a coherent and fully Asymptotic-Preserving framework.

* a third-order IMEX Runge-Kutta time integrator for stiff kinetic equations [9, 44];

* a high-resolution CWENO3 finite-volume reconstruction that provides third-order spatial accuracy [40, 45];

* a conservative velocity discretization that preserves the macroscopic moments of each species, obtained by extending the weighted $L^2$ minimization technique of [11] to the multi-species setting;
* a multithreaded CPU implementation that ensures computational efficiency in high-dimensional phase space;

* the Asymptotic-Preserving IMEX formulation specifically developed in this work, which generalizes the construction of [26] to the multi-species ES-BGK model of [16].

To the best of our knowledge, this is the first work that combines a high-order AP IMEX time discretization with a multi-species ES-BGK collision operator in multiple space dimensions, together with a high-order spatial and velocity discretization and a parallel implementation capable of handling an arbitrary number of species.

The rest of the manuscript is organized as follows. In Section 2 we introduce the AP IMEX formulation based on [26] for the ES-BGK model for gas mixtures proposed by [16], we also prove the Asymptotic-Preserving property of the method, we discuss how to arbitrarily increase the temporal order, and then we present the techniques that we implemented for the conservation of macroscopic moments. In Section 3 we show some numerical tests used to validate the scheme, and in the last Section 4 we summarize the key aspects of this work.

# 2 An asymptotic-Preserving scheme for the ES-BGK equation for gas mixtures

As summarized in [26, 35, 36] a scheme for the Boltzmann equation can be defined AP (Asymptotic-Preserving) if

* it preserves the discrete analogy of the Chapman-Enskog expansion: when holding the mesh size and time step fixed and letting the Knudsen number goes to zero, then one recovers a suitable scheme for the limiting Euler equations;

* implicit collision is treated explicitly, or at least, implemented in a manner that is computationally more efficient than the use of Newton-type solvers for nonlinear algebraic systems.

## 2.1 A first-order IMEX strategy

The following AP-IMEX scheme for this multi-species ES-BGK model is strongly inspired by the one proposed by Filbet and Jin in [26], and we will present it in a three-dimensional velocity space. For each species, a first order IMEX scheme for the Boltzmann equation with this multi-species ES-BGK operator will read as

$$ \begin{cases} \frac{f_p^{n+1} - f_p^n}{\Delta t} + v \cdot \nabla_x f_p^n = \frac{\lambda^{n+1}}{\varepsilon} (G_p^{n+1} - f_p^{n+1}), \\ f_p^0(x, v) = f_p(0, x, v) \end{cases} \qquad (13) $$

The previous equation can be rearranged as

$$ f_p^{n+1} = \frac{\varepsilon}{\varepsilon + \lambda^{n+1} \Delta t} [f_p^n - \Delta t v \cdot \nabla_x f_p^n] + \frac{\lambda^{n+1} \Delta t}{\varepsilon + \lambda^{n+1} \Delta t} G_p^{n+1}. \qquad (14) $$

Integrating each of these last $P$ equations with respect to velocity, and using the conservative properties of the $Q_p$ then it is possible to compute the numerical density of each species at the next timestep from the transport term. Analogously, multiplying each of these last $P$ equations by $m_p v$ and by $m_p |v|^2$, integrating with respect to velocity, summing up all the equations with respect to the species and using the conservative properties of $Q_p$ then it is possible to compute the corresponding total macroscopic quantities of the fluid (average velocity and average temperature). Doing that we
obtain the following equations for $n_p^{n+1}$, $u^{n+1}$ and $E^{n+1}$ (and thus $T^{n+1}$) as follows

$$ egin{cases} n_p^{n+1} = \int (f_p^n - \Delta t \cdot \nabla_x f_p^n) \, dv, \\ \rho^{n+1} = \sum_{p=1}^{P} \int m_p (f_p^n - \Delta t \cdot \nabla_x f_p^n) \, dv = \sum_{p=1}^{P} m_p n_p^{n+1}, \\ \rho^{n+1} u^{n+1} = \sum_{p=1}^{P} \int m_p v (f_p^n - \Delta t \cdot \nabla_x f_p^n) \, dv, \\ E^{n+1} = \sum_{p=1}^{P} \int \frac{m_p}{2} |v|^2 (f_p^n - \Delta t \cdot \nabla_x f_p^n) \, dv, \end{cases} \quad (15) $$

and using the relation $E = \sum_p \int \frac{m_p}{2} |v|^2 f_p \, dv = \frac{1}{2} \rho |u|^2 + \frac{3}{2} n T$ then it is possible to compute $T^{n+1}$.

Unfortunately, this is not sufficient to compute the operators $G_p^{n+1}$. So let's introduce the following quantity

$$ \Sigma^{n+1} := \sum_{p=1}^{P} \int v \otimes v \, m_p f_p^{n+1} \, dv = \rho^{n+1} [\Theta^{n+1} + u^{n+1} \otimes u^{n+1}], \quad (16) $$

where we used equation (12) in this last equality. Let's note that

$$ \sum_{p=1}^{P} \int m_p (v \otimes v) G_p \, dv = \rho (\tau + u \otimes u). \quad (17) $$

So multiplying equation (14) by $m_p v \otimes v$, integrating with respect to velocity, summing up with respect to all species, and using equations (16) and (17) then one gets

$$ \begin{aligned} \Sigma^{n+1} = \frac{\varepsilon}{\varepsilon + \lambda^{n+1} \Delta t} \left( \Sigma^n - \Delta t \sum_{p=1}^{P} \int m_p v \otimes v \, v \cdot \nabla_x f_p^n \, dv \right) & \quad (18) \\ + \frac{\lambda^{n+1} \Delta t}{\varepsilon + \lambda^{n+1} \Delta t} \rho^{n+1} (\tau^{n+1} + u^{n+1} \otimes u^{n+1}) & . \end{aligned} $$

Using equation (11)

$$ \begin{aligned} \Sigma^{n+1} = \frac{\varepsilon}{\varepsilon + \lambda^{n+1} \Delta t} \left( \Sigma^n - \Delta t \sum_{p=1}^{P} \int m_p v \otimes v \, v \cdot \nabla_x f_p^n \, dv \right) & \quad (19) \\ + \frac{\lambda^{n+1} \Delta t}{\varepsilon + \lambda^{n+1} \Delta t} \rho^{n+1} \left( \nu \Theta^{n+1} + (1-\nu) T^{n+1} \frac{n^{n+1}}{\rho^{n+1}} I + u^{n+1} \otimes u^{n+1} \right) & . \end{aligned} $$

The quantity $\Theta^{n+1}$ can be replaced by $\Sigma^{n+1}/\rho^{n+1} - u^{n+1} \otimes u^{n+1}$ (see equation (16)), and by factoring out the therm in $\Sigma^{n+1}$ on the left-hand side the following equation is obtained

$$ \begin{aligned} \Sigma^{n+1} = \frac{\varepsilon}{\varepsilon + \lambda^{n+1} \Delta t (1-\nu)} \left( \Sigma^n - \Delta t \sum_{p=1}^{P} \int m_p v \otimes v \, v \cdot \nabla_x f_p^n \, dv \right) & \quad (20) \\ + \frac{\lambda^{n+1} \Delta t (1-\nu)}{\varepsilon + \lambda^{n+1} \Delta t (1-\nu)} [T^{n+1} n^{n+1} I + \rho^{n+1} u^{n+1} \otimes u^{n+1}] & . \end{aligned} $$

The quantity in round brackets on the right-hand side can be directly computed since the distribution functions $f_p^n$ are known, while the macroscopic moments in square brackets on the right-hand side can be computed explicitly from $f_p^n$ using equation (15). By virtue of equation (20) it is possible to compute $\Sigma^{n+1}$, and then the quantity $\Theta^{n+1}$ by using equation (16). After that $\tau^{n+1}$ and $\tilde{\tau}_p^{n+1}$ are obtained thanks to equation (11) and in the end $G_p^{n+1}$ (equation (10)).

At this point it is possible to update each distribution function by computing its value at the next timestep by using equation (14). It is important to remark that this IMEX strategy (based on [26]) allows to solve explicitly the equation (13) which is nonlinearly implicit. In the following section we will discuss per AP property of this scheme.
## 2.2 Asymptotic-Preserving property

Now we prove the Asymptotic-Preserving property of the scheme similarly to [26].

**Proposition 1.** Consider the solution $f_p^{n+1}$ computed as indicated in the last section 2.1 then,

1. For all $\varepsilon \to 0$ and $\Delta t > 0$, the distribution functions $f_p^{n+1}$ satisfy

$$0 \le f_p^{n+1}(x, v) \le \max(||f_p^n||_\infty, ||G_p^{n+1}||_\infty).$$

2. For all $\Delta t > 0$ and $f_p^0$, the distribution functions $f_p^n$ converge to the equilibrium Maxwellians $M_p^{n^n, u^n, T^n}$, that is

$$\lim_{\varepsilon \to 0} f_p^n = M_p^{n^n, u^n, T^n},$$

and the scheme gives a first order approximation in time of the compressible Euler system for gas mixtures.

3. Assuming that $||f_p(t^n) - M_p(t^n)|| = O(\varepsilon)$, for $n \ge 2$ and

$$\left\| \frac{U^{n+1} - U^n}{\Delta t} \right\| \le C, \qquad (21)$$

(where $U^n$ are the numerical density of each species $n_p^n$, the total velocity $u^n$ and the total energy $E^n$ of the mixture computed at time $t^n$) the scheme asymptotically becomes a first order in time approximation of the compressible Navier-Stokes for gas mixtures given by

$$\left\{ \begin{array}{l} \displaystyle \frac{n_p^{n+1} - n_p^n}{\Delta t} + \nabla_x \cdot (n_p^n u^n) = -\varepsilon \nabla_x \cdot (J_p^{n-1}), \quad \forall p, \\[2ex] \displaystyle \frac{\rho^{n+1} \Delta t}{u^{n+1}} - \rho^n u^n + \nabla_x \cdot (\rho^n u^n \otimes u^n) + \nabla_x (n_p^n T^n) = -\varepsilon \nabla_x \cdot (\rho^{n-1} \Theta^{n-1}), \\[2ex] \displaystyle \frac{E^{n+1} \Delta t}{\Delta t} + \nabla_x \cdot ([E^n + n_p^n T^n] u^n) = -\varepsilon \nabla_x \cdot (\mathbb{Q}^{n-1} + \rho^{n-1} \Theta^{n-1} u^{n-1}), \end{array} \right. \qquad (22)$$

where $\Theta^n$ is given by equation (12), while $J_p^n$ and $\mathbb{Q}^n$ at time $t^n$ are

$$\left\{ \begin{array}{l} \mathbb{Q}^n = \sum_{p=1}^{P} \int m_p f_p^n \frac{|v - u^n|^2}{2} (v - u^n) \, dv, \\ J_p = \int v f_p^n \, dv, \end{array} \right. \qquad (23)$$

and they will be computed in this section providing their analytical expression.

*Proof* To prove the first property, let's note that $f_p^{n+1}$ is a linear combination of $f_p^n$ and $G_p^{n+1}$ thus we get the first assertion.

To prove the second property let us consider a generic initial data and $f_p^n$ for $n \ge 1$ and let us compute the asymptotic limit of $\Sigma^n$ when $\varepsilon$ goes to zero in equation (20). It yields

$$\Sigma^n = T^n n^n \mathbf{I} + \rho^n u^n \otimes u^n, \qquad (24)$$

and using (16)

$$\rho^n \Theta^n = T^n n^n \mathbf{I}, \qquad (25)$$

and from (11)

$$\tilde{\tau}_p^n = T^n \frac{1}{m_p} \mathbf{I}. \qquad (26)$$

So in the asymptotic limit $\varepsilon \to 0$ the quantities $G_p^n$ become isotropic, which means that they converge to the corresponding equilibrium Maxwellians $M_p^{n^n, u^n, T^n}$. So the solution at zeroth order is obtained
by taking $f_p^n = \mathcal{M}_p^{n_p^n, u^n, T^n}$ in the conservation laws (15)

$$
\left\{
\begin{aligned}
n_p^{n+1} &= \int \left( \mathcal{M}_p^{n_p^n, u^n, T^n} - \Delta t \, v \cdot \nabla_x \mathcal{M}_p^{n_p^n, u^n, T^n} \right) \, dv, \\
\rho^{n+1} &= \sum_{p=1}^{P} \int m_p \left( \mathcal{M}_p^{n_p^n, u^n, T^n} - \Delta t \, v \cdot \nabla_x \mathcal{M}_p^{n_p^n, u^n, T^n} \right) \, dv = \sum_p m_p n_p^{n+1}, \\
\rho^{n+1} u^{n+1} &= \sum_{p=1}^{P} \int m_p v \left( \mathcal{M}_p^{n_p^n, u^n, T^n} - \Delta t \, v \cdot \nabla_x \mathcal{M}_p^{n_p^n, u^n, T^n} \right) \, dv, \\
E^{n+1} &= \sum_{p=1}^{P} \int \frac{m_p}{2} |v|^2 \left( \mathcal{M}_p^{n_p^n, u^n, T^n} - \Delta t \, v \cdot \nabla_x \mathcal{M}_p^{n_p^n, u^n, T^n} \right) \, dv,
\end{aligned}
\right.
\qquad (27)
$$

from which one gets the following approximation of the multi-species Euler equation

$$
\left\{
\begin{aligned}
& \frac{n_p^{n+1} - n_p^n}{\rho^{n+1} u^{n+1}} + \nabla_x \cdot (n_p^n u^n) = 0, \quad \forall p \\
& \frac{\rho^{n+1} u^{n+1} - \rho^n u^n}{\Delta t} + \nabla_x \cdot (\rho^n u^n \otimes u^n) + \nabla_x (n^n T^n) = 0, \\
& \frac{E^{n+1} - E^n}{\Delta t} + \nabla_x \cdot \left[ \left( \frac{1}{2} \rho^n |u^n|^2 + \frac{5}{2} n^n T^n \right) u^n \right] = 0.
\end{aligned}
\right.
\qquad (28)
$$

This concludes the proof of the second assertion.

Now let's prove the third property regarding the asymptotically convergence towards the Navier-Stokes equation, using the Chapman-Enskog expansion, which consists of expanding each distribution function $f_p^n$ as

$$f_p^n = \mathcal{M}_p^n + \varepsilon h_p^n, \qquad (29)$$

from which one gets

$$
\left\{
\begin{aligned}
& \int h_p^n \, dv = 0 \quad \forall p, \\
& \sum_{p=1}^{P} \int m_p v h_p^n \, dv = \sum_{p=1}^{P} \int m_p |v|^2 h_p \, dv = 0.
\end{aligned}
\right.
\qquad (30)
$$

Similarly we also expand

$$\Theta^n = \frac{1}{\rho^n} \sum_{p=1}^{P} \int m_p (v - u^n) \otimes (v - u^n) f_p^n \, dv = \frac{T^n n^n}{\rho^n} I + \varepsilon \Theta_1^n, \qquad (31)$$

and

$$\mathbb{Q}^n = \sum_{p=1}^{P} \int m_p \frac{|v - u^n|^2}{2} (v - u^n) f_p^n \, dv = 0 + \varepsilon \mathbb{Q}_1^n, \qquad (32)$$

where

$$
\left\{
\begin{aligned}
\rho^n \Theta_1^n &= \sum_{p=1}^{P} \int m_p (v-u^n) \otimes (v-u^n) h_p^n \, dv, \\
\mathbb{Q}_1^n &= \sum_{p=1}^{P} \int m_p \frac{|v-u^n|^2}{2} (v-u^n) h_p^n \, dv,
\end{aligned}
\right.
\qquad (33)
$$

where $\Theta_1^n$ is traceless, because of equation (30).

Let us note that

$$J_{p,1} = \int h_p v \, dv = \int h_p (v-u) \, dv, \qquad (34)$$

where in the last equality we use equation (30). We can use these quantities inside the conservations
laws (15), obtaining

$$ egin{cases} \frac{n_p^{n+1} - n_p^n}{\rho^{n+1} u^{n+1}} + \nabla_x \cdot (n_p^n u^n) = -\epsilon \nabla_x \cdot (J_p^n), & \forall p, \\ \frac{\rho^{n+1} \Delta t}{u^{n+1}} + \nabla_x \cdot (\rho^n u^n \otimes u^n) + \nabla_x (n_p^n T^n) = -\epsilon \nabla_x \cdot (\rho^n \Theta_1^n), \\ \frac{E^{n+1} \Delta t}{\Delta t} + \nabla_x \cdot \left[ \left( \frac{1}{2} \rho^n |u^n|^2 + \frac{5}{2} n_p^n T^n \right) u^n \right] = -\epsilon \nabla_x \cdot (\mathbb{Q}_1^n + \rho^n \Theta_1^n u^n), \end{cases} \quad (35) $$

In the same way, we also expand each anisotropic Gaussian

$$ G_p^n = \mathcal{M}_p^n + \epsilon g_p^n. \quad (36) $$

From the definition of $\tau$ (11) and (31) one gets that

$$ \tau^n = T^n \frac{n^n}{\rho^n} \mathrm{I} + \nu \epsilon \Theta_1^n. \quad (37) $$

Using equation (11)

$$ \tilde{\tau}_p^n = \frac{T^n}{m_p} \mathrm{I} + \nu \epsilon \Theta_1^n \frac{\rho^n}{m_p n^n}, \quad (38) $$

whose determinant is equal to $(\frac{T^n}{m_p})^{D_v} + \mathcal{O}(\epsilon^2)$ [26], and the inverse is [26]

$$ [\tilde{\tau}_p^n]^{-1} = \frac{m_p}{T^n} \left[ \mathrm{I} - \frac{\nu \epsilon \rho^n}{n^n T^n} \Theta_1^n \right] + \mathcal{O}(\epsilon^2). \quad (39) $$

The quantity $g_p^n$ can be computed as

$$ g_p^n = \frac{G_p^n - \mathcal{M}_p^n}{\epsilon} = \frac{1}{2} \frac{m_p}{n^n (T^n)^2} \nu \rho^n (v - u^n)^{\tau} \Theta_1^n (v - u^n) \mathcal{M}_p^n. \quad (40) $$

Considering the terms of order zero into the scheme (14), one get an expression for $h_p^n$

$$ h_p^n = g_p^n - \frac{\mathcal{M}_p^n - \mathcal{M}_p^{n-1}}{\lambda^n \Delta t} - \frac{v}{\lambda^n} \nabla_x \mathcal{M}_p^{n-1} + \mathcal{O}(\epsilon). \quad (41) $$

From [7] one gets

$$ (\partial_t + v \nabla_x) \mathcal{M}_p = \mathcal{M}_p \left[ m_p c \cdot \left( \frac{\nabla_x n_p}{\rho_p} - \frac{\nabla_x n}{\rho} \right) + \frac{m_p}{T} \left( c \otimes c - \frac{1}{3} |c|^2 \mathrm{I} \right) : \nabla_x u + \left( \frac{m_p}{2T} |c|^2 - \frac{3}{2} - \frac{m_p n}{\rho} \right) c \cdot \frac{\nabla_x T}{T} \right], \quad (42) $$

where $c = v - u$.

Combining equations (40) and (42) into equation (41), and using assumption (21), it is possible to get an expression for $h_p^n$

$$ \begin{aligned} h_p^n = & -\frac{1}{\lambda^n} \mathcal{M}_p^{n-1} \left[ m_p c^{n-1} \cdot \left( \frac{\nabla_x n_p^{n-1}}{\rho_p^{n-1}} - \frac{\nabla_x n^{n-1}}{\rho^{n-1}} \right) + \frac{m_p}{T^{n-1}} \left( c^{n-1} \otimes c^{n-1} - \frac{1}{3} |c^{n-1}|^2 \mathrm{I} \right) : \nabla_x u^{n-1} \right. \\ & \left. + \left( \frac{m_p}{2T^{n-1}} |c^{n-1}|^2 - \frac{3}{2} - \frac{m_p n^{n-1}}{\rho^{n-1}} \right) c^{n-1} \cdot \frac{\nabla_x T^{n-1}}{T^{n-1}} \right] \\ & + \frac{1}{2} \frac{m_p}{n^n (T^n)^2} \nu \rho^n (c^n)^{\tau} \Theta_1^n (c^n) \mathcal{M}_p^n + \mathcal{O}(\Delta t) + \mathcal{O}(\epsilon). \end{aligned} \quad (43) $$
It's now possible to compute $\Theta_1^n$ and $\mathbb{Q}_1^n$ at the first order of the Chapman-Enskog expansion, by substituting this last expression for $h_p^n$ into equation (33).

$$ \rho^n \Theta_1^n = \sum_p \int m_p (v-u^n) \otimes (v-u^n) h_p^n dv = -\frac{n^{n-1}T^{n-1}}{\lambda^n} \sigma(u^{n-1}) + \nu \Theta_1^n \rho^n + \mathcal{O}(\Delta t) + \mathcal{O}(\varepsilon), \quad (44) $$

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

Rearranging the terms one gets

$$ \rho^n \Theta_1^n = -\frac{n^{n-1}T^{n-1}}{(1-\nu)\lambda^n} \sigma(u^{n-1}) = -\mu\sigma(u^{n-1}) + \mathcal{O}(\Delta t) + \mathcal{O}(\varepsilon), \quad (45) $$

where we used the fact that the viscosity [16] is given by

$$ \mu = \frac{nT}{(1-\nu)\lambda}. \quad (46) $$

Similarly

$$ \mathbb{Q}_1^n = \sum_{p=1}^{P} \int m_p \frac{|v-u^n|^2}{2} (v-u^n) h_p^n dv = \frac{5}{2} T^{n-1} \left( \sum_p J_{p,1}^n \right) - k \nabla_x T^{n-1} + \mathcal{O}(\Delta t) + \mathcal{O}(\varepsilon), \quad (47) $$

where the thermal conductivity $k$, consistently with [16], is

$$ k = \frac{5}{2} \frac{1}{\lambda^n} \left( \sum_p \frac{n_p^{n-1}}{m_p} \right) T^{n-1}, \quad (48) $$

and

$$ J_{p,1}^n = -\frac{n_p^{n-1}}{\lambda^n} \left( \frac{\nabla_x n_p^{n-1}}{\rho_p^{n-1}} - \frac{\nabla_x n^{n-1}}{\rho^{n-1}} \right) T^{n-1} - \frac{n_p^{n-1}}{\lambda^n} \left( \frac{1}{m_p} - \frac{n^{n-1}}{\rho^{n-1}} \right) \nabla_x T^{n-1} + \mathcal{O}(\Delta t) + \mathcal{O}(\varepsilon). \quad (49) $$

The Prandtl Number is related to the coefficient $\nu$ by [16, 26]

$$ \mathrm{Pr} = \frac{5 \mu}{2 k} = \frac{n^{n-1}}{1 - \nu \sum_p \frac{n_p^{n-1}}{m_p}}. \quad (50) $$

## 2.3 High-order IMEX schemes

The IMEX scheme showed in the previous section is first order in time. However it can be extended to an arbitrary higher temporal order using, for example, the IMEX DIRK (Diagonally Implicit Runge Kutta) schemes proposed by [44], that we briefly review in the following (much more details on the topic can be found in [9]).

The approach presented in [44] is adoptable for any stiff systems of differential equation in the form

$$ y' = q(y) + \frac{1}{\varepsilon} g(y), \quad (51) $$

where $y = y(t) \in \mathbb{R}^N$, $q, g: \mathbb{R}^N \to \mathbb{R}^N$. In our case $y$ is one of the distribution functions of the mixture, when considering only the temporal dependence, $q(y)$ is the transport part of the Boltzmann equation, and $g(y)$ is the ES-BGK operator. In order to extend this procedure to the gas mixture, then we have to apply this scheme to each distribution function.

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

$$ \left\{ \begin{array}{l} 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}{\varepsilon} 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}{\varepsilon} g(t_n + c_i \Delta t, Y_i). \end{array} \right. \quad (52) $$
The matrices $\tilde{A} = (\tilde{a}_{ij})$, $\tilde{a}_{ij} = 0$ for $j \ge i$ and $A = (a_{ij})$ are $\Lambda \times \Lambda$ matrices such that the 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}$, and it can be represented by a double tableau in the Butcher notation,

$$ \begin{array}{c|c} \tilde{c} & \tilde{A} \\ \hline & \tilde{\omega}^{\mathrm{T}} \end{array}, \quad \begin{array}{c|c} c & A \\ \hline & \omega^{\mathrm{T}} \end{array}. $$

In our code we implement the third order time IMEX DIRK scheme ARS(2,3,3) [44], which is characterized by the following Butcher tableaus

$$ \begin{array}{c|ccc} 0 & 0 & 0 & 0 \\ \gamma & \gamma & 0 & 0 \\ 1-\gamma & \gamma-1 & 2-2\gamma & 0 \\ \hline 0 & 1/2 & 1/2 & 0 \end{array}, \quad \begin{array}{c|ccc} 0 & 0 & 0 & 0 \\ \gamma & 0 & \gamma & 0 \\ 1-\gamma & 0 & 1-2\gamma & \gamma \\ \hline 0 & 1/2 & 1/2 & 0 \end{array} $$

where $\gamma = (3 + \sqrt{3})/6$.

## 2.4 Conservation of macroscopic quantities

It is well known that the discretization of the velocity space and the introduction of velocity grid introduce a loss of conservation [24]. A natural solution to the problem is to consider a velocity support big enough such that the distribution function is smooth and that the energy outside the domain in negligible. However, it may be considerable to construct a solver where conservation properties are exactly maintained during the simulation.

Several approaches have been proposed in literature, exploiting an Entropy minimization (like in [43]), or a $L^2$-minimization (like in [11, 29]). The Entropy approach allows the construction of a conservative discrete Maxwellian, while $L^2$-minimization can be applied to more general distribution functions. In this paper we implement, describe and extend to a multi-species context, the weighted $L^2$-minimization approach proposed in [11]. In particular, we generalize and extended this method to a multi-species context applying it to each anisotropic Maxwellian $G_p$ computed with our scheme. We will present the method considering an arbitrary quantity of species, and in a spatially homogeneous settings. Generalizing to a spatially non homogeneous settings is straightforward: it is sufficient to apply the following strategy cell-by-cell in the space domain. In the following we will consider a space of dimension 3 in velocity, but of course it can be generalized to an arbitrary amounts of dimensions in the velocity space.

Let $\tilde{G}_p$ be the initial guess for the anisotropic Maxwellian of the $p$-th species obtained with the previous AP-IMEX scheme. Let us consider a different generic weight function $1/h_p$ for each species (we will specify in the following the choice of the weight functions), then we look for the target anisotropic Maxwellians $G_p$ exactly preserving the moments. Let $M$ be the total number of discretization points of the velocity space such that

$$ \tilde{G}_p = \begin{pmatrix} \tilde{G}_p^1 \\ \tilde{G}_p^2 \\ \vdots \\ \tilde{G}_p^M \end{pmatrix}. \qquad (53) $$

Let $Z_p$ be the $(D_V + 2) \times M$ (in this case $D_V = 3$, number of dimensions of the velocity space) integration matrix

$$ Z_p = \begin{pmatrix} w_1(h_p)_1 & w_2(h_p)_2 & \dots & w_M(h_p)_M \\ m_p v_1^x w_1(h_p)_1 & m_p v_2^x w_2(h_p)_2 & \dots & m_p v_M^x w_M(h_p)_M \\ m_p v_1^y w_1(h_p)_1 & m_p v_2^y w_2(h_p)_2 & \dots & m_p v_M^y w_M(h_p)_M \\ m_p v_1^z w_1(h_p)_1 & m_p v_2^z w_2(h_p)_2 & \dots & m_p v_M^z w_M(h_p)_M \\ 0.5 m_p |v_1|^2 w_1(h_p)_1 & 0.5 m_p |v_2|^2 w_2(h_p)_2 & \dots & 0.5 m_p |v_M|^2 w_M(h_p)_M \end{pmatrix}, \qquad (54) $$
where $w_i$ are the integration weights.

The moments of the p-th species are the following

$$a_p = \begin{pmatrix} n_p \\ \rho_p u_p^x \\ \rho_p u_p^y \\ \rho_p u_p^z \\ E_p \end{pmatrix}, \qquad (55)$$

where $n_p$ is the numerical density of the p-th species, $\rho_p$ is its mass density, $u_p^x, u_p^y, u_p^z$ are respectively, the total velocities of the p-th species in the three directions and $E_p$ is its energy.

The goal is to find $G_p$ that minimizes the weighted $L^2$ distance from $\tilde{G}_p \circ \frac{1}{h_p}$ under the constraint that $Z_p G_p - a_p = 0$, where the symbol $\circ$ is the Hadamard product.

Similarly to [11], we introduce the following Lagrangian for each species

$$L_p = \left\| \tilde{G}_p \circ \frac{1}{h_p} - G_p \right\|_2^2 + \lambda_p^T (Z_p G_p - a_p). \qquad (56)$$

From the stationary points of $L_p$ one gets an expression for $\lambda_p$ from which [11]

$$G_p \circ h_p = \tilde{G}_p + Z_p^T (Z_p Z_p^T)^{-1} \left( a_p - Z_p \left( \tilde{G}_p \circ \frac{1}{h_p} \right) \right) \circ h_p. \qquad (57)$$

Then in order to impose the conservation of moments we applied this approach to each anisotropic Maxwellian $\tilde{G}_p$ computed using our scheme. In particular, at each step of the Runge-Kutta scheme, we compute the anisotropic Maxwellians $\tilde{G}_p$. Then we use this $L^2$-minimization strategy, to estimate the anisotropic Maxwellians (the quantity $G_p$ in this last approach) having exactly the desired moments used to construct $\tilde{G}_p$ using our scheme. Similarly to [11], we impose $h_p$ to be equal to the Maxwellian $M_p$ having the moments used to construct $\tilde{G}_p$.

# 3 Numerical tests

In the following, we present several simulations conducted to numerically test and validate our numerical scheme. All simulations have been conducted on a non-dedicated computer equipped with a 12th Gen Intel(R) Core(TM) i7-12700H CPU and 16 GB of RAM, except the Test 2 which has been executed on a dedicated computer equipped with two Intel(R) Xeon(R) Gold 5418Y CPUs and 128 GB of RAM, kindly provided by the Mathematics Department of the University of Ferrara. In order to speed-up the code execution a multithread approach has been implemented for the non homogeneous simulations 2D in space: each thread operates on a subdomain of the entire spatial domain. For the spatially non-homogeneous simulations the variables in the transport term have been reconstructed using a CWENO 3 scheme [40, 45], while in time we use the IMEX scheme presented before extended to the third order using the ARS(2,3,3) scheme [44]. For the sake of simplicity, in all simulations we considered $\nu = -0.5$ and $\lambda = 1$, except where clearly indicated.

In Section 3.1 we perform some numerical experiments to validate the scheme and our numerical implementation. In Section 3.2 we show some interesting numerical applications for fluid simulations 1D and 2D in space, like the classical SOD shock tube benchmark test, or some intriguing phenomena, like the Kelvin-Helmholtz instability and the flow of a fluid around a cylinder.

## 3.1 Validation

In this section we perform some numerical experiments to validate the scheme and our numerical implementation. We perform a first test in order to verify the moments relaxation of two species toward those of equilibrium. In a second numerical experiment we study the convergence rate in time and both in space and time of our code.
### 3.1.1 Homogeneous relaxation test

The aim of this spatially homogeneous test 1D in velocity is to check the correct convergence of the temperature of each species toward the equilibrium one. We have considered two species of masses $m_1$ and $m_2$ described by the following distribution functions $f_1$ and $f_2$ both initialized as sum of two Maxwellians

$$ \left\{ \begin{aligned} f_1(0, v) &= \mathcal{M}_1^{n_a, v_a, T_a} + \mathcal{M}_1^{0.3n_a, 1.5v_a, 1.3T_a}, \\ f_2(0, v) &= \mathcal{M}_2^{n_b, v_b, 0.7T_b} + \mathcal{M}_2^{3n_b, 0.5v_b, 2T_b}, \\ n_a &= 1/m_1, \quad T_a = 1, \quad v_a = 0.5, \\ n_b &= 2/m_2, \quad T_b = 2, \quad v_b = -0.3. \end{aligned} \right. \qquad (58) $$

It is important to remark that the macroscopic quantities $n_a$, $n_b$, $v_a$, $v_b$, $T_a$ and $T_b$ do not correspond to the macroscopic moments of the two species $n_1$, $n_2$, $v_1$, $v_2$, $T_1$ and $T_2$, since both species are sum of two different Maxwellians. Also, the Maxwellians $\mathcal{M}_p$ are computed considering the mass of the p-th species. Thus in the expression of the Maxwellians $\mathcal{M}_1$, the mass $m_1$ is employed; similarly, in the expression of the Maxwellians $\mathcal{M}_2$, the mass $m_2$ is employed.

The simulations have been conducted in two different configurations of mass ratio (MR): configuration 1 with $m_2/m_1 = 1$ and configuration 2 with $m_2/m_1 = 100$. In all simulations and both configurations we considered $m_1 = 1$, $\epsilon = 10^{-3}$ (Knudsen Number), $\Delta t = 10^{-4}$ (time step). The small value of $\Delta t$ is not due to any stability condition, but only in order to better graphically visualize the convergence of the parameters of each species towards the one of the equilibrium Maxwellian. The final time of the integration is equal to $16\epsilon$. We considered a velocity box $[-20, 20]$ for both configurations, discretized with a total of 256 points when $MR = 100$ and 32 points when $MR = 1$. In Figure 1, the temperature and velocity evolution in time of both species are plotted, and it is evident that they are converging toward the corresponding asymptotic value.

The time integration in this simulation has been performed using the third-order IMEX strategy presented above, without considering the transport term in the Boltzmann equation.

### 3.1.2 Convergence test

**Time convergence test** The aim of this test is to check the third-order temporal convergence of the implemented ARS(2,3,3) IMEX scheme. We considered two gases with a given initial condition 1D in space and 1D in velocity, and we evolved the solution using the IMEX scheme presented above, using different timesteps $(\Delta t)_i$ up to a time $T_f = 1$ (where $(\Delta t)_{i+1} = 0.5(\Delta t)_i$). Then we evaluated the error $\delta_{(\Delta t)_i}$ associated to each $(\Delta t)_i$ at the final integration time $T_f$ as following

$$ \delta_{(\Delta t)_i} = \| \rho_{(\Delta t)_i}(T_f) - \rho_{(\Delta t)_{rif}}(T_f) \|_2, \qquad (59) $$

where $\rho_{(\Delta t)_{rif}}$ is the solution computed using the timestep $(\Delta t)_{rif} = 0.5(\Delta t)_{min}$, and $(\Delta t)_{min}$ is the smallest timestep for which we compute the error convergence. We considered two gases of masses $m_2 = 2m_1 = 2$, initialised as

$$ \left\{ \begin{aligned} f_1 &= \mathcal{M}_1^{n_1(x), v(x), T(x)} \\ f_2 &= \mathcal{M}_2^{n_2(x), v(x), T(x)}, \end{aligned} \right. \qquad (60) $$

where $v(x) = 0$, $T(x) = 0.3(1 + \exp(-3|x - L_x/2|^2))$, $n_i(x) = \rho(x)/m_i$ and $\rho(x) = 1 + 0.1 \exp(-3|x - L_x/2|^2)$, where $L_x = 10$, and $x \in [0, L_x]$ discretized using 100 cells. The velocity space $[-10, 10]$ is discretized using 40 points, and we used periodic boundary conditions in the space domain. The convergence curves for different values of $\epsilon$ are reported in Figure 2. They clearly show the third order of convergence of the scheme.

**Space and time convergence test** Lastly we have also checked the global third-order accuracy of the scheme in space and in time. We considered two gases with a given initial condition 1D in space and 1D in velocity, and we evolved the solution using the IMEX scheme presented above, using different spatial discretization steps $(\Delta x)_i$ (obtained by doubling the number of discretization cells $N_{i+1} = 2N_i$). The timestep associated to each run is $(\Delta t)_i = (\Delta x)_i/L_v$ (where $L_v = 10$ is the maximum value of
<figure><img src="image_4.png" /><figcaption>(a) Temperature Evolution; Mass Ratio = 1</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(b) Velocity Evolution; Mass Ratio = 1</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(c) Temperature Evolution; Mass Ratio = 100</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(d) Velocity Evolution; Mass Ratio = 100</figcaption></figure>

Figure 1: Moments evolution for the homogeneous relaxation test. In the first row it is reported the evolution of the temperature and the velocity of each species for a mass ratio equal to 1. The same quantities, but for a mass ratio equal to 100, are reported in the second row. The solid line in each figure represents the asymptotic value associated to the corresponding equilibrium Maxwellian given by equation (5), while the crosses and the open circles represents, respectively, the first and the second species.
<figure><img src="image_8.png" /><figcaption>Figure 2: Convergence rate in time for different values of ε for the time convergence test.</figcaption></figure>

the velocity box [-10, 10] discretized using 40 points). Then we evaluated the error $\delta_{(\Delta x)_i}$ associated to each $(\Delta x)_i$ at the final integration time $T_f = 0.2$ as following

$$ \delta_{(\Delta x)_i} = \|\rho_{(\Delta x)_i}(T_f) - \rho_{(\Delta x)_{rif}}(T_f)\|_2, \qquad (61) $$

where $\rho_{(\Delta x)_{rif}}$ is the solution computed using the timestep $(\Delta x)_{rif} = 0.5(\Delta x)_{min}$, and $(\Delta x)_{min}$ is the smallest step for which we compute the error convergence. We considered two gases of masses $m_2 = 2m_1 = 2$, initialised as

$$ \begin{cases} f_1 = M_1^{n_1(x),v(x),T(x)} \\ f_2 = M_2^{n_2(x),v(x),T(x)}, \end{cases} \qquad (62) $$

where $v(x) = 0$, $T(x) = 0.3(1 + 0.1 \exp(-\Theta|x - L_x/2|^2))$, $n_i(x) = \rho(x)/m_i$ and $\rho(x) = 1 + 0.1 \exp(-\Theta|x - L_x/2|^2)$, where $\Theta = 1000$, $L_x = 1$, and $x \in [0, L_x]$. The first run $i = 1$ is the one obtained using 50 cells to discretize the physical space. We used periodic boundary conditions in the space domain. In Figure 3 the convergence curves for different values of ε are reported. They clearly show the global third order of convergence of the scheme.

## 3.2 Numerical applications

We show some interesting numerical applications for fluid simulations 1D and 2D in space, like the classical SOD shock tube benchmark test, or some intriguing phenomena, like the Kelvin-Helmholtz instability and the flow over a cylinder. Lastly, we perform a comparison between BGK, ES-BGK and full Boltzmann operator. This last test shows that the ES-BGK performs better than the BGK operator in the rarefied regime. At the end of this section we also discuss the big speed-up offered by such kind of IMEX scheme compared to fully explicit integrators.

### 3.2.1 Test 1: A Sod shock tube problem

A classical benchmark test for the 1D Euler equations is the so-called Sod shock tube test problem [50]: its solution consists of a left rarefaction, a contact and a right shock. Here we have considered its extension to the multi-species configuration using similar initial conditions proposed by [2]. We consider two species in a 1D domain $[0, L_x]$ where $L_x = 1$, discretized using $N_x = 200$ cells. In order to initially confine the first gas on the left side of the discontinuity, and the second gas on
<figure><img src="image_9.png" /><figcaption>Figure 3: Convergence rate in space and time for different values of ε for the space and time convergence test.</figcaption></figure>

the right side of the discontinuity we initialise the two distribution functions as Maxwellians whose moments are

$$ \begin{cases} \rho_1 = (1 - \delta)\rho_L, & \rho_2 = \delta\rho_L, & v_1 = v_2 = v_L, & T_1 = T_2 = P_L/(\rho_1/m_1), & \text{if } x < 0.5, \\ \rho_1 = \delta\rho_R, & \rho_2 = (1 - \delta)\rho_R, & v_1 = v_2 = v_R, & T_1 = T_2 = P_R/(\rho_2/m_2), & \text{if } x \ge 0.5, \end{cases} \quad (63) $$

where $\rho_L = 1$, $v_L = 0$, $P_L = 1$, $\rho_R = 2^{-3}$, $v_R = 0$, $P_R = 2^{-5}$.

As in [2], we set $\delta = 10^{-5}$, because its ideal value ($\delta = 0$), would lead to ill-defined macroscopic moments.

Each species has its own mass, thus we introduce the mass $m_1$ associated to the first species, and the mass $m_2$ associated to the second species.

In all simulations we considered $\Xi = 3$ as specific heat ratio in the computation of the analytical solution.

All simulations are executed up to a final integration time equal to 0.15.

We executed two sets of simulations and we compared them with the analytical solution formally obtained when $\epsilon \to 0$. This reference solution is obtained analytically solving the Riemann problem.

* Set A: two simulations with mass ratio equal to 1, and two different values of the Knudsen number, $\epsilon = 10^{-2}$ and $\epsilon = 10^{-6}$. The timestep is equal to $\Delta t = 0.1 dx$. The goal of this test is to check the correct implementation of the integration method and that the scheme is able to detect the correct asymptotic behaviour (Euler equations) when the Knudsen number goes to zero. In both simulations we considered a velocity domain given by $[-L, L]$ with $L = 8$ and discretized with 32 points. The density, velocity and temperature of the mixture at the final time are showed in Figure 4, from which it is evident that the solution computed with the smallest value of the Knudsen number is converging towards the hydrodynamic (analytic) solution ($\epsilon \to 0$) and the correctness of the numerical schemes.

* Set B: three simulations with different mass ratios: $MR = 1$ and $MR = 10$ and $MR = 20$, and with $\epsilon = 10^{-6}$. This is an important test to check that the scheme is properly able to numerically integrate those configurations where the two species have different masses. The velocity space of the simulation with $MR = 1$ is given by $[-L, L]$ with $L = 8$ and it is discretized with 32 points, and the timestep is equal to $\Delta t = 0.1 dx$; while for the simulations with $MR = 10$ and $MR = 20$ we considered $L = 40$ and 200 points, and the timestep is equal to $\Delta t = 0.025 dx$. The
<figure><img src="image_10.png" /><figcaption>Figure 4: Total density, temperature and velocity for Sod Shock Tube problem (Test 1), for different values of Knudsen number at the final time $t = 0.15$ and with unitary mass ratio.</figcaption></figure>

densities, velocities and pressure of the mixture are showed in Figure 5. It is evident that the code is able to correctly capture the correct asymptotic behaviour, even with high mass ratios. The bumps that can be seen in the numerical solutions have been already detected and discussed by other authors, like [2]: they are due to the effect of high mass ratio in proximity of the contact discontinuity between the two species. The choice of the velocity box can affect the magnitude of these bumps [2].

### 3.2.2 Test 2: Kelvin-Helmholtz instability

The Kelvin-Helmholtz instability is a phenomenon where vortices arise at the interface of two fluids moving at different horizontal speeds and with a perturbative vertical sinusoidally spatially distributed velocity.

The initial conditions are inspired from [41, 42] and have been adapted for Gas-Mixtures by Bailo and Rey in [2]. This is an important test to check the implementations of our 2D integration schemes for our kinetic model.

We consider two gases with mass ratio $m_2/m_1 = 2$ in a domain $[0, 1] \times [0, 0.5]$, periodic along x-direction, and discretized using $N_x = 150$ and $N_y = 75$ cells. The velocity domain $[-8, 8]^2$ is discretized using 32 points in both directions. The Knudsen number is equal to $\epsilon = 5 \times 10^{-5}$. The two gases are initially
<figure><img src="image_11.png" /><figcaption>Figure 5: Total Density, Pressure and Velocity in the Sod Shock Tube problem (Test 1), for different values of Mass Ratio ($MR = m_2/m_1$), with a Knudsen number equal to $\epsilon = 10^{-6}$, at the final time $t = 0.15$.</figcaption></figure>
separated along the interface surface placed at $y = L_y/2$ and they are initialized as Maxwellians whose moments correspond to the following values:

$$y > L_y/2: \begin{cases} \rho_1 = \tilde{\rho}_1(1-\delta), \quad v_1 = (0.5, 0.01 \sin(4\pi x)), \quad T_1 = (\tilde{\rho}_1/m_1)^{-1}, \\ \rho_2 = \tilde{\rho}_2\delta, \quad v_2 = (0.5, 0.01 \sin(4\pi x)), \quad T_2 = T_1 = (\tilde{\rho}_1/m_1)^{-1}, \end{cases} \quad (64)$$

$$y \le L_y/2: \begin{cases} \rho_1 = \tilde{\rho}_1\delta, \quad v_1 = (-0.5, 0.01 \sin(4\pi x)), \quad T_1 = T_2 = (\tilde{\rho}_2/m_2)^{-1}, \\ \rho_2 = \tilde{\rho}_2(1-\delta), \quad v_2 = (-0.5, 0.01 \sin(4\pi x)), \quad T_2 = (\tilde{\rho}_2/m_2)^{-1}, \end{cases} \quad (65)$$

where $\delta = 10^{-5}$, $\tilde{\rho}_1 = 1$ and $\tilde{\rho}_2 = 2$.

The densities of the two species for three different times $t = 0.9$, $t = 1.7$ and $t = 2.70$ are plotted in Figure 6. At the final timestep, it is clear the formation of the vortices, consistently with [2].

### 3.2.3 Test 3: Flow around a cylinder

The aim of this test is to verify that the numerical scheme is able to handle complex geometries, like a flow over a round object, similarly to [26, 52]. We consider a 2D domain of size $[0, L_x] \times [0, L_y]$ where $L_x = L_y = 0.8$, and discretized with 66 cells in both directions. The velocity space is a domain equal to $[-L, L]^2$ where $L = 8$ and it is discretized using 32 points in both directions. The Knudsen number is equal to $10^{-6}$ and we consider a mass ratio $MR = m_2/m_1 = 5$ ($m_1 = 1$ and $m_2 = 5$). The round object is modelled considering a round shape centred in the cell of indices $[4L_x/10, N_y/2]$ and whose radius is equal to $R = L_x/10$ in physical coordinates, such that the integrator is never executed inside this round shape. We implement specular boundary conditions [22] on the walls of the object. The distribution functions for both species are initialized like a Maxwellian whose moments outside of the round object are given by

$$\begin{cases} \rho_1 = 0.8, \quad v_1 = (0.8, 0), \quad T_1 = 1, \\ \rho_2 = 1.0, \quad v_2 = (0.8, 0), \quad T_2 = 1 \end{cases}. \quad (66)$$

In Figure 7 they are plotted the densities of the two species for three different times $t = 0.05$ and $t = 0.10$ and $t = 0.15$, from which it is evident that the code is able to handle complex geometries. The timestep is approximately equal to $\Delta t = 7.7 \times 10^{-4}$.

From this test it is also clear that the densities of the two species are different one from each other, due to their different masses and different initial values: thus it is remarkable the importance of using a multi-species framework rather than the single species one.

### 3.2.4 Test 4: Comparison between BGK, ES-BGK and full Boltzmann operator

The aim of this test, inspired by [26], is to show the ability of our multi-species ES-BGK scheme to reproduce the correct behaviour of the full Boltzmann operator.

For simplicity we consider the case of two species with same unitary masses $m_1 = m_2 = m = 1$. In what follows we strongly use this assumption, that it is also really useful, because allows us to check the indifferentiability principle [16] of our scheme: when all the species have the same masses and their crosses section are equal, then the system of equations reduces to a single one by adding the distribution functions.

Here we will follow and summarize [26] in presenting this test performed in a 1D space domain $[-L_x, L_x]$, with $L_x = 1$, and 2D velocity space. We will consider two species of same masses, and the distribution functions initialized as $f_1(t = 0) = f_2(t = 0)$ and evolved using the scheme presented in this paper.

Considering a smooth solution with periodic boundary conditions, in a 1D space domain $[-1, 1]$, from the conservation laws for mass, momentum and energy, one can determine the stationary state of the distribution function $f = f_1 + f_2$, which is the normalized global Maxwellian distribution [26]

$$M_g(v) = \frac{\rho_g}{2\pi T_g} \exp\left(-\frac{|v|^2}{2T_g}\right), \quad (67)$$

where $\rho_g$ and $T_g$ to be determined, under the assumption that

$$\frac{1}{2} \int_{-1}^{1} \int_{\mathbb{R}^2} f_1(t=0) v \, dv \, dx = \frac{1}{2} \int_{-1}^{1} \int_{\mathbb{R}^2} f_2(t=0) v \, dv \, dx = 0. \quad (68)$$
<figure><img src="image_12.png" /><figcaption>Figure 6: Densities of the two species at $t = 0.9$, $t = 1.7$ and $t = 2.7$ for the Kelvin-Helmholtz instability (Test 2).</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>(a) Species 1, $\rho_1$</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(b) Species 2, $\rho_2$</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(c) Species 1, $\rho_1$</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(d) Species 2, $\rho_2$</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(e) Species 1, $\rho_1$</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(f) Species 2, $\rho_2$</figcaption></figure>

Figure 7: Densities evolution of the two species at $t = 0.05$, $t = 0.10$ and $t = 0.15$ for the flow around a circular shape (Test 3).
Under some boundness assumption, it is expected the distribution function $f$ converges to the global Maxwellian $M_g$ as $t \to \infty$. One expects that the if the solution is smooth enough then [26]

$$||f(t) - M_g|| = O(t^{-\infty}), \qquad (69)$$

meaning that the solution converges almost exponentially fast to the global equilibrium. In order to check this convergence behaviour we study the evolution of the following quantity

$$||\rho(t, x) - \rho_g||_{L_x^1} = \int_{-1}^{1} |\rho(t, x) - \rho_g| dx, \qquad (70)$$

where $\rho$ is the total mass density of the distribution function $f$ (where $f = f_1 + f_2$) and

$$\rho_g = \frac{1}{2} \int_{-1}^{1} \int_{\mathbb{R}^2} (f_1(t=0, x, v) + f_2(t=0, x, v)) \, dv \, dx. \qquad (71)$$

In particular we will compare the evolution of this quantity obtained from three different approaches:

* multi-species ES-BGK model presented in this paper by setting $\nu = -1$ in our scheme,

* multi-species BGK model obtained by putting $\nu = 0$ in our scheme,

* spectral full Boltzmann operator for single species [27].

We use the first two integrators to evolve the multi-species system given by $f_1$ and $f_2$, initialized using the same distribution function and having the same unitary mass; and we use the third one to evolve $f$ which is initialized as $f(t=0) = f_1(t=0) + f_2(t=0)$. Indeed thanks to the indifferentiability principle the total density $\rho$ obtained from the multi-species model ($\rho = \rho_1 + \rho_2$) and the total density $\rho$ obtained from the single-species such defined, are the same.

In order to obtain the correct full Boltzmann Prandtl number with the ES-BGK operator we set $\nu = -1$ [26]. By equating the viscosity of the Boltzmann operator with the one of the ES-BGK operator (see equation (46)) one gets a constraint on the parameter $\lambda$ [26], that is $\lambda = \frac{P}{1-\nu} \frac{A_2 3\pi}{\sqrt{2T}}$, where $A_2 = 0.436$. We use this equation to set the value of $\lambda$ for the ES-BGK operator (by setting $\nu = -1$) and for the BGK operator (by setting $\nu = 0$). We initialize the two distribution functions as

$$f_1(t=0, x, v) = f_2(t=0, x, v) = \frac{1 + A_0 \sin(\pi x)}{4(2\pi T_0)} \left[ \exp\left(-\frac{|v-u_0|^2}{2T_0}\right) + \exp\left(-\frac{|v+u_0|^2}{2T_0}\right) \right], \quad (72)$$

where $A_0 = 0.5$, $T_0 = 0.125$ and $u_0 = (1/2, 1/2)$. We update $f_1$ and $f_2$ using the multi-species ES-BGK operator, and the multi-species BGK operator, while we update $f$ (initialized as $f(t=0) = f_1(t=0) + f_2(t=0)$) using the full Boltzmann operator. We consider $L_x = 1$ with a periodic physical space $[-L_x, L_x]$ discretized using 100 cells, and a velocity box $[-8, 8]^2$ discretized using 64 points in both directions.

In all simulations the space reconstruction is performed using the already discussed CWENO3 scheme. We perform the test using two different values of the Knudsen number $\epsilon = 10^{-1}$ and $\epsilon = 10^{-4}$. For the sake of simplicity, the time integration of the ES-BGK and BGK operators is performed using the first order AP IMEX scheme discussed in this paper, while for the full Boltzmann operator we use the first order IMEX AP scheme proposed in [25]. The timestep is equal to $dt = 0.1dx$ in all simulations, and the solutions are evolved up to a final time $t_F = 10$.

The curves $||\rho(t, x)_{ESBGK} - \rho_g||_{L_x^1}$, $||\rho(t, x)_{BGK} - \rho_g||_{L_x^1}$, $||\rho(t, x)_{BOL} - \rho_g||_{L_x^1}$, reported in Figure 8, are obtained respectively using the multi-species ES-BGK operator, the multi-species BGK operator and the single species spectral full Boltzmann operator. An important indicator of the correctness our results, are the curves in Figure 8 obtained with $\epsilon = 0.1$, which are in perfect agreement with those presented in [26]. It is evident that the ES-BGK operator is able to correctly match and capture the behaviour of the curve obtained using the full Boltzmann operator even in the rarefied regime $\epsilon = 0.1$, while the BGK operator only in the collision dominated regime $\epsilon = 10^{-4}$. This test clearly shows the importance of using the ES-BGK operator rather than the BGK operator in the rarefied regime.
<figure><img src="image_19.png" /><figcaption>Figure 8: Evolution of the quantity $||\rho(t, x) - \rho_g||_{L_x^1}^{\text{norm}}$ (Test 4) obtained using three different operators: multi-species ES-BGK, multi-species BGK and single species spectral full Boltzmann, for two different values of Knudsen number $\epsilon$. The multi-species solvers are used to evolve the distribution functions $f_1$ and $f_2$ of two species that have the same mass and the same cross section, while the single-species Boltzmann operator is used to evolve the single distribution $f$ initialized as the sum of the two species $f(t=0) = f_1(t=0) + f_2(t=0)$. Then thanks to the indifferentiability principle it is possible to compare the evolution of the sum of the two distribution functions $f_1(t) + f_2(t)$ with the evolution of $f(t)$ for every time $t$.</figcaption></figure>

### 3.2.5 Discussion on the IMEX scheme computational performance

We now briefly discuss the speed-up offered by the proposed IMEX scheme when compared to a standard fully explicit integrator, using our numerical tests as a representative example. More details can be found in [9]. The first clear advantage of the IMEX scheme is the already discussed Asymptotic- Preserving property, which ensures that, for sufficiently small Knudsen numbers, the method becomes an exact solver for the multi-species Euler and Navier-Stokes limits (see Section 2.2). Moreover, the IMEX structure guarantees that the method can be implemented at the cost of a fully explicit scheme but with a timestep constrained only by the CFL condition associated with the transport term, and not by the value of the Knudsen number. In contrast, in a standard fully explicit solver the timestep must satisfy

$$ \Delta t^{\text{Explicit}} \lesssim \min\{\epsilon, \Delta x/v_{\text{max}}\}, $$

where $\epsilon$ is the Knudsen number, $v_{\text{max}}$ the maximal velocity in the velocity domain, and $\Delta x$ the spatial grid size.

To make this concrete, consider Test 3 (flow around a cylinder). In this simulation we used:

$$ \epsilon = 10^{-6}, \quad \Delta x = 0.0123, \quad v_{\text{max}} = 8. $$

Hence,

$$ \min\{\epsilon, \Delta x/v_{\text{max}}\} = \epsilon = 10^{-6}, $$

so a fully explicit solver would be forced to use a timestep of order $\Delta t^{\text{Explicit}} \approx 10^{-6}$.

In contrast, our IMEX scheme is constrained only by $\Delta x/v_{\text{max}}$, yielding a timestep of approximately

$$ \Delta t^{\text{IMEX}} = 7.7 \times 10^{-4}. $$

Therefore, the effective speed-up is

$$ \frac{\Delta t^{\text{IMEX}}}{\Delta t^{\text{Explicit}}} \approx 770, $$

which is extremely significant.

Of course, the precise value of the speed-up depends on both the Knudsen number and the CFL constraint, which are problem-dependent. Remarkably, in the worst-case scenario (when $\epsilon > \Delta x/v_{\text{max}}$),
the maximal allowable IMEX timestep coincides with that of the explicit scheme. Thus, IMEX schemes
are *always* at least as good as explicit ones, and often dramatically better.

This qualitative discussion highlights the importance of developing IMEX methods for accelerating
stiff kinetic simulations. Finally, let us mention that further reductions in computational time could
be achieved by introducing a local velocity grid adaptation strategy, such as the one proposed in [3].

4 Conclusion

We have developed a high-order Asymptotic-Preserving IMEX scheme for the ES-BGK model of gas mixtures. The construction extends IMEX strategies previously available for single-species BGK models to the multi-species ES-BGK collision operator and is formulated in such a way that the time discretization is explicitly implementable and uniformly stable with respect to the Knudsen number. In the fluid regime, the scheme reduces to a consistent and high-order accurate discretization of the limiting macroscopic equations for the mixture, while in the rarefied regime it remains a genuine kinetic solver.

The proposed method has been embedded in a high-order numerical framework that combines
third-order IMEX Runge-Kutta time integration, CWENO3 spatial reconstruction, conservation of
macroscopic moments in the discrete velocity space, and a multithreaded implementation. The nu-
merical tests confirm the expected accuracy of the scheme in both kinetic and fluid regimes, as well
as a smooth transition between them, and illustrate its capability to handle an arbitrary number of
species without loss of stability.

Several directions for future work are natural. On the modeling side, one can consider the extension of the present framework to more complex collisional models, and to mixtures in plasma physics. On the numerical side, it would be of interest to couple the present IMEX strategy with different high-order spatial discretization and unstructured meshes, and to investigate its performance in multi-scale applications where kinetic and fluid regions coexist in a more complex geometry.

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 acknowledges the support of the Mathematics Department of the University of Ferrara through access to their computational resources. 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". LP also acknowledges the partial support by Fondo Italiano per la Scienza (FIS2023-01334) advanced grant ADAMUS. This work has been written within the activities of GNCS group of INdAM (Italian National Institute of High Mathematics).

References

[1] P. Andries, K. Aoki, and B. Perthame. "A consistent BGK-type model for gas mixtures". In: *J. Stat. Phys.* 106.5-6 (2002), pp. 993–1018.

[2] R. Bailo and T. Rey. "Projective and telescopic projective integration for non-linear kinetic mixtures". In: *J. Comput. Phys.* 458 (2022), p. 111082.

[3] F. Bernard, A. Iollo, and G. Puppo. "A Local Velocity Grid Approach for BGK Equation". In: *Commun. Comput. Phys.* 16 (Feb. 2014).

[4] P. Bhatnagar, E. Gross, and M. Krook. "A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems". In: *Phys. Rev.* 94.3 (1954), pp. 511–525.

[5] M. Bisi, M. Groppi, and G. Spiga. “Kinetic Bhatnagar-Gross-Krook model for fast reactive mixtures and its hydrodynamic limit”. In: *Phys. Rev. E* 81 (3 Mar. 2010), p. 036327.
[6] M. Bisi, W. Boscheri, G. Dimarco, M. Groppi, and G. Martalò. "A new mixed Boltzmann-BGK model for mixtures solved with an IMEX finite volume scheme on unstructured meshes". In: *Appl. Math. Comput.* 433 (2022), p. 127416.

[7] M. Bisi, M. Groppi, and G. Martalò. "Macroscopic equations for inert gas mixtures in different hydrodynamic regimes". In: *J. Phys. Math. Theor.* 54.8 (Feb. 2021), p. 085201.

[8] A. Bobylev, M. Bisi, M. Groppi, G. Spiga, and I. Potapenko. "A general consistent BGK model for gas mixtures". In: *Kinet. Relat. Mod.* 11 (Dec. 2018), pp. 1377-1393.

[9] S. Boscarino, L. Pareschi, and G. Russo. *Implicit-Explicit Methods for Evolutionary Partial Differential Equations*. Philadelphia, PA: Society for Industrial and Applied Mathematics, 2024.

[10] S. Boscarino, S. Y. Cho, M. Groppi, and G. Russo. "BGK models for inert mixtures: Comparison and applications". In: *Kinet. Relat. Mod.* 14.5 (2021), pp. 895-928.

[11] S. Boscarino, S. Y. Cho, and G. Russo. "A local velocity grid conservative semi-Lagrangian schemes for BGK model". In: *J. Comput. Phys.* 460 (Apr. 2022), p. 111178.

[12] S. Boscarino, S.-Y. Cho, G. Russo, and S.-B. Yun. "High Order Conservative Semi-Lagrangian Scheme for the BGK Model of the Boltzmann Equation". In: *Commun. Comput. Phys.* 29 (June 2021), pp. 1-56.

[13] A. F. Bott, S. Cowley, and A. Schekochihin. "Kinetic stability of Chapman-Enskog plasmas". In: *J. Plasma Phys.* 90.2 (2024).

[14] M. Briant. "Stability of global equilibrium for the multi-species Boltzmann equation in $L^\infty$ settings". In: *Discrete Cont. Dyn-a.* 36.12 (2016), pp. 6669-6688.

[15] M. Briant and E. Daus. "The Boltzmann Equation for a Multi-species Mixture Close to Global Equilibrium". In: *Arch. Ration. Mech. An.* 222 (Dec. 2016).

[16] S. Brull. "An Ellipsoidal Statistical Model for gas mixtures". In: *Commun. Math. Sci.* 13 (Jan. 2015), pp. 1-13.

[17] S. Brull, V. Pavan, and J. Schneider. "Derivation of a BGK model for mixtures". In: *Eur. J. Mech. B. Fluids* 33 (2012), pp. 74-86.

[18] S. Brull and J. Schneider. "A new approach of the Ellipsoidal Statistical Model". In: *Continuum Mech. Therm.* 20 (May 2008), pp. 63-74.

[19] S. Brull and J. Schneider. "On the Ellipsoidal Statistical Model for polyatomic gases". In: *Continuum Mech. Therm.* 20 (Jan. 2009), pp. 489-508.

[20] D. Caparello, L. Pareschi, and T. Rey. *Hierarchical dynamic domain decomposition for the multiscale Boltzmann equation*. Preprint, arXiv:2505.03360 [math.NA] (2025). 2025.

[21] C. Cercignani. *The Boltzmann Equation and Its Applications*. Springer, 1987.

[22] C. Cercignani, R. Illner, and M. Pulvirenti. *The Mathematical Theory of Dilute Gases*. Springer, 2013.

[23] M. A. Diaz, M.-H. Chen, and J.-Y. Yang. "High-Order Conservative Asymptotic-Preserving Schemes for Modeling Rarefied Gas Dynamical Flows with Boltzmann-BGK Equation". In: *Commun. Comput. Phys.* 18.4 (2015), pp. 1012-1049.

[24] G. Dimarco and L. Pareschi. "Numerical methods for kinetic equations". In: *Acta Numer.* 23 (2014), pp. 369-520.

[25] F. Filbet and S. Jin. "A class of asymptotic-preserving schemes for kinetic equations and related problems with stiff sources". In: *J. Comput. Phys.* 229.20 (2010), pp. 7625-7648.

[26] F. Filbet and S. Jin. "An Asymptotic Preserving Scheme for the ES-BGK Model of the Boltzmann Equation". In: *J. Sci. Comput.* 46 (Feb. 2011), pp. 204-224.

[27] F. Filbet, C. Mouhot, and L. Pareschi. "Solving the Boltzmann Equation in $N \log 2N$". In: *SIAM J. Sci. Comput.* 28 (Aug. 2006), pp. 1029-1053.

[28] F. Filbet and T. Rey. "A Hierarchy of Hybrid Numerical Methods for Multiscale Kinetic Equations". In: *SIAM J. Sci. Comput.* 37 (May 2015), A1218-A1247.

[29] I. Gamba and S. Tharkabhushanam. "Spectral-Lagrangian methods for Collisional Models of Non-Equilibrium Statistical States". In: *J. Comput. Phys.* 228 (Apr. 2009), pp. 2012-2036.
[30] M. Groppi, S. Monica, and G. Spiga. "A kinetic ellipsoidal BGK model for a binary gas mixture". In: *Europhys. Lett.* 96 (Dec. 2011).

[31] J. Haack, C. Hauck, and M. Murillo. "A Conservative, Entropic Multispecies BGK Model". In: *J. Stat. Phys.* 168 (Aug. 2017).

[32] J. Holway Lowell H. "New Statistical Models for Kinetic Theory: Methods of Construction". In: *Phys. Fluids* 9.9 (Sept. 1966), pp. 1658-1673.

[33] J. Hu and L. Wang. "An asymptotic-preserving scheme for the semiconductor Boltzmann equation toward the energy-transport limit". In: *J. Comput. Phys.* 281 (2015), pp. 806-824.

[34] S. Jaiswal, A. A. Alexeenko, and J. Hu. "A discontinuous Galerkin fast spectral method for the multi-species Boltzmann equation". In: *Computer Methods in Applied Mechanics and Engineering* 352 (2019), pp. 56-84.

[35] S. Jin. "Asymptotic preserving (AP) schemes for multiscale kinetic and hyperbolic equations: a review". In: *Riv. Math. Univ. Parma (N.S.)* 2 (Jan. 2010).

[36] S. Jin. "Efficient Asymptotic-Preserving (AP) Schemes For Some Multiscale Kinetic Equations". In: *SIAM J. Sci. Comput.* 21.2 (1999), pp. 441-454.

[37] S. Jin and Q. Li. "A BGK-penalization-based asymptotic-preserving scheme for the multispecies Boltzmann equation". In: *Numer. Meth. Part. D. E.* 29 (May 2013).

[38] M. Junk. "Maximum entropy for reduced moment problems". In: *Math. Models Method Appl. Sci.* 10 (Oct. 2000), pp. 1001-1025.

[39] C. Levermore. "Moment closure hierarchies for kinetic theories". In: *J. Stat. Phys.* 83 (June 1996), pp. 1021-1065.

[40] D. Levy, G. Puppo, and G. Russo. "Compact Central WENO Schemes for Multidimensional Conservation Laws". In: *SIAM J. Sci. Comput.* 22 (1999), pp. 656-672.

[41] C. P. McNally, W. Lyra, and J.-C. Passy. "A well-posed Kelvin-Helmholtz instability test and comparison". In: *Astrophys. J. Suppl. Ser.* 201.2 (June 2012), p. 18.

[42] W. Melis, T. Rey, and G. Samaey. "Projective and telescopic projective integration for the nonlinear BGK and Boltzmann equations". In: *SMAI J. Comput. Math.* 5 (2019), pp. 53-88.

[43] L. Mieussens. "Discrete velocity model and implicit scheme for the BGK equation of rarefied gas dynamic". In: *Math. Models Method Appl. Sci.* 10 (Nov. 2000).

[44] L. Pareschi and G. Russo. "Implicit-explicit Runge-Kutta schemes for stiff systems of differential equations". In: *Recent Trends in Numerical Analysis*. USA: Nova Science Publishers, Inc., 2000, pp. 269-288.

[45] G. Puppo, M. Semplice, and G. Visconti. "Quinpi: Integrating Conservation Laws with CWENO Implicit Methods". In: *Commun. Appl. Math. Comput.* 5 (Dec. 2023), pp. 343-369.

[46] G. Puppo, T. Rey, and T. Tenna. *Formal derivation of an isentropic two-phase flow model from the multi-species Boltzmann equation*. 2025. arXiv: 2506.13480 [math.AP].

[47] T. Rey and T. Tenna. "The Boltzmann Equation for a Multi-Species Inelastic Mixture". In: *J. Stat. Phys.* 192.11 (2025), p. 151.

[48] J. Schneider. "Entropic approximation in kinetic theory". In: *ESAIM Math. Model. Numer. Anal.* 38.3 (2004), pp. 541-561.

[49] S. Tiwari. "Coupling of the Boltzmann and Euler Equations with Automatic Domain Decomposition". In: *J. Comput. Phys.* 144.2 (1998), pp. 710-726.

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

[51] L. Wu, J. Zhang, Y.-h. Zhang, and J. Reese. "A fast spectral method for the Boltzmann equation for monatomic gas mixtures". In: *J. Comput. Phys.* 298 (June 2015), pp. 602-621.

[52] J. Yang and J. Huang. "Rarefied Flow Computations Using Nonlinear Model Boltzmann Equations". In: *J. Comput. Phys.* 120.2 (1995), pp. 323-339.

[53] Y. Zhang, C. Zhang, X. Song, and Z. Guo. "Implicit Discrete Unified Gas Kinetic Scheme for Steady Flows of Binary Gas Mixtures". In: *Commun. Comput. Phys.* 37 (Jan. 2025), pp. 383-419.