# Diffusion asymptotics of a kinetic model for gas-particle mixtures with energy exchanges

Frédérique Charles, Annamaria Massimini, Francesco Salvarani

► To cite this version:

Frédérique Charles, Annamaria Massimini, Francesco Salvarani. Diffusion asymptotics of a kinetic model for gas-particle mixtures with energy exchanges. 2025. hal-05004987v2

HAL Id: hal-05004987

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

Preprint submitted on 16 Sep 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Diffusion asymptotics of a kinetic model for gas-particle mixtures with energy exchanges

Frédérique Charles¹, Annamaria Massimini², Francesco Salvarani³

¹: Université Grenoble Alpes, LJK, UMR CNRS 5224, 38058 GRENOBLE CEDEX 9, France.
frederique.charles@univ-grenoble-alpes.fr

²: CERMICS, Ecole des Ponts Paris-Tech and Inria MATHERIALS team-project, 6 et 8 av. Blaise Pascal Cité Descartes 77455 Marne-La-Vallée, France. annamaria.massimini@enpc.fr

³: De Vinci Higher Education, De Vinci Research Center, Paris, France & Dipartimento di Matematica "F. Casorati", Università di Pavia, Via Ferrata 5, 27100 Pavia, Italy. francesco.salvarani@unipv.it

## Abstract

In this paper, we propose a non-isothermal kinetic model describing interactions between dust particles and gas molecules. We assume that gas-dust collisions follow a diffuse reflection mechanism at the surface of the dust particles. The surface temperature of the particles is treated as a function of time and space, satisfying a transport-like equation. The key novelty of this model is its ability to derive both the conservation of total energy in the system and an explicit expression for entropy. After adimensionalising the equations, we perform a formal diffusive asymptotic analysis leading to a coupled system of partial differential equations. In particular, we examine the limiting case in which the mass ratio between gas and dust, as well as the Mach and Knudsen numbers, tend to zero.

Keywords: Kinetic theory; Boltzmann equation; Diffusive asymptotics

## 1 Introduction

The study of gas-dust mixtures at the mathematical level is of fundamental importance in several fields, such as astrophysics, environmental science, and industrial processes. At its core, this importance lies in the intricate and often nonlinear interactions between gas and dust particles, which can lead to complex dynamics and phenomena that are critical for understanding the behavior and evolution of various systems.

In astrophysics, for example, gas-dust mixtures play a key role in the formation of stars, planets, and other celestial bodies. In environmental science, understanding the dynamics of gas-dust mixtures is essential for modeling atmospheric phenomena such as pollutant dispersion, dust storms, and cloud formation. These mixtures affect air quality, weather patterns, and climate change. In various industrial processes, such as pharmaceutical manufacturing, mining, and chemical processing, controlling the behavior of gas-dust mixtures is critical for efficiency, safety, and environmental compliance. In reactive gas-dust systems, the presence of solid particles can catalyze chemical reactions or alter the heat transfer properties of the gas. This is particularly relevant in combustion processes, where dust may either inhibit or enhance combustion, depending on the conditions. Mathematical modeling of these processes is essential for the design and optimization of reactors, combustion engines, and pollution control systems.

When dealing with diffusive phenomena, the fundamental model in the framework of multicomponent mixtures is given by the Maxwell-Stefan equations. Introduced in the scientific literature by James Clerk Maxwell and Josef Stefan in the 19th century [27, 29], these equations have become very popular in the study and modeling of diffusion processes in gaseous mixtures. The equations offer a more detailed and accurate description of mass transport in mixtures than classical Fickian diffusion models [21, 22], especially when dealing with high concentration gradients and non-ideal mixtures.

Although widely used in chemical engineering [30], numerical studies of the Maxwell-Stefan system and its variants begun around the end of the 20th century [20, 23]. Starting from [5] and [7], the literature has shown a renewed interest in these equations in the mathematical community.
In particular, a key research direction is the derivation of multicomponent diffusion equations from a system of Boltzmann equations, establishing a bridge between microscopic particle dynamics and macroscopic gas mixture behavior. This derivation relies on the assumption of local equilibrium, meaning that over sufficiently small spatial scales, the system can be approximated by a local Maxwell-Boltzmann distribution. This assumption is valid when the mean free path of gas molecules is much smaller than the characteristic length scale over which macroscopic gradients develop.

This scaling approach has been extensively studied, particularly for non-reactive mixtures [8, 6, 25, 4], and extended to reactive cases in [1, 2, 3] under different frameworks.

In the case of models for gas-particle mixtures, the starting point is again a system of two Boltzmann equations (one for the evolution of gas molecules, and one for the evolution of dust particles). The interaction between the particles and the surrounding gas is described by specific collision operators, which take into account the macroscopic character of the particles compared to the gas molecules. In [12], the collisional mechanism is based on the hypothesis of diffuse reflection (at a fixed and constant temperature) of the molecules on the surface of the particles. The collisions preserve the total momentum, but not the kinetic energy. However, the model has an entropy [14] (different from the classical Boltzmann entropy). In order to take into account the heat exchange between the gas and the dispersed phase, the authors have modified the model of [12] by introducing the surface temperature of the particles as a variable. The model preserves the total energy of the system, but no explicit expression for the entropy has been provided.

Asymptotic models can be derived by letting the mass ratio between a gas molecule and a particle go to 0. They are very interesting for numerical simulation purposes, since the discretization of the kinetic model is very expensive when the mass ratio is weak [13]. In the asymptotic model proposed in [15], the evolution of the particle density is described by the Vlasov equation with a term that can be interpreted as a drag force exerted by the gas on the particles. This term has some similarities with existing models of the drag force, such as the formula of Otterman and Levinet [28] (for a two-phase fluid model), or that of [17], which proposes a correction of the Stokes drag force to extend its validity as the Knudsen number increases. The derivation of macroscopic models for binary mixtures with large mass disparities is also studied in [18], devoted to the study of Grad's epochal relaxation phenomenon. When the gas and the dispersed phase are coupled not only by the drag force but also by the temperature exchange between the gas and the droplets/dust particles, the asymptotic model also includes a term describing the internal energy exchange between the species [16].

In the case of a porous media [31], homogenized kinetic equations to derive effective transport properties (such as permeability, Knudsen diffusivity, and thermal transpiration ratio) have been investigated. The authors of [31] describe the evolution of the gas via a Boltzmann equation, while the dust is fixed in space and assumed to be part of the spatial domain. The interaction of the gas with the solid surface is described by diffuse reflection boundary conditions, where the surface temperature of the dust is considered as a function of time and space. In particular, it satisfies a classical heat equation.

Our work fits within the dynamics of the previously described literature on non-isothermal dust-gas mixtures. As in previous studies, such as [12], we have investigated gas-dust interactions using a diffuse reflection mechanism. However, the novelty of our model lies in the fact that it possesses an explicit entropy structure and, at the same time, it conserves the total energy of the system.

Rather than following the methodology in [16], we draw inspiration from the porous media literature to treat the surface temperature as a macroscopic quantity. This approach leads to a transport-like equation that ensures the conservation of total energy within the system.

Our study has two main objectives. First, we aim to provide a comprehensive description of non-isothermal dust-gas mixtures within this new framework. Second, we seek to perform a formal asymptotic diffusive limit, leading to a coupled system. This analysis offers deeper insight into the coupling mechanisms between gas and dust particles and their thermodynamic behavior.

The structure of the article is the following. After the Introduction, in Section 2 and 3, we introduce the model and describe its properties. In Section 4 we perform the adimensionalization of the system and introduce small parameters. Macroscopic equations for the macroscopic scaled
quantities are derived in Section 5. The equation system obtained by formally making the small parameters tend to zero is presented in Section 6. In particular, we prove that, under a hard-sphere collision kernel, the gas kinetic temperature, the particle density temperature and the particle surface temperature all converge to the same limit function.

## 2 The model

Consider a binary aerosol mixture consisting of solid dust particles or liquid droplets suspended in a monatomic carrier gas. The dust particles are assumed to be identical, spherical (i.e., all having the same radius and mass), and significantly larger and more massive than the gas molecules. Both species are described by appropriate density functions. In what follows, the function $F = F(t, x, v)$ represents the density of particles, whereas $f = f(t, x, w)$ denotes the density of gaseous atoms (for example, a noble gas, such as Helium or Neon).

### 2.1 The Boltzmann system

At the kinetic level, we describe the gas-particle system by the following system of Boltzmann-like equations

$$ \left\{ \begin{aligned} \partial_t F + v \cdot \nabla_x F &= \mathcal{D}(F, f) + \mathcal{P}(F), \\ \partial_t f + w \cdot \nabla_x f &= \mathcal{R}(F, f) + \mathcal{C}(f), \end{aligned} \right. \qquad (1) $$

where the exact form of the four collision operators $\mathcal{P}$, $\mathcal{D}$, $\mathcal{C}$ and $\mathcal{R}$ depends on the precise microscopic behavior of the collisions. It is worth noting that the form of the collision kernel can have a significant influence on the mathematical properties of the system, even though it does not alter the bilinear structure of the collision operator.

### 2.2 Mono-species collision operators

In (1), $\mathcal{C}(f)$ is the classical Boltzmann operator for monoatomic gases, describing interactions between molecules of the same species [11]:

$$ \mathcal{C}(f)(w) = \sigma_{gg} \int_{\mathbb{R}^3} \int_{S^2} [f(w'')f(w_*'') - f(w)f(w_*)] c(w - w_*, n) \, dndw_*, \qquad (2) $$

where $\sigma_{gg}$ is a physical constant of the order of magnitude of $r_g^2$, with $r_g$ the radius of a molecule, and the pre-post collisional velocity rules are given by:

$$ (w'', w_*'') = (w - [(w - w_*) \cdot n]n, w_* + [(w - w_*) \cdot n]n). \qquad (3) $$

In particular, we see from (3) that the momentum and the kinetic energy are preserved during each collision:

$$ \left\{ \begin{aligned} w'' + w_*'' &= w + w_*, \\ |w''|^2 + |w_*''|^2 &= |w|^2 + |w_*|^2. \end{aligned} \right. \qquad (4) $$

The nonnegative function $c$ (collision kernel) in (2) describes the interactions between molecules and can be written as

$$ c(w - w_*, n) = |w - w_*| q \left( |w - w_*|, \frac{w - w_*}{|w - w_*|} \cdot n \right), $$

where $q$ is a dimensionless function. This operator has the following weak formulations: for any test function $\psi \in L^\infty(\mathbb{R}^3)$, we have (see [11]):

$$ \int_{\mathbb{R}^3} \mathcal{C}(f)(w)\psi(w)dw = \sigma_{gg} \int_{\mathbb{R}^3} (\psi(w') - \psi(w))f(w)f(w_*)c(w - w_*, n) \, dw \, dndw_*, \qquad (5) $$
and

$$ \int_{\mathbb{R}^3} \mathcal{C}(f)(w)\psi(w)dw = -\frac{\sigma_{gg}}{2} \int_{\mathbb{R}^3} [\psi(w')-\psi(w)][f(w'')f(w''')-f(w)f(w')]c(w-w_*,n) dwdndw_* \quad (6) $$

We assume that $\mathcal{P}(F)$ is also a Boltzmann operator, describing the elastic collisions between dust particles. The expression of $\mathcal{P}(F)(v)$ is similar to (2), where $\sigma_{gg}$ is replaced by $\sigma_{pp}$ (which has order of magnitude $r_p^2$, where $r_p$ is the radius of the particles), and $q$ is replaced by another dimensionless function, $p$.

The weak formulations of $\mathcal{P}(F)$ are very similar to those in (5) and (6). We recall the properties of $\mathcal{C}(f)$ and $\mathcal{P}(F)(v)$ concerning the conservation of mass, momentum, and kinetic energy, which follow directly from (5) and (4):

$$ \int_{\mathbb{R}^3} \mathcal{C}(f)(w) \begin{pmatrix} 1 \\ w \\ |w|^2/2 \end{pmatrix} dw = \int_{\mathbb{R}^3} \mathcal{P}(f)(v) \begin{pmatrix} 1 \\ v \\ |v|^2/2 \end{pmatrix} dv = \begin{pmatrix} 0 \\ 0 \\ 0 \end{pmatrix}. \quad (7) $$

## 2.3 Bi-species collision operators

In gas mixtures, the collision mechanism between two molecules of different masses, when neither internal energy exchange nor chemical reactions are involved, can be described by elastic collisions (see, for example, [10] and [13]). In [12], a new model was proposed in which collisions between dust particles and gas molecules are assumed to be inelastic and governed by a diffuse reflection mechanism (at a fixed temperature) on the surface of the dust particles. This collision mechanism accounts for the macroscopic nature of dust particles compared to gas molecules. In a collision between a particle and a molecule, the post-collision relative velocity is drawn from a half-Maxwellian distribution at the temperature of the particle's surface. This assumption implies that a molecule coming into contact with a particle will thermalize with the molecules forming the particle's surface within a negligible time compared to other characteristic time scales. Consequently, the post-collisional relative velocity lies within the half-space defined by the tangent plane to the surface of the particle (see Figure 1).

<figure><img src="image_2.png" /><figcaption>Figure 1: The diffuse reflection mechanism considered in the article.</figcaption></figure>

Collision operators describing this mechanism were derived in [12] for the case of a fixed particle surface temperature $T_p$, assumed to be independent of time and space. However, this model is not entirely satisfactory, as it does not conserve the total energy of the system. To address this issue, an improved model was introduced in [16], incorporating the temperature change of a dust particle during a collision. This was achieved by treating the surface temperature of the dust particles as a variable in the density function $F$. Correspondingly, modified collision operators were derived based on this mechanism. However, deriving an entropy function for this model remains both a challenging task and an open problem.

Here, we propose an alternative model in which the surface temperature of the dust particles is treated as a function of both $t$ and $x$. In this model, the conservation of total energy is ensured
through an additional equation governing this temperature. We assume that at position $x$ and time $t$ the particles have the same surface temperature $T_p(t, x) > 0$, and, as in previous models, the post-collision relative velocity is given by a diffuse reflection at temperature $T_p(t, x)$. Let $m_g \in \mathbb{R}_+^*$ and $m_p \in \mathbb{R}_+^*$ denote the masses of the gas molecules and the particles, respectively. The collision mechanism between dust and gas thus leads to the introduction, for $n \in \mathbb{S}^2$, of the probability density $h_n$ of the post-collision relative velocities, given by

$$h_n(t, x, z) = \frac{1}{2\pi} \frac{m_g^2}{k_B^2 T_p(t, x)^2} (n \cdot z) e^{-\frac{m_g|z|^2}{2k_B T_p(t, x)}} 1_{\{z \cdot n \ge 0\}}, \quad (8)$$

where $k_B$ is the Boltzmann constant. The expression of the moments of $h_n$ are recalled in Appendix A.

The expression of the collision operators derived in [13] can be directly extended to the context of a time- and space-dependent surface temperature of the particles. Thus, the operators $\mathcal{D}$ and $\mathcal{R}$ can be written in the following form (where, for simplicity, we omit the dependence on $t$ and $x$ in $f$ and $F$):

$$\begin{aligned} \mathcal{D}(F, f)(t, x, v) = & \sigma_{gp} \left( \frac{m_g + m_p}{m_g} \right)^3 \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} \int_{\mathbb{S}^2} f(w^{\diamond}) F(v^{\diamond}) h_n \left( t, x, \frac{m_g + m_p}{m_g} (v_B - v) \right) \zeta(v^{\diamond} - w^{\diamond}, n) dndv^{\diamond} dw^{\diamond} \\ & - \sigma_{gp} \int_{\mathbb{R}^3} \int_{\mathbb{S}^2} f(w) F(v) \zeta(v - w, n) dndw, \end{aligned} \quad (9)$$

and

$$\begin{aligned} \mathcal{R}(F, f)(t, x, w) = & \sigma_{gp} \left( \frac{m_g + m_p}{m_p} \right)^3 \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} \int_{\mathbb{S}^2} f(w^{\diamond}) F(v^{\diamond}) h_n \left( t, x, \frac{m_g + m_p}{m_p} (v_B - w) \right) \zeta(w^{\diamond} - v^{\diamond}, n) dndv^{\diamond} dw^{\diamond} \\ & - \sigma_{gp} \int_{\mathbb{R}^3} \int_{\mathbb{S}^2} f(w) F(v) \zeta(w - v, n) dndv, \end{aligned} \quad (10)$$

where the collision kernel $\sigma_{gp}$ is of order of magnitude $(r_p + r_g)^2$. The quantity

$$v_B = \frac{m_g w^{\diamond} + m_p v^{\diamond}}{m_g + m_p} \quad (11)$$

is the velocity of the barycenter (which is a collision invariant with the assumptions made) and $\zeta(\cdot, \cdot)$ encodes the interspecies collision kernel. For the general case, this collision kernel has the following form:

$$\zeta(v - w, n) = b \left( |v - w|, \frac{(v - w) \cdot n}{|v - w|} \right), \quad (12)$$

where $b$ is a given function.

We aim to keep our analysis as general as possible regarding the choice of the collision kernel. However, obtaining explicit results requires specifying the collision kernel. In what follows, we primarily use the classical hard-sphere collision kernel, which models particle interactions as billiard ball collisions:

$$\zeta(v - w, n) = [(v - w) \cdot n]_+, \quad (13)$$

where we denote $s_+ = s1_{\{s \ge 0\}}$.

Another popular choice is the Maxwellian collision kernel, a special case where the collision kernel depends only on the deviation angle and not on the relative velocity:

$$\zeta(v - w, n) = b \left( \frac{(v - w) \cdot n}{|v - w|} \right). \quad (14)$$

This simplification enables explicit calculations that, in some cases, closely match physical observations [32]. Some results for this class of collision kernels are discussed in Appendix E.
In the following, we define

$$ \Sigma(v-w) := \int_{S^2} \zeta(v-w, n) \mathrm{d}n, \qquad (15) $$

the integral of the collision kernel $\zeta(\cdot, n)$ with respect to $n$.

We first establish the following weak formulations of the operators $\mathcal{D}(F, f)$ and $\mathcal{R}(F, f)$.

**Proposition 1.** Let $\varphi \in L_{\mathrm{loc}}^{\infty}(\mathbb{R}^3)$ and $\psi \in L_{\mathrm{loc}}^{\infty}(\mathbb{R}^3)$ be two test functions, and $\mathcal{D}(F, f)$ and $\mathcal{R}(F, f)$ the operators defined in (9) and (10). Then

$$ \int_{\mathbb{R}^3} \mathcal{D}(F, f)\varphi(v)dv = \sigma_{gp} \iiint\int (\varphi(v') - \varphi(v)) F(v)f(w)h_n(z)\zeta(v-w, n) \mathrm{d}ndz\mathrm{d}w\mathrm{d}v, \quad (16) $$

where

$$ v' = \frac{m_g w + m_p v}{m_g + m_p} - \frac{m_g}{m_g + m_p} z, \qquad (17) $$

and

$$ \int_{\mathbb{R}^3} \mathcal{R}(F, f)(w)\psi(w)dw = \sigma_{gp} \iiint\int (\psi(w') - \psi(w)) F(v)f(w)h_n(z)\zeta(v-w, n) \mathrm{d}ndz\mathrm{d}w\mathrm{d}v, \quad (18) $$

where

$$ w' = \frac{m_g w + m_p v}{m_g + m_p} + \frac{m_p}{m_g + m_p} z. \qquad (19) $$

*Proof.* Thanks to (9) and to the change of variable

$$ v \to z = \frac{m_g + m_p}{m_g} (v_B - v) $$

in the first integral, where $v_B$ is defined by (11), with Jacobian $((m_g + m_p)/m_g)^3$, we have:

$$ \begin{aligned} & \int_{\mathbb{R}^3} \mathcal{D}(F, f)(t, x, v)\varphi(v)dv \\ &= \sigma_{gp} \left(\frac{m_g + m_p}{m_g}\right)^3 \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} \int_{S^2} f(w^{\diamond}) F(v^{\diamond}) h_n\left(t, x, \frac{m_g + m_p}{m_g}(v_B - v)\right) \zeta(v^{\diamond} - w^{\diamond}, n)\varphi(v) \mathrm{d}ndv^{\diamond} dw^{\diamond} dv \\ & - \sigma_{gp} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} \int_{S^2} f(w) F(v) \zeta(v - w, n)\varphi(v) \mathrm{d}ndw\mathrm{d}v \\ &= \sigma_{gp} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} \int_{S^2} f(w^{\diamond}) F(v^{\diamond}) h_n(z) \zeta(v^{\diamond} - w^{\diamond}, n)\varphi(v') \mathrm{d}ndv^{\diamond} dw^{\diamond} dz \\ & - \sigma_{gp} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} \int_{S^2} f(w) F(v) \zeta(v - w, n) h_n(t, x, z) \varphi(v) \mathrm{d}z\mathrm{d}vdz\mathrm{d}ndw, \end{aligned} $$

where we used the fact that

$$ \int h_n(t, x, z) \mathrm{d}z = 1 $$

for the second integral. We then obtain (16). Moreover, the change of variable

$$ (w, n) \to \left( z = -\frac{m_g + m_p}{m_p} (v_B - w), -n \right) $$
in the first integral, and $n \to -n$ in the second integral, yield

$$
\begin{align*}
& \int_{\mathbb{R}^3} \mathcal{R}(F, f)(t, x, w)\psi(w)dw \\
&= \sigma_{gp} \left(\frac{m_g + m_p}{m_p}\right)^3 \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} \int_{S^2} f(w^{\diamond}) F(v^{\diamond}) h_n\left(\frac{m_g + m_p}{m_p}(v_B - w)\right) \varsigma(w^{\diamond} - v^{\diamond}, n)\psi(w) dndv^{\diamond}dw^{\diamond} \\
&\quad - \sigma_{gp} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} \int_{S^2} f(w) F(v) \varsigma(w - v, n)\psi(w) dndvdw \\
&= \sigma_{gp} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} \int_{S^2} f(w^{\diamond}) F(v^{\diamond}) h_{-n}(t, x, -z) \varsigma(w^{\diamond} - v^{\diamond}, -n)\psi(w) dndv^{\diamond}dw^{\diamond} \\
&\quad - \sigma_{gp} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} \int_{S^2} f(w) F(v) \varsigma(w - v, -n)\psi(w) dndvdw
\end{align*}
$$

Expression (18) is obtained by observing in (8) that $h_{-n}(t, x, -z) = h_n(t, x, z)$ and $\varsigma(w - v, -n) = \varsigma(v - w, n)$. $\square$

For deriving another weak formulation of the collision operators, the following auxiliary Lemma is crucial, and it is also essential for obtaining their strong forms.

**Lemma 1.** Let $\varphi \in L_{\text{loc}}^{\infty}(\mathbb{R}^3)$ and $\psi \in L_{\text{loc}}^{\infty}(\mathbb{R}^3)$ be two test functions, and $\mathcal{D}(F, f)$ and $\mathcal{R}(F, f)$ the operators defined in (9) and (10). Then

$$
\int_{\mathbb{R}^3} \mathcal{D}(F, f)\varphi(v)dv = \sigma_{gp} \iiint\int \varphi(v) [F(v')f(w')h_n(v-w)\varsigma(z,n) - F(v)f(w)h_n(z)\varsigma(v-w,n)] dndzdwdv, \quad (20)
$$

and

$$
\int_{\mathbb{R}^3} \mathcal{R}(F, f)(w)\psi(w)dw = \sigma_{gp} \iiint \psi(w) [F(v')f(w')h_n(v-w)\varsigma(z,n) - F(v)f(w)h_n(z)\varsigma(v-w,n)] dndzdwdv, \quad (21)
$$

where $v'$ and $w'$ are given by (17) and (19), respectively. Moreover, the following equalities hold:

$$
\begin{align}
& \iiint\int \varphi(v')F(v')f(w')h_n(v-w)\varsigma(z,n) \, dndz\,dw\,dv \tag{22} \\
& \qquad = \iiint\int \varphi(v)F(v)f(w)h_n(z)\varsigma(v-w,n) \, dndz\,dw\,dv, \nonumber
\end{align}
$$

and

$$
\begin{align}
& \iiint\int \psi(w')F(v')f(w')h_n(v-w)\varsigma(z,n) \, dndz\,dw\,dv \tag{23} \\
& \qquad = \iiint\int \psi(w)F(v)f(w)h_n(z)\varsigma(v-w,n) \, dndz\,dw\,dv. \nonumber
\end{align}
$$

*Proof.* In the gain term of (16), we perform the change of variables $(v, w, z, n) \to (v', w', u = w - v, -n)$, where $v'$ and $w'$ are given by (17) and (19) respectively. This transformation is an involution, since $w' - v' = z$, and consequently has a Jacobian equal to 1. Hence, we obtain:

$$
\begin{align*}
& \iiint\int \varphi(v')F(v)f(w)h_n(z)\varsigma(v-w,n) \, dndz\,dw\,dv \\
&= \iiint\int \varphi(v)F(v')f(w')h_{-n}(w-v)\varsigma(-z,-n) \, dndz\,dw\,dv.
\end{align*}
$$

From (12), we see that $\varsigma(-z, -n) = \varsigma(z, n)$, and, from (8), that $h_{-n}(w - v) = h_n(v - w)$. This allows us to obtain (20). We then apply the change of variables $(v, w, z, n) \to (v', w', u = w - v, -n)$
in the gain term of (18):

$$ 
\begin{aligned}
& \iiint \psi(w')F(v)f(w)h_n(z)\zeta(v-w, -n) \, dndzdwdv \\
&= \iiint \psi(w)F(v')f(w')h_{-n}(w-v)\zeta(-z, -n) \, dndzdwdv \\
&= \iiint \psi(w)F(v')f(w')h_n(v-w)\zeta(z, n) \, dndzdwdv
\end{aligned}
$$

to get (21). The same change of variables, together with the same observations gives (22) and (23) as well. □

**Corollary 1.** Let $\varphi \in L_{\text{loc}}^{\infty}(\mathbb{R}^3)$ and $\psi \in L_{\text{loc}}^{\infty}(\mathbb{R}^3)$ be two test functions, and $\mathcal{D}(F, f)$ and $\mathcal{R}(F, f)$ the operators defined in (9) and (10). Then

$$ 
\begin{aligned}
& \int_{\mathbb{R}^3} \mathcal{D}(F, f)\varphi(v)dv + \int_{\mathbb{R}^3} \mathcal{R}(F, f)\psi(w)dw \\
&= -\frac{\sigma_{gp}}{2} \iiint \left(\varphi(v') + \psi(w') - \varphi(v) - \psi(w)\right) \\
& \qquad \times [F(v')f(w')h_n(v-w)\zeta(z, n) - F(v)f(w)h_n(z)\zeta(v-w, n)] dndzdwdv,
\end{aligned}
\quad (24)
$$

where $v'$ and $w'$ are given by (17) and (19) respectively.

*Proof.* The proof is an immediate consequence of Lemma 1 together with (16)–(21). Indeed, we consider separately $\int_{\mathbb{R}^3} \mathcal{D}(F, f)\varphi(v)dv$ and $\int_{\mathbb{R}^3} \mathcal{R}(F, f)\psi(w)dw$. For the first term, we split the integral into two parts, rewriting one with (16) and the other with (20). Analogously, for the second term we apply (18) and (21). Summing the resulting equalities and then using (22)–(23) yields (24). □

In [13], strong forms of the operators $\mathcal{D}$ and $\mathcal{R}$ were established for the hard-sphere collision kernel. Here we extend these expressions to the case of a general collision kernel $\zeta(v-w, n)$.

**Proposition 2.** The operators $\mathcal{D}(F, f)$ and $\mathcal{R}(F, f)$, given by (9) and (10), can be expressed as follows (where we omit the $t$ and $x$ variables in order to simplify the notation):

$$ 
\mathcal{D}(F, f)(v) = \sigma_{gp} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} \int_{S^2} [f(w')F(v')h_n(v-w)\zeta(z, n) - f(w)F(v)h_n(z)\zeta(v-w, n)] \, dzdndw, \quad (25)
$$

and

$$ 
\mathcal{R}(F, f)(w) = \sigma_{gp} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} \int_{S^2} [F(v')f(w')h_n(v-w)\zeta(z, n) - F(v)f(w)h_n(z)\zeta(v-w, n)] \, dzdndv, \quad (26)
$$

where $v'$ and $w'$ are given by (17) and (19) respectively.

*Proof.* It is a direct consequence of equations (20) and (21), respectively, derived in Lemma 1. □

**Remark 1.** Note that the total momentum is conserved during a collision. If $v'$ and $w'$ are given by (17) and (19), then

$$ m_p v' + m_g w' = m_p v + m_g w. \quad (27) $$

However, in contrast to the classical collision mechanism found in Boltzmann operators, the total kinetic energy is not conserved during a collision:

$$ m_p|v'|^2 + m_g|w'|^2 \neq m_p|v|^2 + m_g|w|^2. $$
Nevertheless, the difference between the kinetic energy after the collision and before the collision can be expressed in terms of the square of the relative velocities before and after the collision, as established in the following Lemma.

**Lemma 2.** Let $v'$ and $w'$ the post-collision velocities, given by (17) and (19) respectively. Then we have the following relation

$$m_p|v'|^2 - m_p|v|^2 + m_g|w'|^2 - m_g|w|^2 = \frac{m_p m_g}{m_p + m_g} (|z|^2 - |v-w|^2), \quad (28)$$

*Proof.* If we replace $(m_p v + m_g w)/(m_p + m_g)$ by $v_B$ in (17) and (19), we get:

$$\begin{aligned} m_p|v'|^2 + m_g|w'|^2 &= m_p \left|v_B - \frac{m_g}{m_p + m_g}z\right|^2 + m_g \left|v_B + \frac{m_p}{m_g + m_p}z\right|^2 \\ &= (m_g + m_p)|v_B|^2 + \frac{m_p m_g}{m_p + m_g}|z|^2. \end{aligned} \quad (29)$$

Using the involutive character of transformation $(v, w, z) \rightarrow (v', w', w - v)$, and the conservation of the momentum (27), we deduce from (29) that

$$m_p|v|^2 + m_g|w|^2 = (m_p + m_g)|v_B|^2 + \frac{m_p m_g}{m_p + m_g}|w - v|^2.$$

By subtracting the latter relation and relation (29), we obtain (28). □

We introduce the ratio

$$\mu_{gp} = \frac{m_p + m_g}{m_p}, \quad (30)$$

and derive an alternative weak form of the operators for the hard-sphere collision kernel. This formulation will be utilized to establish the entropy inequality in Section 3.

**Proposition 3.** In the case of a hard-sphere collision kernel, that is $\zeta(v-w, n) = [(v-w) \cdot n]_+$, we have the following weak form, for $\varphi \in L_{\text{loc}}^{\infty}(\mathbb{R}^3)$ and $\psi \in L_{\text{loc}}^{\infty}(\mathbb{R}^3)$ test functions:

$$\begin{aligned} & \int_{\mathbb{R}^3} \mathcal{D}(F, f)\varphi(v)dv + \int_{\mathbb{R}^3} \mathcal{R}(F, f)\psi(w)dw \\ &= -\frac{\sigma_{gp}}{2} \iiint \int \left[\varphi(v') + \psi(w') - \varphi(v) - \psi(w)\right] [(v-w) \cdot n]_+ h_n(z) e^{-\frac{\mu_{gp}}{2k_B T_p} (m_p|v'|^2 + m_g|w'|^2)} \\ & \times \left[ F(v')e^{\frac{\mu_{gp}}{2k_B T_p} m_p|v'|^2} f(w')e^{\frac{\mu_{gp}}{2k_B T_p} m_g|w'|^2} - F(v) e^{\frac{\mu_{gp}}{2k_B T_p} m_p|v|^2} f(w) e^{\frac{\mu_{gp}}{2k_B T_p} m_g|w|^2} \right] dndzdwdv. \end{aligned} \quad (31)$$

*Proof.* From (24) and the definition (8) of $h_n$, we get

$$\begin{aligned} & \int_{\mathbb{R}^3} \mathcal{D}(F, f)\varphi(v)dv \\ &= -\frac{\sigma_{gp}}{2} \iiint \int (\varphi(v') - \varphi(v)) [(v-w) \cdot n]_+ h_n(z) \\ & \times [F(v')f(w')e^{-\frac{m_g}{2k_B T_p} |v-w|^2} e^{\frac{m_g}{2k_B T_p} |z|^2} - F(v)f(w)] dndzdwdv, \\ &= -\frac{\sigma_{gp}}{2} \iiint \int (\varphi(v') - \varphi(v)) [(v-w) \cdot n]_+ h_n(z) \\ & \times [F(v')f(w')e^{\frac{1}{2k_B T_p} \frac{m_g+m_p}{m_p} (m_p|v'|^2 - m_p|v|^2 + m_g|w'|^2 - m_g|w|^2)} - F(v)f(w)] dndzdwdv, \end{aligned}$$

thanks to (28). The proof for the part $\int_{\mathbb{R}^3} \mathcal{R}(F, f)\psi(w)dw$ follows similarly. Finally, by using the definition of $\mu_{gp}$ given in Equation (30), we deduce (31). □
**Remark 2.** The fact that the term $\int_{\mathbb{R}^3} \mathcal{D}(F, f)\varphi(v)dv + \int_{\mathbb{R}^3} \mathcal{R}(F, f)\psi(w)dw$ can be put into the form

$$\iiint \int [\varphi(v') + \psi(w') - \varphi(v) - \psi(w)] [G(v')H(w') - G(v)H(w)] r(v, w, z, n)dndzdwdv,$$

with some non-negative functions $r$, $G$, $H$ – which is the case of expression (31) – seems to be specific to the case of the hard-sphere collision kernel. For general collision kernels, this simple factorization does not hold because the symmetry between pre- and post-collisional velocities is lost. More precisely, using the relation on the velocities (28) and that $z = w' - v'$, one obtains

$$F(v')f(w')h_n(v-w)\zeta(z,n) - F(v)f(w)h_n(z)\zeta(v-w,n) \\ = \frac{1}{2\pi} \frac{m_g^2}{k_B^2 T_p^2} \left[ F(v')e^{\frac{\mu_{gp}}{2k_B T_p} m_p |v'|^2} f(w')e^{\frac{\mu_{gp}}{2k_B T_p} m_g |w'|^2} [(v-w)\cdot n]_+ \zeta(z,n) \right. \\ \left. - F(v)e^{\frac{\mu_{gp}}{2k_B T_p} m_p |v'|^2} f(w)e^{\frac{\mu_{gp}}{2k_B T_p} m_g |w'|^2} [z\cdot n]_+ \zeta(v-w,n) \right] e^{-\frac{m_g|z|^2}{2k_B T_p}} e^{-\frac{\mu_{gp}}{2k_B T_p} (m_p|v|^2 + m_g|w|^2)}.$$

In the case of a hard-sphere collision kernel, one has

$$[(v-w) \cdot n]_+ \zeta(z, n) = [z \cdot n]_+ \zeta(v-w, n),$$

which allows to factor this term and then write the entire expression as

$$[G(v')H(w') - G(v)H(w)] r(v, w, z, n).$$

This is no longer the case for a different collision kernel. This property will be useful for establishing an inequality on the entropy dissipation term (see Subsection 3.2).

## 2.4 Surface temperature evolution

In order to model the evolution of the surface temperature of the particles, we first introduce a key quantity associated with the system: its total energy. At time $t$ and position $x$, the energy of the system is defined as follows:

$$E_T(t, x) = \frac{m_p}{2} \int_{\mathbb{R}^3} F(t, x, v) |v|^2 \, dv + \frac{m_g}{2} \int_{\mathbb{R}^3} f(t, x, w) |w|^2 \, dw + m_p c_p T_p(t, x) \int_{\mathbb{R}^3} F(t, x, v) \, dv, \quad (32)$$

which correspond to the sum of the kinetic energies of each species and the internal energy of the particles, with $c_p$ being the particle mass heat capacity. We are interested in the time evolution of the total energy.

**Proposition 4.** Let $(F, f)$ satisfy the system (1) in $\mathbb{R}^3$, with $\mathcal{D}(F, f)$ and $\mathcal{R}(F, f)$ given by (9) and (10), and $\lim_{|x| \to +\infty} F(t, x, v) = \lim_{|x| \to +\infty} f(t, x, w) = 0$. If $(F, f, T_p)$ are sufficiently smooth, then, at a formal level, the time evolution of the total energy of the system satisfies

$$\frac{d}{dt} \int_{\mathbb{R}^3} E_T(t, x) \, dx = m_p \int_{\mathbb{R}^3} \left[ c_p \left( \partial_t T_p(t, x) \int_{\mathbb{R}^3} F(t, x, v) \, dv + \nabla_x T_p(t, x) \cdot \int_{\mathbb{R}^3} F(t, x, v) v \, dv \right) \right. \\ \left. + 2k_b \Upsilon_1(F, f)(t, x) T_p(t, x) - \frac{1}{2} m_g \Upsilon_2(F, f)(t, x) \right] \, dx, \quad (33)$$

where

$$\Upsilon_1(F, f)(t, x) := \sigma_{gp} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} F(t, x, v) f(t, x, w) \Sigma(v-w) dv dw, \quad (34)$$

$$\Upsilon_2(F, f)(t, x) := \sigma_{gp} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} F(t, x, v) f(t, x, w) \Sigma(v-w) |v-w|^2 dv dw. \quad (35)$$
Proof. Let us write

$$ \partial_t E_T(t, x) = \partial_t E_K(t, x) + \partial_t E_I(t, x), $$

where $E_K$ and $E_I$ denote, respectively, the kinetic energies of both species and the internal energy of the particles:

$$ \begin{aligned} \partial_t E_K(t, x) &= \frac{m_p}{2} \partial_t \int_{\mathbb{R}^3} F(t, x, v) |v|^2 \, dv + \frac{m_g}{2} \partial_t \int_{\mathbb{R}^3} f(t, x, w) |w|^2 \, dw, \\ \partial_t E_I(t, x) &= m_p c_p \partial_t T_p(t, x) \int_{\mathbb{R}^3} F(t, x, v) \, dv + m_p c_p T_p(t, x) \partial_t \int_{\mathbb{R}^3} F(t, x, v) \, dv. \end{aligned} $$

Since 1 is a collisional invariant of $\mathcal{P}(F)$ and $\mathcal{D}(F, f)$, $F(t, x, v)$ satisfies the conservation equation

$$ \partial_t \int_{\mathbb{R}^3} F(t, x, v) \, dv + \nabla_x \cdot \left( \int_{\mathbb{R}^3} F(t, x, v) v \, dv \right) = 0. $$

Then

$$ \begin{aligned} \int_{\mathbb{R}^3} T_p(t, x) \partial_t \int_{\mathbb{R}^3} F(t, x, v) \, dv \, dx &= - \int_{\mathbb{R}^3} T_p(t, x) \nabla_x \cdot \left( \int_{\mathbb{R}^3} F(t, x, v) v \, dv \right) \, dx \\ &= \int_{\mathbb{R}^3} \nabla_x T_p(t, x) \cdot \int_{\mathbb{R}^3} F(t, x, v) v \, dv \, dx, \end{aligned} $$

and

$$ \int_{\mathbb{R}^3} \partial_t E_I(t, x) dx = m_p c_p \int_{\mathbb{R}^3} \left( \partial_t T_p(t, x) \int_{\mathbb{R}^3} F(t, x, v) \, dv + \nabla_x T_p(t, x) \cdot \int_{\mathbb{R}^3} F(t, x, v) v \, dv \right) dx. $$

Moreover, since $|v|^2$ and $|w|^2$ are collisional invariants of $\mathcal{P}(F)$ and $\mathcal{C}(f)$, respectively (see Equation (7)), and thanks to the weak formulations (16) and (18) of $\mathcal{D}(F, f)$ and $\mathcal{R}(F, f)$, we deduce that

$$ \begin{aligned} \partial_t E_K(t, x) &= \frac{m_p}{2} \int_{\mathbb{R}^3} (-v \cdot \nabla_x F + \mathcal{D}(F, f) + \mathcal{P}(F)) |v|^2 \, dv \\ &\quad + \frac{m_g}{2} \int_{\mathbb{R}^3} (-w \cdot \nabla_x f + \mathcal{R}(F, f) + \mathcal{C}(f)) |w|^2 \, dw \\ &= -\frac{m_p}{2} \nabla_x \cdot \left( \int_{\mathbb{R}^3} vF|v|^2 \, dv \right) - \frac{m_g}{2} \nabla_x \cdot \left( \int_{\mathbb{R}^3} wF|w|^2 \, dw \right) \\ &\quad + \frac{1}{2} \sigma_{gp} \iiint F(t, x, v) f(t, x, w) h_n(t, x, z) \zeta(v-w, n) \\ &\quad [m_p|v'|^2 - m_p|v|^2 + m_g|w'|^2 - m_g|w|^2] \, dndzdvdw. \end{aligned} $$

Thanks to the hypotheses on the tails of $F$ and $f$, we have

$$ \int_{\mathbb{R}^3} \nabla_x \cdot \left( \int_{\mathbb{R}^3} vF|v|^2 \, dv \right) \, dx = 0 \quad \text{and} \quad \int_{\mathbb{R}^3} \nabla_x \cdot \left( \int_{\mathbb{R}^3} wF|w|^2 \, dw \right) \, dx = 0. $$

Moreover, using (28) together with the following properties of $h_n$

$$ \int_{\mathbb{R}^3} h_n(t, x, z) dz = 1, \quad \int_{\mathbb{R}^3} h_n(t, x, z) |z|^2 dz = \frac{4k_b}{m_g} T_p(t, x), $$

(see Equations (110) and (112) for $h_n$, with $c = \sqrt{\frac{m_g}{k_b T_p(t, x)}}$ in Appendix A), we can compute the second part of the term $\partial_t E_K(t, x)$:

$$ \begin{aligned} & \frac{1}{2} \sigma_{gp} \iiint F(v) f(w) h_n(z) \zeta(v-w, n) [m_p|v'|^2 - m_p|v|^2 + m_g|w'|^2 - m_g|w|^2] \, dndzdvdw \\ &= \frac{m_p}{m_p + m_g} \sigma_{gp} \iint F(t, x, v) f(t, x, w) \Sigma(v-w) \left( 2k_b T_p(t, x) - \frac{1}{2} m_g |v-w|^2 \right) \, dvdw. \end{aligned} $$

We finally get (33).

□
The evolution of the surface temperature of the particles $T_p(t, x)$ is then described by the following PDE, which ensures the conservation of the total energy of the system :

$$c_p \left( \partial_t T_p(t, x) \int_{\mathbb{R}^3} F(t, x, v) \, dv + \nabla_x T_p(t, x) \cdot \int_{\mathbb{R}^3} vF(t, x, v) \, dv \right) + \frac{2k_b}{m_p + m_g} \Upsilon_1(F, f)(t, x)T_p(t, x) = \frac{1}{2} \frac{m_g}{m_p + m_g} \Upsilon_2(F, f)(t, x), \quad (36)$$

where $\Upsilon_1(F, f)$ and $\Upsilon_2(F, f)$ are defined by (34) and (35). We immediately see that the system conserves the total energy if $T_p$ satisfies (36):

**Corollary 2.** Let $(F, f, T_p)$ satisfy the system (1) and (36) in $\mathbb{R}^3$, with $\mathcal{D}(F, f)$ and $\mathcal{R}(F, f)$ given by (9) and (10), and $\lim_{|x| \to +\infty} F(t, x, v) = \lim_{|x| \to +\infty} f(t, x, w) = 0$. If $(F, f, T_p)$ are sufficiently smooth, then, at a formal level, the time evolution of the total energy of the system satisfies

$$\frac{d}{dt} \int_{\mathbb{R}^3} E_T(t, x) \, dx = 0. \quad (37)$$

# 3 Collision invariants, entropy and equilibrium states

## 3.1 Collision invariants

We recall that property (7) of conservation of mass, momentum and energy holds for operators $\mathcal{C}$ and $\mathcal{P}$. For operators $\mathcal{D}$ and $\mathcal{R}$, we deduce the conservation of mass and momentum from Corollary 1.

**Corollary 3.** The operators $\mathcal{D}$ and $\mathcal{R}$ defined by (9) and (10) satisfy

$$\int_{\mathbb{R}^3} \mathcal{R}(F, f)(w)dw = \int_{\mathbb{R}^3} \mathcal{D}(F, f)(v)dv = 0, \quad (38)$$

and

$$\eta \int_{\mathbb{R}^3} \mathcal{R}(F, f)(w)wdw + \int_{\mathbb{R}^3} \mathcal{D}(F, f)(v)vdv = 0. \quad (39)$$

*Proof.* The starting point is given by the weak formulations (16) and (18) of the operators $\mathcal{D}$ and $\mathcal{R}$, respectively. Thanks to the pre-post interaction rule (17) and (19), by taking $\varphi = \psi = 1$ and $\varphi(v) = v$ and $\psi(w) = \eta w$, the thesis of the Corollary follows immediately. $\square$

**Remark 3.** The collisions between a gas molecule and a dust particle do not conserve the kinetic energy, as a consequence of Remark 1.

## 3.2 Entropy decay

A key quantity associated with kinetic collisional equations is the entropy. As recalled in [32], for the classical Boltzmann equation in the whole space

$$\partial_t f(t, x, w) + w \cdot \nabla_x f(t, x, w) = \mathcal{C}(f)(t, x, w), \quad t \ge 0, x \in \mathbb{R}^3, w \in \mathbb{R}^3 \quad (40)$$

with $\mathcal{C}(f)$ given by (2), the entropy is defined by

$$H_B(f)(t) = \int_{\mathbb{R}^3} h(t, x)dx, \quad \text{with} \quad h(t, x) = \int_{\mathbb{R}^3} f(t, w) \ln f(t, w)dw.$$

For $f$ a smooth solution of (40) such that $\lim_{|x| \to +\infty} f(t, x, w) = 0$, we have :

$$\frac{dH_B(f)}{dt}(t) = \int_{\mathbb{R}^3} D_B(f, f)(t, x)dx$$
where $D_B(f, f)$ is the entropy dissipation function, defined by

$$ 
\begin{aligned}
D_B(f, f)(t, x) &= \int_{\mathbb{R}^3} \mathcal{C}(f)(w) \ln f(w) dw \\
&= -\frac{\sigma_{gg}}{4} \iiint [f'f'_* - ff_*] [\ln(f'f'_*) - \ln(ff_*)] c(w - w_*, n) dndw_* dw.
\end{aligned}
$$

The algebraic inequality $(y_1 - y_2)(\ln(y_1) - \ln(y_2)) \ge 0$ for all $(y_1, y_2) \in \mathbb{R}^{+*} \times \mathbb{R}^{+*}$ implies that $D_B(f, f)(t, x) \le 0$. The Boltzmann entropy is, therefore, a non-increasing function of time.

We define here an entropy function for the system (1). Let

$$ H(F, f)(t) = \int h_1(t, x) dx + \int h_2(t, x) dx, \qquad (41) $$

with

$$ h_1(t, x) = \int f(t, x, w) \ln f(t, x, w) dw + \int F(t, x, v) \ln F(t, x, v) dv, $$

and

$$ h_2(t, x) = -\frac{m_g + m_p}{k_B} c_p \ln T_p(t, x) \int F(t, x, v) dv. $$

Here $h = h_1 + h_2$ define the *entropy density*. In the following Proposition, we show that the functional $H(F, f)$ has the property that characterizes it as an entropy for System (1) in the case of hard-sphere binary interactions.

**Proposition 5.** Let $(f, F)$ a smooth solution of (1), with

$$ \lim_{|x| \to +\infty} F(t, x, v) = \lim_{|x| \to +\infty} f(t, x, w) = 0, $$

and assume that the collision kernels of $\mathcal{R}(F, f)$ and $\mathcal{D}(F, f)$ are hard-sphere collision kernels given by (13). Then the functional $H(F, f)$, defined in (41), does not increase with time.

To show this result, we start by defining the *entropy dissipation* functional:

$$ D(F, f)(t, x) = D^d(F, f)(t, x) + D^e(F, f)(t, x), \qquad (42) $$

where

$$ 
\begin{aligned}
D^d(F, f)(t, x) &= \int \mathcal{R}(F, f)(t, x, w) \ln \left( f(t, x, w) e^{\frac{\mu_{gp}}{2k_B T_p(t,x)} m_g |w|^2} \right) dw \\
&\quad + \int \mathcal{D}(F, f)(t, x, v) \ln \left( F(t, x, v) e^{\frac{\mu_{gp}}{2k_B T_p(t,x)} m_p |v|^2} \right) dv, \\
&= \int \mathcal{R}(F, f)(t, x, w) \ln f(t, x, w) dw + \int \mathcal{D}(F, f)(t, x, v) \ln F(t, x, v) dv \\
&\quad + \frac{\mu_{gp}}{2k_B T_p(t, x)} \left( m_g \int \mathcal{R}(F, f)(t, x, w) |w|^2 dw + m_p \int \mathcal{D}(F, f)(t, x, v) |v|^2 dv \right)
\end{aligned}
\qquad (43) $$

is the entropy dissipation associated to the diffuse reflection mechanism, and

$$ D^e(F, f)(t, x) = \int \mathcal{C}(f)(t, x, w) \ln f(t, x, w) dw + \int \mathcal{P}(F)(t, x, v) \ln F(t, x, v) dv \qquad (44) $$

is the entropy dissipation associated with the elastic collisions.

**Lemma 3.** Assume that $\zeta$ is given by (13). The entropy dissipation (42) is such that

$$ D(F, f)(t, x) \le 0, \quad \text{for all } (t, x) \in \mathbb{R}^+ \times \mathbb{R}^3. $$
Proof. To show that the dissipation functional $D(F, f)$ is negative, we study the sign of $D^e(F, f)$ and $D^d(F, f)$ separately. The negativity of the dissipation functional $D^e(F, f)$ associated with the Boltzmann operators is well-known. By choosing $\psi(w) = \ln(f(w))$ and $\varphi(v) = \ln(F(v))$ in the weak formulation (6), we obtain:

$$ D^e(F, f)(t, x) = -\frac{\sigma_{gg}}{4} \iiint [f'f'_* - ff_*] [\ln(f'f'_*) - \ln(ff_*)] c(w - w_*, n) dndw_*dw \\ -\frac{\sigma_{pp}}{4} \iiint [F'F'_* - FF_*] [\ln(F'F'_*) - \ln(FF_*)] p(v - v_*, n) dndv_*dv. $$

The algebraic inequality $(y_1 - y_2)(\ln(y_1) - \ln(y_2)) \ge 0$ then shows that $D^e(F, f)(t, x) \le 0$. Then, by using

$$ \varphi(t, x, v) = \ln \left( F(t, x, v) e^{\frac{\mu_{gp}}{2k_B T_p(t,x)} m_p |v|^2} \right) \quad \text{and} \quad \psi(t, x, w) = \ln \left( f(t, x, w) e^{\frac{\mu_{gp}}{2k_B T_p(t,x)} m_g |w|^2} \right) $$

in the weak formulation (31), we get:

$$ D^d(F, f)(t, x) = -\frac{\sigma_{gp}}{2} \iiint \left[ (v-w) \cdot n \right]_+ h_n(z) e^{-\frac{\mu_{gp}}{2k_B T_p} \left(m_p|v|^2 + m_g|w|^2\right)} \\ \times \left[ \ln(F(v')e^{\frac{\mu_{gp}}{2k_B T_p} m_p |v'|^2} f(w')e^{\frac{\mu_{gp}}{2k_B T_p} m_g |w'|^2}) - \ln(F(v)e^{\frac{\mu_{gp}}{2k_B T_p} m_p |v|^2} f(w)e^{\frac{\mu_{gp}}{2k_B T_p} m_g |w|^2}) \right] \\ \times \left[ F(v')e^{\frac{\mu_{gp}}{2k_B T_p} m_p |v'|^2} f(w')e^{\frac{\mu_{gp}}{2k_B T_p} m_g |w'|^2} - F(v)e^{\frac{\mu_{gp}}{2k_B T_p} m_p |v|^2} f(w)e^{\frac{\mu_{gp}}{2k_B T_p} m_g |w|^2} \right] dndzdwdv \le 0. $$

□

**Remark 4.** The weak formulation (31), which is written as

$$ \iiint \left[ \varphi(v') + \psi(w') - \varphi(v) - \psi(w) \right] \left[ G(v')H(w') - G(v)H(w) \right] r(v, w, z, n) dndzdwdv, $$

with some non-negative function $G$, $H$, $r$, and which is specific to the hard-sphere collision kernel (see Remark 2), is the key point for obtaining the inequality $D^d(F, f)(t, x) \le 0$. This suggests that the proof of the entropy dissipation inequality might be closely tied to the specific structure of the hard-sphere collision kernel. Whether similar inequalities hold for other kernels remains an open question, potentially requiring alternative approaches.

*Proof of Proposition 5.* We start considering $h_1(t, x)$, and we will omit the dependence of $f$ and $F$ on $(t, x)$ for simplifying the notation:

$$ \begin{aligned} \partial_t h_1(t, x) &= \int \partial_t f(w) (\ln f(w) + 1) dw + \int \partial_t F(v) (\ln F(v) + 1) dv \\ &= -\nabla_x \cdot \int w f(w) \ln f(w) dw - \nabla_x \cdot \int v F(v) \ln F(v) dv \\ &\quad + \int (\mathcal{R}(F, f)(w) + \mathcal{C}(f)(w)) \ln f(w) dw + \int (\mathcal{D}(F, f)(v) + \mathcal{P}(F)(v)) \ln F(v) dv, \end{aligned} $$

where we made use of the Boltzmann system (1) and Corollary 3. Using the hypothesis on the tails of both $F$ and $f$, we have

$$ \int \nabla_x \cdot \int w f(w) \ln f(w) dw dx = 0 \quad \text{and} \quad \int \nabla_x \cdot \int v F(v) \ln F(v) dv dx = 0. $$
Moreover, using the second line of (43), we deduce the following equality:

$$ \int \partial_t h_1(t, x)dx = \int D(F, f)(t, x)dx - \frac{\mu_{gp}}{k_B} \int \frac{1}{T_p(t, x)} \left( \frac{m_g}{2} \int \mathcal{R}(F, f)|w|^2 dw + \frac{m_p}{2} \int \mathcal{D}(F, f)|v|^2 dv \right) dx. $$

If now we use the weak formulations (16) and (18) for $\mathcal{D}(F, f)$ and $\mathcal{R}(F, f)$ respectively, we can proceed in the same way as in the proof of Proposition 2. It follows, thanks to (28) and properties (110) and (112) for $h_n$ (with $c = \sqrt{\frac{m_g}{k_B T_p(t, x)}}$; see Appendix A), that

$$ \begin{aligned} & \frac{m_g}{2} \int \mathcal{R}(F, f)|w|^2 dw + \frac{m_p}{2} \int \mathcal{D}(F, f)|v|^2 dv \\ &= \frac{\sigma_{gp}}{2} \iiint \int F(t, x, v)f(t, x, w)h_n(t, x, z)\zeta(v-w, n) [m_p|v'|^2 - m_p|v|^2 + m_g|w'|^2 - m_g|w|^2] \, dndzdvdw \\ &= \frac{\sigma_{gp}}{2} \frac{m_p m_g}{m_p + m_g} \iiint \int F(t, x, v)f(t, x, w)h_n(t, x, z)\zeta(v-w, n) [|z|^2 - |v-w|^2] \, dndzdvdw \\ &= \frac{\sigma_{gp}}{2} \frac{m_p m_g}{m_p + m_g} \left( \frac{4k_B T_p(t, x)}{m_g} \Upsilon_1(F, f)(t, x) - \Upsilon_2(F, f)(t, x) \right), \end{aligned} $$

where $\Upsilon_1(F, f)$ and $\Upsilon_2(F, f)$ are defined by (34) and (35) respectively. Then, thanks to (30) and the surface temperature equation (36), we obtain

$$ \begin{aligned} & \int \partial_t h_1(t, x)dx \\ &= \int D(F, f)(t, x)dx - \frac{\mu_{gp}}{2k_B} m_p \int \frac{1}{T_p(t, x)} \left( \frac{4k_B T_p(t, x)}{m_p + m_g} \Upsilon_1(F, f)(t, x) - \frac{m_g}{m_p + m_g} \Upsilon_2(F, f)(t, x) \right) dx \\ &= \int D(F, f)(t, x)dx + \frac{(m_g + m_p)}{k_B} c_p \int \frac{1}{T_p(t, x)} \left( \partial_t T_p(t, x) \int_{\mathbb{R}^3} F(t, x, v)dv + \nabla_x T_p(t, x) \int F(v)vdv \right) dx \\ &= \int D(F, f)(t, x)dx + \frac{(m_g + m_p)}{k_B} c_p \int \left( \partial_t (\ln(T_p(t, x))) \int F(v)dv + \nabla_x (\ln(T_p(t, x))) \cdot \int F(v)vdv \right) dx. \end{aligned} $$

Moreover, since 1 is a collision invariant of $\mathcal{D}(F, f)$ and $\mathcal{P}(F, f)$, from Equation (1), we deduce

$$ \int \partial_t F(t, x, v)dv + \int \nabla_x \cdot (F(t, x, v)v)dv = 0, $$

and thus

$$ \begin{aligned} & \int \partial_t h_2(t, x)dx \\ &= -\frac{m_g + m_p}{k_B} c_p \int \partial_t \left( \ln T_p(t, x) \int F(t, x, v) \right) dv dx \\ &= -\frac{m_g + m_p}{k_B} c_p \int \left( \partial_t (\ln T_p(t, x)) \int F(t, x, v) dv - \ln T_p(t, x) \int \nabla_x \cdot (F(t, x, v)v) dv \right) dx \\ &= -\frac{m_g + m_p}{k_B} c_p \left( \int \partial_t (\ln T_p(t, x)) \int F(t, x, v) dv dx + \int \nabla_x (\ln T_p(t, x)) \cdot \int F(t, x, v) v dv dx \right), \end{aligned} $$

from the divergence theorem. We finally obtain

$$ \int \partial_t h_1(t, x)dx + \int \partial_t h_2(t, x)dx = \int D(F, f)(t, x)dx, $$

and so, by Lemma 3, the thesis.

□
### 3.3 Equilibrium states

We first recall the case of the Boltzmann operator $C(f)$. We have the following Theorem [32]:

**Theorem 1** (Boltzmann's H-Theorem). Let $f: \mathbb{R}^3 \to \mathbb{R}^+$ be a smooth positive solution of the functional equation such that

$$ \forall (w, w_*, n) \in \mathbb{R}^3 \times \mathbb{R}^3 \times \mathbb{S}^2, \quad f(w'')f(w_*') = f(w)f(w_*), $$

where $w''$ and $w_*'$ are given by (3). Then there exists $\rho_g, u_g, \theta_g > 0$ such that $f = \mathcal{M}[\rho_g, u_g, \theta_g]$, where

$$ \mathcal{M}[\rho_g, u_g, \theta_g](w) := \frac{\rho_g m_g^{3/2}}{(2\pi k_B \theta_g)^{3/2}} e^{-\frac{m_g |w-u_g|^2}{2k_B \theta_g}}. \quad (45) $$

Moreover, relations (4) implies that

$$ |w'' - u_g|^2 + |w_*' - u_g|^2 = |w - u_g|^2 + |w_* - u_g|^2 $$

for any $u_g \in \mathbb{R}^3$, and then the reciprocal is true. It follows that, in case of equilibrium,

$$ D_B(f) = \int_{\mathbb{R}^3} \ln f(w) C(f)(w) dw = 0 \quad \Leftrightarrow \quad f = \mathcal{M}[\rho_g, u_g, \theta_g], $$

with

$$ \begin{aligned} \rho_g(t, x) &:= \int_{\mathbb{R}^3} f(t, x, w) dw, \quad u_g(t, x) := \frac{1}{\rho_g(t, x)} \int_{\mathbb{R}^3} f(t, x, w) w dw, \\ \theta_g(t, x) &:= \frac{m_g}{3k_B \rho_g(t, x)} \int_{\mathbb{R}^3} f(t, x, w) |w - u_g(t, x)|^2 dw. \end{aligned} \quad (46) $$

We also have

$$ \int_{\mathbb{R}^3} \ln F(v) \mathcal{P}(F)(v) dv = 0 \quad \Leftrightarrow \quad F = \mathcal{M}[\rho_p, u_p, \theta_p], $$

with

$$ \begin{aligned} \rho_p(t, x) &:= \int_{\mathbb{R}^3} F(t, x, v) dv, \quad u_p(t, x) := \frac{1}{\rho_p(t, x)} \int_{\mathbb{R}^3} F(t, x, v) v dv, \\ \theta_p(t, x) &:= \frac{m_p}{3k_B \rho_p(t, x)} \int_{\mathbb{R}^3} F(t, x, v) |v - u_p(t, x)|^2 dv. \end{aligned} \quad (47) $$

We now extend this result to the case of operators $\mathcal{D}(F, f)$ and $\mathcal{R}(F, f)$.

**Lemma 4.** Let $F$ and $f$ in $C^1(\mathbb{R}^3) \cap L^1(\mathbb{R}^3)$, such that

$$ \forall (v, w, z) \in \mathbb{R}^3 \times \mathbb{R}^3 \times \mathbb{R}^3, \quad F(v')f(w')e^{c(m_p|v'|^2+m_g|w'|^2)} = F(v)f(w)e^{c(m_p|v|^2+m_g|w|^2)}, \quad (48) $$

where $v'$ and $w'$ are given by (17) and (19), where $c$ is a given constant. Then $F$ and $f$ are Maxwellian states. More precisely, there exist $u \in \mathbb{R}^3$, $N_f > 0$, $N_F > 0$ such that

$$ \begin{cases} F(v) = N_F e^{-c(m_p|v-u|^2)}, \\ f(w) = N_F e^{-c(m_g|w-u|^2)}. \end{cases} \quad (49) $$

Moreover, the converse is also true. That is, if $F$ and $f$ are given by (49) for a given $u \in \mathbb{R}^3$, $c \in \mathbb{R}$, then relation (48) holds.

*Proof.* The proof follows the same idea as in [33], but is adapted to the particular collision mechanism that we consider here. We multiply Equation (48) by $h_n(z)$ and integrate over the parameters $z \in \mathbb{R}^3$ and $n \in \mathbb{S}$. This yields:

$$ \begin{aligned} & \int_{\mathbb{S}^2} \int_{\mathbb{R}^3} F(v) f(w) e^{c(m_p|v|^2+m_g|w|^2)} h_n(z) dz dn \\ &= \int_{\mathbb{S}^2} \int_{\mathbb{R}^3} F(v') f(w') e^{c(m_p|v'|^2+m_g|w'|^2)} h_n(z) dz dn. \end{aligned} \quad (50) $$
Let

$$v_B = \frac{m_p v + m_g w}{m_p + m_g}.$$

Thanks to relations (17) and (19), we have

$$m_p |v'|^2 + m_g |w'|^2 = (m_p + m_g) |v_B|^2 + \frac{m_p m_g}{m_p + m_g} |z|^2.$$

Then, by using the expression (8) for $h_n(z)$ and its property (110) with $c = \sqrt{\frac{m_g}{k_b T_p(t,x)}}$ (Appendix A), we arrive at the rewritten equality:

$$\begin{aligned} & 4\pi F(v)f(w)e^{c(m_p|v|^2+m_g|w|^2)} \\ &= \int_{\mathbb{R}^3} \int_{S^2} F\left(v_B - \frac{m_g}{m_p+m_g}z\right) f\left(v_B + \frac{m_p}{m_p+m_g}z\right) e^{c\left((m_p+m_g)|v_B|^2+\frac{m_p m_g}{m_p+m_g}|z|^2\right)} h_n(z) \mathrm{d}n \mathrm{d}z \\ &= \frac{1}{2} \frac{m_g^2}{k_B^2 T_p^2} \int_{\mathbb{R}^3} F\left(v_B - \frac{m_g}{m_p+m_g}z\right) f\left(v_B + \frac{m_p}{m_p+m_g}z\right) e^{c\left((m_p+m_g)|v_B|^2+\frac{m_p m_g}{m_p+m_g}|z|^2\right)} e^{-\frac{m_g}{2k_B T_p} |z|^2} |z| \mathrm{d}z \mathrm{d}z, \end{aligned}$$

where the last relation arises from the fact that

$$\int_{S^2} [n \cdot z]_+ \mathrm{d}n = \pi |z|.$$

It follows that

$$\Psi(v_B) := F(v)f(w)e^{c(m_p|v|^2+m_g|w|^2)}$$

is a (smooth) function of the vector $v_B$ only. Then, if we denote $E = m_p|v|^2 + m_g|w|^2$, we have

$$\ln F(v) + \ln f(w) = \ln \Psi(v_B) - cE.$$

This allows to obtain that

$$\begin{aligned} \nabla_v \ln F(v) &= \frac{m_p}{m_g+m_p} \nabla(\ln \Psi) \left( \frac{m_p v + m_g w}{m_p+m_g} \right) - 2c m_p v, \\ \nabla_w \ln f(w) &= \frac{m_g}{m_g+m_p} \nabla(\ln \Psi) \left( \frac{m_p v + m_g w}{m_p+m_g} \right) - 2c m_g w, \end{aligned} \quad (51)$$

and then that $m_g \nabla_v \ln F(v) - m_p \nabla_w \ln f(w)$ is colinear to $v-w$:

$$m_g \nabla_v \ln F(v) - m_p \nabla_w \ln f(w) = -2c m_p m_g (v-w), \quad \forall (v,w) \in \mathbb{R}^3 \times \mathbb{R}^3. \quad (52)$$

We can hence conclude that there exists $u \in \mathbb{R}^3$ such that

$$\left\{ \begin{array}{l} \nabla_v \ln F(v) = -2c m_p (v-u), \\ \nabla_w \ln f(w) = -2c m_g (w-u), \end{array} \right. \quad (53)$$

and thus that there exists also $\lambda_1 \in \mathbb{R}$, $\lambda_2 \in \mathbb{R}$ such that

$$\left\{ \begin{array}{l} F(v) = e^{-c(m_p|v-u|^2+\lambda_1)}, \\ f(w) = e^{-c(m_g|w-u|^2+\lambda_2)}, \end{array} \right. \quad (54)$$

which proves (49). Finally, if (49) is satisfied, the conservation of momentum $m_p v' + m_g w' = m_p v + m_g w$ implies

$$-m_p|v' - u|^2 - m_g|w' - u|^2 + m_p|v'|^2 + m_g|w'|^2 = -m_p|v - u|^2 - m_g|w - u|^2 + m_p|v|^2 + m_g|w|^2$$
for any $c \in \mathbb{R}^3$, and (48) holds. $\square$

This leads to the following (formal) characterization of the equilibrium states of System (1).
**Proposition 6.** Let $f$ and $F$ smooth functions that verify (1), with $\zeta$ given by (13). Then,

$$D(F, f)(t, x) = 0, \quad \text{for any } (t, x) \in \mathbb{R}^{+} \times \mathbb{R}^{3},$$

if and only if there exists $u^{eq} = u^{eq}(t, x) \in \mathbb{R}^{3}$, $N_{f} = N_{f}(t, x) \ge 0$, and $N_{F} = N_{F}(t, x) \ge 0$ such that

$$\left\{ \begin{array}{l} f(t, x, w) = N_{f}(t, x) \exp \left( -\frac{\mu_{gp}}{2k_{B}T_{p}(t, x)} m_{g} |w - u^{eq}(t, x)|^{2} \right) =: f^{eq}(t, x, w), \\ F(t, x, v) = N_{F}(t, x) \exp \left( -\frac{\mu_{gp}}{2k_{B}T_{p}(t, x)} m_{p} |v - u^{eq}(t, x)|^{2} \right) =: F^{eq}(t, x, v), \end{array} \right.$$

where

$$\mu_{gp} = \frac{m_{p} + m_{g}}{m_{p}}.$$

In particular, we have

$$\theta_{g}(t, x) = \theta_{p}(t, x) = \frac{T_{p}(t, x)}{\mu_{gp}},$$

where $\theta_g$ and $\theta_p$ are defined by (46) and (47), and the macroscopic velocities of $f$ and $F$ are equal:

$$u_{g}(t, x) = u_{p}(t, x) = u^{eq}(t, x).$$

*Proof.* Thanks to Lemma 3, we have

$$\begin{array}{l} D(F, f)(t, x) = 0 \\ \Leftrightarrow D^e(F, f)(t, x) = 0 \quad \text{and} \quad D^d(F, f)(t, x) = 0 \\ \Leftrightarrow \text{for any } w \in \mathbb{R}^3, n \in \mathbb{S}^2, v \in \mathbb{R}^3, z \in \mathbb{R}^3 \\ \left\{ \begin{array}{l} f(w'')f(w''_{*}) = f(w)f(w_{*}), \\ F(v'')F(v''_{*}) = F(v)F(v_{*}), \\ F(v')f(w')e^{\frac{\mu_{gp}}{2k_{B}T_{p}}(m_{p}|v'|^{2}+m_{g}|w'|^{2})} = F(v)f(w)e^{\frac{\mu_{gp}}{2k_{B}T_{p}}(m_{p}|v'|^{2}+m_{g}|w'|^{2})}. \end{array} \right. $$

From Theorem 1, we deduce that the two first relations are satisfied if and only if $f$ and $F$ are Maxwellian states, and, from Lemma 4, that the last line gives an additional condition on the macroscopic velocities and temperatures. $\square$

**Remark 5.** We underline the importance of treating the temperature $T_p(t, x)$ as a space-time dependent unknown rather than as a microscopic independent variable, as done in [16]. This choice enables us to derive an expression for the system’s entropy (41), which plays a crucial role in determining the system’s equilibria and in the analytical study of the model itself. The latter will be explored in a future work.

# 4 Adimensionalization

Macroscopic quantities such as density, momentum, and energy, for both gases and dust, can be obtained by integrating the distribution functions $f$ and $F$ with respect to velocities $v$ and $w$, respectively. These quantities are called 'moments', and the equations describing their evolution in space and time can be derived by directly integrating the first equation of System (1) with respect to $w$, after multiplying by $1, w, |w|^2$, and the second equation of System (1) with respect to $v$, after multiplying by $1, v, |v|^2$.

However, the resulting system involves integrals of $f$ and $F$, which cannot be expressed in terms of the macroscopic quantities of mass, momentum, and energy, unless $f$ and $F$ have a particular explicit expression. Thus, the resulting system is generally not closed. To resolve this, it is necessary to specify a form of the distribution functions in terms of the independent variables. This cannot be done in a fully general manner; instead, it requires defining a specific regime in which the approximate form of the distribution functions remains valid. Consequently, the equations must be nondimensionalized using an appropriate scaling.
## 4.1 Scaling units and main hypothesis

We present here the adopted scaling and the main hypotheses we made to perform the diffusion asymptotics.

**Scaling of variables** In Table 1, we introduce the scaling units. We then make the following

<table><caption>Table 1: Scaling units</caption><thead><tr><th>Scaling unit</th><th>Meaning</th></tr></thead><tbody><tr><td><i>L</i></td><td>macroscopic length scale</td></tr><tr><td><i>t</i><sup>o</sup></td><td>macroscopic time scale</td></tr><tr><td><i>V</i><sup>o</sup></td><td>reference velocity of particles</td></tr><tr><td><i>W</i><sup>o</sup></td><td>reference velocity of molecules</td></tr><tr><td><i>V</i><sub>rel</sub></td><td>scaling of the relative velocity</td></tr><tr><td><i>N</i><sub>g</sub></td><td>order of magnitude of the number density in molecules</td></tr><tr><td><i>N</i><sub>p</sub></td><td>order of magnitude of the number density in particles</td></tr><tr><td><i>T</i><sup>o</sup></td><td>reference temperature</td></tr></tbody></table>

change of space, time, and velocities coordinates:

$$ \hat{x} = \frac{x}{L}, \quad \hat{t} = \frac{t}{t^{\circ}}, \quad \hat{v} = \frac{v}{V^{\circ}}, \quad \hat{w} = \frac{w}{W^{\circ}}, \quad \hat{z} = \frac{z}{V_{\text{rel}}}, $$

and we define the dimensionless unknowns

$$ \hat{F}(\hat{t}, \hat{x}, \hat{v}) = \frac{F(t, x, v)}{F^{\circ}}, \quad \text{with } F^{\circ} = \frac{N_p}{(V^{\circ})^3}, $$

$$ \hat{f}(\hat{t}, \hat{x}, \hat{w}) = \frac{f(t, x, w)}{f^{\circ}}, \quad \text{with } f^{\circ} = \frac{N_g}{(W^{\circ})^3}, $$

$$ \hat{T}_p(\hat{t}, \hat{x}) = \frac{T_p(t, x)}{T^{\circ}}. $$

**Small parameters** We introduce the following small parameters:

* $\eta$, the mass ratio between a molecule and a particle,

$$ \eta := \frac{m_g}{m_p}, \qquad (55) $$

* $\nu$, the ratio of the molecular radius of the gas to the particle radius of the dust,

$$ \nu := \frac{r_g}{r_p}, \qquad (56) $$

* $\delta$, the Knudsen number related to the gas

$$ \delta := \frac{1}{N_g \sigma_{gg} L}, $$

* $\gamma$, the Mach number of the gas

$$ \gamma := \frac{L}{W^{\circ} t^{\circ}}. $$
**Remark 6.** Note that the Knudsen number of the dust, $\delta_p$, depends on the Knudsen number of the gas and the ratio of their respective radii:

$$ \delta_p := \frac{1}{N_p \sigma_{pp} L} \sim \delta \nu^2. \qquad (57) $$

Similarly, the Mach number of the dust, $\gamma_p$, depends on the Mach number of the gas and on the ratio of their respective masses:

$$ \gamma_p := \frac{L}{V^\circ t^\circ} = \frac{\gamma}{\sqrt{\eta}}. \qquad (58) $$

Moreover, we have:

$$ \frac{\sigma_{gg}}{\sigma_{pp}} \sim \nu^2, \quad \frac{\sigma_{gg}}{\sigma_{gp}} \sim \frac{\nu^2}{(1 + \nu)^2} \sim \nu^2. $$

In particular, the previous relationships imply that $\sigma_{gp}/\sigma_{pp} = O(1)$.

## Hypotheses on the orders of magnitude

**Hypothesis 1.** We assume that

$$ \eta \ll 1, \quad \nu \ll 1, \quad \delta \ll 1, \quad \gamma \ll 1. $$

**Hypothesis 2** (Order of magnitude of the kinetic temperatures). We assume that the velocity scales $V^\circ$ and $W^\circ$ of dust and gas are given by their thermal velocities, and that the order of magnitude of the kinetic temperatures of gas and dust is $T^\circ$:

$$ V^\circ = \sqrt{\frac{k_B T^\circ}{m_p}}, \quad W^\circ = \sqrt{\frac{k_B T^\circ}{m_g}}. $$

Consequently, these velocities depend only on the masses, and the relation

$$ \frac{V^\circ}{W^\circ} \sim \sqrt{\eta} $$

holds, where $\eta$ is the mass ratio (55).

**Remark 7.** Hypothesis 2 explains the introduction of two different scales for the dust and gas velocities, since - from the latter - it follows that

$$ V^\circ = \sqrt{\eta} W^\circ \ll W^\circ. $$

In particular, it follows that

$$ V_{\text{rel}} \sim W^\circ. $$

**Hypothesis 3** (Order of magnitude of the number densities). We assume that

$$ \frac{N_p}{N_g} \sim 1. $$

This means that the ratio of the orders of magnitude of the particle number density and the gas molecule number density remains constant, regardless of the parameters $\eta, \nu, \delta, \gamma$.
## 4.2 Adimensionalized Boltzmann system

By considering the physical dimension [9], the collision kernels $c(w - w_*, n)$, $p(v - v_*, n)$ and $\zeta(v - w, n)$ are scaled as follows:

$$ 
\begin{aligned}
\zeta(v - w, n) &= W^{\circ}\zeta(\sqrt{\eta}\hat{v} - \hat{w}, n), \\
c(w - w_*, n) &= W^{\circ}\hat{c}(\hat{w} - \hat{w}_*), \\
p(v - v_*, n) &= V^{\circ}\hat{p}(\hat{v} - \hat{v}_*).
\end{aligned}
 $$

The dimensionless operators are then derived using the following formulas:

$$ 
\begin{aligned}
\mathcal{D}(F, f) &= \sigma_{gp} F^{\circ} f^{\circ} (W^{\circ})^4 \hat{\mathcal{D}}(\hat{F}, \hat{f}), \\
\mathcal{P}(F) &= \sigma_{pp} (F^{\circ})^2 (V^{\circ})^4 \hat{\mathcal{P}}(\hat{F}), \\
\mathcal{R}(F, f) &= \sigma_{gp} F^{\circ} f^{\circ} W^{\circ} (V^{\circ})^3 \hat{\mathcal{R}}(\hat{F}, \hat{f}), \\
\mathcal{C}(f) &= \sigma_{gg} (f^{\circ})^2 (W^{\circ})^4 \hat{\mathcal{C}}(\hat{f}),
\end{aligned}
 $$

where $\hat{\mathcal{D}}(\hat{F}, \hat{f})$, $\hat{\mathcal{P}}(\hat{F})$, $\hat{\mathcal{R}}(\hat{F}, \hat{f})$, and $\hat{\mathcal{C}}(\hat{f})$ are the dimensionless version of $\mathcal{C}(f)$, $\mathcal{P}(F)$, $\mathcal{D}(F, f)$ and $\mathcal{R}(F, f)$, given by:

$$ \hat{\mathcal{C}}(\hat{f})(\hat{w}) = \int_{\mathbb{R}^3} \int_{S^2} [\hat{f}(\hat{w}')\hat{f}(\hat{w}_*)' - \hat{f}(\hat{w})\hat{f}(\hat{w}_*)] \hat{c}(\hat{w} - \hat{w}_*, n) \, dnd\hat{w}_*, \quad (59) $$

$$ \hat{\mathcal{P}}(\hat{F})(\hat{v}) = \int_{\mathbb{R}^3} \int_{S^2} [\hat{F}(\hat{v}')\hat{F}(\hat{v}_*)' - \hat{F}(\hat{v})\hat{F}(\hat{v}_*)] \hat{p}(\hat{v} - \hat{v}_*, n) \, dnd\hat{v}_*, \quad (60) $$

$$ \begin{aligned} \hat{\mathcal{D}}(\hat{F}, \hat{f})(\hat{t}, \hat{x}, \hat{v}) = & \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} \int_{S^2} [\hat{f}(\hat{w}') \hat{F}(\hat{v}') h_n(\hat{t}, \hat{x}, \sqrt{\eta}\hat{v} - \hat{w}) \zeta(\hat{z}, n) - \hat{f}(\hat{w}) \hat{F}(\hat{v}) \hat{h}_n(\hat{t}, \hat{x}, \hat{z}) \zeta(\sqrt{\eta}\hat{v} - \hat{w}, n)] \, d\hat{z} dnd\hat{w}, \end{aligned} \quad (61) $$

and

$$ \begin{aligned} \hat{\mathcal{R}}(\hat{F}, \hat{f})(\hat{t}, \hat{x}, \hat{w}) = & \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} \int_{S^2} [\hat{f}(\hat{w}') \hat{F}(\hat{v}') h_n(\hat{t}, \hat{x}, \sqrt{\eta}\hat{v} - \hat{w}) \zeta(\hat{z}, n) - \hat{f}(\hat{w}) \hat{F}(\hat{v}) \hat{h}_n(\hat{t}, \hat{x}, \hat{z}) \zeta(\sqrt{\eta}\hat{v} - \hat{w}, n)] \, d\hat{z} dnd\hat{v}, \end{aligned} \quad (62) $$

with

$$ \left\{ \begin{aligned} \hat{v}' &= \frac{\hat{v} + \sqrt{\eta \hat{w}}}{1 + \eta} - \frac{\sqrt{\eta}}{1 + \eta} \hat{z}, \\ \hat{w}' &= \frac{\sqrt{\eta \hat{v}} + \eta \hat{w}}{1 + \eta} + \frac{1}{1 + \eta} \hat{z}, \end{aligned} \right. \quad (63) $$

and

$$ \hat{h}_n(\hat{t}, \hat{x}, \hat{z}) = \frac{1}{2\pi \hat{T}_p(\hat{t}, \hat{x})^2} (n \cdot \hat{z}) e^{-\frac{|\hat{z}|^2}{2\hat{T}_p^2}} 1_{\{\hat{z} \cdot n \ge 0\}}. $$

For the collision operators $\hat{\mathcal{D}}$ and $\hat{\mathcal{R}}$, we have made the dependence on $\hat{t}$ and $\hat{x}$ explicit, although we have expressed it only in $\hat{T}_p$ and $\hat{h}_n$, and not in $\hat{F}$ and $\hat{f}$. This choice is just to remind the reader of the dependence of the surface temperature on the space-time variables themselves. System (1) is then written as:

$$ \left\{ \begin{aligned} \partial_{\hat{t}} \hat{F} + \frac{V^{\circ} t^{\circ}}{L} \hat{v} \cdot \nabla_{\hat{x}} \hat{F} &= N_g \sigma_{gp} W^{\circ} t^{\circ} \hat{\mathcal{D}}(\hat{F}, \hat{f}) + N_p \sigma_{pp} V^{\circ} t^{\circ} \hat{\mathcal{P}}(\hat{F}), \\ \partial_{\hat{t}} \hat{f} + \frac{W^{\circ} t^{\circ}}{L} \hat{w} \cdot \nabla_{\hat{x}} \hat{f} &= N_p \sigma_{gp} W^{\circ} t^{\circ} \hat{\mathcal{R}}(\hat{F}, \hat{f}) + N_g \sigma_{gg} W^{\circ} t^{\circ} \hat{\mathcal{C}}(\hat{f}). \end{aligned} \right. $$
## 4.3 Adimensionalized surface temperature's equation

Expressions $\Upsilon_1(F, f)$ and $\Upsilon_2(F, f)$ defined by (34) and (35) are scaled in the following way

$$ \Upsilon_1(F, f) = N_p N_g \sigma_{gp} W^\circ \hat{\Upsilon}_1(\hat{F}, \hat{f}), \quad \Upsilon_2(F, f) = N_p N_g \sigma_{gp} (W^\circ)^3 \hat{\Upsilon}_2(\hat{F}, \hat{f}), $$

with

$$ \left\{ \begin{aligned} \hat{\Upsilon}_1(\hat{F}, \hat{f}) &= \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} \hat{F}(\hat{v}) \hat{f}(\hat{w}) \hat{\Sigma}(\sqrt{\eta}\hat{v} - \hat{w}) d\hat{v} d\hat{w}, \\ \hat{\Upsilon}_2(\hat{F}, \hat{f}) &= \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} \hat{F}(\hat{v}) \hat{f}(\hat{w}) \hat{\Sigma}(\sqrt{\eta}\hat{v} - \hat{w}) |\sqrt{\eta}\hat{v} - \hat{w}|^2 d\hat{v} d\hat{w}. \end{aligned} \right. $$

Therefore, the expression for (36) showed the scaled quantities is

$$ \begin{aligned} & 2c_p T^\circ N_p \left[ \frac{1}{t^\circ} \partial_t \hat{T}_p \int_{\mathbb{R}^3} \hat{F}(\hat{v}) d\hat{v} + \frac{V^\circ}{L} \nabla_{\hat{x}} \hat{T}_p \cdot \int_{\mathbb{R}^3} \hat{F}(\hat{v}) \hat{v} d\hat{v} \right] + \frac{4k_B}{(m_g + m_p)} N_p N_g \sigma_{gp} W^\circ T^\circ \hat{T}_p \hat{\Upsilon}_1(\hat{F}, \hat{f}) \\ &= \frac{m_g}{(m_g + m_p)} N_p N_g \sigma_{gp} (W^\circ)^3 \hat{\Upsilon}_2(\hat{F}, \hat{f}). \end{aligned} $$

## 4.4 Summary of the dimensionless equations

We now remove all the hats from the dimensionless formulas and rewrite the scaling units in terms of the small parameters $\eta$, $\mu$, $\delta$, and $\gamma$. We thus end up with the following coupled integro-differential system:

$$ \left\{ \begin{aligned} \frac{\gamma}{\sqrt{\eta}} \partial_t F + v \cdot \nabla_x F &= \frac{1}{\delta \nu^2 \sqrt{\eta}} \mathcal{D}(F, f) + \frac{1}{\delta \nu^2} \mathcal{P}(F), \\ \gamma \partial_t f + w \cdot \nabla_x f &= \frac{1}{\delta \nu^2} \mathcal{R}(F, f) + \frac{1}{\delta} \mathcal{C}(f), \end{aligned} \right. \qquad (64) $$

where $\mathcal{D}$, $\mathcal{P}$, $\mathcal{R}$, and $\mathcal{C}$ are the dimensionless operators (61), (60), (62), and (59), and the non-linear transport equation for $T_p$:

$$ \begin{aligned} & \tau \left( \gamma \partial_t T_p(t, x) \int F(t, x) dv + \sqrt{\eta} \nabla_x T_p(t, x) \cdot \int v F(t, x) dv \right) + \frac{4}{\delta \nu^2} \frac{1}{1 + \eta} \Upsilon_1(F, f)(t, x) T_p(t, x) \\ &= \frac{1}{\delta \nu^2} \frac{1}{1 + \eta} \Upsilon_2(F, f)(t, x), \end{aligned} \qquad (65) $$

with

$$ \Upsilon_1(F, f)(t, x) = \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} F(v) f(w) \Sigma(\eta v - w) dv dw, \qquad (66) $$

and

$$ \Upsilon_2(F, f)(t, x) = \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} F(v) f(w) \Sigma(\eta v - w) |\eta v - w|^2 dv dw. \qquad (67) $$

The constant

$$ \tau := \frac{2m_p c_p}{k_B}, \qquad (68) $$

thanks to the Hypothesis 2, this can be interpreted as the ratio between the order of magnitude of the internal energy of the dust ($m_p c_p T^\circ$) and that of the kinetic temperature of the gas, i.e. $\frac{1}{2} m_g (W^\circ)^2$. It depends only on the gas and the material of which the particles are made (but not on the size or number of particles) and has no unit of measurement. Moreover, thanks to (63), the relation (39) become after scaling

$$ \sqrt{\eta} \int \mathcal{R}(F, f)(w) w dv + \int \mathcal{D}(F, f)(v) v dv = 0. \qquad (69) $$
## 4.5 Dimensionless equilibrium states and entropy dissipation

With this scaling, and considering that $\hat{u}^{eq} = W^{\circ}u^{eq}$, we obtain the following forms for entropy dissipation and equilibrium, whose expressions can be derived directly from Lemma 3:

**Lemma 5.** The re-scaled equilibrium states have the following form:

$$ egin{cases} f^{\text{eq}}(t, x, w) = N_f(t, x) \exp \left( -\frac{(1+\eta)}{2T_p(t, x)} |w - u^{\text{eq}}(t, x)|^2 \right), \\ F^{\text{eq}}(t, x, v) = N_F(t, x) \exp \left( -\frac{(1+\eta)}{2T_p(t, x)} |v - \frac{1}{\sqrt{\eta}} u^{\text{eq}}(t, x)|^2 \right). \end{cases} $$

# 5 Diffusive scaling and balance laws

This section is devoted to the study of the diffusion limit of the system (64)-(65), i.e. when the small parameters $\delta$ and $\gamma$, representing the Knudsen number and the Mach number of the gas, respectively, tend to 0. In addition, we also simultaneously consider the asymptotic of the mass and size ratios between the two species to 0, i.e. $\eta \to 0$ and $\nu \to 0$.

After having introduced $\epsilon = (\gamma, \nu, \delta, \eta)$, we consider $(f^\epsilon, F^\epsilon, T_p^\epsilon)$ a smooth solution of (64)-(65). The macroscopic densities of $f^\epsilon$ and $F^\epsilon$ are defined by:

$$ \rho_F^\epsilon(t, x) = \int F^\epsilon(t, x, v)dv, \quad \rho_f^\epsilon(t, x) = \int f^\epsilon(t, x, w)dw. \quad (70) $$

Since we are interested in diffusive behavior, we make the following additional hypothesis, as in [8].

**Hypothesis 4.** The bulk velocities of the gas and particles are small and approach zero as the respective Mach numbers, $\gamma$ and $\gamma_p$, approach zero. In the limit of vanishing Mach number, each bulk velocity is at least of order one.

Hypothesis 4 leads to the definition of $u_F^\epsilon$ and $u_f^\epsilon$, which is consistent with the different speed scaling of the two species:

$$ u_F^\epsilon(t, x) = \frac{\sqrt{\eta}}{\gamma \rho_F^\epsilon(t, x)} \int F^\epsilon(t, x, v) v dv, \quad u_f^\epsilon(t, x) = \frac{1}{\gamma \rho_f^\epsilon(t, x)} \int f^\epsilon(t, x, w) w dw. \quad (71) $$

**Remark 8.** Note that $u_F^\epsilon$ and $u_f^\epsilon$ do not correspond to the macroscopic velocities as they are usually defined, but rather to the first-order term in $\gamma/\sqrt{\eta}$ (respectively $\gamma$) of the asymptotic expansion of these macroscopic velocities.

We also define the kinetic temperature of the species:

$$ \theta_F^\epsilon(t, x) = \frac{1}{3\rho_F^\epsilon} \int F^\epsilon(v) \left| v - \frac{\gamma}{\sqrt{\eta}} u_F^\epsilon \right|^2 dv, \quad \theta_f^\epsilon(t, x) = \frac{1}{3\rho_f^\epsilon} \int f^\epsilon(w) |w - \gamma u_f^\epsilon|^2 dw. \quad (72) $$

According to definitions (70), (71), and (72), the moments of order 0, 1 and 2 for $F^\epsilon$ are as follows:

$$ \int_{\mathbb{R}^3} F^\epsilon(t, x, v) \begin{pmatrix} 1 \\ v \\ |v|^2 \end{pmatrix} dv = \begin{pmatrix} \rho_F^\epsilon \\ \frac{\gamma}{\sqrt{\eta}} \rho_F^\epsilon u_F^\epsilon \\ 3\rho_F^\epsilon \theta_F^\epsilon + \frac{\gamma^2}{\eta} |u_F^\epsilon|^2 \rho_F^\epsilon \end{pmatrix}, \quad (73) $$

and the moments of order 0, 1 and 2 for $f^\epsilon$ are:

$$ \int_{\mathbb{R}^3} f^\epsilon(t, x, w) \begin{pmatrix} 1 \\ w \\ |w|^2 \end{pmatrix} dw = \begin{pmatrix} \rho_f^\epsilon \\ \gamma \rho_f^\epsilon u_f^\epsilon \\ 3 \rho_f^\epsilon \theta_f^\epsilon + \gamma^2 |u_f^\epsilon|^2 \rho_f^\epsilon \end{pmatrix}. \quad (74) $$
## 5.1 Ansatz

The main hypothesis of our work is that the initial conditions ($f^\epsilon(0, \cdot, \cdot)$, $F^\epsilon(0, \cdot, \cdot)$) are local Maxwellian states and that the evolution governed by (64) and (65) preserves the local Maxwellian nature of the distribution functions, with an error which remains small with respect to all the small parameters involved in the computations. Physically, this assumption is equivalent to stating that the mixtures have reached their local mechanical equilibrium. This approximation is justified in the limit as the vector parameter $\epsilon$ approaches the zero vector, as we will see in the next Section. We then have

$$F^{\epsilon}(t, x, v) = \rho_{F}^{\epsilon}(t, x) \left( \frac{1}{2\pi\theta_{F}^{\epsilon}(t, x)} \right)^{\frac{3}{2}} \exp \left( -\frac{1}{2\theta_{F}^{\epsilon}(t, x)} |v - \frac{\gamma}{\sqrt{\eta}} u_{F}^{\epsilon}(t, x)|^2 \right), \quad (75)$$

$$f^{\epsilon}(t, x, w) = \rho_{f}^{\epsilon}(t, x) \left( \frac{1}{2\pi\theta_{f}^{\epsilon}(t, x)} \right)^{\frac{3}{2}} \exp \left( -\frac{1}{2\theta_{f}^{\epsilon}(t, x)} |w - \gamma u_{f}^{\epsilon}(t, x)|^2 \right), \quad (76)$$

where $\rho_F^\epsilon, \rho_f^\epsilon, u_F^\epsilon, u_f^\epsilon, \theta_F^\epsilon, \theta_f^\epsilon$ are defined by (70), (71), (72). The choice of ansatz (75) and (76) is the usual line of attack to treat diffusion limit methods in the context of Boltzmann models, see for example [1, 3, 8, 24, 25].

## 5.2 Balance laws

To obtain a macroscopic description of the gas-dust aerosol mixture, we therefore start with the kinetic description (64), combined with the equation (65) describing the surface temperature of the dust, and integrate them only with respect to the velocity variables $v$ and $w$.

### 5.2.1 Conservation of mass

**Proposition 7.** Let ($f^\epsilon, F^\epsilon$) be a solution of system (64), and let $\rho_F^\epsilon, \rho_f^\epsilon, u_F^\epsilon, u_f^\epsilon$ be defined by (70) and (71). We have the mass balance equations:

$$\left\{ \begin{array}{l} \partial_t \rho_F^\epsilon + \nabla_x \cdot (\rho_F^\epsilon u_F^\epsilon) = 0, \\ \\ \partial_t \rho_f^\epsilon + \nabla_x \cdot (\rho_f^\epsilon u_f^\epsilon) = 0. \end{array} \right. \qquad (77)$$

*Proof.* This result follows directly from integrating the first and second equations of (64) with respect to $v$ and $w$, respectively, since 1 is a collision invariant for all $D, P, R$ and $C$, and we can explicitly compute the momentum of order 0 and 1 of $F^\epsilon$ and $f^\epsilon$. $\square$

### 5.2.2 Balance of momentum

**Proposition 8.** Under the same assumptions as in Proposition 7, we obtain the momentum balance:

$$\left\{ \begin{array}{l} \frac{\gamma^2}{\eta} [\partial_t (\rho_F^\epsilon u_F^\epsilon) + \nabla_x \cdot (u_F^\epsilon \otimes u_F^\epsilon \rho_F^\epsilon)] + \nabla_x (\rho_F^\epsilon \theta_F^\epsilon) = \frac{1}{\delta\nu^2\sqrt{\eta}} \Theta_F^\epsilon, \\ \\ \gamma^2 [\partial_t (\rho_f^\epsilon u_f^\epsilon) + \nabla_x \cdot (u_f^\epsilon \otimes u_f^\epsilon \rho_f^\epsilon)] + \nabla_x (\rho_f^\epsilon \theta_f^\epsilon) = \frac{1}{\delta\nu^2} \Theta_f^\epsilon, \end{array} \right. \qquad (78)$$

where

$$\Theta_F^\epsilon(t, x) := \int_{\mathbb{R}^3} vD(F^\epsilon, f^\epsilon)(t, x, v) \, dv, \quad \Theta_f^\epsilon(t, x) := \int_{\mathbb{R}^3} wR(F^\epsilon, f^\epsilon)(t, x, w) \, dw.$$

In addition, for a general collision kernel, the operators $\Theta_F^\epsilon$ and $\Theta_f^\epsilon$ can be explicitly written as

$$\Theta_F^\epsilon = -\frac{\sqrt{\eta}}{1+\eta} \left[ \begin{aligned} & \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} (\sqrt{\eta}v - w) F^\epsilon(v) f^\epsilon(w) \Sigma(\sqrt{\eta}v - w) dw dv \\ & + \sqrt{\frac{\pi}{2}} \sqrt{T_p} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} F^\epsilon(v) f^\epsilon(w) \int_{S^2} n\zeta(\sqrt{\eta}v - w, n) dndw dv \end{aligned} \right], \quad (79)$$
and

$$ \Theta_f^\epsilon = \frac{1}{1+\eta} \left[ \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} (\sqrt{\eta}v-w) F^\epsilon(v) f^\epsilon(w) \Sigma(\sqrt{\eta}v-w) \mathrm{d}w \mathrm{d}v \right. \\ \left. \qquad + \sqrt{\frac{\pi}{2}} \sqrt{T_p} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} F^\epsilon(v) f^\epsilon(w) \int_{S^2} n_\zeta(\sqrt{\eta}v-w, n) \mathrm{d}n \mathrm{d}w \mathrm{d}v \right]. \quad (80) $$

*Proof.* See Appendix B.1. □

### 5.2.3 Balance of energy

**Proposition 9.** Under the same assumptions as in Proposition 7, the energy balance also holds:

$$ \left\{ \begin{array}{l} \partial_t (3\rho_F^\epsilon \theta_F^\epsilon) + \nabla_x \cdot (5u_F^\epsilon \theta_F^\epsilon \rho_F^\epsilon) + \frac{\gamma^2}{\eta} (\partial_t (|u_F^\epsilon|^2 \rho_F^\epsilon) + \nabla_x \cdot (|u_F^\epsilon|^2 u_F^\epsilon \rho_F^\epsilon)) = \frac{1}{\gamma \delta \nu^2} \Xi_F^\epsilon, \\ \\ \partial_t (3\rho_f^\epsilon \theta_f^\epsilon) + \nabla_x \cdot (5u_f^\epsilon \theta_f^\epsilon \rho_f^\epsilon) + \gamma^2 (\partial_t (|u_f^\epsilon|^2 \rho_f^\epsilon) + \nabla_x \cdot (|u_f^\epsilon|^2 u_f^\epsilon \rho_f^\epsilon)) = \frac{1}{\gamma \delta \nu^2} \Xi_f^\epsilon, \end{array} \right. \quad (81) $$

where

$$ \Xi_F^\epsilon := \int_{\mathbb{R}^3} |v|^2 \mathcal{D}(F^\epsilon, f^\epsilon) \, \mathrm{d}v, \quad \Xi_f^\epsilon := \int_{\mathbb{R}^3} |w|^2 \mathcal{R}(F^\epsilon, f^\epsilon) \, \mathrm{d}w. \quad (82) $$

Additionally, for a general collision kernel, the operators $\Xi_F^\epsilon$ and $\Xi_f^\epsilon$ can be explicitly expressed as

$$ \Xi_F^\epsilon = \frac{1}{(1+\eta)^2} \left[ \begin{aligned} & -(\eta^2 + 2\eta) \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} |v|^2 F^\epsilon(v) f^\epsilon(w) \Sigma(\sqrt{\eta}v - w) \mathrm{d}w \mathrm{d}v \\ & + \eta \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} |w|^2 F^\epsilon(v) f^\epsilon(w) \Sigma(\sqrt{\eta}v - w) \mathrm{d}w \mathrm{d}v \\ & + 4T_p^\epsilon \eta \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} F^\epsilon(v) f^\epsilon(w) \Sigma(\sqrt{\eta}v - w) \mathrm{d}w \mathrm{d}v \\ & + 2\sqrt{\eta} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} v \cdot w F^\epsilon(v) f^\epsilon(w) \Sigma(\sqrt{\eta}v - w) \mathrm{d}w \mathrm{d}v \\ & - \sqrt{2\pi T_p^\epsilon} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} (\sqrt{\eta}v + \eta w) F^\epsilon(v) f^\epsilon(w) \cdot \int_{S^2} n_\zeta(\sqrt{\eta}v - w, n) \mathrm{d}n \mathrm{d}w \mathrm{d}v \end{aligned} \right]. \quad (83) $$

and

$$ \Xi_f^\epsilon = \frac{1}{(1+\eta)^2} \left[ \begin{aligned} & \eta \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} |v|^2 F^\epsilon(v) f^\epsilon(w) \Sigma(\sqrt{\eta}v - w) \mathrm{d}w \mathrm{d}v \\ & - (1+2\eta) \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} |w|^2 F^\epsilon(v) f^\epsilon(w) \Sigma(\sqrt{\eta}v - w) \mathrm{d}w \mathrm{d}v \\ & + 4T_p^\epsilon \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} F^\epsilon(v) f^\epsilon(w) \Sigma(\sqrt{\eta}v - w) \mathrm{d}w \mathrm{d}v \\ & + 2\eta\sqrt{\eta} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} v \cdot w F^\epsilon(v) f^\epsilon(w) \Sigma(\sqrt{\eta}v - w) \mathrm{d}w \mathrm{d}v \\ & + \sqrt{2\pi T_p^\epsilon} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} (\sqrt{\eta}v + \eta w) F^\epsilon(v) f^\epsilon(w) \cdot \int_{S^2} n_\zeta(\sqrt{\eta}v - w, n) \mathrm{d}n \mathrm{d}w \mathrm{d}v \end{aligned} \right]. \quad (84) $$

Moreover, the following relation is satisfied

$$ \Xi_F^\epsilon + \Xi_f^\epsilon = \frac{1}{1+\eta} (4T_p^\epsilon \Upsilon_1(F^\epsilon, f^\epsilon) - \Upsilon_2(F^\epsilon, f^\epsilon)), \quad (85) $$
where $\Upsilon_1(F^\epsilon, f^\epsilon)$ and $\Upsilon_2(F^\epsilon, f^\epsilon)$ are given by

$$ egin{cases} 
\Upsilon_1(F^\epsilon, f^\epsilon) = \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} F^\epsilon(v) f^\epsilon(w) \Sigma(\sqrt{\eta}v - w) \mathrm{dvdw}, \\ 
\Upsilon_2(F^\epsilon, f^\epsilon) = \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} F^\epsilon(v) f^\epsilon(w) \Sigma(\sqrt{\eta}v - w) |\sqrt{\eta}v - w|^2 \mathrm{dvdw}. 
\end{cases} 
\quad (86) $$

*Proof.* See Appendix B.2. □

### 5.2.4 Surface temperature's equation

The following Corollary is a direct consequence of the relation (65) and the definition (71) of $u_F^\epsilon$.

**Corollary 4.** If the functions $F^\epsilon$ and $f^\epsilon$ are given by the ansatz (75) and (76), respectively, then the dust surface temperature $T_p^\epsilon$ satisfies the following equation:

$$ egin{aligned} 
\rho_F^\epsilon \tau \gamma (\partial_t T_p^\epsilon + u_F^\epsilon \cdot \nabla_x T_p^\epsilon) &= -\frac{4}{\delta\nu^2} \frac{1}{(1+\eta)} T_p^\epsilon \Upsilon_1(F^\epsilon, f^\epsilon) + \frac{1}{\delta\nu^2} \frac{1}{(1+\eta)} \Upsilon_2(F^\epsilon, f^\epsilon) \\ 
&= -\frac{1}{\delta\nu^2} (\Xi_F^\epsilon + \Xi_f^\epsilon). 
\end{aligned} 
\quad (87) $$

## 5.3 Hard-sphere collision kernel

To make the above calculations more explicit and give them a physical meaning, we now focus on the case of hard-sphere collisional kernel. This will also be useful for performing the asymptotic diffusive limit later on. In a model with a hard-sphere collision kernel, we assume that a particle moves in a straight line until it collides with another particle. In this case, we hence have the following explicit expression of $\zeta$:

$$ \zeta(v-w,n) := [(v-w) \cdot n]_+ = ((v-w) \cdot n) 1_{\{(v-w) \cdot n \ge 0\}}. $$

Therefore, for $y \in \mathbb{R}^3$,

$$ \int_{S^2} \zeta(y, n) \, dn = \pi|y|, \quad \int_{S^2} n\zeta(y, n) \, dn = \frac{2\pi}{3}y. \quad (88) $$

The proofs of these results can be found, for example, in [13].

**Proposition 10.** For a hard-sphere collision kernel, the expressions of the operators $\Theta_F^\epsilon$ and $\Theta_f^\epsilon$ of Proposition 8 become

$$ \Theta_F^\epsilon = \frac{\sqrt{\eta}}{1+\eta} \pi \left[ \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} (w - \sqrt{\eta}v) F^\epsilon(v) f^\epsilon(w) |\sqrt{\eta}v - w| \mathrm{dwdv} + \frac{\gamma}{3} \sqrt{2\pi T_p^\epsilon} \rho_F^\epsilon \rho_f^\epsilon (u_f^\epsilon - u_F^\epsilon) \right], \quad (89) $$

and

$$ \Theta_f^\epsilon = \frac{1}{1+\eta} \pi \left[ \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} (\sqrt{\eta}v - w) F^\epsilon(v) f^\epsilon(w) |\sqrt{\eta}v - w| \mathrm{dwdv} + \frac{\gamma}{3} \sqrt{2\pi T_p^\epsilon} \rho_F^\epsilon \rho_f^\epsilon (u_F^\epsilon - u_f^\epsilon) \right]. \quad (90) $$

Furthermore, the operators $\Xi_F^\epsilon$ and $\Xi_f^\epsilon$, expressed for a general collision kernel in Proposition 9,
can be written in the hard-sphere case as

$$
\begin{equation}
\Xi_F^\epsilon = \frac{\pi}{(1+\eta)^2} \left[
\begin{aligned}
& (-2\eta - \eta^2) \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} |v|^2 F^\epsilon(v) f^\epsilon(w) |\sqrt{\eta v} - w| \, dw dv \\
& + \eta \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} |w|^2 F^\epsilon(v) f^\epsilon(w) |\sqrt{\eta v} - w| \, dw dv \\
& + 4\eta T_p \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} F^\epsilon(v) f^\epsilon(w) |\sqrt{\eta v} - w| \, dw dv \\
& + 2\sqrt{\eta} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} v \cdot w F^\epsilon(v) f^\epsilon(w) |\sqrt{\eta v} - w| \, dw dv \\
& - \frac{2}{3}\pi \sqrt{2\pi T_p} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} (\sqrt{\eta v} + \eta w) F^\epsilon(v) f^\epsilon(w) \cdot (\sqrt{\eta v} - w) \, dw dv
\end{aligned}
\right],
\end{equation}
$$

and

$$
\begin{equation}
\Xi_f^\epsilon = \frac{\pi}{(1+\eta)^2} \left[
\begin{aligned}
& (-1 - 2\eta) \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} |w|^2 F^\epsilon(v) f^\epsilon(w) |\sqrt{\eta v} - w| \, dw \, dv \\
& + \eta \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} |v|^2 F^\epsilon(v) f^\epsilon(w) |\sqrt{\eta v} - w| \, dw \, dv \\
& + 4T_p \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} F^\epsilon(v) f^\epsilon(w) |\sqrt{\eta v} - w| \, dw \, dv \\
& + 2\eta \sqrt{\eta} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} v \cdot w F^\epsilon(v) f^\epsilon(w) |\sqrt{\eta v} - w| \, dw \, dv \\
& + \frac{2}{3} \sqrt{2\pi T_p \pi} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} (\sqrt{\eta v} + \eta w) F^\epsilon(v) f^\epsilon(w) \cdot (\sqrt{\eta v} - w) \, dw \, dv
\end{aligned}
\right].
\end{equation}
$$

Proof. It is a direct consequence of applying (88) to the general expressions (79), (80), (83), and (84). □

**Proposition 11.** In the case of the hard-sphere collision kernel, the equation for $T_p^\epsilon$ becomes:

$$
\begin{align}
& \rho_F^\epsilon \frac{\tau}{\pi} \gamma \delta \nu^2 (1+\eta) (\partial_t T_p^\epsilon + u_F^\epsilon \cdot \nabla_x T_p^\epsilon) + 4T_p^\epsilon \iint F^\epsilon(v) f^\epsilon(w) |\sqrt{\eta v} - w| dv dw \tag{93} \\
& = \iint F^\epsilon(v) f^\epsilon(w) |\sqrt{\eta v} - w|^3 dv dw. \nonumber
\end{align}
$$

Proof. It is a direct consequence of applying (88) to $\Upsilon_1(F^\epsilon, f^\epsilon)$ and $\Upsilon_2(F^\epsilon, f^\epsilon)$. $\square$

6 Macroscopic equations and formal asymptotic

At this point, we are ready to explore the formal macroscopic limit as all the parameters approach
zero. The underlying idea is to adopt the viewpoint of the particles. As $\eta \to 0$, we interpret this
as the particles perceiving the gas molecules as much faster, yet much lighter and smaller, with
$\nu \to 0$. Specifically, in this limit, we expect the system to be described by a fluid model, as the
Knudsen number of the gas, $\delta$, approaches zero. Additionally, the Mach number $\gamma$ vanishes in the
diffusive asymptotic regime. To summarize, we take the limit:

$$
\epsilon = (\gamma, \nu, \delta, \eta) \to 0.
$$

We further assume that the macroscopic quantities (70), (71), (72) have a limit when $\epsilon \to 0$. We denote

$$
\rho_F(t, x) := \lim_{\epsilon \to 0} \rho_F^\epsilon(t, x), \quad \rho_f(t, x) := \lim_{\epsilon \to 0} \rho_f^\epsilon(t, x), \qquad (94)
$$

$$
u_F(t, x) := \lim_{\epsilon \to 0} u_F^\epsilon(t, x), \quad u_f(t, x) := \lim_{\epsilon \to 0} u_f^\epsilon(t, x), \qquad (95)
$$
$$ \theta_F(t, x) := \lim_{\epsilon \to 0} \theta_F^\epsilon(t, x), \qquad \theta_f(t, x) := \lim_{\epsilon \to 0} \theta_f^\epsilon(t, x). \qquad (96) $$

As a consequence, $f^\epsilon(t, x)$ and $F^\epsilon(t, x)$ have a limit when $\epsilon \to 0$, which we denote $f(t, x)$ and $F(t, x)$, respectively. Moreover, we also assume that the derivatives with respect to $t$ and $x_i$ (for $1 \le i \le d$) of $\rho_F^\epsilon, \rho_f^\epsilon, u_F^\epsilon, u_f^\epsilon, \theta_F^\epsilon, \theta_f^\epsilon$, and $T_p^\epsilon$ are bounded as $\epsilon \to 0$. Since the results presented here are formal, we will exchange limits and derivatives without further justification. We also assume that we have the following link between parameters.

**Hypothesis 5.** It is assumed that the Knudsen number of the dust particles, defined by (57), tends to 0 at the same rate as the Mach number of the gas $\gamma$, i.e.,

$$ \delta\nu^2 \sim \gamma, \qquad (97) $$

and that the Mach number of the dust particles, defined by (58), tends to zero, i.e.

$$ \frac{\gamma}{\sqrt{\eta}} \to 0, \quad \text{when } \epsilon \to 0. \qquad (98) $$

To simplify the link between the parameters $\eta$ and $\gamma$, we choose $\gamma$ as the reference parameter. Then, hypothesis (98) can equivalently be expressed as

$$ \eta \sim \gamma^s, \quad \text{with } s < 2. \qquad (99) $$

## 6.1 General collision kernel

We start by expressing the equations that can be obtained for any collision kernel.

**Lemma 6.** Under the previous assumption about the unknown functions and the Hypothesis (97), we have

$$ 4\Upsilon_1(F^\epsilon, f^\epsilon)T_p^\epsilon - \Upsilon_2(F^\epsilon, f^\epsilon) = \mathcal{O}(\gamma^2), \qquad (100) $$

and thus

$$ T_p^\epsilon(t, x) \to T_p := \frac{\Upsilon_2(F, f)}{4\Upsilon_1(F, f)}, \quad \text{when } \epsilon \to 0. \qquad (101) $$

*Proof.* Since the left side of the system equations (81) is bounded, Hypothesis (97) leads to $\Xi_f^\epsilon + \Xi_f^s = O(\gamma^2)$, and we conclude with (85). □

**Theorem 2.** Let $(\rho_F, u_F, \theta_F)$, $(\rho_f, u_f, \theta_f)$, and $T_p$ defined by (94), (95), (96), and (101). For general collision kernels, the macroscopic variables satisfy the conservation of mass equations

$$ \left\{ \begin{array}{l} \partial_t \rho_F + \nabla_x \cdot (\rho_F u_F) = 0, \\ \partial_t \rho_f + \nabla_x \cdot (\rho_f u_f) = 0, \end{array} \right. \qquad (102) $$

the equation

$$ \nabla_x (\rho_F \theta_F) + \nabla_x (\rho_f \theta_f) = 0, \qquad (103) $$

and the conservation of the total energy

$$ 3\partial_t (\rho_F \theta_F + \rho_f \theta_f) + 5\nabla_x \cdot (u_F \theta_F \rho_F + u_f \theta_f \rho_f) + 2\rho_F \tau (\partial_t T_p + u_F \cdot \nabla_x T_p) = 0. \qquad (104) $$

Moreover,

$$ \left\{ \begin{array}{l} \nabla_x (\rho_F \theta_F) = \Theta_F, \\ \nabla_x (\rho_f \theta_f) = \Theta_f, \end{array} \right. \qquad (105) $$
$$ \textit{and} \qquad \left\{ \begin{array}{l} \partial_t (3\rho_F \theta_F) + \nabla_x \cdot (5u_F \theta_F \rho_F) = \Xi_F, \\ \partial_t (3\rho_f \theta_f) + \nabla_x \cdot (5u_f \theta_f \rho_f) = \Xi_f, \end{array} \right. \qquad (106) $$

$$ \textit{with} \qquad \Theta_F := \lim_{\epsilon \to 0} \frac{1}{\delta\nu^2\sqrt{\eta}} \Theta_F^\epsilon, \qquad \Theta_f := \lim_{\epsilon \to 0} \frac{1}{\delta\nu^2} \Theta_f^\epsilon, $$

$$ \textit{and} \qquad \Xi_F := \lim_{\epsilon \to 0} \frac{1}{\gamma \delta \nu^2} \Xi_F^\epsilon, \qquad \Xi_f := \lim_{\epsilon \to 0} \frac{1}{\gamma \delta \nu^2} \Xi_f^\epsilon. $$

*Proof.* Equations (102) are a direct consequence of taking the limit in the mass balance equations (77), using that the derivatives with respect to $t$ and $x_i$ are bounded as $\epsilon \to 0$. By combining the two equations of the momentum balance (78), we obtain:

$$ \begin{aligned} & \frac{\gamma^2}{\eta} [\partial_t (\rho_F^\epsilon u_F^\epsilon) + \nabla_x \cdot (u_F^\epsilon \otimes u_F^\epsilon \rho_F^\epsilon)] + \gamma^2 [\partial_t (\rho_f^\epsilon u_f^\epsilon) + \nabla_x \cdot (u_f^\epsilon \otimes u_f^\epsilon \rho_f^\epsilon)] + \nabla_x (\rho_F^\epsilon \theta_F^\epsilon) + \nabla_x (\rho_f^\epsilon \theta_f^\epsilon) \\ &= \frac{1}{\delta\nu^2\sqrt{\eta}} \Theta_F^\epsilon + \frac{1}{\delta\nu^2} \Theta_f^\epsilon \\ &= \frac{1}{\delta\nu^2\sqrt{\eta}} \left( \int_{\mathbb{R}^3} v \mathcal{D}(F^\epsilon, f^\epsilon)(t, x, v) \, dv + \sqrt{\eta} \int_{\mathbb{R}^3} w \mathcal{R}(F^\epsilon, f^\epsilon)(t, x, w) \, dw \right) \\ &= 0 \end{aligned} $$

thanks to (69). Assumption (98) allows us to obtain (103) when $\epsilon \to 0$. By combining the two equations of (81) and using (87), we get

$$ \begin{aligned} & \frac{\gamma^2}{\eta} [\partial_t (|u_F^\epsilon|^2 \rho_F^\epsilon) + \nabla_x \cdot (|u_F^\epsilon|^2 u_F^\epsilon \rho_F^\epsilon)] + \gamma^2 [\partial_t (|u_f^\epsilon|^2 \rho_f^\epsilon) + \nabla_x \cdot (|u_f^\epsilon|^2 u_f^\epsilon \rho_f^\epsilon)] \\ & \quad + \partial_t (3\rho_F^\epsilon \theta_F^\epsilon) + \nabla_x \cdot (5u_F^\epsilon \theta_F^\epsilon \rho_F^\epsilon) + \partial_t (3\rho_f^\epsilon \theta_f^\epsilon) + \nabla_x \cdot (5u_f^\epsilon \theta_f^\epsilon \rho_f^\epsilon) \\ &= \frac{1}{\gamma \delta \nu^2} (\Xi_F^\epsilon + \Xi_f^\epsilon) \\ &= -\rho_F^\epsilon \tau (\partial_t T_p^\epsilon + u_F^\epsilon \cdot \nabla_x T_p^\epsilon). \end{aligned} $$

Assumption (98) allows us to obtain (104) when $\epsilon \to 0$. Equations (105) and (106) respectively are a direct consequence of taking the limit in each equation of (78) and (81) respectively, under assumption (98). $\square$

**Remark 9.** Note that we have in particular the following relation

$$ \Xi_f + \Xi_F = -2\tau\rho_F(\partial_t T_p + u_F \cdot \nabla_x T_p). $$

**Remark 10.** The results of Theorem 2 allowed us to obtain 14 equations for 19 unknowns: $T_p$, $\rho_F$, $\rho_f$, $\theta_F$, $\theta_f$, $\Xi_F$, $\Xi_f$ which take real values, and $u_F$, $u_f$, $\Theta_F$, $\Theta_f$ which take values in $\mathbb{R}^3$. It is possible to derive additional equations explicitly by specifying the form of the collision kernel.

## 6.2 Hard-Sphere collision kernel

The choice of a specific kernel allows us to deduce other independent relationships between the unknowns of the limit problem. In the physically relevant case of the hard-sphere dynamics, the results are as follows. To proceed, we make the subsequent additional assumption.

**Hypothesis 6.** We assume

$$ \frac{\eta^3}{\gamma} = o(1), \quad \text{when} \quad \epsilon \to 0. \qquad (107) $$
Taking $\gamma$ as the reference parameter, hypotheses (98) and (107) can be combined in the following relation

$$ \eta = \gamma^s, \quad \text{with } s \in (\frac{1}{3}, 2). \qquad (108) $$

**Proposition 12.** In the case of the hard-sphere collision kernel, under Hypothesis (97) and (107), we have:

$$ \left\{ \begin{aligned} \Theta_F &= \rho_F \rho_f \frac{\sqrt{\pi}}{3\sqrt{2}} (16\sqrt{\theta_f} + 2\pi\sqrt{T_p}) (u_f - u_F), \\ \Theta_f &= \rho_F \rho_f \frac{\sqrt{\pi}}{3\sqrt{2}} (16\sqrt{\theta_f} + 2\pi\sqrt{T_p}) (u_F - u_f). \end{aligned} \right. \qquad (109) $$

Furthermore, in the limit, the following temperature equality holds:

$$ \theta_f = \theta_F = T_p. $$

*Proof.* We begin by applying the asymptotic expansions (123) for $\Theta_F^\epsilon/(\eta\gamma)$ and $\Theta_f^\epsilon/\gamma$, as given in Proposition D.1 of Appendix D. These expansions lead to (109) when $\epsilon \to 0$, under Hypotheses (98) and (107). Next, we utilize the asymptotic expansion (129) for $\Xi_F^\epsilon$ provided in Proposition D.2 of Appendix D, to obtain, according to the reference parameter $\gamma$:

$$ \frac{\Xi_f^\epsilon}{\gamma^2} = \frac{1}{\gamma^2} \frac{16\pi}{\sqrt{2\pi}} \rho_F^\epsilon \rho_f^\epsilon \sqrt{\theta_f^\epsilon} (T_p^\epsilon - \theta_f^\epsilon) + O(\gamma^{s-2}) $$

Since we assume the boundedness of the derivative of the macroscopic quantities in the first equation of (81) and Hypothesis (97), the term $\Xi_f^\epsilon/\gamma^2$ should be bounded. This implies that the term

$$ \frac{16\pi}{\sqrt{2\pi}} \rho_F^\epsilon \rho_f^\epsilon \sqrt{\theta_f^\epsilon} (T_p^\epsilon - \theta_f^\epsilon) $$

must tend to zero when $\epsilon \to 0$, and we get $T_p = \theta_f$. We then use the asymptotic expansion (128) for $\Xi_F^\epsilon$, which gives:

$$ \frac{\Xi_F^\epsilon}{\gamma^2} = \frac{1}{\gamma^{2-s}} \frac{16\pi}{\sqrt{2\pi}} \rho_F^\epsilon \rho_f^\epsilon \sqrt{\theta_f^\epsilon} \left( -2\theta_F^\epsilon + \theta_f^\epsilon + T_p^\epsilon - \frac{\pi^2}{4} \frac{\sqrt{T_p^\epsilon}}{\sqrt{\theta_f^\epsilon}} (\theta_F^\epsilon - \theta_f^\epsilon) \right) + O(\gamma^{2s-2}) $$

Similarly, the boundedness of the derivatives in the second equation of (81), along with Hypothesis (97), leads us to conclude that, since $2-s>0$, the term of order $\gamma^{s-2}$ must tend to zero as $\gamma \to 0$. Given that $T_p = \theta_f$, this implies that $\theta_F = \theta_f$. $\square$

**Remark 11.** Likewise, in the Maxwellian collision kernel case, we can explicitly compute $\Theta_F^\epsilon$, $\Theta_f^\epsilon$, $\Xi_F^\epsilon$, and $\Xi_f^\epsilon$. We get (see Proposition E.3 in Appendix E), under Hypothesis (97), the following result:

$$ \Theta_F = \lim_{\epsilon \to 0} \frac{1}{\delta\nu^2\sqrt{\eta}} \Theta_F^\epsilon = 2\pi \|b\|_1 \rho_f \rho_F (u_f - u_F), $$

$$ \Theta_f = \lim_{\epsilon \to 0} \frac{1}{\delta\nu^2} \Theta_f^\epsilon = 2\pi \|b\|_1 \rho_f \rho_F (u_F - u_f). $$

Moreover, in the limit, we have

$$ \theta_F = \theta_f = \frac{4}{3} T_p. $$

The last relation between temperatures does not align with what one would intuitively expect from physics. It suggests that the Maxwellian collision cross-section is a more theoretical model, somewhat more distant from physical reality compared to the hard-sphere model. It is worth noting, however, that this result is similar to the findings in [13], particularly in the case of a spatially homogeneous system:

$$ \left\{ \begin{aligned} \partial_t F &= \mathcal{D}(F, f), \\ \partial_t f &= \mathcal{R}(F, f) + C(f, f), \end{aligned} \right. $$
where the collision operators $D(F, f)$ and $R(F, f)$ are based on the diffuse reflection mechanism at fixed temperature (that is, $T_p$ constant in space and time). In the case of a Maxwellian collision kernel, the author proved that

$$ \lim_{t \to +\infty} \theta_F(t) = \lim_{t \to +\infty} \theta_f(t) = \frac{4}{3} T_p. $$

**Remark 12.** The asymptotic limit derived in Section 6 does not provide a fully closed system: as noted in Remark 10, the Maxwell-Stefan-type relations obtained from the kinetic model are rank-deficient, so some degrees of freedom remain undetermined. For hard-sphere kernels, additional independent relations can be extracted from the detailed structure of the collision operator (see Hypothesis 6), which reduces but does not completely eliminate this deficiency. In the case of Maxwellian molecules (Remark 11), two degrees of freedom remain because the leading-order asymptotics eliminate higher-order moments, leaving two more constraints unspecified. These missing conditions are usually supplied by macroscopic closures, such as an equation of state linking pressure, density, and temperature. Therefore, specifying the collision kernel alone is not sufficient to close the system; additional macroscopic constraints are required for well-posedness.

This behavior is consistent with the other results on this subject (see, e.g., [1, 2, 3, 8, 25, 24]) and with the need of closure relationships in the macroscopic derivation of the Maxwell-Stefan equations themselves [26].

# 7 Conclusion

In this paper, we have derived a kinetic model to describe a gas-dust aerosol mixture, capable of accounting for the macroscopic nature of the dust relative to the gas and, consequently, the energy exchange between them. At the microscopic level, this interaction is modeled as diffuse reflection-type collisions between the gas and the dust, occurring on the surface of the dust particles. A surface temperature, $T_p = T_p(t, x)$, is introduced as a macroscopic variable, dependent on both time and space.

We formulate a Boltzmann-type system with elastic collision operators for mono-species interactions and an inelastic diffuse reflection operator for bi-species interactions. The system is closed by a transport-like equation (36) for the temperature $T_p$. This approach ensures the conservation of the system's total energy. In the case of a hard-sphere collision kernel, the model additionally possesses an explicit formulation of an entropy functional.

The system has been adimensionalized, and the key dimensionless quantities have been identified. The relationships between the various parameters have been made explicit by taking $\gamma$ as the reference parameter. The following table summarizes these relationships and the meaning of each parameter:

<table><caption>Table 2: Summary of the links between the small parameters. The assumption $s > 1/3$ is used only in the hard-sphere case.</caption><thead><tr><th>Parameter</th><th>Physical meaning</th><th>Scaling with respect to &gamma;</th></tr></thead><tbody><tr><td>&gamma;</td><td>Mach number of the gas</td><td>&mdash;</td></tr><tr><td>&delta;</td><td>Knudsen number of the gas</td><td>&mdash;</td></tr><tr><td>&nu;</td><td>Ratio of molecular radius</td><td>&mdash;</td></tr><tr><td>&eta;</td><td>Mass ratio</td><td>&eta; &sim; &gamma;<sup>s</sup>, s &isin; (<sup>1</sup>&frasl;<sub>3</sub>, 2)</td></tr><tr><td>&gamma;<sub>p</sub></td><td>Mach number of the dust</td><td>&gamma;<sub>p</sub> = &gamma;/&radic;[&eta;] &sim; &gamma;<sup>1-s/2</sup></td></tr><tr><td>&delta;<sub>p</sub></td><td>Knudsen number of dust</td><td>&delta;<sub>p</sub> &sim; &delta;&nu;<sup>2</sup> &sim; &gamma;</td></tr></tbody></table>

By letting these parameters tend to zero, we obtained a coupled system of partial differential equations describing the diffusive limit of the gas-particle mixture. In the hard-sphere setting, we have formally shown that the kinetic temperatures of the species and the surface temperature of the particles converge to the same limit function.
As a perspective, future work may focus on rigorously justifying the diffusion asymptotic and on performing a comprehensive analytical and numerical study of both the Boltzmann and macroscopic systems. Moreover, extending the current results to a multi-species setting would provide a more accurate representation of gas-dust mixtures, which in real-world applications typically involve more than two components.

# A Moments of $h_n$

We summarize in the following some of the calculations used in this paper.

**Lemma A.1.** Given a function

$$h_n(t, x, z) = \frac{c^4}{2\pi} (n \cdot z) e^{-\frac{c^2}{2}|z|^2} 1_{\{z \cdot n \ge 0\}},$$

for $c = c(t, x) \in \mathbb{R}$ constant in $n$ and $z$, it holds that

$$\int_{\mathbb{R}^3} h_n(z) \, dz = 1. \qquad (110)$$

Moreover,

$$\int_{\mathbb{R}^3} zh_n(z) \, dz = \frac{1}{c} \sqrt{\frac{\pi}{2}} n, \qquad (111)$$

and

$$\int_{\mathbb{R}^3} |z|^2 h_n(z) \, dz = \frac{4}{c^2}. \qquad (112)$$

The interested lector can find the proof of these results in [13, Annexe C].

# B Proofs of Propositions 8 and 9

## B.1 Proof of Proposition 8

To deduce (78), we first multiply the first and second equations of (64) by $w_{(\ell)}$ and $v_{(\ell)}$, respectively, for $\ell = 1, 2, 3$, and then we integrate the first equation of (64) respect to $v$ and the second one respect to $w$. The resulting system is:

$$\left\{ \begin{aligned} & \frac{\gamma^2}{\eta} \partial_t (\rho_F^\epsilon(u_F^\epsilon)_{(\ell)}) + \nabla_x \cdot \left( \int_{\mathbb{R}^3} v_{(\ell)} F^\epsilon(v) v \, dv \right) = \frac{1}{\delta\nu^2\sqrt{\eta}} (\Theta_F^\epsilon)_{(\ell)}, \\ & \gamma^2 \partial_t (\rho_f^\epsilon(u_f^\epsilon)_{(\ell)}) + \nabla_x \cdot \left( \int_{\mathbb{R}^3} w_{(\ell)} f^\epsilon(w) w \, dw \right) = \frac{1}{\delta\nu^2} (\Theta_f^\epsilon)_{(\ell)}. \end{aligned} \right.$$

We first study the divergence term appearing in the equation for the dust. For $k = 1, 2, 3$ and $\ell = 1, 2, 3$, thanks to the change of variables $v \mapsto z := v - \gamma u_F^\epsilon / \sqrt{\eta}$, and by noting that the following integrals are zero

$$\left\{ \begin{aligned} & \int_{\mathbb{R}^3} v_{(k)} v_{(\ell)} \exp \left( -\frac{|v|^2}{2\theta_F^\epsilon} \right) \, dv = 0 && \text{for all } k \neq \ell, \\ & \int_{\mathbb{R}^3} v_{(k)} \exp \left( -\frac{|v|^2}{2\theta_F^\epsilon} \right) \, dv = 0 && \text{for all } k, \end{aligned} \right.$$

we have that:

$$\begin{aligned} \int_{\mathbb{R}^3} v_{(\ell)} F^\epsilon(v) v_{(k)} \, dv &= \left(\frac{1}{2\pi\theta_F^\epsilon}\right)^{\frac{3}{2}} \rho_F^\epsilon \int_{\mathbb{R}^3} v_{(\ell)} v_{(k)} \exp\left(-\frac{1}{2\theta_F^\epsilon}\left|v - \frac{\gamma}{\sqrt{\eta}}u_F^\epsilon\right|^2\right) \, dv \\ &= \left(\frac{1}{2\pi\theta_F^\epsilon}\right)^{\frac{3}{2}} \rho_F^\epsilon \int_{\mathbb{R}^3} \left(\frac{\gamma^2}{\eta}(u_F^\epsilon)_{(\ell)}(u_F^\epsilon)_{(k)} + v_{(\ell)}^2 \delta_{k\ell}\right) \exp\left(-\frac{|v|^2}{2\theta_F^\epsilon}\right) \, dv \\ &= \frac{\gamma^2}{\eta}(u_F^\epsilon)_{(\ell)}(u_F^\epsilon)_{(k)}\rho_F^\epsilon + \rho_F^\epsilon\theta_F^\epsilon. \end{aligned}$$
Regarding the divergence term appearing in the equation for the gas, following the same pro-
cedure, we deduce that

$$
\int_{\mathbb{R}^3} w_{(\ell)} f^\epsilon(w) w_{(k)} \, dw = \gamma^2 (u_f^\epsilon)_{(\ell)} (u_f^\epsilon)_{(k)} \rho_f^\epsilon + \rho_f^\epsilon \theta_f^\epsilon.
$$

Moreover, via the weak formulation (16), we can rewrite Θ<sub>F</sub><sup>ε</sup> as

$$
\Theta_F^\epsilon = \int_{\mathbb{R}^3} v \mathcal{D}(F^\epsilon, f^\epsilon) \, dv = \int_{\mathbb{S}^2} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} (v' - v) F^\epsilon(v) f^\epsilon(w) h_n(z) \zeta(\sqrt{\eta}v - w, n) \, dndzdwdv.
$$

Similarly, Θ<sub>f</sub><sup>ε</sup> can be treated in a parallel manner, starting by using the weak formulation (18)

$$
\Theta_f^\epsilon = \int_{\mathbb{R}^3} w \mathcal{R}(F^\epsilon, f^\epsilon) \, dw = \int_{\mathbb{S}^2} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} (w' - w) F^\epsilon(v) f^\epsilon(w) h_n(z) \zeta(\sqrt{\eta}v - w, n) \, dndzdwdv.
$$

Thanks to (63), we have for the rescaled variables:

$$
\left\{
\begin{aligned}
v' - v &= - \frac{\sqrt{\eta}}{1 + \eta} (\sqrt{\eta}v - w + z), \\
w' - w &= \frac{1}{1 + \eta} (\sqrt{\eta}v - w + z).
\end{aligned}
\right.
$$

Hence, the expressions (79) and (80) follows directly from computations (110) and (111) of Lemma A.1,
with

$$
c = \sqrt{\frac{1}{T_p(t, x)}}.
$$

B.2 Proof of Proposition 9

We now integrate the first equation of (64) respect to v, after having multiplied it by |v|^2, and the
second one respect to w, after having multiplied it by |w|^2. The system we obtain has the form:

$$
\left\{
\begin{aligned}
& \frac{\gamma}{\sqrt{\eta}} \partial_t \left( 3\rho_F^\epsilon \theta_F^\epsilon + \frac{\gamma^2}{\eta} |u_F^\epsilon|^2 \rho_F^\epsilon \right) + \nabla_x \cdot \left( \int_{\mathbb{R}^3} |v|^2 F^\epsilon(v) v \, dv \right) = \frac{1}{\delta\nu^2\sqrt{\eta}} \Xi_F^\epsilon, \\
& \gamma \partial_t \left( 3\rho_f^\epsilon \theta_f^\epsilon + \gamma^2 |u_f^\epsilon|^2 \rho_f^\epsilon \right) + \nabla_x \cdot \left( \int_{\mathbb{R}^3} |w|^2 f^\epsilon(w) w \, dw \right) = \frac{1}{\delta\nu^2} \Xi_f^\epsilon,
\end{aligned}
\right.
\quad (113)
$$

As in Proof B.1, we begin by studying the divergence term that appears in the dust equation,
utilizing the definition (75) for $F^\epsilon(v)$. Hence, for every $k = 1, 2, 3$, we perform the change of
variable $z := v - \frac{\gamma}{\sqrt{\eta}} u_F^\epsilon$ and rename $z$ as $v$. Then, using the properties of the Maxwellian, we
deduce that

$$
\begin{align*}
\int_{\mathbb{R}^3} |v|^2 F^\epsilon(v) v_{(k)} \, dv &= \frac{\gamma}{\sqrt{\eta}} (u_F^\epsilon)_{(k)} \rho_F^\epsilon \left( \frac{1}{2\pi \theta_F^\epsilon} \right)^{\frac{3}{2}} \int_{\mathbb{R}^3} |v|^2 \exp\left(-\frac{|v|^2}{2\theta_F^\epsilon}\right) \, dv \\
&\quad + \frac{\gamma^3}{\eta^{3/2}} |u_F^\epsilon|^2 (u_F^\epsilon)_{(k)} \rho_F^\epsilon \left( \frac{1}{2\pi \theta_F^\epsilon} \right)^{\frac{3}{2}} \int_{\mathbb{R}^3} \exp\left(-\frac{|v|^2}{2\theta_F^\epsilon}\right) \, dv \\
&\quad + 2 \frac{\gamma}{\sqrt{\eta}} \sum_{\ell=1}^3 (u_F^\epsilon)_{(\ell)} \rho_F^\epsilon \left( \frac{1}{2\pi \theta_F^\epsilon} \right)^{\frac{3}{2}} \int_{\mathbb{R}^3} (v_{(\ell)})^2 \delta_{k\ell} \exp\left(-\frac{|v|^2}{2\theta_F^\epsilon}\right) \, dv \\
&= 5 \frac{\gamma}{\sqrt{\eta}} (u_F^\epsilon)_{(k)} \theta_F^\epsilon \rho_F^\epsilon + \frac{\gamma^3}{\eta^{3/2}} |u_F^\epsilon|^2 (u_F^\epsilon)_{(k)} \rho_F^\epsilon,
\end{align*}
$$

where the last equality is simply a result of using the expressions (73) for the moments of $F^\epsilon(v)$.
In a similar way, we can express the divergence term appearing in the equation for the gas. We
obtain (81) by multiplying the first equation of (113) by $\sqrt{\eta}/\gamma$ and the second one by $1/\gamma$.
The operators $\Xi_F^\epsilon$ and $\Xi_f^\epsilon$ can be rewritten using the weak formulations (16) and (18), respectively. In particular,

$$ \Xi_F^\epsilon = \int_{\mathbb{R}^3} |v|^2 \mathcal{D}(F^\epsilon, f^\epsilon)(v) \, dv = \int_{\mathbb{S}^2} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} (|v'|^2 - |v|^2) F^\epsilon(v) f^\epsilon(w) h_n(z) \varsigma(\sqrt{\eta}v - w, n) \, dndzdwdv, $$

and

$$ \Xi_f^\epsilon = \int_{\mathbb{R}^3} |w|^2 \mathcal{R}(F^\epsilon, f^\epsilon)(w) \, dw = \int_{\mathbb{S}^2} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} (|w'|^2 - |w|^2) F^\epsilon(v) f^\epsilon(w) h_n(z) \varsigma(\sqrt{\eta}v - w, n) \, dndzdwdv. $$

Using the dimensionless version of relation (28), which is

$$ |v'|^2 + |w'|^2 - |v|^2 - |w|^2 = \frac{1}{1 + \eta} (|z|^2 - |\sqrt{\eta}v - w|^2), $$

we get

$$ \begin{aligned} \Xi_F^\epsilon + \Xi_f^\epsilon &= \frac{1}{1 + \eta} \int_{\mathbb{S}^2} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} (|z|^2 - |\sqrt{\eta}v - w|^2) F^\epsilon(v) f^\epsilon(w) h_n(z) \varsigma(\sqrt{\eta}v - w, n) \, dndzdwdv. \end{aligned} $$

It directly gives (85), using (112) of Lemma A.1, with $c = \sqrt{\frac{1}{T_p(t,x)}}$. Via (63), we know that for the rescaled variables

$$ \left\{ \begin{aligned} |v'|^2 - |v|^2 &= \frac{1}{(1 + \eta)^2} (\eta|w|^2 + \eta|z|^2 + 2\sqrt{\eta}v \cdot w - 2\sqrt{\eta}v \cdot z - 2\eta w \cdot z) + \left(\frac{1}{(1 + \eta)^2} - 1\right) |v|^2, \\ |w'|^2 - |w|^2 &= \frac{1}{(1 + \eta)^2} (\eta|v|^2 + |z|^2 + 2\sqrt{\eta\eta}v \cdot w + 2\sqrt{\eta\eta}v \cdot z + 2\eta w \cdot z) + \left(\frac{\eta^2}{(1 + \eta)^2} - 1\right) |w|^2. \end{aligned} \right. $$

Hence, thanks to computations (110), (111) and (112) of Lemma A.1, with $c = \sqrt{\frac{1}{T_p(t,x)}}$, we can obtain (83) and (84).

### C Functions $q$ and $q_k$, $k \in \{1, 2, 3, 4\}$

Here, we introduce the following functions, which will be used later in Appendix D:

$$ q(a) := \int_{\mathbb{R}^3} (a - y) |a - y| e^{-\frac{|y|^2}{2}} \frac{dy}{(2\pi)^{3/2}}, \qquad (114) $$

and for $k \in \{1, 2, 3, 4\}$,

$$ q_k(a) := \int_{\mathbb{R}^3} |a - y|^k e^{-\frac{|y|^2}{2}} \frac{dy}{(2\pi)^{3/2}}. $$

The expression of $q_2(a)$ and $q_4(a)$ can be explicitly computed as

$$ q_2(a) = \int_{\mathbb{R}^3} |a - y|^2 e^{-\frac{|y|^2}{2}} \frac{dy}{(2\pi)^{3/2}} = |a|^2 + 3, \qquad (115) $$

and

$$ q_4(a) = \int_{\mathbb{R}^3} |a - y|^4 e^{-\frac{|y|^2}{2}} \frac{dy}{(2\pi)^{3/2}} = |a|^4 + 10|a|^2 + 15. \qquad (116) $$

Moreover, it is possible to establish the following result for these functions (see [19], Lemmas A.3 and A.4):
**Lemma C.1.** For all $a \in \mathbb{R}^3$,

$$q(a) = a\bar{q}(|a|), \quad q_1(a) = \tilde{q}_1(|a|),$$

with

$$\bar{q}(r) = \frac{1}{\sqrt{2\pi}} \left\{ 2rI_2(r) + \frac{4}{3}\frac{I_4(r)}{r} - \frac{2}{15}\frac{I_6(r)}{r^3} + \frac{8}{15}r^2J_1(r) + \frac{8}{3}J_3(r) \right\},$$

$$\tilde{q}_1(r) = \frac{1}{\sqrt{2\pi}} \left\{ 2rI_2(r) + \frac{2}{3}\frac{I_4(r)}{r} + 2J_3(r) + \frac{2}{3}r^2J_1(r) \right\},$$

where (for $k \in \mathbb{N}$, $x \in \mathbb{R}$)

$$I_k(x) = \int_0^x t^k e^{-t^2/2} dt, \quad J_k(x) = \int_x^{+\infty} t^k e^{-t^2/2} dt.$$

Moreover,

$$p_3(a) := \int_{\mathbb{R}^3} |y|^2 |a - y| e^{-\frac{|y|^2}{2}} \frac{dy}{(2\pi)^{3/2}} = \tilde{p}_3(|a|),$$

with

$$\tilde{p}_3(r) = \frac{1}{\sqrt{2\pi}} \left\{ 2rI_4(r) + \frac{2}{3}\frac{I_6(r)}{r} + 2J_5(r) + \frac{2}{3}r^2J_3(r) \right\}.$$

The following equalities are also valid:

$$J_1(0) = 1, \quad J_3(0) = 2, \quad J_5(0) = 8. \qquad (117)$$

We refine and clarify these statements as detailed below.

**Corollary C.1.** There exists some functions $\bar{q}, \bar{\bar{q}}, \bar{q}_1, \bar{\bar{q}}_3 \in C^\infty(\mathbb{R})$ such that

$$q(a) = \frac{1}{\sqrt{2\pi}} a \left( \frac{16}{3} + \frac{8}{15} |a|^2 + |a|^4 \bar{q}(|a|) \right), \qquad (118)$$

$$q_1(a) = \frac{1}{\sqrt{2\pi}} \left( 4 + \frac{2}{3} |a|^2 + |a|^4 \bar{\bar{q}}_1(|a|) \right), \qquad (119)$$

$$q_3(a) = \frac{1}{\sqrt{2\pi}} (16 + 8|a|^2 + |a|^4 \bar{\bar{q}}_3(|a|)). \qquad (120)$$

*Proof.* We first note that

$$I_k(x) = \int_0^x t^k e^{-t^2/2} dt = x^{k+1} \int_0^1 u^k e^{-x^2 u^2/2} du = x^{k+1} \tilde{I}_k(x),$$

where $\tilde{I}_k \in C^\infty(\mathbb{R})$, and that

$$J_k(x) = J_k(0) - I_k(x).$$

To show (118), (119), and (120), we use the reformulations of the functions $q$, $q_1$, and $p_3$ from Lemma C.1. We then observe that $\bar{q}, \tilde{q}_1$, and $\tilde{p}_3$ can be rewritten as

$$\bar{q}(r) = \frac{1}{\sqrt{2\pi}} \left\{ \frac{8}{3} J_3(0) + \frac{8}{15} J_1(0) r^2 + r^4 \left( 2 \tilde{I}_2(r) + \frac{4}{3} \tilde{I}_4(r) - \frac{2}{15} \tilde{I}_6(r) - \frac{8}{15} \tilde{I}_1(r) - \frac{8}{3} \tilde{I}_3(r) \right) \right\},$$

$$\tilde{q}_1(r) = \frac{1}{\sqrt{2\pi}} \left\{ 2J_3(0) + \frac{2}{3}J_1(0)r^2 + r^4 \left( 2\tilde{I}_2(r) + \frac{2}{3}\tilde{I}_4(r) - 2\tilde{I}_3(r) - \frac{2}{3}\tilde{I}_1(r) \right) \right\},$$

$$\tilde{p}_3(r) = \frac{1}{\sqrt{2\pi}} \left\{ 2J_5(0) + \frac{2}{3}J_3(0)r^2 + r^6 \left( 2\tilde{I}_4(r) + \frac{2}{3}\tilde{I}_6(r) - 2\tilde{I}_5(r) - \frac{2}{3}\tilde{I}_3(r) \right) \right\}.$$

Finally, we rely $q_3(a)$ to $p_3(a)$, $q(a)$ and $q_1(a)$:

$$q_3(a) = p_3(a) + 2a \cdot q(a) - |a|^2 q_1(a).$$
We deduce that $q_3$ writes $q_3(a) = \tilde{q}_3(|a|)$, where

$$ 
\begin{aligned}
\tilde{q}_3(r) &= \tilde{p}_3(r) + 2r^2\bar{q}(r) - r^2\tilde{q}_1(r) \\
&= \frac{1}{\sqrt{2\pi}} \left\{ 2J_5(0) + 4J_3(0)r^2 + \frac{2}{5}J_1(0)r^4 \right. \\
&\qquad \left. + r^6 \left( -\frac{2}{5}\tilde{I}_1(r) + 2\tilde{I}_2(r) - 4\tilde{I}_3(r) + 4\tilde{I}_4(r) - 2\tilde{I}_5(r) + \frac{2}{5}\tilde{I}_6(r) \right) \right\}.
\end{aligned}
$$

By using the values given by (117), we can conclude. $\square$

## D Computations for the Hard-Sphere Kernel

The goal of this Appendix is to derive the asymptotic expansions of

$$ \frac{\Theta_F^\epsilon}{\sqrt{\eta\gamma}}, \quad \frac{\Theta_f^\epsilon}{\gamma}, \quad \text{and} \quad \frac{(\Xi_F^\epsilon + \Xi_f^\epsilon)}{\gamma^2} $$

with respect to the small parameter $\gamma$, assuming that $\delta\nu^2 \sim \gamma$ (Hypothesis (97)) and $\eta = \gamma^s$ with $s \in (\frac{1}{3}, 2)$ (Hypothesis (108)), in the context of the hard-sphere collision kernel. This procedure involves inserting the expressions for $F^\epsilon$ and $f^\epsilon$ from (75) into the integrals defining $\Theta_F^\epsilon$, $\Theta_f^\epsilon$, $\Xi_F^\epsilon$, and $\Xi_f^\epsilon$, as given in (89), (90), (91), and (92). The resulting expressions are then rewritten in terms of the functions $q$, $q_1$, and $q_3$.

Consider an integral of the form

$$ 
\begin{aligned}
\mathcal{I}^\epsilon &:= \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} G_\epsilon(v, w) F^\epsilon(v) f^\epsilon(w) |w - \sqrt{\eta}v| dw dv \\
&= \frac{1}{(2\pi)^3} \frac{\rho_F^\epsilon}{(\theta_F^\epsilon)^{3/2}} \frac{\rho_f^\epsilon}{(\theta_f^\epsilon)^{3/2}} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} G_\epsilon(v, w) \exp\left(-\frac{|v - \frac{\gamma}{\sqrt{\eta}}u_F^\epsilon|^2}{2\theta_F^\epsilon}\right) \exp\left(-\frac{|w - \gamma u_f^\epsilon|^2}{2\theta_f^\epsilon}\right) |w - \sqrt{\eta}v| dw dv.
\end{aligned}
$$

Different choices of the function $G_\epsilon$ will give different integrals, which appear in the following. Our goal here is to rewrite $G_\epsilon$ in a way that highlights the small parameters. With the change of variables

$$ \xi = \frac{1}{\sqrt{\theta_F^\epsilon}}(v - \frac{\gamma}{\sqrt{\eta}}u_F^\epsilon), \quad \varpi = \frac{1}{\sqrt{\theta_f^\epsilon}}(w - \gamma u_f^\epsilon), $$

we get

$$ \mathcal{I}^\epsilon = \rho_F^\epsilon \rho_f^\epsilon \sqrt{\theta_f^\epsilon} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} G_\epsilon \left( \sqrt{\theta_F^\epsilon} \xi + \frac{\gamma}{\sqrt{\eta}} u_F^\epsilon, \sqrt{\theta_f^\epsilon} \varpi + \gamma u_f^\epsilon \right) |a_\epsilon(\xi) - \varpi| e^{-|\xi|^2/2} e^{-|\varpi|^2/2} \frac{d\xi}{(2\pi)^{3/2}} \frac{d\varpi}{(2\pi)^{3/2}}, \quad (121) $$

with

$$ a^\epsilon(\xi) = \frac{1}{\sqrt{\theta_f^\epsilon}} \left( \sqrt{\theta_F^\epsilon} \sqrt{\eta} \xi + \gamma (u_F^\epsilon - u_f^\epsilon) \right) = \frac{\sqrt{\theta_F^\epsilon}}{\sqrt{\theta_f^\epsilon}} \gamma^{s/2} (\xi - b^\epsilon), \quad (122) $$

$$ \text{where } b^\epsilon = \frac{1}{\sqrt{\theta_F^\epsilon}} \frac{\gamma}{\sqrt{\eta}} (u_f^\epsilon - u_F^\epsilon) = O(\gamma^{1-s/2}) = o(1). $$

Notice that the integrals of the following parts are such that can be expressed in terms of $q$, $q_1$, and $q_3$, as defined in Appendix C.
## D.1 Asymptotic expansions of $\Theta_F^\epsilon$ and $\Theta_f^\epsilon$

**Proposition D.1.** In the case of the hard-sphere kernel, under Hypothesis (108), we have the following asymptotic expansions

$$ \left\{ \begin{aligned} \frac{\Theta_F^\epsilon}{\sqrt{\eta}\gamma} &= \pi \rho_F^\epsilon \rho_f^\epsilon (u_f^\epsilon - u_F^\epsilon) \frac{1}{3\sqrt{2\pi}} \left( 16\sqrt{\theta_f^\epsilon} + 2\pi\sqrt{T_p^\epsilon} \right) + O(\gamma^{3s-1}) + o(1), \\ \frac{\Theta_f^\epsilon}{\gamma} &= \pi \rho_F^\epsilon \rho_f^\epsilon (u_F^\epsilon - u_f^\epsilon) \frac{1}{3\sqrt{2\pi}} \left( 16\sqrt{\theta_f^\epsilon} + 2\pi\sqrt{T_p^\epsilon} \right) + O(\gamma^{3s-1}) + o(1). \end{aligned} \right. \quad (123) $$

where $\Theta_F^\epsilon$ and $\Theta_f^\epsilon$ are defined by (89) and (90).

*Proof.* We can express expressions (89) and (90) of $\Theta_F^\epsilon$ and $\Theta_f^\epsilon$ as

$$ \Theta_F^\epsilon = \frac{\sqrt{\eta}}{1+\eta} \pi \left[ \mathcal{I}_2^\epsilon + \frac{\gamma}{3} \sqrt{2\pi T_p^\epsilon} \rho_F^\epsilon \rho_f^\epsilon (u_f^\epsilon - u_F^\epsilon) \right], $$

and

$$ \Theta_f^\epsilon = \frac{1}{1+\eta} \pi \left[ -\mathcal{I}_2^\epsilon + \frac{\gamma}{3} \sqrt{2\pi T_p^\epsilon} \rho_F^\epsilon \rho_f^\epsilon (u_F^\epsilon - u_f^\epsilon) \right]. $$

where

$$ \mathcal{I}_2^\epsilon := \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} (w - \sqrt{\eta}v) F^\epsilon(v) f^\epsilon(w) |w - \sqrt{\eta}v| dwdv. $$

According to (121), we have

$$ \mathcal{I}_2^\epsilon = -\rho_F^\epsilon \rho_f^\epsilon \theta_f^\epsilon \int q(a^\epsilon(\xi)) e^{-|\xi|^2/2} \frac{d\xi}{(2\pi)^{3/2}} $$

where $q$ is defined by (114). Using (118) and an asymptotic expansion of $\bar{q}$ around 0 we get an asymptotic expansion of $q(a^\epsilon(\xi)):

$$ \begin{aligned} q(a^\epsilon(\xi)) &= \frac{1}{\sqrt{2\pi}} a^\epsilon(\xi) \left( \frac{16}{3} + \frac{8}{15} |a^\epsilon(\xi)|^2 + |a^\epsilon(\xi)|^4 \bar{q}(|0|) + O(|a^\epsilon(\xi)|^5) \right), \\ &= \frac{1}{\sqrt{2\pi}} \frac{\sqrt{\theta_F^\epsilon}}{\sqrt{\theta_f^\epsilon}} \gamma^{s/2} \left( \frac{16}{3}(\xi - b^\epsilon) + \frac{8}{15} \gamma^s \frac{\theta_F^\epsilon}{\theta_f^\epsilon} (\xi - b^\epsilon) |\xi - b^\epsilon|^2 + \bar{q}(|0|) \gamma^{2s} \left( \frac{\theta_F^\epsilon}{\theta_f^\epsilon} \right)^2 (\xi - b^\epsilon) |\xi - b^\epsilon|^4 \right) + O(\gamma^{3s}) \end{aligned} \quad (124) $$

where $b^\epsilon = o(1)$ is defined by (122). Then

$$ \int q(a^\epsilon(\xi)) e^{-|\xi|^2/2} \frac{d\xi}{(2\pi)^{3/2}} = -\frac{1}{\sqrt{2\pi}} \frac{\sqrt{\theta_F^\epsilon}}{\sqrt{\theta_f^\epsilon}} \gamma^{s/2} \left( \frac{16}{3} + \frac{8}{3} \gamma^s \frac{\theta_F^\epsilon}{\theta_f^\epsilon} \right) b^\epsilon + O(\gamma^3) + O(\gamma^{2s+1}) + O(\gamma^{3s}), \quad (125) $$

since

$$ \int (\xi - b^\epsilon) |\xi - b^\epsilon|^2 e^{-|\xi|^2/2} \frac{d\xi}{(2\pi)^{3/2}} = -b^\epsilon (5 + |b^\epsilon|^2), \quad (126) $$

and

$$ \int (\xi - b^\epsilon) |\xi - b^\epsilon|^4 e^{-|\xi|^2/2} \frac{d\xi}{(2\pi)^{3/2}} = b^\epsilon g(b^\epsilon) = O(\gamma^{1-s/2}). \quad (127) $$

Then, by replacing $b^\epsilon$ with $\frac{1}{\sqrt{\theta_F^\epsilon}} \gamma^{1-s/2} (u_f^\epsilon - u_F^\epsilon)$, we get

$$ \frac{\mathcal{I}_2^\epsilon}{\gamma} = \rho_F^\epsilon \rho_f^\epsilon (u_f^\epsilon - u_F^\epsilon) \frac{\sqrt{\theta_f^\epsilon}}{\sqrt{2\pi}} \left( \frac{16}{3} + \frac{8}{3} \gamma^s \frac{\theta_F^\epsilon}{\theta_f^\epsilon} \right) + O(\gamma^{3s-1}) + o(1) $$

and we obtain (123).

□
## D.2 Asymptotic expansions of $\Xi_F^\epsilon$ and $\Xi_f^\epsilon$

**Proposition D.2.** Under Hypothesis (99), we have the following asymptotic expansions:

$$ \Xi_F^\epsilon = \gamma^s \frac{16\pi}{\sqrt{2\pi}} \rho_F^\epsilon \rho_f^\epsilon \sqrt{\theta_f^\epsilon} \left( -2\theta_F^\epsilon + \theta_f^\epsilon + T_p^\epsilon - \frac{\pi^2}{4} \frac{\sqrt{T_p^\epsilon}}{\sqrt{\theta_f^\epsilon}} (\theta_F^\epsilon - \theta_f^\epsilon) \right) + O(\gamma^{2s}) \quad (128) $$

and

$$ \Xi_f^\epsilon = \frac{16\pi}{\sqrt{2\pi}} \rho_F^\epsilon \rho_f^\epsilon \sqrt{\theta_f^\epsilon} (T_p^\epsilon - \theta_f^\epsilon) + O(\gamma^{2s}). \quad (129) $$

*Proof.* We rewrite the expressions (91) and (92) of $\Xi_F^\epsilon$ and $\Xi_f^\epsilon$ according to the following integrals:

$$ \mathcal{I}_1^\epsilon := \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} F^\epsilon(v) f^\epsilon(w) |w - \sqrt{\eta}v| \, dw \, dv. $$

$$ \mathcal{I}_4^\epsilon := \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} |v|^2 F^\epsilon(v) f^\epsilon(w) |w - \sqrt{\eta}v| \, dw \, dv $$

$$ \mathcal{I}_5^\epsilon := \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} |w|^2 F^\epsilon(v) f^\epsilon(w) |w - \sqrt{\eta}v| \, dw \, dv $$

$$ \mathcal{I}_6^\epsilon := \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} (\sqrt{\eta}v + \eta w) \cdot (\sqrt{\eta}v - w) F^\epsilon(v) f^\epsilon(w) \, dw \, dv. $$

$$ \mathcal{I}_7^\epsilon := \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} v \cdot w F^\epsilon(v) f^\epsilon(w) |w - \sqrt{\eta}v| \, dw \, dv. $$

We have:

$$ \Xi_F^\epsilon = \frac{\pi}{(1+\eta)^2} \left[ -\eta(2+\eta)\mathcal{I}_4^\epsilon + \eta\mathcal{I}_5^\epsilon + 4\eta T_p \mathcal{I}_1^\epsilon + 2\sqrt{\eta}\mathcal{I}_7^\epsilon - \frac{2}{3}\pi\sqrt{2\pi T_p \mathcal{I}_6^\epsilon} \right], \quad (130) $$

and

$$ \Xi_f^\epsilon = \frac{\pi}{(1+\eta)^2} \left[ (-1-2\eta)\mathcal{I}_5^\epsilon + \eta\mathcal{I}_4^\epsilon + 4T_p\mathcal{I}_1^\epsilon + 2\eta\sqrt{\eta}\mathcal{I}_7^\epsilon + \frac{2}{3}\pi\sqrt{2\pi T_p\mathcal{I}_6^\epsilon} \right]. \quad (131) $$

Moreover,

$$ \mathcal{I}_5^\epsilon = \mathcal{I}_3^\epsilon - \eta\mathcal{I}_4^\epsilon + 2\sqrt{\eta}\mathcal{I}_7^\epsilon, $$

where

$$ \mathcal{I}_3^\epsilon := \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} F^\epsilon(v) f^\epsilon(w) |w - \sqrt{\eta}v|^3 \, dw \, dv, $$

and then if we express $\Xi_F^\epsilon$ and $\Xi_f^\epsilon$ according to $\mathcal{I}_3^\epsilon$ instead of $\mathcal{I}_5^\epsilon$, we get

$$ \Xi_F^\epsilon = \frac{\pi}{(1+\eta)^2} \left[ -2\eta(1+\eta)\mathcal{I}_4^\epsilon + \eta\mathcal{I}_3^\epsilon + 4\eta T_p \mathcal{I}_1^\epsilon + 2\sqrt{\eta}(1+\eta)\mathcal{I}_7^\epsilon - \frac{2}{3}\pi\sqrt{2\pi T_p \mathcal{I}_6^\epsilon} \right], \quad (132) $$

$$ \Xi_f^\epsilon = \frac{\pi}{(1+\eta)^2} \left[ 2\eta(1+\eta)\mathcal{I}_4^\epsilon - (1+2\eta)\mathcal{I}_3^\epsilon + 4T_p\mathcal{I}_1^\epsilon - 2\sqrt{\eta}(1+\eta)\mathcal{I}_7^\epsilon + \frac{2}{3}\pi\sqrt{2\pi T_p \mathcal{I}_6^\epsilon} \right]. \quad (133) $$

We can compute directly $\mathcal{I}_6^\epsilon$:

$$ \mathcal{I}_6^\epsilon = \rho_F^\epsilon \rho_f^\epsilon \{3\eta(\theta_F^\epsilon - \theta_f^\epsilon) + \gamma^2 (|u_F^\epsilon|^2 - |u_f^\epsilon|^2 + (\eta-1)u_F^\epsilon \cdot u_f^\epsilon)\} \quad (134) $$

and we perform the asymptotic expansion of $\mathcal{I}_1^\epsilon$, $\mathcal{I}_3^\epsilon$, $\eta\mathcal{I}_4^\epsilon$, $\sqrt{\eta}\mathcal{I}_7^\epsilon$, rewritten as it was done for the generic $\mathcal{I}^\epsilon$ in (121), and by replacing $\eta$ with $\gamma^s$. Thanks to the reformulation of $q_1$ and $q_3$ deduced in Corollary C.1, we can perform the following expansion of $q_1$ and $q_3$ around 0

$$ q_1(a) = \frac{1}{\sqrt{2\pi}} \left( 4 + \frac{2}{3} |a|^2 + |a|^4 \bar{q}_1(0) \right) + O(|a|^5), \quad (135) $$

$$ q_3(a) = \frac{1}{\sqrt{2\pi}} (16 + 8|a|^2 + |a|^4 \bar{q}_3(0)) + O(|a|^5). $$
Applied to $a^\epsilon(\xi)$, it gives

$$ 
\begin{aligned}
q_1(a^\epsilon(\xi)) &= \frac{1}{\sqrt{2\pi}} \left( 4 + \frac{2}{3} \gamma^s \frac{\theta_F^\epsilon}{\theta_f^\epsilon} |\xi - b^\epsilon|^2 + \gamma^{2s} \left( \frac{\theta_F^\epsilon}{\theta_f^\epsilon} \right)^2 |\xi - b^\epsilon|^4 \bar{q}_1(0) \right) + O(\gamma^{5s/2}) \\
q_3(a^\epsilon(\xi)) &= \frac{1}{\sqrt{2\pi}} \left( 16 + 8\gamma^s \frac{\theta_F^\epsilon}{\theta_f^\epsilon} |\xi - b^\epsilon|^2 + \gamma^{2s} \left( \frac{\theta_F^\epsilon}{\theta_f^\epsilon} \right)^2 |\xi - b^\epsilon|^4 \bar{q}_3(0) \right) + O(\gamma^{5s/2})
\end{aligned}
\quad (136) $$

**Remark 13.** The rest of the asymptotic expansion is not uniform in the variable $\xi$, but depends on powers of $|\xi - b^\epsilon|$. It is therefore integrable for the measure $e^{-|\xi|^2/2}d\xi$.

Thus, we get

$$ 
\begin{aligned}
\mathcal{I}_1^\epsilon &= \rho_F^\epsilon \rho_f^\epsilon \sqrt{\theta_f^\epsilon} \int_{\mathbb{R}^3} q_1(a^\epsilon(\xi)) e^{-|\xi|^2/2} \frac{d\xi}{(2\pi)^{3/2}} \\
&= \rho_F^\epsilon \rho_f^\epsilon \sqrt{\theta_f^\epsilon} \frac{1}{\sqrt{2\pi}} \left( 4 + \frac{2}{3} \gamma^s \frac{\theta_F^\epsilon}{\theta_f^\epsilon} \int |\xi - b^\epsilon|^2 e^{-|\xi|^2/2} \frac{d\xi}{(2\pi)^{3/2}} \right) + O(\gamma^{2s}) \\
&= \rho_F^\epsilon \rho_f^\epsilon \sqrt{\theta_f^\epsilon} \frac{1}{\sqrt{2\pi}} \left( 4 + 2\gamma^s \frac{\theta_F^\epsilon}{\theta_f^\epsilon} \right) + O(\gamma^{2s}) + O(\gamma^2).
\end{aligned}
\quad (137) $$

Similarly,

$$ 
\begin{aligned}
\mathcal{I}_3^\epsilon &= \rho_F^\epsilon \rho_f^\epsilon \theta_f^\epsilon^{3/2} \int_{\mathbb{R}^3} q_3(a_\epsilon(\xi)) e^{-|\xi|^2/2} \frac{d\xi}{(2\pi)^{3/2}} \\
&= \rho_F^\epsilon \rho_f^\epsilon \theta_f^\epsilon^{3/2} \frac{1}{\sqrt{2\pi}} \left( 16 + 24\gamma^s \frac{\theta_F^\epsilon}{\theta_f^\epsilon} \right) + O(\gamma^{2s}) + O(\gamma^2).
\end{aligned}
\quad (138) $$

Then we use again (136) to get the asymptotic expansion of $\eta \mathcal{I}_4^\epsilon = \gamma^s \mathcal{I}_4^\epsilon$:

$$ 
\begin{aligned}
\frac{\mathcal{I}_4^\epsilon}{\rho_F^\epsilon \rho_f^\epsilon} &= \sqrt{\theta_f^\epsilon} \theta_F^\epsilon \int_{\mathbb{R}^3} |\xi|^2 q_1(a_\epsilon(\xi)) e^{-|\xi|^2/2} \frac{d\xi}{(2\pi)^{3/2}} \\
&\quad + \sqrt{\theta_f^\epsilon} \gamma^{2-s} |u_F^\epsilon|^2 \int_{\mathbb{R}^3} q_1(a_\epsilon(\xi)) e^{-|\xi|^2/2} \frac{d\xi}{(2\pi)^{3/2}} \\
&\quad + 2\sqrt{\theta_f^\epsilon} \sqrt{\theta_F^\epsilon} \gamma^{1-s} u_F^\epsilon \cdot \int_{\mathbb{R}^3} \xi q_1(a_\epsilon(\xi)) e^{-|\xi|^2/2} \frac{d\xi}{(2\pi)^{3/2}} \\
&= \frac{\sqrt{\theta_f^\epsilon} \theta_F^\epsilon}{\sqrt{2\pi}} \left( 12 + \frac{2}{3} \gamma^s \frac{\theta_F^\epsilon}{\theta_f^\epsilon} \left( 15 + 3 \frac{1}{\theta_F^\epsilon} \gamma^{2-s} |u_F^\epsilon - u_F^\epsilon|^2 \right) \right) + O(\gamma^{2s}) \\
&\quad + \sqrt{\theta_f^\epsilon} \gamma^{2-s} |u_F^\epsilon|^2 \frac{1}{\sqrt{2\pi}} \left( 4 + 2\gamma^s \frac{\theta_F^\epsilon}{\theta_f^\epsilon} + \frac{2}{3} \frac{1}{\theta_f^\epsilon} \gamma^2 |u_F^\epsilon - u_F^\epsilon|^2 + O(\gamma^{2s}) \right) \\
&\quad + 2 \frac{1}{\sqrt{2\pi}} \sqrt{\theta_f^\epsilon} \sqrt{\theta_F^\epsilon} \gamma^{1-s} u_F^\epsilon \cdot \left( -\frac{4}{3} \gamma^s \frac{\theta_F^\epsilon}{\theta_f^\epsilon} \frac{1}{\sqrt{\theta_F^\epsilon}} \gamma^{1-s/2} (u_F^\epsilon - u_F^\epsilon) + O(\gamma^{3s/2+1}) + O(\gamma^{5s/2}) \right),
\end{aligned}
$$

and finally

$$ \gamma^s \frac{\mathcal{I}_4^\epsilon}{\rho_F^\epsilon \rho_f^\epsilon} = \frac{12}{\sqrt{2\pi}} \gamma^s \sqrt{\theta_f^\epsilon} \theta_F^\epsilon + O(\gamma^{2+s/2}) + O(\gamma^{2s}). \quad (139) $$

We use again the asymptotic expansions (124) of $q(a^\epsilon)$ and (135) of $q_1(a^\epsilon)$ to expand $\sqrt{\eta \mathcal{I}_7^\epsilon} =$
$\gamma^{s/2} \mathcal{I}_7^\epsilon :$

$$
\begin{align*}
\gamma^{s/2} \frac{\mathcal{I}_7^\epsilon}{\rho_F^\epsilon \rho_f^\epsilon}
&= \sqrt{\theta_f^\epsilon} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} \left( \gamma^{s/2} \sqrt{\theta_F^\epsilon} \xi + \gamma u_F^\epsilon \right) \cdot \left( \sqrt{\theta_f^\epsilon} \varpi + \gamma u_f^\epsilon \right) |a_\epsilon(\xi) - \varpi| e^{-|\xi|^2/2} e^{-|\varpi|^2/2} \frac{d\xi}{(2\pi)^{3/2}} \frac{d\varpi}{(2\pi)^{3/2}} \\
&= \sqrt{\theta_f^\epsilon} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} \left( \gamma^{s/2} \sqrt{\theta_F^\epsilon} (\xi - b^\epsilon) + \gamma u_f^\epsilon \right) \cdot \left( \sqrt{\theta_f^\epsilon} (\varpi - a^\epsilon) + \sqrt{\theta_f^\epsilon} a^\epsilon + \gamma u_f^\epsilon \right) \\
&\qquad \times |a_\epsilon(\xi) - \varpi| e^{-|\xi|^2/2} e^{-|\varpi|^2/2} \frac{d\xi}{(2\pi)^{3/2}} \frac{d\varpi}{(2\pi)^{3/2}} \\
&= -\theta_f^\epsilon \int \left( \gamma^{s/2} \sqrt{\theta_F^\epsilon} (\xi - b^\epsilon) + \gamma u_f^\epsilon \right) q(a^\epsilon) e^{-|\xi|^2/2} \frac{d\xi}{(2\pi)^{3/2}} \\
&\qquad + \sqrt{\theta_f^\epsilon} \int \left( \gamma^2 |u_f^\epsilon|^2 + 2\gamma^{1+s/2} \sqrt{\theta_F^\epsilon} u_f^\epsilon \cdot (\xi - b^\epsilon) + \gamma^s \theta_F |\xi - b^\epsilon|^2 \right) q_1(a^\epsilon) e^{-|\xi|^2/2} \frac{d\xi}{(2\pi)^{3/2}}
\end{align*}
$$

and then

$$
\gamma^{s/2} \mathcal{I}_7^\epsilon = - \frac{4}{\sqrt{2\pi}} \gamma^s \rho_F^\epsilon \rho_f^\epsilon \sqrt{\theta_f^\epsilon} \theta_F^\epsilon + O(\gamma^2) + O(\gamma^{2s}). \quad (140)
$$

Finally, we use (139), (138), (137), (140) and (134) to get (128) and (129).

□

E The case of the Maxwellian kernel

For a Maxwellian molecule case, the collision kernel $\zeta$ depends on $v$, $w$ and $n$ only trough the
deviation angle $\theta \in [0, \pi]$ between $v-w$ and $n$ and it does not depend on $|v-w|$, i.e.

$$
\zeta(v-w,n) := b \left( \frac{v-w}{|v-w|} \cdot n \right) = b(\cos \theta),
$$

where $b : [-1, 1] \rightarrow \mathbb{R}^{+}$, and we assume $b \in L^{1}([-1, 1])$ and $b \ge 0$. Therefore,

$$
\Sigma(v-w) = 2\pi \|b\|_{L^1(-1,1)}, \quad (141)
$$

and

$$
\int_{S^2} n\zeta(v-w,n) \, dn = 0. \quad (142)
$$

The proofs of these results can be found for example in [8].

**Proposition E.1.** For a Maxwellian collision kernel, the expressions of operators $\Theta_F^\epsilon$ and $\Theta_f^\epsilon$ of Proposition 8 become

$$
\left\{
\begin{array}{l}
\Theta_F^\epsilon = 2\pi \|b\|_{L^1} \frac{\gamma\sqrt{\eta}}{1+\eta} \rho_f^\epsilon \rho_F^\epsilon (u_f^\epsilon - u_F^\epsilon), \\
\\
\Theta_f^\epsilon = 2\pi \|b\|_{L^1} \frac{\gamma}{1+\eta} \rho_f^\epsilon \rho_F^\epsilon (u_F^\epsilon - u_f^\epsilon).
\end{array}
\right.
\qquad (143)
$$

Moreover, the operators $\Xi_F^\epsilon$ and $\Xi_f^\epsilon$, expressed for a general collision kernel in Proposition 9, in
the Maxwellian case can be written as

$$
\Xi_F^\epsilon = \frac{2\pi \|b\|_{L^1}}{(1+\eta)^2} \rho_f^\epsilon \rho_F^\epsilon \left[ \eta (4T_p^\epsilon + 3\theta_f^\epsilon - 3\eta\theta_F^\epsilon) - 6\eta\theta_F^\epsilon \right. \\
\left. \qquad + \eta\gamma^2 (|u_f^\epsilon|^2 - |u_F^\epsilon|^2) + 2\gamma^2 (u_F^\epsilon \cdot u_f^\epsilon - |u_F^\epsilon|^2) \right], \quad (144)
$$
and

$$ \Xi_f^\epsilon = \frac{2\pi\|b\|_{L^1}}{(1+\eta)^2} \rho_f^\epsilon \rho_F^\epsilon \left[ \begin{aligned} & -6\eta\theta_f^\epsilon + 4T_p^\epsilon - 3\theta_f^\epsilon + 3\eta\theta_F^\epsilon \\ & +2\eta\gamma^2 (u_F^\epsilon \cdot u_f^\epsilon - |u_f^\epsilon|^2) + \gamma^2 (|u_F^\epsilon|^2 - |u_f^\epsilon|^2) \end{aligned} \right]. \quad (145) $$

*Proof.* The expressions (143) are direct consequences of the applications of (141) and (142) to the general forms (79) and (80) of $\Theta_F^\epsilon$ and $\Theta_f^\epsilon$, while (144) and (145) follow by using (141) and (142) in the general expressions (83) and (84) of $\Xi_F^\epsilon$ and $\Xi_f^\epsilon$. $\square$

**Proposition E.2.** For the Maxwellian collision kernel, $T_p^\epsilon$ solves

$$ \gamma \delta\nu^2(1+\eta)\bar{\tau}\rho_F^\epsilon(\partial_t T_p^\epsilon + u_F^\epsilon \cdot \nabla_x T_p^\epsilon) = \rho_F^\epsilon \rho_f^\epsilon (3\eta\theta_F^\epsilon + 3\theta_f^\epsilon - 4T_p^\epsilon + \gamma^2|u_F^\epsilon - u_f^\epsilon|^2), \quad (146) $$

where

$$ \bar{\tau} := \frac{\tau}{2\pi\|b\|_{L^1}}. $$

*Proof.* Thanks to (141), $\Upsilon_1(F^\epsilon, f^\epsilon)$ and $\Upsilon_2(F^\epsilon, f^\epsilon)$ take the form

$$ \Upsilon_1(F^\epsilon, f^\epsilon) = 2\pi\|b\|_{L^1} \rho_F^\epsilon \rho_f^\epsilon, \quad \Upsilon_2(F^\epsilon, f^\epsilon) = 2\pi\|b\|_{L^1} \rho_F^\epsilon \rho_f^\epsilon (3\eta\theta_F^\epsilon + 3\theta_f^\epsilon + \gamma^2|u_F^\epsilon - u_f^\epsilon|^2). \quad (147) $$

By substituting these expressions into (87), we arrive at the desired conclusion. $\square$

We now take the limit as $\epsilon \to 0$, under the assumptions outlined in Section 6.

**Proposition E.3.** In the Maxwellian collision kernel case, we have, under Hypothesis (97),

$$ \begin{aligned} \Theta_F &= \lim_{\epsilon \to 0} \frac{1}{\delta\nu^2\sqrt{\eta}} \Theta_F^\epsilon = 2\pi\|b\|_{1} \rho_F \rho_F(u_F - u_F), \\ \Theta_f &= \lim_{\epsilon \to 0} \frac{1}{\delta\nu^2} \Theta_f^\epsilon = 2\pi\|b\|_{1} \rho_F \rho_F(u_F - u_F). \end{aligned} \quad (148) $$

Moreover, in the limit,

$$ \theta_F = \theta_f = \frac{4}{3}T_p. $$

*Proof.* First, the system of equations (148) is derived directly from (143). Secondly equation (146) writes

$$ \bar{\tau}(\partial_t T_p^\epsilon + u_F^\epsilon \cdot \nabla_x T_p^\epsilon) = \frac{1}{(1+\eta)} \rho_f^\epsilon \left[ -\frac{1}{\gamma^2} (4T_p^\epsilon - 3\theta_f^\epsilon - 3\eta\theta_F^\epsilon) + |u_f^\epsilon - u_F^\epsilon|^2 \right] \quad (149) $$

Since the left-hand side of the equation (149) is assumed to be bounded, we deduce that

$$ 4T_p^\epsilon - 3\theta_f^\epsilon - 3\eta\theta_F^\epsilon = O(\gamma^2). \quad (150) $$

In particular, we also have $4T_p^\epsilon - 3\theta_f^\epsilon \to 0$, and thanks to the Proposition 6, we get

$$ \theta_f = \frac{4}{3}T_p. $$

Note that this relation can also be obtained directly thanks to (101) and (147). Finally, we rewrite (144) as

$$ \frac{1}{\gamma^2} \Xi_F^\epsilon = \frac{2\pi\|b\|_{L^1}}{(1+\eta)^2} \rho_f^\epsilon \rho_F^\epsilon \left[ \begin{aligned} & \frac{\eta}{\gamma^2} (4T_p^\epsilon - 3\theta_f^\epsilon - 3\eta\theta_F^\epsilon) + \eta (|u_f^\epsilon|^2 - |u_F^\epsilon|^2) \\ & +6\frac{\eta}{\gamma^2} (\theta_f^\epsilon - \theta_F^\epsilon) + 2(u_F^\epsilon \cdot u_f^\epsilon - |u_f^\epsilon|^2) \end{aligned} \right], $$
where we have added and subtracted $3\eta\theta_f^\epsilon$ inside the square brackets. From the boundedness of the derivative of the macroscopic quantities in Equation (81) and Hypothesis (97), the term $\Xi_F^\epsilon/\gamma^2$ is bounded. Then, by using (150) and by replacing $\eta$ by $\gamma^s$, we get

$$ \theta_f^\epsilon - \theta_F^\epsilon = O(\gamma^{2-s}). $$

Hypothesis (99) allows us to conclude the equality $\theta_F = \theta_f$ in the limit. □

**Remark 14.** We can express

$$ \Xi_F = \lim_{\epsilon \to 0} \frac{1}{\gamma \delta \nu^2} \Xi_F^\epsilon = 2\pi \|b\|_{L^1} \rho_f \rho_F \left[ K_2 + 2(u_F \cdot u_f - |u_F|^2) \right] $$

and

$$ \Xi_f = \lim_{\epsilon \to 0} \frac{1}{\gamma \delta \nu^2} \Xi_f^\epsilon = 2\pi \|b\|_{L^1} \rho_f \rho_F \left[ K_1 - K_2 + (|u_F|^2 - |u_f|^2) \right]. $$

where

$$ K_1 := \lim_{\gamma \to 0} \frac{1}{\gamma^2} (4T_p^\epsilon - 3\theta_f^\epsilon - 3\gamma^s\theta_F^\epsilon) \quad \text{and} \quad K_2 := \lim_{\gamma \to 0} 6\gamma^{s-2} (\theta_f^\epsilon - \theta_F^\epsilon), $$

and the equation satisfied by $T_p$ is

$$ \bar{\tau} \left( \partial_t T_p + u_F \cdot \nabla_x T_p^\epsilon \right) = 2\rho_f \left( |u_f - u_F|^2 - K_1 \right). $$

We see that the system is not completely closed, and that two degrees of freedom remains.

**Acknowledgements.** This project has received funding from: The European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No 101034255, the COST Action: CA18232, via the Short Term Scientific Mission E-COST-GRANT-CA18232-a2e82880, the OeAD organization and the Austrian Ministry of Science, Education and Research through the Marietta Blau-grant. AM acknowledges partial support from the Austrian Science Fund (FWF), grant DOI 10.55776/P33010 and 10.55776/F65. This work has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme, ERC Advanced Grant NEUROMORPH, no. 101018153. The support of INDAM-GNFM group is acknowledged.

**Conflicts of interest.** Authors state that there is no conflict of interest.

## References

[1] B. Anwasia, M. Bisi, F. Salvarani, and A. J. Soares. On the Maxwell-Stefan diffusion limit for a reactive mixture of polyatomic gases in non-isothermal setting. *Kinet. Relat. Models*, 13(1):63–95, 2020.

[2] B. Anwasia, P. Gonçalves, and A. J. Soares. On the formal derivation of the reactive Maxwell-Stefan equations from the kinetic theory. *EPL (Europhysics Letters)*, 129(4):40005, 2020.

[3] B. Anwasia, P. Gonçalves, and A. J. Soares. From the simple reacting sphere kinetic model to the reaction-diffusion system of Maxwell-Stefan type. *Commun. Math. Sci.*, 17(2):507–538, 2019.

[4] A. Bondesan and M. Briant. Stability of the Maxwell-Stefan system in the diffusion asymptotics of the Boltzmann multi-species equation. *Comm. Math. Phys.*, 382(1):381–440, 2021.

[5] D. Bothe. On the Maxwell-Stefan approach to multicomponent diffusion. In *Parabolic problems*, volume 80 of *Progr. Nonlinear Differential Equations Appl.*, pages 81–93. Birkhäuser/Springer Basel AG, Basel, 2011.
[6] L. Boudin, B. Grec, and V. Pavan. The Maxwell-Stefan diffusion limit for a kinetic model of mixtures with general cross sections. *Nonlinear Anal.*, 159:40–61, 2017.

[7] L. Boudin, B. Grec, and F. Salvarani. A mathematical and numerical analysis of the Maxwell-Stefan diffusion equations. *Discrete Contin. Dyn. Syst. Ser. B*, 17(5):1427–1440, 2012.

[8] L. Boudin, B. Grec, and F. Salvarani. The Maxwell-Stefan diffusion limit for a kinetic model of mixtures. *Acta Appl. Math.*, 136:79–90, 2015.

[9] C. Cercignani. *The Boltzmann equation and its applications*, volume 67 of *Applied Mathematical Sciences*. Springer-Verlag, New York, 1988.

[10] C. Cercignani. *Multicomponent flow modeling*. Birkhäuser, 1999.

[11] C. Cercignani. *Rarefied gas dynamics: from basic concepts to actual calculations*, volume 21. Cambridge university press, 2000.

[12] F. Charles. Kinetic Modelling and Numerical Simulations using Particle Methods for the Transport of Dust in a Rarefied Gas. In *Proceeding of the 26 th symposium Rarefied gas dynamics*, volume 1084, pages 409–414. American Institute of Physics, 2009.

[13] F. Charles. *Modélisation mathématique et étude numérique d'un aérosol dans un gaz raréfié. Application à la simulation du transport de particules de poussière en cas d'accident de perte de vide dans ITER*. PhD thesis, École Normale Supérieure de Cachan, 2009.

[14] F. Charles. *Développement de modèles et analyse de méthodes numériques pour des problèmes cinétiques issus de la physique des plasmas et des écoulements gaz-particules*. PhD thesis, Sorbonne Université, 2021.

[15] F. Charles, S. Dellacherie, and J. Segré. Kinetic Modelling of the Transport of Dust Particles in a Rarefied Atmosphere. *Mathematical Models and Methods in Applied Sciences*, 22(40), 2012.

[16] F. Charles and L. Desvillettes. From collisional kinetic models to sprays: internal energy exchanges. *Communications in Mathematical Sciences*, 23(1):173–193, 2024.

[17] E. Cunningham. On the velocity of steady fall of spherical particles through fluid medium. *Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character*, 83(563):357–365, 1910.

[18] P. Degond and B. Lucquin-Desreux. The asymptotics of collision operators for two species of particles of disparate masses. *Mathematical Models and Methods in Applied Sciences*, 6(03):405–436, 1996.

[19] L. Desvillettes, F. Golse, and V. Ricci. A formal passage from a system of boltzmann equations for mixtures towards a vlasov-euler system of compressible fluids. *Acta Mathematicae Applicatae Sinica, English Series*, 35:158–173, 2019.

[20] A. Ern and V. Giovangigli. *Multicomponent transport algorithms*, volume 24 of *Lecture Notes in Physics. New Series M: Monographs*. Springer-Verlag, Berlin, 1994.

[21] A. Fick. On liquid diffusion. *Phil. Mag.*, 10(63):30–39, 1855.

[22] A. Fick. Ueber Diffusion. *Poggendorff's Annel Physik*, 94:59–86, 1855.

[23] V. Giovangigli. *Multicomponent flow modeling*. Modeling and Simulation in Science, Engineering and Technology. Birkhäuser Boston, Inc., Boston, MA, 1999.

[24] H. Hutridurga and F. Salvarani. On the maxwell–stefan diffusion limit for a mixture of monatomic gases. *Mathematical Methods in the Applied Sciences*, 40(3):803–813, 2017.
[25] H. Hutridurga and F. Salvarani. On the Maxwell-Stefan diffusion limit for a mixture of monatomic gases. *Math. Methods Appl. Sci.*, 40(3):803-813, 2017.

[26] R. Krishna and J. A. Wesselingh. The Maxwell-Stefan approach to mass transfer. *Chem. Eng. Sci.*, 52(6):861-911, 1997.

[27] J. C. Maxwell. On the dynamical theory of gases. *Phil. Trans. R. Soc.*, 157:49-88, 1866.

[28] B. Otterman and A. S. Levine. Analysis of gas-solid particle flows in shock tubes. *AIAA Journal*, 12(5):579-580, 1974.

[29] J. Stefan. Ueber das Gleichgewicht und die Bewegung insbesondere die Diffusion von Gasgemengen. *Akad. Wiss. Wien*, 63:63-124, 1871.

[30] R. Taylor and R. Krishna. *Multicomponent mass transfer*, volume 2. John Wiley & Sons, 1993.

[31] G. L. Vignoles, P. Charrier, C. Preux, and B. Dubroca. Rarefied pure gas transport in non-isothermal porous media: effective transport properties from homogenization of the kinetic equation. *Transport in Porous Media*, 73(2):211-232, 2008.

[32] C. Villani. A review of mathematical topics in collisional kinetic theory. *Handbook of mathematical fluid dynamics*, 1(71-305):3-8, 2002.

[33] C. Villani. Mathematics of granular materials. *Journal of statistical physics*, 124(2-4):781-822, 2006.