# Discrete *H*-theorem for a finite volume discretization of a nonlinear kinetic system: application to hypocoercivity

Marianne Bessemoulin-Chatard, Tino Laidin, Thomas Rey

► To cite this version:

Marianne Bessemoulin-Chatard, Tino Laidin, Thomas Rey. Discrete *H*-theorem for a finite volume discretization of a nonlinear kinetic system: application to hypocoercivity. 2025. hal-05369620

HAL Id: hal-05369620

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

Preprint submitted on 17 Nov 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# DISCRETE H-THEOREM FOR A FINITE VOLUME DISCRETIZATION OF A NONLINEAR KINETIC SYSTEM: APPLICATION TO HYPOCOERCIVITY

MARIANNE BESSEMOULIN-CHATARD, TINO LAIDIN, AND THOMAS REY

ABSTRACT. In this article, we study the long-time behavior of a finite-volume discretization for a nonlinear kinetic reaction model involving two interacting species. Building upon the seminal work of [Favre, Pirner, Schmeiser, ARMA, 2023], we extend the discrete exponential convergence to equilibrium result established in [Bessemoulin-Chatard, Laidin, Rey, IMAJNA, 2025], which was obtained in a perturbative framework using weighted $L^2$ estimates. The analysis applies to a broader class of exponentially decaying initial data, without requiring proximity to equilibrium, by exploiting the properties of the Boltzmann entropy. The proof relies on the propagation of the initial $L^\infty$ bounds, derived from monotonicity properties of the scheme, allowing controlled linearizations within the nonlinear entropy estimates. Moreover, we show that the time-discrete dissipation inherent to the numerical scheme plays a crucial stabilizing role, providing control over the nonlinear terms.

2020 MATHEMATICS SUBJECT CLASSIFICATION: 82B40, 65M08, 65M12.

## 1. INTRODUCTION

This work investigates the long-time behavior of an implicit-in-time finite-volume discretization for a nonlinear kinetic reaction model describing two interacting species in one-dimensional position-velocity phase space. We focus on the following system introduced in [33]

$$ \partial_t f_1 + v \partial_x f_1 = \chi_1 - \rho_2 f_1, \qquad (1.1) $$

$$ \partial_t f_2 + v \partial_x f_2 = \chi_2 - \rho_1 f_2, \qquad (1.2) $$

where $f_1$ and $f_2$ denote the phase-space densities of two chemical reactants, A and B. These species are produced by the decomposition of a third substance C (whose density is assumed to be constant) with nonnegative velocity profiles $\chi_1$ and $\chi_2$, and can also recombine to form C, thus leaving the system. The functions $f_1$ and $f_2$ depend on time $t \ge 0$, position $x \in \mathbb{T}$ (the one-dimensional torus), and velocity $v \in \mathbb{R}$. The reaction probability depends on the spatial density of the reaction partner, defined by

$$ \rho_k(t, x) := \int_{\mathbb{R}} f_k(t, x, v) \, dv, \quad k = 1, 2. \qquad (1.3) $$

The system (1.1)-(1.2) is complemented with the initial condition

$$ f_1(0, x, v) = f_1^{\text{in}}(x, v) \ge 0, \quad f_2(0, x, v) = f_2^{\text{in}}(x, v) \ge 0. \qquad (1.4) $$

Since the chemical reaction is assumed to be reversible, we impose the compatibility condition $\int_{\mathbb{R}} (\chi_1 - \chi_2) \, dv = 0$. Moreover, we make the following standard assumptions on the

<footnote>(Marianne Bessemoulin-Chatard) NANTES UNIVERSITÉ, CNRS, LABORATOIRE DE MATHÉMATIQUES JEAN LERAY, LMJL, UMR 6629, F-44000 NANTES, FRANCE</footnote>

<footnote>(Tino Laidin) UNIV BREST, CNRS UMR 6205, LABORATOIRE DE MATHÉMATIQUES DE BRETAGNE AT-LANTIQUE, F-29200 BREST, FRANCE</footnote>

<footnote>(Thomas Rey) UNIVERSITÉ CÔTE D'AZUR, CNRS, LJAD, PARC VALROSE, F-06108 NICE, FRANCE
E-mail addresses: marianne.bessemoulin@univ-nantes.fr, tino.laidin@univ-brest.fr, thomas.rey@univ-cotedazur.fr.</footnote>

<footnote>Key words and phrases. Kinetic equations, hypocoercivity, maximum principle, finite volume methods, large time behavior.</footnote>
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Comparison between the exponential decay of the relative Boltzmann entropy $H$ and of the weighted $L^2$ distance to equilibrium for the numerical scheme studied in [6].</figcaption></figure>

velocity profiles $\chi_k$, which are given positive functions of $v$. For $k = 1, 2$, we require

$$ \int_{\mathbb{R}} \chi_k \, dv = 1, \quad \int_{\mathbb{R}} v \chi_k \, dv = 0, \\ D_k := \int_{\mathbb{R}} v^2 \chi_k \, dv < \infty, \quad Q_k := \int_{\mathbb{R}} v^4 \chi_k \, dv < \infty. \tag{1.5} $$

In [25], additional linear relaxation terms of the form $\rho_k\chi_k - f_k$, $k = 1, 2$ were introduced to model thermalization effects relevant in semiconductor physics. While such relaxation terms could also be included within our analysis, the present work focuses on the nonlinear reaction terms in (1.1)–(1.2), which are sufficient to establish exponential convergence toward equilibrium.

Under these assumptions, the model satisfies a natural conservation property reflecting the reversible nature of the chemical reaction:

$$ \frac{d}{dt} \int_{\mathbb{T} \times \mathbb{R}} (f_1 - f_2) \, dv \, dx = 0. $$

This motivates the introduction of the unique constant $\rho^\infty > 0$ such that

$$ \int_{\mathbb{T} \times \mathbb{R}} (f_1^{\text{in}} - f_2^{\text{in}}) \, dv \, dx = |\mathbb{T}| \left( \rho^{\infty} - \frac{1}{\rho^{\infty}} \right), \tag{1.6} $$

and of the corresponding equilibrium state $\mathbf{F}^\infty = (f_1^\infty, f_2^\infty)^\top$, depending only on velocity, defined by

$$ f_1^\infty(x, v) := \rho_1^\infty \chi_1(v), \quad f_2^\infty(x, v) := \rho_2^\infty \chi_2(v), \tag{1.7} $$

with $\rho_1^\infty = \rho^\infty$ and $\rho_2^\infty = 1/\rho^\infty$.

The exponential convergence to equilibrium for the solution of (1.1)–(1.2) was first established in [33]. This result relies on a perturbative approach around equilibrium: a hypocoercivity result for the linearized system is first obtained in a weighted $L^2$ framework, following the method of [20], and then extended to the nonlinear system through the propagation of $L^\infty$ bounds. More recently, the authors of [25] removed the small-perturbation assumption by exploiting the nonlinear Boltzmann entropy of the system, allowing the treatment of more general initial data that satisfy weaker conditions relative to the equilibrium.

Motivated by these developments, we observed that the Boltzmann entropy $H(\mathbf{F})$ of the solution to the numerical scheme studied in [6]—for which exponential convergence to equilibrium was proven in a perturbative framework—also exhibits exponential decay in time (see Figure 1). The goal of the present work is to provide a rigorous proof of this decay by adapting the approach of [25] to the discrete setting. This extends our previous results to initial data that are not necessarily close to equilibrium, requiring only suitable $L^\infty$ bounds relative to the equilibrium profiles.

The long-time dynamics of (1.1)–(1.2) results from the interplay between the transport operator $v \partial_x$ and the reaction term, which induces a mixing effect in phase space. Under spatial periodic boundary conditions, this mixing drives the system toward the position-independent equilibrium $\mathbf{F}^\infty$. The study of the long-time asymptotics in kinetic equations is
a classical topic, dating back to the seminal work [32] on hypoellipticity for linear operators. Early results on the asymptotic behavior of kinetic equations were obtained in [3]. More recently, for a broad class of collision operators [31, 30, 20, 10], it has been shown that, in a suitable norm, solutions converge exponentially fast toward equilibrium: there exist constants $\kappa > 0$ and $C > 1$ such that

$$ \|\mathbf{F}(t) - \mathbf{F}^{\infty}\|_{\mathcal{X}} \le C \|\mathbf{F}^{\text{in}} - \mathbf{F}^{\infty}\|_{\mathcal{X}} e^{-\kappa t}, \quad (1.8) $$

for an appropriate Hilbert space $\mathcal{X}$. This property is known as hypocoercivity, following [36], whereas the special case $C = 1$ corresponds to the classical notion of coercivity.

Building on the general framework of hypocoercivity, the seminal work [20] introduced a robust and versatile method to establish exponential convergence to equilibrium. The core idea is to construct a modified entropy functional, equivalent to a weighted $L^2$-norm, for which exponential decay can be proven. This approach has since been successfully applied in a variety of kinetic contexts, including chemotaxis models [14], the fractional Fokker-Planck equation [9], the Vlasov-Fokker-Planck equation [11, 12], and the Vlasov-Poisson-Fokker-Planck system [1]. Extension to control problems in plasma theory, using refined estimates on the trajectories of the solutions, was also considered in [18].

More recently, in [25], the authors adapted this strategy to the nonlinear framework. Indeed, instead of relying on the evolution of a weighted $L^2$-norm, they construct a modified entropy functional directly from the nonlinear Boltzmann entropy, which can be shown to decay exponentially. Under an $L^\infty$ assumption on the initial data relative to the equilibrium profiles, the logarithmic terms appearing in the Boltzmann entropy can be carefully linearized, thereby quantifying the distance to equilibrium in the standard weighted $L^2$-norm.

The entropy-based and hypocoercive frameworks developed at the continuous level have also inspired a broad range of studies on their discrete counterparts, aiming to understand how numerical schemes reproduce long-time asymptotics. Over the past decades, the long-time behavior of discretizations of macroscopic models has been extensively studied using energy-based techniques. These methods rely on preserving the entropic structure of the continuous equations at the discrete level [28, 16, 13, 4]. Among them, the $\phi$-entropy method has proven to be a particularly robust tool for analyzing asymptotic behavior [26, 17]. The use of $\phi$-entropies was later extended to linear kinetic equations, such as the Fokker-Planck equation [19] and the linear relaxation model [24], demonstrating the effectiveness of nonlinear techniques in deriving quantitative convergence rates toward equilibrium.

These ideas naturally extend to the discretization of kinetic equations, allowing numerical schemes to accurately capture long-time asymptotics. Early developments in this direction focused on the Kolmogorov equation, where collisions are modelled by a Laplacian in velocity. A finite-difference scheme was analyzed in [34], while finite-element approaches were investigated in [27, 21]. These discrete hypocoercivity results led to the design of an efficient stabilization strategy for finite-element methods, as developed in [22]. The original $H^1$-hypocoercivity framework “à la Villani” [36] was first adapted to the discrete setting through finite-difference discretizations of the kinetic Fokker-Planck equation [23]. Subsequently, the robust weighted $L^2$ approach of [20] was extended to various discrete schemes for linear kinetic models. A discrete hypocoercivity property was first established for finite-volume discretizations of the Fokker-Planck and linear BGK equations in [5], and later generalized to the nonlinear model (1.1)-(1.2) within a perturbative regime in [6]. Long-time asymptotics for a finite-volume scheme of the fractional Fokker-Planck equation were also studied in [2]. Finally, this methodology has been applied to collisional Vlasov equations—both with external potentials and for the linearized Vlasov-Poisson system—via spectral decompositions based on Hermite polynomials [8, 7]. We also mention the recent work [29], which employs a fully spectral phase-space decomposition with a linearized BGK operator preserving the first three moments in an unbounded domain with a confining potential. This approach follows the analysis of [15], where convergence toward specific macroscopic modes (both stationary and time-periodic) was rigorously established.

**Outline of the paper.** Section 2 recall the main arguments of the continuous analysis, reformulated in a form suitable for the discrete framework. Section 3 introduces the numerical scheme, together with its discrete maximum principle and an existence result.
The maximum principle follows from a monotonicity property of the chosen upwind numerical fluxes. Finally, Section 4 establishes the exponential convergence of the discrete solution toward equilibrium by adapting the continuous proof. In this section, we focus on the specific challenges introduced by the discretization—such as the treatment of additional terms, the lack of a discrete chain rule, and the impact of time discretization—while the algebraic computations from the continuous setting carry over directly and require only a brief discussion.

## 2. THE CONTINUOUS SETTING

This works builds on and extends our previous contribution [6]. In particular, it relies on the recent work [25], adopting a fully nonlinear approach rather than relying on a linearization of the model. For clarity and readability, we present here a condensed version of the proof from [25], reformulated in a way that facilitates its adaptation to the discrete setting.

First of all, let us recall the global existence result and maximum principle estimates stated in [25, Theorem 1].

**Proposition 1** (Global existence and maximum principle). Under assumptions (1.5) on $\chi_k$, we further assume that the initial data $\mathbf{F}^{in}$ satisfies for all $x \in \mathbb{T}$, $v \in \mathbb{R}$,

$$ \rho_m \chi_1(v) \le f_1^{in}(x, v) \le \rho_M \chi_1(v), \qquad (2.1) $$

$$ \rho_M^{-1} \chi_2(v) \le f_2^{in}(x, v) \le \rho_m^{-1} \chi_2(v), \qquad (2.2) $$

with constants $0 < \rho_m \le \rho_M < \infty$.

Then there exists a unique global solution $\mathbf{F}$ of (1.1)–(1.2) such that for all $x \in \mathbb{T}$, $v \in \mathbb{R}$, $t \ge 0$,

$$ \rho_m \chi_1(v) \le f_1(t, x, v) \le \rho_M \chi_1(v), \qquad (2.3) $$

$$ \rho_M^{-1} \chi_2(v) \le f_2(t, x, v) \le \rho_m^{-1} \chi_2(v). \qquad (2.4) $$

2.1. Notations and functional setting. We now introduce the notations and functional framework that will be used throughout the analysis, providing the tools needed to handle both the continuous and discrete settings.

The decay towards equilibrium will be estimated quantitatively in the following weighted $L^2$ space

$$ \mathcal{X} := L^2(\mathbb{T} \times \mathbb{R}, dx \, dv/f_1^{\infty}) \times L^2(\mathbb{T} \times \mathbb{R}, dx \, dv/f_2^{\infty}), $$

endowed with the scalar product

$$ \langle \mathbf{F}, \mathbf{G} \rangle = \int_{\mathbb{T} \times \mathbb{R}} \left( \frac{f_1 g_1}{f_1^{\infty}} + \frac{f_2 g_2}{f_2^{\infty}} \right) dx \, dv, \quad \forall \mathbf{F} = (f_1, f_2)^{\top}, \mathbf{G} = (g_1, g_2)^{\top} \in \mathcal{X}. $$

The corresponding norm is denoted by $\| \cdot \|$.

We also introduce the space

$$ L^2(\mathbb{T}, dx/\rho_1^{\infty}) \times L^2(\mathbb{T}, dx/\rho_2^{\infty}), $$

equipped with the scalar product $\langle \cdot, \cdot \rangle_x$ and the associated norm $\| \cdot \|_x$.

We define the orthogonal projection onto the velocity equilibrium by

$$ \Pi \mathbf{F} := \begin{pmatrix} \rho_1 \chi_1 \\ \rho_2 \chi_2 \end{pmatrix}, $$

where $\rho_k$ is defined by (1.3).

To study the convergence of $\mathbf{F}(t)$ toward $\mathbf{F}^{\infty}$ as $t$ goes to $\infty$, we introduce the deviation $\tilde{\mathbf{F}}(t) := \mathbf{F}(t) - \mathbf{F}^{\infty}$ and rewrite the system (1.1)–(1.2) as

$$ \partial_t \tilde{f}_1 + v \partial_x \tilde{f}_1 = \chi_1 - \rho_2 f_1, \qquad (2.5) $$

$$ \partial_t \tilde{f}_2 + v \partial_x \tilde{f}_2 = \chi_2 - \rho_1 f_2. \qquad (2.6) $$
## 2.2. Moment dynamics and estimates.
To analyze the long-time behavior of the system, it is useful to consider its macroscopic quantities [20]. These moments capture essential features of the dynamics and play a key role in establishing coercivity properties. We now define these moments and provide corresponding estimates with respect to the norm $|| \cdot ||$. The moments of $\mathbf{F} = (f_1, f_2)^\top$ are defined as

$$ \rho_{\mathbf{F}} := (\rho_1, \rho_2)^{\top}, \quad \text{with } \rho_k = \int_{\mathbb{R}} f_k \, dv, \qquad (2.7) $$

$$ \mathbf{J}_{\mathbf{F}} := (J_1, J_2)^{\top}, \quad \text{with } J_k = \int_{\mathbb{R}} v f_k \, dv, \qquad (2.8) $$

$$ \mathbf{S}_{\mathbf{F}} := (S_1, S_2)^{\top}, \quad \text{with } S_k = \int_{\mathbb{R}} (v^2 - D_k) f_k \, dv, \qquad (2.9) $$

$D_k$ being the second-order moment of $\chi_k$, defined in (1.5). Let us remark that with these notations,

$$ \tilde{\rho} := \rho_{\tilde{\mathbf{F}}} = \rho_{\mathbf{F}} - \rho^{\infty}, \quad \tilde{\mathbf{J}} := \mathbf{J}_{\tilde{\mathbf{F}}} = \mathbf{J}_{\mathbf{F}}, \quad \tilde{\mathbf{S}} := \mathbf{S}_{\tilde{\mathbf{F}}} = \mathbf{S}_{\mathbf{F}}. $$

By multiplying (2.5) and (2.6) by $(1, v)^{\top}$ and integrating with respect to $v$, we get the following moment equations for $\tilde{\mathbf{F}}$:

$$ \partial_t \tilde{\rho} + \partial_x \tilde{\mathbf{J}} = (1 - \rho_1 \rho_2) \begin{pmatrix} 1 \\ 1 \end{pmatrix}, \qquad (2.10) $$

and

$$ \partial_t \tilde{J}_1 + \partial_x \tilde{S}_1 + D_1 \partial_x \tilde{\rho}_1 = -\rho_2 J_1, \qquad (2.11) $$

$$ \partial_t \tilde{J}_2 + \partial_x \tilde{S}_2 + D_2 \partial_x \tilde{\rho}_2 = -\rho_1 J_2. \qquad (2.12) $$

The following lemma provides estimates, which are essential for controlling the macroscopic contributions in the subsequent hypocoercivity analysis.

**Lemma 1 (Moments estimates).** Under assumptions (1.5), the moments satisfy the following estimates

$$ \|\tilde{\rho}\|_x = \|\Pi \tilde{\mathbf{F}}\|, \qquad (2.13) $$

$$ \|\tilde{\mathbf{J}}\|_x = \|\mathbf{J}_{\mathbf{F}}\|_x \le C_1 \|(I - \Pi)\tilde{\mathbf{F}}\|, \qquad (2.14) $$

$$ \|\tilde{\mathbf{S}}\|_x = \|\mathbf{S}_{\mathbf{F}}\|_x \le C_2 \|(I - \Pi)\tilde{\mathbf{F}}\|, \qquad (2.15) $$

where the constants are given by

$$ C_1 = \sqrt{\max(D_1, D_2)}, \qquad (2.16) $$

$$ C_2 = \sqrt{\max(Q_1 - D_1^2, Q_2 - D_2^2)}. \qquad (2.17) $$

*Proof.* We start by proving (2.13). Using the definition of $\Pi$ and the first property in (1.5), one has

$$ \|\Pi \tilde{\mathbf{F}}\|^2 = \int_{\mathbb{T} \times \mathbb{R}} \left[ \frac{(\tilde{\rho}_1 \chi_1)^2}{\chi_1 \rho_1^\infty} + \frac{(\tilde{\rho}_2 \chi_2)^2}{\chi_2 \rho_2^\infty} \right] dx \, dv = \int_{\mathbb{T} \times \mathbb{R}} \left[ \frac{\tilde{\rho}_1^2 \chi_1}{\rho_1^\infty} + \frac{\tilde{\rho}_2^2 \chi_2}{\rho_2^\infty} \right] dx \, dv = \|\tilde{\rho}\|_x^2. $$

Now for (2.14), we have by definition

$$ \|\mathbf{J}_{\mathbf{F}}\|_x^2 = \int_{\mathbb{T}} \sum_{k=1,2} \frac{1}{\rho_k^\infty} \left( \int_{\mathbb{R}} v f_k \, dv \right)^2 \, dx. $$

Using the second property in (1.5), we get

$$ \left( \int_{\mathbb{R}} v f_k \, dv \right)^2 = \left( \int_{\mathbb{R}} v (f_k - \rho_k \chi_k) \frac{\sqrt{\chi_k}}{\sqrt{\chi_k}} \, dv \right)^2 \\ \le \left( \int_{\mathbb{R}} \frac{(f_k - \rho_k \chi_k)^2}{\chi_k} \, dv \right) \left( \int_{\mathbb{R}} v^2 \chi_k \, dv \right), $$

which yields the result using the definition of $D_k$ given in (1.5).
Finally, for (2.15), we have by definition

$$ \|\mathbf{S}_F\|_x^2 = \sum_{k=1,2} \int_{\mathbb{T}} \frac{1}{\rho_k^\infty} \left( \int_{\mathbb{R}} (v^2 - D_k) f_k \, dv \right)^2 dx. $$

Now, using the definition of $D_k$, we have $\int_{\mathbb{R}} (v^2 - D_k) \rho_k \chi_k \, dv = 0$, which gives

$$ \left( \int_{\mathbb{R}} (v^2 - D_k) f_k \, dv \right)^2 = \left( \int_{\mathbb{R}} (v^2 - D_k) (f_k - \rho_k \chi_k) \frac{\sqrt{\chi_k}}{\sqrt{\chi_k}} \, dv \right)^2 \\ \leq \left( \int_{\mathbb{R}} (v^2 - D_k)^2 \chi_k \, dv \right) \left( \int_{\mathbb{R}} \frac{(f_k - \rho_k \chi_k)^2}{\chi_k} \, dv \right), $$

yielding the result using properties (1.5) of $\chi_k$.

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

**2.3. Boltzmann entropy dissipation.** We define the relative Boltzmann entropy of $\mathbf{F} = (f_1, f_2)^\top$ with respect to the equilibrium $\mathbf{F}^\infty = (f_1^\infty, f_2^\infty)^\top$ as

$$ H(\mathbf{F}) := \sum_{k=1,2} \int_{\mathbb{T} \times \mathbb{R}} \left( f_k \left[ \log \left( \frac{f_k}{f_k^\infty} \right) - 1 \right] + f_k^\infty \right) \, dv \, dx. \quad (2.18) $$

The first step of the analysis is to establish the following microscopic coercivity property, which provides dissipation a priori solely in the velocity variable.

**Lemma 2 (Microscopic coercivity).** Let (1.5) hold and let $\mathbf{F} = (f_1, f_2)^\top$ be the solution to the system (1.1)-(1.2) with initial data $\mathbf{F}^{\text{in}} = (f_1^{\text{in}}, f_2^{\text{in}})^\top \in \mathcal{X}$. Then for every $t \ge 0$,

$$ \frac{d}{dt} H(\mathbf{F}) = -\mathcal{D}(\mathbf{F}) \le 0, \quad (2.19) $$

where the dissipation is defined by

$$ \mathcal{D}(\mathbf{F}) := \int_{\mathbb{T} \times \mathbb{R}^2} (f_1 f_2' - \chi_1 \chi_2') \log \left( \frac{f_1 f_2'}{\chi_1 \chi_2'} \right) \, dv \, dv' \, dx. \quad (2.20) $$

*Proof.* Since $\mathbf{F} = (f_1, f_2)^\top$ satisfies (1.1)-(1.2), we have

$$ \begin{aligned} \frac{d}{dt} H(\mathbf{F}) &= \sum_{k=1,2} \int_{\mathbb{T} \times \mathbb{R}} \partial_t f_k \log (f_k/f_k^\infty) \, dv \, dx \\ &= \int_{\mathbb{T} \times \mathbb{R}} (-v \partial_x f_1 + \chi_1 - \rho_2 f_1) \log (f_1/f_1^\infty) \, dv \, dx \\ &\quad + \int_{\mathbb{T} \times \mathbb{R}} (-v \partial_x f_2 + \chi_2 - \rho_1 f_2) \log (f_2/f_2^\infty) \, dv \, dx. \end{aligned} $$

Performing an integration by parts with respect to $x$ and using the periodic boundary conditions, we remark that the transport terms vanish:

$$ \begin{aligned} - \int_{\mathbb{T} \times \mathbb{R}} v \partial_x f_k \log (f_k/f_k^\infty) \, dv \, dx &= \int_{\mathbb{T} \times \mathbb{R}} v f_k \partial_x (\log (f_k/f_k^\infty)) \, dv \, dx \\ &= \int_{\mathbb{T} \times \mathbb{R}} v \partial_x f_k \, dv \, dx = 0. \end{aligned} $$

It remains to treat the relaxation terms. To do so, we double the velocity variables using that $\int_{\mathbb{R}} \chi_k' \, dv' = 1$ and $\rho_k = \int_{\mathbb{R}} f_k' \, dv'$:

$$ \int_{\mathbb{T} \times \mathbb{R}} \left[ (\chi_1 - \rho_2 f_1) \log \left( \frac{f_1}{\rho_1^\infty \chi_1} \right) + (\chi_2 - \rho_1 f_2) \log \left( \frac{f_2}{\rho_2^\infty \chi_2} \right) \right] \, dv \, dx = \\ \int_{\mathbb{T} \times \mathbb{R}^2} \left[ (\chi_1 \chi_2' - f_2' f_1) \log \left( \frac{f_1}{\rho_1^\infty \chi_1} \right) + (\chi_2 \chi_1' - f_1' f_2) \log \left( \frac{f_2}{\rho_2^\infty \chi_2} \right) \right] \, dv' \, dv \, dx. $$

By reversing the integration variables $v$ and $v'$ in the second term, we get

$$ \frac{d}{dt} H(\mathbf{F}) = \int_{\mathbb{T} \times \mathbb{R}^2} (\chi_1 \chi_2' - f_1 f_2') \left[ \log \left( \frac{f_1 f_2'}{\chi_1 \chi_2'} \right) - \log(\rho_1^\infty \rho_2^\infty) \right] \, dv' \, dv \, dx, $$

which concludes the proof using that $\rho_1^\infty \rho_2^\infty = 1$.

<figure><img src="image_4.png" /></figure>
We also require the following technical result to establish exponential convergence toward equilibrium. It connects the dissipation to the nonlinear reaction terms in the moment equations (2.11)-(2.12), as well as to the dissipation in the direction of $I - \Pi$, which corresponds to dissipation in the velocity variable.

**Lemma 3.** Under assumptions (1.5) and (2.1)-(2.2), the dissipation term defined by (2.20) satisfies

$$ \mathcal{D}(\mathbf{F}) \geq C_3 \|1 - \rho_1 \rho_2\|_{L^2(\mathbb{T})}^2 + C_4 \|(I - \Pi)\tilde{\mathbf{F}}\|^2, \quad (2.21) $$

$$ \text{where } C_3 = \rho_m / \rho_M \text{ and } C_4 = (\rho_m / \rho_M)(\rho_1^\infty + \rho_2^\infty)^{-1} \min(\rho_m, \rho_M^{-1})^2. $$

*Proof.* Using Taylor expansion, we have

$$ \log \left( \frac{f_1 f'_2}{\chi_1 \chi'_2} \right) = \frac{1}{\widehat{f_1 f'_2}} (f_1 f'_2 - \chi_1 \chi'_2), $$

where $\widehat{f_1 f'_2} \in (\min(f_1 f'_2, \chi_1 \chi'_2), \max(f_1 f'_2, \chi_1 \chi'_2))$. In particular $\widehat{f_1 f'_2} \le (\rho_M / \rho_m) \chi_1 \chi'_2$ by assumptions (2.1)-(2.2), and it yields

$$ \mathcal{D}(\mathbf{F}) \geq \frac{\rho_m}{\rho_M} \int_{\mathbb{T} \times \mathbb{R}^2} \frac{(f_1 f'_2 - \chi_1 \chi'_2)^2}{\chi_1 \chi'_2} dv' \, dv \, dx =: \frac{\rho_m}{\rho_M} \mathcal{E}(\mathbf{F}). \quad (2.22) $$

Moreover, using that

$$ \rho_1^\infty + \rho_2^\infty = \frac{1}{\rho_1^\infty} + \frac{1}{\rho_2^\infty} $$

and $\int_{\mathbb{R}} \chi'_k dv' = 1$, we remark that

$$ \begin{aligned} (\rho_1^\infty + \rho_2^\infty) \mathcal{E}(\mathbf{F}) &= \int_{\mathbb{T} \times \mathbb{R}} \frac{1}{\rho_1^\infty} \int_{\mathbb{R}} \frac{(f_1 f'_2 - \chi_1 \chi'_2)^2}{\chi_1 \chi'_2} dv' \, dv \, dx + \int_{\mathbb{T} \times \mathbb{R}} \frac{1}{\rho_2^\infty} \int_{\mathbb{R}} \frac{(f'_1 f_2 - \chi'_1 \chi_2)^2}{\chi'_1 \chi_2} dv' \, dv \, dx \\ &= \sum_{\substack{k,l=1,2 \\ k \neq l}} \int_{\mathbb{T} \times \mathbb{R}} \frac{1}{\rho_k^\infty} \chi_k \left[ \left( \int_{\mathbb{R}} \frac{f_k f'_l - \chi_k \chi'_l}{\chi'_l} dv' \right)^{\frac{1}{2}} \left( \int_{\mathbb{R}} \chi'_l dv' \right)^{\frac{1}{2}} \right]^2 . \end{aligned} $$

Then, thanks to the Cauchy-Schwarz inequality, we get

$$ \begin{aligned} (\rho_1^\infty + \rho_2^\infty) \mathcal{E}(\mathbf{F}) &\geq \sum_{\substack{k,l=1,2 \\ k \neq l}} \int_{\mathbb{T} \times \mathbb{R}} \frac{1}{\rho_k^\infty \chi_k} \left[ \int_{\mathbb{R}} (f_k f'_l - \chi_k \chi'_l) dv' \right]^2 dv \, dx \\ &\geq \sum_{\substack{k,l=1,2 \\ k \neq l}} \int_{\mathbb{T} \times \mathbb{R}} \frac{1}{\rho_k^\infty \chi_k} (f_k \rho_l - \chi_k)^2 dv \, dx \\ &\geq \sum_{\substack{k,l=1,2 \\ k \neq l}} \int_{\mathbb{T} \times \mathbb{R}} \frac{1}{\rho_k^\infty \chi_k} [\rho_l(f_k - \rho_k \chi_k) - \chi_k(1 - \rho_k \rho_l)]^2 dv \, dx. \end{aligned} $$

Noticing that the cross term in the expansion of the square satisfies

$$ \int_{\mathbb{T} \times \mathbb{R}} \frac{1}{\rho_k^\infty \chi_k} \rho_l (f_k - \rho_k \chi_k) \chi_k (1 - \rho_k \rho_l) dv \, dx = 0, $$

we obtain

$$ \begin{aligned} (\rho_1^\infty + \rho_2^\infty) \mathcal{E}(\mathbf{F}) &\geq \sum_{\substack{k,l=1,2 \\ k \neq l}} \int_{\mathbb{T} \times \mathbb{R}} \frac{1}{\rho_k^\infty \chi_k} (\rho_l^2 (f_k - \rho_k \chi_k)^2 + \chi_k^2 (1 - \rho_k \rho_l)^2) \, dv \, dx \\ &\geq \sum_{\substack{k,l=1,2 \\ k \neq l}} \int_{\mathbb{T} \times \mathbb{R}} \frac{1}{f_k^\infty} \rho_l^2 (f_k - \rho_k \chi_k)^2 \, dv \, dx + \left( \frac{1}{\rho_1^\infty} + \frac{1}{\rho_2^\infty} \right) \int_{\mathbb{T}} (1 - \rho_1 \rho_2)^2 \, dv \\ &\geq \min(\rho_m, \rho_M^{-1})^2 \sum_{\substack{k,l=1,2 \\ k \neq l}} \int_{\mathbb{T} \times \mathbb{R}} \frac{1}{f_k^\infty} (f_k - \rho_k \chi_k)^2 \, dv \, dx + (\rho_1^\infty + \rho_2^\infty) \|1 - \rho_1 \rho_2\|_{L^2(\mathbb{T})}^2, \end{aligned} $$
where we have used $\rho_1^2 \ge \rho_m^2$ and $\rho_2^2 \ge \rho_M^{-2}$ in the last inequality. Finally, by definition of $I - \Pi$ and since $\mathbf{F}^\infty$ belongs to the nullspace of $I - \Pi$, we get

$$ (\rho_1^\infty + \rho_2^\infty)\mathcal{E}(\mathbf{F}) \ge \min(\rho_m, \rho_M^{-1})^2 \|(I - \Pi)\mathbf{F}\|^2 + (\rho_1^\infty + \rho_2^\infty)\|1 - \rho_1\rho_2\|_{L^2(\mathbb{T})}^2 \\ \ge \min(\rho_m, \rho_M^{-1})^2 \|(I - \Pi)\tilde{\mathbf{F}}\|^2 + (\rho_1^\infty + \rho_2^\infty)\|1 - \rho_1\rho_2\|_{L^2(\mathbb{T})}^2, $$

which concludes the proof by combining this estimate with (2.22). □

## 2.4. Nonlinear hypocoercivity.
With the microscopic and macroscopic estimates established in the previous subsections, we are now in a position to state the main hypocoercivity result given in [25]. This theorem quantifies the exponential convergence of the solution toward equilibrium, combining the velocity dissipation and the macroscopic coercivity into a unified nonlinear estimate.

**Theorem 1.** Under assumptions (1.5), there exist two constants $C \ge 1$ and $\kappa > 0$ such that for all initial data $\mathbf{F}^{in} = (f_1^{in}, f_2^{in})^\top \in \mathcal{X}$ fulfilling (2.1)-(2.2), the solution $\mathbf{F} = (f_1, f_2)^\top$ of (1.1)-(1.2) satisfies

$$ \|\mathbf{F}(t) - \mathbf{F}^\infty\| \le C \|\mathbf{F}^{in} - \mathbf{F}^\infty\| e^{-\kappa t} \quad \forall t \ge 0. \qquad (2.23) $$

In order to prove Theorem 1, we introduce a modified entropy functional, following the approach of [25, 20]:

$$ \Gamma(t) := H(\mathbf{F}) + \delta \langle \tilde{\mathbf{J}}, \partial_x \mathbf{\Phi} \rangle_x, \qquad (2.24) $$

where $\mathbf{\Phi} = (\phi_1, \phi_2)^\top$, with each $\phi_k$ solving the Poisson equation

$$ -\partial_{xx}\phi_k = \tilde{\rho}_k, \quad \int_{\mathbb{T}} \phi_k \, dx = 0, \quad k = 1, 2, \qquad (2.25) $$

and $\delta > 0$ is a small parameter to be fixed later. The additional term in (2.24) is designed to explicitly recover dissipation in the spatial variable, reflecting the mixing effect induced by the transport operator.

We begin with key estimates on $\mathbf{\Phi}$ that will be used to control the additional term in the modified entropy $\Gamma(t)$.

**Lemma 4.** The function $\mathbf{\Phi}$ satisfies for all $t \ge 0$

$$ \| \partial_x \mathbf{\Phi} \|_x \le C_5 \| \Pi \tilde{\mathbf{F}} \|, \qquad (2.26) $$

$$ \| \partial_t \partial_x \mathbf{\Phi} \|_x \le C_1 \|(I - \Pi)\tilde{\mathbf{F}}\| + C_6 \|1 - \rho_1 \rho_2\|_{L^2(\mathbb{T})}, \qquad (2.27) $$

where $C_5$ is the Poincaré constant of $\mathbb{T}$ and $C_6 = C_5 \sqrt{\rho_1^\infty + \rho_2^\infty}$.

**Proof.** The first estimate is obtained by using the Poisson equation (2.25), applying the Poincaré inequality and equality (2.13):

$$ \| \partial_x \mathbf{\Phi} \|_x^2 = \langle \mathbf{\Phi}, -\partial_{xx} \mathbf{\Phi} \rangle_x = \langle \mathbf{\Phi}, \tilde{\rho} \rangle_x \le \| \mathbf{\Phi} \|_x \| \tilde{\rho} \|_x \le C_5 \| \partial_x \mathbf{\Phi} \|_x \| \Pi \tilde{\mathbf{F}} \|. $$

For the second estimate, we differentiate the Poisson equation with respect to time and use the continuity equation (2.10) to get

$$ \| \partial_t \partial_x \mathbf{\Phi} \|_x^2 = \langle \partial_t \mathbf{\Phi}, -\partial_t \partial_{xx} \mathbf{\Phi} \rangle_x = \langle \partial_t \mathbf{\Phi}, \partial_t \tilde{\rho} \rangle_x \\ = \left\langle \partial_t \mathbf{\Phi}, -\partial_t \tilde{\mathbf{J}} + (1 - \rho_1 \rho_2) \begin{pmatrix} 1 \\ 1 \end{pmatrix} \right\rangle_x. $$

An integration by parts on the first term and the Cauchy-Schwarz inequality give

$$ \| \partial_t \partial_x \mathbf{\Phi} \|_x^2 = \langle \partial_t \partial_x \mathbf{\Phi}, \tilde{\mathbf{J}} \rangle_x + \langle \partial_t \mathbf{\Phi}, (1 - \rho_1 \rho_2) \begin{pmatrix} 1 \\ 1 \end{pmatrix} \rangle_x \\ \le \| \partial_t \partial_x \mathbf{\Phi} \|_x \| \tilde{\mathbf{J}} \|_x + \| \partial_t \mathbf{\Phi} \|_x \left\| \begin{pmatrix} 1 - \rho_1 \rho_2 \\ 1 - \rho_1 \rho_2 \end{pmatrix} \right\|_x. \qquad (2.28) $$

Using that $\rho_1^\infty \rho_2^\infty = 1$, we have

$$ \left\| \begin{pmatrix} 1 - \rho_1 \rho_2 \\ 1 - \rho_1 \rho_2 \end{pmatrix} \right\|_x^2 = \left( \frac{1}{\rho_1^\infty} + \frac{1}{\rho_2^\infty} \right) \int_{\mathbb{T}} (1 - \rho_1 \rho_2)^2 \, dx = (\rho_1^\infty + \rho_2^\infty) \|1 - \rho_1 \rho_2\|_{L^2(\mathbb{T})}^2. $$

Using this together with the Poincaré inequality on the second term of (2.28), and the estimate (2.14) on the first term of (2.28) finally gives the result. □
Using the $L^\infty$ bounds (2.3)-(2.4) on $\mathbf{F}$ together with the moment estimates, we show that, for $\delta > 0$ sufficiently small, the modified entropy $\Gamma(t)$ defines a norm equivalent to $\|\cdot\|$.

**Lemma 5 (Equivalent norm).** There is $\delta_1 > 0$ such that for all $\delta \in (0, \delta_1)$, there are positive constants $0 < c_\delta < C_\delta$ such that if $\mathbf{F} \in \mathcal{X}$ satisfies (2.3)-(2.4), one has

$$c_\delta \|\tilde{\mathbf{F}}(t)\|^2 \leq \Gamma(t) \leq C_\delta \|\tilde{\mathbf{F}}(t)\|^2.$$

*Proof.* We first treat the Boltzmann part of $\Gamma$. Applying a Taylor expansion to the function

$$\xi \mapsto \xi \left( \log \left( \frac{\xi}{f_k^\infty} \right) - 1 \right) + f_k^\infty,$$

there exists $\xi_k \in (\min(f_k, f_k^\infty), \max(f_k, f_k^\infty))$ such that

$$f_k \left[ \log \left( \frac{f_k}{f_k^\infty} \right) - 1 \right] + f_k^\infty = \frac{1}{2\xi_k} (f_k - f_k^\infty)^2.$$

Since $f_k^\infty = \chi_k \rho_k^\infty$ and given the bounds (2.3)-(2.4) on $f_k$, it holds

$$\begin{aligned} \frac{\rho_m}{\rho_1^\infty} f_1^\infty &\le \xi_1 \le \frac{\rho_M}{\rho_1^\infty} f_1^\infty, \\ \frac{1}{\rho_M \rho_2^\infty} f_2^\infty &\le \xi_2 \le \frac{1}{\rho_m \rho_2^\infty} f_2^\infty. \end{aligned}$$

We then have

$$\begin{aligned} H(\mathbf{F}) &= \int_{\mathbb{T} \times \mathbb{R}} \left[ \frac{(f_1 - f_1^\infty)^2}{2\xi_1} + \frac{(f_2 - f_2^\infty)^2}{2\xi_2} \right] dv \, dx \\ &\le \int_{\mathbb{T} \times \mathbb{R}} \left[ \frac{1}{2}(f_1 - f_1^\infty)^2 \frac{\rho_1^\infty}{\rho_m f_1^\infty} + \frac{1}{2}(f_2 - f_2^\infty)^2 \frac{\rho_M \rho_2^\infty}{f_2^\infty} \right] dv \, dx \\ &\le \frac{1}{2} \max \left( \frac{\rho_1^\infty}{\rho_m}, \rho_M \rho_2^\infty \right) \|\tilde{\mathbf{F}}\|^2 := C_H \|\tilde{\mathbf{F}}\|^2. \end{aligned}$$

Proceeding in the same way for the lower bound, we obtain

$$H(\mathbf{F}) \ge \frac{1}{2} \min \left( \frac{\rho_1^\infty}{\rho_M}, \rho_M \rho_2^\infty \right) \|\tilde{\mathbf{F}}\|^2 := c_H \|\tilde{\mathbf{F}}\|^2.$$

We now treat the second term of $\Gamma$:

$$\left| \left\langle \tilde{\mathbf{J}}, \partial_x \mathbf{\Phi} \right\rangle_x \right| \le \|\tilde{\mathbf{J}}\|_x \|\partial_x \mathbf{\Phi}\|_x \le C_1 \|(I - \Pi)\tilde{\mathbf{F}}\| C_5 \|\Pi \tilde{\mathbf{F}}\| \le C_1 C_5 \|\tilde{\mathbf{F}}\|^2.$$

This finally gives the expected result with $c_\delta = c_H - \delta C_1 C_5$ and $C_\delta = C_H + \delta C_1 C_5$. One concludes the proof by choosing $\delta > 0$ such that $c_\delta > 0$, that is $\delta < \delta_1 := C_H/(C_1 C_5)$. $\square$

Combining lemmas 3 and 5, one can then establish the proof of Theorem 1.

*Proof of Theorem 1.* Using the microscopic coercivity (2.19), we have

$$\frac{d}{dt} \Gamma(t) = -\mathcal{D}(\mathbf{F}) + \delta \frac{d}{dt} \left\langle \tilde{\mathbf{J}}, \partial_x \mathbf{\Phi} \right\rangle_x.$$

We start by studying the second term. Using the continuity equations (2.11)-(2.12), one has

$$\frac{d}{dt} \left\langle \tilde{\mathbf{J}}, \partial_x \mathbf{\Phi} \right\rangle_x = \left\langle \partial_t \tilde{\mathbf{J}}, \partial_x \mathbf{\Phi} \right\rangle_x + \left\langle \tilde{\mathbf{J}}, \partial_t \partial_x \mathbf{\Phi} \right\rangle_x = T_1 + T_2 + T_3 + T_4,$$

where

$$T_1 := - \left\langle \partial_x \tilde{\mathbf{S}}, \partial_x \mathbf{\Phi} \right\rangle_x,$$

$$T_2 := - \left\langle \left( \frac{D_1}{D_2} \partial_x \tilde{\rho}_1}{\partial_x \tilde{\rho}_2} \right), \partial_x \mathbf{\Phi} \right\rangle_x,$$

$$T_3 := - \left\langle \left( \frac{\rho_2 J_1}{\rho_1 J_2} \right), \partial_x \mathbf{\Phi} \right\rangle_x,$$

$$T_4 := \left\langle \tilde{\mathbf{J}}, \partial_t \partial_x \mathbf{\Phi} \right\rangle_x.$$

The term $T_1$ can be easily controlled with moment estimates (2.13) and (2.15):

$$|T_1| = \left| \left\langle \tilde{\mathbf{S}}, \partial_{xx} \mathbf{\Phi} \right\rangle_x \right| = \left| - \left\langle \tilde{\mathbf{S}}, \tilde{\rho} \right\rangle_x \right| \le \|\tilde{\mathbf{S}}\|_x \|\tilde{\rho}\|_x \le C_2 \|(I - \Pi)\tilde{\mathbf{F}}\| \|\Pi\tilde{\mathbf{F}}\|.$$
The term $T_2$ provides the macroscopic coercivity. Indeed, using an integration by parts and the Poisson equation (2.25), we get

$$ T_2 = - \sum_{k=1,2} \frac{D_k}{\rho_k^\infty} \langle \partial_x \tilde{\rho}_k, \partial_x \phi_k \rangle_{L^2(\mathbf{T})} = \sum_{k=1,2} \frac{D_k}{\rho_k^\infty} \langle \tilde{\rho}_k, \partial_{xx} \phi_k \rangle_{L^2(\mathbf{T})} \\ \leq - \min(D_1, D_2) \|\tilde{\rho}\|_x^2 = - \min(D_1, D_2) \|\Pi \tilde{\mathbf{F}}\|^2. $$

In the terms $T_3$ and $T_4$, nonlinear contributions appear. These will be controlled on the one hand by the Boltzmann entropy dissipation (microscopic coercivity), and on the other hand by the $L^\infty$ bounds on $\mathbf{F}$. Using the Cauchy-Schwarz inequality, we have

$$ |T_3| \leq \left\| \begin{pmatrix} \rho_2 J_1 \\ \rho_1 J_2 \end{pmatrix} \right\|_x \| \partial_x \mathbf{\Phi} \|_x. $$

From the $L^\infty$ bounds on $\mathbf{F}$, we deduce that

$$ \rho_m \leq \rho_1 \leq \rho_M, \quad \rho_M^{-1} \leq \rho_2 \leq \rho_m^{-1}, $$

which yields

$$ \begin{aligned} \left\| \begin{pmatrix} \rho_2 J_1 \\ \rho_1 J_2 \end{pmatrix} \right\|_x^2 &= \frac{1}{\rho_1^\infty} \int_{\mathbf{T}} (\rho_2 J_1)^2 dx + \frac{1}{\rho_2^\infty} \int_{\mathbf{T}} (\rho_1 J_2)^2 dx \\ &\leq \frac{\rho_m^{-2}}{\rho_1^\infty} \int_{\mathbf{T}} J_1^2 dx + \frac{\rho_M^2}{\rho_2^\infty} \int_{\mathbf{T}} J_2^2 dx \\ &\leq \max(\rho_m^{-1}, \rho_M)^2 \|\mathbf{J}\mathbf{F}\|_x^2 \\ &\leq \max(\rho_m^{-1}, \rho_M)^2 C_1^2 \|(I - \Pi)\tilde{\mathbf{F}}\|^2, \end{aligned} $$

and finally, using also (2.26),

$$ |T_3| \leq \max(\rho_m^{-1}, \rho_M) C_1 C_5 \|(I - \Pi)\tilde{\mathbf{F}}\| \|\Pi \tilde{\mathbf{F}}\|. $$

Let us now turn to the estimate of $T_4$. Thanks to the Cauchy-Schwarz inequality and using estimates (2.14) and (2.27), we write

$$ \begin{aligned} |T_4| &\leq \|\tilde{\mathbf{J}}\|_x \|\partial_t \partial_x \mathbf{\Phi}\|_x \\ &\leq C_1 \|(I - \Pi)\tilde{\mathbf{F}}\| \left[ C_1 \|(I - \Pi)\tilde{\mathbf{F}}\| + C_6 \|1 - \rho_1 \rho_2\|_{L^2(\mathbf{T})} \right] \\ &\leq C_1^2 \|(I - \Pi)\tilde{\mathbf{F}}\|^2 + C_1 C_6 \|(I - \Pi)\tilde{\mathbf{F}}\| \|1 - \rho_1 \rho_2\|_{L^2(\mathbf{T})}. \end{aligned} $$

Using the Young inequality on the second term of the right hand side, we get

$$ |T_4| \leq C_1^2 \left( 1 + \frac{C_6^2}{2} \right) \|(I - \Pi)\tilde{\mathbf{F}}\|^2 + \frac{1}{2} \|1 - \rho_1 \rho_2\|_{L^2(\mathbf{T})}^2. $$

If we now recombine all the terms, and use estimate (2.21) of the dissipation $\mathcal{D}(\mathbf{F})$, we obtain

$$ \begin{aligned} \frac{d}{dt} \Gamma(t) \leq & - \left(C_3 - \frac{\delta}{2}\right) \|1 - \rho_1 \rho_2\|_{L^2(\mathbf{T})}^2 - \left(C_4 - \delta C_1^2 \left(1 + \frac{C_6^2}{2}\right)\right) \|(I - \Pi)\tilde{\mathbf{F}}\|^2 \\ & - \delta \min(D_1, D_2) \|\Pi \tilde{\mathbf{F}}\|^2 + \delta \left(C_2 + \max(\rho_m^{-1}, \rho_M) C_1 C_5\right) \|(I - \Pi)\tilde{\mathbf{F}}\| \|\Pi \tilde{\mathbf{F}}\|. \end{aligned} $$

If $\delta > 0$ is chosen sufficiently small, we obtain two constants $K_\delta > 0$ and $\bar{K}_\delta > 0$ such that

$$ \frac{d}{dt} \Gamma(t) + K_\delta \|\tilde{\mathbf{F}}(t)\|^2 + \bar{K}_\delta \|1 - \rho_1 \rho_2\|_{L^2(\mathbf{T})}^2 \leq 0, $$

and we can conclude by using Lemma 5. □

### 3. THE DISCRETE SETTING

In this section, we recall the numerical scheme introduced in [6]. As discussed therein, the use of monotone numerical fluxes is crucial to ensure a discrete maximum principle. In [6], Lax-Friedrichs fluxes were employed, as they are well suited for extending the hypocoercivity framework to the discrete setting. However, these fluxes do not appear to dissipate the Boltzmann entropy—at least no rigorous proof of such a property is available—making them unsuitable for the present analysis. We therefore consider upwind fluxes, which are both monotone and, as will be shown below, compatible with the Boltzmann entropy structure.
**3.1. Notations.** Before introducing the numerical scheme, we fix the notations and describe the underlying mesh.

**3.1.1. Mesh.** Since implementing a numerical scheme on an unbounded domain is not feasible in practice, we first restrict the velocity domain to a bounded and symmetric interval $[-v^*, v^*]$. We consider a discretization of $[-v^*, v^*]$ composed of 2L control volumes, arranged symmetrically around the origin. Let $v_{j+\frac{1}{2}}$ denote the $2L+1$ interface points, indexed by $j \in \mathcal{J} := \{-L, \cdots, L\}$, so that

$$v_{-L+\frac{1}{2}} = -v^*, \quad v_{\frac{1}{2}} = 0, \quad v_{j+\frac{1}{2}} = -v_{j-\frac{1}{2}} \quad \forall j = 0, \cdots, L.$$

For simplicity, we assume that the velocity grid is uniform, namely that each cell $V_j = (v_{j-\frac{1}{2}}, v_{j+\frac{1}{2}})$ has the same length $\Delta v$. Denoting by $v_j$ the midpoint of the cell $V_j$, we also have $v_j = -v_{-j+1}$ for all $j = 1, \cdots, L$.

In the position variable, we consider a uniform discretization of the torus $\mathbb{T}$ into $N$ control volumes

$$\mathcal{X}_i := (x_{i-\frac{1}{2}}, x_{i+\frac{1}{2}}), \quad i \in \mathcal{I} := \mathbb{Z}/N\mathbb{Z},$$

of constant length $\Delta x$. Once again, uniformity of the mesh is assumed for clarity, although the analysis extends to non-uniform meshes. As in [5], we further impose that $N$ is odd, which in particular ensures the validity of a discrete Poincaré inequality on the torus for the centered discrete derivative operator.

Finally, the control volumes in phase space are defined by

$$K_{ij} := \mathcal{X}_i \times V_j, \quad \forall (i, j) \in \mathcal{I} \times \mathcal{J},$$

and the mesh size is given by $\Delta = (\Delta x, \Delta v)$.

In addition, we denote by $\Delta t > 0$ the time step used for the discrete time integration.

**3.1.2. Discrete velocity profiles.** We now define discrete versions of the velocity profiles $\chi_k$. Our goal is to ensure that these discrete functions satisfy discrete counterparts of assumptions (1.5), which play a crucial role and are repeatedly used throughout the continuous analysis.

For $k = 1, 2$, we consider discrete values $(\chi_{k,j})_{j \in \mathcal{J}} \in \mathbb{R}^\mathcal{J}$ that satisfy the following discrete counterparts of assumptions (1.5)

$$\begin{gathered} \chi_{k,j} > 0, \quad \chi_{k,j} = \chi_{k,-j+1} \quad \forall j = 1, \cdots, L, \\ \sum_{j \in \mathcal{J}} \Delta v \chi_{k,j} = 1, \\ 0 < \underline{D}_k \le D_k^\Delta \le \overline{D}_k < \infty, \quad Q_k^\Delta \le \overline{Q}_k < \infty, \end{gathered} \quad (3.1)$$

where

$$D_k^\Delta := \sum_{j \in \mathcal{J}} \Delta v |v_j|^2 \chi_{k,j}, \quad Q_k^\Delta := \sum_{j \in \mathcal{J}} \Delta v |v_j|^4 \chi_{k,j} = \sum_{j \in \mathcal{J}} \Delta v v_j^2 |v_j|^2 \chi_{k,j},$$

and $\underline{D}_k$, $\overline{D}_k$, $\underline{Q}_k$ are universal constants. Typically, we define $\chi_{k,j} = c_{\Delta v} \chi_k(v_j)$ and compute $c_{\Delta v}$ in such a way that the mass of $(\chi_{k,j})_j$ is 1. Note also that the symmetry properties imply

$$\sum_{j \in \mathcal{J}} \Delta v v_j \chi_{k,j} = 0. \quad (3.2)$$

**3.1.3. Discrete functional setting.** To reproduce at the discrete level the functional framework used in the continuous analysis, we now introduce discrete counterparts of the weighted $L^2$ scalar products and of the associated norms for both microscopic and macroscopic quantities.

For discrete microscopic quantities $\mathbf{F} = (f_{1,ij}, f_{2,ij})_{ij}^\top$ and $\mathbf{G} = (g_{1,ij}, g_{2,ij})_{ij}^\top$, we define the discrete weighted $L^2$ scalar product by

$$\langle \mathbf{F}, \mathbf{G} \rangle_{\Delta} := \sum_{(i,j) \in \mathcal{I} \times \mathcal{J}} \Delta x \Delta v \left( \frac{f_{1,ij} g_{1,ij}}{f_{1,j}^{\infty,*}} + \frac{f_{2,ij} g_{2,ij}}{f_{2,j}^{\infty,*}} \right), \quad (3.3)$$
and denote by $\|\cdot\|_{\Delta}$ the associated norm. For macroscopic quantities $\mathbf{U} = (u_{1,i}, u_{2,i})_{i \in \mathcal{I}}^{\top}$ and $\mathbf{W} = (w_{1,i}, w_{2,i})_{i \in \mathcal{I}}^{\top}$, we also define a discrete weighted scalar product in position $\langle \cdot, \cdot \rangle_{\Delta,x}$ as

$$ \langle \mathbf{U}, \mathbf{W} \rangle_{\Delta,x} = \sum_{i \in \mathcal{I}} \Delta x \left( \frac{u_{1,i} w_{1,i}}{\rho_1^{\infty,*}} + \frac{u_{2,i} w_{2,i}}{\rho_2^{\infty,*}} \right) \quad (3.4) $$

with $\|\cdot\|_{\Delta,x}$ the associated norm.

Due to our choices of discretization, we need to define several discrete gradients in position. Given a macroscopic quantity $u = (u_i)_{i \in \mathcal{I}}$, we define

* the discrete centered gradient $D_x^c u \in \mathbb{R}^{\mathcal{I}}$ given by

$$ (D_x^c u)_i = \frac{u_{i+1} - u_{i-1}}{2\Delta x} \quad \forall i \in \mathcal{I}, $$

* the discrete downstream gradient $D_x^- u \in \mathbb{R}^{\mathcal{I}}$ given by

$$ (D_x^- u)_i = \frac{u_i - u_{i-1}}{\Delta x} \quad \forall i \in \mathcal{I}, $$

* the discrete upstream gradient $D_x^+ u \in \mathbb{R}^{\mathcal{I}}$ given by

$$ (D_x^+ u)_i = \frac{u_{i+1} - u_i}{\Delta x} \quad \forall i \in \mathcal{I}. $$

It is straightforward to see that these discrete gradients satisfy the following properties:

$$ \frac{1}{2}(D_x^- + D_x^+) = D_x^c, \quad D_x^+ D_x^- = D_x^- D_x^+. \quad (3.5) $$

By convention, we assume that these discrete differentiation operators act componentwise on vector-valued quantities. In particular, for a discrete quantity $\mathbf{U} = (U_1, U_2)^{\top}$, we set

$$ D_x^* \mathbf{U} = (D_x^* U_1, D_x^* U_2)^{\top}, \quad \text{for } * \in \{c, +, -\}. $$

Using the definition of the discrete gradients and the periodic boundary conditions, we immediately have the following algebraic properties.

**Lemma 6.** For all discrete macroscopic quantities $\mathbf{U}$ and $\mathbf{W}$, it holds

$$ \langle D_x^c \mathbf{U}, \mathbf{W} \rangle_{\Delta,x} = - \langle \mathbf{U}, D_x^c \mathbf{W} \rangle_{\Delta,x}, \quad (3.6) $$

$$ \langle D_x^+ \mathbf{U}, \mathbf{W} \rangle_{\Delta,x} = - \langle \mathbf{U}, D_x^- \mathbf{W} \rangle_{\Delta,x}, \quad (3.7) $$

$$ \langle (D_x^+ D_x^- + D_x^- D_x^+) \mathbf{U}, \mathbf{W} \rangle_{\Delta,x} = -4 \langle D_x^c \mathbf{U}, D_x^c \mathbf{W} \rangle_{\Delta,x}, \quad (3.8) $$

$$ \Delta x \| D_x^c \mathbf{U} \|_{\Delta,x} \le \| \mathbf{U} \|_{\Delta,x}. \quad (3.9) $$

Let us also recall the discrete Poincaré inequality on the torus (see for example [5, Lemma 6] for a proof of this result).

**Lemma 7** (Discrete Poincaré inequality on the torus). Assume that the number of points $N$ used in the spatial discretization of the torus is odd. Then, there is a constant $C_P > 0$ independent of $\Delta x$ such that for all $\mathbf{U} = (u_{1,i}, u_{2,i})_{i \in \mathcal{I}}^{\top}$ satisfying $\sum_{i \in \mathcal{I}} \Delta x u_{k,i} = 0, k = 1, 2,$

$$ \| \mathbf{U} \|_{\Delta,x} \le C_P \| D_x^c \mathbf{U} \|_{\Delta,x}. $$

Note that in the above lemma, the constant $C_P$ approaches $\frac{1}{\pi}$ whereas in the continuous case the Poincaré constant is $\frac{1}{2\pi}$. This discrepancy arises from the specific choice of the centered discrete gradient, and we refer the reader to [5] for a detailed discussion.

**3.2. Definition of the numerical scheme.** The numerical scheme used to approximate (1.1)–(1.2) is based on a fully implicit finite-volume discretization in the phase space. More precisely, at each time step $t^n = n\Delta t$ and for each control volume $K_{ij}$, the solution is approximated by $\mathbf{F}_{ij}^n = (f_{1,ij}^n, f_{2,ij}^n)^{\top}$. The initial datum $\mathbf{F}^{\text{in}} = (f_1^{\text{in}}, f_2^{\text{in}})^{\top}$ is discretized as

$$ f_{k,ij}^{\text{in}} = \frac{1}{\Delta x \Delta v} \int_{K_{ij}} f_k^{\text{in}}(x,v) \, dx \, dv, \quad \forall (i,j) \in \mathcal{I} \times \mathcal{J}, \quad k = 1, 2. $$
Integrating (1.1)–(1.2) over each cell $K_{ij}$ and applying a fully implicit time discretization leads to the following scheme: for all $n \ge 0$, $i \in \mathcal{I}$, $j \in \mathcal{J}$,

$$ \frac{f_{1,ij}^{n+1} - f_{1,ij}^{n}}{\Delta t} + \frac{1}{\Delta x \Delta v} (\mathcal{F}_{i+\frac{1}{2},j}^{n+1} - \mathcal{F}_{i-\frac{1}{2},j}^{n+1}) = \chi_{1,j} - \rho_{2,i}^{n+1} f_{1,ij}^{n+1}, \quad (3.10) $$

$$ \frac{f_{2,ij}^{n+1} - f_{2,ij}^{n}}{\Delta t} + \frac{1}{\Delta x \Delta v} (\mathcal{G}_{i+\frac{1}{2},j}^{n+1} - \mathcal{G}_{i-\frac{1}{2},j}^{n+1}) = \chi_{2,j} - \rho_{1,i}^{n+1} f_{2,ij}^{n+1}, \quad (3.11) $$

with the discrete densities

$$ \rho_{k,i}^{n+1} = \sum_{j \in \mathcal{J}} \Delta v f_{k,ij}^{n+1}, \quad k = 1, 2. \quad (3.12) $$

We consider upwind fluxes, which are defined by

$$ \mathcal{F}_{i+\frac{1}{2},j}^{n} = \Delta v (v_{j}^{+} f_{1,ij}^{n} + v_{j}^{-} f_{1,i+1,j}^{n}), \quad (3.13) $$

$$ \mathcal{G}_{i+\frac{1}{2},j}^{n} = \Delta v (v_{j}^{+} f_{2,ij}^{n} + v_{j}^{-} f_{2,i+1,j}^{n}), \quad (3.14) $$

where

$$ a^{+} = \max(a, 0) = (|a| + a)/2, $$

$$ a^{-} = \min(a, 0) = (|a| - a)/2. $$

Remark that the scheme (3.10)–(3.11) clearly satisfies the discrete mass conservation of the difference $f_1 - f_2$:

$$ \sum_{(i,j) \in \mathcal{I} \times \mathcal{J}} \Delta x \Delta v (f_{1,ij}^{n} - f_{2,ij}^{n}) = \sum_{(i,j) \in \mathcal{I} \times \mathcal{J}} \Delta x \Delta v (f_{1,ij}^{\text{in}} - f_{2,ij}^{\text{in}}) \quad \forall n \ge 0. $$

Then, as in the continuous setting, we define $\rho^{\infty,*} > 0$ as the unique constant such that

$$ M_0 := \sum_{(i,j) \in \mathcal{I} \times \mathcal{J}} \Delta x \Delta v (f_{1,ij}^{\text{in}} - f_{2,ij}^{\text{in}}) = |\mathcal{T}| \left( \rho^{\infty,*} - \frac{1}{\rho^{\infty,*}} \right). \quad (3.15) $$

In particular, we take

$$ \rho^{\infty,*} = \frac{M_0 + \sqrt{M_0^2 + 4|\mathcal{T}|}}{2|\mathcal{T}|}. \quad (3.16) $$

It is clear that $\mathbf{F}_{ij}^{\infty,*} = (f_{1,ij}^{\infty,*}, f_{2,ij}^{\infty,*})^{\top}$ defined by

$$ f_{1,ij}^{\infty,*} = \rho_{1}^{\infty,*} \chi_{1,j}, \quad f_{2,ij}^{\infty,*} = \rho_{2}^{\infty,*} \chi_{2,j} \quad \forall (i,j) \in \mathcal{I} \times \mathcal{J}, \quad (3.17) $$

with $\rho_1^{\infty,*} = \rho_2^{\infty,*}$ and $\rho_1^{\infty,*} = 1/\rho_2^{\infty,*}$, is an equilibrium for the scheme (3.10)–(3.11).

As in [6, Theorem 3], it is possible to establish the existence of a solution to the numerical scheme (3.10)–(3.11), as well as a discrete maximum principle.

**Theorem 2.** Assume that the discrete velocity profiles satisfy (3.1), and that there exist constants $0 < \rho_m \le \rho_M < \infty$ such that for all $i \in \mathcal{I}$, $j \in \mathcal{J}$

$$ \rho_m \chi_{1,j} \le f_{1,ij}^{\text{in}} \le \rho_M \chi_{1,j}, \quad (3.18) $$

$$ \rho_M^{-1} \chi_{2,j} \le f_{2,ij}^{\text{in}} \le \rho_m^{-1} \chi_{2,j}. \quad (3.19) $$

Then the scheme (3.10)–(3.11) admits a solution $(f_{1,ij}^n, f_{2,ij}^n)_{i \in \mathcal{I}, j \in \mathcal{J}, n \ge 0}$ satisfying for all $i \in \mathcal{I}$, $j \in \mathcal{J}$ and $n \ge 0$,

$$ \rho_m \chi_{1,j} \le f_{1,ij}^{n} \le \rho_M \chi_{1,j}, \quad (3.20) $$

$$ \rho_M^{-1} \chi_{2,j} \le f_{2,ij}^{n} \le \rho_m^{-1} \chi_{2,j}. \quad (3.21) $$

The proof follows the approach in [6], with the Lax-Friedrichs fluxes replaced by the upwind ones. It proceeds in two steps: first, existence is shown for a truncated version of the scheme; second, one proves that any solution to this truncated scheme automatically satisfies the $L^\infty$ bounds (3.20)–(3.21), and therefore also solves the original scheme. This result relies on the monotonicity of the fluxes, a property verified by the upwind discretization (3.24)–(3.25). A detailed proof is provided in Appendix A.
Now, since we want to study the convergence of $(\tilde{\mathbf{F}}_{ij}^n)_{n,i,j}$ to $(\tilde{\mathbf{F}}_{ij}^{\infty,*})_{i,j}$, we introduce the deviation $\tilde{\mathbf{F}}_{ij}^n := \mathbf{F}_{ij}^n - \mathbf{F}_{ij}^{\infty,*}$ and rewrite the scheme (3.10)–(3.11) as

$$ \frac{\tilde{f}_{1,ij}^{n+1} - \tilde{f}_{1,ij}^n}{\Delta t} + \frac{1}{\Delta x \Delta v} (\tilde{\mathcal{F}}_{i+\frac{1}{2},j}^{n+1} - \tilde{\mathcal{F}}_{i-\frac{1}{2},j}^{n+1}) = \chi_{1,j} - \rho_{2,i}^{n+1} f_{1,ij}^{n+1}, \quad (3.22) $$

$$ \frac{\tilde{f}_{2,ij}^{n+1} - \tilde{f}_{2,ij}^n}{\Delta t} + \frac{1}{\Delta x \Delta v} (\tilde{\mathcal{G}}_{i+\frac{1}{2},j}^{n+1} - \tilde{\mathcal{G}}_{i-\frac{1}{2},j}^{n+1}) = \chi_{2,j} - \rho_{1,i}^{n+1} f_{2,ij}^{n+1}, \quad (3.23) $$

where the numerical fluxes $\tilde{\mathcal{F}}_{i+\frac{1}{2},j}$ and $\tilde{\mathcal{G}}_{i+\frac{1}{2},j}$ are defined by (3.13) and (3.14), replacing $f_k$ by $\tilde{f}_k$.

In [6], the discrete hypocoercivity of the linearized system associated with (1.1)–(1.2), discretized using Lax–Friedrichs fluxes, was studied. The central-flux-plus-viscosity formulation of these fluxes proved particularly convenient, as it allowed one to clearly distinguish the discrete analogues of the continuous terms from the additional contributions arising solely from numerical viscosity. Following this formalism, the upwind fluxes can be equivalently written as

$$ \tilde{\mathcal{F}}_{i+\frac{1}{2},j}^n = \Delta v \left( v_j \frac{\tilde{f}_{1,ij}^n + \tilde{f}_{1,i+1,j}^n}{2} - \frac{|v_j|}{2} (\tilde{f}_{1,i+1,j}^n - \tilde{f}_{1,ij}^n) \right), \quad (3.24) $$

$$ \tilde{\mathcal{G}}_{i+\frac{1}{2},j}^n = \Delta v \left( v_j \frac{\tilde{f}_{2,ij}^n + \tilde{f}_{2,i+1,j}^n}{2} - \frac{|v_j|}{2} (\tilde{f}_{2,i+1,j}^n - \tilde{f}_{2,ij}^n) \right). \quad (3.25) $$

This representation highlights the decomposition into a centered term and a stabilizing numerical diffusion, in direct analogy with the Lax–Friedrichs formulation.

**3.3. Definition and estimates of the discrete moments.** We now introduce the discrete moments, derive their evolution equations, and establish the corresponding estimates.

For $n \ge 0$, $i \in \mathcal{I}$ and $k = 1, 2$, let

$$ J_{k,i}^n := \sum_{j \in \mathcal{J}} \Delta v v_j f_{k,ij}^n, \quad S_{k,i}^n := \sum_{j \in \mathcal{J}} \Delta v (v_j^2 - D_k^\Delta) f_{k,ij}^n, \quad (3.26) $$

and remember that $\rho_{k,i}^n$ is defined by (3.12).

We further define the following “skewed” discrete moments, which will appear in the discrete moments equations as viscous terms. For $n \ge 0$, $i \in \mathcal{I}$ and $k = 1, 2$,

$$ J_{k,i}^{s,n} := \sum_{j \in \mathcal{J}} \Delta v |v_j| f_{k,ij}^n, \quad S_{k,i}^{s,n} := \sum_{j \in \mathcal{J}} \Delta v v_j |v_j| f_{k,ij}^n. \quad (3.27) $$

In the following, we use the vectorial notations:

$$ \rho_{\mathbf{F},i}^n = \begin{pmatrix} \rho_{1,i}^n \\ \rho_{2,i}^n \end{pmatrix}, \quad \mathbf{J}_{\mathbf{F},i}^n = \begin{pmatrix} J_{1,i}^n \\ J_{2,i}^n \end{pmatrix}, \quad \mathbf{S}_{\mathbf{F},i}^n = \begin{pmatrix} S_{1,i}^n \\ S_{2,i}^n \end{pmatrix}, \quad \mathbf{J}_{\mathbf{F},i}^{s,n} = \begin{pmatrix} J_{1,i}^{s,n} \\ J_{2,i}^{s,n} \end{pmatrix}, \quad \mathbf{S}_{\mathbf{F},i}^{s,n} = \begin{pmatrix} S_{1,i}^{s,n} \\ S_{2,i}^{s,n} \end{pmatrix}. \quad (3.28) $$

To lighten the notations, the moments of the deviation $\tilde{\mathbf{F}} = \mathbf{F} - \mathbf{F}^\infty$ are denoted by

$$ \rho_{\tilde{\mathbf{F}},i}^n := \tilde{\rho}_i^n, \quad \mathbf{J}_{\tilde{\mathbf{F}},i}^n := \tilde{\mathbf{J}}_i^n, \quad \mathbf{S}_{\tilde{\mathbf{F}},i}^n := \tilde{\mathbf{S}}_i^n, \quad \mathbf{J}_{\tilde{\mathbf{F}},i}^{s,n} := \tilde{\mathbf{J}}_i^{s,n}, \quad \mathbf{S}_{\tilde{\mathbf{F}},i}^{s,n} := \tilde{\mathbf{S}}_i^{s,n}. $$

Let us now define the discrete counterpart of the projection $\Pi$ as

$$ (\Pi^\Delta \mathbf{F})_{ij} := \begin{pmatrix} \rho_{1,i} \chi_{1,j} \\ \rho_{2,i} \chi_{2,j} \end{pmatrix}. \quad (3.29) $$

Using the equivalent formulation (3.24)–(3.25) of the upwind fluxes, the discrete moment equations are easily derived by multiplying the scheme (3.22)–(3.23) by $(1, v_j)^\top$ and summing over the velocity index $j$.
**Lemma 8 (Moments schemes).** The discrete moments satisfy the following schemes: for all $n \ge 0$, $i \in \mathcal{I}$,

$$ \frac{\tilde{\rho}_i^{n+1} - \tilde{\rho}_i^n}{\Delta t} + (D_x^c \tilde{\mathbf{J}}^{n+1})_i - \frac{\Delta x}{2} \left( (D_x^+ D_x^- + D_x^- D_x^+) \tilde{\mathbf{J}}^{s,n+1} \right)_i = (1 - \rho_{1,i}^{n+1} \rho_{2,i}^{n+1}) \begin{pmatrix} 1 \\ 1 \end{pmatrix}, \quad (3.30) $$

$$ \frac{\tilde{J}_{1,i}^{n+1} - \tilde{J}_{1,i}^n}{\Delta t} + (D_x^c \tilde{S}_1^{n+1})_i + D_1^\Delta (D_x^c \tilde{\rho}_1^{n+1})_i - \frac{\Delta x}{2} \left( (D_x^+ D_x^- + D_x^- D_x^+) \tilde{S}_1^{s,n+1} \right)_i = -\tilde{\rho}_{2,i}^{n+1} \tilde{J}_{1,i}^{n+1}, \quad (3.31) $$

$$ \frac{\tilde{J}_{2,i}^{n+1} - \tilde{J}_{2,i}^n}{\Delta t} + (D_x^c \tilde{S}_2^{n+1})_i + D_2^\Delta (D_x^c \tilde{\rho}_2^{n+1})_i - \frac{\Delta x}{2} \left( (D_x^+ D_x^- + D_x^- D_x^+) \tilde{S}_2^{s,n+1} \right)_i = -\tilde{\rho}_{1,i}^{n+1} \tilde{J}_{2,i}^{n+1}. \quad (3.32) $$

Similarly as in the continuous setting, the discrete moments of the deviation $\tilde{\mathbf{F}} = \mathbf{F} - \mathbf{F}^\infty$ can be estimated using purely algebraic properties.

**Lemma 9 (Discrete moments estimates).** The discrete moments satisfy the following estimates for all $t \ge 0$:

$$ \| \tilde{\rho}^n \|_{\Delta, x} = \| \Pi^\Delta \tilde{\mathbf{F}}^n \|_{\Delta}, \quad (3.33) $$

$$ \| \tilde{\mathbf{J}}^n \|_{\Delta, x} \le C_1^* \|(I - \Pi^\Delta) \tilde{\mathbf{F}}^n \|_{\Delta}, \quad (3.34) $$

$$ \| \tilde{\mathbf{S}}^n \|_{\Delta, x} \le C_2^* \|(I - \Pi^\Delta) \tilde{\mathbf{F}}^n \|_{\Delta}, \quad (3.35) $$

$$ \| \tilde{\mathbf{J}}^{s,n} \|_{\Delta, x} \le C_1^* \| \tilde{\mathbf{F}}^n \|_{\Delta}, \quad (3.36) $$

$$ \| \tilde{\mathbf{S}}^{s,n} \|_{\Delta, x} \le C_3^* \|(I - \Pi^\Delta) \tilde{\mathbf{F}}^n \|_{\Delta}, \quad (3.37) $$

where the constants are given by:

$$ C_1^* = \sqrt{\max \{\overline{D_1}, \overline{D_2}\}}, \quad (3.38) $$

$$ C_2^* = \sqrt{\max \{(\overline{Q_1} - \overline{D_1}^2), (\overline{Q_2} - \overline{D_2}^2)\}}, \quad (3.39) $$

$$ C_3^* = \sqrt{\max \{\overline{Q_1}, \overline{Q_2}\}}, \quad (3.40) $$

The proof of this lemma is obtained with the same algebraic manipulations as in the proof of Lemma 1, relying on the symmetry properties of the velocity mesh and of the discrete equilibrium profiles (3.1)–(3.2). The only additional step consists in estimating the new skewed moments, which can be handled using similar arguments.

**Proof.** We obtain estimate (3.36) thanks to definition (3.27) and the Cauchy-Schwarz inequality:

$$ \begin{aligned} \| \tilde{\mathbf{J}}^{s,n} \|_{\Delta, x}^2 &= \sum_{k=1,2} \frac{1}{\rho_k^{\infty,*}} \sum_{i \in \mathcal{I}} \Delta x \left( \sum_{j \in \mathcal{J}} \Delta v |v_j| \tilde{f}_{k,ij}^n \frac{\sqrt{\chi_{k,j}}}{\sqrt{\chi_{k,j}}} \right)^2 \\ &\le \sum_{k=1,2} \frac{1}{\rho_k^{\infty,*}} \sum_{i \in \mathcal{I}} \Delta x \left( \sum_{j \in \mathcal{J}} \Delta v v_j^2 \chi_{k,j} \right) \left( \sum_{j \in \mathcal{J}} \Delta v \frac{\left( \tilde{f}_{k,ij}^n \right)^2}{\chi_{k,j}} \right) \\ &\le \max(\overline{D_1}, \overline{D_2}) \| \tilde{\mathbf{F}}^n \|_{\Delta}^2. \end{aligned} $$

Estimate (3.37) is a consequence of the following relation (due to symmetry):

$$ \sum_{j \in \mathcal{J}} \Delta v v_j |v_j| \chi_{k,j} = 0. $$
Indeed, one has:

$$
\begin{align}
\|\tilde{\mathbf{S}}^{s,n}\|_{\Delta,x}^2 &= \sum_{k=1,2} \sum_{i \in \mathcal{I}} \Delta x \frac{1}{\rho_k^{\infty,*}} \left( \sum_{j \in \mathcal{J}} \Delta v v_j |v_j| \tilde{f}_{k,ij}^n \right)^2 \nonumber \\
&= \sum_{k=1,2} \sum_{i \in \mathcal{I}} \Delta x \frac{1}{\rho_k^{\infty,*}} \left( \sum_{j \in \mathcal{J}} \Delta v v_j |v_j| (\tilde{f}_{k,ij}^n - \tilde{\rho}_{k,i} \chi_{k,j}) \right)^2 \tag{3.41} \\
&\le \max \{\overline{Q}_1, \overline{Q}_2\} \|(I - \Pi^\Delta)\tilde{\mathbf{F}}^n\|_{\Delta}^2. \nonumber
\end{align}
$$

□

4. DISCRETE NONLINEAR HYPOCOERCIVITY

In this section, we adapt the proof of Section 2 to the discrete setting. Since algebraic ma-
nipulations directly carry over to the discrete setting, we shall focus on the main challenges
induced by the discretization.

The discrete analog of the relative Boltzmann entropy is defined as

$$
H^n := \sum_{k=1,2} \sum_{i,j} \Delta x \Delta v \left( f_{k,ij}^n \left( \log \left( \frac{f_{k,ij}^n}{f_{k,j}^{\infty,*}} \right) - 1 \right) + f_{k,j}^{\infty,*} \right). \quad (4.1)
$$

First, we introduce the discrete counterpart of Lemma 2.

**Lemma 10 (Discrete microscopic coercivity).** Under assumptions (3.1) and (3.18)-(3.19), let $\mathbf{F} = (\mathbf{F}_{ij})_{ij}$ be solution to the scheme (3.10)-(3.11). Then for every $n \ge 0$,

$$
\frac{H^{n+1} - H^n}{\Delta t} \le -D^{n+1} - \frac{C_4^*}{\Delta t} \|\tilde{\mathbf{F}}^{n+1} - \tilde{\mathbf{F}}^n\|_{\Delta}^2, \quad (4.2)
$$

where

$$
\mathcal{D}^n = \sum_{(i,j,m) \in \mathcal{I} \times \mathcal{J}^2} \Delta x \Delta v^2 (\chi_{1,j} \chi_{2,m} - f_{1,ij}^n f_{2,im}^n) \log \left( \frac{f_{1,ij}^n f_{2,im}^n}{\chi_{1,j} \chi_{2,m}} \right), \quad (4.3)
$$

and

$$
C_4^* = \frac{1}{2} \frac{\rho_m}{\rho_M} \min \left( \frac{\rho_1^{\infty,*}}{\rho_M}, \rho_2^{\infty,*} \rho_m \right).
$$

For the purely algebraic manipulations carried out in the continuous setting, the transpo-
sition to the discrete level is straightforward. The main difference arises from the absence of
a discrete chain rule, which forces us to rely instead on the concavity of the logarithm func-
tion. More precisely, the remainder term in the first-order Taylor expansion of s ↦→ log(s)
has a definite sign, which we exploit to obtain the desired estimate. This yields the second
term in (4.2), corresponding to an additional numerical dissipation. Furthermore, this ex-
tra dissipation originates purely from the time discretization and is of order Δt; hence, it
vanishes in the continuous-time limit Δt → 0.

Proof. By definition, we have

$$
H^{n+1} - H^n = \sum_{k=1,2} \sum_{i,j} \Delta x \Delta v \left[ f_{k,ij}^{n+1} \left( \log \left( \frac{f_{k,ij}^{n+1}}{f_{k,j}^{\infty}} \right) - 1 \right) - f_{k,ij}^{n} \left( \log \left( \frac{f_{k,ij}^{n}}{f_{k,j}^{\infty}} \right) - 1 \right) \right].
$$

For $k = 1, 2$, for all $i \in \mathcal{I}$, $j \in \mathcal{J}$, we can write

$$
\begin{align*}
& f_{k,ij}^{n+1} \left( \log \left( \frac{f_{k,ij}^{n+1}}{f_{k,j}^{\infty}} \right) - 1 \right) - f_{k,ij}^{n} \left( \log \left( \frac{f_{k,ij}^{n}}{f_{k,j}^{\infty}} \right) - 1 \right) = \\
& (f_{k,ij}^{n+1} - f_{k,ij}^{n}) \log \left( \frac{f_{k,ij}^{n+1}}{f_{k,j}^{\infty}} \right) + f_{k,ij}^{n} \left[ \log \left( \frac{f_{k,ij}^{n+1}}{f_{k,j}^{\infty}} \right) - \log \left( \frac{f_{k,ij}^{n}}{f_{k,j}^{\infty}} \right) \right] + f_{k,ij}^{n+1} - f_{k,ij}^{n}.
\end{align*}
$$

Using a Taylor expansion, the second term of the right hand side can be written as

$$
f_{k,ij}^{n} \left( \log(f_{k,ij}^{n+1}) - \log(f_{k,ij}^{n}) \right) = f_{k,ij}^{n} \left[ \frac{1}{f_{k,ij}^{n}} (f_{k,ij}^{n+1} - f_{k,ij}^{n}) - \frac{1}{2(\tilde{f}_{k,ij})^{2}} (f_{k,ij}^{n+1} - f_{k,ij}^{n})^{2} \right],
$$
with $\tilde{f}_{k,ij} \in (\min(f_{k,ij}^n, f_{k,ij}^{n+1}), \max(f_{k,ij}^n, f_{k,ij}^{n+1}))$. We thus have

$$H^{n+1} - H^n = \sum_{k=1,2} \sum_{i,j} \Delta x \Delta v (f_{k,ij}^{n+1} - f_{k,ij}^n) \log \left( \frac{f_{k,ij}^{n+1}}{f_{k,j}^\infty} \right) \\ - \frac{1}{2} \sum_{k=1,2} \sum_{i,j} \Delta x \Delta v \frac{f_{k,ij}^n}{(\tilde{f}_{k,ij})^2} (f_{k,ij}^{n+1} - f_{k,ij}^n)^2. \quad (4.4)$$

The first term is analog to its continuous counterpart, and is rewritten using the numerical scheme (3.10)-(3.11), yielding

$$\sum_{k=1,2} \sum_{i,j} \Delta x \Delta v (f_{k,ij}^{n+1} - f_{k,ij}^n) \log \left( \frac{f_{k,ij}^{n+1}}{f_{k,j}^{\infty,*}} \right) \quad (4.5) \\ = \Delta t \sum_{i,j} \left( -(\mathcal{F}_{i+\frac{1}{2}}^{n+1} - \mathcal{F}_{i-\frac{1}{2}}^{n+1}) + \Delta x \Delta v (\chi_{1,j} - \rho_{2,i}^{n+1} f_{1,ij}^{n+1}) \right) \log \left( \frac{f_{1,ij}^{n+1}}{f_{1,j}^{\infty,*}} \right) \\ + \Delta t \sum_{i,j} \left( -(\mathcal{G}_{i+\frac{1}{2}}^{n+1} - \mathcal{G}_{i-\frac{1}{2}}^{n+1}) + \Delta x \Delta v (\chi_{2,j} - \rho_{1,i} f_{2,ij}^{n+1}) \right) \log \left( \frac{f_{2,ij}^{n+1}}{f_{2,j}^{\infty,*}} \right).$$

The transport terms are handled in the same way for both species. We therefore focus on the first one. Using a discrete integration by parts on $i$, one has

$$A_T := - \sum_{i,j} (\mathcal{F}_{i+\frac{1}{2},j}^{n+1} - \mathcal{F}_{i-\frac{1}{2},j}^{n+1}) \log \left( \frac{f_{1,ij}^{n+1}}{f_{1,j}^{\infty,*}} \right) \\ = \sum_{i,j} \mathcal{F}_{i+\frac{1}{2},j}^{n+1} (\log (f_{1,i+1,j}^{n+1}) - \log (f_{1,ij}^{n+1})) \\ = \sum_{i \in \mathcal{I}, j > 0} \Delta v v_j f_{1,i,j} (\log (f_{1,i+1,j}^{n+1}) - \log (f_{1,ij}^{n+1})) \\ + \sum_{i \in \mathcal{I}, j < 0} \Delta v v_j f_{1,i+1,j}^{n+1} (\log (f_{1,i+1,j}^{n+1}) - \log (f_{1,ij}^{n+1})).$$

Then, by concavity of $s \mapsto \log(s)$, we have

$$A_T \le \sum_{i \in \mathcal{I}, j > 0} \Delta v v_j f_{1,i,j}^{n+1} \frac{f_{1,i+1,j}^{n+1} - f_{1,ij}^{n+1}}{f_{1,ij}^{n+1}} + \sum_{i \in \mathcal{I}, j < 0} \Delta v v_j f_{1,i+1,j}^{n+1} \frac{f_{1,i+1,j}^{n+1} - f_{1,ij}^{n+1}}{f_{1,i+1,j}^{n+1}} \\ \le 0,$$

where we have used the periodic boundary conditions in the physical space to conclude. The reaction terms in (4.5) are then handled in the same way as in the continuous setting, since their treatment involves only algebraic manipulations. One then obtains

$$\sum_{k=1,2} \sum_{i,j} \Delta x \Delta v (f_{k,ij}^{n+1} - f_{k,ij}^n) \log \left( \frac{f_{k,ij}^{n+1}}{f_{k,j}^\infty} \right) \le -\Delta t \mathcal{D}^{n+1}.$$

For the second term of the right-hand side of (4.4), we use the uniform-in-time $L^\infty$ bounds of $f_{k,ij}^n$. These bounds immediately transfer to $\tilde{f}_{k,ij}$, from which we obtain

$$\frac{f_{1,ij}^n}{(f_{1,ij}^*)^2} \ge \frac{\rho_m \chi_1}{(\rho_M \chi_1)^2} = \frac{\rho_m}{\rho_M^2} \frac{1}{\chi_1}, \quad \frac{f_{2,ij}^n}{(f_{2,ij}^*)^2} \ge \frac{\chi_2 \rho_m^2}{\rho_M \chi_2^2} = \frac{\rho_m^2}{\rho_M} \frac{1}{\chi_2},$$

and finally the expected result (4.2). □

By directly adapting the arguments used in the continuous setting to the discrete framework, we obtain the following discrete analogue of Lemma 3. It relates the discrete dissipation to the nonlinear reaction terms appearing in the discrete moment equation (3.30), as well as to the dissipation in the velocity variable.
**Lemma 11.** Under assumptions (3.1) and (3.18)-(3.19), the dissipation term defined by (4.3) satisfies for all $n \ge 0$

$$D^n \ge C_5^* \|1 - \rho_1^n \rho_2^n\|_{L^2(\mathbb{T})}^2 + C_6^* \|(I - \Pi^\Delta)\tilde{\mathbf{F}}^n\|_{\Delta}^2, \quad (4.6)$$

where $C_5^* = \rho_m/\rho_M$ and $C_6^* = (\rho_m/\rho_M)(\rho_1^{\infty,*} + \rho_2^{\infty,*})^{-1} \min(\rho_m, \rho_M^{-1})^2$.

We are now ready to state the full discrete hypocoercivity result.

**Theorem 3.** Under assumptions (3.1) and (3.18)-(3.19), and assuming that the number of points $N$ of the spatial discretization is odd, there exist constants $C \ge 1$ and $\kappa > 0$, which do not depend on the discretization parameters, such that

$$\|\mathbf{F}^n - \mathbf{F}^{\infty,*}\|_{\Delta} \le C \|\mathbf{F}^{\text{in}} - \mathbf{F}^{\infty,*}\|_{\Delta} e^{-\kappa t^n} \quad \forall n \ge 0. \quad (4.7)$$

In order to prove this result, we introduce the discrete modified entropy functional

$$\Gamma^n := H^n + \delta \langle \tilde{\mathbf{J}}^n, D_x^c \mathbf{\Phi}^n \rangle_{\Delta, x} + \frac{\delta}{2\Delta t} \sum_{i \in \mathcal{I}} \Delta x \left( (D_x^c \mathbf{\Phi}^n)_i - (D_x^c \mathbf{\Phi}^{n-1})_i \right)^2, \quad (4.8)$$

where $\mathbf{\Phi}^n = (\phi_1^n, \phi_2^n)^\top$, $(\phi_{k,i})_{i \in \mathcal{I}}$ being solution to the discrete Poisson equation

$$-D_x^c D_x^c \phi_{k,i}^n = \tilde{\rho}_{k,i}, \quad \sum_{i \in \mathcal{I}} \Delta x \phi_{k,i}^n = 0, \quad k = 1, 2, \quad (4.9)$$

and $\delta > 0$ is a small parameter to be chosen later. For an odd number of points $N$, existence and uniqueness of $\mathbf{\Phi}^n$ satisfying (4.9) is guaranteed (see [5]). The third term in (4.8) arises from the choice of a fully implicit time discretization and is consistent with zero in the limit $\Delta t \to 0$. This term, similarly to the additional numerical dissipation identified in Lemma 10, serves to compensate for the lack of a discrete chain rule, which necessitates to control the corresponding remainder terms.

Analogously to the continuous setting, it is straightforward to show that the discrete modified entropy $\Gamma^n$ defines a norm equivalent to $\|\tilde{\mathbf{F}}^n\|_{\Delta}^2$, as stated in the following lemma from [5].

**Lemma 12 (Equivalent norm).** There is $\delta^* > 0$ such that for all $\delta \in (0, \delta^*)$, there are positive constants $0 < c_\delta < C_\delta$ such that if $\mathbf{F}^n$ satisfies (2.3)-(2.4), one has

$$c_\delta \|\tilde{\mathbf{F}}^n\|^2 \le \Gamma^n \le C_\delta \|\tilde{\mathbf{F}}^n\|^2.$$

We also have the following estimates on $(\mathbf{\Phi}_i^n)_{i \in \mathcal{I}}$. They are obtained as in the continuous setting, using discrete integration by parts and the Poincaré inequality recalled in Lemma 7, although the upwind fluxes give rise to additional terms that must be controlled.

**Lemma 13.** One has for all $n \ge 0$

$$\|D_x^c \mathbf{\Phi}^n\|_{\Delta, x} \le C_P \|\Pi^\Delta \tilde{\mathbf{F}}^n\|_{\Delta}, \quad (4.10)$$

$$\|D_x^c \mathbf{\Phi}^{n+1} - D_x^c \mathbf{\Phi}^n\|_{\Delta, x} \le \Delta t (C_1^* \|(I - \Pi^\Delta)\tilde{\mathbf{F}}^{n+1}\|_{\Delta} \\ + C_7^* \|1 - \rho_1^{n+1} \rho_2^{n+1}\|_{L^2(\mathbb{T})} + C_1^* \|\tilde{\mathbf{F}}^{n+1}\|_{\Delta}), \quad (4.11)$$

where $C_P$ is the discrete Poincaré constant of Lemma 7, $C_1^*$ is defined by (3.38) and

$$C_7^* := C_P \sqrt{\rho_1^{\infty,*} + \rho_2^{\infty,*}}. \quad (4.12)$$

*Proof.* Estimate (4.10) can be proven exactly in the same way as (2.26), by using the discrete Poisson equation (4.9), applying the discrete Poincaré inequality (Lemma 7) and equality (3.33).
Regarding estimate (4.11), we obtain it by taking the difference of the scheme (4.9) at times n and n + 1 and using Lemma 8:

$$
\begin{align*}
\| D_x^c \mathbf{\Phi}^{n+1} - D_x^c \mathbf{\Phi}^n \|_{\Delta,x}^2 &= \langle \mathbf{\Phi}^{n+1} - \mathbf{\Phi}^n, -D_x^c D_x^c (\mathbf{\Phi}^{n+1} - \mathbf{\Phi}^n) \rangle_{\Delta,x} \\
&= \langle \mathbf{\Phi}^{n+1} - \mathbf{\Phi}^n, \tilde{\rho}^{n+1} - \tilde{\rho}^n \rangle_{\Delta,x} \\
&= \Delta t \left\langle \mathbf{\Phi}^{n+1} - \mathbf{\Phi}^n, -D_x^c \tilde{\mathbf{J}}^{n+1} + (1 - \rho_1^{n+1} \rho_2^{n+1}) \begin{pmatrix} 1 \\ 1 \end{pmatrix} \right\rangle_{\Delta,x} \\
&\quad + \Delta t \left\langle \mathbf{\Phi}^{n+1} - \mathbf{\Phi}^n, \frac{\Delta x}{2} (D_x^+ D_x^- + D_x^- D_x^+) \tilde{\mathbf{J}}^{s,n+1} \right\rangle_{\Delta,x}.
\end{align*}
$$

Using the discrete integration by parts (3.6) and (3.8), we get

$$
\| D_x^c \mathbf{\Phi}^{n+1} - D_x^c \mathbf{\Phi}^n \|_{\Delta,x}^2 = \Delta t \langle D_x^c \mathbf{\Phi}^{n+1} - D_x^c \mathbf{\Phi}^n, \tilde{\mathbf{J}}^{n+1} - 2\Delta x D_x^c \tilde{\mathbf{J}}^{s,n+1} \rangle_{\Delta,x} \\
+ \Delta t \left\langle \mathbf{\Phi}^{n+1} - \mathbf{\Phi}^n, (1 - \rho_1^{n+1} \rho_2^{n+1}) \begin{pmatrix} 1 \\ 1 \end{pmatrix} \right\rangle_{\Delta,x}.
$$

Now, thanks to the Cauchy-Schwarz inequality, estimate (3.9) and the discrete Poincaré inequality, we obtain

$$
\begin{align*}
\| D_x^c \mathbf{\Phi}^{n+1} - D_x^c \mathbf{\Phi}^n \|_{\Delta,x}^2 &\le \Delta t \| D_x^c \mathbf{\Phi}^{n+1} - D_x^c \mathbf{\Phi}^n \|_{\Delta,x} \left[ \| \tilde{\mathbf{J}}^{n+1} \|_{\Delta,x} + 2\Delta x \| D_x^c \tilde{\mathbf{J}}^{s,n+1} \|_{\Delta,x} \right] \\
&\quad + \Delta t \| \mathbf{\Phi}^{n+1} - \mathbf{\Phi}^n \|_{\Delta,x} \sqrt{\rho_1^{\infty,*} + \rho_2^{\infty,*}} \| 1 - \rho_1^{n+1} \rho_2^{n+1} \|_{L^2(\mathbb{T})} \\
&\le \Delta t \| D_x^c \mathbf{\Phi}^{n+1} - D_x^c \mathbf{\Phi}^n \|_{\Delta,x} \left( \| \tilde{\mathbf{J}}^{n+1} \|_{\Delta,x} + 2 \| \tilde{\mathbf{J}}^{s,n+1} \|_{\Delta,x} \right. \\
&\qquad \left. + C_7^* \| 1 - \rho_1^{n+1} \rho_2^{n+1} \|_{L^2(\mathbb{T})} \right),
\end{align*}
$$

which concludes the proof thanks to estimates (3.34) and (3.36). □

*Proof of Theorem 3.* We adapt the proof of Theorem 1 to the discrete setting. For the sake of clarity, the terms are numbered as follows: discrete counterparts of the continuous terms retain the same numbering (for instance, $T_1$ becomes $T_1^\Delta$), while additional terms arising from the time and space discretizations are denoted by $R_1^\Delta$, $R_2^\Delta$, and so on.

Using the microscopic coercivity (4.2), we have

$$
\Gamma^{n+1} - \Gamma^n \leq -\Delta t \mathcal{D}^{n+1} - C_4^* \|\tilde{\mathbf{F}}^{n+1} - \tilde{\mathbf{F}}^n\|_{\Delta}^2 + \delta \left( \langle \tilde{\mathbf{J}}^{n+1}, D_x^c \mathbf{\Phi}^{n+1} \rangle_{\Delta,x} - \langle \tilde{\mathbf{J}}^n, D_x^c \mathbf{\Phi}^n \rangle_{\Delta,x} \right) \\
+ \frac{\delta}{2\Delta t} \sum_{i \in \mathcal{I}} \Delta x \left[ \left( (D_x^c \mathbf{\Phi}^{n+1})_i - (D_x^c \mathbf{\Phi}^n)_i \right)^2 - \left( (D_x^c \mathbf{\Phi}^n)_i - (D_x^c \mathbf{\Phi}^{n-1})_i \right)^2 \right]. \quad (4.13)
$$

Let us start with the term corresponding to the discrete counterpart of $\frac{d}{dt}\langle \tilde{\mathbf{J}}, \partial_x \mathbf{\Phi} \rangle_x$:

$$
\langle \tilde{\mathbf{J}}^{n+1}, D_x^c \mathbf{\Phi}^{n+1} \rangle_{\Delta,x} - \langle \tilde{\mathbf{J}}^n, D_x^c \mathbf{\Phi}^n \rangle_{\Delta,x} = \langle \tilde{\mathbf{J}}^{n+1} - \tilde{\mathbf{J}}^n, D_x^c \mathbf{\Phi}^{n+1} \rangle_{\Delta,x} + T_4^\Delta + R_1^\Delta,
$$

where

$$
T_4^\Delta = \langle \tilde{\mathbf{J}}^{n+1}, D_x^c \mathbf{\Phi}^{n+1} - D_x^c \mathbf{\Phi}^n \rangle_{\Delta,x},
$$

$$
R_1^\Delta = \langle \tilde{\mathbf{J}}^n, D_x^c \mathbf{\Phi}^{n+1} \rangle_{\Delta,x} - \langle \tilde{\mathbf{J}}^{n+1}, D_x^c \mathbf{\Phi}^{n+1} \rangle_{\Delta,x} + \langle \tilde{\mathbf{J}}^{n+1}, D_x^c \mathbf{\Phi}^n \rangle_{\Delta,x} - \langle \tilde{\mathbf{J}}^n, D_x^c \mathbf{\Phi}^n \rangle_{\Delta,x}.
$$

Using the scheme (3.31)–(3.32) for $\tilde{\mathbf{J}}$, we can rewrite

$$
\langle \tilde{\mathbf{J}}^{n+1} - \tilde{\mathbf{J}}^n, D_x^c \mathbf{\Phi}^{n+1} \rangle_{\Delta,x} = T_1^\Delta + T_2^\Delta + T_3^\Delta + R_2^\Delta,
$$
where

$$
\begin{align*}
T_1^\Delta &= -\Delta t \langle D_x^c \tilde{S}^{n+1}, D_x^c \Phi^{n+1} \rangle_{\Delta,x}, \\
T_2^\Delta &= -\Delta t \langle \left( \frac{D_1^\Delta D_x^c \tilde{\rho}_1^{n+1}}{D_2^\Delta D_x^c \tilde{\rho}_2^{n+1}} \right), D_x^c \Phi^{n+1} \rangle_{\Delta,x}, \\
T_3^\Delta &= -\Delta t \langle \left( \frac{\tilde{\rho}_2^{n+1} \tilde{J}_1^{n+1}}{\tilde{\rho}_2^{n+1} \tilde{J}_1^{n+1}} \right), D_x^c \Phi^{n+1} \rangle_{\Delta,x}, \\
R_2^\Delta &= \frac{\Delta t \Delta x}{2} \langle (D_x^+ D_x^- + D_x^- D_x^+) \bar{\mathbf{S}}_{\tilde{\mathbf{F}}}^{n+1}, D_x^c \Phi^{n+1} \rangle_{\Delta,x}.
\end{align*}
$$

Let us first deal with the residual term $R_1^\Delta$. Reorganizing the terms and integrating by parts, we have

$$
\begin{align*}
R_1^\Delta &= \langle \tilde{\mathbf{J}}^n - \tilde{\mathbf{J}}^{n+1}, D_x^c \mathbf{\Phi}^{n+1} - D_x^c \mathbf{\Phi}^n \rangle_{\Delta,x} \\
&= \langle D_x^c \tilde{\mathbf{J}}^{n+1} - D_x^c \tilde{\mathbf{J}}^n, \mathbf{\Phi}^{n+1} - \mathbf{\Phi}^n \rangle_{\Delta,x}.
\end{align*}
$$

The next step is to replace the discrete space derivatives of $\tilde{\mathbf{J}}$ using the scheme (3.30), which yields $R_1^\Delta = R_{11}^\Delta + R_{12}^\Delta + R_{13}^\Delta$, with

$$
\begin{align*}
R_{11}^{\Delta} &= -\frac{1}{\Delta t} \sum_{i \in \mathcal{I}} \Delta x (\tilde{\rho}_{i}^{n+1} - 2\tilde{\rho}_{i}^{n} + \tilde{\rho}_{i}^{n-1}) (\mathbf{\Phi}_{i}^{n+1} - \mathbf{\Phi}_{i}^{n}), \\
R_{12}^{\Delta} &= \frac{\Delta x}{2} \sum_{i \in \mathcal{I}} \Delta x \left[ (D_{x}^{+}D_{x}^{-} + D_{x}^{-}D_{x}^{+})(\tilde{\mathbf{J}}^{s,n+1} - \tilde{\mathbf{J}}^{s,n}) \right]_{i} (\mathbf{\Phi}_{i}^{n+1} - \mathbf{\Phi}_{i}^{n}), \\
R_{13}^{\Delta} &= \sum_{i \in \mathcal{I}} \Delta x (\rho_{1,i}^{n+1}\rho_{2,i}^{n+1} - \rho_{1,i}^{n}\rho_{2,i}^{n}) (\mathbf{\Phi}_{i}^{n+1} - \mathbf{\Phi}_{i}^{n}).
\end{align*}
$$

As in [6], we combine the term $R_{11}^\Delta$ with the last line of (4.13) using the auxiliary equation (4.9), an integration by parts, and the identity $-a(a-b) = -(a^2-b^2)/2 - (a-b)^2/2$, with $a = (D_x^c(\mathbf{\Phi}^{n+1} - \mathbf{\Phi}_i^n))_i$ and $b = (D_x^c(\mathbf{\Phi}^n - \mathbf{\Phi}_i^{n-1}))_i$, therefore yielding

$$
\frac{1}{2\Delta t} \sum_{i \in \mathcal{I}} \Delta x \left[ \left( (D_x^c \mathbf{\Phi}^{n+1})_i - (D_x^c \mathbf{\Phi}^n)_i \right)^2 - \left( (D_x^c \mathbf{\Phi}^n)_i - (D_x^c \mathbf{\Phi}^{n-1})_i \right)^2 \right] + R_{11}^\Delta \le 0.
$$

Let us now turn to the estimates of $R_{12}^\Delta$ and $R_{13}^\Delta$. These terms will be bounded in terms of $\|\tilde{\mathbf{F}}^{n+1} - \tilde{\mathbf{F}}^n\|_{\Delta}^2$, which will ultimately be controlled by the additional dissipation term appearing in the microscopic coercivity estimate (4.2).

First of all, we apply two integration by parts and use the auxiliary scheme (4.9) to obtain

$$
\begin{align*}
R_{12}^{\Delta} &= 2\Delta x \sum_{i \in \mathcal{I}} \Delta x (\tilde{\mathbf{J}}_{i}^{s,n+1} - \tilde{\mathbf{J}}_{i}^{s,n}) (D_{x}^{c} D_{x}^{c} (\mathbf{\Phi}^{n+1} - \mathbf{\Phi}^{n}))_{i} \\
&= -2\Delta x \sum_{i \in \mathcal{I}} \Delta x (\tilde{\mathbf{J}}_{i}^{s,n+1} - \tilde{\mathbf{J}}_{i}^{s,n}) (\tilde{\rho}_{i}^{n+1} - \tilde{\rho}_{i}^{n}).
\end{align*}
$$

Using the moments estimates (3.33) and (3.36), and assuming without loss of generality that $\Delta x < 1$, we finally get

$$
|R_{12}^{\Delta}| \leq 2\Delta x C_1^* \|\tilde{\mathbf{F}}^{n+1} - \tilde{\mathbf{F}}^n\|_{\Delta}^2 \leq 2C_1^* \|\tilde{\mathbf{F}}^{n+1} - \tilde{\mathbf{F}}^n\|_{\Delta}^2. \quad (4.14)
$$

Using the Cauchy-Schwarz inequality, we obtain

$$
|R_{13}^{\Delta}| \le \| \rho_1^{n+1} \rho_2^{n+1} - \rho_1^n \rho_2^n \|_{L^2(\mathbf{T})} \| \mathbf{\Phi}^{n+1} - \mathbf{\Phi}^n \|_{\Delta, x}. \quad (4.15)
$$

We treat the first factor by linearizing and using the $L^\infty$ bounds:

$$
\begin{align*}
\| \rho_1^{n+1} \rho_2^{n+1} - \rho_1^n \rho_2^n \|_{L^2(\mathbf{T})}^2 &= \| \rho_1^{n+1} (\rho_2^{n+1} - \rho_2^n) + \rho_2^n (\rho_1^{n+1} - \rho_1^n) \|_{L^2(\mathbf{T})}^2 \\
&\le 2 \left( \rho_M^2 \| \rho_2^{n+1} - \rho_2^n \|_{L^2(\mathbf{T})}^2 + \rho_m^{-2} \| \rho_1^{n+1} - \rho_1^n \|_{L^2(\mathbf{T})}^2 \right) \\
&\le 2 \max(\rho_M^2 \rho_2^{n,*}, \rho_m^{-2} \rho_1^{n,*}) \| \rho_{\mathbf{F}}^{n+1} - \rho_{\mathbf{F}}^n \|_{\Delta, x}^2,
\end{align*}
$$
from which we deduce

$$ \|\rho_1^{n+1} \rho_2^{n+1} - \rho_1^n \rho_2^n\|_{L^2(\mathbb{T})} \le C_8^* \|\tilde{\mathbf{F}}^{n+1} - \tilde{\mathbf{F}}^n\|_{\Delta}, $$

with $C_8^* = 2 \max(\rho_M^2 \rho_2^{\infty,*}, \rho_m^{-2} \rho_1^{\infty,*})$. Now, for the second factor of (4.15), remark that by the Poincaré inequality and an integration by parts,

$$ \begin{aligned} \|\mathbf{\Phi}^{n+1} - \mathbf{\Phi}^n\|_{\Delta,x}^2 &\le C_P^2 \|D_x^c(\mathbf{\Phi}^{n+1} - \mathbf{\Phi}^n)\|_{\Delta,x}^2 \\ &\le -C_P^2 \langle \mathbf{\Phi}^{n+1} - \mathbf{\Phi}^n, D_x^c D_x^c(\mathbf{\Phi}^{n+1} - \mathbf{\Phi}^n) \rangle_{\Delta,x}. \end{aligned} $$

Using now the auxiliary scheme (4.9) and the Cauchy-Schwarz inequality, we get

$$ \begin{aligned} \|\mathbf{\Phi}^{n+1} - \mathbf{\Phi}^n\|_{\Delta,x}^2 &\le C_P^2 \langle \mathbf{\Phi}^{n+1} - \mathbf{\Phi}^n, \tilde{\rho}^{n+1} - \tilde{\rho}^n \rangle_{\Delta,x} \\ &\le C_P^2 \|\mathbf{\Phi}^{n+1} - \mathbf{\Phi}^n\|_{\Delta,x} \|\tilde{\rho}^{n+1} - \tilde{\rho}^n\|_{\Delta,x}, \end{aligned} $$

from which we obtain

$$ \|\mathbf{\Phi}^{n+1} - \mathbf{\Phi}^n\|_{\Delta,x} \le C_P^2 \|\tilde{\rho}^{n+1} - \tilde{\rho}^n\|_{\Delta,x} \le C_P^2 \|\tilde{\mathbf{F}}^{n+1} - \tilde{\mathbf{F}}^n\|_{\Delta}. $$

Finally, we get

$$ |R_{13}^{\Delta}| \le C_8^* C_P^2 \|\tilde{\mathbf{F}}^{n+1} - \tilde{\mathbf{F}}^n\|_{\Delta}^2. \quad (4.16) $$

The rest of the proof follows the same steps as in the continuous setting, with the additional task of controlling the numerical viscosity. The terms $T_k^{\Delta}$, for $k = 1, 2, 3$ correspond respectively to the discrete counterparts of the continuous terms $T_k$, and are handled in exactly the same manner, yielding

$$ \begin{aligned} T_1^{\Delta} &\le \Delta t C_2^* \|(I - \Pi^{\Delta}) \tilde{\mathbf{F}}^{n+1}\|_{\Delta} \|\Pi^{\Delta} \tilde{\mathbf{F}}^{n+1}\|_{\Delta}, \\ T_2^{\Delta} &\le -\Delta t \min(\underline{D}_1, \underline{D}_2) \|\Pi^{\Delta} \tilde{\mathbf{F}}^{n+1}\|_{\Delta}^2, \\ T_3^{\Delta} &\le \Delta t C_9^* \|(I - \Pi^{\Delta}) \tilde{\mathbf{F}}^{n+1}\|_{\Delta} \|\Pi^{\Delta} \tilde{\mathbf{F}}^{n+1}\|_{\Delta}, \end{aligned} $$

with $C_9^* = \max(\rho_m^{-1}, \rho_M) C_1^* C_P$.

The term $R_2^{\Delta}$ arises from the numerical viscosity introduced by the upwind fluxes. It can be estimated by combining (3.9) with the discrete moment estimates (3.33) and (3.37), which yields

$$ R_2^{\Delta} \le 2 \Delta t C_3^* \|(I - \Pi^{\Delta}) \tilde{\mathbf{F}}^{n+1}\|_{\Delta} \|\Pi^{\Delta} \tilde{\mathbf{F}}^{n+1}\|_{\Delta}. $$

Gathering all these estimates, we conclude that

$$ \langle \tilde{\mathbf{J}}^{n+1} - \tilde{\mathbf{J}}^n, D_x^c \mathbf{\Phi}^{n+1} \rangle_{\Delta,x} \le -\Delta t \min(\underline{D}_1, \underline{D}_2) \|\Pi^{\Delta} \tilde{\mathbf{F}}^{n+1}\|_{\Delta}^2 \\ + \Delta t (C_2^* + C_9^* + 2C_3^*) \|(I - \Pi^{\Delta}) \tilde{\mathbf{F}}^{n+1}\|_{\Delta} \|\Pi^{\Delta} \tilde{\mathbf{F}}^{n+1}\|_{\Delta}. $$

Finally, the term $T_4^{\Delta}$ is the discrete counterpart of the term $T_4$. Using the Cauchy-Schwarz inequality and estimates (3.34) and (4.11), we have

$$ T_4^{\Delta} \le \Delta t (C_1^*)^2 \|(I - \Pi^{\Delta}) \tilde{\mathbf{F}}^{n+1}\|_{\Delta}^2 + \Delta t C_1^* C_7^* \|(I - \Pi^{\Delta}) \tilde{\mathbf{F}}^{n+1}\|_{\Delta} \|1 - \rho_1^{n+1} \rho_2^{n+1}\|_{L^2(\mathbb{T})} \\ + \Delta t (C_1^*)^2 \|(I - \Pi^{\Delta}) \tilde{\mathbf{F}}^{n+1}\|_{\Delta} \|\tilde{\mathbf{F}}^{n+1}\|_{\Delta}. $$

Then, using Young's inequality, we have for $\eta > 0$ (which will be chosen conveniently later) some constants $C_{\eta} > 0$ and $C_{10}^*$ such that

$$ T_4^{\Delta} \le \Delta t C_{\eta} \|(I - \Pi^{\Delta}) \tilde{\mathbf{F}}^{n+1}\|_{\Delta}^2 + \Delta t C_{10}^* \|1 - \rho_1^{n+1} \rho_2^{n+1}\|_{L^2(\mathbb{T})}^2 + \Delta t \frac{\eta}{2} \|\tilde{\mathbf{F}}^{n+1}\|_{\Delta}^2. $$

By collecting the estimates of all the terms, we obtain

$$ \begin{aligned} \Gamma^{n+1} - \Gamma^n \le & -\Delta t \mathcal{D}^{n+1} - C_4^* \|\tilde{\mathbf{F}}^{n+1} - \tilde{\mathbf{F}}^n\|_{\Delta}^2 - \Delta t \delta \min(\underline{D}_1, \underline{D}_2) \|\Pi^{\Delta} \tilde{\mathbf{F}}^{n+1}\|_{\Delta}^2 \\ & + \Delta t \delta (C_2^* + C_9^* + 2C_3^*) \|(I - \Pi^{\Delta}) \tilde{\mathbf{F}}^{n+1}\|_{\Delta} \|\Pi^{\Delta} \tilde{\mathbf{F}}^{n+1}\|_{\Delta} \\ & + \Delta t \delta C_{\eta} \|(I - \Pi^{\Delta}) \tilde{\mathbf{F}}^{n+1}\|_{\Delta}^2 + \Delta t \delta C_{10}^* \|1 - \rho_1^{n+1} \rho_2^{n+1}\|_{L^2(\mathbb{T})}^2 \\ & + \delta \frac{\eta}{2} \|\tilde{\mathbf{F}}^{n+1}\|_{\Delta}^2 + 2\delta C_1^* \|\tilde{\mathbf{F}}^{n+1} - \tilde{\mathbf{F}}^n\|_{\Delta}^2 + \delta C_8^* C_P^2 \|\tilde{\mathbf{F}}^{n+1} - \tilde{\mathbf{F}}^n\|_{\Delta}^2. \end{aligned} $$

Now, using that

$$ \|\tilde{\mathbf{F}}^{n+1}\|_{\Delta}^2 = \|\Pi^{\Delta} \tilde{\mathbf{F}}^{n+1}\|_{\Delta}^2 + \|(I - \Pi^{\Delta}) \tilde{\mathbf{F}}^{n+1}\|_{\Delta}^2, $$
and Lemma 11, we finally get

$$
\begin{aligned}
\Gamma^{n+1} - \Gamma^n \leq & - (C_4^* - \delta(2C_1^* + C_8^*C_P^2)) \|\tilde{\mathbf{F}}^{n+1} - \tilde{\mathbf{F}}^n\|_{\Delta}^2 \\
& - \delta\Delta t \left(\min(\underline{D}_1, \underline{D}_2) - \frac{\eta}{2}\right) \|\Pi^{\Delta}\tilde{\mathbf{F}}^{n+1}\|_{\Delta}^2 \\
& + \delta\Delta t (C_2^* + C_9^* + 2C_3^*) \|(I - \Pi^{\Delta})\tilde{\mathbf{F}}^{n+1}\|_{\Delta} \|\Pi^{\Delta}\tilde{\mathbf{F}}^{n+1}\|_{\Delta} \\
& - \Delta t \left(C_6^* - \delta\left(C_{\eta} + \frac{\eta}{2}\right)\right) \|(I - \Pi^{\Delta})\tilde{\mathbf{F}}^{n+1}\|_{\Delta}^2 \\
& - \Delta t (C_5^* - \delta C_{10}^*) \|1 - \rho_1^{n+1} - \rho_2^{n+1}\|_{L^2(\mathbb{T})}^2.
\end{aligned}
$$

First, we fix $\eta > 0$ in such a way that the second term of the right hand side is nonpositive, namely $\eta \le 2\min(\underline{D}_1, \underline{D}_2)$. Second, we also assume that $0 < \delta < \delta_0^*$, with

$$ \delta_0^* = \min \left( \frac{C_4^*}{2C_1^* + C_8^*C_P^2}, \frac{C_5^*}{C_{10}^*}, \frac{C_9^*}{C_\eta + \eta/2} \right), $$

ensuring that the first and the last two terms are nonpositive.

With such an $\eta$ and $\delta$, we now have

$$
\frac{\Gamma^{n+1} - \Gamma^n}{\Delta t} \leq -\delta\alpha_1 \|\Pi^\Delta \tilde{\mathbf{F}}^{n+1}\|_{\Delta}^2 - (C_6^* - \delta\alpha_2) \|(I - \Pi^\Delta) \tilde{\mathbf{F}}^{n+1}\|_{\Delta}^2 \\
+ \delta\alpha_3 \|(I - \Pi^\Delta) \tilde{\mathbf{F}}^{n+1}\|_{\Delta} \|\Pi^\Delta \tilde{\mathbf{F}}^{n+1}\|_{\Delta},
$$

with $\alpha_1 = \min(\underline{D}_1, \underline{D}_2) - \eta/2$, $\alpha_2 = C_\eta + \eta/2$, and $\alpha_3 = C_2^* + C_9^* + 2C_3^*$. If we choose $\delta < \delta_1$ with

$$ \delta_1 < \frac{C_6^* \alpha_1}{\alpha_3^2 + \alpha_1 \alpha_2}, $$

we then get

$$ \frac{\Gamma^{n+1} - \Gamma^n}{\Delta t} \leq - \frac{\delta \alpha_1}{2} \|\Pi^\Delta \tilde{\mathbf{F}}^{n+1}\|_{\Delta}^2 - \frac{1}{2} (C_6^* - \delta \alpha_2) \|(I - \Pi^\Delta) \tilde{\mathbf{F}}^{n+1}\|_{\Delta}^2, $$

from which we conclude the proof using Lemma 12. $\square$

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

ACKNOWLEDGMENTS

TR received funding from the European Union's Horizon Europe research and innovation program under the Marie Skłodowska-Curie Doctoral Network DataHyking (Grant No. 101072546), and from UniCAJEDI Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR-15-IDEX-01.

TL received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No 865711) and from the European Union's Horizon Europe research and innovation program under the Marie Skłodowska-Curie Doctoral Network DataHyking (Grant No. 101072546).

MBC benefits from the support of the French government "Investissements d'Avenir" program integrated to France 2030, under the reference ANR-11-LABX-0020-01, and ANR Project Cookie (ANR-25-CE40-5565-02).

TR would also like to thank Isabelle Tristani for fruitful discussions about the continuous part of this work.

REFERENCES

[1] ADDALA, L., DOLBEAULT, J., LI, X., AND TAYEB, M. L. $L^2$-hypocoercivity and large time asymptotics of the linearized Vlasov-Poisson-Fokker-Planck system. *J. Stat. Phys.* **184**, 1 (2021), Paper No. 4, 34.

[2] AYI, N., HERDA, M., HIVERT, H., AND TRISTANI, I. On a structure-preserving numerical method for fractional Fokker-Planck equations. *Math. Comput.* **92**, 340 (2023), 635–693.

[3] BENSOUSSAN, A., LIONS, J.-L., AND PAPANICOLAU, G. C. Boundary Layers and Homogenization of Transport Processes. *Publ. Res. I. Math. Sci.* (1979).

[4] BESSEMOULIN-CHATARD, M., AND FILBET, F. A finite volume scheme for nonlinear degenerate parabolic equations. *SIAM J. Sci. Comput.* **34**, 5 (2012), B559–B583.

[5] BESSEMOULIN-CHATARD, M., HERDA, M., AND REY, T. Hypocoercivity and diffusion limit of a finite volume scheme for linear kinetic equations. *Math. Comput.* **89**, 323 (2020), 1093–1133.

[6] BESSEMOULIN-CHATARD, M., LAIDIN, T., AND REY, T. Discrete hypocoercivity for a nonlinear kinetic reaction model. *IMA J. Numer. Anal.* **45**, 4 (2025), 2204–2247.

[7] BLAUSTEIN, A., AND FILBET, F. A structure and asymptotic preserving scheme for the Vlasov-Poisson-Fokker-Planck model. *J. Comput. Phys.* **498** (2024), 112693.
[8] BLAUSTEIN, A., AND FILBET, F. On a discrete framework of hypocoercivity for kinetic equations. *Math. Comput.* **93**, 345 (2024), 163–202.

[9] BOUIN, E., DOLBEAULT, J., LAFLECHE, L., AND SCHMEISER, C. Fractional Hypocoercivity. *Commun. Math. Phys.* **390** (2022), 1369–1411.

[10] BOUIN, E., DOLBEAULT, J., MISCHLER, S., MOUHOT, C., AND SCHMEISER, C. Hypocoercivity without confinement. *Pure Appl. Anal.* **2**, 2 (2020), 203–232.

[11] BOUIN, E., DOLBEAULT, J., AND ZIVIANI, L. $L^2$ hypocoercivity methods for kinetic Fokker-Planck equations with factorised Gibbs states. In *Kolmogorov operators and their applications*. Singapore: Springer, 2024, pp. 23–56.

[12] BOUIN, E., AND ZIVIANI, L. Convergence to a non-explicit steady state in non-factorized kinetic Fokker-Planck equations with (very) weak velocity confinements. Preprint, arXiv:2510.12331, 2025.

[13] BURGER, M., CARRILLO, J. A., AND WOLFRAM, M.-T. A mixed finite element method for nonlinear diffusion equations. *Kinet. Relat. Mod.* **3**(1), 1 (2010), 59–83.

[14] CALVEZ, V., RAOUL, G., AND SCHMEISER, C. Confinement by biased velocity jumps: Aggregation of escherichia coli. *Kinet. Relat. Mod.* **8**, 4 (2015), 651–666.

[15] CARRAPATOSO, K., DOLBEAULT, J., HÉRAU, F., MISCHLER, S., MOUHOT, C., AND SCHMEISER, C. Special macroscopic modes and hypocoercivity. *J. Eur. Math. Soc.* (2024).

[16] CHAINAIS-HILLAIRET, C., AND FILBET, F. Asymptotic behaviour of a finite-volume scheme for the transient drift-diffusion model. *IMA J. Numer. Anal.* **27**, 4 (02 2007), 689–716.

[17] CHAINAIS-HILLAIRET, C., AND HERDA, M. Large-time behaviour of a family of finite volume schemes for boundary-driven convection-diffusion equations. *IMA J. Numer. Anal.* **40**, 4 (Oct. 2020), 2473–2505.

[18] DIETERT, H., HÉRAU, F., HUTRIDURGA, H., AND MOUHOT, C. Trajectorial hypocoercivity and application to control theory. In *Séminaire Laurent Schwartz—Équations aux dérivées partielles et applications. Année 2021–2022*. Inst. Hautes Études Sci., Bures-sur-Yvette, [2021–2022] ©2021–2022, pp. Exp. No. VIII, 10.

[19] DOLBEAULT, J., AND LI, X. Phi-Entropies: convexity, coercivity and hypocoercivity for Fokker-Planck and kinetic Fokker-Planck equations. *Math. Models Method Appl. Sci.* **28**, 13 (2018), 2637–2666.

[20] DOLBEAULT, J., MOUHOT, C., AND SCHMEISER, C. Hypocoercivity for linear kinetic equations conserving mass. *Trans. Amer. Math. Soc.* **367**, 6 (2015), 3807–3828.

[21] DONG, Z., GEORGOULIS, E. H., AND HERBERT, P. J. Asymptotic numerical hypocoercivity of the space-time discontinuous Galerkin method for Kolmogorov equation. Preprint, arXiv:2505.19222, 2025.

[22] DONG, Z., GEORGOULIS, E. H., AND HERBERT, P. J. A hypocoercivity-exploiting stabilized finite element method for Kolmogorov equation. *SIAM J. Numer. Anal.* **63**, 3 (2025), 1105–1127.

[23] DUJARDIN, G., HÉRAU, F., AND LAFITTE-GODILLON, P. Coercivity, hypocoercivity, exponential time decay and simulations for discrete Fokker-Planck equations. *Numer. Math.* **144** (2018).

[24] EVANS, J. Hypocoercivity in phi-entropy for the linear relaxation boltzmann equation on the torus. *SIAM J. Math. Anal.* **53**, 2 (2021), 1357–1378.

[25] FAVRE, G., PIRNER, M., AND SCHMEISER, C. Hypocoercivity and reaction-diffusion limit for a nonlinear generation-recombination model. *Arch. Ration. Mech. Anal.* **247**, 4 (2023), 15. Id/No 72.

[26] FILBET, F., AND HERDA, M. A finite volume scheme for boundary-driven convection-diffusion equations with relative entropy structure. *Numer. Math.* (2017).

[27] GEORGOULIS, E. H. Hypocoercivity-compatible finite element methods for the long-time computation of Kolmogorov's equation. *SIAM J. Numer. Anal.* **59**, 1 (2021), 173–194.

[28] GOSSE, L., AND TOSCANI, G. Identification of Asymptotic Decay to Self-Similarity for One-Dimensional Filtration Equations. *SIAM J. Numer. Anal.* **43**, 6 (2006), 2590–2606.

[29] GROSSE, B. Fully spectral scheme for the linear BGK equation on the whole space. Preprint, arXiv:2502.14396, 2025.

[30] HÉRAU, F. Introduction to hypocoercive methods and applications for simple linear inhomogeneous kinetic models. In *Lectures on the Analysis of Nonlinear Partial Differential Equations*, no. MLM 5 in Morning Side Lectures in Mathematics Series. International Press, 2017, pp. 119–147.

[31] HÉRAU, F., AND NIER, F. Isotropic hypoellipticity and trend to the equilibrium for the Fokker-Planck equation with high degree potential. *Arch. Rational Mech. Anal.* **171**, 2 (2004), 151–218.

[32] HÖRMANDER, L. Hypoelliptic second order differential equations. *Acta Math.* **119** (1967), 147–171.

[33] NEUMANN, L., AND SCHMEISER, C. A kinetic reaction model: decay to equilibrium and macroscopic limit. *Kinet. Relat. Mod.* **9**, 3 (2016), 571–585.

[34] PORRETTA, A., AND ZUAZUA, E. Numerical hypocoercivity for the Kolmogorov equation. *Math. Comput.* **86**, 303 (2017), 97–119.

[35] TEMAM, R. *Navier-Stokes equations. Theory and numerical analysis. 3rd (rev.) ed*, vol. 2 of *Stud. Math. Appl.* Elsevier, Amsterdam, 1984.

[36] VILLANI, C. Hypocoercivity. *Mem. Amer. Math. Soc.* (Sept. 2009).
APPENDIX A. PROOF OF THE EXISTENCE AND MAXIMUM PRINCIPLE

To prove Theorem 2, we follow along the lines of [6, Theorem 3], adapted to the upwind fluxes. Let us first introduce the following truncated version of the scheme (3.10)–(3.11):

$$
\frac{f_{1,ij}^{n+1} - f_{1,ij}^{n}}{\Delta t} + \frac{1}{\Delta x \Delta v} (\mathcal{F}_{i+\frac{1}{2},j}^{n+1} - \mathcal{F}_{i-\frac{1}{2},j}^{n+1}) = \chi_{1,j} - \bar{\rho}_{2,i}^{n+1} \bar{f}_{1,ij}^{n+1}, \quad (A.1)
$$

$$
\frac{f_{2,ij}^{n+1} - f_{2,ij}^{n}}{\Delta t} + \frac{1}{\Delta x \Delta v} (\mathcal{G}_{i+\frac{1}{2},j}^{n+1} - \mathcal{G}_{i-\frac{1}{2},j}^{n+1}) = \chi_{2,j} - \bar{\rho}_{1,i}^{n+1} \bar{f}_{2,ij}^{n+1}, \quad (A.2)
$$

with the upwind fluxes (3.13)–(3.14), and where the truncated quantities are defined by

$$
\bar{f}_{1,ij} := \begin{cases} \rho_m \chi_{1,j} & \text{if } f_{1,ij} \le \rho_m \chi_{1,j}, \\ f_{1,ij} & \text{if } \rho_m \chi_{1,j} \le f_{1,ij} \le \rho_M \chi_{1,j}, \\ \rho_M \chi_{1,j} & \text{if } f_{1,ij} \ge \rho_M \chi_{1,j}, \end{cases} \quad (A.3)
$$

and

$$
\bar{f}_{2,ij} := \begin{cases} \rho_M^{-1} \chi_{2,j} & \text{if } f_{2,ij} \le \rho_M^{-1} \chi_{2,j}, \\ f_{2,ij} & \text{if } \rho_M^{-1} \chi_{2,j} \le f_{2,ij} \le \rho_m^{-1} \chi_{2,j}, \\ \rho_m^{-1} \chi_{2,j} & \text{if } f_{2,ij} \ge \rho_m^{-1} \chi_{2,j}. \end{cases} \quad (A.4)
$$

The corresponding truncated densities are then given by

$$
\bar{\rho}_{1,i} = \sum_{j \in \mathcal{J}} \Delta v \bar{f}_{1,ij}, \quad \bar{\rho}_{2,i} = \sum_{j \in \mathcal{J}} \Delta v \bar{f}_{2,ij} \quad \forall i \in \mathcal{I}.
$$

The proof of the existence of solutions to (A.1)–(A.2) relies on the following result which
is proven for example in [35, Lemma 1.4, Chapter II].

**Lemma 14.** Let $X$ be a finite dimensional Hilbert space with scalar product $\langle \cdot, \cdot \rangle_X$ and associated norm $\| \cdot \|_X$. Let $P: X \to X$ be a continuous mapping such that $\langle P(\xi), \xi \rangle_X > 0$ for all $\xi \in X$ such that $\|\xi\|_X = k$ for some fixed $k > 0$. Then there exists $\xi_0 \in X$ such that $\|\xi_0\|_X \le k$ and $P(\xi) = 0$.

We apply this result with $X = (\mathbb{R}^{(2L+1)N})^2$ where for $\mathbf{F} = (f_1, f_2)^\top \in X$ and $\mathbf{G} = (g_1, g_2)^\top \in X$, we set

$$
\langle \mathbf{F}, \mathbf{G} \rangle_X = \langle f_1, g_1 \rangle + \langle f_2, g_2 \rangle,
$$

$\langle \cdot, \cdot \rangle$ being the classical Euclidian dot product on $\mathbb{R}^{(2L+1)N}$. Then, let us define

$$
P : \begin{pmatrix} f_1 \\ f_2 \end{pmatrix} \mapsto \begin{pmatrix} \Delta x \Delta v (f_{1,ij} - f_{1,ij}^n) + \Delta t (\mathcal{F}_{i+\frac{1}{2},j} - \mathcal{F}_{i-\frac{1}{2},j}) - \Delta x \Delta v \Delta t (\chi_{1,j} - \bar{\rho}_{2,i} \bar{f}_{1,ij}) \\ \Delta x \Delta v (f_{2,ij} - f_{2,ij}^n) + \Delta t (\mathcal{G}_{i+\frac{1}{2},j} - \mathcal{G}_{i-\frac{1}{2},j}) - \Delta x \Delta v \Delta t (\chi_{1,j} - \bar{\rho}_{1,i} \bar{f}_{2,ij}) \end{pmatrix},
$$

where $\mathcal{F}$ and $\mathcal{G}$ are the numerical fluxes (3.24)–(3.25) computed with $f_{1,ij}$ and $f_{2,ij}$. We
shall show that

$$
\langle P(\mathbf{F}), \mathbf{F} \rangle_X > 0, \quad \text{for all } \mathbf{F} = (f_1, f_2)^\top \text{ such that } \|\mathbf{F}\|_X = k,
$$

where $k$ is taken large enough. Indeed, this scalar product splits into 3 terms,

$$
\langle P(\mathbf{F}), \mathbf{F} \rangle_X = A_1 + A_2 + A_3,
$$

where

$$
A_1 = \sum_{(i,j) \in \mathcal{I} \times \mathcal{J}} \Delta x \Delta v \left( (f_{1,ij} - f_{1,ij}^n) f_{1,ij} + (f_{2,ij} - f_{2,ij}^n) f_{2,ij} \right),
$$

$$
A_2 = \Delta t \sum_{(i,j) \in \mathcal{I} \times \mathcal{J}} \left( (\mathcal{F}_{i+\frac{1}{2},j} - \mathcal{F}_{i-\frac{1}{2},j}) f_{1,ij} + (\mathcal{G}_{i+\frac{1}{2},j} - \mathcal{G}_{i-\frac{1}{2},j}) f_{2,ij} \right),
$$

$$
A_3 = -\Delta t \sum_{(i,j) \in \mathcal{I} \times \mathcal{J}} \Delta x \Delta v \left( (\chi_{1,j} - \bar{\rho}_{2,i} \bar{f}_{1,ij}) f_{1,ij} + (\chi_{2,j} - \bar{\rho}_{1,i} \bar{f}_{2,ij}) f_{2,ij} \right).
$$

Using the relation $a(a-b) \ge (a^2-b^2)/2$ one gets

$$
A_1 \ge \frac{1}{2} ( \| \mathbf{F} \|_X^2 - \| \mathbf{F}^n \|_X^2 ). \quad (A.5)
$$
Focusing on the first species, by definition of the numerical fluxes, one can use similar computations as in the proof of Lemma 10. In particular, using (3.13), a discrete integration by parts and the periodic boundary conditions, one obtains

$$
\begin{align}
\sum_{i,j} (\mathcal{F}_{i+\frac{1}{2},j} - \mathcal{F}_{i-\frac{1}{2},j}) f_{1,ij} &= -\sum_{i,j} \mathcal{F}_{i+\frac{1}{2},j} (f_{1,i+1,j} - f_{1,ij}) \nonumber \\
&= -\sum_{i \in \mathcal{I}, j > 0} \Delta v v_j f_{1,ij} (f_{1,i+1,j} - f_{1,i,j}) \nonumber \\
&\quad - \sum_{i \in \mathcal{I}, j < 0} \Delta v v_j f_{1,i+1,j} (f_{1,i+1,j} - f_{1,ij}) \tag{A.6} \\
&\geq \frac{1}{2} \sum_{i,j} \Delta v v_j ((f_{1,ij})^2 - (f_{1,i+1,j})^2) = 0. \nonumber
\end{align}
$$

Now, by assumptions (3.1), the velocity profiles are bounded:

$$
\exists C_{\chi} > 0 \quad \text{such that} \quad \forall j \in \mathcal{J}, \quad 0 \le \chi_{1,j}, \chi_{2,j} \le C_{\chi}.
$$

Then, by definition of the truncated quantities, there exists a constant $C > 0$ such that

$$
|\chi_{1,j} - \bar{\rho}_{2,i} \bar{f}_{1,ij}| \le C \quad \text{and} \quad |\chi_{2,j} - \bar{\rho}_{1,i} \bar{f}_{2,ij}| \le C.
$$

Applying the Cauchy-Schwarz inequality on $A_3$, one gets

$$
A_3 \geq -\Delta t C \|\mathbf{F}\|_{X}. \qquad (A.7)
$$

Combining (A.5), (A.6) and (A.7) yields the estimate

$$
\langle P(\mathbf{F}), \mathbf{F} \rangle_X \geq \frac{1}{2} (\|\mathbf{F}\|_X^2 - \|\mathbf{F}^n\|_X^2) - \Delta t C \|\mathbf{F}\|_X.
$$

The right hand side is a second order polynomial in $\|\mathbf{F}\|_X$ with a positive leading coefficient.
Since $\|\mathbf{F}^n\|_X$ is a constant in this context, there exists $k > 0$ such that if $\|\mathbf{F}\|_X \ge k$ then

$$
\langle P(\mathbf{F}), \mathbf{F} \rangle_X > 0.
$$

Finally, applying Lemma 14, one obtains existence of $\mathbf{F}^{n+1}$ such that $P(\mathbf{F}^{n+1}) = 0$, therefore ensuring existence of a solution to the truncated scheme (A.1)-(A.2).

**Lemma 15.** If $(f_{1,ij}^n, f_{2,ij}^n)_{i \in \mathcal{I}, j \in \mathcal{J}}^\top$ satisfies the estimates (3.20)–(3.21), then any solution $(f_{1,ij}^{n+1}, f_{2,ij}^{n+1})_{i \in \mathcal{I}, j \in \mathcal{J}}^\top$ to the truncated scheme (A.1)–(A.2) also satisfies these estimates.

Proof. Let us focus on proving that any solution to the nonlinear truncated scheme (A.1)–
(A.2) satisfies

$$
f_{1,ij}^{n+1} \geq \rho_m \chi_{1,j}.
$$

We start by setting $(i, j) \in \mathcal{I} \times \mathcal{J}$ such that

$$
f_{1,ij}^{n+1} - \rho_m \chi_{1,j} = \min_{(k,l) \in \mathcal{I} \times \mathcal{J}} (f_{1,kl}^{n+1} - \rho_m \chi_{1,l}). \quad (A.8)
$$

Our aim is now to show that this quantity is nonnegative. To this end, we multiply the
equation in (A.1) corresponding to the fixed pair (i, j) by

$$
(f_{1,ij}^{n+1} - \rho_m \chi_{1,j})^{-} = \min(0, f_{1,ij}^{n+1} - \rho_m \chi_{1,j}) \le 0. \quad (A.9)
$$

It yields an expression of the form $B_1 = B_2 + B_3$ where

$$
B_1 = \Delta x \Delta v (f_{1,ij}^{n+1} - f_{1,ij}^n) (f_{1,ij}^{n+1} - \rho_m \chi_{1,j})^{-},
$$

$$
B_2 = -\Delta t (\mathcal{F}_{i+\frac{1}{2},j}^{n+1} - \mathcal{F}_{i-\frac{1}{2},j}^{n+1}) (f_{1,ij}^{n+1} - \rho_m \chi_{1,j})^{-},
$$

$$
B_3 = \Delta t \Delta x \Delta v (\chi_{1,j} - \bar{\rho}_{2,i} \bar{f}_{1,ij}^{n+1}) (f_{1,ij}^{n+1} - \rho_m \chi_{1,j})^{-}.
$$

Starting with $B_1$, we add and subtract $\rho_m \chi_{1,j}$ to obtain

$$
B_1 = \Delta x \Delta v \left( (f_{1,ij}^{n+1} - \rho_m \chi_{1,j}) - (f_{1,ij}^n - \rho_m \chi_{1,j}) \right) (f_{1,ij}^{n+1} - \rho_m \chi_{1,j})^{-}
$$

Then, under the hypothesis that $f_{1,ij}^n$ satisfies the bounds (3.20), and using the definition (A.9) of the negative part,

$$
B_1 \geq \Delta x \Delta v (f_{1,ij}^{n+1} - \rho_m \chi_{1,j}) (f_{1,ij}^{n+1} - \rho_m \chi_{1,j})^{-} \geq 0. \quad (A.10)
$$

Next, we turn our attention to $B_3$. We need to consider two cases:
* If $f_{1,ij}^{n+1} \ge \rho_m \chi_{1,j}$ then $B_3 = 0$.

* Else, by definition of the truncated quantities (A.3) and (A.4), one has $\bar{f}_{1,ij}^{n+1} = \rho_m \chi_{1,j}$ and since the discrete velocity profiles have unit mass, $\bar{\rho}_{2,i}^{n+1} \le \rho_m^{-1}$. Therefore,

$$ \chi_{1,j} - \bar{\rho}_{2,i}^{n+1} \bar{f}_{1,ij}^{n+1} \ge \chi_{1,j} - \frac{\rho_m}{\rho_m} \chi_{1,j} = 0. $$

As a result, one gets

$$ B_3 \le 0. \qquad (A.11) $$

Finally, the sign of $B_2$ relies on the monotonicity of the numerical flux. More precisely, let us consider a monotone numerical flux in a general two-point approximation form:

$$ \mathcal{F}_{i+\frac{1}{2},j}^{n+1} = \varphi_j(f_{1,ij}^{n+1}, f_{1,i+1,j}^{n+1}), $$

where $a \mapsto \varphi_j(a, \cdot)$ is assumed to be nondecreasing and $b \mapsto \varphi_j(\cdot, b)$ is assumed to be nonincreasing. Then, the balance of fluxes at cell $\mathcal{X}_i$ rewrites as

$$ \begin{aligned} \mathcal{F}_{i+\frac{1}{2},j}^{n+1} - \mathcal{F}_{i-\frac{1}{2},j}^{n+1} = & \varphi_j(f_{1,ij}^{n+1}, f_{1,i+1,j}^{n+1}) - \varphi_j(f_{1,ij}^{n+1}, f_{1,i,j}^{n+1}) \\ & + \varphi_j(f_{1,ij}^{n+1}, f_{1,ij}^{n+1}) - \varphi_j(f_{1,i-1,j}^{n+1}, f_{1,ij}^{n+1}). \end{aligned} \qquad (A.12) $$

From our choice (A.8) of the pair $(i, j)$, we deduce that

$$ f_{1,ij}^{n+1} - \rho_m \chi_{1,j} \le f_{1,kl}^{n+1} - \rho_m \chi_{1,l}, \quad \forall (k,l) \in \mathcal{I} \times \mathcal{J}, $$

therefore in particular, for every $k \in \mathcal{I}$, $f_{1,ij}^{n+1} \le f_{1,kj}^{n+1}$. Consequently, due to the monotonicity of the function $\varphi_j$, it yields

$$ B_2 \le 0. \qquad (A.13) $$

Gathering (A.10), (A.11) and (A.13) into $B_1 = B_2 + B_3$, and using the monotonicity of the upwind fluxes, we get that $B_1 \le 0$, and consequently that

$$ f_{1,i,j}^{n+1} - \rho_m \chi_{1,j} \ge 0. $$

Lastly, since we chose $(i, j)$ satisfying (A.8),

$$ f_{1,i,j}^{n+1} \ge \rho_m \chi_{1,j}, \quad \forall (i,j) \in \mathcal{I} \times \mathcal{J}. $$

The remaining bounds can then be obtained following the same steps. $\square$

*Proof of Theorem 2.* We proceed by induction on $n$. The case $n = 0$ is satisfied by assumption. Suppose now that there exists $(f_1^n, f_2^n)^\top$ satisfying the bounds (3.20) and (3.21). Then, we showed the existence of $(f_1^{n+1}, f_2^{n+1})^\top$ solution to the truncated scheme (A.1)-(A.2), which satisfies (3.20) and (3.21) according to Lemma 15. But then, we have

$$ \begin{aligned} \bar{f}_1^{n+1} &= f_1^{n+1}, & \bar{f}_2^{n+1} &= f_2^{n+1}, \\ \bar{\rho}_1^{n+1} &= \rho_1^{n+1}, & \bar{\rho}_2^{n+1} &= \rho_2^{n+1}, \end{aligned} $$

meaning that $(f_1^{n+1}, f_2^{n+1})^\top$ is indeed a solution to the original nonlinear scheme (3.10)-(3.10), which concludes the proof. $\square$