# Discrete stability estimates for the pressureless Euler-Poisson-Boltzmann equations in the Quasi-Neutral limit

Mehdi Badsi, Nicolas Crouseilles

► To cite this version:

Mehdi Badsi, Nicolas Crouseilles. Discrete stability estimates for the pressureless Euler-Poisson-Boltzmann equations in the Quasi-Neutral limit. 2025. hal-05014885

HAL Id: hal-05014885

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

Preprint submitted on 1 Apr 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 stability estimates for the pressureless Euler-Poisson-Boltzmann equations in the Quasi-Neutral limit

Mehdi Badsi*

Nicolas Crouseilles‡

## Abstract

We propose and study a fully implicit finite volume scheme for the pressureless Euler-Poisson-Boltzmann equations on the one dimensional torus. Especially, we design a consistent and dissipative discretization of the force term which yields an unconditional energy decay. In addition, we establish a discrete analogue of the modulated energy estimate around constant states with a small velocity. Numerical experiments are carried to illustrate our theoretical results and to assess the accuracy of our scheme. A test case of the literature is also illustrated.

**keywords:** pressureless Euler-Poisson-Boltzmann, finite-volume, implicit-scheme, stability estimates, modulated energy, plasma, quasi-neutral limit

**MSCcodes:** 65M12, 82D10, 65J15.

## 1 Introduction

We consider a simplified model of a uni-dimensional plasma in which ions are cold and electrons have reached a thermodynamical equilibrium. The macroscopic density of the electrons is thus assumed to obey the Maxwell-Boltzmann law [17]. We model this plasma using a fluid approach where at time $t \in \mathbb{R}^+$ and at position $x \in \mathbb{T} = \mathbb{R}/\mathbb{Z}$, the unknowns are $\rho_{\varepsilon}(t, x) \ge 0$, $u_{\varepsilon}(t, x) \in \mathbb{R}$, $-\phi_{\varepsilon}(t, x) \in \mathbb{R}$ which stand respectively for the ions density, the ions mean velocity and the electrostatic potential. In dimensionless unit, they are assumed to satisfy the pressureless Euler-Poisson-Boltzmann equations posed on $(0, T] \times \mathbb{T}$:

$$ egin{cases} 
\partial_t \rho_{\varepsilon} + \partial_x (\rho_{\varepsilon} u_{\varepsilon}) = 0, \\
\partial_t (\rho_{\varepsilon} u_{\varepsilon}) + \partial_x (\rho_{\varepsilon} u_{\varepsilon}^2) = \rho_{\varepsilon} \partial_x \phi_{\varepsilon}, \\
\varepsilon^2 \partial_{xx} \phi_{\varepsilon} + e^{-\phi_{\varepsilon}} = \rho_{\varepsilon},
\end{cases} 
\quad (1) $$

where $T > 0$ is a time horizon and $\varepsilon > 0$ is a physical parameter called the Debye length. The system (1) is supplemented with an initial condition

$$ \rho_{\varepsilon}(0, x) = \rho_{\varepsilon}^{\text{ini}}(x), \quad u_{\varepsilon}(0, x) = u_{\varepsilon}^{\text{ini}}(x), \quad (2) $$

where $\rho_{\varepsilon}^{\text{ini}}: x \in \mathbb{T} \to \mathbb{R}^+$ and $u_{\varepsilon}^{\text{ini}}: x \in \mathbb{T} \to \mathbb{R}$ are given functions. The existence of local in time strong solutions to (1) has been rigorously established in [16] in the case of the whole space. The proof can be adapted to the case of the torus. More precisely, provided $(\rho_{\varepsilon}^{\text{ini}}, u_{\varepsilon}^{\text{ini}}) \in H^s(\mathbb{T}) \times H^{s+1}(\mathbb{T})$ with $s > \frac{1}{2}$ and $\text{ess inf}_{x \in \mathbb{T}} \rho_{\varepsilon}^{\text{ini}} > 0$ there exists $T_{\varepsilon} > 0$ and a unique strong solution $(\rho_{\varepsilon}, u_{\varepsilon}, \phi_{\varepsilon})$ with the regularity $(\rho_{\varepsilon}, u_{\varepsilon}) \in C^0([0, T_{\varepsilon}]; H^s(\mathbb{T}) \times H^{s+1}(\mathbb{T})) \cap C^1([0, T_{\varepsilon}]; H^{s-1}(\mathbb{T}) \times H^s(\mathbb{T}))$ and $\phi_{\varepsilon} \in C^0([0, T_{\varepsilon}]; H^{s+2}(\mathbb{T})) \cap C^1([0, T_{\varepsilon}]; H^{s+1}(\mathbb{T}))$ and such that $\text{ess inf}_{x \in \mathbb{T}} \rho_{\varepsilon}(t, x) > 0$ for $t \in [0, T_{\varepsilon}]$. The study of the quasi-neutral limit

<footnote>*Nantes Université, Inria centre de l'université de Rennes (Mingus team) and Laboratoire de Mathématiques Jean Leray UMR CNRS 6310, 2 Chemin de la Houssinière BP 92208, 44322 Nantes Cedex 3, France. (mehdi.badsi@univ-nantes.fr).</footnote>

<footnote>†Univ Rennes, Inria centre de l'université de Rennes (Mingus team) and IRMAR UMR CNRS 6625, F-35042 Rennes, France & ENS Rennes, France. (nicolas.crouseilles@inria.fr).</footnote>
$\epsilon \to 0$ for the strong solutions has been carried by Pu and Guo in [19, 20]. Note that in the case of the Euler-Poisson-Boltzmann equations, the quasi-neutral limit has been studied by Cordier and Grenier in [5]. When $\epsilon \to 0$ we formally expect the solution $(\rho_{\epsilon}, u_{\epsilon}, \phi_{\epsilon})$ to converge towards a solution to the isothermal Euler equations:

$$ \left\{ \begin{array}{l} \partial_t \rho_0 + \partial_x (\rho_0 u_0) = 0, \\ \partial_t (\rho_0 u_0) + \partial_x (\rho_0 u_0^2 + e^{-\phi_0}) = 0, \\ \rho_0 = e^{-\phi_0}. \end{array} \right. \qquad (3) $$

The system (3) is hyperbolic symmetrizable. Thus, the existence of local in time strong solutions is an application of the Kato-Lax-Friedrichs theory for symmetric hyperbolic systems [15, 21]. One difficulty in the study of the quasi-neutral limit (see [20]) consists in establishing thanks to high order energy estimates that the local time of existence $T_{\epsilon}$ for (1) does not shrink when $\epsilon \to 0$ (that is $\liminf_{\epsilon \to 0} T_{\epsilon} > 0$) so that there exists a time $T > 0$ on which both system (1) and (3) live.

As far as the numerical approximation of solutions to (1) is concerned, several works propose finite-difference or finite-volume schemes for plasma fluid models, which are linearly stable in the asymptotic $\epsilon \to 0$ and formally converge as $\epsilon \to 0$ see for instance [7, 1, 6]. Precisely on the pressureless Euler-Poisson-Boltzmann equations (1), a recent work is [2] where the authors study a semi-implicit finite volume scheme based on the so-called staggered discretization studied in [14, 9] for the compressible Euler equations and the compressible Stokes equations. The authors prove a discrete energy estimate for their scheme using a stabilization term. This stabilization term formally vanishes when the time step of the scheme tends to zero but is formally inconsistent when the time step is fixed and the mesh size tends to zero. A proof of convergence (up to a subsequence) of the discrete scheme in the limit $\epsilon \to 0$ is given thanks to a finite dimensional argument which, in fine, boils down to an application of the Bolzano-Weierstrass theorem. The question of the convergence rate when $\epsilon \to 0$ and its dependance with respect to the dimension of the problem is so far an open problem.

Our main focus in this work is twofold: firstly, we prove discrete stability estimates for a fully-implicit finite-volume scheme for (1) which as a by product enables us to prove existence for the scheme. Our formalism also uses a staggered discretization which enables us to establish a kinetic energy balance somehow similar to Lemma 3.1 in [14]. What is new with respect to the existing literature, is a consistent space discretization of the force term in the momentum equation of (1) which is compatible with the discrete continuity equation and leads to an unconditional energy decay. Secondly, we propose a discrete analogue of the modulated energy approach [12, 3, 19, 20] to establish non linear stability for the constant solutions of the system (1). Especially, it provides a discrete quantitative stability estimate even when $\epsilon \to 0$ for well-prepared initial data. The consistency analysis (at fixed $\epsilon$) of our scheme based on standard assumptions of the literature is to a certain extent classical, and is as a matter of fact, omitted.

The plan of this work is as follows. In Section (1.1) we establish the conservation properties of (1). In Section (1.2) we recall the modulated energy estimates. Section (2.1) defines the numerical scheme to approximate the solutions to (1). We establish the discrete energy estimates in (2.2) and prove the existence of the scheme. Then in section (2.3), we establish the stability of constant states with a small velocity using the discrete modulated energy. Eventually, we illustrate our results and discuss the numerical accuracy of our scheme in Section (3).

## 1.1 Conservation properties

In the sequel we consider strong solutions to (1) and (3) which are both defined on $[0, T]$ and such that their respective density is a positive function on $[0, T] \times \mathbb{T}$. We first establish the conservations of the system (1).

**Proposition 1.** (Conservations) Let $(\rho_{\epsilon}, u_{\epsilon}, \phi_{\epsilon})$ be a strong solution to (1) on $[0, T]$ with $\rho_{\epsilon} > 0$. Then we have for $t \in [0, T]$:

$$ \int_{\mathbb{T}} \rho_{\epsilon}(t, x)dx = \int_{\mathbb{T}} \rho_{\epsilon}(0, x)dx, \qquad (4) $$

$$ \int_{\mathbb{T}} (\rho_{\epsilon} u_{\epsilon})(t, x)dx = \int_{\mathbb{T}} (\rho_{\epsilon} u_{\epsilon})(0, x)dx, \qquad (5) $$

$$ \mathcal{H}(t) = \mathcal{H}(0), \qquad (6) $$
where $\mathcal{H}$ is the total energy given by

$$ \mathcal{H}(t) = \int_{\mathbb{T}} \rho_{\varepsilon} \frac{u_{\varepsilon}^2}{2} dx + \int_{\mathbb{T}} \frac{\varepsilon^2}{2} |\partial_x \phi_{\varepsilon}|^2 dx + \int_{\mathbb{T}} h(\phi_{\varepsilon}) dx, \quad t \in [0, T], \qquad (7) $$

where

$$ h : s \in \mathbb{R} \mapsto -(s+1)e^{-s}. \qquad (8) $$

*Proof.* The conservation of mass (4) is readily obtained by integration in space of the continuity equation using the periodicity. As for the total momentum conservation (5), we integrate in space the momentum equation and use the periodicity. It yields using the Poisson equation

$$ \int_{\mathbb{T}} \partial_t (\rho_{\varepsilon} u_{\varepsilon})(t, x) dx = \int_{\mathbb{T}} \rho_{\varepsilon} \partial_x \phi_{\varepsilon} = \int_{\mathbb{T}} (\varepsilon^2 \partial_{xx} \phi_{\varepsilon} + e^{-\phi_{\varepsilon}}) \partial_x \phi_{\varepsilon} dx = \int_{\mathbb{T}} \partial_x \left( \frac{\varepsilon^2}{2} |\partial_x \phi_{\varepsilon}|^2 - e^{-\phi_{\varepsilon}} \right) dx = 0. $$

So we get $\frac{d}{dt} \int_{\mathbb{T}} (\rho_{\varepsilon} u_{\varepsilon})(t, x) dx = 0$ for $t \in [0, T]$ and thus (5). We prove the energy conservation (6). We multiply the momentum equation by $u_{\varepsilon}$ to get $u_{\varepsilon} \cdot (\partial_t(\rho_{\varepsilon}u_{\varepsilon}) + \partial_x(\rho_{\varepsilon}u_{\varepsilon}^2)) = \rho_{\varepsilon}u_{\varepsilon} \cdot \partial_x \phi_{\varepsilon}$. Then we re-write the first term as a total derivative plus a residual term. We have,

$$ \partial_t(\rho_{\varepsilon} u_{\varepsilon}^2) + \partial_x(\rho_{\varepsilon} u_{\varepsilon}^3) = \rho_{\varepsilon} u_{\varepsilon} (\partial_t u_{\varepsilon} + \partial_x \frac{u_{\varepsilon}^2}{2}) + \rho_{\varepsilon} u_{\varepsilon} \partial_x \phi_{\varepsilon}. \qquad (9) $$

Thanks to the momentum and the continuity equation, we have $\partial_t u_{\varepsilon} + \partial_x \frac{u_{\varepsilon}^2}{2} = \partial_x \phi_{\varepsilon}$. Plugging this relation in (9) we get

$$ \partial_t(\rho_{\varepsilon} \frac{u_{\varepsilon}^2}{2}) + \partial_x(\rho_{\varepsilon} \frac{u_{\varepsilon}^3}{2}) = \rho_{\varepsilon} u_{\varepsilon} \partial_x \phi_{\varepsilon}. \qquad (10) $$

We then integrate in space (10) and use the periodicity to get

$$ \frac{d}{dt} \int_{\mathbb{T}} \rho_{\varepsilon} \frac{u_{\varepsilon}^2}{2} = \int_{\mathbb{T}} \rho_{\varepsilon} u_{\varepsilon} \partial_x \phi_{\varepsilon} dx = - \int_{\mathbb{T}} \partial_x (\rho_{\varepsilon} u_{\varepsilon}) \phi_{\varepsilon} dx = \int_{\mathbb{T}} (\partial_t \rho_{\varepsilon}) \phi_{\varepsilon} dx \qquad (11) $$

where we used the continuity equation for the last equality. Besides, the Poisson equation gives

$$ \int_{\mathbb{T}} (\partial_t \rho_{\varepsilon}) \phi_{\varepsilon} dx = \int_{\mathbb{T}} [\varepsilon^2 \partial_t (\partial_{xx} \phi_{\varepsilon}) \phi_{\varepsilon} + \partial_t (e^{-\phi_{\varepsilon}}) \phi_{\varepsilon}] dx = \int_{\mathbb{T}} [\varepsilon^2 \partial_x (\partial_t \partial_x \phi_{\varepsilon}) \phi_{\varepsilon} - (\partial_t \phi_{\varepsilon}) e^{-\phi_{\varepsilon}} \phi_{\varepsilon}] dx. \qquad (12) $$

Using an integration by parts for the first term and the definition of the function $h$ we eventually obtain

$$ \int_{\mathbb{T}} (\partial_t \rho_{\varepsilon}) \phi_{\varepsilon} dx = - \frac{d}{dt} \int_{\mathbb{T}} \frac{\varepsilon^2}{2} |\partial_x \phi_{\varepsilon}|^2 dx - \frac{d}{dt} \int_{\mathbb{T}} h(\phi_{\varepsilon}) dx. \qquad (13) $$

Gathering (13) with (11) yields $\frac{d}{dt} \mathcal{H}(t) = 0$ for $t \in [0, T]$ and thus (6). $\square$

## 1.2 The modulated energy estimate

Following [13], for a strong solution $(\rho_{\varepsilon}, u_{\varepsilon}, \phi_{\varepsilon})$ to (1) and a strong solution $(\rho_0, u_0, \phi_0)$ to (3) both defined on $[0, T]$, we define the modulated energy around $(\rho_0, u_0, \phi_0)$ at time $t \in [0, T]$ by

$$ \mathcal{E}(t) := \int_{\mathbb{T}} \rho_{\varepsilon} \frac{(u_{\varepsilon} - u_0)^2}{2} dx + \int_{\mathbb{T}} \frac{\varepsilon^2}{2} |\partial_x \phi_{\varepsilon}|^2 dx + \int_{\mathbb{T}} \left[ \tilde{h}(e^{-\phi_{\varepsilon}}) - \left( \tilde{h}(\rho_0) + \tilde{h}'(\rho_0)(e^{-\phi_{\varepsilon}} - \rho_0) \right) \right] dx, \quad (14) $$

where $\tilde{h} : \mathbb{R}_{\star}^{+} \to \mathbb{R}$ is the Boltzmann entropy function given by

$$ \forall \psi > 0, \tilde{h}(\psi) := h(-\log(\psi)) = \psi \log(\psi) - \psi. \qquad (15) $$

We point out that $\mathcal{E}$ is a non negative functional since it is the sum of three non negative functionals. The fact that the last term is non negative is due to the fact that the function $\tilde{h}$ is convex. Provided $\|\rho_{\varepsilon}\|_{L_{t,x}^{\infty}}$
and $\|\frac{1}{\rho_{\varepsilon}}\|_{L^\infty_{t,x}}$ are uniformly bounded in $\varepsilon$, we have $\mathcal{E}(t) \gtrsim \|u_{\varepsilon}(t)-u_0(t)\|_{L^2(\mathbb{T})}^2+\|e^{-\phi_{\varepsilon}(t)}-\rho_0(t)\|_{L^2(\mathbb{T})}$ where the constant in the inequality is independent of $\varepsilon$. Thus $\mathcal{E}(t)$ yields a control of the distance in $L^2(\mathbb{T})$ between $(\rho_{\varepsilon}, u_{\varepsilon}, \phi_{\varepsilon})$ and $(\rho_0, u_0, \phi_0)$ at time $t$. Expanding the terms of (14), we have the decomposition for $t \in [0, T]$:

$$ \mathcal{E}(t) = \mathcal{H}(t) + \mathcal{E}_{\text{kin}}(t) - \mathcal{E}_{\text{int}}(t), \qquad (16) $$

with $\mathcal{H}(t)$ given by (7) and

$$ \mathcal{E}_{\text{kin}}(t) = \int_{\mathbb{T}} \rho_{\varepsilon} \left( \frac{u_0^2}{2} - u_{\varepsilon} u_0 \right) dx, \qquad (17) $$

$$ \mathcal{E}_{\text{int}}(t) = \int_{\mathbb{T}} (\tilde{h}(\rho_0) + \tilde{h}'(\rho_0)(e^{-\phi_{\varepsilon}} - \rho_0)) dx. \qquad (18) $$

A simple identity that will be used in the modulated energy estimate is the following.

**Lemma 1.** Let $(\rho_{\varepsilon}, u_{\varepsilon}, \phi_{\varepsilon})$ be a strong solution to (1) on $[0, T]$ with $\rho_{\varepsilon} > 0$. Then for every $\psi \in C^1([0, T] \times \mathbb{T})$:

$$ \rho_{\varepsilon}(\partial_t + u_{\varepsilon}\partial_x)(\psi) = \partial_t(\rho_{\varepsilon}\psi) + \partial_x(\rho_{\varepsilon}u_{\varepsilon}\psi). \qquad (19) $$

*Proof.* A direct computation yields

$$ \rho_{\varepsilon}(\partial_t + u_{\varepsilon}\partial_x)(\psi) = \rho_{\varepsilon}\partial_t\psi + \rho_{\varepsilon}u_{\varepsilon}\partial_x\psi = \partial_t(\rho_{\varepsilon}\psi) + \partial_x(\rho_{\varepsilon}u_{\varepsilon}\psi) - \psi(\partial_t\rho_{\varepsilon} + \partial_x(\rho_{\varepsilon}u_{\varepsilon})). $$

It yields the claim thanks to the continuity equation.

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

The main quantitative stability estimate is stated in the following proposition.

**Proposition 2.** (Modulated energy estimate) Let $(\rho_{\varepsilon}, u_{\varepsilon}, \phi_{\varepsilon})$ be a strong solution to (1) on $[0, T]$ with $\rho_{\varepsilon} > 0$ and $(\rho_0, u_0, \phi_0)$ be a strong solution to (3) on $[0, T]$ with $\rho_0 > 0$. Then we have for $t \in [0, T]$,

$$ \mathcal{E}(t) \le \mathcal{E}(0)e^{2\|\partial_x u_0\|_{L^\infty_{t,x}} t} + \varepsilon^2 \int_0^t \mathcal{L}(\tau)e^{2(t-\tau)\|\partial_x u_0\|_{L^\infty_{t,x}}} d\tau. \qquad (20) $$

$$ \text{where } \mathcal{L}(\tau) = \int_{\mathbb{T}} \partial_{xx}\phi_{\varepsilon} \log \rho_0 dx - \int_{\mathbb{T}} \partial_{xx}\phi_{\varepsilon} u_0 \partial_x \log \rho_0 dx. $$

*Proof.* Thanks to the energy conservation (6) and (16), we have for $t \in [0, T]$,

$$ \frac{d}{dt} \mathcal{E}(t) = \frac{d}{dt} \mathcal{E}_{\text{kin}}(t) - \frac{d}{dt} \mathcal{E}_{\text{int}}(t). \qquad (21) $$

We shall now estimate each term. For the first term we have,

$$ \frac{d}{dt} \mathcal{E}_{\text{kin}}(t) = \int_{\mathbb{T}} \partial_t \rho_{\varepsilon} \left( \frac{u_0^2}{2} - u_{\varepsilon} u_0 \right) dx + \int_{\mathbb{T}} \rho_{\varepsilon} \partial_t \left( \frac{u_0^2}{2} - u_{\varepsilon} u_0 \right) dx. $$

Using the continuity equation $\partial_t \rho_{\varepsilon} + \partial_x(\rho_{\varepsilon}u_{\varepsilon}) = 0$ and an integration by parts we obtain

$$ \frac{d}{dt} \mathcal{E}_{\text{kin}}(t) = \int_{\mathbb{T}} \rho_{\varepsilon} u_{\varepsilon} \partial_x \left( \frac{u_0^2}{2} - u_{\varepsilon} u_0 \right) dx + \int_{\mathbb{T}} \rho_{\varepsilon} \partial_t \left( \frac{u_0^2}{2} - u_{\varepsilon} u_0 \right) dx = \int_{\mathbb{T}} \rho_{\varepsilon} (\partial_t + u_{\varepsilon} \partial_x) \left( \frac{u_0^2}{2} - u_{\varepsilon} u_0 \right) dx. $$

We now apply Lemma (1) with the function $\psi = \frac{u_0^2}{2} - u_{\varepsilon}u_0$. Thus,

$$ \frac{d}{dt} \mathcal{E}_{\text{kin}}(t) = \int_{\mathbb{T}} \partial_t \left( \rho_{\varepsilon} \left( \frac{u_0^2}{2} - u_{\varepsilon} u_0 \right) \right) + \partial_x \left( \rho_{\varepsilon} u_{\varepsilon} \left( \frac{u_0^2}{2} - u_{\varepsilon} u_0 \right) \right) dx. $$

We set for ease $I_1(t) = \int_{\mathbb{T}} \partial_t (\rho_{\varepsilon} \frac{u_0^2}{2}) + \partial_x (\rho_{\varepsilon} u_{\varepsilon} \frac{u_0^2}{2}) dx$, and $I_2(t) = - \int_{\mathbb{T}} \partial_t (\rho_{\varepsilon} u_{\varepsilon} u_0) + \partial_x (\rho_{\varepsilon} u_{\varepsilon}^2 u_0) dx$. We begin to treat $I_1$. Multiplying the continuity equation of (1) by $\frac{u_0^2}{2}$ we obtain

$$ \partial_t (\rho_{\varepsilon} \frac{u_0^2}{2}) + \partial_x (\rho_{\varepsilon} u_{\varepsilon} \frac{u_0^2}{2}) = \rho_{\varepsilon} (\partial_t + u_{\varepsilon} \partial_x) \left( \frac{u_0^2}{2} \right). $$
Therefore,

$$I_1(t) = \int_{\mathbb{T}} \rho_{\varepsilon}(\partial_t + u_{\varepsilon}\partial_x)\left(\frac{u_0^2}{2}\right)dx = \int_{\mathbb{T}} \rho_{\varepsilon}u_0(\partial_t + u_{\varepsilon}\partial_x)(u_0)dx$$

As for $I_2$, we multiply the momentum equation of (1) by $u_0$ to get

$$\partial_t(\rho_{\varepsilon}u_{\varepsilon}u_0) + \partial_x(\rho_{\varepsilon}u_{\varepsilon}^2u_0) = \rho_{\varepsilon}u_0\partial_x\phi_{\varepsilon} + \rho_{\varepsilon}u_{\varepsilon}(\partial_t + u_{\varepsilon}\partial_x)(u_0).$$

Therefore,

$$I_2(t) = - \int_{\mathbb{T}} [\rho_{\varepsilon}u_0\partial_x\phi_{\varepsilon} + \rho_{\varepsilon}u_{\varepsilon}(\partial_t + u_{\varepsilon}\partial_x)(u_0)]dx.$$

Combining $I_1$ and $I_2$ we glean,

$$\begin{aligned}
\frac{d}{dt}\mathcal{E}_{kin}(t) &= \int_{\mathbb{T}} [\rho_{\varepsilon}u_0(\partial_t + u_{\varepsilon}\partial_x)(u_0) - \rho_{\varepsilon}u_0\partial_x\phi_{\varepsilon} - \rho_{\varepsilon}u_{\varepsilon}(\partial_t + u_{\varepsilon}\partial_x)(u_0)]dx \\
&= \int_{\mathbb{T}} [\rho_{\varepsilon}(u_0 - u_{\varepsilon})(\partial_t + u_{\varepsilon}\partial_x)(u_0) - \rho_{\varepsilon}u_0\partial_x\phi_{\varepsilon}]dx \\
&= \int_{\mathbb{T}} \rho_{\varepsilon}(u_0 - u_{\varepsilon})(\partial_t + u_0\partial_x)(u_0)dx - \int_{\mathbb{T}} \rho_{\varepsilon}(u_0 - u_{\varepsilon})^2\partial_x u_0 dx - \int_{\mathbb{T}} \rho_{\varepsilon}\partial_x\phi_{\varepsilon}u_0dx.
\end{aligned}$$

Using the Poisson equation, we have for the last term

$$\begin{aligned}
\int_{\mathbb{T}} \rho_{\varepsilon}\partial_x\phi_{\varepsilon}u_0dx &= \int_{\mathbb{T}} (\varepsilon^2\partial_{xx}\phi_{\varepsilon} + e^{-\phi_{\varepsilon}})\partial_x\phi_{\varepsilon}u_0dx = \int_{\mathbb{T}} \partial_x\left(\frac{\varepsilon^2}{2}|\partial_x\phi_{\varepsilon}|^2 - e^{-\phi_{\varepsilon}}\right)u_0dx \\
&= -\int_{\mathbb{T}} \left(\frac{\varepsilon^2}{2}|\partial_x\phi_{\varepsilon}|^2 - e^{-\phi_{\varepsilon}}\right)\partial_x u_0 dx.
\end{aligned}$$

Therefore,

$$\begin{aligned}
\frac{d}{dt}\mathcal{E}_{kin}(t) &= \int_{\mathbb{T}} \left(\frac{\varepsilon^2}{2}|\partial_x\phi_{\varepsilon}|^2 - e^{-\phi_{\varepsilon}}\right)\partial_x u_0 dx - \int_{\mathbb{T}} \rho_{\varepsilon}(u_0 - u_{\varepsilon})^2 \partial_x u_0 dx \quad (22) \\
&\quad + \int_{\mathbb{T}} \rho_{\varepsilon}(u_0 - u_{\varepsilon})(\partial_t + u_0\partial_x)(u_0)dx.
\end{aligned}$$

As for $\mathcal{E}_{int}$, we have

$$\begin{aligned}
-\frac{d}{dt}\mathcal{E}_{int}(t) &= \frac{d}{dt}\int_{\mathbb{T}}(\rho_0 - e^{-\phi_{\varepsilon}}\log \rho_0)dx \\
&= \int_{\mathbb{T}}\partial_t\rho_0dx - \int_{\mathbb{T}}\left[(\partial_t e^{-\phi_{\varepsilon}})\log \rho_0 + e^{-\phi_{\varepsilon}}\frac{\partial_t\rho_0}{\rho_0}\right]dx \\
&= \int_{\mathbb{T}}\left(-\frac{e^{-\phi_{\varepsilon}}}{\rho_0} + 1\right)\partial_t\rho_0dx + \int_{\mathbb{T}}(\varepsilon^2\partial_{txx}\phi_{\varepsilon} - \partial_t\rho_{\varepsilon})\log \rho_0dx \\
&= \int_{\mathbb{T}}\left(-\frac{e^{-\phi_{\varepsilon}}}{\rho_0} + 1\right)\partial_t\rho_0dx + \int_{\mathbb{T}}\varepsilon^2\partial_{txx}\phi_{\varepsilon}\log \rho_0dx + \int_{\mathbb{T}}\partial_x(\rho_{\varepsilon}u_{\varepsilon})\log \rho_0dx, \\
&= \int_{\mathbb{T}}\frac{e^{-\phi_{\varepsilon}}}{\rho_0}(\rho_0\partial_x u_0 + u_0\partial_x\rho_0)dx + \int_{\mathbb{T}}\varepsilon^2\partial_{txx}\phi_{\varepsilon}\log \rho_0dx + \int_{\mathbb{T}}\partial_x(\rho_{\varepsilon}u_{\varepsilon})\log \rho_0dx, \\
&= \int_{\mathbb{T}}e^{-\phi_{\varepsilon}}\partial_x u_0 dx + \int_{\mathbb{T}}e^{-\phi_{\varepsilon}}u_0\partial_x\log \rho_0 dx + \int_{\mathbb{T}}\varepsilon^2\partial_{txx}\phi_{\varepsilon}\log \rho_0 dx + \int_{\mathbb{T}}\partial_x(\rho_{\varepsilon}u_{\varepsilon})\log \rho_0 dx.
\end{aligned}$$

where we have used that $\int_{\mathbb{T}}\partial_t\rho_0dx = 0$.

We then observe that

$$\begin{aligned}
\int_{\mathbb{T}} \rho_{\varepsilon}(u_0 - u_{\varepsilon})\partial_x \log \rho_0 dx &= \int_{\mathbb{T}} \rho_{\varepsilon}u_0 \partial_x \log \rho_0 dx - \int_{\mathbb{T}} (\rho_{\varepsilon}u_{\varepsilon})\partial_x \log \rho_0 dx \\
&= \int_{\mathbb{T}} e^{-\phi_{\varepsilon}} u_0 \partial_x \log \rho_0 dx + \varepsilon^2 \int_{\mathbb{T}} \partial_{xx}\phi_{\varepsilon} u_0 \partial_x \log \rho_0 dx + \int_{\mathbb{T}} \partial_x (\rho_{\varepsilon}u_{\varepsilon}) \log \rho_0 dx,
\end{aligned}$$
so that (23) rewrites

$$ -\frac{d}{dt} \mathcal{E}_{\text{int}}(t) = \int_{\mathbb{T}} e^{-\phi_{\varepsilon}} \partial_x u_0 dx + \int_{\mathbb{T}} \rho_{\varepsilon}(u_0 - u_{\varepsilon}) \partial_x \log \rho_0 dx \\ + \int_{\mathbb{T}} \varepsilon^2 \partial_{txx} \phi_{\varepsilon} \log \rho_0 dx - \varepsilon^2 \int_{\mathbb{T}} \partial_{xx} \phi_{\varepsilon} u_0 \partial_x \log \rho_0 dx. \quad (24) $$

Gathering the equalities (22) and (24), we get

$$ \frac{d}{dt} \mathcal{E}(t) = \int_{\mathbb{T}} \left( \frac{\varepsilon^2}{2} |\partial_x \phi_{\varepsilon}|^2 - e^{-\phi_{\varepsilon}} \right) \partial_x u_0 dx - \int_{\mathbb{T}} \rho_{\varepsilon} (u_0 - u_{\varepsilon})^2 \partial_x u_0 dx \quad (25) $$

$$ + \int_{\mathbb{T}} \rho_{\varepsilon}(u_0 - u_{\varepsilon})(\partial_t + u_0 \partial_x)(u_0) dx \quad (26) $$

$$ + \int_{\mathbb{T}} e^{-\phi_{\varepsilon}} \partial_x u_0 dx + \int_{\mathbb{T}} \rho_{\varepsilon}(u_0 - u_{\varepsilon}) \partial_x \log \rho_0 dx \quad (27) $$

$$ + \int_{\mathbb{T}} \varepsilon^2 \partial_{txx} \phi_{\varepsilon} \log \rho_0 dx - \varepsilon^2 \int_{\mathbb{T}} \partial_{xx} \phi_{\varepsilon} u_0 \partial_x \log \rho_0 dx \quad (28) $$

$$ = \int_{\mathbb{T}} \rho_{\varepsilon}(u_0 - u_{\varepsilon}) \left[ \partial_t u_0 + u_0 \partial_x u_0 + \partial_x \log \rho_0 \right] dx - \int_{\mathbb{T}} \rho_{\varepsilon}(u_0 - u_{\varepsilon})^2 \partial_x u_0 dx \quad (29) $$

$$ + \int_{\mathbb{T}} \frac{\varepsilon^2}{2} |\partial_x \phi_{\varepsilon}|^2 \partial_x u_0 dx + \int_{\mathbb{T}} \varepsilon^2 \partial_{txx} \phi_{\varepsilon} \log \rho_0 dx - \varepsilon^2 \int_{\mathbb{T}} \partial_{xx} \phi_{\varepsilon} u_0 \partial_x \log \rho_0 dx. \quad (30) $$

Moreover, $u_0$ satisfies $\partial_t u_0 + u_0 \partial_x u_0 + \partial_x \log \rho_0 = 0$ so that we eventually obtain

$$ \frac{d}{dt} \mathcal{E}(t) = \int_{\mathbb{T}} \left[ -\rho_{\varepsilon}(u_0 - u_{\varepsilon})^2 + \frac{\varepsilon^2}{2} |\partial_x \phi_{\varepsilon}|^2 \right] \partial_x u_0 dx + \varepsilon^2 \int_{\mathbb{T}} \left[ \partial_{txx} \phi_{\varepsilon} \log \rho_0 - \partial_{xx} \phi_{\varepsilon} u_0 \partial_x \log \rho_0 \right] dx. $$

Hence, we deduce the following inequality

$$ \frac{d}{dt} \mathcal{E}(t) \leq \varepsilon^2 \mathcal{L}(t) + 2\mathcal{E}(t) \|\partial_x u_0\|_{L_{i,x}^{\infty}} $$

with $\mathcal{L}(t) = \int_{\mathbb{T}} [\partial_{txx} \phi_{\varepsilon} \log \rho_0 - \partial_{xx} \phi_{\varepsilon} u_0 \partial_x \log \rho_0] dx$ and a Grönwall lemma enables us to conclude. $\square$

Note that in [13], it is explained how to bound $\mathcal{L}(t)$ with respect to $\varepsilon$.

## 2 Discretization

Let $\Delta x = \frac{1}{N+1}$ where $N \in \mathbb{N}^*$ is fixed. We consider a uniform grid defined by the sequence of points $(x_i = i\Delta x)_{i \in \mathbb{Z}}$. Since we work on $\mathbb{T}$, we shall identify two points of the same grid according to the equivalence relation defined on $\mathbb{R}$ by

$$ \forall x, y \in \mathbb{R}, \quad x \equiv y \mod \mathbb{Z} \Leftrightarrow x - y \in \mathbb{Z}. \quad (31) $$

It yields in particular an identification of the torus $\mathbb{T}$ with the unit interval $[0, 1)$ through the unique decomposition of a real number:

$$ \forall x \in \mathbb{R}, \quad x = \lfloor x \rfloor + \{x\} \quad (32) $$

where $\lfloor \cdot \rfloor$ denotes the integer part function and $\{ \cdot \}$ denotes the fractional part function. Especially,

$$ \forall x \in \mathbb{R}, \quad x = \lfloor x \rfloor + \{x\} \mod \mathbb{Z}. \quad (33) $$

The relation (33) applied to the grid points $(i\Delta x)_{i \in \mathbb{Z}}$ yields

$$ \forall i \in \mathbb{Z}, \quad \exists! i^* \in \{0, \dots, N\} \quad x_i = x_{i^*} \mod \mathbb{Z}, \quad (34) $$

where $i^*$ is the remainder of the Euclidean division of $i$ by $(N+1)$. We then identify the quotient space $\mathbb{Z}/(N+1)\mathbb{Z}$ with the first $N+1$ non negative integers. So, in the following we shall systematically identify
an integer with its unique representation in $\{0, \dots, N\}$. To approximate the solutions of (1), we use a finite volume approach where $\mathbb{T}$ is covered by a union of non empty disjoints intervals of size $\Delta x$. We then define two meshes

$$ \mathcal{T} := \bigcup_{i=0}^{N} C_i, \quad C_i = \left[x_i - \frac{\Delta x}{2}, x_i + \frac{\Delta x}{2}\right), \quad \mathcal{T}^* := \bigcup_{i=0}^{N} C_{i+\frac{1}{2}}, \quad C_{i+\frac{1}{2}} = \left[x_{i+\frac{1}{2}} - \frac{\Delta x}{2}, x_{i+\frac{1}{2}} + \frac{\Delta x}{2}\right), \quad (35) $$

where $x_{i+\frac{1}{2}} = x_i + \frac{\Delta x}{2}$, $\mathcal{T}$ is called the primal mesh and $\mathcal{T}^*$ is called the dual mesh. We then consider two spaces of piecewise constant functions on $\mathbb{T}$:

$$ X(\mathcal{T}) = \{v \in L_{\text{loc}}^{1}(\mathbb{T}; \mathbb{R}), \forall i \in \{0, \dots, N\}, v(x) = v_i := \frac{1}{\Delta x} \int_{C_i} v(x)dx \text{ if } x \in C_i\}, \quad (36) $$

$$ X(\mathcal{T}^*) = \{v \in L_{\text{loc}}^{1}(\mathbb{T}; \mathbb{R}), \forall i \in \{0, \dots, N\}, v(x) = v_{i+\frac{1}{2}} := \frac{1}{\Delta x} \int_{C_{i+\frac{1}{2}}} v(x)dx \text{ if } x \in C_{i+\frac{1}{2}}\}. \quad (37) $$

A natural space to estimate the solution to the Poisson equation is $H^1(\mathbb{T})$ endowed with its canonical norm. We consider the discrete analogue for functions in $X(\mathcal{T})$ by introducing the discrete gradient defined by

$$ \forall \varphi \in X(\mathcal{T}), \quad (\delta\varphi) \in X(\mathcal{T}^*) \text{ and } (\delta\varphi)_{i+\frac{1}{2}} = \frac{\varphi_{i+1} - \varphi_i}{\Delta x}, \quad i \in \{0, \dots, N\}. \quad (38) $$

The discrete $H^1$ semi-norm is defined by

$$ \forall \varphi \in X(\mathcal{T}), \quad |\varphi|_{H^1(\mathbb{T})} = \left( \sum_{i=0}^{N} \left| \frac{\varphi_{i+1} - \varphi_i}{\Delta x} \right|^2 \Delta x \right)^{\frac{1}{2}}. \quad (39) $$

The discrete Laplacian is defined for functions in $X(\mathcal{T})$ by

$$ \forall \varphi \in X(\mathcal{T}), \quad (\Delta\varphi) \in X(\mathcal{T}) \text{ and } (\Delta\varphi)_i = \frac{\varphi_{i+1} - 2\varphi_i + \varphi_{i-1}}{\Delta x^2}, \quad i \in \{0, \dots, N\}. \quad (40) $$

We will use routinely a discrete analogue of the integration by parts formulas. More precisely, we have.

**Lemma 2.** (Discrete integration by parts) It holds,

$$ \forall (v, p) \in X(\mathcal{T}^*) \times X(\mathcal{T}), \quad \sum_{i=0}^{N} v_{i+1/2}(p_{i+1} - p_i)\Delta x = - \sum_{i=0}^{N} (v_{i+1/2} - v_{i-1/2})p_i\Delta x, \quad (41) $$

$$ \forall (\varphi, \psi) \in X(\mathcal{T})^2, \quad \sum_{i=0}^{N} (\Delta\varphi)_i \psi_i \Delta x = - \sum_{i=0}^{N} (\delta\varphi)_{i+\frac{1}{2}} (\delta\psi)_{i+\frac{1}{2}} \Delta x. \quad (42) $$

*Proof.* Both (41) and (42) are obtained using a translation of indices and using the periodicity. We only prove (41). Using a translation of indices, we get

$$ \begin{aligned} \sum_{i=0}^{N} v_{i+1/2}(p_{i+1} - p_i)\Delta x &= \sum_{i=0}^{N} v_{i+1/2}p_{i+1}\Delta x - \sum_{i=0}^{N} v_{i+1/2}p_i\Delta x = \sum_{i=1}^{N+1} v_{i-1/2}p_i\Delta x - \sum_{i=0}^{N} v_{i+1/2}p_i\Delta x \\ &= -\sum_{i=0}^{N} (v_{i+1/2} - v_{i-1/2})p_i\Delta x - v_{-1/2}p_0 + v_{N+1/2}p_{N+1}. \end{aligned} $$

By periodicity, we have $v_{-1/2} = v_{N+1/2}$ and $p_0 = p_{N+1}$ which gives the result. $\square$

As we are mainly concerned with discrete analogues of (6) and (20), we introduce a discrete analogue of the energy functional and of its modulated version. We thus define the discrete total energy for $(\rho, u, \phi) \in X(\mathcal{T}) \times X(\mathcal{T}^*) \times X(\mathcal{T})$ by

$$ \mathcal{H}(\rho, u, \phi) = \sum_{i=0}^{N} \left( \rho_{i+\frac{1}{2}} \frac{u_{i+\frac{1}{2}}^2}{2} + \frac{\varepsilon^2}{2} \left| \frac{\phi_{i+1} - \phi_i}{\Delta x} \right|^2 + h(\phi_i) \right) \Delta x. \quad (43) $$
The discrete modulated version around a constant state $(\bar{u}, \bar{\phi}) \in X(\mathcal{T}^*) \times X(\mathcal{T})$ is given by

$$ \begin{aligned} \mathcal{E}(\rho, u, \phi | \bar{u}, \bar{\phi}) = & \sum_{i=0}^{N} \rho_{i+\frac{1}{2}} \frac{(u_{i+\frac{1}{2}} - \bar{u})^2}{2} \Delta x + \sum_{i=0}^{N} \frac{\varepsilon^2}{2} \left| \frac{\phi_{i+1} - \phi_i}{\Delta x} \right|^2 \Delta x \\ & + \sum_{i=0}^{N} \left[ \tilde{h}(e^{-\phi_i}) - \left( \tilde{h}(e^{-\bar{\phi}}) + \tilde{h}'(e^{-\bar{\phi}})(e^{-\phi_i} - e^{-\bar{\phi}}) \right) \right] \Delta x. \end{aligned} \quad (44) $$

As in (16), expanding the first term we have the decomposition

$$ \mathcal{E}(\rho, u, \phi | \bar{u}, \bar{\phi}) = \mathcal{H}(\rho, u, \phi) + \mathcal{E}_{kin}(\rho, u | \bar{u}) - \mathcal{E}_{int}(\phi | \bar{\phi}) \quad (45) $$

where $\mathcal{H}$ is given in (43) and

$$ \mathcal{E}_{kin}(\rho, u | \bar{u}) = \sum_{i=0}^{N} \rho_{i+\frac{1}{2}} \left( \frac{\bar{u}^2}{2} - u_{i+\frac{1}{2}} \bar{u} \right) \Delta x, \quad (46) $$

$$ \mathcal{E}_{int}(\phi | \bar{\phi}) = \sum_{i=0}^{N} \left( \tilde{h}(e^{-\bar{\phi}}) + \tilde{h}'(e^{-\bar{\phi}})(e^{-\phi_i} - e^{-\bar{\phi}}) \right) \Delta x. \quad (47) $$

## 2.1 The time implicit scheme

We fix $N_T \in \mathbb{N}^*$ and set $\Delta t = \frac{T}{N_T}$. We consider a uniform in time discretization $(t_n)_{n=0,...,N_T} = (n\Delta t)_{n=0,...,N_T}$. For each $n = 0, ..., N_T$, we consider $(\rho^n, u^n, \phi^n) \in X(\mathcal{T}) \times X(\mathcal{T}^*) \times X(\mathcal{T})$ defined by induction for $n \in \{0, ..., N_T - 1\}$:

$$ \left\{ \begin{array}{l} \frac{\rho_i^{n+1} - \rho_i^n}{\Delta t} + \frac{1}{\Delta x} (\mathcal{F}_{i+1/2}^{n+1} - \mathcal{F}_{i-1/2}^{n+1}) = 0, \\ \frac{\rho_{i+1/2}^{n+1} u_{i+1/2}^{n+1} - \rho_i^n u_{i+1/2}^{n+1}}{\Delta t} + \frac{1}{\Delta x} (\mathcal{F}_{i+1}^{n+1} u_{i+1}^{n+1} - \mathcal{F}_i^{n+1} u_i^{n+1}) = \tilde{\rho}_{i+1/2}^{n+1} \delta(\phi^{n+1})_{i+1/2}, \\ \varepsilon^2 \Delta t (\phi^{n+1})_i + e^{-\phi_i^{n+1}} = \rho_i^{n+1}, \end{array} \right. \quad (48) $$

where $i \in \{0, ..., N\}$. The system (48) is supplemented with an initial condition $\rho^0 \in X(\mathcal{T})$ and $u^0 \in X(\mathcal{T}^*)$ such that

$$ \rho_i^0 = \frac{1}{\Delta x} \int_{C_i} \rho_{\varepsilon}^{\text{ini}}(x) dx, \quad u_{i+1/2}^0 = \frac{1}{\Delta x} \int_{C_{i+1/2}} u_{\varepsilon}^{\text{ini}}(x) dx. \quad (49) $$

For the sake of conciseness in the notation we have voluntarily discarded the dependence with respect to $\varepsilon$ of the discrete solution to (48). In (48), the flux of mass at the interface $x_{i+1/2}$ is defined by

$$ \mathcal{F}_{i+\frac{1}{2}}^{n+1} = G(\rho_i^{n+1}, \rho_{i+1}^{n+1}, u_{i+\frac{1}{2}}^{n+1}) \quad (50) $$

where $G: \mathbb{R}^3 \to \mathbb{R}$ is the function defined by

$$ \forall (s, t, u) \in \mathbb{R}^3, \quad G(s, t, u) = sg(u) - t(g(u) - u) \quad (51) $$

where $g: \mathbb{R} \to \mathbb{R}$ is some arbitrary function that satisfies the following assumptions:

$$ g \in \text{Lip}(\mathbb{R}), \quad (52) $$

$$ \forall u \in \mathbb{R}, \quad g(u) \geq \max(u, 0), \quad (53) $$

$$ g \text{ is differentiable at } 0. \quad (54) $$

Under the assumptions (52), (53) and (54), the function $G$ is a Lipschitz continuous, non decreasing in its first variable and non increasing in its second variable. Moreover, it verifies the usual consistency relation in the finite volume sense:

$$ \forall (\rho, u) \in \mathbb{R}^2, \quad G(\rho, \rho, u) = \rho u. \quad (55) $$
The function $g$ could be thought as some regularization of the function $u \mapsto \max(u, 0)$. The density and its flux of mass at the interface $x_{i+\frac{1}{2}}$ are defined by:

$$ \rho_{i+\frac{1}{2}}^{n+1} = \frac{\rho_{i+1}^{n+1} + \rho_i^{n+1}}{2}, \quad \mathcal{F}_i^{n+1} = \frac{\mathcal{F}_{i+\frac{1}{2}}^{n+1} + \mathcal{F}_{i-\frac{1}{2}}^{n+1}}{2}. \qquad (56) $$

An elementary consequence of these two definitions and of the discrete continuity equation is:

$$ \exists i \in \{0, \dots, N\}, \forall j \in \{i, i+1\} : \frac{\rho_j^{n+1} - \rho_j^n}{\Delta t} + \frac{1}{\Delta x}(\mathcal{F}_{j+1/2}^{n+1} - \mathcal{F}_{j-1/2}^{n+1}) = 0 \\ \implies \frac{\rho_{i+1/2}^{n+1} - \rho_{i+1/2}^n}{\Delta t} + \frac{1}{\Delta x}(\mathcal{F}_{i+1}^{n+1} - \mathcal{F}_{i-1}^{n+1}) = 0. $$

It embodies the fact that if the discrete continuity equation is satisfied on the cells $C_i$ and $C_{i+1}$ for some $i$ then it has a dual analogue on the cell $C_{i+\frac{1}{2}}$ which is in between. The velocity at the interface $x_i$ is defined by:

$$ u_i^{n+1} = \begin{cases} u_{i-\frac{1}{2}}^{n+1} & \text{if } \mathcal{F}_i^{n+1} \ge 0, \\ u_{i+\frac{1}{2}}^{n+1} & \text{if } \mathcal{F}_i^{n+1} < 0. \end{cases} \qquad (57) $$

An originality of this work is a definition of the density in the forcing term of discrete momentum equation given by:

$$ \tilde{\rho}_{i+\frac{1}{2}}^{n+1} = \begin{cases} \frac{G(\rho_i^{n+1}, \rho_{i+1}^{n+1}, u_{i+\frac{1}{2}}^{n+1}) - G(\rho_i^{n+1}, \rho_{i+1}^{n+1}, 0)}{u_{i+\frac{1}{2}}^{n+1} - 0}, & \text{if } u_{i+\frac{1}{2}}^{n+1} \neq 0, \\ \rho_{i+1}^{n+1} - (\rho_{i+1}^{n+1} - \rho_i^{n+1})g'(0), & \text{if } u_{i+\frac{1}{2}}^{n+1} = 0, \end{cases} \qquad (58) $$

which yields the unconditional dissipativity of the force term in the sense that

$$ \forall n \in \{0, \dots, N_T - 1\}, \quad \sum_{i=0}^{N} \tilde{\rho}_{i+\frac{1}{2}}^{n+1} u_{i+\frac{1}{2}}^{n+1} \delta(\phi^{n+1})_{i+\frac{1}{2}} \Delta x + \frac{\varepsilon^2}{2} |\phi^{n+1}|_{H^1(\mathbb{T})}^2 + \sum_{i=0}^{N} h(\phi_i^{n+1}) \Delta x \quad (59) \\ \le \frac{\varepsilon^2}{2} |\phi^n|^2_{H^1(\mathbb{T})} + \sum_{i=0}^{N} h(\phi_i^n) \Delta x. $$

Note that the function $g$ is differentiable at the origin, so the definition (58) makes $\tilde{\rho}_{i+\frac{1}{2}}^{n+1}$ a continuous function of each of its arguments. Indeed we prove.

**Lemma 3.** (Continuity of $\tilde{\rho}$) The function $\tilde{\rho} : \mathbb{R}^3 \to \mathbb{R}$ given by

$$ \tilde{\rho}(s, t, u) = \begin{cases} \frac{G(s, t, u) - G(s, t, 0)}{u-0} & (s, t, u) \in \mathbb{R}^3 \setminus F, \\ t - (t-s)g'(0) & (s, t, u) \in F, \end{cases} \qquad (60) $$

where $F = \{(s, t, 0) : (s, t) \in \mathbb{R}^2\}$ is continuous on $\mathbb{R}^3$.

*Proof.* Since $g$ is a Lipschitz continuous on $\mathbb{R}$, the function $G$ is continuous on $\mathbb{R}^3$ as a sum and product of continuous functions. Thus, $\tilde{\rho}$ is a continuous function in the open set $\mathbb{R}^3 \setminus F$. Then remark that for $(s, t, u) \in \mathbb{R}^3 \setminus F$ we have

$$ \tilde{\rho}(s, t, u) = t - (t - s)\hat{g}(u) \qquad (61) $$

where $\hat{g}(u) = \begin{cases} \frac{g(u) - g(0)}{u} & \text{if } u \neq 0, \\ g'(0) & \text{if } u = 0. \end{cases}$ Observe that $\hat{g}$ is a continuous function on $\mathbb{R}$ since $g$ is continuous on $\mathbb{R}^*$ and it is differentiable at the origin. For $(s, t, u) \in \mathbb{R}^3 \setminus F$ we have

$$ |\tilde{\rho}(s, t, u) - (t - (t - s))g'(0)| = |t - s| |\hat{g}(u) - g'(0)|. \qquad (62) $$

By continuity of $\hat{g}$ we deduce that $\tilde{\rho}(s, t, u) \to \tilde{\rho}(s, t, 0)$ as $(s, t, u) \to (s, t, 0)$. $\square$
The previous lemma is important in view of the existence theory for (48) since we shall invoke the Brouwer fixed-point theorem which requires the continuity of an appropriate map. We mention that the main idea behind (58) is the need of compatibility between the discrete continuity equation and the Poisson equation to get discrete energy estimates. More precisely, our definition (58) enables us to reproduce a discrete version of the following computation:

$$\int_{\mathbb{T}} \rho_{\varepsilon} u_{\varepsilon} \partial_x \phi_{\varepsilon} dx = - \int_{\mathbb{T}} \partial_x (\rho_{\varepsilon} u_{\varepsilon}) \phi_{\varepsilon} dx = \int_{\mathbb{T}} \partial_t \rho_{\varepsilon} \phi_{\varepsilon} dx.$$

Then, using the discrete Poisson equation, the implicitness in time yields some expected dissipation and it turns out that the additional spatial part also brings some dissipation since $g(0) \ge 0$. We also mention that the forcing term in the momentum equation is a priori not written as a gradient in space on the contrary to what is proposed in [7]. We nevertheless highlight the fact for the pressure-less Euler-Poisson equation, the theory of weak solutions is not fully understood. In our setting, we always consider strong solutions. Besides note that for $\varepsilon > 0$ fixed, for each $t \in [0, T]$, $\phi_{\varepsilon}(t)$ gains two spatial derivatives thanks to the standard elliptic regularity theory for the Poisson equation. Thus, the product $\rho_{\varepsilon}(t)\partial_x\phi_{\varepsilon}(t)$ is always a function even with $\rho_{\varepsilon}(t) \in L^1(\mathbb{T})$. Of course, what is more challenging is the case $\varepsilon \to 0$ since the Poisson equation becomes algebraic at the limit and there is no gain of regularity for $\phi_0$. The study of the convergence as $\varepsilon \to 0$ towards a weak entropic solution to (3) is to the best of our knowledge an open question. We may anyway, at the discrete level, always consider consistency of the scheme (48) for strong solutions. By the way, our modulated energy estimates in the limit $\varepsilon \to 0$ holds only for constant solutions. Last but not the least, note that the definition (58) does not yield straightforwardly $\tilde{\rho}_{i+\frac{1}{2}}^{n+1} \ge 0$ if $\rho^{n+1} > 0$ on $\mathbb{T}$. Nevertheless, we do observe that it can be re-written for $i \in \{0, \dots, N\}$ as

$$\tilde{\rho}_{i+\frac{1}{2}}^{n+1} = \rho_{i+1}^{n+1}(1 - \hat{g}(u_{i+\frac{1}{2}}^{n+1})) + \rho_{i}^{n+1}\hat{g}(u_{i+\frac{1}{2}}^{n+1}) \quad (63)$$

where $\hat{g}(u) = \begin{cases} \frac{g(u) - g(0)}{u} & \text{if } u \ne 0, \\ g'(0) & \text{if } u = 0. \end{cases}$ As a matter of fact, if $g$ is non decreasing and its Lipschitz constant is such that $\text{Lip}(g) \le 1$ then $\tilde{\rho}_{i+\frac{1}{2}}^{n+1}$ given by (63) is a convex combination of $\rho_{i+1}^{n+1}$ and $\rho_i^{n+1}$ so it preserves the $L^{\infty}$ boundedness provided $\rho^{n+1}$ is. Since this assumption is not needed in the discrete stability estimates, we shall not consider it as granted.

The first main result of this work is.

**Theorem 4.** (Unconditional Energy stability) Let $N_T$ and $N \in \mathbb{N}$ be two positive integers and $(\rho^0, u^0) \in X(\mathcal{T}) \times X(\mathcal{T}^*)$ being given in (49) such that $\rho^0 > 0$. Assume $\phi^0 \in X(\mathcal{T})$ verifies the discrete Poisson equation initially. Then, there exists a solution $(\rho^n, u^n, \phi^n)_{n=0,\dots,N_T} \subset X(\mathcal{T}) \times X(\mathcal{T}^*) \times X(\mathcal{T})$ to the scheme (48). In addition, this solution satisfies for $0 \le n \le N_T$ the estimates:

$$\rho^n > 0 \text{ on } \mathbb{T}, \quad (64)$$

$$\|\rho^n\|_{L^1(\mathbb{T})} = \|\rho^0\|_{L^1(\mathbb{T})}, \quad (65)$$

$$\mathcal{H}(\rho^n, u^n, \phi^n) + \Delta t \sum_{k=1}^{n-1} \tau^n = \mathcal{H}(\rho^0, u^0, \phi^0) \quad (66)$$

where $\tau^n \ge 0$ is given in (83) and $\phi^n$ satisfies the elliptic estimates (85)-(90) for $n \in \{0, \dots, N_T\}$.

Theorem 4 states an unconditional energy decay for the scheme (48). The second main result of this work is.

**Theorem 5.** (Modulated energy estimates around constant states with a small velocity) Let $N_T$ and $N$ be two positive integers and let $(\rho^n, u^n, \phi^n)_{n=0,\dots,N_T} \subset X(\mathcal{T}) \times X(\mathcal{T}^*) \times X(\mathcal{T})$ a solution to (48) such that $\rho^n > 0$ for all $n \in \{0, \dots, N_T\}$. Assume $\phi^0 \in X(\mathcal{T})$ verifies the discrete Poisson equation initially. Let $(\bar{u}, \bar{\phi}) \in X(\mathcal{T}^*) \times X(\mathcal{T})$ a constant state. Then the solution satisfies for $0 \le n \le N_T$ the following modulated energy estimates (recalling the definition (43) of $\mathcal{E}(\rho, u, \phi|\bar{u}, \bar{\phi})$):

a) If $\bar{u} = 0$ then

$$\forall n \in \{0, \dots, N_T\}, \quad \mathcal{E}(\rho^n, u^n, \phi^n|0, \bar{\phi}) \le \mathcal{E}(\rho^0, u^0, \phi^0|0, \bar{\phi}). \quad (67)$$
b) If $\bar{u} \neq 0$ and $|\bar{u}| \le \frac{g(0)}{\mathrm{Lip}(g)}$, provided there exists a constant $K \ge 0$ which is such that $K = \underset{\varepsilon \to 0}{\int} O(1)$, $K = \underset{\Delta x \to 0}{\int} O(1)$ and $\rho^n$ verifies:

$$ \forall n \in \{0, \dots, N_T\}, \quad \|\rho^n\|_{L^\infty(\mathbb{T})} \le K, \qquad (68) $$

then if $\Delta x = \underset{\varepsilon \to 0}{\int} O(\varepsilon)$ we have under the following CFL condition,

$$ |\bar{u}| \frac{\Delta t}{\Delta x} \left( 8\mathrm{Lip}(g) + 4 + \frac{2\Delta x^2}{\varepsilon^2} K \right) < 1 \qquad (69) $$

that

$$ \forall n \in \{0, \dots, N_T\}, \quad \mathcal{E}(\rho^n, u^n, \phi^n|\bar{u}, \bar{\phi}) \le a^n \mathcal{E}(\rho^0, u^0, \phi^0|\bar{u}, \bar{\phi}) \qquad (70) $$

with

$$ a = \frac{1}{1 - |\bar{u}| \frac{\Delta t}{\Delta x} \left( 8\mathrm{Lip}(g) + 4 + \frac{2\Delta x^2}{\varepsilon^2} K \right)}. \qquad (71) $$

Theorem 5 states an unconditional stability estimates around the constant state $(0, \bar{\phi})$. This result goes beyond the linearized studies [8, 7, 10] around constant states. As for the stability around an arbitrary constant state $(\bar{u}, \bar{\phi})$, we have not been able to address the problem in its full generality. It may look paradoxical regarding the case of a null velocity since the system (1) is Galilean invariant. Nevertheless, the discretization (48) is a priori not Galilean invariant, except is some special cases when the scheme is equivalent to a Lagrangian discretization of the continuity equation. The scheme (48) is fundamentally Eulerian. We are able to deal with the case where $\bar{u}$ is small enough. In such a case we require the density to be uniformly bounded with respect to $\varepsilon$ and to $\Delta x$. This assumption was already used in the literature [14, 2]. Provided the mesh size is of order $\varepsilon$ and a hyperbolic type CFL condition is verified, we are able to prove at most exponential growth of the modulated energy with a rate of increase which is bounded uniformly in $\varepsilon$. This is to compare with its continuous analogue (20) which is different in nature. The difficulty comes from the fact that the force term is not conservative at the discrete level:

$$ \forall n \in \{0, \dots, N_T - 1\}, \quad \sum_{i=0}^{N} \tilde{\rho}_{i+\frac{1}{2}}^{n+1} \delta(\phi^{n+1})_{i+\frac{1}{2}} \Delta x = 0 + \text{non zero residual terms.} \qquad (72) $$

The residual terms need to be controlled in the worst case by the modulated energy at step $n + 1$, this is precisely where the smallness condition on $\bar{u}$ and the CFL condition comes from. We mention that the assumption that $\Delta x$ is of order $\varepsilon$ seems us not so natural regarding what is expected for the spatial regularity of the perturbation. Indeed, in [20], it is proven that for strong solutions, we have $\frac{1}{\varepsilon}\|\rho_\varepsilon - \rho_0\|_{L^\infty H^s}$ is bounded uniformly in $\varepsilon$ for $s'$ large enough. So the spatial fluctuation and its high order spatial derivative are of order $\varepsilon$. We thus believe that the restriction on the mesh size is technical. We shall investigate this question in the numerical section. Theorem 5 is not fully satisfactory, it is nevertheless up to our knowledge, the first non linear discrete modulated energy estimate for the pressureless Euler-Poisson-Boltzmann equations. The rest of this section is devoted to the proofs of Theorem 4 and Theorem 5. In all this section $N_T$ and $N$ are fixed positive integers.

## 2.2 A priori estimates and existence

We begin with a discrete analogue of the renormalized continuity equation.

**Lemma 6.** (Discrete renormalized continuity equation) Let $(\rho^n, u^n, \phi^n)_{n=0, \dots, N_T} \subset X(\mathcal{T}) \times X(\mathcal{T}^*) \times X(\mathcal{T})$ be a solution to (48) with $\rho^n > 0$ for all $n \in \{0, \dots, N_T\}$. Let $\psi \in C^2(\mathbb{R}_+^*)$ a convex function. Then we have for $n \in \{0, \dots, N_T - 1\}$ and $i \in \{0, \dots, N\}$:

$$ \frac{\psi(\rho_i^{n+1}) - \psi(\rho_i^n)}{\Delta t} + \frac{\mathcal{G}_{i+\frac{1}{2}}^{n+1} - \mathcal{G}_{i-\frac{1}{2}}^{n+1}}{\Delta x} - \rho_i^{n+1} u_{i+\frac{1}{2}}^{n+1} \frac{\psi'(\rho_{i+1}^{n+1}) - \psi'(\rho_i^{n+1})}{\Delta x} = R_i^{n+1} + \mathcal{D}_{i+\frac{1}{2}}^{n+1} \qquad (73) $$
with

$$G_{i+\frac{1}{2}}^{n+1} = G^{n+1}(\rho_i^{n+1}, \rho_{i+1}^{n+1}, u_{i+\frac{1}{2}}^{n+1})\psi'(\rho_{i+1}^{n+1}) \quad (74)$$

$$D_{i+\frac{1}{2}}^{n+1} = \frac{1}{\Delta x} (G^{n+1}(\rho_i^{n+1}, \rho_{i+1}^{n+1}, u_{i+\frac{1}{2}}^{n+1}) - G^{n+1}(\rho_i^{n+1}, \rho_i^{n+1}, u_{i+\frac{1}{2}}^{n+1})) (\psi'(\rho_{i+1}^{n+1}) - \psi'(\rho_i^{n+1})) \le 0, \quad (75)$$

$$R_i^{n+1} = -\psi''(\xi_i^{n+1}) \frac{(\rho_i^n - \rho_i^{n+1})^2}{2\Delta t} \le 0, \quad (76)$$

$$\text{where } \xi_i^{n+1} \in (\min(\rho_i^{n+1}, \rho_i^n), \max(\rho_i^{n+1}, \rho_i^n)).$$

*Proof.* Let $n \in \{0, \dots, N_T - 1\}$ and $i \in \{0, \dots, N\}$. We multiply the continuity equation by $\psi'(\rho_i^{n+1})$ to get

$$\psi'(\rho_i^{n+1}) \frac{\rho_i^{n+1} - \rho_i^n}{\Delta t} + \frac{1}{\Delta x} (\mathcal{F}_{i+1/2}^{n+1} - \mathcal{F}_{i-1/2}^{n+1}) \psi'(\rho_i^{n+1}) = 0.$$

Let set

$$A_i^{n+1} = \psi'(\rho_i^{n+1}) \frac{\rho_i^{n+1} - \rho_i^n}{\Delta t} \quad \text{and} \quad B_i^{n+1} = \frac{1}{\Delta x} (\mathcal{F}_{i+1/2}^{n+1} - \mathcal{F}_{i-1/2}^{n+1}) \psi'(\rho_i^{n+1}).$$

Let us first focus on the temporal part $A_i^{n+1}$. Since $\psi \in C^2(\mathbb{R}_+)$ and $\rho_i^{n+1} > 0$ everywhere in $\mathbb{T}$, a Taylor-Lagrange expansion of $\psi(\rho_i^n)$ around $\rho_i^{n+1}$ yields the existence of $\xi_i^{n+1} \in (\min(\rho_i^{n+1}, \rho_i^n), \max(\rho_i^{n+1}, \rho_i^n))$ such that

$$\psi(\rho_i^n) = \psi(\rho_i^{n+1}) + \psi'(\rho_i^{n+1})(\rho_i^n - \rho_i^{n+1}) + \psi''(\xi_i^{n+1}) \frac{(\rho_i^n - \rho_i^{n+1})^2}{2},$$

so that

$$A_i^{n+1} = \psi'(\rho_i^{n+1}) \frac{\rho_i^{n+1} - \rho_i^n}{\Delta t} = \frac{\psi(\rho_i^{n+1}) - \psi(\rho_i^n)}{\Delta t} - R_i^{n+1},$$

with $R_i^{n+1}$ given by (76).

Let us now focus on the flux part $B_i^{n+1}$. We omit the time dependence at this step since all the quantities are defined at the discrete time $t_{n+1}$. We also set for ease $G_{i+\frac{1}{2}}(s,t) = G(s,t, u_{i+\frac{1}{2}}^{n+1})$. We have,

$$\begin{align*}
\Delta x B_i^{n+1} &= G_{i+1/2}(\rho_i, \rho_{i+1})\psi'(\rho_i) - G_{i-1/2}(\rho_{i-1}, \rho_i)\psi'(\rho_i) \\
&= G_{i+1/2}(\rho_i, \rho_{i+1})\psi'(\rho_{i+1}) - G_{i-1/2}(\rho_{i-1}, \rho_i)\psi'(\rho_i) + G_{i+1/2}(\rho_i, \rho_{i+1})(\psi'(\rho_i) - \psi'(\rho_{i+1})) \\
&= G_{i+1/2} - G_{i-1/2} - (G_{i+1/2}(\rho_i, \rho_{i+1}) - G_{i+1/2}(\rho_i, \rho_i))(\psi'(\rho_{i+1}) - \psi'(\rho_i)) \\
&\quad - G_{i+1/2}(\rho_i, \rho_i)(\psi'(\rho_{i+1}) - \psi'(\rho_i)) \\
&= G_{i+1/2} - G_{i-1/2} - \Delta x \mathcal{D}_{i+\frac{1}{2}} - \rho_i u_{i+1/2}(\psi'(\rho_{i+1}) - \psi'(\rho_i))
\end{align*}$$

where $\mathcal{G}_{i+1/2}$ and $\mathcal{D}_{i+\frac{1}{2}}$ are defined in (74) and (75). Note that $\mathcal{D}_{i+\frac{1}{2}}$ is non positive since $G$ is non increasing in its second variable and $\psi$ is convex. Hence, summing $A_i^{n+1}$ and $B_i^{n+1}$ together, we finally obtain

$$\frac{\psi(\rho_i^{n+1}) - \psi(\rho_i^n)}{\Delta t} + \frac{\mathcal{G}_{i+\frac{1}{2}}^{n+1} - \mathcal{G}_{i-\frac{1}{2}}^{n+1}}{\Delta x} - \rho_i^{n+1} u_{i+\frac{1}{2}} \frac{\psi'(\rho_{i+1}^{n+1}) - \psi'(\rho_i^{n+1})}{\Delta x} = R_i^{n+1} + \mathcal{D}_{i+\frac{1}{2}}^{n+1}.$$

□

We now give a discrete analogue of a renormalized momentum equation.

**Lemma 7.** Let $(\rho^n, u^n, \phi^n)_{n=0,\dots,N_T} \subset X(\mathcal{T}) \times X(\mathcal{T}^*) \times X(\mathcal{T})$ be a solution to (48) with $\rho^n > 0$ for all $n \in \{0, \dots, N_T\}$. Let $\psi \in C^2(\mathbb{R})$ a convex function. Then we have for $n \in \{0, \dots, N_T - 1\}$ and $i \in \{0, \dots, N\}$:

$$\frac{\rho_{i+\frac{1}{2}}^{n+1} \psi(u_{i+\frac{1}{2}}^{n+1}) - \rho_{i+\frac{1}{2}}^{n} \psi(u_{i+\frac{1}{2}}^{n})}{\Delta t} + \frac{\mathcal{F}_{i+1}^{n+1} \psi(u_{i+1}^{n+1}) - \mathcal{F}_{i}^{n+1} \psi(u_{i}^{n+1})}{\Delta x} = \tilde{\rho}_{i+\frac{1}{2}}^{n+1} \delta(\phi^{n+1})_{i+\frac{1}{2}} \psi'(u_{i+1}^{n+1}) - \mathcal{S}_{i}^{n+1}, \quad (77)$$
with

$$S_i^{n+1} = -\mathcal{F}_{i+1}^{n+1} \frac{(u_{i+1}^{n+1} - u_{i+\frac{1}{2}}^{n+1})^2}{2\Delta x} \psi''(\beta_{i+\frac{1}{2}}^{n+1}) + \mathcal{F}_i^{n+1} \frac{(u_i^{n+1} - u_{i+\frac{1}{2}}^{n+1})^2}{2\Delta x} \psi''(\gamma_{i+\frac{1}{2}}^{n+1}) + \rho_{i+\frac{1}{2}}^n \frac{(u_{i+\frac{1}{2}}^n - u_{i+\frac{1}{2}}^{n+1})^2}{2\Delta t} \psi''(\alpha_{i+\frac{1}{2}}^{n+1}), \quad (78)$$

where $S_i^{n+1}$ is non negative because of the definition of the velocities at the interfaces (57) and the fact that $\psi$ is a convex function.

*Proof.* Let $i \in \{0, \dots, N\}$. We multiply the discrete momentum equation by $\psi'(u_{i+\frac{1}{2}}^{n+1})$ to get

$$A_{i+\frac{1}{2}}^{n+1} + B_{i+\frac{1}{2}}^{n+1} = P_{i+\frac{1}{2}}^{n+1}.$$

with $P_{i+\frac{1}{2}}^{n+1} = \tilde{\rho}_{i+\frac{1}{2}}^{n+1} \delta(\phi^{n+1})_{i+\frac{1}{2}} \psi'(u_{i+\frac{1}{2}}^{n+1})$ and

$$A_{i+\frac{1}{2}}^{n+1} = \frac{\rho_{i+\frac{1}{2}}^{n+1} u_{i+\frac{1}{2}}^{n+1} - \rho_{i+\frac{1}{2}}^{n} u_{i+\frac{1}{2}}^{n}}{\Delta t} \psi'(u_{i+\frac{1}{2}}^{n+1}), \quad B_{i+\frac{1}{2}}^{n+1} = \frac{\mathcal{F}_{i+1}^{n+1} u_{i+1}^{n+1} - \mathcal{F}_i^{n+1} u_i^{n+1}}{\Delta x} \psi'(u_{i+\frac{1}{2}}^{n+1}).$$

We shall now reformulate $A_{i+\frac{1}{2}}^{n+1}$ and $B_{i+\frac{1}{2}}^{n+1}$. We have

$$\begin{aligned} \Delta t A_{i+\frac{1}{2}}^{n+1} &= \rho_{i+\frac{1}{2}}^{n+1} u_{i+\frac{1}{2}}^{n+1} \psi'(u_{i+\frac{1}{2}}^{n+1}) - \rho_{i+\frac{1}{2}}^{n} (u_{i+\frac{1}{2}}^{n} - u_{i+\frac{1}{2}}^{n+1}) \psi'(u_{i+\frac{1}{2}}^{n+1}) - \rho_{i+\frac{1}{2}}^{n} u_{i+\frac{1}{2}}^{n+1} \psi'(u_{i+\frac{1}{2}}^{n+1}) \\ &= u_{i+\frac{1}{2}}^{n+1} \psi'(u_{i+\frac{1}{2}}^{n+1}) (\rho_{i+\frac{1}{2}}^{n+1} - \rho_{i+\frac{1}{2}}^{n}) - \rho_{i+\frac{1}{2}}^{n} (u_{i+\frac{1}{2}}^{n} - u_{i+\frac{1}{2}}^{n+1}) \psi'(u_{i+\frac{1}{2}}^{n+1}). \end{aligned}$$

A Taylor-Lagrange expansion of $\psi(u_{i+\frac{1}{2}}^n)$ around $u_{i+\frac{1}{2}}^{n+1}$ yields the existence of $\alpha_{i+\frac{1}{2}}^{n+1} \in (\min(u_{i+\frac{1}{2}}^n, u_{i+\frac{1}{2}}^{n+1}), \max(u_{i+\frac{1}{2}}^n, u_{i+\frac{1}{2}}^{n+1}))$ such that

$$\psi(u_{i+\frac{1}{2}}^n) = \psi(u_{i+\frac{1}{2}}^{n+1}) + \psi'(u_{i+\frac{1}{2}}^{n+1})(u_{i+\frac{1}{2}}^n - u_{i+\frac{1}{2}}^{n+1}) + \psi''(\alpha_{i+\frac{1}{2}}^{n+1}) \frac{(u_{i+\frac{1}{2}}^n - u_{i+\frac{1}{2}}^{n+1})^2}{2}.$$

Inserting this expansion in $A_{i+\frac{1}{2}}^{n+1}$ leads to

$$\begin{aligned} \Delta t A_{i+\frac{1}{2}}^{n+1} &= u_{i+\frac{1}{2}}^{n+1} \psi'(u_{i+\frac{1}{2}}^{n+1}) (\rho_{i+\frac{1}{2}}^{n+1} - \rho_{i+\frac{1}{2}}^{n}) - \rho_{i+\frac{1}{2}}^{n} \left[ \psi(u_{i+\frac{1}{2}}^{n}) - \psi(u_{i+\frac{1}{2}}^{n+1}) - \frac{(u_{i+\frac{1}{2}}^{n} - u_{i+\frac{1}{2}}^{n+1})^2}{2} \psi''(\alpha_{i+\frac{1}{2}}^{n+1}) \right] \\ &= \rho_{i+\frac{1}{2}}^{n} \psi(u_{i+\frac{1}{2}}^{n+1}) - \rho_{i+\frac{1}{2}}^{n} \psi(u_{i+\frac{1}{2}}^{n}) + u_{i+\frac{1}{2}}^{n+1} \psi'(u_{i+\frac{1}{2}}^{n+1}) (\rho_{i+\frac{1}{2}}^{n+1} - \rho_{i+\frac{1}{2}}^{n}) + \rho_{i+\frac{1}{2}}^{n} \frac{(u_{i+\frac{1}{2}}^{n} - u_{i+\frac{1}{2}}^{n+1})^2}{2} \psi''(\alpha_{i+\frac{1}{2}}^{n+1}). \end{aligned}$$

The continuity equation on the dual mesh writes:

$$\rho_{i+\frac{1}{2}}^{n+1} + \frac{\Delta t}{\Delta x} (\mathcal{F}_{i+1}^{n+1} - \mathcal{F}_i^{n+1}) = \rho_{i+\frac{1}{2}}^{n},$$

hence,

$$A_{i+\frac{1}{2}}^{n+1} = \frac{\rho_{i+\frac{1}{2}}^{n+1} \psi(u_{i+\frac{1}{2}}^{n+1}) - \rho_{i+\frac{1}{2}}^{n} \psi(u_{i+\frac{1}{2}}^{n})}{\Delta t} + \frac{\mathcal{F}_{i+1}^{n+1} - \mathcal{F}_i^{n+1}}{\Delta x} \psi(u_{i+\frac{1}{2}}^{n+1}) + u_{i+\frac{1}{2}}^{n+1} \psi'(u_{i+\frac{1}{2}}^{n+1}) \frac{\rho_{i+\frac{1}{2}}^{n+1} - \rho_{i+\frac{1}{2}}^{n}}{\Delta t} + \rho_{i+\frac{1}{2}}^{n} \frac{(u_{i+\frac{1}{2}}^{n} - u_{i+\frac{1}{2}}^{n+1})^2}{2\Delta t} \psi''(\alpha_{i+\frac{1}{2}}^{n+1}). \quad (79)$$

We now deal with $B_{i+\frac{1}{2}}^{n+1}$ which we rewrite as

$$\begin{aligned} \Delta x B_{i+\frac{1}{2}}^{n+1} = \mathcal{F}_{i+1}^{n+1} (u_{i+1}^{n+1} - u_{i+\frac{1}{2}}^{n+1}) \psi'(u_{i+\frac{1}{2}}^{n+1}) - \mathcal{F}_i^{n+1} (u_i^{n+1} - u_{i+\frac{1}{2}}^{n+1}) \psi'(u_{i+\frac{1}{2}}^{n+1}) \\ + u_{i+\frac{1}{2}}^{n+1} \psi'(u_{i+\frac{1}{2}}^{n+1}) (\mathcal{F}_{i+1}^{n+1} - \mathcal{F}_i^{n+1}). \end{aligned}$$
The Taylor-Lagrange expansions of $\psi$ yield the existence of $\beta_{i+\frac{1}{2}}^{n+1} \in (\min(u_{i+1}^{n+1}, u_{i+\frac{1}{2}}^{n+1}), \max(u_{i+1}^{n+1}, u_{i+\frac{1}{2}}^{n+1}))$ and $\gamma_{i+\frac{1}{2}}^{n+1} \in (\min(u_i^{n+1}, u_{i+\frac{1}{2}}^{n+1}), \max(u_i^{n+1}, u_{i+\frac{1}{2}}^{n+1}))$ such that

$$ \psi(u_{i+1}^{n+1}) = \psi(u_{i+\frac{1}{2}}^{n+1}) + (u_{i+1}^{n+1} - u_{i+\frac{1}{2}}^{n+1})\psi'(u_{i+\frac{1}{2}}^{n+1}) + \frac{(u_{i+1}^{n+1} - u_{i+\frac{1}{2}}^{n+1})^2}{2}\psi''(\beta_{i+\frac{1}{2}}^{n+1}), $$

$$ \psi(u_i^{n+1}) = \psi(u_{i+\frac{1}{2}}^{n+1}) + (u_i^{n+1} - u_{i+\frac{1}{2}}^{n+1})\psi'(u_{i+\frac{1}{2}}^{n+1}) + \frac{(u_i^{n+1} - u_{i+\frac{1}{2}}^{n+1})^2}{2}\psi''(\gamma_{i+\frac{1}{2}}^{n+1}). $$

Inserting these expressions in $B_{i+\frac{1}{2}}^{n+1}$ leads to

$$ \begin{aligned} \Delta x B_{i+\frac{1}{2}}^{n+1} = {}& \mathcal{F}_{i+1}^{n+1} [\psi(u_{i+1}^{n+1}) - \psi(u_{i+\frac{1}{2}}^{n+1}) - \frac{(u_{i+1}^{n+1} - u_{i+\frac{1}{2}}^{n+1})^2}{2} \psi''(\beta_{i+\frac{1}{2}}^{n+1})] \\ & - \mathcal{F}_i^{n+1} [\psi(u_i^{n+1}) - \psi(u_{i+\frac{1}{2}}^{n+1}) - \frac{(u_i^{n+1} - u_{i+\frac{1}{2}}^{n+1})^2}{2} \psi''(\gamma_{i+\frac{1}{2}}^{n+1})] + u_{i+\frac{1}{2}}^{n+1} \psi'(u_{i+\frac{1}{2}}^{n+1}) (\mathcal{F}_{i+1}^{n+1} - \mathcal{F}_i^{n+1}). \end{aligned} $$

Rearranging the terms we eventually obtain

$$ \begin{aligned} \Delta x B_{i+\frac{1}{2}}^{n+1} = {}& \mathcal{F}_{i+1}^{n+1} \psi(u_{i+1}^{n+1}) - \mathcal{F}_i^{n+1} \psi(u_i^{n+1}) - \psi(u_{i+\frac{1}{2}}^{n+1}) (\mathcal{F}_{i+1}^{n+1} - \mathcal{F}_i^{n+1}) + u_{i+\frac{1}{2}}^{n+1} \psi'(u_{i+\frac{1}{2}}^{n+1}) (\mathcal{F}_{i+1}^{n+1} - \mathcal{F}_i^{n+1}) \\ & - \mathcal{F}_{i+1}^{n+1} \frac{(u_{i+1}^{n+1} - u_{i+\frac{1}{2}}^{n+1})^2}{2} \psi''(\beta_{i+\frac{1}{2}}^{n+1}) + \mathcal{F}_i^{n+1} \frac{(u_i^{n+1} - u_{i+\frac{1}{2}}^{n+1})^2}{2} \psi''(\gamma_{i+\frac{1}{2}}^{n+1}). \end{aligned} \quad (80) $$

Using eventually the continuity equation on the dual mesh

$$ \frac{\rho_{i+\frac{1}{2}}^{n+1} - \rho_{i+\frac{1}{2}}^{n}}{\Delta t} + \frac{\mathcal{F}_{i+1}^{n+1} - \mathcal{F}_{i}^{n+1}}{\Delta x} = 0, $$

we get

$$ \begin{aligned} \Delta x B_{i+\frac{1}{2}}^{n+1} = {}& \mathcal{F}_{i+1}^{n+1} \psi(u_{i+1}^{n+1}) - \mathcal{F}_i^{n+1} \psi(u_i^{n+1}) - \psi(u_{i+\frac{1}{2}}^{n+1}) (\mathcal{F}_{i+1}^{n+1} - \mathcal{F}_i^{n+1}) - u_{i+\frac{1}{2}}^{n+1} \psi'(u_{i+\frac{1}{2}}^{n+1}) (\rho_{i+\frac{1}{2}}^{n+1} - \rho_{i+\frac{1}{2}}^{n}) \frac{\Delta x}{\Delta t} \\ & - \mathcal{F}_{i+1}^{n+1} \frac{(u_{i+1}^{n+1} - u_{i+\frac{1}{2}}^{n+1})^2}{2} \psi''(\beta_{i+\frac{1}{2}}^{n+1}) + \mathcal{F}_i^{n+1} \frac{(u_i^{n+1} - u_{i+\frac{1}{2}}^{n+1})^2}{2} \psi''(\gamma_{i+\frac{1}{2}}^{n+1}). \end{aligned} $$

Finally, we glean $A_{i+\frac{1}{2}}^{n+1}$ and $B_{i+\frac{1}{2}}^{n+1}$

$$ \begin{aligned} A_{i+\frac{1}{2}}^{n+1} + B_{i+\frac{1}{2}}^{n+1} = {}& \frac{\rho_{i+\frac{1}{2}}^{n+1} \psi(u_{i+\frac{1}{2}}^{n+1}) - \rho_{i+\frac{1}{2}}^{n} \psi(u_{i+\frac{1}{2}}^{n})}{\Delta t} + \frac{\mathcal{F}_{i+1}^{n+1} \psi(u_{i+1}^{n+1}) - \mathcal{F}_i^{n+1} \psi(u_i^{n+1})}{\Delta x} \\ & - \mathcal{F}_{i+1}^{n+1} \frac{(u_{i+1}^{n+1} - u_{i+\frac{1}{2}}^{n+1})^2}{2\Delta x} \psi''(\beta_{i+\frac{1}{2}}^{n+1}) + \mathcal{F}_i^{n+1} \frac{(u_i^{n+1} - u_{i+\frac{1}{2}}^{n+1})^2}{2\Delta x} \psi''(\gamma_{i+\frac{1}{2}}^{n+1}) + \rho_{i+\frac{1}{2}}^{n} \frac{(u_{i+\frac{1}{2}}^{n} - u_{i+\frac{1}{2}}^{n+1})^2}{2\Delta t} \psi''(\alpha_{i+\frac{1}{2}}^{n+1}). \end{aligned} \quad (81) $$

We obtained the expected result. $\square$

A consequence of Lemmas 6, 7 and the definition (58) is that the discrete energy functional (43) decays along the solutions of the scheme (48). More precisely we have.

**Proposition 3.** (Discrete energy decay) Let $(\rho^n, u^n, \phi^n)_{n=0,...,N_T} \subset X(\mathcal{T}) \times X(\mathcal{T}^*) \times X(\mathcal{T})$ a solution to (48) with $\rho^n > 0$ for all $n \in \{0,...,N_T\}$. Assume $\phi^0 \in X(\mathcal{T})$ satisfies the discrete Poisson equation initially. Then we have

$$ \forall n \in \{0, \dots, N_T - 1\} \quad \mathcal{H}(\rho^{n+1}, u^{n+1}, \phi^{n+1}) - \mathcal{H}(\rho^n, u^n, \phi^n) = -\Delta t \tau^{n+1}, \quad (82) $$
where $\tau^{n+1} \ge 0$ is given by

$$ \tau^{n+1} = \sum_{i=0}^{N} (S_i^{n+1} + e^{-\xi_i^n} \frac{(\phi_i^{n+1} - \phi_i^n)^2}{2\Delta t}) \Delta x + \frac{\varepsilon^2}{2\Delta t} \sum_{i=0}^{N} \left( \frac{\phi_{i+1}^{n+1} - \phi_i^{n+1}}{\Delta x} - \frac{\phi_{i+1}^n - \phi_i^n}{\Delta x} \right)^2 \Delta x \quad (83) \\ + \frac{g(0)}{\Delta x} \sum_{i=0}^{N} \varepsilon^2 |\delta(\phi^{n+1})_{i+\frac{1}{2}} - \delta(\phi^{n+1})_{i-\frac{1}{2}}|^2 \Delta x + \frac{g(0)}{\Delta x} \sum_{i=0}^{N} |(e^{-\phi_i^{n+1}} - e^{-\phi_{i-1}^{n+1}})(\phi_i^{n+1} - \phi_{i-1}^{n+1})| \Delta x $$

where $S_i^{n+1} \ge 0$ is given (78) with $\psi: s \in \mathbb{R} \mapsto \frac{s^2}{2}$ and for all $i \in \{0, \dots, N\}$, $\xi_i^n \in (\min(\phi_i^n, \phi_i^{n+1}), \max(\phi_i^n, \phi_i^{n+1}))$.

*Proof.* Let $n \in \{0, \dots, N_T - 1\}$. Observe on the one hand that applying Lemma 7 with the function $\psi: s \in \mathbb{R} \mapsto \frac{s^2}{2}$ yields for all $i \in \{0, \dots, N\}$,

$$ \frac{\rho_{i+\frac{1}{2}}^{n+1} (u_{i+\frac{1}{2}}^{n+1})^2 - \rho_{i+\frac{1}{2}}^n (u_{i+\frac{1}{2}}^n)^2}{2\Delta t} + \frac{\mathcal{F}_{i+1}^{n+1} (u_{i+1}^{n+1})^2 - \mathcal{F}_i^{n+1} (u_i^{n+1})^2}{2\Delta x} = \tilde{\rho}_{i+\frac{1}{2}}^{n+1} u_{i+\frac{1}{2}}^{n+1} \delta(\phi^{n+1})_{i+\frac{1}{2}} - S_i^{n+1} \quad (84) $$

where $S_i^{n+1} \ge 0$ is given by (78). Let us treat the first term of the right hand side. Using the definition (58) we have

$$ D := \sum_{i=0}^{N} \rho_{i+\frac{1}{2}}^{n+1} u_{i+\frac{1}{2}}^{n+1} \delta(\phi^{n+1})_{i+\frac{1}{2}} \Delta x = \sum_{i=0}^{N} \mathcal{F}_{i+\frac{1}{2}}^{n+1} \delta(\phi^{n+1})_{i+\frac{1}{2}} \Delta x - \sum_{i=0}^{N} G(\rho_i^{n+1}, \rho_{i+1}^{n+1}, 0) \delta(\phi^{n+1})_{i+\frac{1}{2}} \Delta x. $$

Let us set

$$ D_1 = \sum_{i=0}^{N} \mathcal{F}_{i+\frac{1}{2}}^{n+1} \delta(\phi^{n+1})_{i+\frac{1}{2}} \Delta x, $$

and

$$ D_2 = - \sum_{i=0}^{N} G(\rho_i^{n+1}, \rho_{i+1}^{n+1}, 0) \delta(\phi^{n+1})_{i+\frac{1}{2}} \Delta x = g(0) \sum_{i=0}^{N} (\rho_{i+1}^{n+1} - \rho_i^{n+1}) \delta(\phi^{n+1})_{i+\frac{1}{2}} \Delta x. $$

As for the first term, using the discrete integration by parts (41) of Lemma 2 and the continuity equation, we have

$$ D_1 = - \sum_{i=0}^{N} \frac{\mathcal{F}_{i+\frac{1}{2}}^{n+1} - \mathcal{F}_{i-\frac{1}{2}}^{n+1}}{\Delta x} \phi_i^{n+1} \Delta x = \sum_{i=0}^{N} \frac{\rho_i^{n+1} - \rho_i^n}{\Delta t} \phi_i^{n+1} \Delta x. $$

Then, using the discrete Poisson equation (which also holds initially), we obtain

$$ \begin{aligned} & \sum_{i=0}^{N} \frac{\rho_i^{n+1} - \rho_i^n}{\Delta t} \phi_i^{n+1} \Delta x = \frac{\varepsilon^2}{\Delta t} \sum_{i=0}^{N} \left( \Delta(\phi^{n+1})_i - \Delta(\phi^n)_i \right) \phi_i^{n+1} \Delta x + \frac{1}{\Delta t} \sum_{i=0}^{N} \left( e^{-\phi_i^{n+1}} - e^{-\phi_i^n} \right) \phi_i^{n+1} \Delta x \\ & = - \frac{\varepsilon^2}{\Delta t} \sum_{i=0}^{N} \left( \delta(\phi^{n+1})_{i+\frac{1}{2}} - \delta(\phi^n)_{i+\frac{1}{2}} \right) \delta(\phi^{n+1})_{i+\frac{1}{2}} \Delta x + \frac{1}{\Delta t} \sum_{i=0}^{N} \left( e^{-\phi_i^{n+1}} - e^{-\phi_i^n} \right) \phi_i^{n+1} \Delta x, \end{aligned} $$

where we used the discrete integration by parts (42). We now use the identity $-a(a-b) = -a^2/2 + b^2/2 - (a-b)^2/2$ for $a, b \in \mathbb{R}$, so the first term re-writes

$$ \begin{aligned} & - \frac{\varepsilon^2}{\Delta t} \sum_{i=0}^{N} \left( \delta(\phi^{n+1})_{i+\frac{1}{2}} - \delta(\phi^n)_{i+\frac{1}{2}} \right) (\delta\phi^{n+1})_{i+\frac{1}{2}} \Delta x \\ & = - \frac{\varepsilon^2}{2\Delta t} \sum_{i=0}^{N} |\delta(\phi^{n+1})_{i+\frac{1}{2}}|^2 \Delta x + \frac{\varepsilon^2}{2\Delta t} \sum_{i=0}^{N} |\delta(\phi^n)_{i+\frac{1}{2}}|^2 \Delta x - \frac{\varepsilon^2}{2\Delta t} \sum_{i=0}^{N} \left( \delta(\phi^{n+1})_{i+\frac{1}{2}} - \delta(\phi^n)_{i+\frac{1}{2}} \right)^2 \Delta x. \end{aligned} $$
As for the second term we have

$$
\begin{align*}
& \frac{1}{\Delta t} \sum_{i=0}^{N} (e^{-\phi_i^{n+1}} - e^{-\phi_i^n})(\phi_i^{n+1} + 1 - \phi_i^n - 1 + \phi_i^n)\Delta x \\
&= -\frac{1}{\Delta t} \sum_{i=0}^{N} h(\phi_i^{n+1})\Delta x + \frac{1}{\Delta t} \sum_{i=0}^{N} h(\phi_i^n)\Delta x - \frac{1}{\Delta t} \sum_{i=0}^{N} (e^{-\phi_i^{n+1}} - e^{-\phi_i^n} + e^{-\phi_i^n}(\phi_i^{n+1} - \phi_i^n))\Delta x \\
&= -\frac{1}{\Delta t} \sum_{i=0}^{N} h(\phi_i^{n+1})\Delta x + \frac{1}{\Delta t} \sum_{i=0}^{N} h(\phi_i^n)\Delta x - \frac{1}{\Delta t} \sum_{i=0}^{N} e^{-\xi_i^n} \frac{(\phi_i^{n+1} - \phi_i^n)^2}{2}\Delta x,
\end{align*}
$$

where Taylor-Lagrange expansion has been performed with $\xi_i^n \in (\min(\phi_i^n, \phi_i^{n+1}), \max(\phi_i^n, \phi_i^{n+1}))$. Gathering the terms together we obtain

$$
\begin{align*}
\Delta t D_1 = & -\sum_{i=0}^{N} \frac{\varepsilon^2}{2} |\delta(\phi^{n+1})_{i+\frac{1}{2}}|^2 \Delta x + \sum_{i=0}^{N} \frac{\varepsilon^2}{2} |\delta(\phi^n)_{i+\frac{1}{2}}|^2 \Delta x \\
& -\sum_{i=0}^{N} h(\phi_i^{n+1}) \Delta x + \sum_{i=0}^{N} h(\phi_i^n) \Delta x \\
& -\frac{\varepsilon^2}{2} \sum_{i=0}^{N} \left( \delta(\phi^{n+1})_{i+\frac{1}{2}} - \delta(\phi^n)_{i+\frac{1}{2}} \right)^2 \Delta x - \sum_{i=0}^{N} e^{-\xi_i^n} \frac{(\phi_i^{n+1} - \phi_i^n)^2}{2} \Delta x.
\end{align*}
$$

We now treat the second term $D_2$:

$$
D_2 = g(0) \sum_{i=0}^{N} (\rho_{i+1}^{n+1} - \rho_i^{n+1}) \delta(\phi^{n+1})_{i+\frac{1}{2}} \Delta x = - \frac{g(0)}{\Delta x} \sum_{i=0}^{N} \rho_i^{n+1} (\phi_{i+1}^{n+1} - 2\phi_i^{n+1} + \phi_{i-1}^{n+1}) \Delta x,
$$

where we used a discrete integration by parts to get the last equality. Then using the discrete Poisson equation we get

$$
D_2 = - \frac{g(0)}{\Delta x} \sum_{i=0}^{N} \varepsilon^2 |\delta(\phi^{n+1})_{i+\frac{1}{2}} - \delta(\phi^{n+1})_{i-\frac{1}{2}}|^2 \Delta x - \frac{g(0)}{\Delta x} \sum_{i=0}^{N} e^{-\phi_i^{n+1}} (\phi_{i+1}^{n+1} - \phi_i^{n+1} - (\phi_i^{n+1} - \phi_{i-1}^{n+1})) \Delta x.
$$

Since $g(0) \ge 0$ the first sum is non positive. The second sum is also non positive since after using a change of indices we have

$$
- \sum_{i=0}^{N} e^{-\phi_i^{n+1}} (\phi_{i+1}^{n+1} - \phi_i^{n+1} - (\phi_i^{n+1} - \phi_{i-1}^{n+1})) \Delta x = \sum_{i=0}^{N} (e^{-\phi_i^{n+1}} - e^{-\phi_{i-1}^{n+1}}) (\phi_i^{n+1} - \phi_{i-1}^{n+1}) \Delta x.
$$

Since the function $s \mapsto e^{-s}$ is decreasing we thus get $\sum_{i=0}^{N} (e^{-\phi_i^{n+1}} - e^{-\phi_{i-1}^{n+1}})(\phi_i^{n+1} - \phi_{i-1}^{n+1})\Delta x \le 0$.
Thus, it yields for $D_2$

$$
D_2 = - \frac{g(0)}{\Delta x} \sum_{i=0}^{N} \varepsilon^2 |\delta(\phi^{n+1})_{i+\frac{1}{2}} - \delta(\phi^{n+1})_{i-\frac{1}{2}}|^2 \Delta x - \frac{g(0)}{\Delta x} \sum_{i=0}^{N} |(e^{-\phi_i^{n+1}} - e^{-\phi_{i-1}^{n+1}})(\phi_i^{n+1} - \phi_{i-1}^{n+1})| \Delta x.
$$

So we eventually obtain,

$$
\Delta t D := \Delta t \sum_{i=0}^{N} \rho_{i+\frac{1}{2}}^{n+1} u_{i+\frac{1}{2}}^{n+1} \delta(\phi^{n+1})_{i+\frac{1}{2}} \Delta x = \Delta t D_1 + \Delta t D_2,
$$

which enables to recover the terms in (83). Summing (84) over $i \in \{0, \dots, N\}$ then leads to the discrete kinetic energy balance. $\square$

We now establish discrete elliptic estimates for the discrete Poisson equation.
**Proposition 4.** (Discrete elliptic estimates) Let $(\rho^n, u^n, \phi^n)_{n=0,...,N_T} \subset X(\mathcal{T}) \times X(\mathcal{T}^*) \times X(\mathcal{T})$ a solution to (48) with $\rho^n > 0$ for all $n \in \{0,...,N_T\}$. Assume $\phi^0 \in X(\mathcal{T})$ satisfies the discrete Poisson equation initially. Then the discrete potential satisfies for all $n \in \{0,...,N_T\}$ the following estimates:

$$ \forall p \in [1, +\infty), \quad \left\| e^{-\phi^n} \right\|_{L^p(\mathcal{T})} \le \|\rho^n\|_{L^p(\mathcal{T})}, \qquad (85) $$

$$ \forall p \in [1, +\infty), \quad \varepsilon^2 \sum_{i=0}^{N} \left| \frac{\phi_{i+1}^n - \phi_i^n}{\Delta x} \frac{e^{-(p-1)\phi_{i+1}^n} - e^{-(p-1)\phi_i^n}}{\Delta x} \right| \Delta x \le \|\rho^n\|_{L^p(\mathcal{T})}^p, \qquad (86) $$

$$ \forall i \in \{0, \dots, N\}, \quad \min_{i \in \{0, \dots, N\}} \log(\rho_i^n) \le -\phi_i^n \le \max_{i \in \{0, \dots, N\}} \log(\rho_i^n), \qquad (87) $$

$$ \min_{i \in \{0, \dots, N\}} \log(\rho_i^n) \le \langle -\phi^n \rangle \le \max_{i \in \{0, \dots, N\}} \log(\rho_i^n), \qquad (88) $$

$$ \alpha^n \left\| \phi^n - \langle \phi^n \rangle \right\|_{L^2(\mathcal{T})} \le \|\rho^n\|_{L^2(\mathcal{T})} + \|\rho^n\|_{L^1(\mathcal{T})}, \qquad (89) $$

$$ \varepsilon \left| \phi^n - \langle \phi^n \rangle \right|_{H^1(\mathcal{T})} \le \frac{1}{\sqrt{\alpha^n}} \left( \| \rho^n \|_{L^2(\mathcal{T})} + \| \rho^n \|_{L^1(\mathcal{T})} \right), \qquad (90) $$

with $\langle \phi^n \rangle = \sum_{i=0}^{N} \phi_i^n \Delta x$, $\alpha^n = e^{\min_{i \in \{0, \dots, N\}} \log(\rho_i^n)}$ and where the inequality (85) is an equality for $p=1$.

*Proof.* We prove (85). Fix $p \in [1, +\infty)$ and $n \in \{0, \dots, N_T\}$. We set $\psi = -\phi^n$. We multiply the discrete Poisson equation by $e^{(p-1)\psi_i}$ for $i \in \{0, \dots, N\}$ and sum over $i \in \{0, \dots, N\}$. We get after a discrete integration by parts

$$ \sum_{i=0}^{N} \varepsilon^2 \frac{\psi_{i+1} - \psi_i}{\Delta x} \frac{e^{(p-1)\psi_{i+1}} - e^{(p-1)\psi_i}}{\Delta x} \Delta x + \sum_{i=0}^{N} e^{p\psi_i} \Delta x = \sum_{i=0}^{N} \rho_i^n e^{(p-1)\psi_i} \Delta x. \qquad (91) $$

Since $p \ge 1$ and the exponential function is increasing, the first sum is non negative. As a result we obtain on the one hand

$$ \sum_{i=0}^{N} |e^{\psi_i}|^p \Delta x \le \sum_{i=0}^{N} \rho_i^n e^{(p-1)\psi_i} \Delta x \le \left( \sum_{i=0}^{N} |\rho_i^n|^p \Delta x \right)^{\frac{1}{p}} \left( \sum_{i=0}^{N} |e^{\psi_i}|^p \Delta x \right)^{\frac{p-1}{p}}, $$

where for the last inequality we used the Hölder inequality in the duality $(\ell^p(\mathbb{R}^{N+1}), \ell^{\frac{p}{p-1}}(\mathbb{R}^{N+1}))$. It yields (85) after a simplification. The case $p=1$ yields the equality in (85) directly from (91). Using again (91) we obtain on the other hand,

$$ \sum_{i=0}^{N} \varepsilon^2 \frac{\psi_{i+1} - \psi_i}{\Delta x} \frac{e^{(p-1)\psi_{i+1}} - e^{(p-1)\psi_i}}{\Delta x} \Delta x \le \|\rho^n\|_{L^p(\mathcal{T})} \|e^{-\phi^n}\|_{L^p(\mathcal{T})}^{p-1} \le \|\rho^n\|_{L^p(\mathcal{T})}^p. $$

It proves (86). We prove the maximum principle (87). We set $M = \max_{i \in \{0, \dots, N\}} \log(\rho_i^n)$ and $\psi = -\phi^n$. Then for all $i \in \{0, \dots, N\}$:

$$ -\varepsilon^2 \frac{\psi_{i+1} - 2\psi_i + \psi_{i-1}}{\Delta x^2} + e^{\psi_i} - e^M = \rho_i^n - e^M. \qquad (92) $$

By definition of $M$ we have $\rho_i^n - e^M \le 0$ for all $i \in \{0, \dots, N\}$. Let us now prove the claim by contradiction. Assume that there exists $j \in \{0, \dots, N\}$ such that $\psi_j > M$. Then let $i^* \in \{0, \dots, N\}$ such that $\psi_{i^*} = \max_{i \in \{0, \dots, N\}} \psi_i$. Then we have $i \in \{0, \dots, N\}$

$$ \psi_{i^*} > M \implies e^{\psi_{i^*}} > e^M. $$

In addition,

$$ -\varepsilon^2 \frac{\psi_{i^*+1} - 2\psi_{i^*} + \psi_{i^*-1}}{\Delta x^2} \ge 0 $$

so we deduce from (92) applied at $i = i^*$ that $\rho_{i^*} > e^M$ which is the expected contradiction. A similar reasoning with $m = \min_{i \in \{0, \dots, N\}} \log(\rho_i^n)$ yields the lower bound in (87). The estimate (88) is an immediate
consequence of (87) since $\sum_{i=0}^{N} \Delta x = 1$. We now prove the $L^2(\mathbb{T})$ estimate (89). We set here $\psi = -\phi^n + \langle \phi^n \rangle$.

Then we have for $i \in \{0, \dots, N\}$:

$$ -\varepsilon^2 \frac{\psi_{i+1} - 2\psi_i + \psi_{i-1}}{\Delta x^2} + e^{-\phi_i^n} - e^{\langle \phi^n \rangle} = \rho_i^n - e^{\langle \phi^n \rangle}. \quad (93) $$

We multiply (93) by $\psi_i$ for all $i \in \{0, \dots, N\}$ and sum. Then after a discrete integration by parts we have

$$ \sum_{i=0}^{N} \varepsilon^2 \left| \frac{\psi_{i+1} - \psi_i}{\Delta x} \right|^2 \Delta x + \sum_{i=0}^{N} (e^{-\phi_i^n} - e^{\langle \phi^n \rangle}) (-\phi_i^n + \langle \phi^n \rangle) \Delta x = \sum_{i=0}^{N} (\rho_i^n - e^{\langle \phi^n \rangle}) (-\phi_i^n + \langle \phi^n \rangle) \Delta x. \quad (94) $$

Using the mean value theorem, we have for each $i \in \{0, \dots, N\}$ the existence of $\xi_i^n \in (\min(-\phi_i^n, \langle \phi^n \rangle), \max(-\phi_i^n, \langle \phi^n \rangle))$ such that

$$ (e^{-\phi_i^n} - e^{\langle \phi^n \rangle}) (-\phi_i^n + \langle \phi^n \rangle) = e^{\xi_i^n} (-\phi_i^n + \langle \phi^n \rangle)^2. $$

Because of the bounds (87), (88), we deduce that $e^{\xi_i^n} \ge e^{\min \log(\rho_i^n)}$ for $i \in \{0, \dots, N\}$, $\alpha^n$ so that this lower bound is uniform with respect to $i$. We thus obtain

$$ \sum_{i=0}^{N} \varepsilon^2 \left| \frac{\psi_{i+1} - \psi_i}{\Delta x} \right|^2 \Delta x + \alpha^n \sum_{i=0}^{N} \left| -\phi_i^n + \langle \phi^n \rangle \right|^2 \Delta x \le \sum_{i=0}^{N} (\rho_i^n - e^{\langle \phi^n \rangle}) (-\phi_i^n + \langle \phi^n \rangle) \Delta x. \quad (95) $$

Besides, using a Cauchy-Schwarz inequality for the right hand side we deduce the first bound

$$ \alpha^n \|\phi^n - \langle \phi^n \rangle\|_{L^2(\mathbb{T})} \le \|\rho^n - e^{\langle \phi^n \rangle}\|_{L^2(\mathbb{T})} \le \|\rho^n\|_{L^2(\mathbb{T})} + e^{\langle \phi^n \rangle} \quad (96) $$

$$ \le \|\rho^n\|_{L^2(\mathbb{T})} + \|e^{\phi^n}\|_{L^1(\mathbb{T})} \le \|\rho^n\|_{L^2(\mathbb{T})} + \|\rho^n\|_{L^1(\mathbb{T})}. \quad (97) $$

where we used Jensen's inequality and the bound (85). We now prove the last estimate (90). From (95), we have (recalling $\psi_i = -\phi_i^n + \langle \phi^n \rangle$)

$$ \varepsilon^2 |\phi^n - \langle \phi^n \rangle|^2_{H^1(\mathbb{T})} \le \|\rho^n - e^{\langle \phi^n \rangle}\|_{L^2(\mathbb{T})} \|\phi^n - \langle \phi^n \rangle\|_{L^2(\mathbb{T})} \quad (98) $$

$$ \le (\|\rho^n\|_{L^2(\mathbb{T})} + \|\rho^n\|_{L^1(\mathbb{T})}) \|\phi^n - \langle \phi^n \rangle\|_{L^2(\mathbb{T})} \quad (99) $$

$$ \le \frac{1}{\alpha^n} (\|\rho^n\|_{L^2(\mathbb{T})} + \|\rho^n\|_{L^1(\mathbb{T})})^2. \quad (100) $$

where we used (89) for the last inequality. It yields (90). $\square$

For the existence proof, we will need a weaker $H^1$ estimate than (90) for the discrete potential which does not involve the constant $\alpha^n$ given in Proposition 4. In this respect, we will use the discrete analogue of the Poincaré-Wirtinger inequality in $L^2(\mathbb{T})$ that we prove right after.

**Lemma 8.** (Discrete Poincaré-Wirtinger inequality)

$$ \forall u \in X(\mathcal{T}), \quad \left\|u - \langle u \rangle\right\|_{L^2(\mathbb{T})}^2 \le \frac{1}{3} |u|_{H^1(\mathbb{T})}^2, \quad (101) $$

where $\langle u \rangle = \sum_{i=0}^{N} u_i \Delta x$.

*Proof.* The proof mimicks the one in the continuous case. Let $u \in X(\mathcal{T})$. We have

$$ \begin{aligned} \sum_{i=0}^{N} |u_i - \langle u \rangle|^2 \Delta x &= \sum_{i=0}^{N} |u_i - \sum_{j=0}^{N} u_j \Delta x|^2 \Delta x \\ &\le \sum_{i=0}^{N} \left| \sum_{j=0}^{N} (u_i - u_j) \Delta x \right|^2 \Delta x \le \sum_{i=0}^{N} \sum_{j=0}^{N} (u_i - u_j)^2 \Delta x \Delta x \end{aligned} \quad (102) $$
where the last inequality is obtained thanks to Jensen's inequality. Then for $i, j \in \{0, \dots, N\}$ we have, using a Cauchy-Schwarz inequality, that if $i > j$

$$ \left| \frac{u_i - u_j}{\Delta x} \right|^2 = \left| \sum_{k=j}^{i-1} \frac{u_{k+1} - u_k}{\Delta x} \right|^2 \le (i-j) \sum_{k=j}^{i-1} \left| \frac{u_{k+1} - u_k}{\Delta x} \right|^2 \le (i-j) \sum_{k=0}^{N} \left| \frac{u_{k+1} - u_k}{\Delta x} \right|^2. $$

while if $i \le j$

$$ \left| \frac{u_j - u_i}{\Delta x} \right|^2 = \left| \sum_{k=i}^{j-1} \frac{u_{k+1} - u_k}{\Delta x} \right|^2 \le (j-i) \sum_{k=i}^{j-1} \left| \frac{u_{k+1} - u_k}{\Delta x} \right|^2 \le (j-i) \sum_{k=0}^{N} \left| \frac{u_{k+1} - u_k}{\Delta x} \right|^2. $$

Plugging this inequality in (102) we obtain

$$ \sum_{i=0}^{N} \left| u_i - \langle u \rangle \right|^2 \Delta x \le \left( \sum_{k=0}^{N} \left| \frac{u_{k+1} - u_k}{\Delta x} \right|^2 \Delta x \right) \sum_{i=0}^{N} \sum_{j=0}^{N} |(i-j) \Delta x| \Delta x \Delta x. $$

Then a direct computation yields

$$ \sum_{i=0}^{N} \sum_{j=0}^{N} |i - j| \Delta x^3 = \sum_{i=0}^{N} i(i+1) \Delta x^3 = \frac{1}{3} \frac{N(N+1)(N+2)}{(N+1)^3} = \frac{1}{3} \left(1 - \frac{2}{(N+1)^3}\right) $$

which enables to conclude the proof.

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

We are ready to prove Theorem 4.

*Proof.* Consider the assumption of Theorem 4. The existence part is done by an induction argument. So we assume that for a fixed integer $n \in \{0, \dots, N_T - 1\}$ we have been able to construct a solution $(\rho^n, u^n, \phi^n) \in X(\mathcal{T}) \times X(\mathcal{T}^*) \times X(\mathcal{T})$ with $\rho^n > 0$. We want to prove the existence of a solution at step $n+1$. We shall apply the Brouwer fixed-point theorem. Consider the ball of radius $M$ centered at $u^n$

$$ B_M := \{ u \in X(\mathcal{T}^*) : \|u - u^n\|_{L^2(\mathbb{T})} \le M \} $$

where $M > 0$ is to be fixed later. Note that $B_M$ is closed for the $L^2$-topology and that $X(\mathcal{T}^*)$ is a finite dimensional space.

**Definition of a Mapping.** We consider $T: B_M \to X(\mathcal{T}^*)$ which is defined on $B_M$ in three steps. For $u \in B_M$:

* Step 1: Compute $\bar{\rho}(u)$ that solves for $i \in \{0, \dots, N\}$:

$$ \frac{\bar{\rho}_i - \rho_i^n}{\Delta t} + \frac{\mathcal{F}_{i+\frac{1}{2}}(u) - \mathcal{F}_{i-\frac{1}{2}}(u)}{\Delta x} = 0 \qquad (103) $$

where

$$ \mathcal{F}_{i+\frac{1}{2}}(u) = G(\bar{\rho}_i, \bar{\rho}_{i+1}, u_{i+\frac{1}{2}}). \qquad (104) $$

Since the flux (51) is linear in its two first arguments, the equation on $\bar{\rho}$ can be written under the form of a linear system $L(u)\bar{\rho} = \rho^n$ where $L$ is a M-matrix of size $N+1$ given for $i \in \{0, \dots, N\}$ by:

$$ L_{i,i} = 1 + \frac{\Delta t}{\Delta x} (g(u_{i+\frac{1}{2}}) + g(u_{i-\frac{1}{2}}) - u_{i-\frac{1}{2}}), $$

$$ L_{i,i+1} = -\frac{\Delta t}{\Delta x} (g(u_{i+\frac{1}{2}}) - u_{i+\frac{1}{2}}) \mathbf{1}_{i+1 \le N}, $$

$$ L_{i,i-1} = -\frac{\Delta t}{\Delta x} g(u_{i-\frac{1}{2}}) \mathbf{1}_{i-1 \ge 0}. $$
So $\bar{\rho}$ is uniquely defined and since $\rho^n > 0$ we have $\bar{\rho} > 0$. The fact $L$ is a M-matrix comes from $L$ has positive diagonal terms ($L_{i,i} > 0$) and non-positive off-diagonal terms ($L_{i,j} \le 0$ for $i \ne j$) and is strictly diagonally dominant with respect to their columns. For the latter argument, we indeed have $L_{i,i} > \sum_{j \ne i} |L_{j,i}| = |L_{i-1,i}| + |L_{i+1,i}|$ since

$$1 + \frac{\Delta t}{\Delta x} (g(u_{i+\frac{1}{2}}) + g(u_{i-\frac{1}{2}}) - u_{i-\frac{1}{2}}) > \frac{\Delta t}{\Delta x} (g(u_{i-\frac{1}{2}}) - u_{i-\frac{1}{2}} + g(u_{i+\frac{1}{2}})).$$

See [4, 11].

* Step 2: Compute $\varphi = \varphi(\bar{\rho}(u))$ which solves the non linear discrete Poisson equation for $i \in \{0, \dots, N\}$:

$$\varepsilon^2(\Delta\varphi)_i + e^{-\varphi_i} = \bar{\rho}_i.$$

Existence and uniqueness for this non linear equation is classical and can be proven for example by minimization of a strictly convex functional. Mimicking exactly the computation to obtain the estimate (95), we get using a Cauchy-Schwarz inequality combined with the Poincaré-Wirtinger inequality (101),

$$|\varphi|_{H^1(\mathbb{T})} \le \frac{1}{\varepsilon^2\sqrt{3}} \|\bar{\rho} - e^{\langle\varphi\rangle}\|_{L^2(\mathbb{T})},$$

and using (96), we obtain

$$|\varphi|_{H^1(\mathbb{T})} \le \frac{1}{\varepsilon^2\sqrt{3}} \|\bar{\rho} - e^{\langle\varphi\rangle}\|_{L^2(\mathbb{T})} \le \frac{1}{\varepsilon^2\sqrt{3}} (\|\bar{\rho}\|_{L^2(\mathbb{T})} + \|\bar{\rho}\|_{L^1(\mathbb{T})}) \le C \left( \Delta x, \|\rho^n\|_{L^1(\mathbb{T})}, \frac{1}{\varepsilon^2} \right), \quad (105)$$

where $C(\Delta x, \|\rho^n\|_{L^1(\mathbb{T})}, \frac{1}{\varepsilon^2}) > 0$ is constant that depends only on $\Delta x$ and $\|\rho^n\|_{L^1(\mathbb{T})}$. This constant is obtained thanks to the equivalence of norms in finite dimension, the positivity of $\bar{\rho}$ and the conservation of mass given by (103).

* Step 3: Compute $v = v(\bar{\rho}(u), \varphi(\bar{\rho}(u)), u)$ which solves for $i \in \{0, \dots, N\}$:

$$\frac{\bar{\rho}_{i+\frac{1}{2}} v_{i+\frac{1}{2}} - \rho_{i+\frac{1}{2}}^n u_{i+\frac{1}{2}}^n}{\Delta t} + \frac{Q_{i+1}(u) v_{i+1} - Q_i(u) v_i}{\Delta x} = \bar{\rho}_{i+\frac{1}{2}} \delta(\varphi)_{i+\frac{1}{2}}, \quad (106)$$

where

$$\bar{\rho}_{i+\frac{1}{2}} = \frac{\bar{\rho}_i + \bar{\rho}_{i+1}}{2}, \quad Q_i(u) = \frac{\mathcal{F}_{i+\frac{1}{2}}(u) + \mathcal{F}_{i-\frac{1}{2}}(u)}{2}, \quad v_i = \begin{cases} v_{i-\frac{1}{2}} & \text{if } Q_i(u) \ge 0, \\ v_{i+\frac{1}{2}} & \text{if } Q_i(u) < 0 \end{cases} \quad (107)$$

and

$$\bar{\rho}_{i+\frac{1}{2}} = \begin{cases} \frac{G(\bar{\rho}_i, \bar{\rho}_{i+1}, u_{i+\frac{1}{2}}) - G(\bar{\rho}_i, \bar{\rho}_{i+1}, 0)}{u_{i+\frac{1}{2}}} & \text{if } u_{i+\frac{1}{2}} \ne 0, \\ \bar{\rho}_{i+1} - (\bar{\rho}_{i+1} - \bar{\rho}_i)g'(0) & \text{if } u_{i+\frac{1}{2}} = 0. \end{cases} \quad (108)$$

Note that $v$ solves a linear system which is invertible.

Step 1, Step 2 and Step 3 are well defined, so is $T$ on $B_M$. $T$ is moreover continuous on $B_M$ notably because the flux given in (104) is continuous, the flux part in (124) is also continuous and the forcing term in (124) has been designed in such a way that (108) is in particular continuous with respect to $u$.

**Stability of $B_M$**. We want to prove that $T(B_M) \subset B_M$ for a well-chosen $M > 0$. We first mimick the energy estimates as in Proposition 3. Note that because the right hand side in (106) is now explicit, we do not have the energy decay. We have instead, after multiplying (106) by $v_{i+1/2}$

$$\sum_{i=0}^{N} \rho_{i+\frac{1}{2}}^{n} \frac{(v_{i+\frac{1}{2}} - u_{i+\frac{1}{2}}^{n})^{2}}{2} \Delta x + \mathcal{H}(\bar{\rho}, v, \varphi) \le \mathcal{H}(\rho^{n}, u^{n}, \phi^{n}) + \Delta t \sum_{i=0}^{N} \bar{\rho}_{i+\frac{1}{2}} (v_{i+\frac{1}{2}} - u_{i+\frac{1}{2}}) \delta(\varphi)_{i+\frac{1}{2}} \Delta x. \quad (109)$$
The first term comes from the implicitness of the discrete time derivative (its equivalent, is the third term in (78)). We need to estimate the residual term $P = \Delta t \sum_{i=0}^{N} \tilde{\rho}_{i+\frac{1}{2}}(v_{i+\frac{1}{2}} - u_{i+\frac{1}{2}})(\delta\phi)_{i+\frac{1}{2}}\Delta x$. Using the definition of $G$ given in (51) and the fact that $g$ is Lipschitz continuous, we have for $i \in \{0, ..., N\}$:

$$|\tilde{\rho}_{i+\frac{1}{2}}| \le (1 + 2\mathrm{Lip}(g)) \|\bar{\rho}\|_{L^{\infty}(\mathbb{T})} \le \frac{(1 + 2\mathrm{Lip}(g))}{\Delta x} \|\bar{\rho}\|_{L^1(\mathbb{T})} \le \frac{(1 + 2\mathrm{Lip}(g))}{\Delta x} \|\rho^n\|_{L^1(\mathbb{T})}, \quad (110)$$

where we used the equivalence of norms in finite dimension, the conservation of the positivity and the total mass given by (103). Using a Hölder inequality and the elliptic estimate (105) we obtain

$$|P| \le \Delta t C' \left( \Delta x, \mathrm{Lip}(g), \|\rho^n\|_{L^1(\mathbb{T})}, \frac{1}{\varepsilon^2} \right) \|v - u\|_{L^2(\mathbb{T})} \quad (111)$$

where $C'(\Delta x, \mathrm{Lip}(g), \|\rho^n\|_{L^1(\mathbb{T})}, \frac{1}{\varepsilon^2}) > 0$ is a constant that depends only on $\Delta x, \varepsilon, \mathrm{Lip}(g)$ and $\|\rho^n\|_{L^1(\mathbb{T})}$. We set $C^n := C'(\Delta x, \mathrm{Lip}(g), \|\rho^n\|_{L^1(\mathbb{T})}, \frac{1}{\varepsilon^2})$ in the remaining part of the proof. Using a triangular inequality and the fact that $u \in B_M$, we have the following estimate for $P$,

$$|P| \le \Delta t C^n \|v - u^n\|_{L^2(\mathbb{T})} + \Delta t C^n M. \quad (112)$$

Besides, we observe that the discrete energy functional is bounded below,

$$\mathcal{H}(\bar{\rho}, v, \varphi) \ge \bar{h}, \quad \text{where} \quad \bar{h} = \min_{\mathbb{R}} h > -\infty. \quad (113)$$

Combining (112) and (113) with (109), we obtain that

$$\frac{\min_{i \in \{0, \dots, N\}} \rho_{i+\frac{1}{2}}^n}{2} \|v - u^n\|_{L^2(\mathbb{T})}^2 - \Delta t C^n \|v - u^n\|_{L^2(\mathbb{T})} - (\Delta t C^n M + \mathcal{H}(\rho^n, u^n, \phi^n) - \bar{h}) \le 0, \quad (114)$$

where we recall that $\mathcal{H}(\rho^n, u^n, \phi^n) - \bar{h} \ge 0$ since $\sum_{i=0}^N (h(\phi_i^n) - \bar{h})\Delta x \ge 0$ by definition of $\bar{h}$. Note that (114) is a polynomial of second degree in $\|v - u^n\|$ and the inequality (114) tells us that this polynomial is non positive on $\mathbb{R}^+$. It has two roots of opposite sign and the non negative root is given by

$$X = \frac{\Delta t C^n + \sqrt{(\Delta t C^n)^2 + 2(\Delta t C^n M + \mathcal{H}(\rho^n, u^n, \phi^n) - \bar{h}) \min_{i \in \{0, \dots, N\}} \rho_{i+\frac{1}{2}}^n}}{\min_{i \in \{0, \dots, N\}} \rho_{i+\frac{1}{2}}^n}.$$

Note that that the denominator is fixed while the numerator behaves as $O(\sqrt{M})$ as $M \to +\infty$. Therefore, we claim that there exists $M > 0$ large enough, which possibly depends on $(\rho^n, u^n, \phi^n, \Delta x, \Delta t, \varepsilon)$ such that $X \le M$. Thus (114) implies that $\|v - u^n\|_{L^2(\mathbb{T})} \le X \le M$. We have proven that $T(B_M) \subset B_M$. The Brouwer-fixed point theorem thus applies. By induction we deduce the existence of a solution for $n \in \{0, \dots, N_T\}$. The estimate (64) is a consequence of the positivity of the scheme. The estimate (65) is obtained by summation of the discrete continuity equation and the positivity of $\rho^n$. The energy decay (66) is just a direct application of Proposition 3. Finally the discrete elliptic estimates are a consequence of Proposition 4. $\square$

## 2.3 Discrete modulated energy estimate

This section is devoted to the proof of Theorem 5. The first step consists in establishing the discrete evolution law for the modulated energy (44).

**Lemma 9.** (Evolution of the discrete modulated energy) Let $(\rho^n, u^n, \phi^n)_{n=0, \dots, N_T} \subset X(\mathcal{T}) \times X(\mathcal{T}^*) \times X(\mathcal{T})$ a solution to (48) such that $\rho^n > 0$. Assume $\phi^0 \in X(\mathcal{T})$ verifies the discrete Poisson equation initially. Let $(\bar{u}, \bar{\phi}) \in X(\mathcal{T}^*) \times X(\mathcal{T})$ a constant state. Then the solution satisfies for $0 \le n \le N_T - 1$,

$$\mathcal{E}(\rho^{n+1}, u^{n+1}, \phi^{n+1}|\bar{u}, \bar{\phi}) = \mathcal{E}(\rho^n, u^n, \phi^n|\bar{u}, \bar{\phi}) - \Delta t \tau^{n+1} - \Delta t \bar{u} \sum_{i=0}^{N} \tilde{\rho}_{i+\frac{1}{2}}^{n+1} \delta(\phi^{n+1})_{i+\frac{1}{2}} \Delta x, \quad (115)$$

where we recall that $\tau^{n+1} \ge 0$ is given by (83).
Proof. Let $0 \le n \le N_T - 1$. We set $\bar{\rho} = e^{-\bar{\phi}}$. Thanks to the decomposition (45) we have

$$
\begin{align*}
\mathcal{E}(\rho^{n+1}, u^{n+1}, \phi|\bar{u}, \bar{\phi}) &= \mathcal{H}(\rho^{n+1}, u^{n+1}, \phi^{n+1}) + \mathcal{E}_{\text{kin}}(\rho^{n+1}, u^{n+1}|\bar{u}) - \mathcal{E}_{\text{int}}(\phi^{n+1}|\bar{\phi}) \\
&= \mathcal{H}(\rho^n, u^n, \phi^n) - \Delta t \tau^{n+1} + \mathcal{E}_{\text{kin}}(\rho^{n+1}, u^{n+1}|\bar{u}) - \mathcal{E}_{\text{int}}(\phi^{n+1}|\bar{\phi}),
\end{align*}
$$

where we used the energy decay (66). Besides, a direct computation gives

$$
\mathcal{E}_{\text{int}}(\phi^{n+1}|\bar{\phi}) = \sum_{i=0}^{N} \left( \tilde{h}(\bar{\rho}) + \tilde{h}'(\bar{\rho}) (e^{-\phi_i^{n+1}} - \bar{\rho}) \right) \Delta x.
$$

Recall that $\bar{\rho}$ is constant. So, using the discrete Poisson we have $e^{-\phi_i^{n+1}} - \bar{\rho} = \rho_i^{n+1} - \bar{\rho} - \varepsilon^2 \Delta(\phi^{n+1})_i$.
So with the periodicity,

$$
\mathcal{E}_{\text{int}}(\phi^{n+1}|\bar{\phi}) = \sum_{i=0}^{N} \left( \tilde{h}(\bar{\rho}) + \tilde{h}'(\bar{\rho}) (\rho_i^{n+1} - \bar{\rho}) \right) \Delta x = \sum_{i=0}^{N} \left( \tilde{h}(\bar{\rho}) + \tilde{h}'(\bar{\rho}) (\rho_i^n - \bar{\rho}) \right) \Delta x = \mathcal{E}_{\text{int}}(\phi^n|\bar{\phi}),
$$

where the last equality is obtained thanks to the mass conservation. We compute the evolution of the
modulated kinetic energy. A direct computation yields

$$
\frac{\mathcal{E}_{\text{kin}}(\rho^{n+1}, u^{n+1}|\bar{u}) - \mathcal{E}_{\text{kin}}(\rho^n, u^n|\bar{u})}{\Delta t} = \sum_{i=0}^{N} \frac{\bar{u}^2}{2} \left( \frac{\rho_{i+\frac{1}{2}}^{n+1} - \rho_{i+\frac{1}{2}}^n}{\Delta t} \right) \Delta x - \sum_{i=0}^{N} \bar{u} \left( \frac{\rho_{i+\frac{1}{2}}^{n+1} u_{i+\frac{1}{2}}^{n+1} - \rho_{i+\frac{1}{2}}^n u_{i+\frac{1}{2}}^n}{\Delta t} \right) \Delta x.
$$

Using the discrete continuity equation and the discrete momentum we obtain

$$
\frac{\mathcal{E}_{\text{kin}}(\rho^{n+1}, u^{n+1}|\bar{u}) - \mathcal{E}_{\text{kin}}(\rho^n, u^n|\bar{u})}{\Delta t} = -\bar{u} \sum_{i=0}^{N} \bar{\rho}_{i+\frac{1}{2}}^{n+1} \delta(\phi^{n+1})_{i+\frac{1}{2}} \Delta x.
$$

Then gathering all the terms together, we obtain the expected equality.

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

To obtain a close estimate for the equation (115), we need to control the term

$$
W^{n+1} = \sum_{i=0}^{N} \bar{\rho}_{i+\frac{1}{2}}^{n+1} \delta(\phi^{n+1})_{i+\frac{1}{2}} \Delta x. \qquad (116)
$$

In this regard we have.

**Lemma 10.** (Control on W) Let $(\rho^n, u^n, \phi^n)_{n=0,...,N_T} \subset X(\mathcal{T}) \times X(\mathcal{T}^*) \times X(\mathcal{T})$ a solution to (48) such that $\rho^n > 0$. Assume $\phi^0 \in X(\mathcal{T})$ verifies the discrete Poisson equation initially. Then we have for $n \in \{0, ..., N_T - 1\}$,

$$
W^{n+1} = W_1^{n+1} + W_2^{n+1}, \tag{117}
$$

$$
W_1^{n+1} = - \sum_{i=0}^{N} (\rho_{i+1}^{n+1} - \rho_i^{n+1}) \hat{g}(u_{i+\frac{1}{2}}^{n+1}) \delta(\phi^{n+1})_{i+\frac{1}{2}} \Delta x, \quad (118)
$$

$$
W_2^{n+1} = \sum_{i=0}^{N} \rho_{i+1}^{n+1} \delta(\phi^{n+1})_{i+\frac{1}{2}} \Delta x, \qquad (119)
$$

where $\hat{g}(u) = \begin{cases} \frac{g(u) - g(0)}{u} & \text{if } u \neq 0, \\ g'(0) & \text{if } u = 0. \end{cases}$ Moreover, we have the following estimates

$$
|W_1^{n+1}| \le \frac{8\text{Lip}(g)}{\Delta x} \mathcal{E}(\rho^{n+1}, u^{n+1}, \phi^{n+1}|\bar{u}, \bar{\phi}) + \frac{\text{Lip}(g)}{\Delta x} \sum_{i=0}^{N} \left| (\phi_{i+1}^n - \phi_i^n)(e^{\phi_{i+1}^n} - e^{-\phi_i^n}) \right| \Delta x, \quad (120)
$$

$$
|W_2^{n+1}| \le \frac{4}{\Delta x} \mathcal{E}(\rho^{n+1}, u^{n+1}, \phi^{n+1}|\bar{u}, \bar{\phi}) + \frac{2\Delta x}{\varepsilon^2} \|\rho^{n+1}\|_{L^\infty(\mathbb{T})} \mathcal{E}(\rho^{n+1}, u^{n+1}, \phi^{n+1}|\bar{u}, \bar{\phi}). \quad (121)
$$
Proof. Let $0 \le n \le N_T - 1$. Using the definitions (58), (51) and the definition of $\hat{g}$, $W^{n+1}$ can be written as

$$W^{n+1} = - \sum_{i=0}^{N} (\rho_{i+1}^{n+1} - \rho_i^{n+1}) \hat{g}(u_{i+\frac{1}{2}}^{n+1}) \delta(\phi^{n+1})_{i+\frac{1}{2}} \Delta x + \sum_{i=0}^{N} \rho_{i+1}^{n+1} \delta(\phi^{n+1})_{i+\frac{1}{2}} \Delta x$$

Let observe that $\hat{g}$ is a continuous and bounded function with $\|\hat{g}\|_{L^\infty(\mathbb{R})} \le \text{Lip}(g)$. It shows the decomposition (117). We now estimate each term separately. Using the discrete Poisson equation, we decompose the first term as

$$W_1^{n+1} = W_{1,1}^{n+1} + W_{1,2}^{n+1},$$

$$W_{1,1}^{n+1} = - \sum_{i=0}^{N} \varepsilon^2 (\Delta(\phi^{n+1})_{i+1} - \Delta(\phi^{n+1})_i) \delta(\phi^{n+1})_{i+\frac{1}{2}} \hat{g}(u_{i+\frac{1}{2}}^{n+1}) \Delta x,$$

$$W_{1,2}^{n+1} = - \sum_{i=0}^{N} \left( e^{-\phi_{i+1}^{n+1}} - e^{-\phi_i^{n+1}} \right) \delta(\phi^{n+1})_{i+\frac{1}{2}} \hat{g}(u_{i+\frac{1}{2}}^{n+1}) \Delta x.$$

To estimate $W_{1,1}^{n+1}$, we use a discrete integration by parts then we get

$$W_{1,1}^{n+1} = \sum_{i=0}^{N} \varepsilon^2 \Delta(\phi^{n+1})_i (\hat{g}(u_{i+\frac{1}{2}}^{n+1}) \delta(\phi^{n+1})_{i+\frac{1}{2}} - \hat{g}(u_{i-\frac{1}{2}}^{n+1}) \delta(\phi^{n+1})_{i-\frac{1}{2}}) \Delta x \\ = \frac{1}{\Delta x} \sum_{i=0}^{N} \varepsilon^2 (\delta(\phi^{n+1})_{i+\frac{1}{2}} - \delta(\phi^{n+1})_{i-\frac{1}{2}}) (\hat{g}(u_{i+\frac{1}{2}}^{n+1}) \delta(\phi^{n+1})_{i+\frac{1}{2}} - \hat{g}(u_{i-\frac{1}{2}}^{n+1}) \delta(\phi^{n+1})_{i-\frac{1}{2}}) \Delta x.$$

Expanding the product, using a Young inequality and a translation of indices, we obtain

$$|W_{1,1}^{n+1}| \le \frac{8\text{Lip}(g)}{\Delta x} \sum_{i=0}^{N} \frac{\varepsilon^2}{2} |\delta(\phi^{n+1})_{i+\frac{1}{2}}|^2 \Delta x \le \frac{8\text{Lip}(g)}{\Delta x} \mathcal{E}(\rho^{n+1}, \phi^{n+1}, u^{n+1} |\bar{u}, \bar{\phi}). \quad (122)$$

As for $W_{1,2}^{n+1}$, we obtain readily

$$|W_{1,2}^{n+1}| \le \frac{\text{Lip}(g)}{\Delta x} \sum_{i=0}^{N} |(\phi_{i+1}^{n} - \phi_{i}^{n})(e^{\phi_{i+1}^{n}} - e^{-\phi_{i}^{n}})| \Delta x.$$

We thus infer the estimate (120).

We now treat the term $W_2^{n+1}$. We first use the discrete Poisson equation to get

$$W_2^{n+1} = \sum_{i=0}^{N} \varepsilon^2 \Delta(\phi^{n+1})_{i+1} \delta(\phi^{n+1})_{i+\frac{1}{2}} \Delta x + \sum_{i=0}^{N} e^{-\phi_{i+1}^{n+1}} \delta(\phi^{n+1})_{i+\frac{1}{2}} \Delta x.$$

We set

$$W_{2,1}^{n+1} = \sum_{i=0}^{N} \varepsilon^2 \Delta(\phi^{n+1})_{i+1} \delta(\phi^{n+1})_{i+\frac{1}{2}} \Delta x \quad \text{and} \quad W_{2,2}^{n+1} = \sum_{i=0}^{N} e^{-\phi_{i+1}^{n+1}} \delta(\phi^{n+1})_{i+\frac{1}{2}} \Delta x.$$

We are going to bound each term separately. We have,

$$W_{2,1}^{n+1} = \sum_{i=0}^{N} \frac{\varepsilon^2}{\Delta x} \left( \delta(\phi^{n+1})_{i+\frac{3}{2}} - \delta(\phi^{n+1})_{i+\frac{1}{2}} \right) \delta(\phi^{n+1})_{i+\frac{1}{2}} \Delta x.$$

Using a Young inequality and a translation of indices we obtain the following bound for $W_{2,1}^{n+1}$,

$$|W_{2,1}^{n+1}| \le \frac{2}{\Delta x} \sum_{i=0}^{N} \varepsilon^2 |\delta(\phi^{n+1})_{i+\frac{1}{2}}|^2 \Delta x \le \frac{4}{\Delta x} \mathcal{E}(\rho^{n+1}, \phi^{n+1}, u^{n+1} |\bar{u}, \bar{\phi}).$$
As for the term $W_{2,2}^{n+1}$, we use a Taylor-Lagrange expansion: for each $i \in \{0, \dots, N\}$ there exists $\zeta_i^{n+1} \in (\min(\phi_{i+1}^{n+1}, \phi_i^{n+1}), \max(\phi_{i+1}^{n+1}, \phi_i^{n+1}))$

$$e^{-\phi_i^{n+1}} = e^{-\phi_{i+1}^{n+1}} - e^{-\phi_{i+1}^{n+1}}(\phi_{i+1}^{n+1} - \phi_i^{n+1}) + e^{-\zeta_i^{n+1}} \frac{(\phi_{i+1}^{n+1} - \phi_i^{n+1})^2}{2}.$$

Using the periodicity, we obtain after summation that

$$W_{2,2}^{n+1} = -\frac{1}{\Delta x} \sum_{i=0}^{N} e^{-\zeta_i^{n+1}} \frac{(\phi_{i+1}^{n+1} - \phi_i^{n+1})^2}{2} \Delta x = -\Delta x \sum_{i=0}^{N} e^{-\zeta_i^{n+1}} |\delta(\phi^{n+1})_{i+\frac{1}{2}}|^2 \Delta x.$$

Using the maximum principle (87) we obtain the bound $e^{-\xi_i^{n+1}} \le \|\rho^{n+1}\|_{L^\infty(\mathbb{T})}$ for $i \in \{0, \dots, N\}$, so that we deduce

$$|W_{2,2}^{n+1}| \le \frac{2\Delta x}{\varepsilon^2} \|\rho^{n+1}\|_{L^\infty(\mathbb{T})} \mathcal{E}(\rho^{n+1}, \phi^{n+1}, u^{n+1}|\bar{u}, \bar{\phi}).$$

So that eventually,

$$|W_2^{n+1}| \le \frac{4}{\Delta x} \mathcal{E}(\rho^{n+1}, \phi^{n+1}, u^{n+1}|\bar{u}, \bar{\phi}) + \frac{2\Delta x}{\varepsilon^2} \|\rho^{n+1}\|_{L^\infty(\mathbb{T})} \mathcal{E}(\rho^{n+1}, \phi^{n+1}, u^{n+1}|\bar{u}, \bar{\phi})$$

which is the expected estimate (121)

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

We are ready to prove Theorem 5.

*Proof.* Combining the estimate of Lemma 10, with the evolution of the modulated energy (115), we obtain the following following closed estimate for the modulated energy:

$$\begin{aligned} \forall n \in \{0, \dots, N_T - 1\}, \quad & \mathcal{E}(\rho^{n+1}, u^{n+1}, \phi^{n+1}|\bar{u}, \bar{\phi}) \le \mathcal{E}(\rho^n, u^n, \phi^n|\bar{u}, \bar{\phi}) \tag{123} \\ & + |\bar{u}| \frac{\Delta t}{\Delta x} (8\mathrm{Lip}(g) + 4 + \frac{2\Delta x^2}{\varepsilon^2} \|\rho^{n+1}\|_{L^\infty(\mathbb{T})}) \mathcal{E}(\rho^{n+1}, u^{n+1}, \phi^{n+1}|\bar{u}, \bar{\phi}) \\ & - \Delta t \tau^{n+1} + \frac{\Delta t |\bar{u}| \mathrm{Lip}(g)}{\Delta x} \sum_{i=0}^{N} \left| (\phi_{i+1}^n - \phi_i^n)(e^{\phi_{i+1}^n} - e^{-\phi_i^n}) \right| \Delta x. \end{aligned}$$

Now remark that in $\Delta t \tau^{n+1}$ (83) we have the term

$$\frac{\Delta t g(0)}{\Delta x} \sum_{i=0}^{N} \left| (e^{-\phi_i^{n+1}} - e^{-\phi_{i-1}^{n+1}})(\phi_i^{n+1} - \phi_{i-1}^{n+1}) \right| \Delta x.$$

So we see that the last term in (123) can be absorbed by $-\Delta t \tau^{n+1}$ if $|\bar{u}| \le \frac{g(0)}{\mathrm{Lip}(g)}$. The case $\bar{u} = 0$ is trivial. It completes the proof of Theorem 5.

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

# 3 Numerical experiment

This section is dedicated to some numerical experiments. First, we describe how we solve the non linear scheme. Then, a test illustrating the theorems proven above (namely Theorem 4 and Theorem 5) is performed and we explore the accuracy of our scheme. Finally, a test taken from [18, 7] is considered.

## 3.1 Description of the non linear solver

Since we shall use a Newton method, we assume that the function $g$ verifies, in addition to (52)-(54), $g \in C^1(\mathbb{R})$ and that it is twice differentiable near the origin. Suppose now that for a fixed integer $n \in \{0, \dots, N_T - 1\}$ we have constructed the solution $(\rho^n, u^n, \phi^n) \in X(\mathcal{T}) \times X(\mathcal{T}^*) \times X(\mathcal{T})$. We seek a solution of the scheme (48) at step $n+1$ as a fixed point of a certain map just exactly as in the proof of existence. More precisely, we look for $u^{n+1} \in X(\mathcal{T}^*)$ which solves

$$\mathbf{T}^n(u^{n+1}) = u^{n+1}$$

where $\mathbf{T}^n : X(\mathcal{T}^*) \to X(\mathcal{T}^*)$ is defined for $u \in X(\mathcal{T}^*)$ in three steps:
• Step 1: we solve the linear continuity equation of unknown $\bar{\rho} \equiv \bar{\rho}(u)$ which solves (103).

• Step 2: we solve the non linear Poisson equation of unknown $\varphi \equiv \varphi(\bar{\rho}) \in X(\mathcal{T})$ which verifies for $i \in \{0, ..., N\}$:

$$ \varepsilon^2(\Delta\varphi)_i + e^{-\varphi_i} = \bar{\rho}_i. $$

We use a Newton-method with an error bound fixed to the zero machine which is in our case $10^{-15}$.

• Step 3: we solve the non linear momentum equation for $v \equiv v(\bar{\rho}, u) \in X(\mathcal{T}^*)$ which verifies for $i \in \{0, ..., N\}$:

$$ \frac{\bar{\rho}_{i+\frac{1}{2}} v_{i+\frac{1}{2}} - \rho_{i+\frac{1}{2}}^n u_{i+\frac{1}{2}}^n}{\Delta t} + \frac{Q_{i+1}(u)v_{i+1} - Q_i(u)v_i}{\Delta x} = \bar{\rho}(v)_{i+\frac{1}{2}}(\delta\varphi)_{i+\frac{1}{2}}, \quad (124) $$

where

$$ \bar{\rho}_{i+\frac{1}{2}} = \frac{\bar{\rho}_i + \bar{\rho}_{i+1}}{2}, \quad Q_i(u) = \frac{\mathcal{F}_{i+\frac{1}{2}}(u) + \mathcal{F}_{i-\frac{1}{2}}(u)}{2}, \quad v_i = \begin{cases} v_{i-\frac{1}{2}} & \text{if } Q_i(u) \ge 0, \\ v_{i+\frac{1}{2}} & \text{if } Q_i(u) < 0, \end{cases} \quad (125) $$

and

$$ \bar{\rho}(v)_{i+\frac{1}{2}} = \begin{cases} \frac{G(\bar{\rho}_i, \bar{\rho}_{i+1}, v_{i+\frac{1}{2}}) - G(\bar{\rho}_i, \bar{\rho}_{i+1}, 0)}{v_{i+\frac{1}{2}}} & \text{if } v_{i+\frac{1}{2}} \neq 0, \\ \bar{\rho}_{i+1} - (\bar{\rho}_{i+1} - \bar{\rho}_i)g'(0) & \text{if } v_{i+\frac{1}{2}} = 0. \end{cases} \quad (126) $$

We also use a Newton-method with an error bound fixed to the zero machine.

Once these three steps are accomplished we consider that we have computed (an approximation of) $\mathbf{T}^n(u)$ for a given $u \in X(\mathcal{T}^*)$. We thus use a Picard-iteration scheme which consists in the sequence $(u_k^n)_{k \in \mathbb{N}}$ defined by induction as follows:

$$ \begin{cases} u_0^n = u^n, \\ u_{k+1}^n = \mathbf{T}^n(u_k^n), \quad k \in \mathbb{N}. \end{cases} \quad (127) $$

If the sequence $(u_k^n)_{k \in \mathbb{N}}$ converges to some $u_*^n \in X(\mathcal{T}^*)$ then, since $\mathbf{T}^n$ is a continuous map, the limit verifies $\mathbf{T}^n(u_*^n) = u_*^n$ which exactly means that $u_*^n$ is a solution of the non linear system (48). Our stopping criterion for the algorithm is (provided $u_k^n \neq 0$ for each $k$ and the sequence converges to some $u_*^n \neq 0$):

• Compute:

$$ N_* := \inf \left\{ k \in \mathbb{N} : \frac{\|u_{k+1}^n - u_k^n\|_{L^\infty(\mathcal{T})}}{\|u_k^n\|_{L^\infty(\mathcal{T})}} \le 10^{-7} \right\}. $$

Since the sequence is assumed to converge towards a non zero limit this number is well-defined.

• Update the approximate solution by setting:

$$ u^{n+1} = u_{N_*}^n, \quad (128) $$

$$ \rho^{n+1} = \bar{\rho}(u_{N_*}), \quad (129) $$

$$ \phi^{n+1} = \varphi(\rho^{n+1}). \quad (130) $$

Of course, this is an approximation of a fixed point up to the threshold error.

In the following numerical experiment we fix the function $g$ to be given by

$$ g(u) = \begin{cases} u & \text{if } u \ge \Delta x, \\ \frac{(u+\Delta x)^2}{4\Delta x} & \text{if } -\Delta x < u < \Delta x, \\ 0 & \text{if } u \le -\Delta x. \end{cases} \quad (131) $$
## 3.2 Non linear stability around constant states

We consider a constant state of the form

$$ \begin{cases} \bar{u} \in \mathbb{R}, \\ \bar{\phi} = 0, \\ \bar{\rho} = e^{-\bar{\phi}} = 1. \end{cases} \qquad (132) $$

This constant state is clearly a stationary solution of both (1) and (3). We consider a fluctuation around the constant state $(\bar{\rho}, \bar{u})$ in the form

$$ \rho_{\varepsilon}^{\text{ini}}(x) - \bar{\rho} = \frac{\varepsilon^s}{2} \sin(2\pi x \lfloor \varepsilon^{-1} \rfloor), \quad u_{\varepsilon}^{\text{ini}}(x) - \bar{u} = \varepsilon \sin(2\pi x), \quad x \in [0, 1], \qquad (133) $$

where $\varepsilon \in (0, 1]$ and $s \ge 0$. In particular, we see that the fluctuation around $\bar{\rho}$ oscillates at the spatial scale $\varepsilon$. Observe besides that for $0 < s' < s$ we have

$$ \| \rho_{\varepsilon}^{\text{ini}} - \bar{\rho} \|_{H^{s'}(\mathbb{T})} \to 0 \text{ as } \varepsilon \to 0. \qquad (134) $$

In [20], it is shown that provided $s$ is large enough (134) is propagated on $[0, T]$ for a Sobolev exponent which is smaller than $s'$. The modulated energy estimate (20) enables to show the convergence in $L^2(\mathbb{T})$ on $[0, T]$ provided the initial data is such that $\mathcal{E}(0) \to 0$ as $\varepsilon \to 0$ (see [13]). Its discrete analogue is (67) for constant states with $\bar{u} = 0$. The initial data is discretized in a finite volume manner, that is

$$ \rho_i^0 - \bar{\rho} = \frac{\varepsilon^s}{4\pi\Delta x[\varepsilon^{-1}]} \left( \cos\left(2\pi x_{i-\frac{1}{2}} \lfloor \varepsilon^{-1} \rfloor\right) - \cos\left(2\pi x_{i+\frac{1}{2}} \lfloor \varepsilon^{-1} \rfloor\right) \right), \quad i \in \{0, \dots, N\} \qquad (135) $$

$$ u_i^0 - \bar{u} = \frac{\varepsilon}{2\pi\Delta x} \left( \cos(2\pi x_{i-\frac{1}{2}}) - \cos(2\pi x_{i+\frac{1}{2}}) \right), \quad i \in \{0, \dots, N\}. \qquad (136) $$

The initial potential $-\phi^0$ satisfies the discrete nonlinear Poisson equation:

$$ \varepsilon^2(\Delta\phi^0)_i + e^{-\phi^0} = \rho_i^0, \quad i \in \{0, \dots, N\}. $$

### 3.2.1 Convergence as $\varepsilon \to 0$ with a fixed mesh-size

a) A well-prepared data on a coarse mesh. The mesh size is $\Delta x = 10^{-2}$ and the time step is $\Delta t = \frac{1}{2}\Delta x$. The final time is $T = 1000\Delta t$. In Table 1, we report the values of the modulated energy at initial and final time for different values of $\varepsilon$ for an initial data of the form (135)-(136) with $s = 1$.

<table><caption>Table 1: Modulated energy at final and initial time for different values of $\varepsilon$ for an initial data of the form (135)-(136) with $s = 1$.</caption><thead><tr><th>ε</th><th>E(ρ<sup>N<sub>T</sub></sup>, u<sup>N<sub>T</sub></sup>, φ<sup>N<sub>T</sub></sup>|ū, φ̄)</th><th>E(ρ<sup>0</sup>, u<sup>0</sup>, φ<sup>0</sup>|ū, φ̄)</th></tr></thead><tbody><tr><td>0.1</td><td>0.00147172</td><td>0.0225411</td></tr><tr><td>0.01</td><td>1.61869e-05</td><td>0.000224353</td></tr><tr><td>0.001</td><td>1.61757e-07</td><td>2.24352e-06</td></tr><tr><td>0.0001</td><td>1.61787e-09</td><td>2.24348e-08</td></tr></tbody></table>

We observe that whatever the value of $\varepsilon$ is, the modulated energy at final time is lower than the modulated energy at initial time. It is an expected behavior of our scheme. Moreover, we see that when $\varepsilon$ decreases to zero, the modulated energy also decreases towards zero. We measure a convergence rate towards zero in $\varepsilon$ which is $O(\varepsilon^2)$. It is exactly the same rate as the rate of decrease towards zero of the modulated energy at initial time. During the simulation, we have checked the total energy decay, the mass conservation and the conservation of positivity of the density. Note that the time step and the mesh size are fixed and completely independent of $\varepsilon$. These results are in good agreement with Theorem 4 and the item a) of Theorem 5. We have performed the same test with $\bar{u} \in \{-4, -2, 2, 4\}$ and we have obtained comparable results. These results illustrate the unconditional stability of our scheme. However, the results must be interpreted with care since it is only a rough illustration of the convergence as $\varepsilon \to 0$ on a coarse mesh. We do not claim that our scheme is accurate when $\varepsilon \to 0$.
<figure><img src="image_7.png" /><figcaption>Figure 1: Left: time evolution of the modulated energy on $[0, T]$ with $\varepsilon = 10^{-1}$ for three values of $\Delta x$. Right: density $\rho(T = 0.2, \cdot)$ with $\varepsilon = 10^{-1}$ for three values of $\Delta x$.</figcaption></figure>

b) Evaluation of the numerical dissipation rate when $\varepsilon$ is fixed. We quantify the numerical dissipation when $\varepsilon$ is fixed but $\Delta x$ and $\Delta t$ tend both to zero. We expect that for smooth fluctuations, the numerical dissipation rate should tend towards zero since in the continuous case the energy of smooth solutions is conserved. The numerical parameters are: $\varepsilon = 10^{-1}$, $\Delta t = \frac{1}{2}\Delta x$. The final time is $T = 0.2$. The initial data is of the form (135)-(136) with $s = 1$. In Table 2, we report the numerical dissipation rate defined by

$$ \tau(\Delta x) = \frac{\log(\mathcal{E}(\rho^{N_T}, u^{N_T}, \phi^{N_T}|\bar{u}, \bar{\phi})) - \log(\mathcal{E}(\rho^0, u^0, \phi^0|\bar{u}, \bar{\phi}))}{T} \quad (137) $$

We see that with the CFL condition $\Delta t = \frac{1}{2}\Delta x$ the numerical dissipation rate is of order one in $\Delta x$.

<table><caption>Table 2: Numerical dissipation rate for $\varepsilon = 10^{-1}$ for three values of $\Delta x$.</caption><thead><tr><th>Δx</th><th>τ(Δx)</th></tr></thead><tbody><tr><td>0.01</td><td>-0.796237</td></tr><tr><td>0.005</td><td>-0.416073</td></tr><tr><td>0.0025</td><td>-0.222488</td></tr></tbody></table>

In Figure 1, we still consider $\Delta x \in \{0.01, 0.005, 0.0025\}$ and plot for each $\Delta x$ the time evolution of the modulated energy on $[0, T]$ (left part) and the density $\rho(T, \cdot)$ (right part).

c) A well-prepared data on a fine mesh. Here, we still consider the initial data of the form (135)-(136) with $s = 1$. The mesh size is $\Delta x = 10^{-3}$ and the time step is $\Delta t = \frac{1}{2}\Delta x$. The final time is $T = 20\Delta t$. We report in Table 3 the values of the modulated energy at initial and final time for different values of $\varepsilon$. We expect the scheme to converge in $L^\infty([0; T]; L^2(\mathbb{T}))$ but a priori not in $L^\infty([0; T]; H^1(\mathbb{T}))$ since there is a loss of one power of $\varepsilon$ when we differentiate (133). We see that the modulated energy

<table><caption>Table 3: Modulated energy at final and initial time for different values of $\varepsilon$ for an initial data of the form (135)-(136) with $s = 1$.</caption><thead><tr><th>ε</th><th>E(ρ<sup>N<sub>T</sub></sup>, u<sup>N<sub>T</sub></sup>, φ<sup>N<sub>T</sub></sup>|ū, φ̄)</th><th>E(ρ<sup>0</sup>, u<sup>0</sup>, φ<sup>0</sup>|ū, φ̄)</th></tr></thead><tbody><tr><td>0.1</td><td>0.0225142</td><td>0.022534</td></tr><tr><td>0.05</td><td>0.00562696</td><td>0.00563343</td></tr><tr><td>0.025</td><td>0.00140584</td><td>0.00140829</td></tr><tr><td>0.0125</td><td>0.000353224</td><td>0.00035338</td></tr></tbody></table>

still decays as $\varepsilon$ decays towards zero. The order of convergence is $O(\varepsilon^2)$. In Figures 2, we represent the initial density and the final density on the refined mesh for the values of $\varepsilon$ given in Table 3. We see that the spatial oscillations are still present and it seems that there is also an oscillatory behavior in time. To investigate the time oscillations, we plot in Figure 3 the time evolution of $\varepsilon^{-2}\|\rho(t, \cdot) - \bar{\rho}\|_{L^2(\mathbb{T})}$
<figure><img src="image_8.png" /><figcaption>Figure 2: Plots of the initial density $\rho(0, \cdot)$ and final density $\rho(20\Delta t, \cdot)$ for $\varepsilon \in \{0.1, 0.05, 0.025\}$ on the fine mesh: $\Delta x = 10^{-3}$ and $\Delta t = \frac{1}{2}\Delta x$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 3: Time evolution of $t \in [0, 5] \mapsto \|\rho(t, \cdot) - \bar{\rho}\|_{L^2}/\varepsilon^2$ for $\varepsilon \in \{0.01, 0.05, 0.1\}$. $\Delta x = 1/200$ and $\Delta t = \frac{1}{2}\Delta x$.</figcaption></figure>

(the rescaling by $\varepsilon^2$ is needed to get comparable amplitudes). The oscillation period does not seem to depend strongly on $\varepsilon$. This behavior has already been observed in [7] thanks to a linear stability analysis. The main reason is that the zero order term in the Poisson equation tends to stabilize the high spatial frequency mode for the electric potential.

### 3.3 An ill-prepared data on a fine mesh.

We now consider the following initial data

$$ \rho_{\varepsilon}^{\text{ini}}(x) - \bar{\rho} = \frac{1}{2} \sin(2\pi x \lfloor \varepsilon^{-1} \rfloor), \quad u_{\varepsilon}^{\text{ini}}(x) - \bar{u} = \sin(2\pi x), \quad x \in [0, 1], \qquad (138) $$

which is actually not well prepared. The mesh size is $\Delta x = 10^{-3}$, the time step is $\Delta t = \frac{1}{2}\Delta x$ and the final time is $T = 20\Delta t$. We report on Table 4 the values of the modulated energy at initial and final time for three values of $\varepsilon$. In figures 4, we represent the density at initial and final time for three different values of $\varepsilon$.
<table><caption>Table 4: Modulated energy at final and initial time for different values of ε for an ill prepared initial data of the form (138).</caption><thead><tr><th>ε</th><th>E(ρ<sup>N<sub>T</sub></sup>, u<sup>N<sub>T</sub></sup>, φ<sup>N<sub>T</sub></sup>|ū, φ̄)</th><th>E(ρ<sup>0</sup>, u<sup>0</sup>, φ<sup>0</sup>|ū, φ̄)</th></tr></thead><tbody><tr><td>0.1</td><td>2.24793</td><td>2.2534</td></tr><tr><td>0.05</td><td>2.24739</td><td>2.25337</td></tr><tr><td>0.025</td><td>2.2465</td><td>2.25326</td></tr></tbody></table>

<figure><img src="image_10.png" /><figcaption>Figure 4: Plots of the initial density ρ(0, x) and final density ρ(20Δt, x) for ε ∈ {0.1, 0.05, 0.025} on a fine mesh: Δx = 10⁻³ and Δt = 1/2Δx.</figcaption></figure>

## 3.4 Five-branch solution

This test case is inspired from [18, 7, 2] but here we consider (because of the periodic boundary conditions), the following initial condition

$$ \rho^{\text{ini}}(x) = \begin{cases} 0.1 + \exp\left(\frac{0.1}{(x-3\pi/4)(x-5\pi/4)}\right) & \text{if } \frac{3\pi}{4} < x < \frac{5\pi}{4}, \\ 0.1 & \text{if } x \in [0, 2\pi] \setminus (\frac{3\pi}{4}, \frac{5\pi}{4}). \end{cases} $$

$$ u^{\text{ini}}(x) = \sin^3(x), \quad x \in [0, 2\pi]. $$

Note that $\rho^{\text{ini}} \in C_c^{\infty}((0, 2\pi))$ so it has a smooth $2\pi$-periodic extension. The electric potential $\phi^{\text{ini}}$ is computed from the nonlinear Poisson equation. The space domain is $2\pi\mathbb{T} \equiv [0, 2\pi)$ discretized with a mesh size $\Delta x = 2\pi/400$ and the time step is $\Delta t = \Delta x/2$. We set the final time $T = 0.5$. We run the method presented before with the limit scheme ($\epsilon = 0$) which consists in replacing the Poisson equation by the algebraic relation $\phi^{n+1} = -\log(\rho^{n+1})$. Moreover, we also present some results obtained by a numerical method for the limit model (3) based on the same spatial and temporal discretization. The main difference lies in the fact the discretization we use for (3) is conservative which is not the case for the limit scheme ($\epsilon = 0$). The numerical parameters are the same for the three solvers.

In Figures 5 and 6, we plot the space dependency of the density and velocity at time $t = 0.5$ for $\epsilon \in \{10^{-4}, 10^{-2}, 0\}$. The case $\epsilon = 0$ corresponds to the limit scheme. One can observe the scheme for $\epsilon = 10^{-4}$ and $\epsilon = 0$ are almost indistinguishable whereas oscillations are present for $\epsilon = 10^{-2}$. In Figure 6, the density is plotted at time $t = 0.5$ for $\epsilon \in \{10^{-4}, 0\}$ as before but we also add the result obtained by the scheme for (3). In the region where the solution is smooth, the three curves are very similar whereas some differences can be observed around the discontinuity $x \approx 1.9$. Indeed, the propagation of the discontinuity seems to be different due to the non conservative treatment of the term $\rho \partial_x \phi$ (see the inset in Figure 6).

**Acknowledgement.** MB acknowledges the financial support from Inria Rennes. He is also grateful to Institut de Recherche Mathématiques de Rennes and Ecole Normale Supérieure de Rennes for the material support during the elaboration of this work. This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (grant agreement No 101052200 EUROfusion). Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Commission. The work has been supported by the French Federation for Magnetic Fusion Studies (FR-FCM).
<figure><img src="image_11.png" /><figcaption>Figure 5: Five-branch test: comparison of the scheme for $\epsilon = 10^{-2}, 10^{-4}$ and the asymptotic scheme $\epsilon = 0$. Left: density $\rho(t = 0.5, \cdot)$. Right: velocity $u(t = 0.5, \cdot)$. $\Delta x = 2\pi/400$ and $\Delta t = \frac{1}{2}\Delta x$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 6: Five-branch test: comparison of the scheme for $\epsilon = 10^{-4}$, the asymptotic scheme $\epsilon = 0$ and the limit model. Density $\rho(t = 0.5, x)$. The inset is a zoom around $x = 1.9$. $\Delta x = 2\pi/400$ and $\Delta t = \frac{1}{2}\Delta x$.</figcaption></figure>
References

[1] A. Alvarez Laguna, T. Pichard, T. Magin, P. Chabert, A. Bourdon, and M. Massot. An asymptotic preserving well-balanced scheme for the isothermal fluid equations in low-temperature plasmas at low-pressure. *Journal of Computational Physics*, 419:109634, 2020.

[2] K. Arun and R. Ghorai. An asymptotic preserving scheme for the Euler-Poisson-Boltzmann system in the quasineutral limit. *Computers & Mathematics with Applications*, 185:1–28, 2025.

[3] D. Bresch, M. Gisclon, and I. Lacroix-Violet. On Navier-Stokes-Korteweg and Euler-Korteweg systems: Application to quantum fluids models. *Archive for Rational Mechanics and Analysis*, 233:975–1025, 2019.

[4] C. Chainais-Hillairet, B. Merlet, and A. Zurek. Convergence of a finite volume scheme for a parabolic system with a free boundary modeling concrete carbonation. *ESAIM: M2AN*, (52):457–480, 2018.

[5] S. Cordier and E. Grenier. Quasineutral limit of an Euler-Poisson system arising from plasma physics. *Comm. Partial Differential Equations*, 25:1099–1113, 2000.

[6] P. Crispel, P. Degond, and M.-H. Vignal. An asymptotic preserving scheme for the two-fluid Euler-Poisson model in the quasineutral limit. *J. Comput. Phys.*, 223:208–234, 2007.

[7] P. Degond, H. Liu, D. Savelief, and M.-H. Vignal. Approximation of the Euler-Poisson-Boltzmann model in the quasineutral limit. *J. Sci. Comput.*, 51:59–86, 2012.

[8] P. Degond, J.-G. Liu, and M.-H. Vignal. Analysis of an asymptotic preserving scheme for the Euler-Poisson system in the quasineutral limit. *SIAM J. Numer. Anal.*, 46:1298–1322, 2008.

[9] R. Eymard, T. Gallouët, R. Herbin, and J.-C. Latché. Convergence of the MAC scheme for the compressible Stokes equations. *SIAM J. Numer. Anal.*, 48:2218–2246, 2010.

[10] S. Fabre. Stability analysis of the Euler-Poisson system. *J. Comput. Phys.*, 101:445–451, 1992.

[11] L. Gastaldo, R. Herbin, and J.-C. Latché. A discretization of the phase mass balance in fractional step algorithms for the drift-flux model. *IMA Journal of Numerical Analysis*, 31(1):116–146, 09 2009.

[12] D. Gérard-Varet, D. Han-Kwan, and F. Rousset. Quasineutral limit of the Euler-Poisson system for ions in a domain with boundaries. *Indiana Univ. Math. J.*, 2:359–402, 2013.

[13] D. Han-Kwan. Quasineutral limit of the Vlasov-Poisson system with massless electrons. *Communications in Partial Differential Equations*, 36(8):1385–1425, 2011.

[14] R. Herbin, W. Kheriji, and J.-C. Latché. On some implicit and semi-implicit staggered schemes for the Shallow-Water and Euler equations. *Mathematical Modelling and Numerical Analysis*, 48:1807–1857, 48.

[15] T. Kato. The cauchy problem for quasi-linear symmetric hyperbolic systems. *Archive for Rational Mechanics and Analysis*, 58:180–205, 1995.

[16] D. Lannes, F. Linares, and J.-C. Saut. The Cauchy Problem for the Euler-Poisson System and Derivation of the Zakharov-Kuznetsov Equation, pages 181–213. Springer New York, New York, NY, 2013.

[17] M. Lieberman and A. Lichtenberg. *Principles of plasma discharge and materials processing*. New-York,Wiley, 1994.

[18] H. Liu and Z. Wang. A field space-based level set method for computing multi-valued solutions to 1d euler-poisson equations. *J. Comput. Phys.*, (25):591–614, 2007.

[19] X. Pu. Quasineutral limit of the pressureless Euler-Poisson equation. *Appl. Math. Lett.*, 30:33–37, 2014.

[20] X. Pu and C. Guo. Quasineutral limit of the pressureless Euler-Poisson equation for ions. *Quarterly of Applied Mathematics*, 74(2):245–273, 2016.
[21] S. Benzoni-Gavage and D. Serre. *Multi-dimensional Hyperbolic Partial Differential Equations*. Oxford Mathematical Monographs, Oxford University Press, 2007.