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

HAL
open science

# Formal derivation of an isentropic two-phase flow model from the multi-species Boltzmann equation

Gabriella Puppo, Thomas Rey, Tommaso Tenna

► To cite this version:

Gabriella Puppo, Thomas Rey, Tommaso Tenna. Formal derivation of an isentropic two-phase flow model from the multi-species Boltzmann equation. 2025. hal-05115202

HAL Id: hal-05115202

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

Preprint submitted on 16 Jun 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
FORMAL DERIVATION OF AN ISENTROPIC TWO-PHASE
FLOW MODEL FROM THE MULTI-SPECIES BOLTZMANN
EQUATION

GABRIELLA PUPPO, THOMAS REY, AND TOMMASO TENNA

ABSTRACT. Starting from the multi-species Boltzmann equation for a gas mixture, we propose the formal derivation of the isentropic two-phase flow model introduced in [48]. We examine the asymptotic limit as the Knudsen numbers approach zero, in a regime characterized by resonant intra-species collisions, where interactions between particles of the same species dominate. This specific regime leads to a multi-velocity and multi-pressure hydrodynamic model, enabling the explicit computation of the coefficients for the two-phase macroscopic model. Our derivation also accounts for the inclusion of the evolution of the volume fraction, which is a key variable in many macroscopic multiphase models.

2020 MATHEMATICS SUBJECT CLASSIFICATION: 82B40, 76T10, 76P05.

# 1. INTRODUCTION

Multiphase flow models are widely diffused for the description of fluids with dif-
ferent physical or thermodynamic properties, usually separated by interfaces which
are assumed to be diffuse. In these diffuse interface models, the fluids have different
macroscopic properties, but are allowed to mix across the interface zone. One of the
first models of this type was introduced by Baer and Nunziato in [5] and it is still the
basis for the development of new strategies to describe compressible two-phase flows.
After this pioneering work, many other models were proposed to approximate multi-
phase fluids, introducing some slight modifications to the original system regarding
closure relations and relaxation parameters, see [34, 35, 45, 50, 51]. Their derivation
is mainly based on phenomenological assumptions, which affect the validity of the
model in more general situations. An alternative approach to the Baer-Nunziato
model is the hyperbolic multi-phase model introduced by Romenski, Drikakis and
Toro in [47, 48], based on the theory of symmetric hyperbolic thermodynamically
compatible systems [38]. Despite the fact that this model has been validated from a
numerical point of view (see [44] for its approximation in all Mach number regimes
and [52] for a complete analysis on the Riemann problem in the barotropic case), the
derivation has an intrinsic phenomenological nature and it contains some coefficients
which cannot be explicitly determined. The purpose of the present paper is a formal
derivation of macroscopic equations for multiphase flows from Boltzmann equations
for mixtures of rarefied gases.

The behavior of rarefied gases is usually modeled through kinetic theory, in or-
der to reproduce non-equilibrium phenomena. The seminal Boltzmann Equation

<footnote>(Gabriella Puppo) DIPARTIMENTO DI MATEMATICA, SAPIENZA UNIVERSITÀ DI ROMA, P.LE ALDO MORO 5, 00185 ROME, ITALY</footnote>

<footnote>(Thomas Rey) UNIVERSITÉ CÔTE D'AZUR, CNRS, LJAD, PARC VALROSE, F-06108 NICE, FRANCE</footnote>

<footnote>(Tommaso Tenna) UNIVERSITÉ CÔTE D'AZUR, CNRS, LJAD, PARC VALROSE, F-06108 NICE, FRANCE - DIPARTIMENTO DI MATEMATICA, SAPIENZA UNIVERSITÀ DI ROMA, P.LE ALDO MORO 5, 00185 ROME, ITALY</footnote>

<footnote>Key words and phrases. Hydrodynamic limit; Boltzmann equation; Multi-species mixture; Two-phase flow; Multi-velocity models; Multi-pressure models.</footnote>
represents the cornerstone to describe the dynamics of a gas which is not in ther-
modynamical equilibrium, composed by a very large number of identical particles
[26, 27].

The necessity of studying systems involving gases composed of molecules of dif-
ferent types leads to kinetic models for mixtures. Indeed, realistic kinetic approx-
imations of such situations can be performed only through multi-species models,
in which interactions and exchanges between different components must be proper-
ly described. Whereas the mathematical properties of the Boltzmann equation
for a single-species gas are well known, the analytical study of multi-species kinetic
models is still not complete. Most of the literature concerns theory for Bhatnagar-
Gross-Krook (BGK) models for gas mixtures [3, 14, 17, 41, 43]. To the best of
our knowledge, the first derivation of the Boltzmann equation for binary gas mix-
tures is due to Chapman and Cowling [29]. More recently, Briant and coauthors
[21, 22] have investigated the existence, uniqueness, positivity and exponential trend
to equilibrium for the linearized version of the full non-linear Boltzmann equation
in a perturbative isotropic $L_v^1 L_x^\infty$ setting, namely $L^1$ in velocity space and $L^\infty$ in
the physical space. A spectral gap for the linearized operator around global equi-
librium has been proven in [22], then extended in [15] to the case of the linearized
operator around non-equilibrium states and in [16], providing the explicit spectral
gap for the linearized operator modeling reactive mixtures. In another recent paper,
Gamba and Pavić-Čolić [37] have given an existence and uniqueness result for the
spatially homogeneous multi-species Boltzmann equations with binary interactions
in the Banach space represented by the $L^1$ space weighted polynomially in velocity.
Properties of such solutions have been later analyzed in [31]. Further investigations
on the mixing of monatomic and polyatomic gases can be found in [2].

The formal hydrodynamic limit for kinetic mixtures has been investigated for different regimes, both for the BGK model [10], for a mixed BGK-Boltzmann model [9, 11] and for the Boltzmann model [6, 8, 12, 13, 18, 19]. Let us mention also the work of Briant and Grec [23] on the rigorous derivation of the Fick cross-diffusion system from the multi-species Boltzmann equation.

All hydrodynamic limits of the Boltzmann equation correspond to situations where the frequency of collisions is rapidly increasing. This means that the collision mechanism governed by the collision operator holds on a time scale which is much smaller than the observation time scale and we have an almost instantaneous convergence of the distribution function to thermodynamic equilibrium. We refer to the book of Saint-Raymond for a complete review of hydrodynamic limits of the Boltzmann equation in the single-species framework [49]. In any case, the fluid limit provides a set of equations for the macroscopic fields, obtained as moments of the kinetic distribution function. This approach leads to a hierarchy of equations coupling lower order moments to those of higher order, which could be truncated to obtain the desired order of approximation [10]. In a very recent paper [30] the author has derived a multi-velocity and multi-temperature macroscopic model without phase separation obtained within extended thermodynamics from Boltzmann-like equations, for mixtures of monatomic and polyatomic gases.

The aim of this paper is to propose a formal derivation of the model proposed by Romenski, Drikakis and Toro [48] using kinetic theory. In particular, after a brief introduction to the multi-species Boltzmann equation, we discuss a hydrodynamic limit for this model in the regime of resonant intra-species collisions [36], namely when the interactions between particles of the same species represent the fast process of the dynamics. This regime allows fast convergence towards local “equilibria”, in which each species reaches its own velocity, whereas inter-species interactions are characterized by a much slower process. The novelty of the work lies in the justification of the different scalings, in the introduction of the notion of “volume fraction”
starting from the kinetic formulation and in the use of the full Boltzmann collision operator. Moreover, deriving a constitutive law from the microscopic dynamics at the interface of two gases, we are able to recover the original model and to propose explicit expressions for the coefficients appearing in the macroscopic model.

## 2. THE MULTI-SPECIES BOLTZMANN EQUATION

The derivation of the multi-species Boltzmann equation is mostly formal, unlike that of the single-species case. Indeed, it is possible to consider some modeling assumptions, extending the microscopic dynamics of the single-species case to gas mixtures, in which different species are taken into account.

Let us consider a mixture of $M$ species, each described by a distribution function $f_p = f_p(t, x, v)$, where the mass of a molecule of component $p$ is denoted by $m_p$. Given $x \in \Omega \in \mathbb{R}^d$, $v \in \mathbb{R}^d$, the distribution functions $f_p$ evolve according to the so-called multi-species Boltzmann equation, which, in absence of external forces, takes the following form [22]

$$ (2.1) \qquad \begin{cases} \partial_t f_p + v \cdot \nabla_x f_p = \frac{1}{\epsilon_p} Q_p(f), & p = 1, \dots, M, \\ Q_p(f) = \sum_{q=1}^{M} Q^{pq}(f_p, f_q), \\ f_p(0, x, v) = f_{p,0}(x, v), \end{cases} $$

where $f = (f_1, \dots, f_M)$. The parameter $\epsilon_p > 0$ is the Knudsen number, that is the ratio between the mean free path of particles before each collision and the length scale of observation. Therefore, $\epsilon_p$ governs the frequency of collisions for the species $p$: taking a small Knudsen number corresponds to more frequent collisions. This particular scaling of the Boltzmann equation is usually known as hyperbolic scaling. The time evolution of the distribution function due to particle collisions in the Boltzmann equation is given by a sum of collision operators $Q^{pq}$, referred to as self-collision operators (or intra-species collision operators) if $p = q$ and cross-collision operators (or inter-species collision operators) if $p \neq q$. The microscopic dynamics is fundamental to understand the formulation of the Boltzmann equation and to derive an explicit formulation for the collision operator in (2.1). We need the following reasonable assumptions:

(i) The particles interact via binary collisions. This means that collisions between 3 or more particles can be neglected.

(ii) These binary collisions are localized in space and time.

(iii) Collisions preserve mass for each species, total momentum and total kinetic energy, namely for two particles of species $p$ and $q$ involved in the collision it holds

$$ \begin{cases} m_p v' + m_q v'_* = m_p v + m_q v_*, \\ m_p |v'|^2 + m_q |v'_*|^2 = m_p |v|^2 + m_q |v_*|^2. \end{cases} $$

We avoid any superscripts $v^p$ in the notations and we simply refer to the pre-collisional velocities as $(v, v_*)$ and to the post-collisional ones as $(v', v'_*)$.

In particular, the conservation of total momentum and total energy yields the following parametrizations for the post-collisional velocities $(v', v'_*)$ in terms of the pre-collisional ones $(v, v_*)$:

* The $\omega$-representation

$$ (2.2) \qquad \begin{aligned} v' &= v - \frac{2m_q}{m_p+m_q}((v-v_*)\cdot\omega)\omega, \\ v'_* &= v_* + \frac{2m_p}{m_p+m_q}((v-v_*)\cdot\omega)\omega, \end{aligned} $$
<figure><img src="image_3.png" /><figcaption>FIGURE 1. Geometrical configuration of elastic collisions.</figcaption></figure>

where $\omega \in \mathbb{S}^{d-1}$.

• The $\sigma$-representation

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

where $\sigma \in \mathbb{S}^{d-1}$.

The $\sigma$-representation is related to the $\omega$-one through the following relation

$$ (2.4) \qquad ((v - v_*) \cdot \omega) \omega = \frac{1}{2}(v - v_* - |v - v_*|\sigma). $$

Defining $\mathcal{V}$ as the velocity of the center of mass of the particles $p$ and $q$:

$$ (2.5) \qquad \mathcal{V} = \frac{m_p}{m_p + m_q} v + \frac{m_q}{m_p + m_q} v_*, $$

the laws of elastic collisions gives

$$ (2.6) \qquad v = \mathcal{V} + \frac{m_q}{m_p + m_q} (v - v_*) \qquad v_* = \mathcal{V} - \frac{m_p}{m_p + m_q} (v - v_*). $$

and we can also rewrite the post-collisional velocities as

$$ (2.7) \qquad v' = \mathcal{V} + \frac{m_q}{m_p + m_q} (v' - v_*) \qquad v'_* = \mathcal{V} - \frac{m_p}{m_p + m_q} (v' - v'_*). $$

Hence, the points $v, v_*, v'$ and $v'_*$ belong to the plane defined by $\mathcal{V}$ and $\text{Span}(\sigma, v - v_*)$ and we have the geometric configuration in Fig.1. The angle $\theta$ depends on $\sigma$ according to the relation

$$ (2.8) \qquad \cos \theta := \frac{v - v_*}{|v - v_*|} \cdot \sigma. $$
Using the microscopic hypotheses (i)-(ii)-(iii), it is possible to derive the explicit expression for the strong form of the collision operators as

$$ (2.9) \quad Q^{pq}(f_p, f_q) = \int_{\mathbb{R}^d} \int_{S^{d-1}} \mathbf{B}_{pq}(|v - v_*|, \theta) [f_q(v'_*) f_p(v') - f_q(v_*) f_p(v)] d\Omega dv_* $$

The collision kernels $\mathbf{B}_{pq}$ are non-negative functions, which contain all the physical information about particles interactions, we refer to [26, 27] for a detailed analysis on this subject.

In the following, we consider a collision kernel $\mathbf{B}_{pq}$ which depends only on $|v - v_*|$ and $\cos \theta$. For a general spherically symmetrical intermolecular potential $\Phi(r)$, the deflection angle $\theta$ is given by [7, 28]:

$$ (2.10) \quad \theta(|u|) = \pi - 2 \int_{0}^{W_1} \left[ 1 - W^2 - \frac{4\Phi(r)}{m|u|^2} \right]^{-1/2} dW, $$

where $W = d_p \cos(\theta/2)/r$, $W_1$ is the positive root of the term in the brackets, $d_p$ is the molecular diameter and $r$ is the distance between two molecules. For the case of $(k-1)$th inverse-power laws, namely for a potential of the form $\Phi(r) \propto 1/(k-1) r^{k-1}$, the relevant simplification is that $\mathbf{B}_{pq}(|v - v_*|, \theta)$ becomes the product of a function of $\theta$ alone by a fractional power of $|v - v_*|$. Thus, in the case of inverse-power law potentials, it is possible to rewrite the collision kernel as [7]

$$ (2.11) \quad \mathbf{B}_{pq}(|v - v_*|, \theta) = b_{pq}(\theta)|v - v_*|^\gamma, \quad \gamma = \frac{k - (2d - 1)}{k - 1}, $$

where $b_{pq}(\theta)$ is a nonelementary function of $\theta$. In particular, in $d=3$ dimensions of space, we obtain $\gamma = (k-5)/(k-1)$. Thus, the choice $k=5$ corresponds to the particular case in which the collision kernel does not depend on the relative velocity, but only on the deviation angle. Such particles are called Maxwellian molecules and $\mathbf{B}_{pq}(|v - v_*|, \theta) = b_{pq}(\cos \theta)$. The pseudo-Maxwellian case is given by a constant value of $b_\gamma(\cos \theta)$. For numerical purpose a widely diffused model is the variable hard sphere (VHS) model introduced by Bird [7], corresponding to the choice $b_{pq}(\cos \theta) = C_\gamma$, where $C_\gamma > 0$. The case $\gamma = 0$ is the aforementioned pseudo-Maxwellian molecules case, whereas $\gamma = 1$ gives the hard-spheres molecules case.

The hypothesis of constant $\gamma$ for different species could be weakened, for instance by considering a collision kernel (2.11) with different $\gamma_p$ according to the species $p$, as done in [2], but it does not affect the main conclusions of this work.

### 2.1. Weak Form.
In this section we investigate the weak form of the multi-species Boltzmann equation, following the clear derivation presented in [19] (see also [26] for the single-species case). The operators $Q^{pp}$ are the classical Boltzmann collision operators. Starting from the definition of velocities obtained in (2.3) and using the microscopic hypotheses, we derive the weak form in the $\sigma$-representation, considering a $C^\infty$ test function $\psi: \mathbb{R}^d \to \mathbb{R}$. We assume that the cross section $\mathbf{B}_{pq}$ satisfies the micro-reversibility assumption $\mathbf{B}_{pq}(v, v_*, \sigma) = \mathbf{B}_{pq}(v_*, v, \sigma)$. We obtain

$$ (2.12) \quad \int_{\mathbb{R}^d} Q^{pp}(f_p, f_p)(v) \psi(v) dv = \frac{1}{2} \int_{\mathbb{R}^d \times \mathbb{R}^d \times S^{d-1}} f_{p*} f_p (\psi' + \psi'_* - \psi - \psi_*) \mathbf{B}_{pp}(|v - v_*|, \cos \theta) d\sigma dv dv_*, $$

where we have introduced the shorthands $f'_p = f_p(v')$, $f_{p*} = f_p(v_*)$, $f'_{p*} = f_p(v'_*)$, $\psi' = \psi(v')$, $\psi_* = \psi(v_*)$ and $\psi'_* = \psi(v'_*)$.

In an analogous way, it is possible to characterize the bi-species collision operators, under the assumption of a collision kernel independent from the species. Therefore,
considering the definition (2.3) with $p \neq q$, in the same way of the mono-species framework, we obtain the following formulation of $Q^{pq}$ in the $\sigma$-representation

$$ (2.13) \quad \int_{\mathbb{R}^d} Q^{pq}(f_p, f_q)(v)\psi(v) dv = \int_{\mathbb{R}^d \times \mathbb{R}^d \times \mathbb{S}^{d-1}} f_p f_{q*}(\psi' - \psi) \mathbf{B}_{pq}(|v - v_*|, \cos \theta) d\sigma \, dv \, dv_* $$

In addition, if we consider both the bi-species collision operators, we can give the following characterization:

$$ (2.14) \quad \int_{\mathbb{R}^d} Q^{pq}(f_p, f_q)(v)\psi(v) dv + \int_{\mathbb{R}^d} Q^{qp}(f_q, f_p)(v)\phi(v) dv = \int_{\mathbb{R}^d \times \mathbb{R}^d \times \mathbb{S}^{d-1}} f_p f_{q*}(\psi' + \phi_*' - \psi - \phi_*) \mathbf{B}_{pq}(|v - v_*|, \cos \theta) d\sigma \, dv \, dv_* $$

for any $C^\infty$ test functions $\psi, \phi: \mathbb{R}^d \to \mathbb{R}$.

Using the weak form, one can easily observe that the conservation properties of the Boltzmann equation are extended to the multi-species cases. The mass of each species is conserved, since

$$ (2.15) \quad \int_{\mathbb{R}^d} Q_p(f, g)(v) dv = 0, $$

for any functions $f(v)$ and $g(v)$, namely the constant function $v \mapsto 1$ is a collisional invariant for this operator. On the other hand, the functions $v$ and $|v|^2$ are not collision invariants, but it holds

$$ (2.16) \quad \int_{\mathbb{R}^d} v Q^{pq}(f,g)(v) dv = - \int_{\mathbb{R}^d} v Q^{qp}(g,f)(v) dv, $$

and

$$ (2.17) \quad \int_{\mathbb{R}^d} |v|^2 Q^{pq}(f,g)(v) dv = - \int_{\mathbb{R}^d} |v|^2 Q^{qp}(g,f)(v) dv. $$

In particular, the choice of $\psi(v) = m_i v$ and $\phi(v) = m_j v$ in (2.14) leads to the conservation of the total momentum, whereas $\psi(v) = m_i |v|^2$ and $\phi(v) = m_j |v|^2$ leads to the conservation of the total kinetic energy of the mixture.

## 2.2. Entropy dissipation.
The large time behavior of the multi-species Boltzmann equation (2.1) is here briefly described. In this subsection we restrict our analysis to the case of a bi-species mixture for the sake of simplicity, but it can be easily extended to the case of M-species mixtures. In particular, we will consider the space homogeneous setting

$$ (2.18) \quad \left\{ \begin{array}{l} \partial_t f_1 = Q^{11}(f_1, f_1) + Q^{12}(f_1, f_2), \\ \partial_t f_2 = Q^{21}(f_2, f_1) + Q^{22}(f_2, f_2). \end{array} \right. $$

**Definition 2.1 (Multi-species Entropy).** Let $(f_1, f_2)$ be a solution to the homogeneous Boltzmann equation (2.18). The multi-species entropy is described by the so-called Boltzmann's H-functional

$$ (2.19) \quad H = \int_{\mathbb{R}^d} f_1 \log f_1(v) \, dv + \int_{\mathbb{R}^d} f_2 \log f_2(v) \, dv. $$

**Remark 2.1.** This definition can be easily extended to the case of $M \ge 2$ species by defining

$$ (2.20) \quad H = \sum_{i=1}^{M} \int_{\mathbb{R}^d} f_p \log f_p(v) \, dv. $$
Let us now recall the following result on entropy dissipation, proven in [32], here explicitly detailed for a non-reactive mixture of two species.

**Proposition 2.2.** Let us assume strictly positive collision kernels $\mathbf{B}_{pq}$ and let us consider a positive solution $(f_1, f_2)$ to (2.1). Then, the entropy is dissipated, namely

$$ (2.21) \qquad \frac{dH}{dt} \le 0. $$

Moreover, we have $\frac{dH}{dt} = 0$ if and only if the equilibrium $(f_1, f_2)$ is given by a pair of global Maxwellian distributions $(M_1, M_2)$ with common average velocity $\bar{v}$ and temperature $T_{eq}$:

$$ (2.22) \qquad \mathcal{M}_p^{n_p, \bar{v}, T}(t, x, v) = n_p \left( \frac{m_p}{2\pi T_{eq}} \right)^{\frac{d}{2}} \exp \left( -\frac{m_p}{2T_{eq}} |\bar{v} - v|^2 \right), \quad p = 1, 2. $$

*Sketch of the proof.* Let us consider the expression

$$ (2.23) \qquad \frac{dH}{dt} = \int_{\mathbb{R}^d} Q^{11}(f_1, f_1) \log f_1(v) dv + \int_{\mathbb{R}^d} Q^{12}(f_1, f_2) \log f_1(v) dv + \\ \qquad + \int_{\mathbb{R}^d} Q^{21}(f_2, f_1) \log f_2(v) dv + \int_{\mathbb{R}^d} Q^{22}(f_2, f_2) \log f_2(v) dv. $$

We consider the weak formulation of the collision operators (2.12)-(2.14) with $\phi = \log f_1$ and $\psi = \log f_2$:

$$ \begin{aligned} \frac{dH}{dt} &= \frac{1}{2} \int_{\mathbb{R}^d \times \mathbb{R}^d \times \mathbb{S}^{d-1}} f_{1*} f_1 \log \left( \frac{f_1' f_{1*}'}{f_1 f_{1*}} \right) \mathbf{B}_{11} d\sigma \, dv \, dv_* \\ &\quad - \frac{1}{2} \int_{\mathbb{R}^d \times \mathbb{R}^d \times \mathbb{S}^{d-1}} [f_1' f_{2*}' - f_1 f_{2*}] \log \left( \frac{f_1' f_{2*}'}{f_1 f_{2*}} \right) \mathbf{B}_{12} d\sigma \, dv \, dv_* \\ &\quad + \frac{1}{2} \int_{\mathbb{R}^d \times \mathbb{R}^d \times \mathbb{S}^{d-1}} f_{2*} f_2 \log \left( \frac{f_2' f_{2*}'}{f_2 f_{2*}} \right) \mathbf{B}_{22} d\sigma \, dv \, dv_* \\ &= \frac{1}{2} \sum_{p,q \in \{1,2\}} \int_{\mathbb{R}^d \times \mathbb{R}^d \times \mathbb{S}^{d-1}} f_{p*} f_q \left[ \log \left( \frac{f_p' f_{q*}'}{f_p f_{q*}} \right) - \frac{f_p' f_{q*}'}{f_p f_{q*}} + 1 \right] \mathbf{B}_{pq} d\sigma \, dv \, dv_* \\ &\quad + \frac{1}{2} \sum_{p,q \in \{1,2\}} \int_{\mathbb{R}^d \times \mathbb{R}^d \times \mathbb{S}^{d-1}} (f_{q*}' f_p' - f_{q*} f_p) \mathbf{B}_{pq} d\sigma \, dv \, dv_* \le 0. \end{aligned} $$

Indeed, since we are assuming strictly positive cross sections, $\log x - x + 1 \le 0$ and the collisional transformation $(v, v_*) \to (v', v_*')$ is an involution, all terms on the right hand side of the equation above are non-positive and this proves the entropy dissipation.

To prove the second part, we first observe that, since we are considering a sum of non-positive contributions, $\frac{dH}{dt} = 0$ implies that all these terms must be zero. The terms corresponding to intra-species collisions (namely depending only on $f_p$) vanish if there exist $n_p \ge 0$, $\bar{v}_p \in \mathbb{R}^d$ and $T_p > 0$ such that

$$ (2.24) \qquad f_p(v) = n_p \left( \frac{m_p}{2\pi T_p} \right)^{\frac{d}{2}} \exp \left( -\frac{m_p}{2T_p} |\bar{v}_p - v|^2 \right), \quad p = 1, 2. $$

To conclude the proof, in [32] the authors prove that $v_1 = v_2 = \bar{v}$ and $T_1 = T_2 = T_{eq}$. $\square$

**Remark 2.3.** Let us underline that the local Maxwellians $\mathcal{M}_p(t, x, v)$ defined in (2.24) are maximum of the Boltzmann entropy in the mono-species setting, but they are not local equilibria for the gas mixture. Nevertheless, they play a fundamental role in the derivation of hydrodynamic limits for different regimes.
Let us consider the inhomogeneous system with a hyperbolic scaling $\epsilon$, under the assumption that the relaxation velocity is the same for each collision operator, namely

$$ (2.25) \qquad \left\{ \begin{aligned} \partial_t f_1 + v \cdot \nabla_x f_1 &= \frac{1}{\epsilon} Q^{11}(f_1, f_1) + \frac{1}{\epsilon} Q^{12}(f_1, f_2), \\ \partial_t f_2 + v \cdot \nabla_x f_2 &= \frac{1}{\epsilon} Q^{21}(f_2, f_1) + \frac{1}{\epsilon} Q^{22}(f_2, f_2). \end{aligned} \right. $$

In this framework, the equilibrium $(f_1, f_2)$ is reached in the limit $\epsilon \to 0$ and the distribution functions converge towards Maxwellian functions with common average velocity and temperature, as shown in Proposition 2.2.

In mixtures of gases with disparate masses the exchange of momentum between particles of the same species is usually a faster process and the relaxation velocity is not the same for each collision operator. In this case, the inhomogeneous system is scaled as follows

$$ (2.26) \qquad \left\{ \begin{aligned} \partial_t f_1 + v \cdot \nabla_x f_1 &= \frac{1}{\epsilon} Q^{11}(f_1, f_1) + \frac{1}{\kappa} Q^{12}(f_1, f_2), \\ \partial_t f_2 + v \cdot \nabla_x f_2 &= \frac{1}{\kappa} Q^{21}(f_2, f_1) + \frac{1}{\epsilon} Q^{22}(f_2, f_2), \end{aligned} \right. $$

where $\epsilon > 0$ and $\kappa > 0$ may not have the same order of magnitude. Despite the minimum of the entropy (the global equilibrium) is still a pair of global Maxwellian $(M_1, M_2)$, the dynamics involves processes at two different time scales and the convergence no longer occurs on the same timescale. Let us assume that the two species are segregated, as in Fig. 2. For the sake of illustration, we can think of a setting in which droplets (or bubbles) of light gas are separated by a heavy gas. The separation between the two species leads to a configuration in which inter-species collisions can occur only along the separation layers (interface). On the contrary, intra-species collisions can occur inside the volume (in Fig. 2, dashed volumes for one species and white for the other species). This corresponds to the regime $\kappa = O(L/T)$ where $L$ and $T$ are the respective characteristic length and time scales, and $\epsilon \ll \kappa$, in which intra-species collisions are much more frequent than inter-species ones. In this case, for $\epsilon \to 0$, each distribution function $f_p$ converges towards the local Maxwellian (2.24) with local average velocity $\bar{u}_p$ and local temperature $T_p$, before reaching a global equilibrium. The process which leads to common average velocities and common temperatures is much slower and we assume it can be disregarded for finite times.

This regime is physically realistic for mixtures of heavy and light gases, usually referred to as $\epsilon$-mixtures, as described in [10, 36].

## 2.3. Macroscopic Quantities.
Let us define the macroscopic fields for each species as suitable moments of the distribution functions $f_p$. In particular, the number density $n_p$, the mass density $\rho_p$, the average macroscopic velocity $\bar{u}_p$, the partial pressure $P_p$ and the temperature $T_p$ for each species are computed as

$$ (2.27) \qquad n_p = \int_{\mathbb{R}^d} f_p(v) dv, \quad \rho_p \bar{u}_p = \int_{\mathbb{R}^d} m_p v f_p dv, \quad T_p = \frac{m_p}{d n_p} \int_{\mathbb{R}^d} |v - \bar{u}_p|^2 f_p dv, $$

with $\rho_p = m_p n_p$ and $P_p = n_p T_p$. We can also define the total moments for the mixture, given by

$$ (2.28) \qquad \begin{aligned} n &= \sum_p n_p, \quad \rho = \sum_p \rho_p, \quad \rho \bar{u} = \sum_p \rho \bar{u}_p, \\ P &= \sum_p P_p, \quad T = \sum_p \frac{m_p}{d n_p} \int_{\mathbb{R}^d} |v - \bar{u}|^2 f_p dv. \end{aligned} $$
<figure><img src="image_4.png" /><figcaption>FIGURE 2. Diffused interface gases in a control volume.</figcaption></figure>

All these quantities are valid for situations in which the two species are not well separated. When we deal with gases in which species are segregated, intra-species collisions become more frequent and different regimes must be taken into account. For this purpose, we introduce a new quantity $\alpha_p$, describing the volume fraction occupied by the species $p$. Heuristically, let us consider a control volume $V$ and let $\chi_p(x, t)$ be the probability of finding a particle of the species $p$ in $(x, t)$. Then, the volume fraction $\alpha_p(x, t)$ occupied by the phase $p$ is defined as

$$ (2.29) \qquad \alpha_p(x,t) = \frac{1}{|V|} \int_V \chi_p(x,t) dV. $$

In this setting, $\sum_p \alpha_p = 1$ and the situation is depicted in Fig. 2 for a mixture of two species.

**Remark 2.4.** The role of the volume fraction for each species is fundamental for models of separated gases [5, 47]. Indeed, in the passage to the macroscopic level, the system could be described using volume fractions that smoothly vary from one region to the other in the space scale of the control volume $V$.

### 3. HYDRODYNAMIC LIMIT FOR MULTI-SPECIES BOLTZMANN EQUATION

3.1. **The single-velocity case.** In this section we describe the formal derivation of the Euler equations from the multi-species Boltzmann equation, under the assumption that the Knudsen number is approximately the same for all species. In the limit of small Knudsen numbers, that is $\epsilon \to 0$, the frequency of collisions increases to infinity and the solution can be completely described by its local hydrodynamic fields. To do so, we recall the multi-species Boltzmann equation rescaled according to the hyperbolic scaling, where all the collision operators scale with the same rate

$$ (3.1) \qquad \partial_t f_p + v \cdot \nabla_x f_p = \frac{1}{\epsilon} Q_p(f) \qquad p = 1, \dots, M $$

and we apply the Chapman-Enskog expansion as in [10]. Formally, we can consider an expansion of the distribution functions in term of $\epsilon$ as

$$ (3.2) \qquad f_p = f_p^0 + \epsilon f_p^1 + \mathcal{O}(\epsilon^2), $$

obtaining expansions also for the other quantities of interest. In particular the number densities $n_p$, mass velocities $u_p$ and temperatures $T_p$ can be expanded as

$$ (3.3) \qquad n_p = n_p^0 + \epsilon n_p^1 + \mathcal{O}(\epsilon^2), \quad \bar{u}_p = \bar{u}_p^0 + \epsilon \bar{u}_p^1 + \mathcal{O}(\epsilon^2), \quad T_p = T_p^0 + \epsilon T_p^1 + \mathcal{O}(\epsilon^2). $$
Let us consider the weak form of the Boltzmann equation truncated at the zeroth order of expansion

$$ (3.4) \qquad \frac{\partial}{\partial t} \int_{\mathbb{R}^d} \varphi(v)(f_p^0 + \varepsilon f_p^1)dv + \nabla_x \cdot \int_{\mathbb{R}^d} \varphi(v)v(f_p^0 + \varepsilon f_p^1)dv \\ = \frac{1}{\varepsilon} \sum_{q=1}^{N} \int_{\mathbb{R}^d} \varphi(v)Q^{pq}(f_p^0 + \varepsilon f_p^1, f_q^0 + \varepsilon f_q^1)dv, $$

where $\varphi(v)$ is a $C^\infty$ test function.

**Remark 3.1.** In the Chapman-Enskog expansion, $f_p^1$ designates a quantity that depends smoothly on the macroscopic quantities and any finite number of its derivatives with respect to the $(x, v)$-variables (but not on time). In particular, $f_p^1$ satisfies

$$ (3.5) \qquad \int_{\mathbb{R}^d} f_p^1 \begin{pmatrix} 1 \\ v \\ \frac{1}{2}|v|^2 \end{pmatrix} dv = 0. $$

The limit for $\varepsilon \to 0$ in (3.1) implies naturally a vanishing leading order collision operator, thus $f_i^0$ are the Boltzmann equilibria (2.22). The 0-th order term of equation (3.4) provides the multi-species (or multi-component) Euler equations for macroscopic quantities and the first order correction corresponds to the multi-species (or multi-component) Navier-Stokes system [4]. In other words, we can claim that formally the solution to the multi-species Euler equations and the solution to the multi-species Navier-Stokes equations are close to the moments of the solution $f$ to the multi-species Boltzmann equation respectively at order $O(\varepsilon)$ and $O(\varepsilon^2)$.

The multi-species Euler system can be directly recovered with a different approach, considering the collision invariants as test functions in the weak formulation

$$ \frac{\partial}{\partial t} \int_{\mathbb{R}^d} \varphi(v) f_p dv + \nabla_x \cdot \int_{\mathbb{R}^d} \varphi(v) v f_p dv = \sum_{q=1}^{N} \int_{\mathbb{R}^d} \varphi(v) Q^{pq}(f_p, f_q) dv. $$

Let us consider $\varphi(v) \equiv 1$. Since it corresponds to a collision invariant of the Boltzmann operator, the right-hand side vanishes and it leads to the equation

$$ (3.6) \qquad \frac{\partial n_p}{\partial t} + \nabla_x \cdot (n_p \bar{u}) = 0, \quad p = 1, \dots, M, $$

where we have used $\bar{u}_p = \bar{u}$ for all $p$. Multiplying each equation for the respective mass $m_p$ we obtain

$$ (3.7) \qquad \frac{\partial \rho_p}{\partial t} + \nabla_x \cdot (\rho_p \bar{u}) = 0, \quad p = 1, \dots, M. $$

Finally, summing over the index $p$, the continuity equation for the total density holds

$$ (3.8) \qquad \frac{\partial \rho}{\partial t} + \nabla_x \cdot (\rho \bar{u}) = 0. $$

Now consider $\varphi(v) = m_p v$ and sum over the index $p$, in order to obtain

$$ (3.9) \qquad \frac{\partial}{\partial t}(\rho \bar{u}) + \nabla_x \cdot \left[ \sum_{p=0}^{P} \left( \rho_p \bar{u} \otimes \bar{u} + n_p T_{eq} \right) \right] = 0, $$

where the single-species temperature $T_p = T_{eq}$. Let us consider $\varphi(v) = m_p|v|^2$ and sum again over the index $p$. Thus, we obtain

$$ (3.10) \qquad \frac{\partial}{\partial t} \left( \frac{1}{2} \rho |\bar{u}|^2 + \frac{d}{2} n T_{eq} \right) + \nabla_x \cdot \left[ \left( \frac{1}{2} \rho |\bar{u}|^2 + \frac{d+2}{2} n T_{eq} \right) \bar{u} \right] = 0. $$
Thus, in the limit $\epsilon \to 0$, we have formally obtained the so-called multi-species Euler limit:

$$ (3.11) \qquad \left\{ \begin{aligned} & \frac{\partial \rho_p}{\partial t} + \nabla_x \cdot (\rho_p \bar{u}) = 0, \quad p = 1, \dots, M, \\ & \frac{\partial}{\partial t} (\rho \bar{u}) + \nabla_x \cdot \left[ \sum_{p=1}^{N} \rho_p (\bar{u} \otimes \bar{u} + n_p T_{eq}) \right] = 0, \\ & \frac{\partial}{\partial t} \left( \frac{1}{2} \rho |\bar{u}|^2 + \frac{d}{2} n T_{eq} \right) + \nabla_x \cdot \left[ \left( \frac{1}{2} \rho |\bar{u}|^2 + \frac{d+2}{2} n T_{eq} \right) \bar{u} \right] = 0, \end{aligned} \right. $$

with a single velocity and a single temperature for the whole mixture.

3.2. **The multi-velocity case.** Now, we are interested in describing the formal hydrodynamic limit in a different regime, in order to recover the isentropic two-phase flow model proposed in [47], which in the one-dimensional case reads as

$$ \begin{align} (3.12\text{a}) \qquad & \left\{ \begin{aligned} & \frac{\partial \rho}{\partial t} + \frac{\partial}{\partial x}(\rho \bar{u}) = 0, \\ & \frac{\partial}{\partial t}(\alpha_1 \rho) + \frac{\partial}{\partial x}(\alpha_1 \rho \bar{u}) = -\frac{1}{\tau}(P_2 - P_1), \end{aligned} \right. \notag \\ (3.12\text{b}) \qquad & \left\{ \begin{aligned} & \frac{\partial}{\partial t}(\alpha_1 \rho_1) + \frac{\partial}{\partial x}(\alpha_1 \rho_1 \bar{u}_1) = 0, \\ & \frac{\partial}{\partial t}(\rho \bar{u}) + \frac{\partial}{\partial x}(\alpha_1 \rho_1 \bar{u}_1^2 + \alpha_1 P_1 + \alpha_2 \rho_2 \bar{u}_2^2 + \alpha_2 P_2) = 0, \end{aligned} \right. \notag \\ (3.12\text{c}) \qquad & \left\{ \begin{aligned} & \frac{\partial}{\partial t}(\bar{u}_1 - \bar{u}_2) + \frac{\partial}{\partial x} \left( \frac{1}{2} \bar{u}_1^2 + h_1 - \frac{1}{2} \bar{u}_2^2 - h_2 \right) = -\zeta(\bar{u}_1 - \bar{u}_2), \end{aligned} \right. \notag \\ (3.12\text{d}) \qquad & \left\{ \begin{aligned} & \frac{\partial}{\partial t}(\rho \bar{u}) + \frac{\partial}{\partial x}(\alpha_1 \rho_1 \bar{u}_1^2 + \alpha_1 P_1 + \alpha_2 \rho_2 \bar{u}_2^2 + \alpha_2 P_2) = 0, \\ & \frac{\partial}{\partial t}(\bar{u}_1 - \bar{u}_2) + \frac{\partial}{\partial x} \left( \frac{1}{2} \bar{u}_1^2 + h_1 - \frac{1}{2} \bar{u}_2^2 - h_2 \right) = -\zeta(\bar{u}_1 - \bar{u}_2), \end{aligned} \right. \notag \\ (3.12\text{e}) \qquad & \left\{ \begin{aligned} & \frac{\partial}{\partial t}(\rho \bar{u}) = 0, \\ & \frac{\partial}{\partial t}(\alpha_1 \rho_1) = -\frac{1}{\tau}(P_2 - P_1), \end{aligned} \right. \end{align} $$

where $\alpha_p \in (0,1)$ is the volume fraction defined in (2.29), $\tau$ and $\zeta$ are opportune constants, and the partial pressures $P_p$ will be described later. In this framework, the total density $\rho$ is defined as

$$ (3.13) \qquad \rho = \alpha_1 \rho_1 + \alpha_2 \rho_2, $$

and the mixture velocity $\bar{u}$ as

$$ (3.14) \qquad \bar{u} = \kappa_1 \bar{u}_1 + \kappa_2 \bar{u}_2, \qquad \kappa_p = \frac{\alpha_p \rho_p}{\rho}. $$

Compared to the model of two-phase mixtures proposed by Baer and Nunziato in [5], the main advantage of this alternative model is the possibility to write the equations in a conservative form.

Equations (3.12a) and (3.12c) describe the conservation of the total mass and the partial mass, respectively. Equation (3.12b) consists of the balance of the evolution of the volume fraction with respect to a pressure relaxation term. Further, Equation (3.12d) gives the conservation of the total momentum $\rho \bar{u}$ and Equation (3.12e) the balance of the relative velocity $\bar{u}_1 - \bar{u}_2$ with respect to a friction term. Since the system (3.12a)-(3.12e) is isentropic, we set $P_p = P(\rho_p) = c \rho_p^\gamma$. This closure law for the pressure terms implies that an equation for the second moment is not required, differently from (3.11).

To obtain the hyperbolic two-phase flow model as a limit system from a kinetic model, we take into account the case of a gas mixture as in Fig. 2. We recall that in this setting, self-collisions are much more frequent than collisions among particles of different species, namely $\kappa = 1$ and $\epsilon \ll 1$ in (2.26), as described in Section
2.3. For this purpose, we then consider the following rescaling of the multi-species Boltzmann equation

$$ (3.15) \qquad \left\{ \begin{array}{l} \partial_t f_p + v \cdot \nabla_x f_p = \frac{1}{\varepsilon} Q^{pp}(f_p, f_p) + \sum_{\substack{q=1 \\ q \neq p}}^{M} Q^{pq}(f_p, f_q), \\ f_p(0, x, v) = f_{p,0}(x, v). \end{array} \right. $$

The weak form of the Boltzmann equation in this scaling becomes:

$$ (3.16) \qquad \begin{aligned} & \frac{\partial}{\partial t} \int_{\mathbb{R}^d} \varphi(v) f_p dv + \nabla_x \cdot \int_{\mathbb{R}^d} \varphi(v) v f_p dv \\ & = \frac{1}{\varepsilon} \int_{\mathbb{R}^d} \varphi(v) Q^{pp}(f_p, f_p) dv + \sum_{\substack{q=1 \\ q \neq p}}^{M} \int_{\mathbb{R}^d} \varphi(v) Q^{pq}(f_p, f_q) dv. \end{aligned} $$

Next, we average the weak form over the control volume V:

$$ \int_V \chi_p \frac{\partial}{\partial t} \int_{\mathbb{R}^d} \varphi(v) f_p dv \, dV + \int_V \chi_p \nabla_x \cdot \int_{\mathbb{R}^d} \varphi(v) v f_p dv \, dV \\ = \frac{1}{\varepsilon} \int_V \chi_p \int_{\mathbb{R}^d} \varphi(v) Q^{pp}(f_p, f_p) dv \, dV + \sum_{\substack{q=1 \\ q \neq p}}^{M} \int_V \chi_p \int_{\mathbb{R}^d} \varphi(v) Q^{pq}(f_p, f_q) dv \, dV, $$

where $\chi_p$ is the probability to find a particle of species p in $(x, t)$. Each $\chi_p$ is associated to the corresponding volume fraction $\alpha_p$, as introduced in (2.29). Averaging procedure are widely used in the derivation of multi-phase models from the corresponding single-phase formulations (e.g. from multi-component Euler systems), see for instance [33, 1]. In this work, the spatial averaging procedure is instead directly applied to the kinetic equation, following the same rules presented in [33]. Indeed, the averaging procedure over a control volume V commutes with time and space derivatives, and the following calculation rules hold for any function g:

* Gauss rule

$$ \int_V \chi_p \nabla_x g \, dV = \int_V \nabla_x (\chi_p g) \, dV - \int_{\partial V} g_{\text{interface}} \nabla_x \chi_p \, dV, $$

* Leibniz rule

$$ \int_V \chi_p \partial_t g \, dV = \int_V \partial_t (\chi_p g) \, dV - \int_{\partial V} g_{\text{interface}} \partial_t \chi_p \, dV, $$

where $g_{\text{interface}}$ is the value of g at $\partial V$. The boundary $\partial V$ of the domain includes both the external physical boundary and all interfaces separating the phases, where state variables may exhibit discontinuities. Let us then define macroscopic fields as moments of the distribution functions averaged over a unit control volume V (avoiding the scaling factor $|V|^{-1}$ for the sake of simplicity) according to the quantity $\chi_p$, namely

$$ (3.17) \qquad \int_V \int_{\mathbb{R}^d} \chi_p m_p f_p(v) dv \, dV = \alpha_p \rho_p, \qquad \int_V \int_{\mathbb{R}^d} \chi_p m_p v f_p(v) dv \, dV = \alpha_p \rho_p \bar{v}_p, $$

where we recall that $\alpha_p$ stands for the volume fraction. This choice allows us to give a kinetic interpretation of the macroscopic quantities involved in (3.12a)-(3.12e). The probability $\chi_p$ evolves according to the equation

$$ (3.18) \qquad \partial_t \chi_p + \bar{u}_{\text{interface}} \cdot \nabla_x \chi_p = 0, $$

where $\bar{u}_{\text{interface}}$ is the interface velocity between the two phases.
<figure><img src="image_5.png" /><figcaption>FIGURE 3. Interface between two different species.</figcaption></figure>

**Remark 3.2.** The spatial average procedure is fundamental to obtain equations for different phases, describing “non-local” processes which happen at the scale of the control volume.

In order to derive the equation for the fraction $\alpha_1$, we have to determine an expression for the transfer of volume in the dynamics. We will show that the volume transfer is caused by the difference in pressures between the two gases, which in the literature is usually assumed as a phenomenological hypothesis [5, 50, 35].

Let us consider a zoom of the setting of Fig. 2 on the thin layer which separates the volume occupied by one species from the volume occupied by the other, as shown in Fig. 3. We now focus on the displacement of the interface separating the two species. Following the idea of Bresch, Burtea and Lagoutière in [20], at this scale the volume fraction variation can be reinterpreted as a one-dimensional shift orthogonal to the separating interface. Denoting with $x(t)$ the average distance of particles from the interface of separation between the two species, let us consider $\eta_1(t) := x_1(t) - x(t)$ and $\eta_2(t) := x(t) - x_2(t)$ where $x_1$ and $x_2$ are the local centers of mass of the species 1 and 2 respectively. We indeed perform a spatial average over a thin layer which is much smaller than $V$, so we refer to the “centers of mass” of the molecules within this tiny volume.

Assuming that $\eta_p \ll \text{diameter}(V)$, let us define

$$ (3.19) \qquad \alpha_1(t) = \frac{\eta_1(t)}{\eta_1(t) + \eta_2(t)}. $$

Let us introduce the material derivative as

$$ (3.20) \qquad \frac{D}{Dt} := \frac{\partial}{\partial t} + \bar{u} \cdot \nabla_x, $$

where we recall that $\bar{u}$ is the mixture velocity. Since we are considering segregated gases, the mixture velocity satisfies $\bar{u} = \bar{u}_p$ in the region occupied only by species p. In the reference frame of the interface, the material derivative of the volume fraction $\alpha_1$ is given by

$$ (3.21) \qquad \frac{D}{Dt} \alpha_1(t) = \frac{D_t \eta_1(\eta_1 + \eta_2) - \eta_1(D_t \eta_1 + D_t \eta_2)}{(\eta_1 + \eta_2)^2}, $$
where we have omitted the explicit dependence of $\eta_p$ on time. As we considered a small shift with respect to the interface, we have

$$ (3.22) \qquad D_t \eta_p(t) = \bar{u}_p(t, x_p(t)) - \bar{u}(t, x(t)). $$

Therefore, the expression of the derivative of the volume fraction becomes

$$ (3.23) \qquad \frac{D}{Dt} \alpha_1(t) = \frac{(\bar{u}_1(t, x_1(t)) - \bar{u}(t, x(t)))(\eta_1 + \eta_2) - \eta_1(\bar{u}_1(t, x_1(t)) - \bar{u}_2(t, x_2(t)))}{(\eta_1 + \eta_2)^2} \\ = \frac{\eta_1(\bar{u}_2(t, x_2(t)) - \bar{u}(t, x(t))) + \eta_2(\bar{u}_1(t, x_1(t)) - \bar{u}(t, x(t)))}{(\eta_1 + \eta_2)^2}. $$

Let us assume that across the diffuse interface there is a “boundary layer” really separating the two gases, which provides the existence of a continuous stress tensor at the interface: namely given a coefficient $\lambda$ representing the kinematic viscosity, the quantity

$$ (3.24) \qquad P(t, x(t)) - \sqrt{\varepsilon} \lambda \partial_x \bar{u}(t, x(t)) $$

is continuous at the interface. The derivation of (3.24) at the boundary will be detailed in the separate Section 3.2.1 for the reader's convenience.

Approximating the derivative at first order with respect to $\eta_1$ and $\eta_2$ respectively, we have

$$ (3.25) \qquad P(t, x_1(t)) - \sqrt{\varepsilon} \lambda \frac{\bar{u}_1(t, x_1(t)) - \bar{u}(t, x(t))}{\eta_1(t)} \\ = P(t, x_2(t)) - \sqrt{\varepsilon} \lambda \frac{\bar{u}(t, x(t)) - \bar{u}_2(t, x_2(t))}{\eta_2(t)} + \mathcal{O}(\eta_1 + \eta_2). $$

Finally, substituting this expression in (3.23) and considering the limit for $\eta_1 + \eta_2 \to 0$, we obtain

$$ (3.26) \qquad P_1 - P_2 = \lambda \sqrt{\varepsilon} \frac{(\eta_1 + \eta_2)^2}{\eta_1 \eta_2} D_t \alpha_1(t). $$

Now, let us compute

$$ (3.27) \qquad \frac{\partial}{\partial t}(\alpha_1 \rho) + \nabla_x (\alpha_1 \rho \bar{u}) = \alpha_1 \frac{\partial \rho}{\partial t} + \rho \frac{\partial \alpha_1}{\partial t} + \alpha_1 \nabla_x (\rho \bar{u}) + (\rho \bar{u}) \cdot \nabla_x (\alpha_1). $$

Recalling the definition of the material derivative and using Equations (3.8)-(3.26), we obtain

$$ (3.28) \qquad \frac{\partial}{\partial t}(\alpha_1 \rho) + \nabla_x \cdot (\alpha_1 \rho \bar{u}) = \rho \frac{D}{Dt}(\alpha_1) = -\frac{1}{\tau}(P_2 - P_1), $$

where $\tau := \frac{\sqrt{\varepsilon} \lambda (\eta_1 + \eta_2)^2}{\rho \eta_1 \eta_2}$ is the relaxation rate.

Thus, the pressure relaxation process drives a volume variation which evolves according to (3.28). In particular, we observe that the pressure relaxation implies a volume transfer towards the gas with highest pressure, according to the relaxation parameter $\tau > 0$.

**Remark 3.3.** The macroscopic velocity $\bar{u}$ for each point $(x, t)$ is given by (3.14), obtained from the kinetic equation as (2.28). Inside a region occupied only by one species p, the macroscopic velocity automatically reduces to the single-species macroscopic velocity $\bar{u}_p$.

Then, since we have conservation of the single-species mass and conservation of total momentum at the kinetic level, the derivation of (3.12c)-(3.12d) is obtained taking respectively $\varphi(v) = 1$, $m_p v$, averaging in space and summing over the index p, using Gauss and Leibniz rules and (3.28). We observe that, despite setting the
number of species M to 2, the validity of these equations remains independent from the number of species. The equation for total mass conservation (3.12a) is a straightforward consequence of the conservation of single-species mass. Finally, we would like to recover the equation for the relative velocity. Considering the test function $\varphi(v) = m_p v$ in the weak formulation (3.16), we obtain

$$ (3.29) \qquad \frac{\partial}{\partial t}(\rho_p \bar{u}_p) + \nabla_x \cdot (\rho_p \bar{u}_p \otimes \bar{u}_p + P_p) = \sum_{q \neq p} \mathbf{R}_{p,q}, $$

where $\mathbf{R}_{p,q}$ is the momentum exchange rate. Following the idea presented in [13], we shall compute explicitly $\mathbf{R}_{p,q}$ in the case of Maxwellian molecules and hard spheres molecules. To this aim, in Section A we compute the product of two local Maxwellians $M_p(v) M_q(w)$, whose expression depends on the center of mass $\mathcal{V}$ and relative velocities $|v - w|$. Then, assuming pseudo-Maxwellian molecules (constant kernel $\mathbf{B}$ independent from the species) or hard-spheres molecules (kernel of the form $\mathbf{B} = |u|$ independent from the species), the explicit computation described in A leads to the following expression for the momentum exchange rate:

$$ (3.30) \qquad \mathbf{R}_{p,q} = -\xi_{p,q} \rho_p \rho_q (\bar{u}_p - \bar{u}_q), \quad \xi_{p,q} \in \mathbb{R}, $$

where the coefficients $\xi_{p,q}$ satisfy the symmetric property $\xi_{p,q} = \xi_{q,p}$. We first notice that the computation can be done for general distribution function $f_p$ and $f_q$, if we are considering pseudo-Maxwellian molecules. Indeed, for this particular choice, the interaction kernel does not depend on the deviation angle nor the relative velocity and the exchange rate can be written using (2.13) and the expression for the moments as

$$ (3.31) \qquad \mathbf{R}_{p,q} = \int_{\mathbb{R}^d \times \mathbb{R}^d \times \mathbb{S}^{d-1}} m_p m_q f_p(v) f_q(v_*) (v - v_*) d\sigma dv dv_* = \rho_p \rho_q (\bar{u}_p - \bar{u}_q). $$

For general interaction particle, the presence of the interaction kernel make the computation really insidious and their explicit expressions are given in the Appendix A.

Now, we can rewrite explicitly the equation (3.29) and subtract the equation for the other species, obtaining

$$ (3.32) \qquad \frac{1}{\rho_1} \frac{\partial}{\partial t}(\rho_1 \bar{u}_1) - \frac{1}{\rho_2} \frac{\partial}{\partial t}(\rho_2 \bar{u}_2) \\ \qquad \qquad + \frac{1}{\rho_1} \nabla_x \cdot (\rho_1 \bar{u}_1 \otimes \bar{u}_1 + P_1) - \frac{1}{\rho_2} \nabla_x \cdot (\rho_2 \bar{u}_2 \otimes \bar{u}_2 + P_2) \\ \qquad \qquad = \frac{\mathbf{R}_{1,2}}{\rho_1} - \frac{\mathbf{R}_{2,1}}{\rho_2} = -\xi (\rho_2 + \rho_1) (\bar{u}_1 - \bar{u}_2). $$

We can rewrite the pressure term in terms of the internal energy and the phase enthalpy as

$$ (3.33) \qquad P_p = \rho_p (h_p - e_p), $$

where $h_p$ is the phase enthalpy and $e_p$ is the internal energy. Assuming a polytropic gas for which $P(\rho) = \rho^\gamma$, the following relation holds

$$ (3.34) \qquad \partial_x P_p = \rho_p \partial_x h_p. $$

Indeed, for an ideal gas it holds $P = \rho e(\gamma - 1)$

$$ (3.35) \qquad \frac{1}{\rho} \nabla_x P = \frac{\gamma}{\rho} \rho^{\gamma-1} \nabla_x \rho, \qquad \partial_x h = \nabla_x \left( \frac{P}{\rho} + e \right) = \gamma \rho^{\gamma-2} \nabla_x \rho, $$

which leads to (3.34).

Therefore, using (3.30) and (3.33) and defining the quantity

$$ (3.36) \qquad \zeta := \xi (\rho_1 + \rho_2) $$
we recover the following equation in the one-dimensional case

$$ (3.37) \qquad \frac{\partial}{\partial t}(\bar{u}_1 - \bar{u}_2) + \frac{\partial}{\partial x}\left(\frac{1}{2}\bar{u}_1^2 - \frac{1}{2}\bar{u}_2^2 + h_1 - h_2\right) = -\zeta(\bar{u}_1 - \bar{u}_2). $$

Therefore, through a formal hydrodynamic limit of the multi-species Boltzmann equation when $\epsilon \to 0$ we have obtained a model similar to the two-phase flow model introduced in [47].

### 3.2.1. Derivation of the continuous stress tensor at the interface.
The idea is to obtain an expression for the stress tensor at the boundary separating two phases, inspired by the the compressible boundary layer equations derived from the Boltzmann equation in [25] (see also [40] for the incompressible case).

Let us focus on the thin boundary layer across the interface of Fig. 3 and let us introduce a Prandtl scaling. From Prandtl theory, we know that if the Knudsen number is $O(\epsilon)$ and the fluid is compressible, then the thickness of such a region can be assumed $O(\sqrt{\epsilon})$. For the sake of simplicity, we can assume that the interface is located at $x_D = 0$ and the species 1 moves in the region described by the half space $x' = (x_1, \dots, x_{D-1})$ with $x_D < 0$. Let us rescale the independent variables by obtaining:

$$ (3.38) \qquad x' = x, \qquad x'_D = \frac{x_D}{\sqrt{\epsilon}}, \qquad \text{and} \qquad t' = t. $$

The rescaled Boltzmann equation reads as

$$ (3.39) \qquad \partial_t f_p + v' \cdot \nabla_{x'} f_p + \frac{1}{\sqrt{\epsilon}} v_D \partial_{x_D} f_p = \frac{1}{\epsilon} Q_{pp}(f_p, f_p) + Q_{pq}(f_p, f_q), $$

where the advection term in the D-direction describes the effect of the boundary layer at the interface. Since the process is “local” across the interface, no spatial averaging procedure is performed.

By using the Chapman-Enskog expansion introduced before in the weak formulation of (3.39) and matching different powers of $\epsilon$, the truncation at the 0-th order level needs $Q_{pp} = 0$ (which implies $f_p$ must be a local Maxwellian) and

$$ (3.40) \qquad \frac{1}{\sqrt{\epsilon}} \int_{\mathbb{R}^d} \varphi(v) v_D \partial_{x_D} \left( f_p^0 + \epsilon f_p^1 \right) dv = 0. $$

Since (3.40) is the only $1/\sqrt{\epsilon}$ term, this implies

$$ (3.41) \qquad \int_{\mathbb{R}^d} \varphi(v) v_D \partial_{x_D} \left( f_p^0 + \epsilon f_p^1 \right) dv = 0. $$

Let us first take $\varphi(v) = m_p$ in (3.40) and sum over the species $p$, from which using (3.5) we have

$$ (3.42) \qquad \partial_{x_D} \left( \sum_p \rho_p \bar{u}_p^D \right) = 0, $$

where $u^D$ represents the component of $u$ in the D-direction.

Then, taking $\varphi(v) = m_p v_D$ in (3.40) and summing over the species $p$, we obtain

$$ (3.43) \qquad \partial_{x_D} \left( \sum_p \rho_p (u_p^D)^2 + \sum_p P_p - \sum_p \sqrt{\epsilon} \lambda \partial_{x_D} \bar{u}_p^D - \sum_p \epsilon \lambda \mathcal{D}_{x'}(\bar{u}_p) \right) = 0, $$

where $\lambda$ is the proper viscosity coefficient. The components depending on $\sqrt{\epsilon}$ and $\epsilon$ are obtained as corrections in $\epsilon$ in the Chapman-Enskog expansion. Indeed, as shown in [39, 49], the conservation law at order 0 and the Hilbert's Lemma leads to a first-order flux correction of the form $\epsilon \operatorname{div}_x(\lambda \mathcal{D}(\bar{u}_p))$, where

$$ (3.44) \qquad \mathcal{D}(u) = \frac{1}{2} (\nabla_x u + (\nabla_x u)^T) - \frac{1}{3} (\nabla_x \cdot u) \mathrm{Id}, $$
namely $D(\bar{u}_p)$ is the traceless part of the deformation tensor of $\bar{u}_p$. In the one dimensional case it reduces to $D(u) = \partial_x u$.

Let us focus on the expression (3.43) at the interface $x_D = 0$. We recall that the two species are separated, therefore we can assume without loss of generality that the half-plane $x_D < 0$ is occupied by the species 1 and the half-plane $x_D > 0$ is occupied by the species 2. This implies that the partial momentum $\rho_2 u_2$ and the partial pressure $P_2$, related to the species 2, are null in the half-plane $x_D < 0$ and viceversa. The first equation (3.42) implies

$$ (3.45) \qquad \rho_1(t, x(t)) \bar{u}_1(t, x(t)) = k = \rho_2(t, \tilde{x}(t)) \bar{u}_2(t, \tilde{x}(t)), $$

where $x_D < 0$ and $\tilde{x}_D > 0$. In other words the interface is a contact discontinuity separating the two gases.

By using this relation in the second equation (3.43) and truncating at first order in $\epsilon$, we obtain

$$ (3.46) \qquad \begin{aligned} & k(\bar{u}_1(t, x(t)) + P_1(t, x(t))) - \sqrt{\epsilon} \lambda \frac{\partial \bar{u}_1(t, x(t))}{\partial x_D} \\ & = k \bar{u}_2(t, \tilde{x}(t)) + P_2(t, \tilde{x}(t)) - \sqrt{\epsilon} \lambda \frac{\partial \bar{u}_2(t, \tilde{x}(t))}{\partial x_D}, \end{aligned} $$

where $x_D < 0$ and $\tilde{x}_D > 0$. By imposing continuity conditions of the normal velocity at the interface, we finally obtain the continuity of the stress tensor at the boundary at first order in $\epsilon$, namely

$$ (3.47) \qquad P_1 - P_2 = \sqrt{\epsilon} \lambda \partial_{x_D} \bar{u}_1(t, x(t)) - \sqrt{\epsilon} \lambda \partial_{x_D} \bar{u}_2(t, x(t)). $$

### 3.2.2. Comparison with the Baer-Nunziato model.
Baer-Nunziato-type models [5, 1], differently to the model of two-phase compressible flow with two velocities and two pressures proposed in [48], are not in a conservative form. From a modeling point of view, the difference between these two approaches relies in the physical description of the phenomena: the conservative model has been proposed to deal with processes in which the thermal state of phases are almost in equilibrium state, while Baer-Nunziato-type models can take into account thermal non-equilibrium processes. To this aim, in this subsection we propose a comparison with the isentropic model derived above, by considering Baer-Nunziato models in which thermal processes are ignored and the phases behavior is isentropic. As observed in [48, 52], the one-dimensional versions of the models are very similar, if properly rewritten. The main difference is the definition of an interfacial pressure, which provides different results in the modelling of the phenomena. The isentropic version of a simplified Baer-Nunziato model reads as follows

$$ (3.48a) \qquad \frac{\partial}{\partial t}(\alpha_1) + u_{\text{interface}} \partial_x(\alpha_1) = -\frac{1}{\tau}(P_2 - P_1), $$

$$ (3.48b) \qquad \frac{\partial}{\partial t}(\alpha_1 \rho_1) + \partial_x(\alpha_1 \rho_1 \bar{u}_1) = 0, $$

$$ (3.48c) \qquad \frac{\partial}{\partial t}(\alpha_1 \rho_1 \bar{u}_1) + \partial_x(\alpha_1 \rho_1 \bar{u}_1^2 + \alpha_1 P_1) = P_{\text{interface}} \partial_x \alpha_1 + \zeta(\bar{u}_2 - \bar{u}_1), $$

$$ (3.48d) \qquad \frac{\partial}{\partial t}(\alpha_2 \rho_2) + \partial_x(\alpha_2 \rho_2 \bar{u}_2) = 0, $$

$$ (3.48e) \qquad \frac{\partial}{\partial t}(\alpha_2 \rho_2 \bar{u}_2) + \partial_x(\alpha_2 \rho_2 \bar{u}_2^2 + \alpha_2 P_2) = P_{\text{interface}} \partial_x \alpha_2 - \zeta(\bar{u}_2 - \bar{u}_1), $$
where $u_{\text{interface}}$ and $P_{\text{interface}}$ are the interfacial velocity and the interfacial pressure respectively. The definition of interfacial velocity is the same for both models

$$ (3.49) \qquad u_{\text{interface}} = \bar{u} = \frac{\alpha_1 \rho_1 \bar{u}_1 + \alpha_2 \rho_2 \bar{u}_2}{\alpha_1 \rho_1 + \alpha_2 \rho_2}, $$

since the first equation (3.48a) could be easily derived from (3.12b)-(3.12c). Equations (3.48b)-(3.48d) are derived from the multi-species Boltzmann equation (3.15), as already shown for the conservative model.

The main difference lies in the equation for the balance of the single-species momentum. Indeed, by taking $\varphi(v) = m_p v$ in (3.16), followed by a spatial averaging procedure, we recover Equations (3.48c)-(3.48e), where $\tilde{\zeta}$ now is given by

$$ (3.50) \qquad \tilde{\zeta} := \frac{\xi}{\alpha_1 \alpha_2} \rho. $$

We conclude that both models could be derived from the multi-species Boltzmann equation (3.15), by considering different closure assumptions and different quantities of interest.

## 4. CONCLUSIONS

In this paper we have formally derived the hyperbolic model introduced in [48] to reproduce the dynamics of two-phase flows, starting from a kinetic model for gas mixtures. We have first investigated the properties of the multi-species Boltzmann equations and we have given a brief description of the behavior of the limit solution. We have then examined the hydrodynamic limit of the kinetic model under the regime of resonant intra-species collisions, aiming to derive different equilibrium velocities for the two species. This choice allows to recover a description for two-phase flows at the macroscopic level. We have kept an isentropic closure assumption on the pressure term, in order to recover the original model. This assumption could be weakened in future works, leading to different macroscopic equations for multi-phase flows, see for instance [53].

Several new directions could be pursued. The idea is to develop an asymptotic-preserving numerical scheme for the stiff kinetic model (3.15). This will enable numerical simulations, to provide a comparison between the behavior of the kinetic solution and that of the macroscopic equations (3.12a)-(3.12e) in the limit $\varepsilon \to 0$. Further, we could extend this strategy to more complex kinetic models, taking into account Boltzmann equations with different type of collisions, such as mixtures of gases undergoing inelastic collisions [46].

## ACKNOWLEDGEMENT

GP and TT would like to thank Michael Dumbser for the fruitful discussion on this topic. TT would also like to thank Andrea Bondesan for his help.

The authors received funding from the European Union's Horizon Europe research and innovation program under the Marie Skłodowska-Curie Doctoral Network DataHyking (Grant No. 101072546). TR and TT are also supported by the French government, through the UniCAJEDI Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR-15-IDEX-01. 46). GP is also supported by MUR (Ministry of University and Research) under the PRIN-2022 project on "High order structure-preserving semi-implicit schemes for hyperbolic equations" (number 2022JH87B4).

GP and TT are members of the INdAM Research National Group of Scientific Computing (INDAM-GNCS).
**Data Availability.** Data sharing not applicable to this article as no datasets were generated or analysed during the current study.

**Conflict of interest.** No conflict of interest is extant in the present work.

## APPENDIX A. COMPUTATION OF THE MOMENTUM EXCHANGE RATE FOR HARD-SPHERES MOLECULES

We compute here the momentum exchange rate $\mathbf{R}_{p,q}$, describing the change in momentum for the species $p$ due to the interactions with the other species $q$. We partially follow [24] and the more recent paper [42]. We assume a collision kernel $\mathbf{B}$ independent from the species. Since we are interested in an approximation of the exchange rate in the limit for $\varepsilon \to 0$, it is reasonable to consider $f_p$ and $f_q$ as two Maxwellian distribution functions. In particular we can write the exchange rate in the following form

$$ (A.1) \quad \mathbf{R}_{p,q} = \int_{\mathbb{R}^{d} \times \mathbb{R}^{d} \times \mathbb{S}^{d-1}} m_p m_q \mathcal{M}_p(v) \mathcal{M}_q(v_*) (v - v_*) \mathbf{B}(|v - v_*|, \cos \theta) d\sigma dv dv_* $$

Let us consider the product of two local Maxwellian

$$ (A.2) \quad \mathcal{M}_p(v) \mathcal{M}_q(v_*) = n_p \left(\frac{m_p}{2\pi T}\right)^{\frac{d}{2}} \exp\left(-\frac{m_p}{2T}|\bar{u}_p - v|^2\right) n_q \left(\frac{m_q}{2\pi T}\right)^{\frac{d}{2}} \exp\left(-\frac{m_q}{2T}|\bar{u}_q - v_*|^2\right). $$

We introduce the following relative center of mass with respect to the velocities at the equilibrium. We can define

$$ (A.3) \quad \mathcal{W} = \frac{m_p \bar{u}_p + m_q \bar{u}_q}{m_p + m_q} - \mathcal{V}, $$

where $\mathcal{V}$ is the center of mass defined in (2.5). Now, we can rewrite

$$ \bar{u}_p - v = \mathcal{W} + \frac{m_q}{m_p + m_q} [(\bar{u}_p - v) - (\bar{u}_q - v_*)], $$

$$ \bar{u}_q - v_* = \mathcal{W} + \frac{m_p}{m_p + m_q} [(\bar{u}_q - v_*) - (\bar{u}_p - v)]. $$

Using this new variable, we can rewrite the product (A.2) as

$$ (A.4) \quad \mathcal{M}_p(v) \mathcal{M}_q(v_*) = K_{pq} \exp \left( - \frac{|\mathcal{W}|^2}{\mu_{pq}} - \frac{|(v - v_*) - (\bar{u}_p - \bar{u}_q)|^2}{\gamma_{pq}} \right), $$

where

$$ (A.5) \quad K_{pq} = n_p n_q \frac{(m_p m_q)^{\frac{d}{2}}}{(2\pi T)^d}, \quad \mu_{pq} = \frac{2T}{(m_p + m_q)}, \quad \gamma_{pq} = \frac{2(m_p + m_q)T}{m_p m_q}. $$

Now we can perform the integral with respect the new variables $\mathcal{W}$ and $w := v - v_*$, observing that the Jacobian of the transformation is unitary, which implies $dv dv_* = dw d\mathcal{W}$. Hence, integrating firstly over $d\mathcal{W}$:

$$ (A.6) \quad \mathbf{R}_{p,q} = \pi^{\frac{d}{2}} \mu_{pq}^d K_{pq} m_p m_q \int_{\mathbb{R}^{d} \times \mathbb{S}^{d-1}} \mathbf{B}(|w|, \cos \theta) w \exp \left( - \frac{|w - (\bar{u}_p - \bar{u}_q)|^2}{\gamma_{pq}} \right) d\sigma dw $$

We now restrict our analysis to the hard-spheres case for $d = 3$, where the interaction kernel is

$$ (A.7) \quad \mathbf{B}(|w|, \cos \theta) = |w|, $$

from which we obtain

$$ (A.8) \quad \mathbf{R}_{p,q} = 2\pi^{\frac{d}{2}} \mu_{pq}^3 K_{pq} m_p m_q \int_{\mathbb{R}^3} |w| w \exp \left( - \frac{|w - (\bar{u}_p - \bar{u}_q)|^2}{\gamma_{pq}} \right) dw. $$
Let us rewrite $w$ in spherical coordinates

$$ (A.9) \qquad w = |w|(\sin \varphi \cos \psi, \sin \varphi \sin \psi, \cos \varphi). $$

Then,

$$ \begin{aligned} I_{p,q} &:= \int_{\mathbb{R}^3} |w| w \exp \left( - \frac{|w - (\bar{u}_p - \bar{u}_q)|^2}{\gamma_{pq}} \right) dw \\ &= \int_0^\pi \int_0^{2\pi} \int_0^\infty |w|^3 \hat{w} \exp \left( - \frac{|w - (\bar{u}_p - \bar{u}_q)|^2}{\gamma_{pq}} \right) \sin \varphi \, d\varphi \, d\psi \, d|w| \\ &= 2\pi \hat{e}_3 \int_0^\pi \int_0^\infty |w|^4 \exp \left( - \frac{|w - (\bar{u}_p - \bar{u}_q)|^2}{\gamma_{pq}} \right) \cos \varphi \sin \varphi \, d\varphi \, d|w|, \end{aligned} $$

where $\hat{w} = w/|w|$ and $\hat{e}_3$ could be chosen in the direction of $\bar{u}_p - \bar{u}_q$.

Introducing the first change of variable

$$ (A.10) \qquad \alpha = (\bar{u}_p - \bar{u}_q) \cos \varphi, \quad \beta = |w|, \quad \nu = \bar{u}_p - \bar{u}_q, $$

the integral becomes

$$ (A.11) \qquad I_{p,q} = \frac{2\pi}{|\nu|^3} \nu \exp\left(-\frac{\nu^2}{\gamma_{pq}}\right) \int_{-\nu}^{\nu} \int_{0}^{\infty} \beta^4 \alpha \exp\left(-\frac{\beta^2}{\gamma_{pq}} - \frac{2\beta\alpha}{\gamma_{pq}}\right) d\alpha d\beta. $$

Now we introduce another change of variable

$$ (A.12) \qquad \eta = \beta + \alpha, $$

in order to rewrite $I_{p,q}$ as

$$ I_{p,q} = \frac{2\pi}{|\nu|^3} \nu \exp\left(-\frac{\nu^2}{\gamma_{pq}}\right) \int_{-\nu}^{\nu} \int_{0}^{\infty} \alpha(\eta - \alpha)^4 \exp\left(-\frac{\eta^2}{\gamma_{pq}} + \frac{\alpha^2}{\gamma_{pq}}\right) d\alpha d\eta. $$

Let us consider the rescaled variables

$$ \eta = \frac{\eta}{\sqrt{\gamma_{pq}}}, \quad \alpha = \frac{\alpha}{\sqrt{\gamma_{pq}}}. $$

Now, avoiding the label in $\gamma_{pq}$, we have

$$ I_{p,q} = \frac{2\pi}{|\nu|^3} \gamma^{11/2} \nu \exp\left(-\frac{\nu^2}{\gamma}\right) \int_{-\frac{\nu}{\gamma}}^{\frac{\nu}{\gamma}} \int_{0}^{\infty} \alpha(\eta - \alpha)^4 \exp(-\eta^2 + \alpha^2) d\alpha d\eta. $$

We first perform the integration in $d\eta$,

$$ \begin{aligned} \int_{\alpha}^{\infty} (\eta - \alpha)^4 \exp(-\eta^2) d\eta &= \frac{1}{8} \sqrt{\pi} (3 + 12\alpha^2 + 4\alpha^4) \operatorname{erfc}(\alpha) - \frac{1}{4} \alpha \exp(-\alpha^2) (5 + 2\alpha^2), \end{aligned} $$

where $\operatorname{erfc}$ is the complementary error function.

The final expression for $I_{p,q}$ is

$$ \begin{aligned} I_{p,q} &= \frac{\pi}{4|\nu|^3} \gamma^{11/2} \nu e^{\left(-\frac{\nu^2}{\gamma}\right)} \left(-4\left(\frac{\nu}{\gamma}\right)^3 + 2\left(\frac{\nu}{\gamma}\right) - \sqrt{\pi}\left(4\left(\frac{\nu}{\gamma}\right)^4 + 4\left(\frac{\nu}{\gamma}\right)^2 - 1\right) e^{\left(\frac{\nu}{\gamma}\right)^2} \operatorname{erf}\left(\frac{\nu}{\gamma}\right)\right) \\ &= \pi \nu \left[e^{\left(-\frac{\nu^2}{\gamma}\right)} \left(\gamma^{5/2} + \frac{\gamma^{9/2}}{2\nu^2}\right) - \sqrt{\pi} \left(\nu\gamma^{3/2} + \frac{\gamma^{7/2}}{\nu} - \frac{\gamma^{11/2}}{4\nu^3}\right) \operatorname{erf}\left(\frac{\nu}{\gamma}\right)\right]. \end{aligned} $$

We can conclude that

$$ \mathbf{R}_{p,q} = \rho_p \rho_q (\bar{u}_p - \bar{u}_q) \Psi_{p,q}, $$
where $\Psi_{p,q} = \frac{(m_p m_q)^{\frac{d}{2}}}{\sqrt{\pi}(m_p + m_q)} \frac{I_{p,q}}{(\bar{u}_p - \bar{u}_q)}$.

REFERENCES

[1] ABGRALL, R., AND SAUREL, R. Discrete equations for physical and numerical compressible multiphase mixtures. *J. Comput. Phys.* **186**, 2 (2003), 361–396.

[2] ALONSO, R. J., ČOLIĆ, M., AND GAMBA, I. M. The Cauchy problem for Boltzmann bi-linear systems: the mixing of monatomic and polyatomic gases. *J. Stat. Phys.* **191**, 1 (2024), 50.

[3] ANDRIES, P., AOKI, K., AND PERTHAME, B. A consistent BGK-type model for gas mixtures. *J. Stat. Phys.* **106**, 5–6 (2002), 993–1018.

[4] ARNAULT, P., AND GUISSET, S. Chapman–Enskog derivation of multicomponent Navier–Stokes equations. *Physics of Plasmas* **29**, 9 (2022).

[5] BAER, M. R., AND NUNZIATO, J. W. A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials. *Int. J. Multiphase Flow* **12** (1986), 861–889.

[6] BIANCA, C., AND DOGBE, C. On the Boltzmann gas mixture equation: linking the kinetic and fluid regimes. *Commun. Nonlinear Sci. Numer. Simul.* **29**, 1–3 (2015), 240–256.

[7] BIRD, G. A. *Molecular Gas Dynamics And The Direct Simulation Of Gas Flows*. Oxford University Press, 05 1994.

[8] BISI, M., AND DESVILLETTE'S, L. Formal passage from kinetic theory to incompressible Navier-Stokes equations for a mixture of gases. *ESAIM, Math. Model. Numer. Anal.* **48**, 4 (2014), 1171–1197.

[9] BISI, M., GROPP, M., LUCCHI, E., AND MARTALÒ, G. A mixed Boltzmann-BGK model for inert gas mixtures. *Kinet. Relat. Models* **17**, 5 (2024), 674–696.

[10] BISI, M., GROPP, M., AND MARTALÒ, G. Macroscopic equations for inert gas mixtures in different hydrodynamic regimes. *J. Phys. A, Math. Theor.* **54**, 8 (2021), 23.

[11] BISI, M., GROPP, M., AND MARTALÒ, G. A space-dependent Boltzmann-BGK model for gas mixtures and its hydrodynamic limits. Preprint, arXiv:2505.19735 [math-ph] (2025), 2025.

[12] BISI, M., MARTALÒ, G., AND SPIGA, G. Multi-temperature Euler hydrodynamics for a reacting gas from a kinetic approach to rarefied mixtures with resonant collisions. *Europhysics Letters* **95**, 5 (2011), 55002.

[13] BISI, M., MARTALÒ, G., AND SPIGA, G. Multi-temperature hydrodynamic limit from kinetic theory in a mixture of rarefied gases. *Acta Appl. Math.* **122**, 1 (2012), 37–51.

[14] BOBYLEV, A. V., BISI, M., GROPP, M., SPIGA, G., AND POTAPENKO, I. F. A general consistent BGK model for gas mixtures. *Kinet. Relat. Models* **11**, 6 (2018), 1377–1393.

[15] BONDESAN, A., BOUDIN, L., BRIANT, M., AND GREC, B. Stability of the spectral gap for the Boltzmann multi-species operator linearized around non-equilibrium Maxwell distributions. *Commun. Pure Appl. Anal.* **19**, 5 (2020), 2549–2573.

[16] BONDESAN, A., AND TANG, B. Q. Explicit spectral gap estimates for the linearized Boltzmann operator modeling reactive gaseous mixtures. Preprint, arXiv:2410.15055 [math.AP] (2024), 2024.

[17] BOSCARINO, S., CHO, S. Y., GROPP, M., AND RUSSO, G. BGK models for inert mixtures: comparison and applications. *Kinet. Relat. Models* **14**, 5 (2021), 895–928.

[18] BOUDIN, L., GREC, B., AND PAVAN, V. The Maxwell-Stefan diffusion limit for a kinetic model of mixtures with general cross sections. *Nonlinear Anal., Theory Methods Appl., Ser. A, Theory Methods* **159** (2017), 40–61.

[19] BOUDIN, L., GREC, B., AND SALVARANI, F. The Maxwell-Stefan diffusion limit for a kinetic model of mixtures. *Acta Appl. Math.* **136**, 1 (2015), 79–90.

[20] BRESCH, D., BURTEA, C., AND LAGOUTIÈRE, F. Mathematical justification of a compressible bi-fluid system with different pressure laws: a semi-discrete approach and numerical illustrations. *J. Comput. Phys.* **490** (2023), 39.

[21] BRIANT, M. Stability of global equilibrium for the multi-species Boltzmann equation in $L^\infty$ settings. *Discrete and Continuous Dynamical Systems* **36**, 12 (2016), 6669–6688.

[22] BRIANT, M., AND DAUS, E. S. The Boltzmann equation for a multi-species mixture close to global equilibrium. *Archive for Rational Mechanics and Analysis* **222** (2016), 1367–1443.

[23] BRIANT, M., AND GREC, B. Rigorous derivation of the Fick cross-diffusion system from the multi-species Boltzmann equation in the diffusive scaling. *Asymptotic Anal.* **135**, 1-2 (2023), 55–80.

[24] BURGERS, J. M. Flow equations for composite gases, vol. 11 of *Appl. Math. Sci.* New York etc.: Springer-Verlag, 1969.
[25] ČASTULÍK, P. The compressible boundary layer and the Boltzmann equation. *Appl. Math. Lett.* **15**, 1 (2002), 17–23.

[26] CERCIGNANI, C. *The Boltzmann equation and its applications*, vol. 67 of *Appl. Math. Sci.* New York etc.: Springer-Verlag, 1988.

[27] CERCIGNANI, C., ILLNER, R., AND PULVIRENTI, M. *The mathematical theory of dilute gases*, vol. 106 of *Appl. Math. Sci.* New York, NY: Springer-Verlag, 1994.

[28] CHAPMAN, S., AND COWLING, T. G. The mathematical theory of non-uniform gases. Cambridge: Cambridge University Press. xxiii, 404 p. (1939)., 1939.

[29] CHAPMAN, S., AND COWLING, T. G. The mathematical theory of non-uniform gases. An account of the kinetic theory of viscosity, thermal conduction, and diffusion in gases. Cambridge: Cambridge University Press; New York: The Macmillan Company. xxiv, 404 p. (1939)., 1939.

[30] ČOLIĆ, M. Explicit model for a mixture of Eulerian fluids based on kinetic theory. *Physics of Fluids* **36**, 6 (06 2024), 067134.

[31] DE LA CANAL, E., GAMBA, I. M., AND PAVIĆ-ČOLIĆ, M. Propagation of $L_{\beta}^p$-norm, $1 < p \le \infty$, for the system of Boltzmann equations for monatomic gas mixtures. Preprint, arXiv:2001.09204 [math.AP] (2020), 2020.

[32] DESVILLETTES, L., MONACO, R., AND SALVARANI, F. A kinetic model allowing to obtain the energy law of polytropic gases in the presence of chemical reactions. *European Journal of Mechanics - B/Fluids* **24**, 2 (2005), 219–236.

[33] DREW, D. A., AND PASSMAN, S. L. *Theory of multicomponent fluids*, vol. 135 of *Appl. Math. Sci.* New York, NY: Springer, 1998.

[34] DUMBSER, M. A simple two-phase method for the simulation of complex free surface flows. *Computer Methods in Applied Mechanics and Engineering* **200**, 9 (2011), 1204–1219.

[35] FLÄTTEN, T., AND LUND, H. Relaxation two-phase flow models and the subcharacteristic condition. *Math. Models Methods Appl. Sci.* **21**, 12 (2011), 2379–2407.

[36] GALKIN, V. S., AND MAKASHEV, N. K. Kinetic derivation of the gas-dynamic equations for multicomponent mixtures of light and heavy particles. *Fluid Dyn.* **29**, 1 (1994), 140–155.

[37] GAMBA, I. M., AND PAVIĆ-ČOLIĆ, M. On existence and uniqueness to homogeneous Boltzmann flows of monatomic gas mixtures. *Arch. Ration. Mech. Anal.* **235**, 1 (2020), 723–781.

[38] GODUNOV, S. K., AND ROMENSKY, E. I. Thermodynamics, conservation laws, and symmetric forms of differential equations in mechanics of continuous media. In *Computational fluid dynamics review 1995*. Chichester: John Wiley & Sons, 1995, pp. 19–30.

[39] GOLSE, F. The Boltzmann equation and its hydrodynamic limits. In *Handbook of differential equations: Evolutionary equations. Vol. II*. Amsterdam: Elsevier/North-Holland, 2005, pp. 159–301.

[40] GOLSE, F. From the Boltzmann equation to hydrodynamic equations in thin layers. *Boll. Unione Mat. Ital.* **(9)** 4, 2 (2011), 163–186.

[41] HAACK, J. R., HAUCK, C. D., AND MURILLO, M. S. A conservative, entropic multispecies BGK model. *J. Stat. Phys.* **168**, 4 (2017), 826–856.

[42] HUNANA, P., PASSOT, T., KHOMENKO, E., MARTÍNEZ-GÓMEZ, D., COLLADOS, M., TENERANI, A., ZANK, G. P., MANEVA, Y., GOLDSTEIN, M. L., AND WEBB, G. M. Generalized Fluid Models of the Braginskii Type. *The Astrophysical Journal Supplement Series* **260**, 2 (2022), 26.

[43] KLINGENBERG, C., PIRNER, M., AND PUPPO, G. A consistent kinetic model for a two-component mixture of polyatomic molecules. *Commun. Math. Sci.* **17**, 1 (2019), 149–173.

[44] LUKÁČOVÁ-MEDVIDOVÁ, M., PUPPO, G., AND THOMANN, A. An all Mach number finite volume method for isentropic two-phase flow. *J. Numer. Math.* **31**, 3 (2023), 175–204.

[45] LUND, H. A hierarchy of relaxation models for two-phase flow. *SIAM J. Appl. Math.* **72**, 6 (2012), 1713–1741.

[46] REY, T., AND TENNA, T. The Boltzmann equation for a multi-species inelastic mixture. Preprint, arXiv:2410.17027 [math.AP] (2024), 2024.

[47] ROMENSKI, E., DRIKAKIS, D., AND TORO, E. Conservative models and numerical methods for compressible two-phase flow. *J. Sci. Comput.* **42**, 1 (2010), 68–95.

[48] ROMENSKI, E., AND TORO, E. F. Compressible two-phase flows: two-pressure models and numerical methods. *Comput. Fluid Dyn. J* **13** (2004), 403–416.

[49] SAINT-RAYMOND, L. *Hydrodynamic limits of the Boltzmann equation*, vol. 1971 of *Lect. Notes Math*. Berlin: Springer, 2009.

[50] SAUREL, R., AND ABGRALL, R. A multiphase Godunov method for compressible multifluid and multiphase flows. *J. Comput. Phys.* **150**, 2 (1999), 425–467.

[51] SAUREL, R., AND LEMETAYER, O. A multiphase model for compressible flows with interfaces, shocks, detonation waves and cavitation. *J. Fluid Mech.* **431** (2001), 239–271.
[52] THEIN, F., ROMENSKI, E., AND DUMBSER, M. Exact and numerical solutions of the Riemann problem for a conservative model of compressible two-phase flows. *J. Sci. Comput.* **93**, 3 (2022), 60.

[53] THOMANN, A., AND DUMBSER, M. Thermodynamically compatible discretization of a compressible two-fluid model with two entropy inequalities. *J. Sci. Comput.* **97**, 1 (2023), 33.