[2]⟨#1,#2⟩ ()^+ ()^- [2]∥#2∥_L^#1 [3]∥#3∥_L^#1(#2) [2]∥#1∥_L^*#2^#1
[3]∥#2∥_L^*(#1)#3^#1 ref eqref Cref

# Projective and Telescopic Projective Integration for Non-Linear Kinetic
Mixtures

Rafael Bailo  Univ. Lille, CNRS, Inria, UMR 8524 - Laboratoire Paul Painlevé
Mathematical Institute, University of Oxford  Thomas Rey  Univ. Lille, CNRS,
Inria, UMR 8524 - Laboratoire Paul Painlevé

###### Abstract

We propose fully explicit projective integration and telescopic projective
integration schemes for the multispecies Boltzmann and Bhatnagar-Gross-Krook
(BGK) equations. The methods employ a sequence of small forward-Euler steps,
intercalated with large extrapolation steps. The telescopic approach repeats
said extrapolations as the basis for an even larger step. This hierarchy
renders the computational complexity of the method essentially independent of
the stiffness of the problem, which permits the efficient solution of
equations in the hyperbolic scaling with very small Knudsen numbers. We
validate the schemes on a range of scenarios, demonstrating its prowess in
dealing with extreme mass ratios, fluid instabilities, and other complex
phenomena.

AMS Subject Classification —  65M08; 65M12; 76P05; 82C40.

Keywords —  Multispecies gas; kinetic mixture; Boltzmann equation; BGK model;
projective integration; Sod Tube; Kelvin-Helmholtz; Richtmyer-Meshkov.

AP

    asymptotic preserving
BGK

    Bhatnagar-Gross-Krook
PInt

    projective integration
TPInt

    telescopic projective integration

##  1 Introduction

Mixtures of rarefied gases are found in a wide variety of systems, ranging
from the re-entry of an interplanetary probe in the upper atmosphere [21] to
microscale flows in pumps which use no moving parts (viz. the Knudsen
compressor, [37]). Such complex gases cannot be described by classical fluid
models, such as the compressible Euler or the Navier-Stokes systems, because
of their non-equilibrium behaviour, induced by their rarefaction. Kinetic
models, such as the seminal Boltzmann equation, are thus favoured to describe
these systems because they are able to reflect the non-equilibrium character
of the gases, retaining information about the microscopic many-particle
dynamics while avoiding the sheer complexity of the microscopic approach.
Furthermore, actual gases are usually a mixture of different chemical species;
for instance, the chemistry of the upper atmosphere is made of up to 20
different species (mostly recombinations of O2, CO2, H2, and CH4). The
realistic simulations of such systems must involve multispecies kinetic
models.

Whereas the mathematical properties of the classical Boltzmann equation for a
single species gas are well known by now (crucially, its derivation from
Newtonian dynamics was addressed in [28]), many questions remain open for the
multispecies case. The most recent theoretical results on the topic come from
the series of papers [9, 10]; they prove the existence, uniqueness,
positivity, and exponential trend to equilibrium for the full non-linear
multispecies Boltzmann equation in a perturbative, polynomially-weighted, and
isotropic $L^{1}_{v}L^{\infty}_{x}$ setting. The case of the unperturbed
setting remains mostly open.

At the numerical level, the most advanced deterministic methods use a Fourier
approach to evaluate the full collision operator: a fast spectral algorithm
was recently introduced in [61] for the multispecies case. Nevertheless, the
method was not asymptotic preserving (AP), namely, stable in the small
relaxation parameter limit; the only recent paper on this subject can be found
in [7], where the Maxwell-Stefan limit for a multispecies gas is investigated
numerically via a moment method (which does not compute the full operator). In
this work, we shall introduce a new family of numerical integrators for full
kinetic multispecies models that are able to deal with a large range of values
of the relaxation parameters uniformly on the numerical parameters.

Our work will rely on projective integration (PInt), a robust and fully
explicit method that allows for the time integration of (two-scale) stiff
systems with arbitrary order of accuracy in time. The method was first
proposed in [30] for stiff systems of ordinary differential equations with a
clear gap in their eigenvalue spectrum. In such problems, the fast modes,
corresponding to the Jacobian eigenvalues with large negative real parts,
decay quickly; it is the slow modes, related to eigenvalues of smaller
magnitude, that are of practical interest. In this regime, PInt permits a
stable yet explicit integration by combining small and large steps. The
integrator performs a few small (inner) steps of an explicit method, using a
step size $\delta t$ , until the transients corresponding to the fast modes
have decayed; subsequently, the solution is projected (extrapolated) forward
in time over a large (outer) time step of size ${\Delta t\gg\delta t}$.

PI was analysed for kinetic equations with a diffusive scaling in [41]. An
arbitrary order version, based on Runge-Kutta methods, has been proposed
recently in [39], where it was also analysed for kinetic equations in the
advection-diffusion limit. In [40], the scheme was used to construct a
explicit, flexible, arbitrary order method for general non-linear hyperbolic
conservation laws, based on their relaxation to a kinetic equation.
Alternative approaches to obtain higher-order PInt schemes have been proposed
in [42, 52]. These methods align with recent research efforts on numerical
methods for multiscale simulation [25, 29].

For problems exhibiting more than two time scales, telescopic projective
integration (TPInt) was proposed in [31]. In these methods, the projective
idea is applied recursively. Starting from an inner integrator at the fastest
time scale, a PInt method is constructed with a time step that corresponds to
the second-fastest scale. This PInt method is then considered as the inner
integrator of yet another PInt method at a coarser level. By repeating as
required, TPInt methods construct a hierarchy of projective levels, each using
the previous one as its inner step. This idea was explored for linear kinetic
equations in [49]. These methods turn out to have a computational cost that is
essentially independent of the stiffness of the collision operator. This
property was used with great success in the series of papers [47, 48] to
develop PInt and TPInt methods for the full non-linear Bhatnagar-Gross-Krook
(BGK) and Boltzmann equations of single-species rarefied gas dynamics.

PInt methods are not asymptotic preserving (AP) methods as such, because the
schemes cannot be evaluated explicitly at $\varepsilon=0$ to obtain a
classical numerical scheme for the limiting equation. Nevertheless, PInt and
TPInt methods share important features with AP methods. In particular, their
computational cost does (in many cases) not depend on the stiffness of the
problem. To be precise, it was shown in [49] for linear kinetic equations that
the number of inner time steps at each level of the telescopic hierarchy is
independent of the small-scale parameter $\varepsilon$, as is the step size of
the outermost integrator. The only parameter in the method which may depend on
$\varepsilon$ is the number of levels in the telescopic hierarchy. For systems
in which the spectrum of the collision operator falls apart into a set of
clearly separate clusters (each corresponding to a specific time scale), the
number of levels equals the number of spectral clusters. In this regime, the
computational cost is completely independent from $\varepsilon$. When the
collision operator comprises a continuum of time scales, the number of TPInt
levels increases logarithmically with $\varepsilon$.

The linearisation of the multispecies Boltzmann operator around its equilibria
has very similar properties to the classical single species Boltzmann
operator, as shown in [20, 10]. Among others, its spectrum is well separated,
with slow modes close to the origin, fast modes at a distance of order
$\varepsilon^{-1}$ left of the imaginary axis, and an essential part even
farther away from the axis [24, 51]; this information was the basis for the
development of the TPInt method in [48]. The spectral gap estimates in the
multispecies setting are also very similar to the single species case. In this
paper, we will exploit this structure to develop a TPInt method for the
multispecies Boltzmann equation.

The rest of the work is organised as follows. In Section 2 we recall the
elements of kinetic theory before presenting the models under consideration:
the multispecies Boltzmann and BGK equations. Section 3 develops the numerical
schemes which we shall use to study these models, with emphasis on the novelty
of this paper: the use of projective integration and telescopic projective
integration to construct uniformly accurate schemes for the multispecies BGK
model. In Section 4 we compute the approximate spectrum of the numerical
scheme and discuss the strategy for the choice of the parameters of the
numerical methods. Finally, Section 5 presents a variety of numerical
experiments, both in $1+1$ and $2+2$ dimensions of the phase space, which
demonstrate how the schemes can handle extreme mass ratios, fluid
instabilities, and other complex phenomena.

##  2 On Kinetic Equations and Multiple Species Models

The cornerstone of kinetic theory is the Boltzmann equation, which describes
the evolution in time of the distribution of the particles of a rarefied gas.
Each particle is subject to ballistic motion, travelling in a straight line at
a given velocity, which may only change when it collides with another
particle. These collisions are assumed to take place instantaneously and to
conserve momentum and kinetic energy. Moreover, the gas is assumed to be
dilute enough so that collisions between three or more particles can be
neglected.

The distribution of particles, $f^{\varepsilon}(t,\mathbf{x},\mathbf{v})$, is
a non-negative function which describes for every time $t$ the likelihood of
finding a particle in the infinitesimal configuration
$(\mathop{}\\!\mathrm{d}\mathbf{x},\mathop{}\\!\mathrm{d}\mathbf{v})$. In this
work, $f^{\varepsilon}$ is understood in the number sense, meaning it is not a
probability measure. To be precise111All the integrals in the text are taken
over the full position domain, the full velocity domain, or the entire domain
of Equation 2.4, unless otherwise stated:
$\displaystyle\int\mathop{}\\!\mathrm{d}\mathbf{x}\equiv\int_{\Omega}\mathop{}\\!\mathrm{d}\mathbf{x},\quad\int\mathop{}\\!\mathrm{d}\mathbf{v}\equiv\int_{{\mathbb{R}}^{{D_{\mathbf{v}}}}}\mathop{}\\!\mathrm{d}\mathbf{v},\quad\iint\mathop{}\\!\mathrm{d}\mathbf{x}\mathop{}\\!\mathrm{d}\mathbf{v}\equiv\iint_{\Omega\times{{\mathbb{R}}^{{D_{\mathbf{v}}}}}}\mathop{}\\!\mathrm{d}\mathbf{x}\mathop{}\\!\mathrm{d}\mathbf{v}.$
(2.1) , $\iint
f^{\varepsilon}(t,\mathbf{x},\mathbf{v})\mathop{}\\!\mathrm{d}\mathbf{x}\mathop{}\\!\mathrm{d}\mathbf{v}=N$
for all times $t\geq 0$, where $N$ is the number of particles in the gas. All
particles are assumed identical, each with mass $m>0$.

The distribution $f^{\varepsilon}$ evolves according to the kinetic equation

$\displaystyle\left\\{\begin{array}[]{{>{\displaystyle\vphantom{\frac{1}{1}}}l}}\partial_{t}f^{\varepsilon}+\mathbf{v}\cdot\nabla_{\mathbf{x}}f^{\varepsilon}=\frac{1}{\varepsilon}\mathcal{Q}[f^{\varepsilon}],\quad\mathbf{x}\in\Omega\subseteq{{\mathbb{R}}^{{D_{\mathbf{x}}}}},\,\mathbf{v}\in{{\mathbb{R}}^{{D_{\mathbf{v}}}}},\,t>0,\\\
f^{\varepsilon}(0,\mathbf{x},\mathbf{v})=f_{0}(\mathbf{x},\mathbf{v}),\end{array}\right.$
(2.4)

for some given Lipschitz domain $\Omega$, which has to be supplemented with
boundary conditions. The most simple physically relevant boundary conditions
are the so-called Maxwell boundary conditions, where a fraction
$1-\alpha\in[0,1]$ of the particles is specularly reflected by a wall, whereas
the remaining fraction $\alpha$ is thermalized, leaving the wall (moving at a
velocity $u_{w}$) at a temperature $T_{w}>0$. Such conditions can be written
as follows: given the unit outer normal ${n}_{\mathbf{x}}$ for all
$\mathbf{x}\in\partial\Omega$, for an incoming velocity $\mathbf{v}$ (namely
$\mathbf{v}\cdot{n}_{\mathbf{x}}\leq 0$), we set
$f(t,\mathbf{x},\mathbf{v})=\alpha\,\mathcal{R}f(t,\mathbf{x},\mathbf{v})\,+\,(1-\alpha)\,\mathcal{M}f(t,\mathbf{x},\mathbf{v}),$
with

$\left\\{\begin{aligned}
&\mathcal{R}f(t,\mathbf{x},\mathbf{v})\,\,=\,\,f(t,\mathbf{x},\mathbf{v}\,-\,2({n}_{\mathbf{x}}\cdot\mathbf{v})\,{n}_{\mathbf{x}}),\\\
&\mathcal{M}f(t,\mathbf{x},\mathbf{v})\,\,=\,\,\frac{\mu(t,\mathbf{x})}{(2\pi\,T_{w})^{\mathop{}\\!\mathrm{d}\mathbf{v}/2}}\exp\left(-\frac{|\mathbf{v}-u_{w}|^{2}}{2\,T_{w}}\right).\end{aligned}\right.$

In this last expression, $\mu$ insures global mass conservation. Nevertheless,
for the sake of simplicity, we shall only consider either periodic or influx
boundary conditions. This equation couples a transport operator, which models
the displacement of particles, with a collision operator $\mathcal{Q}$, which
describes the changes in momentum due to collisions. In this work,
${D_{\mathbf{x}}}={D_{\mathbf{v}}}=1$ or $2$.

The frequency of collisions is governed by the dimensionless Knudsen number
$\varepsilon$, defined as the ratio between the mean free path of particles
and the length scale of observation; broadly, a measure of the average
distance that a particle can travel before colliding with another one. The
value of $\varepsilon$ distinguishes the kinetic regime, where collisions are
rare, from the hydrodynamic regime, where collisions dominate the dynamics;
the former corresponds to larger values of $\varepsilon$, and the latter, to
smaller values.

The seminal Boltzmann collision operator $\mathcal{Q}$ is a quadratic operator
local in $(t,\mathbf{x})$. The time and position act only as parameters in
$\mathcal{Q}$ and therefore will be omitted in its description. It is given,
for ${D_{\mathbf{v}}}\geq 2$, by

$\mathcal{Q}[f](\mathbf{v})=\iint_{{{\mathbb{R}}^{{D_{\mathbf{v}}}}}\times\mathbb{S}^{{D_{\mathbf{v}}}-1}}B(|\mathbf{v}-\mathbf{v}_{*}|,\cos\theta)\,\left(f^{\prime}f^{\prime}_{*}-ff_{*}\right)\mathop{}\\!\mathrm{d}\mathbf{v}_{*}\mathop{}\\!\mathrm{d}\sigma,$
(2.5)

where we have used the shorthand $f=f(\mathbf{v})$, $f_{*}=f(\mathbf{v}_{*})$,
$f^{{}^{\prime}}=f(\mathbf{v}^{\prime})$,
$f_{*}^{{}^{\prime}}=f(\mathbf{v}_{*}^{{}^{\prime}})$. The velocities of the
colliding pairs, $(\mathbf{v},\mathbf{v}_{*})$ and
$(\mathbf{v}^{\prime},\mathbf{v}^{\prime}_{*})$, can be parametrized as

$\mathbf{v}^{\prime}=\frac{\mathbf{v}+\mathbf{v}_{*}}{2}+\frac{|\mathbf{v}-\mathbf{v}_{*}|}{2}\sigma,\qquad\mathbf{v}^{\prime}_{*}=\frac{\mathbf{v}+\mathbf{v}^{*}}{2}-\frac{|\mathbf{v}-\mathbf{v}_{*}|}{2}\sigma.$
(2.6)

The collision kernel $B$ is a non-negative function which, by physical
arguments of invariance, may only depend on $|\mathbf{v}-\mathbf{v}_{*}|$ and
$\cos\theta={\hat{g}}\cdot\sigma$ (where
${\hat{g}}=(\mathbf{v}-\mathbf{v}_{*})/|\mathbf{v}-\mathbf{v}_{*}|$). It
characterises the details of the binary interactions, and has the form

$B(|\mathbf{v}-\mathbf{v}_{*}|,\cos\theta)=|\mathbf{v}-\mathbf{v}_{*}|\,\Phi(|\mathbf{v}-\mathbf{v}_{*}|,\cos\theta).$
(2.7)

The scattering cross-section $\Phi$, in the case of inverse $k$th power forces
between particles, can be written as

$\Phi(|\mathbf{v}-\mathbf{v}_{\ast}|,\cos\theta)=b_{\alpha}(\cos\theta)\,|\mathbf{v}-\mathbf{v}_{\ast}|^{\alpha-1},$

with $\alpha=(k-5)/(k-1)$. The special situation $k=5$ gives the so-called
Maxwell pseudo-molecules model with
$B(|\mathbf{v}-\mathbf{v}_{*}|,\cos\theta)=b_{0}(\cos\theta).$ For the Maxwell
case the collision kernel is independent of the relative velocity. For
numerical purposes, a widely used model is the variable hard sphere (VHS)
model introduced by Bird [5]. The model corresponds to
$b_{\alpha}(\cos\theta)=C_{\alpha}$, where $C_{\alpha}$ is a positive
constant, and hence
$\Phi(|\mathbf{v}-\mathbf{v}_{\ast}|,\cos\theta)=C_{\alpha}|\mathbf{v}-\mathbf{v}_{\ast}|^{\alpha-1}.$
For further details on the physical background and derivation of the Boltzmann
equation, we refer to [16, 59]

Equation 2.4 is rich with structural properties. In terms of conservation, its
weak form reveals that any conserved quantity of the dynamics corresponds to
an associated collision invariant of $\mathcal{Q}$: a function
$\varphi(\mathbf{v})$ such that

$\int\mathcal{Q}[f](\mathbf{v})\,\varphi(\mathbf{v})\mathop{}\\!\mathrm{d}\mathbf{v}=0$
(2.8)

for any function $f(\mathbf{v})$. For the Boltzmann operator, the space of
collision invariants is given by
$\operatorname{span}\\{1,\mathbf{v},|\mathbf{v}|^{2}\\}$ [15]. In consequence,
the mass, momentum, and kinetic energy of the distribution, respectively
$\iint
mf^{\varepsilon}(t,\mathbf{x},\mathbf{v})\mathop{}\\!\mathrm{d}\mathbf{x}\mathop{}\\!\mathrm{d}\mathbf{v}$,
$\iint
m\mathbf{v}f^{\varepsilon}(t,\mathbf{x},\mathbf{v})\mathop{}\\!\mathrm{d}\mathbf{x}\mathop{}\\!\mathrm{d}\mathbf{v}$,
and
$\iint\frac{m}{2}|\mathbf{v}|^{2}f^{\varepsilon}(t,\mathbf{x},\mathbf{v})\mathop{}\\!\mathrm{d}\mathbf{x}\mathop{}\\!\mathrm{d}\mathbf{v}$,
are conserved.

In terms of dynamics, the relaxation properties of Eq. 2.4 are well-
understood. The celebrated H-theorem characterises the dissipation of the
Boltzmann entropy222 As is usual, we choose units such that Boltzmann’s
constant equals one, $k_{B}=1$, throughout the text. ,

$\displaystyle\mathcal{H}=\int
f^{\varepsilon}(\mathbf{v})\log(f^{\varepsilon}(\mathbf{v}))\mathop{}\\!\mathrm{d}\mathbf{v}.$
(2.9)

At the local level, the theorem states

$\displaystyle\mathcal{S}=\int\mathcal{Q}[f^{\varepsilon}](\mathbf{v})\log(f^{\varepsilon}(\mathbf{v}))\mathop{}\\!\mathrm{d}\mathbf{v}\leq
0,$ (2.10)

where $\frac{\partial\mathcal{H}}{\partial
t}+\nabla_{\mathbf{x}}\cdot\mathcal{J}=\frac{\mathcal{S}}{\varepsilon}$ and
$\mathcal{J}=\int\mathbf{v}f^{\varepsilon}(\mathbf{v})\log(f^{\varepsilon}(\mathbf{v}))\mathop{}\\!\mathrm{d}\mathbf{v}$,
see, for instance, [14, 16]. Furthermore, the equality in (2.10) is only
achieved for states in the kernel of $\mathcal{Q}$. Such equilibrium states
are always Maxwellian distributions:

$\displaystyle
f^{\varepsilon}=\mathcal{M}^{n,\bar{\mathbf{v}},T}(t,\mathbf{x},\mathbf{v})\coloneqq
n(t,\mathbf{x})\,\left(\frac{m}{2\pi
T}\right)^{{D_{\mathbf{v}}}/2}\exp\left(-\frac{m(\bar{\mathbf{v}}-\mathbf{v})^{2}}{2T}\right)$
(2.11)

for particles of mass $m$. The number density $n$, the mass density $\rho$,
the average velocity $\bar{\mathbf{v}}$, the temperature $T$, and the pressure
$\mathcal{P}$ are (local) moments computed from the distribution which depend
only on time and position; they are given by

$\displaystyle n=\int
f^{\varepsilon}\mathop{}\\!\mathrm{d}\mathbf{v},\quad\rho=mn,\quad\rho\bar{\mathbf{v}}=\int
m\mathbf{v}f^{\varepsilon}\mathop{}\\!\mathrm{d}\mathbf{v},\quad\frac{{D_{\mathbf{v}}}}{2}nT=\int\frac{m}{2}|\bar{\mathbf{v}}-\mathbf{v}|^{2}f^{\varepsilon}\mathop{}\\!\mathrm{d}\mathbf{v},$
(2.12)

where $\rho=mn$ and $\mathcal{P}=nT$. Note, in particular, that the Maxwellian
$\mathcal{M}^{n,\bar{\mathbf{v}},T}$ has moments $n,\bar{\mathbf{v}},T$.

###  2.1 The BGK Model

The collision operator (2.5) has a complicated structure from both the
analytical and the numerical perspectives. Because of this, many works have
proposed simpler operators which attempt to capture some or all of the
structural properties of $\mathcal{Q}$. A ubiquitous simplification is the
model of Bhatnagar-Gross-Krook [4]; since the overall effect of the Boltzmann
operator is to drive $f^{\varepsilon}$ towards the corresponding Maxwellian,
they propose a reduced operator which makes the relaxation explicit:

$\displaystyle\mathcal{Q}^{\text{BGK}}[f^{\varepsilon}]=\nu(\mathcal{M}[f^{\varepsilon}]-f^{\varepsilon}),$
(2.13)

where $\mathcal{M}[f^{\varepsilon}]$ is the Maxwellian whose moments
$n,\bar{\mathbf{v}},T$ are the ones of $f^{\varepsilon}$, and where
$\nu(t,\mathbf{x})$ is a positive collision rate to be determined. This
simplification can be derived from $\mathcal{Q}$ by assuming that the
distribution $f^{\varepsilon}$ is already close to equilibrium. The BGK
equation is thus defined as the evolution law (2.4) together with the operator
(2.13).

Despite its apparent simplicity, the BGK model captures many structural
aspects of the Boltzmann collision operator. Because
$\mathcal{M}[f^{\varepsilon}]$ and $f^{\varepsilon}$ share the first three
moments, the invariants of $\mathcal{Q}$ are also invariants of
$\mathcal{Q}^{\text{BGK}}$ whenever the collision rate does not depend on
$\mathbf{v}$. Therefore, the new dynamics conserve the mass, momentum, and
kinetic energy of $f^{\varepsilon}$ as well. The dynamical structure of
$\mathcal{Q}$ is also preserved by the BGK model because the H-theorem still
holds:

$\displaystyle\mathcal{S}=\int\mathcal{Q}^{\text{BGK}}[f^{\varepsilon}](\mathbf{v})\log(f^{\varepsilon}(\mathbf{v}))\mathop{}\\!\mathrm{d}\mathbf{v}\leq
0,$ (2.14)

see, for instance, [56]. Once again, the equality is only achieved by the
equilibrium states, which are the Maxwellians by construction.

###  2.2 The Multispecies Boltzmann Equation

A limitation of the Boltzmann equation (2.4) is the assumption that all gas
particles are identical. However, it is possible to extend this equation to
the case of a mixture of gases by the physical arguments used in the single
species case [8]. As stated in the introduction, the derivation of the
multispecies Boltzmann equation is mostly formal, unlike that of the single
species case, which has been clearly established in a variety of works such as
[28].

We will consider a mixture of $P$ species, each described by a
distribution333The $\varepsilon$ notation is dropped in the interest of
simplicity. Henceforth, $f$ stands for $f^{\varepsilon}$.
$f_{p}(t,\mathbf{x},\mathbf{v})$, a time-dependent non-negative function as
before. Each distribution is understood in the number sense:

$\displaystyle\iint
f_{p}(t,\mathbf{x},\mathbf{v})\mathop{}\\!\mathrm{d}\mathbf{x}\mathop{}\\!\mathrm{d}\mathbf{v}=N_{p}$
(2.15)

for all times $t\geq 0$, where $N_{p}$ is the number of particles in the
$p\textsuperscript{th}$ species. All particles of that species are assumed
identical, each with mass $m_{p}$.

The distributions $f_{p}$ evolve according to the so-called multispecies
Boltzmann equation, given by

$\displaystyle\left\\{\begin{array}[]{{>{\displaystyle\vphantom{\frac{1}{1}}}l}}\partial_{t}f_{p}+\mathbf{v}\cdot\nabla_{\mathbf{x}}f_{p}=\frac{1}{\varepsilon}\mathcal{Q}_{p}[\boldsymbol{f}],\quad\mathbf{x}\in\Omega\subseteq{{\mathbb{R}}^{{D_{\mathbf{x}}}}},\,\mathbf{v}\in{{\mathbb{R}}^{{D_{\mathbf{v}}}}},\,t>0,\\\
\mathcal{Q}_{p}[\boldsymbol{f}]=\sum_{q=1}^{P}\mathcal{Q}_{p,\,q}[f_{p},f_{q}],\\\
f_{p}(0,\mathbf{x},\mathbf{v})=f_{p,\,0}(\mathbf{x},\mathbf{v}),\end{array}\right.$
(2.19)

posed with appropriate boundary conditions. The changes of momentum of the
$p\textsuperscript{th}$ species are now governed by a sum of collision
operators, one for each species in the gas. These multispecies collision
operators $\mathcal{Q}_{p,\,q}$ are similar to the classical Boltzmann
operator (2.5); namely, they are given by

$\mathcal{Q}_{p,\,q}[f_{p},f_{q}]=\iint_{{{\mathbb{R}}^{{D_{\mathbf{v}}}}}\times\mathbb{S}^{{D_{\mathbf{v}}}-1}}B_{p,\,q}(|\mathbf{v}-\mathbf{v}_{*}|,\cos\theta)\left[f_{p}\left(\mathbf{v}^{\prime}\right)f_{q}\left(\mathbf{v}^{\prime}_{*}\right)-f_{p}(\mathbf{v})\,f_{q}(\mathbf{v}_{*})\right]\mathop{}\\!\mathrm{d}\mathbf{v}_{*}\mathop{}\\!\mathrm{d}\sigma,$
(2.20)

where inter-species collisions are given by

$\displaystyle\begin{split}&\mathbf{v}^{\prime}\,=\,\frac{1}{m_{p}+m_{q}}\left(m_{p}\mathbf{v}+m_{q}\mathbf{v}_{*}+m_{q}|\mathbf{v}-\mathbf{v}_{*}|\sigma\right),\\\
&\mathbf{v}^{\prime}_{*}\,=\,\frac{1}{m_{p}+m_{q}}\left(m_{p}\mathbf{v}+m_{q}\mathbf{v}_{*}-m_{p}|\mathbf{v}-\mathbf{v}_{*}|\sigma\right).\end{split}$
(2.21)

Note that in the intra-species collision case $p=q$, (2.20) corresponds to the
classical Boltzmann operator (2.5). We shall also assume for the sake of
simplicity that the inter-species collision kernel $B_{p,\,q}$ is independent
on $p$ and $q$, and equal to the variable hard spheres kernel (2.7).

The conservation properties of the Boltzmann equation persist in the
multispecies case at the level of the mixture. The mass of each species,
$\iint
m_{p}f_{p}(t,\mathbf{x},\mathbf{v})\mathop{}\\!\mathrm{d}\mathbf{x}\mathop{}\\!\mathrm{d}\mathbf{v},$
is conserved because $1$ remains an invariant of each operator:

$\displaystyle\int\mathcal{Q}_{p,\,q}[f,g](\mathbf{v})\mathop{}\\!\mathrm{d}\mathbf{v}=0$
(2.22)

for any functions $f(\mathbf{v})$ and $g(\mathbf{v})$. However, the functions
$\mathbf{v}$ and $|\mathbf{v}|^{2}$ are only invariants in the sense that

$\displaystyle\int\mathbf{v}\mathcal{Q}_{p,\,q}[f,g](\mathbf{v})\mathop{}\\!\mathrm{d}\mathbf{v}=-\int\mathbf{v}\mathcal{Q}_{q,\,p}[g,f](\mathbf{v})\mathop{}\\!\mathrm{d}\mathbf{v}$
(2.23)

and

$\displaystyle\int|\mathbf{v}|^{2}\mathcal{Q}_{p,\,q}[f,g](\mathbf{v})\mathop{}\\!\mathrm{d}\mathbf{v}=-\int|\mathbf{v}|^{2}\mathcal{Q}_{q,\,p}[g,f](\mathbf{v})\mathop{}\\!\mathrm{d}\mathbf{v}.$
(2.24)

As such, it is only the total momentum and the total kinetic energy of the
gas, respectively

$\displaystyle\sum\limits_{p=1}^{P}\iint
m_{p}\mathbf{v}f_{p}(t,\mathbf{x},\mathbf{v})\mathop{}\\!\mathrm{d}\mathbf{x}\mathop{}\\!\mathrm{d}\mathbf{v},\text{
and
}\,\sum\limits_{p=1}^{P}\iint\frac{m_{p}}{2}|\mathbf{v}|^{2}f_{p}(t,\mathbf{x},\mathbf{v})\mathop{}\\!\mathrm{d}\mathbf{x}\mathop{}\\!\mathrm{d}\mathbf{v},$
(2.25)

that are conserved.

The H-theorem also applies to the multispecies case, for the total Boltzmann
entropy and total dissipation:

$\displaystyle\mathcal{H}=\sum\limits_{p=1}^{P}\int
f_{p}(\mathbf{v})\log(f_{p}(\mathbf{v}))\mathop{}\\!\mathrm{d}\mathbf{v},\quad\mathcal{S}=\sum\limits_{p=1}^{P}\int\mathcal{Q}_{p}[f_{p}](\mathbf{v})\log(f_{p}(\mathbf{v}))\mathop{}\\!\mathrm{d}\mathbf{v}\leq
0,$ (2.26)

see [17]. The equilibrium distributions are all Maxwellians with common
average velocity, $\bar{\mathbf{v}}_{\text{eq}}$, and temperature,
$T_{\text{eq}}$:

$\displaystyle
f_{p}=\mathcal{M}^{n,\bar{\mathbf{v}}_{\text{eq}},T_{\text{eq}}}(t,\mathbf{x},\mathbf{v})=n_{p}(t,\mathbf{x})\,\left(\frac{m_{p}}{2\pi
T_{\text{eq}}}\right)^{{D_{\mathbf{v}}}/2}\exp\left(-\frac{m_{p}(\bar{\mathbf{v}}_{\text{eq}}-\mathbf{v})^{2}}{2T_{\text{eq}}}\right).$
(2.27)

Note $n_{p}$ and the rest of individual moments are defined just as in Eq.
2.12:

$\displaystyle\begin{split}n_{p}=\int
f_{p}\mathop{}\\!\mathrm{d}\mathbf{v},\quad\rho_{p}\bar{\mathbf{v}}_{p}=\int
m_{p}\mathbf{v}f_{p}\mathop{}\\!\mathrm{d}\mathbf{v},\quad\frac{{D_{\mathbf{v}}}}{2}n_{p}T_{p}=\int\frac{m_{p}}{2}(\bar{\mathbf{v}}_{p}-\mathbf{v})^{2}f_{p}\mathop{}\\!\mathrm{d}\mathbf{v},\end{split}$
(2.28)

where $\rho_{p}=m_{p}n_{p}$ and $\mathcal{P}_{p}=n_{p}T_{p}$. The total
moments for the mixture are given by $n=\sum n_{p}$, $\rho=\sum\rho_{p}$,
$\rho\bar{\mathbf{v}}=\sum\rho_{p}\bar{\mathbf{v}}_{p}$,
$\mathcal{P}=\sum\mathcal{P}_{p}$, and

$\displaystyle\begin{split}\frac{{D_{\mathbf{v}}}}{2}nT=\sum\limits_{p=1}^{P}\int\frac{m_{p}}{2}(\bar{\mathbf{v}}-\mathbf{v})^{2}f_{p}\mathop{}\\!\mathrm{d}\mathbf{v}.\end{split}$
(2.29)

###  2.3 A Multispecies BGK Model

Although the multispecies Boltzmann equation described in Section 2.2 shares
many features with the single species Boltzmann equation, the fact that
conservations are only global in space and all species, e.g. (2.24), makes its
simplification to a relaxation operator such as the BGK model (2.13) much
harder. Many different approaches exists to derive such a multispecies BGK
operator, but very few are satisfactory in regards to the macroscopic
properties of the resulting relaxation operator. The interested reader can
find more details about these various approaches in the recent papers [3, 11].

This work will follow along the lines of [33]. The formulation of the
multispecies BGK operator requires the definition of the mixture Maxwellians:

$\displaystyle\mathcal{M}_{p,\,q}[f_{p},f_{q}](t,\mathbf{x},\mathbf{v})=n_{p}(t,\mathbf{x})\left(\frac{m_{p}}{2\pi
T_{p,\,q}}\right)^{{{D_{\mathbf{v}}}}/{2}}\exp\left(-\frac{m_{p}(\bar{\mathbf{v}}_{p,\,q}-\mathbf{v})^{2}}{2T_{p,\,q}}\right),$
(2.30)

where $\bar{\mathbf{v}}_{p,\,q}$ and $T_{p,\,q}$ are mixture moments to be
determined. It is assumed that the interaction of the $p\textsuperscript{th}$
and the $q\textsuperscript{th}$ species drives $f_{p}$ towards
$\mathcal{M}_{p,\,q}$, as well as $f_{q}$ towards $\mathcal{M}_{q,\,p}$; thus
we require $\bar{\mathbf{v}}_{p,\,q}=\bar{\mathbf{v}}_{q,\,p}$,
$T_{p,\,q}=T_{q,\,p}$. Assuming the distributions are close to equilibrium
permits the simplification of $\mathcal{Q}_{p,\,q}$ to obtain

$\displaystyle\mathcal{Q}^{\text{BGK}}_{p,\,q}[f_{p},f_{q}]=\nu_{p,\,q}(\mathcal{M}_{p,\,q}[f_{p},f_{q}]-f_{p}),$
(2.31)

where $\nu_{p,\,q}$ are positive collision rates444 From the derivation of the
multispecies BGK model,
$\nu_{p,\,q}(\mathbf{v})=\iint_{{{\mathbb{R}}^{{D_{\mathbf{v}}}}}\times\mathbb{S}^{{D_{\mathbf{v}}}-1}}B_{p,\,q}(|\mathbf{v}-\mathbf{v}_{*}|,\cos\theta)\mathcal{M}_{q,\,p}(\mathbf{v})\mathop{}\\!\mathrm{d}\mathbf{v}_{*}\mathop{}\\!\mathrm{d}\sigma,$
(2.32) see [33]. However, as done in the reference, we assume the collision
frequencies are given constants, $\nu_{p,\,q}\equiv 1$. . The multispecies BGK
model is thus defined as the evolution law (2.19) together with the operator
(2.31).

In order to recover the structural properties of the Boltzmann equation, the
mixture moments have to be chosen appropriately. Considering the interaction
of the $p\textsuperscript{th}$ species with itself, we deduce
$\bar{\mathbf{v}}_{p,\,p}=\bar{\mathbf{v}}_{p}$ and $T_{p,\,p}=T_{p}$. To
define $\bar{\mathbf{v}}_{p,\,q}$, it is enough to impose the conservation of
momentum in the form of Eq. 2.23. This readily yields

$\displaystyle\bar{\mathbf{v}}_{p,\,q}=\frac{\rho_{p}\nu_{p,\,q}\bar{\mathbf{v}}_{p}+\rho_{q}\nu_{q,\,p}\bar{\mathbf{v}}_{q}}{\rho_{p}\nu_{p,\,q}+\rho_{q}\nu_{q,\,p}}.$
(2.33)

This expression is reminiscent to the microscopic parametrisation of the
postcollisional velocities as a function of the precollisional ones (2.21).
Similarly, requiring Eq. 2.24 as an expression of conservation of kinetic
energy produces

$\displaystyle
T_{p,\,q}=\frac{n_{p}\nu_{p,\,q}T_{p}+n_{q}\nu_{q,\,p}T_{q}}{n_{p}\nu_{p,\,q}+n_{q}\nu_{q,\,p}}+\frac{\rho_{p}\nu_{p,\,q}(\bar{\mathbf{v}}_{p}^{2}-\bar{\mathbf{v}}_{p,\,q}^{2})+\rho_{q}\nu_{q,\,p}(\bar{\mathbf{v}}_{q}^{2}-\bar{\mathbf{v}}_{q,\,p}^{2})}{{D_{\mathbf{v}}}(n_{p}\nu_{p,\,q}+n_{q}\nu_{q,\,p})}.$
(2.34)

This expression guarantees $T_{p,\,q}\geq 0$, see [33, Appendix 9]. These
choices ensure, by construction, the conservation of total momentum and total
kinetic energy of the gas.

This model also captures the dissipative structure of the multispecies
Boltzmann equations. It is shown in [33] that the collision operator (2.31)
verifies the H-theorem:

$\displaystyle\mathcal{S}=\sum\limits_{p=1}^{P}\int\mathcal{Q}^{\text{BGK}}_{p}[f_{p}](\mathbf{v})\log(f_{p}(\mathbf{v}))\mathop{}\\!\mathrm{d}\mathbf{v}\leq
0.$ (2.35)

###  2.4 Hydrodynamic Limits

We will briefly present here the zeroth order hydrodynamic limit of the
multispecies BGK model (2.31), following along the lines of the recent paper
[6]. Taking the limit $\varepsilon\to 0$ in that kinetic equation will
formally project the distributions $f_{p}$ towards the equilibrium
distributions (2.27) whose moments solve a system of conservation laws. This
system is obtained by multiplying the kinetic multispecies equation (2.19) by
the powers of $\mathbf{v}$, integrating in velocity, and summing the equations
on the individual moments in order to have only ${D_{\mathbf{v}}}$ equations
for the total momentum.

Under the Maxwellian closure, we obtain the following system of
${D_{\mathbf{v}}}+P+1$ conservation law, the so-called multispecies Euler
limit of (2.19):

$\left\\{\begin{array}[]{{>{\displaystyle\vphantom{\frac{1}{1}}}l}}\partial_{t}n_{p}+\nabla_{\mathbf{x}}\cdot\left(n_{p}\bar{\mathbf{v}}\right)=0,\quad
p\in\\{1,\cdots,P\\},\\\\[5.0pt]
\partial_{t}\left(\rho\,\bar{\mathbf{v}}\right)+\nabla_{\mathbf{x}}\cdot\left(\rho\,\bar{\mathbf{v}}\otimes\bar{\mathbf{v}}\right)+\nabla_{\mathbf{x}}\left(n\,T\right)=\mathbf{0},\\\\[5.0pt]
\partial_{t}\left(\frac{1}{2}\rho\,|\bar{\mathbf{v}}|^{2}+\frac{{D_{\mathbf{v}}}}{2}n\,T\right)+\nabla_{\mathbf{x}}\cdot\left[\left(\frac{1}{2}\rho\,|\bar{\mathbf{v}}|^{2}+\frac{{D_{\mathbf{v}}}+2}{2}n\,T\right)\bar{\mathbf{v}}\right]=0.\end{array}\right.$
(2.36)

The interested reader can consult [6] for the terms corresponding to the next
order in the expansion, the multispecies Navier-Stokes limit of (2.19).

##  3 Numerical Schemes

This section is devoted to the development of efficient numerical methods for
the multispecies BGK model, given by Eqs. 2.19 and 2.31. The numerical
simulation of the model is hindered by the presence of two radically different
scales. The transport term classically imposes a CFL condition on the mesh
size: $\Delta t=\mathcal{O}(\Delta x)$ is needed for stability. However, in
the regimes where $\varepsilon\ll 1$, the collision term dominates the
dynamics, further restricting the time step to $\Delta
t=\mathcal{O}(\varepsilon)$. This requirement is particularly pernicious in
view of the large cost per step (due to the computation of non-local terms),
and it would render a naïve scheme inefficient. The development of schemes for
kinetic and hyperbolic equations which can handle multiple scales without such
penalties (sometimes called asymptotic-preserving schemes) is a very active
area of research, see the recent reviews [36, 23, 34].

In order to construct an efficient numerical scheme for Equation 2.19, we
resort to projective and telescopic projective integration. First introduced
in [30], projective integration provides an efficient framework for the
numerical solution of differential equations with two distinct timescales,
characterised by two separate clusters of eigenvalues. Telescopic projective
integration was introduced by [31] in order to extend the method to problems
with many scales, or without a clear spectral separation. These methods are
not asymptotic preserving in the classical sense [35] because they cannot be
evaluated at $\varepsilon=0$ to recover a scheme for the asymptotic equation;
nevertheless, they are able to solve stiff problems efficiently, and have been
successfully adapted to some kinetic equations: radiative transfer [41], the
linearised BGK equation with multiple relaxation times [49], and both the full
BGK model and the Boltzmann equation [47, 48]. Higher-order approaches have
also been developed [39, 40].

In the remainder of this section we describe the numerical method. First we
introduce the phase space discretisation, using a discrete velocity method for
the collision operator, coupled with a finite-volume scheme for the transport.
Then we present the projective and telescopic schemes.

###  3.1 Phase Space Discretisation

##### Velocity

The discretisation of Equation 2.19 can be conducted in a number of ways; see
[23] for a survey. Here we shall employ a discrete velocity approach as in
[53, 12, 50], where the velocity space is approximated by a finite grid. The
choice of grid is non-trivial, and it can affect the entropic properties of
the numerical scheme [1]; for simplicity, we choose a Cartesian grid.

The velocity space, ${\mathbb{R}}^{D_{\mathbf{v}}}$, is first restricted to a
bounded domain $(-L_{v},L_{v})^{D_{\mathbf{v}}}$. Each dimension is
discretised into $2N_{v}+1$ points, separated by a distance $\Delta
v=2L_{v}/(2N_{v}+1)$; the $j\textsuperscript{th}$ point is given by
$v_{j}=j\Delta v$, where
$j\in\mathcal{N}_{v}\coloneqq\\{-N_{v},\cdots,N_{v}\\}$. The discretised
velocity space becomes
$\mathcal{V}=\\{\mathbf{v}_{\boldsymbol{j}}=(v_{j_{1}},\cdots,v_{j_{D_{\mathbf{v}}}})\\}$,
where
${\boldsymbol{j}}=(j_{1},\cdots,j_{D_{\mathbf{v}}})\in\mathcal{N}_{v}^{D_{\mathbf{v}}}$
is a multi-index. Equation 2.19 thus becomes a system of transport equations:

$\displaystyle\left\\{\begin{array}[]{{>{\displaystyle\vphantom{\frac{1}{1}}}l}}\partial_{t}f_{p,\,{\boldsymbol{j}}}+\mathbf{v}_{{\boldsymbol{j}}}\cdot\nabla_{\mathbf{x}}f_{p,\,{\boldsymbol{j}}}=\frac{1}{\varepsilon}\mathcal{Q}_{p,\,{\boldsymbol{j}}}[f_{p}],\quad\mathbf{x}\in\Omega\subseteq{{\mathbb{R}}^{{D_{\mathbf{x}}}}},\,{\boldsymbol{j}}\in\mathcal{N}_{v}^{D_{\mathbf{v}}},\,t>0,\\\
\mathcal{Q}_{p,\,{\boldsymbol{j}}}[f_{p}]=\sum_{q=1}^{P}\mathcal{Q}_{p,\,q,\,{\boldsymbol{j}}}[f_{p},f_{q}],\end{array}\right.$
(3.3)

where $f_{p,\,{\boldsymbol{j}}}$ approximates
$f_{p}(\mathbf{v}_{\boldsymbol{j}})$, and
$\mathcal{Q}_{p,\,q,\,{\boldsymbol{j}}}[f_{p},f_{q}]$ approximates
$\mathcal{Q}_{p,\,q}[f_{p},f_{q}](\mathbf{v}_{\boldsymbol{j}})$. The initial
datum $f_{p,\,{\boldsymbol{j}}}(0,\mathbf{x})$ is prescribed by evaluating the
continuous datum, $f_{p,\,0}$, at $\mathbf{v}_{\boldsymbol{j}}$. Note that the
evolution of $f_{p,\,{\boldsymbol{j}}}$ depends on
$\mathcal{Q}_{p,\,{\boldsymbol{j}}}$, which itself depends on the entire
distribution of every species, coupling all the equations in the system.

The evaluation of the collision operator will require the approximation of the
moments of the distribution based on the discretised velocity grid. The semi-
discrete individual moments are555 As with the integrals, the sums in this
text are taken over the full velocity domain, unless otherwise stated:
$\displaystyle\sum\Delta\mathbf{v}\equiv\sum_{{\boldsymbol{j}}\in\mathcal{N}_{v}^{D_{\mathbf{v}}}}\Delta\mathbf{v},$
(3.6) where $\mathcal{N}_{v}\coloneqq\\{-N_{v},\cdots,N_{v}\\}$, see the
discussion around Eq. 3.3.

$\displaystyle n^{h}_{p}=\sum
f_{p,\,{\boldsymbol{j}}}\Delta\mathbf{v},\quad\rho^{h}_{p}\bar{\mathbf{v}}^{h}_{p}=\sum
m_{p}\mathbf{v}_{\boldsymbol{j}}f_{p,\,{\boldsymbol{j}}}\Delta
v,\quad\frac{{D_{\mathbf{v}}}}{2}n^{h}_{p}T^{h}_{p}=\sum\frac{m_{p}}{2}(\bar{\mathbf{v}}^{h}_{p}-\mathbf{v}_{\boldsymbol{j}})^{2}f_{p,\,{\boldsymbol{j}}}\Delta
v,$ (3.4)

where $\Delta\mathbf{v}=\Delta v^{D_{\mathbf{v}}}$, and
$\rho^{h}_{p}=m_{p}n^{h}_{p}$ as before. The approximate mixture moments are
computed as in Eqs. 2.33 and 2.34. The collision operator is, therefore,

$\displaystyle\mathcal{Q}^{\text{BGK}}_{p,\,q,\,{\boldsymbol{j}}}[f_{p},f_{q}]=\nu_{p,\,q}(\mathcal{M}_{p,\,q,\,{\boldsymbol{j}}}[f_{p},f_{q}]-f_{p,\,{\boldsymbol{j}}}),$
(3.5)

where
$\mathcal{M}_{p,\,q,\,{\boldsymbol{j}}}[f_{p},f_{q}]=\mathcal{M}_{p,\,q}[f_{p},f_{q}](\mathbf{v}_{\boldsymbol{j}})$,
and $\mathcal{M}_{p,\,q}[f_{p},f_{q}]$ is the Maxwellian (2.30) generated by
the approximate moments.

##### Space

The spatial discretisation of Equation 3.3 will be guided by the treatment of
the transport term alone. Linear transport terms such as
$\mathbf{v}_{{\boldsymbol{j}}}\cdot\nabla_{\mathbf{x}}f_{p,\,{\boldsymbol{j}}}$
can be discretised in a number of ways, including semi-Lagrangian methods [19,
18, 22] and finite-volume methods [43, 57, 26, 44]. Here we shall prescribe a
first-order finite-volume scheme, but this can be replaced with a higher-
resolution method (a WENO scheme [45], for instance) independently of the
choices for the discretisations of velocity and time.

The spatial domain, $\Omega$, is assumed to be rectangular, of the form
$(0,L_{x,1})\times\cdots\times(0,L_{x,{D_{\mathbf{x}}}})$. Each dimension $d$
is discretised into $N_{x,d}$ volumes of size $\Delta x_{d}=L_{x,d}/N_{x,d}$:
every one, a cell $C_{i}=(x_{i-1/2},x_{i+1/2})$ with centre $x_{i}$, where
$x_{i}=i\Delta x$, for $i\in\mathcal{N}_{x}\coloneqq\\{1,\cdots,N_{x}\\}$. The
discretised spatial domain becomes
$\mathcal{X}=\\{\boldsymbol{C}_{\boldsymbol{i}}=C_{i_{1}}\times\cdots\times
C_{i_{D_{\mathbf{x}}}}\\}$, where
${\boldsymbol{i}}=(i_{1},\cdots,i_{D_{\mathbf{x}}})\in\mathcal{N}_{x}^{D_{\mathbf{x}}}$
is a multi-index. Equation 3.3 thus becomes a system of equations:

$\displaystyle\left\\{\begin{array}[]{{>{\displaystyle\vphantom{\frac{1}{1}}}l}}\partial_{t}f_{p,\,{\boldsymbol{i}},\,{\boldsymbol{j}}}+\Phi_{p,\,{\boldsymbol{i}},\,{\boldsymbol{j}}}[f_{p}]=\frac{1}{\varepsilon}\mathcal{Q}_{p,\,{\boldsymbol{i}},\,{\boldsymbol{j}}}[f_{p}],\quad{\boldsymbol{i}}\in\mathcal{N}_{x}^{D_{\mathbf{x}}},\,{\boldsymbol{j}}\in\mathcal{N}_{v}^{D_{\mathbf{v}}},\,t>0,\\\
\mathcal{Q}_{p,\,{\boldsymbol{i}},\,{\boldsymbol{j}}}[f_{p}]=\sum_{q=1}^{P}\mathcal{Q}_{p,\,q,\,{\boldsymbol{i}},\,{\boldsymbol{j}}}[f_{p},f_{q}],\end{array}\right.$
(3.9)

where $f_{p,\,{\boldsymbol{i}},\,{\boldsymbol{j}}}$ is a finite-volume
approximation of $f_{p,\,{\boldsymbol{j}}}$,
$f_{p,\,{\boldsymbol{i}},\,{\boldsymbol{j}}}(t)\simeq\frac{1}{|\boldsymbol{C}_{\boldsymbol{i}}|}\int_{\boldsymbol{C}_{\boldsymbol{i}}}f_{p,\,{\boldsymbol{j}}}(t,\mathbf{x})\mathop{}\\!\mathrm{d}\mathbf{x},$
with equality for the initial datum.

The term $\Phi_{p,\,{\boldsymbol{i}},\,{\boldsymbol{j}}}$, expresses the net
flux across the boundary of $\boldsymbol{C}_{\boldsymbol{i}}$. In the one-
dimensional setting, ${D_{\mathbf{x}}}=1$, it is written in conservation form
as

$\displaystyle\Phi_{p,\,i,\,j}[f_{p}]=\frac{F_{p,\,i+1/2,\,j}-F_{p,\,i-1/2,\,j}}{\Delta
x},$ (3.10)

where

$\displaystyle
F_{p,\,i+1/2,\,j}=f_{p,\,i,\,j}(v_{j})^{+}+f_{p,\,i+1,\,j}(v_{j})^{-}$ (3.11)

is the upwind flux discretisation. Respectively, $(v_{j})^{+}$ and
$(v_{j})^{-}$ are the positive and negative parts of $v_{j}$;
$(s)^{+}=(-s)^{-}=\max\\{s,0\\}$ for any $s\in{\mathbb{R}}$. The
generalisation to higher dimensions is immediate, see any of [43, 57, 26, 44].

The discretisation in space replaces the collision operator (3.5) with

$\displaystyle\mathcal{Q}^{\text{BGK}}_{p,\,q,\,{\boldsymbol{i}},\,{\boldsymbol{j}}}[f_{p},f_{q}]=\nu_{p,\,q}(\mathcal{M}_{p,\,q,\,{\boldsymbol{i}},\,{\boldsymbol{j}}}[f_{p},f_{q}]-f_{p,\,{\boldsymbol{i}},\,{\boldsymbol{j}}}),$
(3.12)

where $\mathcal{M}_{p,\,q,\,{\boldsymbol{i}},\,{\boldsymbol{j}}}[f_{p},f_{q}]$
is simply the Maxwellian computed from the discretisations of $f_{p},f_{q}$ in
velocity and space. The moments given in (3.4) are easily discretised; for
example,

$\displaystyle n_{p}(t,\mathbf{x})=\sum
f_{p,\,{\boldsymbol{j}}}(t,\mathbf{x})\Delta\mathbf{v}\quad\text{becomes}\quad
n_{p,\,{\boldsymbol{i}}}(t)=\sum
f_{p,\,{\boldsymbol{i}},\,{\boldsymbol{j}}}(t)\Delta\mathbf{v}.$ (3.13)

###  3.2 Projective Integration

Following the phase space discretisation, the multispecies BGK equation leads
to the semi-discrete system (3.9). This can be summarised as

$\displaystyle\partial_{t}f_{p,\,{\boldsymbol{i}},\,{\boldsymbol{j}}}=\mathcal{D}^{\varepsilon}[f_{p,\,{\boldsymbol{i}},\,{\boldsymbol{j}}}],$
(3.14)

where $p$, ${\boldsymbol{i}}$, and ${\boldsymbol{j}}$ are, respectively, the
species, spatial, and velocity indices, and where $\mathcal{D}^{\varepsilon}$
is an operator which comprises both discrete transport and collisions:

$\displaystyle\mathcal{D}^{\varepsilon}[f_{p,\,{\boldsymbol{i}},\,{\boldsymbol{j}}}]\coloneqq-\Phi_{p,\,{\boldsymbol{i}},\,{\boldsymbol{j}}}[f_{p}]+\frac{1}{\varepsilon}\mathcal{Q}_{p,\,{\boldsymbol{i}},\,{\boldsymbol{j}}}[f_{p}].$
(3.15)

Projective integration essentially consists of two stages: an inner step,
where an elementary time-stepping method (the inner integrator) is used to
solve (3.14) over a short interval; and a projective step, where the previous
information is used to approximate the time derivative of $f$ and extrapolate
forward (through the outer integrator) over a large time interval.

We define step sizes $\Delta t_{0}$ and $\Delta t_{1}$, respectively for the
inner and outer levels. The solution
$f_{p,\,{\boldsymbol{i}},\,{\boldsymbol{j}}}(t)$ is approximated by
$f_{p,\,{\boldsymbol{i}},\,{\boldsymbol{j}}}^{n,\,k}$ at the time $t=n\Delta
t_{1}+k\Delta t_{0}$. The evolution is computed as follows:

##### Inner Integrator

The inner step is computed through the forward Euler method with step size
$\Delta t_{0}$:

$\displaystyle f^{n,\,k+1}=\operatorname{FE}_{\Delta
t_{0}}[f^{n,\,k}]\coloneqq f^{n,\,k}+\Delta
t_{0}\mathcal{D}^{\varepsilon}[f^{n,\,k}],$ (3.16)

where we have omitted the $p$, ${\boldsymbol{i}}$, and ${\boldsymbol{j}}$
indices for simplicity. The inner step size will be chosen as $\Delta
t_{0}=\mathcal{O}(\varepsilon)$ to ensure the stability of (3.16), just as it
would be in a classical scheme.

While higher-order inner integrators could easily be constructed, it is shown
in [49] that forward Euler is the choice which gives the projective scheme the
best stability properties. If high-order accuracy is desired, this can be
achieved at the outer level.

##### Outer Integrator

The outer step computes the update $f^{n+1,\,0}$ from $f^{n,\,0}$, which can
be simply be understood as the steps $f^{n+1}$ and $f^{n}$ of a time
discretisation $t=n\Delta t_{1}$. The outer step size will be chosen as
$\Delta t_{1}=\mathcal{O}(\Delta x)$.

In order to compute the update, the method first computes a sequence of $K+1$
inner updates from $f^{n,\,0}$, using the inner integrator. The last two
updates, $f^{n,\,K}$ and $f^{n,\,K+1}$, are used to project the solution
forward:

$\displaystyle f^{n+1}=\operatorname{PFE}_{\Delta t_{1}}[f^{n}]\coloneqq
f^{n,\,K+1}+(\Delta t_{1}-(K+1)\Delta
t_{0})\frac{f^{n,\,K+1}-f^{n,\,K}}{\Delta t_{0}};$ (3.17)

the relative size of the extrapolation is $M=\frac{\Delta t_{1}}{\Delta
t_{0}}-(K+1).$

Since the extrapolation is performed with a first-order approximation of the
derivative, this is known as the projective forward Euler method. Higher order
methods can be constructed by using more points to approximate the time
derivative, in the vein of Runge-Kutta methods, as discussed in [30].

A schematic summary of the method can be found in Algorithm 1. The exact
choice of $K$, $\Delta t_{0}$, and $\Delta t_{1}$ will be discussed in Section
4.

Data: $N$, $K$, $\Delta t_{1}$, $\Delta t_{0}$

$f^{0}\leftarrow f(0)$ ;

for  _$0\leq n\leq N-1$_ do

$f^{n,\,0}\leftarrow f^{n}$ ;

for  _$0\leq k\leq K$_ do

$f^{n,\,k+1}\leftarrow\operatorname{FE}_{\Delta t_{0}}[f^{n,\,k}]$ ;

end for

$f^{n+1}\leftarrow\operatorname{PFE}_{\Delta t_{1}}[f^{n}]$ ;

end for

Algorithm 1 Projective Integration

##### Efficiency of the Method

As discussed in [31], the efficiency of projective integration can be measured
by comparing the number of inner steps required to integrate over a time
interval directly to that required by the projective method. Assuming a
homogeneous cost per step and negligible overhead arising from the projection
step (which is reasonable as the projective step does not involve the
computation of moments of $f$), projective integration reduces the overall
computational cost of a simulation by a factor of

$\displaystyle S=\frac{\Delta t_{1}\left(\Delta
t_{0}\right)^{-1}}{K+1}=\frac{M+K+1}{K+1}.$ (3.18)

###  3.3 Telescopic Integration

The stability properties of projective integration are well-understood: the
method is suited for problems with exactly two timescales which are spectrally
separate. There are, however, many multi-scale problems of interest whose
spectra lack the sufficient structure: for instance, the linearised and the
full BGK models, as soon as the collision rate $\nu$ depends on space or on
the distribution $f$ [49, 48, 38]. This is also the case for the multi-species
BGK model (2.19, 2.31), as will be shown in Section 4.

Telescopic integration extends the projective method to include such problems.
The method nests several levels: an innermost level, consisting, as before, of
an elementary integrator; and a number of outer levels, each performing a
projection step based on the iterations of the level below it. We briefly
discuss an approach with two projective levels, which will prove sufficient
for our needs, though the method can easily be extended.

We define step sizes $\Delta t_{0}$, $\Delta t_{1}$, and $\Delta t_{2}$,
respectively for the inner, middle, and outer levels. The solution
$f_{p,\,{\boldsymbol{i}},\,{\boldsymbol{j}}}(t)$ is approximated by
$f_{p,\,{\boldsymbol{i}},\,{\boldsymbol{j}}}^{n,\,k,\,l}$ at the time
$t=n\Delta t_{2}+k\Delta t_{1}+l\Delta t_{0}$. The evolution is computed as
follows:

The step of the outer integrator is given by

$\displaystyle f^{n+1}=\operatorname{PFE}_{\Delta t_{2}}[f^{n}]\coloneqq
f^{n,\,K_{1}+1}+(\Delta t_{2}-(K_{1}+1)\Delta
t_{1})\frac{f^{n,\,K_{1}+1}-f^{n,\,K_{1}}}{\Delta t_{1}};$ (3.19)

the relative size of the extrapolation is $M_{1}=\frac{\Delta t_{2}}{\Delta
t_{1}}-(K_{1}+1).$ In order to evaluate the update, $K_{1}+1$ steps of the
middle integrator have to be performed, to find $f^{n,\,K_{1}}$ and
$f^{n,\,K_{1}+1}$ from $f^{n}$. Each of those steps is given by

$\displaystyle f^{n,\,k+1}=\operatorname{PFE}_{\Delta
t_{1}}[f^{n,\,k}]\coloneqq f^{n,\,k,\,K_{0}+1}+(\Delta t_{1}-(K_{0}+1)\Delta
t_{0})\frac{f^{n,\,k,\,K_{0}+1}-f^{n,\,k,\,K_{0}}}{\Delta t_{0}},$ (3.20)

extrapolations of size $M_{0}=\frac{\Delta t_{1}}{\Delta t_{0}}-(K_{0}+1).$
For each step to be performed, $K_{0}+1$ steps of the inner integrator have to
be computed, to find $f^{n,\,k,\,K_{0}}$ and $f^{n,\,k,\,K_{0}+1}$ from
$f^{n,\,k}$. These inner steps are, once again, given by the forward Euler
method

$\displaystyle f^{n,\,k,\,l+1}=\operatorname{FE}_{\Delta
t_{0}}[f^{n,\,k,\,l}]\coloneqq f^{n,\,k,\,l}+\Delta
t_{0}\mathcal{D}^{\varepsilon}[f^{n,\,k,\,l}],$ (3.21)

and can be computed directly. Since the outermost extrapolation is once again
performed with a first-order approximation of the derivative, this is known as
the (two-level) telescopic projective forward Euler method.

A schematic summary of the method is presented in Algorithm 1. Once again, the
choice of parameters $K_{0}$, $K_{1}$, $\Delta t_{0}$, $\Delta t_{1}$, and
$\Delta t_{2}$ will be discussed in Section 4.

Data: $N$, $K_{0}$, $K_{1}$, $\Delta t_{2}$, $\Delta t_{1}$, $\Delta t_{0}$

$f^{0}\leftarrow f(0)$ ;

for  _$0\leq n\leq N-1$_ do

$f^{n,\,0}\leftarrow f^{n}$ ;

for  _$0\leq k\leq K_{0}$_ do

$f^{n,\,k,\,0}\leftarrow f^{n,\,k}$ ;

for  _$0\leq l\leq K_{1}$_ do

$f^{n,\,k,\,l+1}\leftarrow\operatorname{FE}_{\Delta t_{0}}[f^{n,\,k,\,l}]$ ;

end for

$f^{n,\,k+1}\leftarrow\operatorname{PFE}_{\Delta t_{1}}[f^{n,\,k}]$ ;

end for

$f^{n+1}\leftarrow\operatorname{PFE}_{\Delta t_{2}}[f^{n}]$ ;

end for

Algorithm 2 Telescopic Integration (two levels)

##### Efficiency of the Method

The efficiency of telescopic projective integration can be measured as in the
previous section. The $n\textsuperscript{th}$ level reduces the computational
cost of a simulation by a factor of $S_{n}$, and the overall efficiency factor
is $S$:

$\displaystyle S_{n}=\frac{\Delta t_{n+1}\left(\Delta
t_{n}\right)^{-1}}{K_{n}+1}=\frac{M_{n}+K_{n}+1}{K_{n}+1},\quad
S=\prod_{n=1}^{N}S_{n}.$ (3.22)

##  4 Spectra & Stability

Figure 1:  Trajectories in solution space for systems with Knudsen number in
the hydrodynamic regime.  Figure 2:  Sketch of projective integration: the
method performs steps of size $\mathcal{O}(1)$, but introduces intermediate
steps of size $\mathcal{O}(\varepsilon)$ to ensure the overall stability of
the numerical solution.  Figure 3:  Asymptotic stability region for the
projective forward Euler integrator.  Figure 4:  Asymptotic stability region
for the telescopic projective forward Euler integrator.  Figure 5:
Approximate spectra for mixture Sod Tube problem of Section 5.1: various mass
ratios. Top: fast modes, $\mathcal{O}(\varepsilon)$ (split across two
regions). Bottom left: middle modes. Bottom right: slow modes,
$\mathcal{O}(1)$.

The purpose of this section is to discuss the choice of integration
parameters, $\Delta t_{i}$ and $K_{i}$. We shall recall the stability analysis
of [30, 31], and illustrate it in the context of kinetic equations, and the
multispecies BGK model in particular.

The suitability of projective integration as a numerical method for kinetic
equations with very small Knudsen numbers can be understood intuitively by
considering their hydrodynamic asymptotics. As referenced in Section 2.4,
Equations 2.4 and 2.19 can be studied in the limit $\varepsilon\rightarrow 0$,
and their limiting behaviour can be characterised by systems of equations for
the evolution of their moments. In the space of solutions, as the collision
operator drives the state towards local equilibrium, trajectories rapidly
arrive at a hypersurface comprised only of Maxwellian distributions, see Fig.
1 for a diagram. Once on this surface, the only possible evolution is within,
which corresponds to the evolution of the distribution’s moments as governed
by the corresponding macroscopic system. For small values of $\varepsilon$,
the arrival to and the evolution within the Maxwellian surface occur in very
different time scales.

As discussed in Section 3, resolving both of these time scales with a typical
scheme proves very costly, since a typical integrator will require steps of
size $\Delta t=\mathcal{O}(\varepsilon)$ for stability. However, the two-scale
nature of the problem can be exploited to escape this restriction: the
integrator is permitted large steps (whose size are limited only by a
hyperbolic CFL-like condition), as long as they alternate with smaller steps
which ensure the return to the equilibrium hypersurface, as sketched in Fig.
2.

This intuitive portrayal can be formalised by performing error stability
analysis on the projective integration and telescopic integration methods. An
asymptotic stability result is presented in [30]: for large $M$, the stability
region for the projective Euler method approaches two disjoint disks on the
$\Delta t_{0}\lambda$ complex plane; $\mathcal{D}\left(-1,M^{-1/K}\right)$ and
$\mathcal{D}\left(-M^{-1},M^{-1}\right)$, where by
$\mathcal{D}\left(c,r\right)$ we mean the disk centred at the complex number
$c$ with radius $r$ (see Fig. 3). The two-disk stability region is perfectly
suited for a two-cluster spectrum, which often lies behind the fast-slow
dynamics described above. Once the spectrum is known, the parameters $\Delta
t_{0}$, $K$, and $M$ can be chosen so that the scaled eigenmodes, $\Delta
t_{0}\lambda$, lie within the two stable disks: the fast modes will be
contained in the left disk, and the slow modes, in the right disk. As a
curiosity, note that the right disk corresponds to the stability region of a
standard forward Euler method with step $\Delta t_{0}M$; this is indeed the
step size of the projection (3.17).

In more complex settings, a two-disk stability region might not suffice. For
instance, the BGK model with a density-dependent collision frequency, the full
Boltzmann equation, and the multispecies BGK model which we consider in this
work all lead to more complicated spectra. Fortunately, telescopic projective
methods also lead to richer stability regions, and can be used for our
purposes. In the same vein as [30], we find that the two-level telescopic
projective Euler method comprises three disks when $M_{0}$ and $M_{1}$ are
large: $\mathcal{D}\left(-1,M_{0}^{-1/K_{0}}M_{1}^{-1/K_{1}}\right)$,
$\mathcal{D}\left(-M_{0}^{-1},M_{0}^{-1}M_{1}^{-1/K_{1}}\right)$, and
$\mathcal{D}\left(-M_{0}^{-1}M_{1}^{-1},M_{0}^{-1}M_{1}^{-1}\right)$ (see Fig.
4). As before, the parameters $\Delta t_{0}$, $K_{0}$, $K_{1}$, $M_{0}$, and
$M_{1}$ have to be chosen so that the scaled eigenmodes lie within the
stability region. This three-disk configuration is featured in the spectra of
the numerical experiments discussed in Sections 5.2, 5.3 and 5.4.

###  4.1 Estimating Spectra

The parameter choice strategy presented above relies entirely on accurate
spectral information for the problems in question. However, finding the
spectrum of kinetic equations, or even establishing partial information (such
as bounds, number of clusters, or spectral gaps) is a difficult problem. For
the single species Boltzmann equation, we refer to the classic paper of [24],
its ${L^{1}}$ extension [51], or its ${L^{2}}$ extension [32]. The
multispecies theory is sparse, though we highlight the recent works [10, 20].

Our approach in this work will be to estimate the spectra of the problems
numerically. We will consider the semi-discrete problem (3.14), evaluate the
Jacobian of the operator $\mathcal{D}^{\varepsilon}$ through finite
differencing, and compute the eigenvalues of the resulting matrix as an
approximate spectrum for the problem. This approach proves extremely
successful, as it permits finding optimal integration parameters for all of
the numerical experiments of Section 5. Furthermore, exploiting the mild
dependence of the spectra in the dimension of the problem [24], we are able to
use the spectra of one-dimensional problems in order to find suitable
parameters for problems in higher dimensions.

We present in Fig. 5 the spectra used in the mixture Sod Tube problem of
Section 5.1 for mass ratios ranging from $m_{1}/m_{2}=25$ to
$m_{1}/m_{2}=100$. The slow modes and middle modes, shown in the bottom right
and bottom left plots respectively, are seen not to vary significantly across
the different ratios, and therefore have similar stability properties.

However, this is not the case for the fastest modes. The fast cluster is shown
across two plot panels in Fig. 5: the top left panel comprises the majority of
fast modes, which vary slightly; the top right panel shows the rightmost fast
modes, which vary greatly as the mass ratio increases. Across the different
ratios, the width of the fast cluster grows by an order of magnitude. This
spreading of the fast cluster, and the associated reduction in the gap, is
consistent with the results in [10].

In practice, the widening of the cluster leads to different stability
requirements: in the corresponding numerical test, the parameters used in the
$m_{1}/m_{2}=1$ case also work for the $m_{1}/m_{2}=100$ case, with the
exception of $K_{1}$; the step number of the second level is increased from
$6$ to $14$. This logarithmic scaling of the computational cost with the
stiffness of the problem has been reported in [47] in the single-species
setting.

##  5 Numerical Experiments

In this section, we apply our schemes to various physical scenarios in order
to showcase the robustness and versatility of the methods. We will consider
one example with ${D_{\mathbf{x}}}={D_{\mathbf{v}}}=1$, the mixture Sod Tube,
and three examples with ${D_{\mathbf{x}}}={D_{\mathbf{v}}}=2$: a shock-bubble
interaction, a Kelvin-Helmholtz instability, and a Richtmyer-Meshkov
instability. The examples demonstrate the ability of the schemes to deal with
complex scenarios, as well as their superior efficiency.

###  5.1 Mixture Sod Tube & Extreme Mass Ratios (1D / 1D)

Figure 6:  Mixture Sod Tube problem: numerical solutions with decreasing
Knudsen number and unit mass ratio. The solutions approach the correct
$\varepsilon\rightarrow 0$ limit.  Figure 7:  Mixture Sod Tube problem:
extreme mass ratios. Left: higher mass ratios yield slower hydrodynamic
convergence. Right: high mass ratios require larger velocity domains
$(-L_{v},L_{v})$ in order to capture the correct behaviour.

Our first experiment is a generalisation of the Sod Tube problem [55]. The Sod
Tube is a generic Riemann problem for the Euler system (2.36) (with a single
species). A discontinuous initial datum is prescribed, which immediately
develops a range of hyperbolic phenomena: a rarefaction wave, a shock wave,
and a contact discontinuity. Nevertheless, the solution is analytically
tractable, yielding an ideal test case for numerical schemes.

We shall employ here a setting twice removed from the classical Sod Tube.
Firstly, we entertain an analogue of the Riemann problem which involves two
initially separate gases, whose solution can be constructed in the same
fashion as that of the classical problem; see the appendix for the derivation.
Secondly, rather than solving a multifluid Euler system, we solve the
corresponding problem for the multispecies BGK model, and compare its moments
to the aforementioned analytic solution, exploiting the limiting behaviour
discussed in Section 2.4.

####  5.1.1 Mixture Sod Tube

We will consider two gases in a one-dimensional domain $\Omega=(0,1)$. The
mass ratio $m_{2}/m_{1}$ and Knudsen number will vary through the examples.
The initial configuration will emulate a Riemann problem with left state
$(\rho_{L},v_{L},\mathcal{P}_{L})=(1,0,1)$ and right state
$(\rho_{R},v_{R},\mathcal{P}_{R})=(2^{-3},0,2^{-5})$, where the left state is
entirely comprised of the first gas, and the right state of the second. To
that end, we prescribe their initial distributions as Maxwellians with moments

$\left\\{\begin{array}[]{lllll}\displaystyle\rho_{1}=(1-\delta)\rho_{L},&\rho_{2}=\delta\rho_{L},&\bar{v}_{1}=\bar{v}_{2}=v_{L},&\mathcal{P}_{1}=\mathcal{P}_{2}=\mathcal{P}_{L},&\text{if
}x\leq 0.5;\\\
\displaystyle\rho_{1}=\delta\rho_{R},&\rho_{2}=(1-\delta)\rho_{R},&\bar{v}_{1}=\bar{v}_{2}=v_{R},&\mathcal{P}_{1}=\mathcal{P}_{2}=\mathcal{P}_{R},&\text{if
}x>0.5.\end{array}\right.$ (5.1)

Ideally, we would set $\delta=0$, but this leads to an ill-defined
temperature, so we let $\delta=10^{-5}$. This pattern will also be used in
later experiments whenever a density is zero. We henceforth refer to this
problem as the mixture Sod Tube problem.

To begin, we verify the behaviour of this problem in the hydrodynamic limit.
We will compute numerically the solution with datum Eq. 5.1 with mass ratio
$m_{2}/m_{1}=1$, for three different values of the Knudsen number:
$\varepsilon=10^{-1}$, $10^{-2}$, and $10^{-6}$. The larger values of epsilon
do not pose significant stiffness, so a direct integration method can be used;
for $\varepsilon=10^{-6}$, we resort to a telescopic two-level method.

The solution is computed over the time interval $t\in(0,0.15)$. The domain
$\Omega$ is discretised with $\Delta x=2^{-10}$, and the velocity space is set
as $(-20,20)$, with $\Delta v=2^{-4}$. For the direct method employed on the
larger values of $\varepsilon$, we let $\Delta
t=$1.53\text{\times}{10}^{-05}$$. In the telescopic method, we choose $\Delta
t_{0}=$5\text{\times}{10}^{-07}$$, $\Delta t_{1}=$2\text{\times}{10}^{-06}$$,
and $\Delta t_{2}=$6.1\text{\times}{10}^{-05}$$, and step numbers $K_{0}=1$
and $K_{1}=6$. We impose no-flux boundary conditions.

Figure 6 shows the numerical solutions superimposed on the analytical limiting
solution (computed through the procedure detailed in the appendix). As the
Knudsen number decreases, the moments of the solution approach the correct
hydrodynamic limit.

####  5.1.2 A CPU Benchmark

$\varepsilon$ | CPU Time ($\mathrm{s}$) | Improvement factor
---|---|---
|  Direct |  Telescopic |  Real |  Theoretical
$1\text{\times}{10}^{-5}$ |  $54$ |  $47$ |  $1.148\,936\,17$ |  $3.487\,723\,21$
$1\text{\times}{10}^{-6}$ |  $508$ |  $47$ |  $10.808\,510\,638$ |  $34.877\,232\,142\,857\,146$
$1\text{\times}{10}^{-7}$ |  $5660$ |  $48$ |  $117.916\,666\,667$ |  $348.772\,321\,428\,571\,44$
Table 1:  CPU Time (rounded to nearest second) and improvement factor for the
benchmark of Section 5.1.2

We now revisit the problem (5.1) to conduct a CPU benchmark. We will solve the
problem for $\varepsilon=10^{-5}$, $\varepsilon=10^{-6}$, and
$\varepsilon=10^{-7}$ with the telescopic method as well as the direct method,
and compare the

The domain $\Omega$ is discretised with $\Delta x=2^{-8}$, and the velocity
space is set as $(-20,20)$, with $\Delta v=2^{-4}$. For the direct method, we
let $\Delta t=\varepsilon/2$. In the telescopic method, we choose $\Delta
t_{0}=\varepsilon/2$, $\Delta t_{1}=2\varepsilon$, and $\Delta t_{2}=\Delta
x/16$, and step numbers $K_{0}=1$ and $K_{1}=6$. We impose no-flux boundary
conditions.

Table 1 shows the CPU time of each simulation, the real improvement factor
(the ratio of CPU times of the direct method to the telescopic method), and
the theoretical improvement factor as defined in Eq. 3.22. The theoretical
factor is overly optimistic, as it assumes that the projection steps have
negligible computational cost; nevertheless, the telescopic method is faster
for $\varepsilon=10^{-5}$ and vastly superior for smaller values of the
Knudsen number.

####  5.1.3 Extreme Mass Ratios

It is of great interest to attempt the numerical solution of problems with
extremely large mass ratios. Indeed, a simple mixture of Ar and He exhibits a
mass ratio of $10$, which can be a problem for some numerical methods, as
noted in [61]. Their method can deal with mass ratios up to $35$; however,
ratios twice as large can easily be found in other scenarios, such as mixtures
of H2 and Xe.

We will demonstrate the behaviour of our scheme in the context of a
hydrodynamic limit under an extreme mass ratio. Conveniently, the asymptotic
behaviour of the mixture Sod Tube problem remains unchanged if we assume that
the gases have different molecular masses; therefore, it remains a suitable
validation case.

We solve the problem for mass ratios $m_{2}/m_{1}=5$, $30$, and $100$, in the
hydrodynamic regime and compare the solutions. We let $\varepsilon=10^{-6}$.
The spatial discretisation is done as before. The velocity space is set as
$(-60,60)$, with $\Delta v=2^{-4}$; such large domain is unnecessary for the
smaller mass ratios, but will be required for $m_{2}/m_{1}=100$. We again
choose $\Delta t_{0}=$5\text{\times}{10}^{-07}$$, $\Delta
t_{1}=$2\text{\times}{10}^{-06}$$, and $\Delta
t_{2}=$6.1\text{\times}{10}^{-05}$$, though this time we set $K_{0}=1$ and
$K_{1}=14$; again, the large step number is only required for the larger mass
ratios. In order to justify our velocity discretisation, we will also solve
the $m_{2}/m_{1}=100$ case with varying velocity spaces $(-L_{v},L_{v})$, for
$L_{v}=40$, $60$, and $80$, keeping the rest of the parameters fixed.

In all cases we deal with the large velocity supports directly. The recent
work [11] has used a rescaling velocity approach reminiscent of [27] to
overcome the same issue, but their strategy remains limited to mass ratios up
to $20$.

Figure 7 shows the numerical solutions, again superimposed on the limiting
solution. The effects of the extreme mass ratios can be seen at point of
contact discontinuity (the boundary between the two gases), magnified in the
figure. The left column shows the solutions with various mass ratios; the
overall hydrodynamic limit is captured well. However, the interfacial effects
are more pronounced as the mass ratio increases, and will require a smaller
Knudsen number before they become imperceptible. The right column shows the
effect of the choice of velocity domain in the solution: $L_{v}=40$ leads to
widespread error, whereas $L_{v}=60$ recovers the correct behaviour, and is in
fact indistinguishable from $L_{v}=80$.

###  5.2 Shock-Bubble Interaction (2D / 2D)

Physical | Phase space discretisation | Time discretisation
---|---|---
$m_{1}$ |  $1$ |  $\Delta x$ |  $0.025$ |  |  $K_{0}$ |  $1$ |
$m_{2}$ |  $5$ |  $\Delta y$ |  $0.025$ |  ($=\Delta x$) |  $K_{1}$ |  $6$ |
$\varepsilon$ |  $1\text{\times}{10}^{-5}$ |  $\Delta v$ |  $0.25$ |  |  $\Delta t_{0}$ |  $0.5\text{\times}{10}^{-5}$ |  ($=\varepsilon/2$)
|  |  |  |  |  $\Delta t_{1}$ |  $2\text{\times}{10}^{-5}$ |  ($=2\varepsilon$)
|  |  |  |  |  $\Delta t_{2}$ |  $0.001\,25$ |  ($=\Delta x/20$)
Table 2:  Parameters for the shock-bubble interaction test of Section 5.2, see
Figs. 8, 9 and 10 Figure 8:  Shock-bubble interaction: evolution of the first
(lighter) species. Left: datum. Middle: shock interacts with bubble. Right:
pressure waves arise as a result of the shock-bubble interaction. Simulation
parameters can be found in Table 2.  Figure 9:  Shock-bubble interaction:
evolution of the densities. Left: shock interacts with bubble. Right: bubble
is dragged and deformed by the stream. Simulation parameters can be found in
Table 2.  Figure 10:  Shock-bubble interaction: detail of the bubble (second
species). Left: the bubble has been compressed by the shock. Right: the bubble
is displaced by the stream and deforms into a “C shape”. Simulation parameters
can be found in Table 2.

We investigate the interaction between a travelling shock and a smooth
stationary bubble. This is a multispecies adaptation of a one-species test
case proposed in [58]; the original test has been used to validate numerical
schemes for the BGK model [13], including projective integration [48].

We will consider a mixture of two gases with mass ratio $m_{2}/m_{1}=5$ in a
rectangular domain $\Omega=(-1.5,3)\times(-1.5,1.5)$. The Knudsen number is
chosen in the hydrodynamic regime, $\varepsilon=$1\text{\times}{10}^{-5}$$.
The initial configuration of the first (lighter) gas is a normal shock wave of
Mach number $\mbox{\text{Ma}}=2$ propagating in the positive $x$-direction.
Its initial distribution, $f_{1}(0,\mathbf{x},\mathbf{v})$, is chosen as the
Maxwellian corresponding to the following Riemann datum:

$\left\\{\begin{array}[]{lllll}\displaystyle\rho_{1}=2,&\bar{v}_{1,x}=1.414,&\bar{v}_{1,y}=0,&T_{1}=2.5,&\text{if
}x\leq-1;\\\
\displaystyle\rho_{1}=1,&\bar{v}_{1,x}=0,&\bar{v}_{1,y}=0,&T_{1}=1,&\text{if
}x>-1.\end{array}\right.$ (5.2)

The second (heavier) gas is initially stationary: a Maxwellian with density

$\displaystyle\rho_{2}(0,\mathbf{x})=\exp\left\\{-16\left|\mathbf{x}-\mathbf{x}_{0}\right|^{2}\right\\}.$
(5.3)

The temperatures of both gases around the bubble are chosen equal to each
other and in such a way that there is unit total pressure:

$\displaystyle
T_{1}(0,\mathbf{x})=T_{2}(0,\mathbf{x})=\left(n_{1}(0,\mathbf{x})+n_{2}(0,\mathbf{x})\right)^{-1},$
(5.4)

which is consistent with the right state in Eq. 5.2.

Numerically, the domain $\Omega$ is discretised with $\Delta x=\Delta
y=$0.025$$ ($180\times 120$ cells). The velocity space is set as
$(-12,12)^{2}$, with $\Delta v=$0.25$$ ($96^{2}$ cells). In total, the phase
space is discretised with $199,065,600$ cells. We impose outflow boundary
conditions.

The solution corresponding to Eqs. 5.2, 5.3 and 5.4 is computed over the time
$t\in(0,1.5)$. We employ the telescopic two-level method, with step sizes
$\Delta t_{0}=$0.5\text{\times}{10}^{-5}$$, $\Delta
t_{1}=$2\text{\times}{10}^{-5}$$, and $\Delta t_{2}=$0.001\,25$$, and step
numbers $K_{0}=1$ and $K_{1}=6$. The inner steps follow the pattern $4\Delta
t_{0}=2\varepsilon=\Delta t_{1}$. The outermost step is restricted by the
stability of the transport scheme, rather than the projective integration. The
use of telescopic projective integration decreases the computational cost by a
factor of $S\simeq 18$, as defined in Eq. 3.22. This factor, considering the
fine and high-dimensional mesh employed here, is extremely beneficial,
reducing the simulation time from a matter of days to a matter of hours!

Figure 8 shows the evolution of the first (lighter) gas at several times. The
shock travels to meet the bubble; when they meet, the central ($y\simeq 0$)
portion of the shock is slowed down. As the shock traverses the bubble, two
pressure waves are formed: one travelling upstream, appears as a reflection
from the initial shock-bubble interaction; the second, travelling downstream,
arises as the shock surrounds the obstacle and both “arms” meet behind the
bubble.

Figure 9 compares the densities of either gas, and Fig. 10 presents a detailed
view of the bubble. The bubble is compressed by the shock when it first
arrives. As the shock passes, the bubble is displaced downstream and deforms
into a “C shape”.

A summary of the simulation parameters can be found in Table 2.

###  5.3 Kelvin-Helmholtz Instability (2D / 2D)

Physical | Phase space discretisation | Time discretisation
---|---|---
$m_{1}$ |  $1$ |  $\Delta x$ |  $0.007\,812\,5$ |  ($=2^{-7}$) |  $K_{0}$ |  $2$ |
$m_{2}$ |  $5$ |  $\Delta y$ |  $0.007\,812\,5$ |  ($=\Delta x$) |  $K_{1}$ |  $4$ |
$\varepsilon$ |  $1\text{\times}{10}^{-5}$ |  $\Delta v$ |  $0.5$ |  |  $\Delta t_{0}$ |  $0.5\text{\times}{10}^{-5}$ |  ($=\varepsilon/2$)
|  |  |  |  |  $\Delta t_{1}$ |  $2\text{\times}{10}^{-5}$ |  ($=2\varepsilon$)
|  |  |  |  |  $\Delta t_{2}$ |  $0.000\,390\,625$ |  ($=\Delta x/20$)
Table 3:  Parameters for the Kelvin-Helmholtz instability test of Section 5.3.
Figure 11:  Kelvin-Helmholtz instability: evolution in time. Left: undulations
appear on the gas interface. Middle: clear vortices have developed. Right:
vortices begin to smear and merge. Simulation parameters can be found in Table
3.  Figure 12:  Kelvin-Helmholtz instability: modulus of the total velocity
$\bar{\mathbf{v}}$ and stream lines at peak vorticity. Left: whole domain.
Right: detail of a vortex. Simulation parameters can be found in Table 3.

We turn our attention to the Kelvin-Helmholtz instability. This is a well-
known phenomenon where vortices arise at the interface of two fluids with
different density moving at different speeds. This test was also used to
validate projective integration for the single-species BGK model in [48], with
an initial configuration drawn from [46]. Here we shall present a two-species
analogue.

We consider two gases with mass ratio $m_{2}/m_{1}=5$ in a square domain
$\Omega=(-0.5,0.5)^{2}$, periodic along the $x$-direction. The Knudsen number
is chosen in the hydrodynamic regime, $\varepsilon=$1\text{\times}{10}^{-5}$$;
higher values do not yield sufficiently defined vortices. The gases are
initially separate: the first (lighter) gas occupies the $y\geq 0$ region,
with density $\rho_{1}=1$ and horizontal velocity $\bar{v}_{1,x}=0.5$; the
second (heavier) gas occupies the $y<0$ region, with density $\rho_{2}=2$ and
horizontal velocity $\bar{v}_{2,x}=-0.5$. Their temperatures are chosen to
ensure there is equal unit pressure across the boundary: $T_{1}=n_{1}^{-1}$,
$T_{2}=n_{2}^{-1}$. Both gases are given a small vertical velocity to induce
vorticity: $\bar{v}_{1,y}=\bar{v}_{2,y}=$1\text{\times}{10}^{-2}$\sin(4\pi
x).$ The initial distribution of the gases are the Maxwellians corresponding
to these moments.

Numerically, the domain $\Omega$ is discretised with $\Delta x=\Delta
y=$0.007\,812\,5$$ ($128^{2}$ cells). The velocity space is set as
$(-8,8)^{2}$, with $\Delta v=$0.5$$ ($32^{2}$ cells). In total, the phase
space is discretised with $16,777,216$ cells. We impose outflow boundary
conditions along the non-periodic boundaries.

The solution is computed over the time $t\in(0,3)$. We employ the telescopic
two-level method, with step sizes $\Delta t_{0}=$0.5\text{\times}{10}^{-5}$$,
$\Delta t_{1}=$2\text{\times}{10}^{-5}$$, and $\Delta
t_{2}=$0.000\,390\,625$$, and inner steps $K_{0}=2$ and $K_{1}=4$. Again the
inner steps follow the pattern $4\Delta t_{0}=2\varepsilon=\Delta t_{1}$, and
the outermost step is restricted by the hyperbolic CFL condition. The use of
telescopic projective integration decreases the computational cost by a factor
of $S\simeq 5$, as defined in Eq. 3.22.

Figure 11 shows the evolution of the densities of each gas as well as the
total pressure. The small initial vertical perturbation results in the
undulation of the interface. The relative difference in horizontal velocity
rapidly causes the formation of clearly defined vortices. As the shear stress,
the vortices are smeared horizontally.

Figure 12 studies the vortex structure at the time corresponding to the middle
column of Fig. 11, where vorticity is at a peak. The figure displays the
modulus of the total velocity $\bar{\mathbf{v}}$; the streamlines are
superimposed.

A summary of the simulation parameters can be found in Table 3.

###  5.4 Richtmyer-Meshkov Instability (2D / 2D)

Physical | Phase space discretisation | Time discretisation
---|---|---
$m_{1}$ |  $1$ |  $\Delta x$ |  $0.000\,25$ |  |  $K_{0}$ |  $1$ |
$m_{2}$ |  $5$ |  $\Delta y$ |  $0.000\,25$ |  ($=\Delta x$) |  $K_{1}$ |  $6$ |
$\varepsilon$ |  $1\text{\times}{10}^{-6}$ |  $\Delta v$ |  $0.25$ |  |  $\Delta t_{0}$ |  $0.5\text{\times}{10}^{-6}$ |  ($=\varepsilon/2$)
|  |  |  |  |  $\Delta t_{1}$ |  $2\text{\times}{10}^{-6}$ |  ($=2\varepsilon$)
|  |  |  |  |  $\Delta t_{2}$ |  $6.25\text{\times}{10}^{-5}$ |  ($=\Delta x/4$)
Table 4:  Parameters for the Richtmyer-Meshkov instability test of Section
5.4.  Figure 13:  Richtmyer-Meshkov instability: density of the second species
and modulus of the total momentum. Left: after the initial shock interaction,
vortices arise. Right: the interface is deformed into a mushroom-like shape.
Simulation parameters can be found in Table 4.  Figure 14:  Richtmyer-Meshkov
instability: modulus of the total velocity $\bar{\mathbf{v}}$ and stream lines
at the final time. Left: whole domain. Right: detail of the lower vortex.
Simulation parameters can be found in Table 4.  Figure 15:  Richtmyer-Meshkov
instability: detail of the gas interface (density of the second species).
Left: vortices appear. Right: mushroom-like interface is formed. Simulation
parameters can be found in Table 4.

To conclude this section, we study the Richtmyer-Meshkov instability. This is
a phenomenon which takes place when the perturbed interface between a thin and
a dense gas is momentarily accelerated by a passing shock. The misalignment of
the pressure gradient (due to the shock) and the density gradient induces
vorticity, which leads to the formation of a mushroom-like interface666This
instability is not to be confused with the related Rayleigh-Taylor
instability, which arises on a perturbed interface when the dense gas is
accelerated continuously by the thin gas, rather than by a shock. More details
as well as experimental images can be found in [60].. A more detailed physical
background, as well as experimental images, can be found in [2]. This
instability has been used to validate a lattice-Boltzmann method in [54].

We shall consider two gases with mass ratio $m_{2}/m_{1}=5$ in a rectangular
domain $\Omega=(-0.5,0.5)\times(0,0.5)$. The Knudsen number is here chosen as
$\varepsilon=$1\text{\times}{10}^{-6}$$ in order to correctly resolve the
instability. The gases are initially separate: the first (lighter) gas
occupies the $x\leq b(y)$ region, with a perturbed boundary given by
$x=-$1\text{\times}{10}^{-2}$\sin(20\pi y),$ whereas the second (heavier) gas
lies in the rest of the domain. The initial configuration of the first gas is
a normal shock wave of Mach number $\mbox{\text{Ma}}=1.21$ (inspired by [54])
propagating in the positive $x$-direction. Its initial distribution is chosen
according to the following Riemann datum:

$\left\\{\begin{array}[]{lllll}\displaystyle\rho_{1}=1.268,&\bar{v}_{1,x}=0.256,&\bar{v}_{1,y}=0,&p_{1}=0.809,&\text{if
}x\leq s_{0};\\\
\displaystyle\rho_{1}=1,&\bar{v}_{1,x}=0,&\bar{v}_{1,y}=0,&p_{1}=0.5,&\text{if
}x>s_{0};\end{array}\right.$ (5.5)

$s_{0}$, the initial position of the shock, is a negative constant. The second
gas is initially stationary, with density $\rho_{2}=5$, and pressure
$p_{2}=0.5$ to ensure the boundary is initially not forced. The initial
distribution of the gases are the Maxwellians corresponding to these moments.
For convenience, in order to avoid the interface leaving the domain, all the
horizontal velocities are decreased by $0.07$.

Numerically, the domain $\Omega$ is discretised with $\Delta x=\Delta
y=$0.000\,25$$ ($400\times 200$ cells). The velocity space is set as
$(-4,4)^{2}$, with $\Delta v=$0.25$$ ($32^{2}$ cells). In total, the phase
space is discretised with $81,920,000$ cells. We impose outflow boundary
conditions.

The solution is computed over the time $t\in(-0.02,1.0)$. We choose the
negative initial time and let $s_{0}=$0.0242$$ so that the shock crosses the
$x=0$ line exactly at $t=0$. We employ the telescopic two-level method, with
step sizes $\Delta t_{0}=$0.5\text{\times}{10}^{-6}$$, $\Delta
t_{1}=$2\text{\times}{10}^{-6}$$, and $\Delta
t_{2}=$6.25\text{\times}{10}^{-5}$$, and inner steps $K_{0}=1$ and $K_{1}=6$.
Once more, the inner steps follow the pattern $4\Delta
t_{0}=2\varepsilon=\Delta t_{1}$ and the outermost step is restricted by the
hyperbolic CFL condition. The use of telescopic projective integration
decreases the computational cost by a factor of $S\simeq 9$, as defined in Eq.
3.22.

Figure 13 shows the density of the heavier gas and the modulus of the total
momentum of the system at two different times. The gas interface is seen
deforming as it develops a mushroom-like shape. Vortices are visible on either
side of the perturbation, and much of the system’s momentum is found at the
interface.

Figure 14 studies the vortices at the final simulation time. The figure shows
the modulus of the total velocity with streamlines, detailing the high speed
along the lateral gas interface, and offering a magnified view of the lower
vortex.

Figure 15 offers a detailed view of the gas interface. The density of the
heavier gas has been recoloured to omit the extrema, highlighting only the
density gradient. The vortex formation and the deformation of the interface
are particularly clear here.

A summary of the simulation parameters can be found in Table 4.

##### Acknowledgments

TR was partially funded by Labex CEMPI (ANR-11-LABX-0007-01) and ANR Project
MoHyCon (ANR-17-CE40-0027-01). RB was partially funded by Labex CEMPI
(ANR-11-LABX-0007-01) and by the Advanced Grant Nonlocal-CPD (Nonlocal PDEs
for Complex Particle Dynamics: Phase Transitions, Patterns and
Synchronization) of the European Research Council Executive Agency (ERC) under
the European Union’s Horizon 2020 research and innovation programme (grant
agreement No. 883363).

## Appendix: Solving the Mixture Sod Tube Problem

Figure 16:  Structure of the solution to the Sod Tube problem

We recall, for the sake of a complete exposition, the solution to the
classical Sod Tube problem and the mixture Sod Tube problem, as can be found
in the literature of hyperbolic problems [57]. In the usual setting, the
initial Riemann datum is given by the left ($\rho_{L}$, $v_{L}$,
$\mathcal{P}_{L}$) and right ($\rho_{R}$, $v_{R}$, $\mathcal{P}_{R}$) states:

$\left\\{\begin{array}[]{llll}\displaystyle\rho=\rho_{L},&v=v_{L},&\mathcal{P}=\mathcal{P}_{L},&\text{if
}x\leq 0;\\\
\displaystyle\rho=\rho_{R},&v=v_{R},&\mathcal{P}=\mathcal{P}_{R},&\text{if
}x>0.\end{array}\right.$ (A.1)

As the datum evolves, typically we observe two non-linear waves enveloping a
new central region; see Fig. 16 for a diagram. This region comprises two
states, separated by a contact discontinuity. Both have equal speed and
pressure, $v_{C}$, $\mathcal{P}_{C}$, but different densities, $\rho_{C,L}$,
and $\rho_{C,R}$. We will assume here
$\mathcal{P}_{L}\geq\mathcal{P}_{C}\geq\mathcal{P}_{R}$, meaning the leftmost
wave is a rarefaction, and the rightmost wave is a shock.

The solution strategy involves first finding values $v_{C}$, $\mathcal{P}_{C}$
which are consistent with the rarefaction/shock wave structure on either side.
Given a guess for $\mathcal{P}_{C}$, the value of $v_{C}$ can be computed from
the left state using the Riemann invariants (quantities which are constant
across the rarefaction). The precise relation is

$\displaystyle
v_{\text{rare}}(\mathcal{P}_{C})=v_{L}+2c_{L}\left[1-\left(\frac{\mathcal{P}_{C}}{\mathcal{P}_{L}}\right)^{\frac{\gamma-1}{2\gamma}}\right]\left(\gamma-1\right)^{-1},$
(A.2)

where $c_{L}$ is the speed of sound on the left state, computed via
$c=\sqrt{\gamma\mathcal{P}/\rho}$, and where $\gamma$ is the adiabatic
exponent; for an ideal monatomic gas, it is given by $\gamma=1+2/D$, where $D$
is the dimension. Similarly, the value of $v_{C}$ can be computed from the
right state using the Rankine-Hugoniot conditions, which relate the quantities
on either side of the shock. In this case, the expression is

$\displaystyle
v_{\text{shock}}(\mathcal{P}_{C})=v_{R}+2c_{R}\left(1-\frac{\mathcal{P}_{C}}{\mathcal{P}_{R}}\right)\left[2\gamma\left(\gamma-1+\frac{\mathcal{P}_{C}(\gamma+1)}{\mathcal{P}_{R}}\right)\right]^{-\frac{1}{2}}.$
(A.3)

A suitable value of $\mathcal{P}_{C}$ can thus be found by solving
$v_{\text{rare}}(\mathcal{P}_{C})=v_{\text{shock}}(\mathcal{P}_{C})$. This
determines the value of $v_{C}$ also. Having established these, the densities
on either side of the contact discontinuity are found as

$\displaystyle\rho_{C,L}=\rho_{L}\left(\frac{\mathcal{P}_{C}}{\mathcal{P}_{L}}\right)^{\frac{1}{\gamma}},\quad\rho_{C,R}=\rho_{R}\left(\frac{\mathcal{P}_{C}}{\mathcal{P}_{R}}+\frac{\gamma-1}{\gamma+1}\right)\left(\frac{\mathcal{P}_{C}(\gamma-1)}{\mathcal{P}_{R}(\gamma+1)}+1\right)^{-1}.$
(A.4)

A detailed derivation of these relations can be found in [57], for instance.

We now proceed to establish the spatial structure of the solution. The
rarefaction, contact discontinuity, and shock wave all radiate from the
initial point of discontinuity. The head and the tail of the rarefaction
radiate at speeds $v_{L}-c_{L}$ and $v_{C}-c_{C}$ respectively. The contact
discontinuity travels at speed $v_{C}$. Lastly, the shock propagates at speed
$s=\frac{\rho_{R}v_{R}-\rho_{C,R}v_{C}}{\rho_{R}-\rho_{C,R}}.$ These four
speeds and their associated trajectories partition space into five intervals.
The solution is constant inside each interval, except inside the rarefaction
wave, where it takes the form:

$\begin{array}[]{ll}v_{\text{rare}}(\xi)=\left[(\gamma-1)v_{L}+2(c_{l}+\xi)\right]\left(\gamma+1\right)^{-1},&c_{\text{rare}}(\xi)=v_{\text{rare}}(\xi)-\xi,\\\
\rho_{\text{rare}}(\xi)=\left[\rho_{L}^{\gamma}c_{\text{rare}}^{2}(\xi)\left(\gamma\mathcal{P}_{L}\right)^{-1}\right]^{\frac{1}{\gamma-1}},&\mathcal{P}_{\text{rare}}(\xi)=\rho_{\text{rare}}^{\gamma}(\xi)\mathcal{P}_{L}\rho_{L}^{-\gamma},\end{array}$
(A.5)

where $\xi=x/t$. In summary, at time $t$ and position $x$, the solution takes
the form

$\left\\{\begin{array}[]{llll}\displaystyle\rho=\rho_{L},&v=v_{L},&\mathcal{P}=\mathcal{P}_{L},&\text{if
}\xi\leq v_{L}-c_{L};\\\
\displaystyle\rho=\rho_{\text{rare}}(\xi),&v=v_{\text{rare}}(\xi),&\mathcal{P}=\mathcal{P}_{\text{rare}}(\xi),&\text{if
}v_{L}-c_{L}\leq\xi\leq v_{C}-c_{C};\\\
\displaystyle\rho=\rho_{C,L},&v=v_{C},&\mathcal{P}=\mathcal{P}_{C},&\text{if
}v_{C}-c_{C}\leq\xi\leq v_{C};\\\
\displaystyle\rho=\rho_{C,R},&v=v_{C},&\mathcal{P}=\mathcal{P}_{C},&\text{if
}v_{C}\leq\xi\leq s;\\\
\displaystyle\rho=\rho_{R},&v=v_{R},&\mathcal{P}=\mathcal{P}_{R},&\text{if
}s\leq\xi.\end{array}\right.$ (A.6)

A two fluid Sod Tube problem can be constructed by considering two initially
separate fluids, in contact at the point $x=0$; the left fluid with state
($\rho_{L}$, $v_{L}$, $\mathcal{P}_{L}$), and the right fluid with state
($\rho_{R}$, $v_{R}$, $\mathcal{P}_{R}$). This datum generates a solution
profile which is identical to the one described above, and the two fluids
simply remain unmixed, separated by the contact discontinuity [44]. This is
true even if the fluids have different adiabatic exponents, in which case the
solution is calculated by using $\gamma_{L}$ and $\gamma_{R}$ respectively
left and right of the contact discontinuity.

## References

  * [1] A. Alaia, S. Pieraccini, and G. Puppo.  Velocity Discretization in Numerical Schemes for BGK Equations.  In S. Benzoni-Gavage and D. Serre, editors, Hyperbolic Problems: Theory, Numerics, Applications, pages 857–864. Springer Berlin Heidelberg, Berlin, Heidelberg, 2008.
  * [2] R. Aure and J. W. Jacobs.  Particle image velocimetry study of the shock-induced single mode Richtmyer–Meshkov instability.  Shock Waves, 18(3):161–167, Aug. 2008.
  * [3] F. Bernard, A. Iollo, and G. Puppo.  BGK Polyatomic Model for Rarefied Flows.  J. Sci. Comput., 78(3):1893–1916, Mar. 2019.
  * [4] P. L. Bhatnagar, E. P. Gross, and M. Krook.  A Model for Collision Processes in Gases. I. Small Amplitude Processes in Charged and Neutral One-Component Systems.  Phys. Rev., 94(3):511–525, May 1954.
  * [5] G. A. Bird.  Molecular Gas Dynamics and the Direct Simulation of Gas Flows.  Number 42 in Oxford Engineering Science Series. Clarendon Press ; Oxford University Press, Oxford : New York, 1994.
  * [6] M. Bisi, A. V. Bobylev, M. Groppi, and G. Spiga.  Hydrodynamic equations from a BGK model for inert gas mixtures.  In 31st International Symposium on Rarefied Gas Dynamics, Glasgow, UK, 2019.
  * [7] A. Bondesan, L. Boudin, and B. Grec.  A numerical scheme for a kinetic model for mixtures in the diffusive limit using the moment method.  Numer. Meth. Part. Differ. Equat., 35(3):1184–1205, 2019.
  * [8] M. Briant.  étude des gaz multi-espèces proches de l’équilibre.  Séminaire L. Schwartz, pages 1–17, 2015.
  * [9] M. Briant.  Stability of global equilibrium for the multi-species Boltzmann equation in $L^{\infty}$ settings.  DCDS-A, 36:6669–6688, Oct. 2016.
  * [10] M. Briant and E. S. Daus.  The Boltzmann Equation for a Multi-species Mixture Close to Global Equilibrium.  Arch. Ration. Mech. Anal., 222(3):1367–1443, Dec. 2016.
  * [11] S. Brull and C. Prigent.  Local Discrete Velocity Grids for Multi-Species Rarefied Flow Simulations.  Commun. Comput. Phys., 28(4):1274–1304, June 2020.
  * [12] C. Buet.  A discrete-velocity scheme for the Boltzmann operator of rarefied gas dynamics.  Transport Theor. Stat., 25(1):33–60, Jan. 1996.
  * [13] Z. Cai and R. Li.  Numerical Regularized Moment Method of Arbitrary Order for Boltzmann-BGK Equation.  SIAM J. Sci. Comput., 32(5):2875–2907, Jan. 2010.
  * [14] C. Cercignani.  The Boltzmann Equation and Its Applications, volume 67 of Applied Mathematical Sciences.  Springer New York, New York, NY, 1988.
  * [15] C. Cercignani.  Are there more than five linearly-independent collision invariants for the Boltzmann equation?  J. Stat. Phys., 58(5-6):817–823, Mar. 1990.
  * [16] C. Cercignani, R. Illner, and M. Pulvirenti.  The Mathematical Theory of Dilute Gases, volume 106 of Applied Mathematical Sciences.  Springer New York, New York, NY, 1994.
  * [17] S. Chapman and T. G. Cowling.  The Mathematical Theory of Non-Uniform Gases: An Account of the Kinetic Theory of Viscosity, Thermal Conduction, and Diffusion in Gases.  Cambridge Mathematical Library. Cambridge University Press, Cambridge ; New York, 3rd ed edition, 1990.
  * [18] N. Crouseilles, M. Mehrenberger, and E. Sonnendrücker.  Conservative semi-Lagrangian schemes for Vlasov equations.  J. Comput. Phys., 229(6):1927–1953, Mar. 2010.
  * [19] N. Crouseilles, T. Respaud, and E. Sonnendrücker.  A forward semi-Lagrangian method for the numerical solution of the Vlasov equation.  Comput. Phys. Commun., 180(10):1730–1745, Oct. 2009.
  * [20] E. S. Daus, A. Jüngel, C. Mouhot, and N. Zamponi.  Hypocoercivity for a Linearized Multispecies Boltzmann System.  SIAM J. Math. Anal., 48(1):538–568, Jan. 2016.
  * [21] R. Dillman and J. Corliss.  Overview of the Mars sample return Earth entry vehicle.  In Sixth International Planetary Probe Workshop, 2008.
  * [22] G. Dimarco, R. Loubère, J. Narski, and T. Rey.  An efficient numerical method for solving the Boltzmann equation in multidimensions.  J. Comput. Phys., 353:46–81, Jan. 2018.
  * [23] G. Dimarco and L. Pareschi.  Numerical methods for kinetic equations.  Acta Numer., 23:369–520, May 2014.
  * [24] R. S. Ellis and M. A. Pinsky.  The first and second fluid approximations to the linearized Boltzmann equation.  J. Math. Pures Appl., 54(9):125–156, 1975.
  * [25] B. Engquist, X. Li, W. Ren, E. Vanden-Eijnden, and E. Weinan.  Heterogeneous multiscale methods: A review.  Commun. Comput. Phys., 2(3):367–450, 2007.
  * [26] R. Eymard, T. Gallouët, and R. Herbin.  Finite volume methods.  Handbook of numerical analysis, 7:713–1018, 2000.
  * [27] F. Filbet and T. Rey.  A rescaling velocity method for dissipative kinetic equations. Applications to granular media.  J. Comput. Phys., 248:177–199, Sept. 2013.
  * [28] I. Gallagher, L. Saint-Raymond, and B. Texier.  From Newton to Boltzmann: Hard Spheres and Short-Range Potentials.  Zurich Lectures in Advanced Mathematics. European Mathematical Society, Zürich, Switzerland, 2013.
  * [29] C. W. Gear, J. M. Hyman, P. G. Kevrekidid, I. G. Kevrekidis, O. Runborg, and C. Theodoropoulos.  Equation-Free, Coarse-Grained Multiscale Computation: Enabling Mocroscopic Simulators to Perform System-Level Analysis.  Commun. Math. Sci., 1(4):715–762, Dec. 2003.
  * [30] C. W. Gear and I. G. Kevrekidis.  Projective Methods for Stiff Differential Equations: Problems with Gaps in Their Eigenvalue Spectrum.  SIAM J. Sci. Comput., 24(4):1091–1106, Jan. 2003.
  * [31] C. W. Gear and I. G. Kevrekidis.  Telescopic projective methods for parabolic differential equations.  J. Comput. Phys., 187(1):95–109, May 2003.
  * [32] P. Gervais.  Spectral study of the linearized Boltzmann operator in $L^{2}$ spaces with polynomial and Gaussian weights.  arXiv:2010.10339 [math], Oct. 2020.
  * [33] J. R. Haack, C. D. Hauck, and M. S. Murillo.  A Conservative, Entropic Multispecies BGK Model.  J. Stat. Phys., 168(4):826–856, Aug. 2017.
  * [34] J. Hu, S. Jin, and Q. Li.  Asymptotic-Preserving Schemes for Multiscale Hyperbolic and Kinetic Equations.  In Handbook of Numerical Analysis, volume 18, pages 103–129. Elsevier, 2017.
  * [35] S. Jin.  Efficient Asymptotic-Preserving (AP) Schemes For Some Multiscale Kinetic Equations.  SIAM J. Sci. Comput., 21(2):441–454, Jan. 1999.
  * [36] S. Jin.  Asymptotic preserving (AP) schemes for multiscale kinetic and hyperbolic equations: A review.  Riv. Mat. Univ. Parma, 3:177–216, 2012.
  * [37] T. Kataoka, M. Tsutahara, K. Ogawa, Y. Yamamoto, M. Shoji, and Y. Sakai.  Knudsen Pump and Its Possibility of Application to Satellite Control, 2004.
  * [38] J. Koellermeier and G. Samaey.  Projective Integration for Moment Models of the BGK Equation.  In V. V. Krzhizhanovskaya, G. Závodszky, M. H. Lees, J. J. Dongarra, P. M. A. Sloot, S. Brissos, and J. Teixeira, editors, Computational Science – ICCS 2020, volume 12142, pages 321–333. Springer International Publishing, Cham, 2020.
  * [39] P. Lafitte, A. Lejon, and G. Samaey.  A High-Order Asymptotic-Preserving Scheme for Kinetic Equations Using Projective Integration.  SIAM J. Numer. Anal., 54(1):1–33, Jan. 2016.
  * [40] P. Lafitte, W. Melis, and G. Samaey.  A high-order relaxation method with projective integration for solving nonlinear systems of hyperbolic conservation laws.  J. Comput. Phys., 340:1–25, July 2017.
  * [41] P. Lafitte and G. Samaey.  Asymptotic-preserving Projective Integration Schemes for Kinetic Equations in the Diffusion Limit.  SIAM J. Sci. Comput., 34(2):A579–A602, Jan. 2012.
  * [42] S. L. Lee and C. W. Gear.  Second-order accurate projective integrators for multiscale problems.  J. Comput. Appl. Math., 201(1):258–274, Apr. 2007.
  * [43] R. J. LeVeque.  Numerical Methods for Conservation Laws, volume 57.  Birkhäuser Basel, Basel, 1990.
  * [44] R. J. LeVeque.  Finite Volume Methods for Hyperbolic Problems.  Cambridge University Press, Cambridge, 2002.
  * [45] X.-D. Liu, S. Osher, and T. Chan.  Weighted Essentially Non-oscillatory Schemes.  J. Comput. Phys., 115(1):200–212, Nov. 1994.
  * [46] C. P. McNally, W. Lyra, and J.-C. Passy.  A Well-Posed Kelvin-Helmholtz Instability Test and Comparison.  ApJS, 201(2):18, Aug. 2012.
  * [47] W. Melis, T. Rey, and G. Samaey.  Projective Integration for Nonlinear BGK Kinetic Equations.  In C. Cancès and P. Omnes, editors, Finite Volumes for Complex Applications VIII - Hyperbolic, Elliptic and Parabolic Problems, volume 200, pages 145–153. Springer International Publishing, Cham, 2017.
  * [48] W. Melis, T. Rey, and G. Samaey.  Projective and telescopic projective integration for the nonlinear BGK and Boltzmann equations.  SMAI J. Comput. Math., 5:53–88, 2019.
  * [49] W. Melis and G. Samaey.  Telescopic Projective Integration for Linear Kinetic Equations with Multiple Relaxation Times.  J. Sci. Comput., 76(2):697–726, Aug. 2018.
  * [50] L. Mieussens.  Discrete Velocity Model and Implicit Scheme for the BGK Equation of Rarefied Gas Dynamics.  Math. Models Methods Appl. Sci., 10(08):1121–1149, Nov. 2000.
  * [51] T. Rey.  A Spectral Study of the Linearized Boltzmann Equation for Diffusively Excited Granular Media.  arXiv:1310.7234, Oct. 2013.
  * [52] R. Rico-Mart_́inez, C. Gear, and I. G. Kevrekidis.  Coarse projective kMC integration: Forward/reverse initial and boundary value problems.  J. Comput. Phys., 196(2):474–489, May 2004.
  * [53] F. Rogier and J. Schneider.  A direct method for solving the Boltzmann equation.  Transport Theor. Stat., 23(1-3):313–338, Jan. 1994.
  * [54] G. N. Sashi Kumar and N. K. Maheshwari.  Viscous multi-species lattice Boltzmann solver for simulating shock-wave structure.  Comput. Fluids, 203:104539, May 2020.
  * [55] G. A. Sod.  A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws.  J. Comput. Phys., 27(1):1–31, Apr. 1978.
  * [56] H. Struchtrup.  Macroscopic Transport Equations for Rarefied Gas Flows: Approximation Methods in Kinetic Theory.  Interaction of Mechanics and Mathematics. Springer-Verlag, Berlin Heidelberg, 2005.
  * [57] E. F. Toro.  Riemann Solvers and Numerical Methods for Fluid Dynamics.  Springer Berlin Heidelberg, Berlin, Heidelberg, 1999.
  * [58] M. Torrilhon.  Two-Dimensional Bulk Microflow Simulations Based on Regularized Grad’s 13-Moment Equations.  Multiscale Model. Simul., 5(3):695–728, Jan. 2006.
  * [59] C. Villani.  A Review of Mathematical Topics in Collisional Kinetic Theory.  In Handbook of Mathematical Fluid Dynamics, volume 1, pages 71–74. Elsevier, 2002.
  * [60] J. T. Waddell, C. E. Niederhaus, and J. W. Jacobs.  Experimental study of Rayleigh–Taylor instability: Low Atwood number liquid systems with single-mode initial perturbations.  Physics of Fluids, 13(5):1263–1273, May 2001.
  * [61] L. Wu, J. Zhang, J. M. Reese, and Y. Zhang.  A fast spectral method for the Boltzmann equation for monatomic gas mixtures.  J. Comput. Phys., 298:602–621, Oct. 2015.

  *[BGK]: Bhatnagar-Gross-Krook
  *[AP]: asymptotic preserving
  *[PInt]: projective integration
  *[TPInt]: telescopic projective integration