# Convergence and long-time behavior of finite volumes for a generalized Poisson-Nernst-Planck system with cross-diffusion and size exclusion

Clément Cancès, Maxime Herda, Annamaria Massimini

► To cite this version:

Clément Cancès, Maxime Herda, Annamaria Massimini. Convergence and long-time behavior of finite volumes for a generalized Poisson-Nernst-Planck system with cross-diffusion and size exclusion. 2024. [hal-04777272v2](https://hal.science/hal-04777272v2)

HAL Id: hal-04777272

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

Preprint submitted on 29 Aug 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
# CONVERGENCE AND LONG-TIME BEHAVIOR OF FINITE VOLUMES FOR A GENERALIZED POISSON-NERNST-PLANCK SYSTEM WITH CROSS-DIFFUSION AND SIZE EXCLUSION

CLÉMENT CANCÈS, MAXIME HERDA, AND ANNAMARIA MASSIMINI

ABSTRACT. We present a finite volume scheme for modeling the diffusion of charged particles, specifically ions, in constrained geometries using a degenerate Poisson-Nernst-Planck system with size exclusion yielding cross-diffusion. Our method utilizes a two-point flux approximation and is part of the exponentially fitted scheme framework. The scheme is shown to be thermodynamically consistent, as it ensures the decay of some discrete version of the free energy. Classical numerical analysis results – existence of a discrete solution, convergence of the scheme as the grid size and the time step go to 0 – follow. We also investigate the long-time behavior of the scheme, both from a theoretical and numerical point of view. Numerical simulations confirm our findings, but also point out some possibly very slow convergence towards equilibrium of the system under consideration.

## 1. PRESENTATION OF THE PROBLEM

1.1. **The continuous generalized Poisson-Nernst-Planck model.** Motivated by the transfer of ions in confined geometries, Burger et al. introduced in [5] a model accounting for cross-diffusion and size-exclusion effects. The model [6] further incorporated self-consistent electric interaction. In this model, I species, the volume fractions of which being denoted by $U = (u_i)_{1\le i\le I}$, are subject to diffusion as well as to electric forces induced by a self-consistent electrostatic potential. Denote by $\Omega \subset \mathbb{R}^d$ a bounded connected polyhedral domain, then the conservation of the volume occupied by the species i is written

$$ (1) \qquad \partial_t u_i + \nabla \cdot F_i = 0, \quad i = 1, \dots, I, $$

with the flux of the species i being (formally) given by

$$ (2) \qquad F_i = -D_i (u_0 \nabla u_i - u_i \nabla u_0 + u_0 u_i z_i \nabla \phi) = -D_i u_i u_0 \nabla \left( \log \left( \frac{u_i}{u_0} \right) + z_i \phi \right). $$

In the above expression, $D_i > 0$ denotes the diffusion coefficient of the species i. The quantity

$$ (3) \qquad u_0 = 1 - \sum_{i=1}^{I} u_i $$

shall be thought of as the volume fraction of available space for the ions, possibly occupied by a motile and electro-neutral solvent. The quantity $u_0$ is then required to remain non-negative, leading to size exclusion for the other species $u_i$, $i = 1, \dots, I$. Denoting by $z_i \in \mathbb{Z}$ the charge of species i and by $\lambda > 0$ the (scaled) Debye length, then the electrostatic potential $\phi$ solves the Poisson equation

$$ (4) \qquad -\lambda^2 \Delta \phi = \sum_{i=1}^{I} z_i u_i + f $$

for some prescribed background charge density $f$. We consider boundary conditions of mixed type for the electric potential. More precisely, we assume that the boundary $\partial\Omega$ of the domain can

<footnote>Date: August 29, 2025.</footnote>

<footnote n="2010">2010 Mathematics Subject Classification. 65M08, 65M12, 35K51.</footnote>

<footnote>Key words and phrases. Drift-diffusion, cross-diffusion, exponential fitting, free energy decay.</footnote>
be split into an insulator part $\Gamma^N$, which is a relatively open $C^1$ part of the boundary, and its complement $\Gamma^D$ on which a Dirichlet boundary condition is imposed:

$$ (5) \qquad \nabla\phi \cdot n = 0 \quad \text{on } \Gamma^N \quad \text{and} \quad \phi = \phi^D \text{ on } \Gamma^D, $$

where $n$ is the unit normal to $\partial\Omega$. Throughout this paper, we will assume that $f \in L^\infty(\Omega)$ and that $\phi^D$ is the trace of an $L^\infty \cap H^1(\Omega)$ function (which we also denote by $\phi^D$). Neither $f$ nor $\phi^D$ depend on time. Boundary conditions of various types can be considered for the conservation laws (1)–(2), for instance a Robin type boundary condition modeling electrochemical reaction thanks to a Butler–Volmer type formula, see for instance [8], or boundary conditions of mixed Dirichlet–Neumann type as in [31]. In the presentation of the scheme, we assume for simplicity that the system is isolated, in the sense that

$$ (6) \qquad F_i \cdot n = 0 \quad \text{on } \partial\Omega, \quad i = 1, \dots, I. $$

The system is finally complemented with initial conditions $u_i(t=0) = u_i^0$ with

$$ (7) \qquad u_i^0 \ge 0 \quad \text{and} \quad \int_{\Omega} u_i^0 > 0 \quad \text{for} \quad i=0, \dots, I \quad \text{and} \quad \sum_{i=0}^{I} u_i^0 = 1. $$

We then denote by

$$ \mathcal{A} = \left\{ U = (u_i)_{1 \le i \le I} \in (\mathbb{R}_{\ge 0})^I \;\middle|\; \sum_{i=1}^{I} u_i \le 1 \right\} $$

the set in which the volume fractions have to take their values.

## 1.2. Entropy structure of the model.
Let us now describe the entropy (or formal gradient flow) structure of the model. Introduce the Slotboom variables $w_i = \frac{u_i}{u_0}e^{z_i\phi}$, then the fluxes (2) rewrite as

$$ (8) \qquad F_i = -D_i u_0^2 e^{-z_i \phi} \nabla w_i, \quad i = 1, \dots, I. $$

If the so-called electrochemical potentials defined by $\mu_i = \log w_i = \log \frac{u_i}{u_0} + z_i\phi$ are regular enough, say in $H^1 \cap L^\infty(\Omega)$, then $\nabla \log w_i \cdot \nabla w_i = 4|\nabla\sqrt{w_i}|^2$. Therefore, multiplying (1) by $\mu_i$, integrating over $\Omega$ and summing over $i = 1, \dots, I$ then yields

$$ (9) \qquad \frac{d}{dt} \mathcal{H} + 4 \int_{\Omega} \sum_{i=1}^{I} D_i u_0^2 e^{-z_i \phi} |\nabla \sqrt{w_i}|^2 = 0, $$

where, denoting the mixing (neg)entropy density function $H: \mathcal{A} \to [-\log(I+1), +\infty)$ by

$$ (10) \qquad H(U) = u_0 \log(u_0) + \sum_{i=1}^{I} u_i \log(u_i), $$

with $u_0$ seen as a function of $U$, the free energy $\mathcal{H}$ is given by

$$ (11) \qquad \mathcal{H} = \int_{\Omega} H(U) + \frac{\lambda^2}{2} \int_{\Omega} |\nabla\phi|^2 - \lambda^2 \int_{\Gamma^D} \phi^D \nabla\phi \cdot n. $$

Assume that the $u_i$ are positive for $i = 0, \dots, I$ (as proved in the discrete case later on), then the second term in (9) is well-defined and non-negative. As a consequence, the free energy decays along-time, as a manifestation of the second principle of thermodynamics. Observe that $\mathcal{H}$ need not be non-negative but may be bounded uniformly from below by a constant depending only on $\lambda, f$, and $\phi^D$.
1.3. **Weak solution.** As we aim to prove rigorously the convergence of the scheme to be presented in Section 2, we need a proper notion of solution for the continuous problem (1)–(7). The volume fractions we seek are such that

$$ (12) \qquad 0 \le u_i \le 1 \text{ for } i = 0, \dots, I \text{ and a.e. } (t, x) \in \mathbb{R}_{\ge 0} \times \Omega. $$

The electric potential $\phi$ then solves the Poisson equation (4) with a bounded right-hand side and boundary condition, so it satisfies

$$ (13) \qquad \|\phi\|_{L^{\infty}(\mathbb{R}_{\ge 0} \times \Omega)} + \|\phi - \phi^D\|_{L^{\infty}(\mathbb{R}_{\ge 0}; V)} \le C(\lambda, (z_i)_i, \Omega, f), $$

where

$$ V = \{v \in H^1(\Omega) \mid v|_{\mathbb{R}^D} = 0\} $$

is equipped with the $H^1(\Omega)$ semi-norm $\|\nabla v\|_{L^2}$, which is a norm thanks to the Poincaré inequality. We have emphasized in the right-hand side of (13) the dependence of the bound on the data, especially on the Debye length $\lambda$. The estimate (9) is the other main estimate on which our study will rely. It provides enough control to define a proper notion of weak solutions, but at the price of a suitable reformulation of the fluxes (2). Let us first remark that, since the free energy $\mathcal{H}$ is bounded (see Lemma 3.7), integrating (9) over $t \in \mathbb{R}_{\ge 0}$ provides

$$ (14) \qquad \iint_{\mathbb{R}_{\ge 0} \times \Omega} \sum_{i=1}^{I} D_i u_0^2 e^{-z_i \phi} |\nabla \sqrt{w_i}|^2 \le C. $$

As already noticed in [31], the above inequality yields a control in $L_{\text{loc}}^2(\mathbb{R}_{\ge 0}; H^1(\Omega))$ on $\sqrt{u_0}$ and on $u_0$ itself, as well as control on some product terms. More precisely, in the numerical analysis exposed hereafter, we derive a $L_{\text{loc}}^2(\mathbb{R}_{\ge 0}; H^1(\Omega))$ estimate on $u_i\sqrt{u_0}$, which, together with (12) provides a $L_{\text{loc}}^2(\mathbb{R}_{\ge 0}; H^1(\Omega))$ control on $u_i u_0$ too. Hence, all the terms in the following expression of the fluxes

$$ (15) \qquad F_i = -D_i (\nabla(u_0 u_i) - 4u_i \sqrt{u_0} \nabla \sqrt{u_0} + u_i u_0 z_i \nabla \phi), \quad 1 \le i \le I, $$

have a clear mathematical sense, motivating the following notion of weak solution.

**Definition 1.1** (Weak solution). $(U, \phi)$ is said to be a weak solution to (1)–(7) if

* $u_i \in L^\infty(\mathbb{R}_{\ge 0} \times \Omega)$ with $U \in \mathcal{A}$ a.e. in $\mathbb{R}_{\ge 0} \times \Omega$, with moreover $\sqrt{u_0}$ and $u_i\sqrt{u_0}$ belonging to $L_{\text{loc}}^2(\mathbb{R}_{\ge 0}; H^1(\Omega))$ for $u_0$ defined as in (3);

* $\phi \in L^\infty(\mathbb{R}_{\ge 0} \times \Omega)$ with $\phi - \phi^D \in L^\infty(\mathbb{R}_{\ge 0}; V)$ satisfies

$$ (16) \qquad \int_{\Omega} \lambda^2 \nabla \phi(t, x) \cdot \nabla \psi(x) dx = \int_{\Omega} \left( \sum_{i=1}^{I} z_i u_i(t, x) + f(x) \right) \psi(x) dx $$

for all $\psi \in C_c^1(\Omega \cup \Gamma^N)$ and a.e. $t \ge 0$;

* for all $\varphi \in C_c^1(\mathbb{R}_{\ge 0} \times \bar{\Omega})$ and all $i = 1, \dots, I$, there holds

$$ (17) \qquad \iint_{\mathbb{R}_{\ge 0} \times \Omega} u_i \partial_t \varphi + \int_{\Omega} u_i^0 \varphi(0, \cdot) - \iint_{\mathbb{R}_{\ge 0} \times \Omega} D_i (\nabla(u_0 u_i) - 4u_i \sqrt{u_0} \nabla \sqrt{u_0} + u_i u_0 z_i \nabla \phi) \cdot \nabla \varphi = 0. $$

1.4. **Equilibrium and modified Poisson–Boltzmann equation.** As the model has a gradient flow structure and is equipped with no-flux boundary conditions, its solution converges towards a steady state $(u^\infty, \phi^\infty)$ which correspond to vanishing fluxes $F_i^\infty$ for all $i = 1, \dots, I$. It is suggested in [6] that the long-time limit profile $(u^\infty, \phi^\infty)$ should correspond to constant in space electrochemical potentials $\mu^\infty = (\mu_i^\infty)_{1 \le i \le I}$. Define the smooth functions $v_i : \mathbb{R} \times \mathbb{R}^I \to (0, 1)$ by

$$ (18) \qquad v_i(y, \boldsymbol{\xi}) = \frac{e^{\xi_i - z_i y}}{1 + \sum_{j=1}^{I} e^{\xi_j - z_j y}} \quad \text{for } 1 \le i \le I, \quad \text{and} \quad v_0(y, \boldsymbol{\xi}) = 1 - \sum_{j=1}^{I} v_j(y, \boldsymbol{\xi}), $$

where $\boldsymbol{\xi} = (\xi_j)_{1 \le j \le I} \in \mathbb{R}^I$. Then the volume fractions $U^\infty = (u_i^\infty)_{1 \le i \le I}$ satisfy

$$ (19) \qquad u_i^\infty(x) = v_i(\phi^\infty(x), \boldsymbol{\mu}^\infty), \quad i = 0, \dots, I, x \in \Omega. $$
Plugging these expressions in the Poisson equation (4) provides what we call the *modified Poisson-Boltzmann equation*, i.e.

$$ (20) \qquad -\lambda^2 \Delta \phi^\infty + r(\phi^\infty, \boldsymbol{\mu}^\infty) = f \quad \text{in } \Omega, $$

complemented with the same boundary conditions (5). In (20), $r$ is the smooth function from $\mathbb{R} \times \mathbb{R}^I$ to $\mathbb{R}$ defined by

$$ (21) \qquad r(y, \boldsymbol{\xi}) = - \frac{\sum_{i=1}^{I} z_i e^{\xi_i - z_i y}}{1 + \sum_{i=1}^{I} e^{\xi_i - z_i y}} = - \sum_{i=1}^{I} z_i v_i(y, \boldsymbol{\xi}). $$

To close the system, we still need to determine the real constants $\boldsymbol{\mu}^\infty$. Their values are calculated to guarantee the conservation-of-total-mass constraints

$$ (22) \qquad \int_{\Omega} u_i^\infty = \int_{\Omega} u_i^0. $$

By setting $z_0 = 0$, then one readily checks that since the ($v_i$) by convexity

$$ \partial_y r(y, \boldsymbol{\xi}) = \sum_{i=0}^{I} z_i^2 v_i(y, \boldsymbol{\xi}) - \left( \sum_{i=0}^{I} z_i v_i(y, \boldsymbol{\xi}) \right)^2 \geq 0, $$

the last inequality being a consequence of the convexity of $y \mapsto y^2$ as $\sum_{i=0}^{I} v_i(y, \boldsymbol{\xi}) = 1$. Therefore, $r$ is non-decreasing w.r.t. its first variable $y$. The monotonicity and boundedness of $r$ allow us to conclude, via classical monotonicity arguments, the well-posedness of the problem (20)&(5). Going one step further allows us to check that the solution $(\phi^\infty, \boldsymbol{\mu}^\infty)$ is a minimizer of the functional $\Psi : (\phi^D + V) \times \mathbb{R}^I \to \mathbb{R}$ defined by

$$ (23) \qquad \Psi(y, \boldsymbol{\xi}) = \int_{\Omega} \left[ \frac{\lambda^2}{2} |\nabla y|^2 + \log \left( 1 + \sum_{i=1}^{I} e^{\xi_i - z_i y} \right) - fy - \sum_{i=1}^{I} \xi_i u_i^0 \right]. $$

One readily checks that $\Psi$ is strictly convex, and thus the minimiser $(\phi^\infty, \boldsymbol{\mu}^\infty)$ is unique. Moreover, one easily verifies that minimizers for $\Psi$ are steady states of the system (1)-(7). Despite hints in these directions, the uniqueness of the steady state is neither established in [6], nor, up to our knowledge, in more recent contributions. Furthermore, no quantitative convergence result has been established so far. Non-quantitative results were however recently presented in [19] in the absence of electrical interactions, i.e. $z_i = 0$ for all $i$. This gap is partially filled in the discrete case. We rigorously establish in Section 5 the uniqueness of the steady states, which are minimisers of some discrete counterpart of the functional $\Psi$, cf (35), and the (non-quantitative) convergence towards such a steady profile as time goes to $+\infty$.

**1.5. Goal and positioning of the paper.** The goal of this paper is to provide a developed mathematical study of a numerical scheme that was proposed in our previous contribution [11]. It enters the broad family of the two-point flux approximation (TPFA) finite volume schemes [26], which has been applied to several cross-diffusion problems in recent years (see e.g. [10, 15, 36, 30, 13, 37, 9]). Our scheme mixes ideas from the so-called square-root approximation (SQRA) finite volume scheme [12] (see also [38, 33]) with some features of the Scharfetter-Gummel scheme [39, 18], especially the use of the Bernoulli function (28) in the definition of the fluxes. Note that a purely SQRA approach was also proposed in [11], but since it is overpassed by the mixed approach both in terms of accuracy and of robustness in the small Debye length (or quasi-neutral) limit, we do not develop the analysis here. Let us however stress that our whole study extends to this full-SQRA scheme at the price of very minor modifications sketched in [11].

As highlighted in [11], our scheme fulfils key properties. First, it is second order accurate in space, in opposition to the scheme proposed in [7] in which upwinding leads to a mere first order consistency. Even though our scheme does not use the so-called entropy variables (the electrochemical potentials in our framework) as unknowns in opposition to the finite element scheme presented in [32], our scheme satisfies a discrete energy dissipation inequality, which should be thought of as a discrete counterpart to (9), and which is key for its numerical analysis. Keeping
conservative variables as primary variables yields well-behaved non-linear solvers (we used Newton-Raphson), allowing us to consider larger time steps than for the finite element methods of [32] or for the finite volume method of [2]. Together with the stability to be used for proving the convergence of the scheme, the discrete energy dissipation inequality implies that the scheme preserves exactly the steady states of the continuous model in the sense that, at the discrete level, fluxes vanish when concentrations are of the form (19). The unique steady state of our scheme is then characterized only by the electric potential which solves a discrete counterpart of the Poisson-Boltzmann equation (20), which can be computed directly in an efficient way.

In the present contribution, we first rigorously establish the convergence of the scheme when the discretization parameters (mesh size and time step) tend to 0. As in the continuous setting [31], the proof is based on compactness arguments. It strongly relies on the discrete counterpart of the energy dissipation estimate (9), following the boundedness-by-entropy method [35], and more precisely a discrete version of it [16, 37]. Uniqueness of the weak solutions is moreover established in [31] under moderate regularity assumptions thanks to the uniqueness technique of Gajewski [27, 28].

We also investigate the long-time behaviour of the numerical scheme. We show in particular that the discrete solution to our scheme for the evolution problem converges as time goes to infinity towards the unique solution to the discrete steady problem. As the proof relies on LaSalle's invariance principle, no quantitative convergence rate is derived, similarly to what is proposed in [19] in the continuous setting but in absence of electrical interactions. By establishing this result, we rigorously establish in the discrete setting the long-time behavior conjectured in [6] and recalled in Section 1.4. Numerical simulations suggest that the convergence rate strongly depends on the Debye length, but also on the initial profile $U^0$ (in opposition to the classical Poisson-Nernst-Planck problem without volume filling [3]).

## 2. THE FINITE VOLUME SCHEME AND MAIN RESULTS

First, we introduce the time discretization and the spatial mesh of the domain $\Omega$.

### 2.1. Space and time discretizations.
The mesh will be assumed to be admissible in the sense of [26], namely it fulfils the so-called orthogonality condition, which is usual for two-point flux approximation finite volumes.

Let $\mathcal{T}$ denote a family of non-empty, disjoint, convex, open, and polygonal control volumes $K \in \mathcal{T}$, whose Lebesgue measure is denoted by $m_K$. We also assume that control volumes partition the domain, meaning that $\bar{\Omega} = \bigcup_{K \in \mathcal{T}} \bar{K}$. Further, we call $\mathcal{E}$ a family of edges/faces, where each $\sigma \in \mathcal{E}$ is a closed subset of $\bar{\Omega}$ contained in a hyperplane of $\mathbb{R}^d$. Each $\sigma$ has a strictly positive $(d-1)$-dimensional Hausdorff (or Lebesgue) measure, denoted by $m_\sigma$. We use the abbreviation $K|L = \partial K \cap \partial L$ for the intersection between two distinct control volumes which is either empty or reduces to a face contained in $\mathcal{E}$. The subset of all interior faces is denoted by

$$ \mathcal{E}_{\text{int}} = \{\sigma \in \mathcal{E} \mid \sigma = K|L \text{ for some } K, L \in \mathcal{T}\}. $$

For any $K \in \mathcal{T}$, we assume that there exists a subset $\mathcal{E}_K$ of distinct elements of $\mathcal{E}$ such that the boundary of a control volume can be described by $\partial K = \bigcup_{\sigma \in \mathcal{E}_K} \sigma$ and, consequently, $\mathcal{E} = \bigcup_{K \in \mathcal{T}} \mathcal{E}_K$. Additionally, we assume that boundary edges $\mathcal{E}_{\text{ext}} = \mathcal{E} \setminus \mathcal{E}_{\text{int}}$ are either subsets of $\Gamma^D$ or $\Gamma^N$. To each control volume $K \in \mathcal{T}$ we assign a cell center $x_K \in K$ which satisfies the orthogonality condition: If $K, L$ share a face $\sigma = K|L$, then the vector $x_L - x_K$ is orthogonal to $\sigma = K|L$. The triplet $(\mathcal{T}, \mathcal{E}, \{x_K\}_{K \in \mathcal{T}})$ is called an admissible mesh.

A classical way to construct such admissible mesh is to place some seeds $(x_K)_{K \in \mathcal{T}}$ in the domain $\Omega$ and to define the corresponding Voronoï tessellation. As in practice the condition $x_K \in K$ can be relaxed, but still with the constraint $x_L - x_K$ points outward with respect to $K$, another popular choice to build admissible meshes from a simplicial mesh satisfying the Delaunay condition is to choose $x_K$ as the center of the circumcircle of the simplex $K$.

We introduce the notation $d_\sigma$ for the Euclidean distance between $x_K$ and $x_L$ if $\sigma = K|L$ or between $x_K$ and the affine hyperplane spanned by $\sigma$ if $\sigma \subset \partial\Omega$. We also denote by $d_{K\sigma} =$
$\text{dist}(x_K, \sigma)$, so that $d_\sigma = d_{K\sigma} + d_{L\sigma}$ if $\sigma = K|L \in \mathcal{E}_{\text{int}}$ and $d_\sigma = d_{K\sigma}$ if $\sigma \in \mathcal{E}_K \cap \mathcal{E}_{\text{ext}}$. The *transmittivity* of the edge $\sigma \in \mathcal{E}$ is defined by $a_\sigma = \frac{m_\sigma}{d_\sigma}$. The size of the mesh is

$$h_{\mathcal{T}} = \max_{K \in \mathcal{T}} \text{diam}(K)$$

where $\text{diam}(K)$ denotes the diameter of the cell $K$. The regularity of the mesh is defined by

$$\zeta_{\mathcal{T}} = \max_{K \in \mathcal{T}} \left( \text{card} \ \mathcal{E}_K ; \max_{\sigma \in \mathcal{E}_K} \frac{\text{diam}(K)}{d_{K\sigma}} \right).$$

For the time discretization, we decompose the time interval $\mathbb{R}_{\ge 0} := [0, +\infty)$ into an unbounded increasing sequence $(t^n)_{n \ge 0}$ with $t^0 = 0$ and possibly non-uniform time steps

$$\tau^n = t^n - t^{n-1} > 0, \quad n > 0.$$

We finally introduce $\Delta t = \sup_{n \in \mathbb{N} \setminus \{0\}} \tau^n$, which we assume to be finite. As we consider the unbounded time interval $\mathbb{R}_{\ge 0}$, we furthermore have to assume that $\inf_{n \ge 1} \tau^n > 0$.

## 2.2. The finite volume schemes.
We are now in position to define the finite volume scheme. Let us start with the discretization of the Poisson equation (4)–(5), which relies on a classical two-point flux approximation

$$(24) \qquad \lambda^2 \sum_{\sigma \in \mathcal{E}_K} a_\sigma (\phi_K^n - \phi_{K\sigma}^n) = m_K \left( f_K + \sum_{i=1}^I z_i u_{i,K}^n \right), \quad K \in \mathcal{T},$$

where $f_K$ is (possibly an approximation of) the mean value of $f$ on the cell $K$, and where

$$\phi_{K\sigma}^n = \begin{cases} \phi_L^n & \text{if } \sigma = K|L \in \mathcal{E}_{\text{int}}, \\ \phi_K^n & \text{if } \sigma \subset \Gamma^N, \\ \phi_\sigma^D = \frac{1}{m_\sigma} \int_\sigma \phi^D & \text{if } \sigma \subset \Gamma^D. \end{cases}$$

The equation (1) is discretized using a backward Euler method in time and finite volumes in space, leading to

$$(25) \qquad \frac{u_{i,K}^n - u_{i,K}^{n-1}}{\tau^n} m_K + \sum_{\sigma \in \mathcal{E}_K} F_{i,K\sigma}^n = 0, \quad i = 1, \dots, I, \quad K \in \mathcal{T}.$$

In accordance with (6), we set $F_{i,K\sigma}^n = 0$ if $\sigma \subset \partial\Omega$. For $\sigma = K|L$ an internal edge, then we define

$$(26) \qquad F_{i,K\sigma}^n = a_\sigma D_i (u_{i,K}^n u_{0,L}^n \mathfrak{B}(z_i(\phi_L^n - \phi_K^n)) - u_{i,L}^n u_{0,K}^n \mathfrak{B}(z_i(\phi_K^n - \phi_L^n))),$$

with

$$(27) \qquad u_{0,K}^n = 1 - \sum_{i=1}^I u_{i,K}^n, \quad K \in \mathcal{T}.$$

Formula (26) involves a function $\mathfrak{B} \in C^1(\mathbb{R}; \mathbb{R})$ which is (strictly) positive and satisfies $\mathfrak{B}(0) = 1$ and $\mathfrak{B}'(0) = -1/2$. The SQRA scheme studied in [11] corresponds to the choice $\mathfrak{B}(y) = e^{-y/2}$. Even if our analysis extends to this case without any particular difficulty, for the sake of readability we will focus here on the scheme referred to as the SG scheme in [11], for which $\mathfrak{B}$ is the Bernoulli function

$$(28) \qquad \mathfrak{B}(y) = \frac{y}{e^y - 1},$$

and which shows a better behavior (accuracy and robustness) with respect to the SQRA scheme. We refer to [11] for a comparison of the two approaches.

We would like to emphasize that the construction of the scheme which we will refer to as the SG scheme is not based on the original idea of [39]. Rather, we take advantage of the free-energy diminishing character of the SG scheme highlighted in [18] and exploited in [4, 40].

In order to close the system, it remains to define the discrete counterpart of $u^0$ as follows:

$$(29) \qquad u_{i,K}^0 = \frac{1}{m_K} \int_K u_i^0, \quad K \in \mathcal{T}, \quad i = 0, \dots, I.$$
Then we infer from (7) that

$$ (30) \qquad \sum_{i=0}^{I} u_{i,K}^{0} = 1 \quad \text{for all } K \in \mathcal{T}, \quad \text{and} \quad \sum_{K \in \mathcal{T}} u_{i,K}^{0} m_{K} = \int_{\Omega} u_{i}^{0} > 0 \quad \text{for } i = 0, \dots, I. $$

In what follows, we denote by $U_K^n = (u_{i,K}^n)_{i=1,\dots,I}$ for $K \in \mathcal{T}$ and $n \ge 0$.

The consistency of the discrete fluxes (26) with the continuous ones (2) might not look completely obvious. It follows from the identity

$$ (31) \qquad F_{i,K\sigma}^n = a_\sigma D_i u_{0,K}^n u_{0,L}^n \mathfrak{M}(e^{-z_i \phi_K^n}, e^{-z_i \phi_L^n}) (w_{i,K}^n - w_{i,L}^n) $$

where $w_{i,K}^n = \frac{u_{i,K}^n}{u_{0,K}^n} e^{z_i \phi_K^n}$ are the discrete counterpart of the Slotboom variables and, by assumption, $U_K^n \in (0,1)^I$. Later, in Proposition 3.2, this bound will be rigorously established. The mean function $\mathfrak{M}$ appearing in (31) is defined by

$$ \mathfrak{M}(a,b) = \frac{\log(1/a) - \log(1/b)}{1/a - 1/b} \quad \text{for } a,b > 0 \text{ with } a \neq b, \text{ and } \mathfrak{M}(a,a) = a \text{ for } a > 0. $$

We refer to [34] for an extensive discussion on the influence of the choice of the Stolarsky mean $\mathfrak{M}$ on the scheme behavior. The consistency of the formula (31) with the expression (8) of the continuous flux is now clear, and assuming the existence of a regular solution to the continuous problem, inserting it in the scheme and performing Taylor expansions as in [12] shows the second order accuracy in space of our scheme on uniform Cartesian grids.

**2.3. Main results and organization of the paper.** We provide here simple presentations of our main results, which will be detailed and proven in the next sections. As already mentioned, the goals of this contribution are twofold. It first aims to complement the elements of analysis sketched in the contribution [11], by showing existence of discrete solutions and convergence of the scheme. Secondly it intends to investigate the long-time behavior of the scheme.

Our first main result concerns the well-posedness of the scheme.

**Theorem 2.1.** Given an admissible mesh $(\mathcal{T}, \mathcal{E}, (x_K)_{K \in \mathcal{T}})$ of $\Omega$ and a sequence of time steps $(\tau^n)_{n \ge 0}$ as in Section 2.1, then there exists (at least) one solution to the scheme (24)-(29) which satisfies

$$ 0 < u_{i,K}^n < 1 \quad \forall i = 0, \dots, I, K \in \mathcal{T}, n \ge 1. $$

Moreover, the discrete free energy $\mathcal{H}_{\mathcal{T}}^n$ defined later on in (41) is decaying along the time iterations

$$ (32) \qquad \mathcal{H}_{\mathcal{T}}^n + \tau^n \mathcal{D}_{\mathcal{T}}^n \le \mathcal{H}_{\mathcal{T}}^{n-1}, \quad n \ge 1, $$

for some dissipation rate $\mathcal{D}_{\mathcal{T}}^n \ge 0$ vanishing if and only if $((U_K^n)_{K \in \mathcal{T}}, (\phi_K^n)_{K \in \mathcal{T}})$ is the stationary solution to the scheme.

Theorem 2.1 states that, at fixed mesh and time step, the numerical scheme admits a positive solution and preserves the $L^\infty$-bounds. Moreover, the inequality (32) is a discrete counterpart to (9). The discrete version $\mathcal{H}_{\mathcal{T}}^n$ of the free energy functional (11) is rigorously defined in formula (41). Schemes that verify (32) in addition to more classical properties such as mass conservation or positivity preservation are often referred to in the literature as thermodynamically consistent or structure-preserving. Our scheme, therefore, falls into this class of schemes.

Theorem 2.1 allows to define the so-called approximate solution to the problem, which consists of functions $U_{\mathcal{T},\tau} = (u_{i,\mathcal{T},\tau})_{1 \le i \le I}$ and $\phi_{\mathcal{T},\tau}$, piecewise constant in both space and time, defined by

$$ (33) \qquad u_{i,\mathcal{T},\tau}(t,x) = u_{i,K}^n \quad \text{and} \quad \phi_{\mathcal{T},\tau}(t,x) = \phi_K^n \quad \text{if } (t,x) \in K \times (t^{n-1}, t^n]. $$

Here again, we reconstruct $u_{0,\mathcal{T},\tau}$ from $U_{\mathcal{T},\tau}$ by setting

$$ (34) \qquad u_{0,\mathcal{T},\tau} = 1 - \sum_{i=1}^{I} u_{i,\mathcal{T},\tau}. $$
Now, let $(\mathcal{T}_\ell, \mathcal{E}_\ell, \{x_K\}_{K \in \mathcal{T}_\ell})_{\ell \ge 1}$ and $(\tau_\ell)_{\ell \ge 1} = ((\tau_\ell^n)_{n \ge 1})_{\ell \ge 1}$ be respectively a sequence of admissible meshes and a sequence of time steps, in the sense of Section 2.1, such that

$$ \lim_{\ell \to \infty} h_{\mathcal{T}_\ell} = \lim_{\ell \to \infty} \Delta t_\ell = 0, \quad \Delta t_\ell = \max_n \tau_\ell^n, $$

while the mesh regularity factor $\zeta_{\mathcal{T}_\ell}$ remains uniformly bounded w.r.t. $\ell$, i.e. $\zeta_{\mathcal{T}_\ell} \le \zeta_\star < +\infty$. The convergence of our scheme can be stated in a very simplified way as follows:

**Theorem 2.2.** There exists a weak solution $(U, \phi)$ in the sense of Definition 1.1 such that, up to the extraction of a subsequence,

$$ \phi_{\mathcal{T}_\ell, \tau_\ell} \underset{\ell \to +\infty}{\longrightarrow} \phi \quad \text{and} \quad u_{0, \mathcal{T}_\ell, \tau_\ell} \underset{\ell \to +\infty}{\longrightarrow} u_0 \quad \text{in } L_{\text{loc}}^p(\mathbb{R}_{\ge 0}; L^p(\Omega)) \quad \forall p \in [1, +\infty), $$

and

$$ U_{\mathcal{T}_\ell, \tau_\ell} \underset{\ell \to +\infty}{\longrightarrow} U \quad \text{in the } L^\infty(\mathbb{R}_{\ge 0} \times \Omega)^I \text{ weak-\star sense.} $$

More convergence properties have to be established to prove Theorem 2.2. They are detailed in Section 4, complementing the preliminary presentation of [11].

Our last main theoretical result is about the long-time asymptotic of the scheme for a fixed mesh and sequence of time steps. We denote by $\Psi_\mathcal{T}: \mathbb{R}^T \times \mathbb{R}^I \to \mathbb{R}$ the functional defined by

$$ (35) \quad \Psi_{\mathcal{T}}(y_{\mathcal{T}}, \xi) = \frac{\lambda^2}{2} \sum_{\sigma \in \mathcal{E}} a_{\sigma}(y_K - y_{K\sigma})^2 + \sum_{K \in \mathcal{T}} m_K \log \left( 1 + \sum_{i=1}^I e^{\xi_i - z_i y_K} \right) \\ \qquad \qquad \qquad \qquad \qquad - \sum_{K \in \mathcal{T}} m_K \left[ f_K y_K + \sum_{i=1}^I \xi_i u_{i,K}^0 \right]. $$

In (35), we made a slight abuse of notation as we denote by $y_\mathcal{T} = (y_K)_{K \in \mathcal{T}}$ the element of $\mathbb{R}^T$, and the mirror values $y_{K\sigma}$ are defined by

$$ y_{K\sigma} = \begin{cases} y_L & \text{if } \sigma = K | L \in \mathcal{E}_{\text{int}}, \\ y_K & \text{if } \sigma \subset \Gamma^N, \\ \phi_\sigma^D & \text{if } \sigma \subset \Gamma^D. \end{cases} $$

**Theorem 2.3.** The functional $\Psi_\mathcal{T}$ is strictly convex and admits a unique minimizer denoted by $((\phi_K^\infty)_{K \in \mathcal{T}}, (\mu_{i,\mathcal{T}}^\infty)_{1 \le i \le I})$. Moreover, one has

$$ \phi_K^n \underset{n \to +\infty}{\longrightarrow} \phi_K^\infty \quad \forall K \in \mathcal{T}, $$

and

$$ u_{i,K}^n \underset{n \to +\infty}{\longrightarrow} v_i(\phi_K^\infty, \mu_\mathcal{T}^\infty) \quad \forall K \in \mathcal{T}, \ 0 \le i \le I. $$

We deduce from this theorem an efficient way to compute the long-time behavior of the scheme, which consists in seeking for $((\phi_K^\infty)_{K \in \mathcal{T}}, (\mu_{i,\mathcal{T}}^\infty)_{1 \le i \le I})$ belonging to argmin $\Psi_\mathcal{T}$ and to a posteriori reconstruct the volume fraction profiles thanks to the functions $v_i$ defined by (18). We also remark that the uniqueness of the long-time limit is established in the discrete setting, while the question remains open in the continuous setting. Conjecturing the uniqueness of such an asymptotic limit for the continuous model, which then has to be the unique minimizer of $\Psi$ defined by (23), then Theorem 2.3 illustrates in yet another aspect the structure-preserving nature of our scheme.

The rest of the paper is organized as follows. Section 3 is devoted to the proof of Theorem 2.1, as well as to some post-treatment of the energy / energy dissipation inequality (32) in order to derive some uniform bounds on quantities to be used in the convergence analysis carried out in Section 4, where the proof of Theorem 2.2 is detailed. Our last main theoretical result Theorem 2.3 is proved in Section 5. The numerical results presented in Section 6 contain some numerical evidence of the convergence of the scheme. They also illustrate that the situation is more complex than what Theorem 2.3 suggests at first glance. Although the approximate volume fractions are (strictly) positive as pointed out in Theorem 2.1, $u_{0,\mathcal{T},\tau}$ can take very small values in some parts
of the domain. The discrete fluxes (26) then take very small values as well, making the evolution extremely slow. As the linear and non-linear systems corresponding to the scheme are merely solved approximately, one furthermore has to pay a particular attention to preserve the positivity at the discrete level.

### 3. UNIFORM A PRIORI BOUNDS AND EXISTENCE OF A DISCRETE SOLUTION

The goal of this section is to prove Theorem 2.1, i.e. to show that the non-linear system corresponding to the scheme (24)-(27) admits at least one solution, and that beyond local conservativity, this solution preserves at the discrete level some key features of the model, namely the positivity of the volume fractions, as well as uniform $L^\infty$ bounds for them, and the decay of the free energy. Furthermore, in this part, precise quantification of the dissipated total free energy supplies enough uniform estimates, in the mesh size $h_T$ and time steps $(\tau_n)_{n \ge 1}$, to perform the convergence analysis in Section 4. The grid $\mathcal{T}$ and the time steps $(\tau^n)_{n \ge 1}$ remain fixed throughout this section.

3.1. Uniform a priori bounds and existence of a discrete solution. The first Lemma regards the discrete counterpart of the conservation of the total mass of the system.

**Lemma 3.1 (Mass conservation).** It holds that

$$ (36) \quad \sum_{K \in \mathcal{T}} u_{i,K}^n m_K = \sum_{K \in \mathcal{T}} u_{i,K}^{n-1} m_K = \sum_{K \in \mathcal{T}} u_{i,K}^0 m_K = \int_{\Omega} u_i^0 > 0, \quad 0 \le i \le I. $$

*Proof.* Since our scheme is locally conservative, i.e., $F_{i,K\sigma}^n + F_{i,L\sigma}^n = 0$ for all $\sigma = K|L \in \mathcal{E}_{\text{int}}$, $i = 1, \dots, I$, then summing (25) over $K$ shows by induction and thanks to (29) the result. $\square$

Since we are interested in discrete solutions with positive volume fractions $u_{i,K}^n$, we perform an eventually harmless modification of the flux formula (26) into

$$ (37) \quad F_{i,K\sigma}^n = a_\sigma D_i \left( (u_{i,K}^n)^+ (u_{0,L}^n)^+ \mathfrak{B} (z_i(\phi_L^n - \phi_K^n)) - (u_{i,L}^n)^+ (u_{0,K}^n)^+ \mathfrak{B} (z_i(\phi_K^n - \phi_L^n)) \right). $$

**Proposition 3.2 (Positivity of volume fractions).** Let $n \ge 1$, and let $(U_K^{n-1})_{K \in \mathcal{T}}$ be such that

$$ (38) \quad u_{i,K}^{n-1} \ge 0, \quad \sum_{i=0}^{I} u_{i,K}^{n-1} = 1 \quad \forall K \in \mathcal{T}, \quad \text{and} \quad \sum_{K \in \mathcal{T}} u_{i,K}^{n-1} m_K > 0. $$

Then any solution $(U_K^n, \phi_K^n)_{K \in \mathcal{T}, n \ge 1}$ to the modified scheme with (37) instead of (26) satisfies $u_{i,K}^n > 0$ for all $i = 0, \dots, I$ and all $K \in \mathcal{T}$.

*Proof.* Let us start by establishing the positivity of $u_{0,K}^n$. Assume for contradiction that there exists a cell $K \in \mathcal{T}$ such that $u_{0,K}^n \le 0$. Then we deduce from formula (37) that $F_{i,K\sigma}^n \ge 0$ for all $\sigma \in \mathcal{E}_K$ and all $i = 1, \dots, I$. Because of (27) and (38), this implies that

$$ 0 \ge u_{0,K}^n = u_{0,K}^{n-1} + \frac{\tau^n}{m_K} \sum_{i=1}^{I} \sum_{\sigma \in \mathcal{E}_K} F_{i,K\sigma}^n \ge 0. $$

In particular, all the fluxes $F_{i,K\sigma}^n$, $i = 1, \dots, I$ and $\sigma \in \mathcal{E}_K$, are equal to 0. In view of formula (37) and of the strict positivity of $\mathfrak{B}$, this implies either that $u_{i,K}^n \le 0$ for all $i$, which yields a contradiction with (27), or that $u_{0,L}^n \le 0$ for all the cells $L$ sharing an edge $\sigma = K|L$ with $K$. Since $\Omega$ is connected, one would obtain that $u_{0,K}^n = 0$ for all $K \in \mathcal{T}$ and thus that $\sum_{K \in \mathcal{T}} u_{0,K}^n m_K = 0$. This contradicts (36), and thus we necessarily have that $u_{0,K}^n > 0$ for all $K \in \mathcal{T}$.

With the positivity of $u_{0,K}^n$, $K \in \mathcal{T}$, at hand, let us focus on the $u_{i,K}^n$ for an arbitrary $i = 1, \dots, I$. Similarly, we assume that there exists some $K \in \mathcal{T}$ such that $u_{i,K}^n \le 0$. Then owing to (37), we infer that $F_{i,K\sigma}^n \le 0$ for all $\sigma \in \mathcal{E}_K$, and then that

$$ 0 \ge u_{i,K}^n = u_{i,K}^{n-1} - \frac{\tau^n}{m_K} \sum_{\sigma \in \mathcal{E}_K} F_{i,K\sigma}^n \ge 0. $$
This leads to $u_{i,K}^n = 0$ and to $F_{i,K\sigma}^n = 0$ for all $\sigma \in \mathcal{E}_K$. Since we already know that $u_{0,K}^n > 0$, we deduce from (37) that $u_{i,L}^n \le 0$ for all cells $L$ sharing a face $\sigma = K|L$ with $K$. As above, this implies as $u_{0,K}^n = 0$ for all $K \in \mathcal{T}$, which contradicts (36). Then $u_{i,K}^n > 0$ for all $K \in \mathcal{T}$, concluding the proof of Proposition 3.2. $\square$

A consequence of the previous proposition is that a solution to the modified scheme with (37) instead of (26) is also a solution to the original scheme (24)-(27).

As we did assume that the background charge density $f$ and thus its discrete counterpart $(f_K)_{K \in \mathcal{T}}$ are uniformly bounded, and that $\phi^D$ belongs to $L^\infty \cap H^{1/2}(\Gamma^D)$, we deduce a uniform bound for the electric potential.

**Proposition 3.3 (Uniform bound for the electric potential).** There exists $C = C(\phi^D, \Omega, \lambda, f, (z_i)_i, \zeta_{\mathcal{T}})$ such that

$$ (39) \qquad \max_{K \in \mathcal{T}} |\phi_K^n| + \left( \sum_{\sigma \in \mathcal{E}} a_\sigma (\phi_K^n - \phi_{K\sigma}^n)^2 \right)^{1/2} \le C, \quad n \ge 1. $$

The proof of this result can be found in [14, Proposition A.1] and is a consequence of the uniform boundedness of the right-hand side of the discrete Poisson equation (24).

These a priori estimates are sufficient to prove the existence of a solution to the scheme. We end up with the following proposition.

**Proposition 3.4 (Existence of solutions).** There exists at least one solution to the numerical scheme (24)-(27) such that $u_{i,K}^n > 0$ for all $i = 0, \dots, I$, for all $K \in \mathcal{T}$ and all $n \ge 1$.

*Proof.* The proof is based on an inductive procedure in time. At each time step $n \ge 1$, we use a topological degree argument to show the existence of a solution $(U_K^n, \phi_K^n)_{K \in \mathcal{T}} \in (0,1)^{I \times \mathcal{T}} \times [-C, C]^\mathcal{T}$ to (24)-(25), with the fluxes defined as in (37). Because of the positivity of the discrete volume fractions stated in Proposition 3.2 and to the equivalence between (37) and (26) for such positive volume fractions, this ensures the existence for the original scheme (24)-(26).

Let $n \ge 1$ be such that $(U_K^{n-1}, \phi_K^{n-1})_{K \in \mathcal{T}} \in [0,1]^{I \times \mathcal{T}} \times [-C, C]^\mathcal{T}$ is given (that is the case for $n = 1$, thanks to (29) and (7)). The idea is to deform our non-linear system continuously until it is transformed into a linear one with known solutions. Therefore, let introduce a parameter $s \in [0, \tau^n]$ and define, for every $K \in \mathcal{T}$,

$$ (40) \qquad \begin{cases} m_K (u_{i,K}^{(s)} - u_{i,K}^{n-1}) + s \sum_{\sigma \in \mathcal{E}_K} F_{i,K\sigma}^{(s)} = 0, & i = 1, \dots, I, \\ \lambda^2 \sum_{\sigma \in \mathcal{E}_K} a_\sigma (\phi_K^{(s)} - \phi_{K\sigma}^{(s)}) = m_K (f_K + \sum_{i=1}^I z_i u_{i,K}^{(s)}), & \end{cases} $$

with $F_{i,K\sigma}^{(s)}$ defined by (37), which corresponds to the original schemes when $s = \tau^n$. For $s = 0$, the two equations can be decoupled and the first one can be rewritten in matrix form

$$ \mathbb{M}(U^{(0)} - U^{n-1}) = 0, $$

where $U^{(s)} = (U_K^{(s)})_{K \in \mathcal{T}}$, for $s \in [0, \tau^n]$, $U^{n-1} = (U_K^{n-1})_{K \in \mathcal{T}}$ and $\mathbb{M} = \text{diag}((m_K)_{K \in \mathcal{T}})$ is a positive definite matrix. Therefore, there exists a unique solution $U^{(0)} = U^{n-1} \in [0,1]^{I \times \mathcal{T}}$. Via the Proposition 3.3, there also exists a unique $\phi^{(s)} = (\phi_K^{(s)})_{K \in \mathcal{T}} \in [-C, C]^\mathcal{T}$. Moreover, thanks to the continuity of the discrete fluxes, the functional

$$ S : \begin{cases} [0, \tau^n] \times [-2, 2]^{I \times \mathcal{T}} \times [-C-1, C+1]^\mathcal{T} \to \mathbb{R}^{\mathcal{T}} \times \mathbb{R}^{\mathcal{T}} \\ (s, U, \phi) \mapsto S(s, U, \phi) \end{cases} $$

whose zeros $(s, U^{(s)}, \phi^{(s)})$ are the solutions of (40), is hence an homotopy. Furthermore, all along it, its zeros $(U^{(s)}, \phi^{(s)})$ remain inside the compact subset $[0, 1]^{I \times \mathcal{T}} \times [-C, C]^\mathcal{T}$.

Thus, the topological degree corresponding to $S(s, U, \phi) = 0$ is equal to one, all along the homotopy, and hence in particular for $s = \tau^n$. That implies the existence of a solution to the scheme (but it does not guarantee uniqueness). $\square$
We note that the proof proposed here is simpler than that found in, for example, [14, Proposition 3.2]. The reason lies in the fact that we did not use an entropy inequality to derive the uniform estimates on volume fractions and, therefore, we do not have to guarantee the validity of it during the homotopy.

**Remark 3.5** (On the uniqueness of solutions). The uniqueness of weak solutions for (1)-(7), as well as for numerical schemes found in the literature, remains largely unresolved. This challenge is inherent to the cross-diffusion structure of the problem. In [31], uniqueness is established under strong assumptions on the diffusion coefficients and valences, specifically $D_i \equiv D$ and $z_i \equiv z$ for all $i = 1, \dots, I$.

For the numerical scheme (25)-(29), by following the argument in [7, Section 3.2, Uniqueness of solutions], uniqueness of discrete solutions can also be shown under the same assumptions: constant diffusion coefficients $D_i \equiv D$ and the absence of a drift term in the fluxes. The assumption of constant diffusion coefficients enables a preliminary proof of uniqueness for the discrete variable $u_0$. In our scheme, since no drift term is present and $\mathfrak{B}(0) = 1$, the function $\mathfrak{B}$ does not appear in the definition of the fluxes. Consequently, the equation for the discrete $u_0$ is significantly simpler than that in [7], and uniqueness of $u_0$ follows straightforwardly. The absence of a drift term also allows the application of Gajewski's entropy method to deduce uniqueness for the remaining components $(u_i)_{1 \le i \le I}$, as in [7].

**3.2. Energy dissipation at the discrete level.** The next proposition is about the thermodynamical consistency of our scheme and the decay of a discrete counterpart of the free energy.

**Proposition 3.6** (Free energy decay). Let $(U_K^n, \phi_K^n)_{K \in \mathcal{T}, n \ge 1}$ be a solution to the scheme (24)-(27) as in Proposition 3.4, then define for $n \ge 0$ the discrete free energy at the $n^{th}$ time step

$$ (41) \qquad \mathcal{H}_{\mathcal{T}}^{n} = \sum_{K \in \mathcal{T}} m_K H(U_K^n) + \frac{\lambda^2}{2} \sum_{\sigma \in \mathcal{E}} a_{\sigma} (\phi_K^n - \phi_{K\sigma}^n)^2 - \lambda^2 \sum_{\sigma \in \mathcal{E}^D} a_{\sigma} \phi_{\sigma}^D (\phi_{\sigma}^D - \phi_K^n), $$

the discrete electrochemical potentials $\mu_{i,K}^n = \log \left( \frac{u_{i,K}^n}{u_{0,K}^n} \right) + z_i \phi_K^n$ of species $i$, and

$$ \mathcal{D}_{\mathcal{T}}^{n} = \sum_{i=1}^{I} \sum_{\sigma \in \mathcal{E}_{\text{int}}} F_{i,K\sigma}^{n} (\mu_{i,K}^{n} - \mu_{i,L}^{n}) $$

the discrete dissipation, which vanishes if and only if $(\mu_{i,K}^n)_{K \in \mathcal{T}}$ is constant for all $i = 1, \dots, I$. Then there holds

$$ (42) \qquad \mathcal{H}_{\mathcal{T}}^{n} + \tau^{n} \mathcal{D}_{\mathcal{T}}^{n} \le \mathcal{H}_{\mathcal{T}}^{n-1}, \quad n \ge 1. $$

*Proof.* As a consequence of the positivity of $u_{0,K}^n$ and of the monotonicity of the exponential function, one then easily infers from reformulation (31) of the discrete fluxes that

$$ (43) \qquad \mathcal{D}_{i,\sigma}^n := F_{i,K\sigma}^n (\mu_{i,K}^n - \mu_{i,L}^n) \ge 0, \quad \forall i = 1, \dots, I, \sigma = K | L \in \mathcal{E}_{\text{int}}, $$

whence the non-negativity of $\mathcal{D}^n$. As $u_{0,K}^n > 0$ for all $K \in \mathcal{T}$ and $n \ge 1$, and as $y \mapsto e^y$ is strictly increasing, one gets that $\mathcal{D}_{i,\sigma}^n = 0$ iff $\mu_{i,K}^n = \mu_{i,L}^n$.

To proceed, we first define by $\mu_{i,K}^n = \log\left(\frac{u_{i,K}^n}{u_{0,K}^n}\right) + z_i \phi_K^n = \log(w_{i,K}^n)$ the electrochemical potential of species $i$, then we multiply the discrete conservation law (25) by $\tau^n \mu_{i,K}^n$ and sum over $i = 1, \dots, I$ and $K \in \mathcal{T}$. Thanks to the local conservativity of the fluxes, i.e., $F_{i,K\sigma}^n + F_{i,L\sigma}^n = 0$ for all $\sigma = K | L \in \mathcal{E}_{\text{int}}$, and to the no-flux boundary condition $F_{i,K\sigma}^n = 0$ for $\sigma \in \mathcal{E}_{\text{ext}}$, $i = 1, \dots, I$, a discrete integration by parts in space allows the flux terms to be rewritten as a sum over interior edges only, namely $\sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K} F_{i,K\sigma}^n \mu_{i,K}^n = \sum_{\sigma \in \mathcal{E}_{\text{int}}} \mathcal{D}_{i,\sigma}^n$ for all $i = 1, \dots, I$. This leads to the identity

$$ (44) \qquad \mathcal{A}_{\mathcal{T}}^{n} + \mathcal{B}_{\mathcal{T}}^{n} + \tau^{n} \mathcal{D}_{\mathcal{T}}^{n} = 0. $$
where we have set

$$
\begin{align*}
\mathcal{A}_{\mathcal{T}}^{n} &= \sum_{i=1}^{I} \sum_{K \in \mathcal{T}} (u_{i,K}^{n} - u_{i,K}^{n-1}) \log \left( \frac{u_{i,K}^{n}}{u_{0,K}^{n}} \right) m_{K} \\
&\stackrel{(27)}{=} \sum_{i=0}^{I} \sum_{K \in \mathcal{T}} (u_{i,K}^{n} - u_{i,K}^{n-1}) \frac{\partial \tilde{H}}{\partial u_{i}} (u_{0,K}^{n}, \dots, u_{I,K}^{n}) m_{K},
\end{align*}
$$

where $\tilde{H} : \{(u_i)_{0 \le i \le I} \mid \sum_{i=0}^I u_i = 1\} \to [-\log(I+1), +\infty)$, defined by $\tilde{H}(u_0, \dots, u_I) := \sum_{i=0}^I \log(u_i)$, and

$$
\begin{align*}
\mathcal{B}_{\mathcal{T}}^{n} &= \sum_{i=1}^{I} \sum_{K \in \mathcal{T}} (u_{i,K}^{n} - u_{i,K}^{n-1}) z_{i} \phi_{K}^{n} m_{K} \\
&\stackrel{(24)}{=} \lambda^{2} \sum_{K \in \mathcal{T}} \phi_{K}^{n} \sum_{\sigma \in \mathcal{E}_{K}} a_{\sigma} (\phi_{K}^{n} - \phi_{K}^{n-1} - (\phi_{K\sigma}^{n} - \phi_{K\sigma}^{n-1})).
\end{align*}
$$

From the convexity of $\tilde{H}$ and the identity $\tilde{H}(u_0, \dots, u_I) = H(U)$, we deduce that

$$
(45) \qquad \mathcal{A}_{\mathcal{T}}^{n} \geq \sum_{K \in \mathcal{T}} (H(U_{K}^{n}) - H(U_{K}^{n-1})) m_{K},
$$

while the term $\mathcal{B}_{\mathcal{T}}^n$ can be rewritten by observing that each interior edge is counted twice in the sum $\sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K}$. The resulting equivalent reformulation is

$$
\begin{align*}
\mathcal{B}_{\mathcal{T}}^{n} &= \lambda^2 \sum_{\sigma \in \mathcal{E}} a_{\sigma} (\phi_{K}^{n} - \phi_{K}^{n-1} - (\phi_{K\sigma}^{n} - \phi_{K\sigma}^{n-1})) (\phi_{K}^{n} - \phi_{K\sigma}^{n}) \\
&\quad + \lambda^2 \sum_{\sigma \in \mathcal{E}^D} a_{\sigma} \phi_{\sigma}^{D} (\phi_{K}^{n} - \phi_{K}^{n-1}),
\end{align*}
$$

which is valid since boundary edges $\sigma \in \mathcal{E}^N$ do not contribute to the first sum, and the contributions from boundary edges $\sigma \in \mathcal{E}^D$ cancel with the second part of the new expression of $\mathcal{B}_{\mathcal{T}}^n$. Then using the elementary convexity inequality $a(a-b) \ge (a^2-b^2)/2$ in the above term and combining the result with (45) in (44) provides the desired result (42). $\square$

Proposition 3.6 allows us to complete the proof of Theorem 2.1, but it also contains important information for proving the convergence of the scheme. It is the purpose of the next section to extract this information.

**3.3. Further uniform estimates on the discrete solution.** To pass to the limit in the scheme and to identify the limit as a weak solution, we need to extract some further uniform estimates, in particular the discrete $L_{\text{loc}}^2(H^1)$ estimates on the discrete counterparts of $u_0$ and $\sqrt{u_i u_0}$. We prove these estimates in Lemma 3.8. As an intermediate result we need a uniform bound on the discrete free energy.

**Lemma 3.7 (Uniform bound for the free energy).** There exists $C > 0$ depending only on $\Omega$, $\phi^D$, $\lambda$, $f$, $(z_i)_i$, and $\zeta_T$ such that, for all $N \ge 1$, there holds $|\mathcal{H}_T^N| \le C$.

*Proof.* Because of the bound $0 \le u_{i,K}^n \le 1$ for all $i$ and $K$, it is clear that the first two contributions of (41) remain uniformly bounded by a constant depending only on $\Omega$. Indeed, for the first term, this is a direct consequence of the definition (10) of $H$, whereas the second term can be controlled thanks to Proposition 3.3. Concerning the last contribution observe that if one defines $\phi_K^D$ and $\phi_\sigma^D$ as the averages of $\phi^D$ on $K \in \mathcal{T}$ and $\sigma \in \mathcal{E}$ respectively, then a suitable reorganization of the sum shows that

$$
\begin{align*}
(46) \quad & \sum_{\sigma \in \mathcal{E}^D} a_{\sigma} \phi_{\sigma}^D (\phi_K^n - \phi_{\sigma}^D) = \sum_{\sigma \in \mathcal{E}} a_{\sigma} (\phi_{\sigma}^D - \phi_K^D) (\phi_K^n - \phi_{K\sigma}^n) + \sum_{K \in \mathcal{T}} \phi_K^D \sum_{\sigma \in \mathcal{E}_K} a_{\sigma} (\phi_K^n - \phi_{K\sigma}^n) \\
& \stackrel{(24)}{=} \sum_{\sigma \in \mathcal{E}} a_{\sigma} (\phi_{\sigma}^D - \phi_K^D) (\phi_K^n - \phi_{K\sigma}^n) + \frac{1}{\lambda^2} \sum_{K \in \mathcal{T}} m_K \phi_K^D \left( f_K + \sum_{i=1}^I z_i u_{i,K}^n \right).
\end{align*}
$$
Using Young's inequality and $\mathcal{E}^D \subset \mathcal{E}$ for the first term in the above right-hand side gives

$$ \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}^D} a_{\sigma}(\phi_{\sigma}^{D} - \phi_{K}^{D})(\phi_{K}^{n} - \phi_{K\sigma}^{n}) \leq \sum_{\sigma \in \mathcal{E}} a_{\sigma}(\phi_{\sigma}^{D} - \phi_{K}^{D})^{2} + \frac{1}{4} \sum_{\sigma \in \mathcal{E}} a_{\sigma}(\phi_{K}^{n} - \phi_{K\sigma}^{n})^{2}. $$

Then Lemma 9.4 of [26] shows that

$$ \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}^D} a_{\sigma}(\phi_{\sigma}^D - \phi_{K}^D)^2 \leq C(\zeta_{\mathcal{T}}) \|\nabla \phi^D\|_{L^2(\Omega)}^2. $$

Then since $f$ and $U$ are bounded, we deduce from (46) that there exists $C$ depending only on $\|f\|_{\infty}$, $\max_i |z_i|$, $\Omega$, $\|\phi^D\|_{H^1}$, $\lambda$, and $\zeta_{\mathcal{T}}$ such that

$$ \sum_{\sigma \in \mathcal{E}^D} a_{\sigma} \phi_{\sigma}^{D} (\phi_{K}^{n} - \phi_{\sigma}^{D}) \leq C. $$

The result of the lemma follows.

□

The next lemma shows that the control of the energy dissipation $\sum_{n \geq 1} \tau^n \mathcal{D}_{\mathcal{T}}^n$ inferred from Proposition 3.6 gives some $L_{\text{loc}}^2(H^1)$ type control on the discrete counterparts of $u_0$, $\sqrt{u_0}$, and $\sqrt{u_i u_0}$.

**Lemma 3.8** (Uniform bounds on $u_0$, $\sqrt{u_0}$, and $\sqrt{u_i u_0}$). There exists $C > 0$ depending only on $\Omega$, $\phi^D$, $\lambda$, $f$, $(z_i)_i$, $(D_i)_i$, and $\zeta_{\mathcal{T}}$ such that, for all $N \geq 1$, there holds

$$ (47) \quad \sum_{n=1}^{N} \tau^n \sum_{i=1}^{I} \sum_{\sigma \in \mathcal{E}_{\text{int}}} a_{\sigma} \left( \sqrt{u_{i,K}^n u_{0,K}^n} - \sqrt{u_{i,L}^n u_{0,L}^n} \right)^2 \\ \qquad + \sum_{n=1}^{N} \tau^n \sum_{\sigma \in \mathcal{E}_{\text{int}}} a_{\sigma} \left( \sqrt{u_{0,K}^n} - \sqrt{u_{0,L}^n} \right)^2 \\ \qquad + \sum_{n=1}^{N} \tau^n \sum_{\sigma \in \mathcal{E}_{\text{int}}} a_{\sigma} (u_{0,K}^n - u_{0,L}^n)^2 \leq C(1+t^N). $$

As a consequence, one also has

$$ (48) \quad \sum_{n=1}^{N} \tau^n \sum_{i=1}^{I} \sum_{\sigma \in \mathcal{E}_{\text{int}}} a_{\sigma} \left( u_{i,K}^n \sqrt{u_{0,K}^n} - u_{i,L}^n \sqrt{u_{0,L}^n} \right)^2 \leq C(1+t^N) $$

and

$$ (49) \quad \sum_{n=1}^{N} \tau^n \sum_{i=1}^{I} \sum_{\sigma \in \mathcal{E}_{\text{int}}} a_{\sigma} (u_{i,K}^n u_{0,K}^n - u_{i,L}^n u_{0,L}^n)^2 \leq C(1+t^N). $$

*Proof.* One gets from the elementary inequality $(a-b)(\log(a) - \log(b)) \geq 4(\sqrt{a} - \sqrt{b})^2$ applied to (43) that

$$ \mathcal{D}_{i,\sigma}^n \geq 4a_{\sigma} D_i \mathfrak{R}(e^{-z_i \phi_K^n}, e^{-z_i \phi_L^n}) \left( \sqrt{u_{i,K}^n u_{0,L}^n} e^{\frac{z_i}{4}(\phi_K^n - \phi_L^n)} - \sqrt{u_{i,L}^n u_{0,K}^n} e^{\frac{z_i}{4}(\phi_K^n - \phi_L^n)} \right)^2 $$

with $\mathfrak{R}(e^{-z_i \phi_K^n}, e^{-z_i \phi_L^n}) = \mathfrak{R}(e^{-z_i \phi_K^n}, e^{-z_i \phi_L^n}) e^{\frac{z_i}{2}(\phi_K^n + \phi_L^n)}$. Thanks to (39) and since $D_i > 0$ for all $i$, there holds

$$ 2D_i \mathfrak{R}(e^{-z_i \phi_K^n}, e^{-z_i \phi_L^n}) \geq \kappa $$

for some $\kappa > 0$ uniform w.r.t. $K$, $i$, and $n$. As a consequence, using furthermore that $(a+b)^2 \geq \frac{1}{2}a^2 - b^2$,

$$ \mathcal{D}_{i,\sigma}^n \geq \kappa a_{\sigma} \cosh^2\left(\frac{z_i}{4}(\phi_K^n - \phi_L^n)\right) \left(\sqrt{u_{i,K}^n u_{0,L}^n} - \sqrt{u_{i,L}^n u_{0,K}^n}\right)^2 \\ - 2\kappa a_{\sigma} \left(\sqrt{u_{i,K}^n u_{0,L}^n} + \sqrt{u_{i,L}^n u_{0,K}^n}\right)^2 \sinh^2\left(\frac{z_i}{4}(\phi_K^n - \phi_L^n)\right). $$
Since $|\phi_K^n| \le C$ owing to (39), one has $\sinh^2(\frac{z_i}{4}(\phi_K^n - \phi_L^n)) \le C(\phi_K^n - \phi_L^n)^2$. Using moreover that $0 < u_{i,K}^n, u_{0,K}^n < 1$ and that $\cosh(a) \ge 1$, one gets that

$$ \mathcal{D}_{i,\sigma}^n \ge a_{\sigma} \kappa \left( \sqrt{u_{i,K}^n u_{0,L}^n} - \sqrt{u_{i,L}^n u_{0,K}^n} \right)^2 - C a_{\sigma} (\phi_K^n - \phi_L^n)^2. $$

Since

$$ (50) \quad \left( \sqrt{u_{i,K}^n u_{0,L}^n} - \sqrt{u_{i,L}^n u_{0,K}^n} \right)^2 = \left( \sqrt{u_{i,K}^n u_{0,K}^n} - \sqrt{u_{i,L}^n u_{0,L}^n} \right)^2 - (u_{i,K}^n - u_{i,L}^n)(u_{0,K}^n - u_{0,L}^n), $$

then summing over $i = 1, \dots, I$ and $\sigma = K|L$ and using (27) leads to

$$ \mathcal{D}_T^n \ge \kappa \sum_{i=1}^{I} \sum_{\sigma \in \mathcal{E}_{\text{int}}} a_{\sigma} \left( \sqrt{u_{i,K}^n u_{0,K}^n} - \sqrt{u_{i,L}^n u_{0,L}^n} \right)^2 + \kappa \sum_{\sigma \in \mathcal{E}_{\text{int}}} a_{\sigma} (u_{0,K}^n - u_{0,L}^n)^2 - C \sum_{\sigma \in \mathcal{E}_{\text{int}}} a_{\sigma} (\phi_K^n - \phi_{K\sigma}^n)^2. $$

Bearing in mind Proposition 3.3, we obtain that

$$ (51) \quad \mathcal{D}_T^n \ge \kappa \sum_{i=1}^{I} \sum_{\sigma \in \mathcal{E}_{\text{int}}} a_{\sigma} \left( \sqrt{u_{i,K}^n u_{0,K}^n} - \sqrt{u_{i,L}^n u_{0,L}^n} \right)^2 + \kappa \sum_{\sigma \in \mathcal{E}_{\text{int}}} a_{\sigma} (u_{0,K}^n - u_{0,L}^n)^2 - C. $$

Moreover, since $\sqrt{ab} \le (a+b)/2$ for all $a, b \ge 0$, there holds

$$ \sum_{i=0}^{I} \sqrt{u_{i,K}^n u_{i,L}^n} \le \sum_{i=0}^{I} \frac{u_{i,K}^n + u_{i,L}^n}{2} = 1, \quad \forall \sigma = K|L \in \mathcal{E}_{\text{int}}. $$

This provides

$$ \sum_{i=1}^{I} \sum_{\sigma \in \mathcal{E}_{\text{int}}} a_{\sigma} \left( \sqrt{u_{i,K}^n u_{0,K}^n} - \sqrt{u_{i,L}^n u_{0,L}^n} \right)^2 \\ \ge \sum_{\sigma \in \mathcal{E}_{\text{int}}} a_{\sigma} \left( (1-u_{0,K}^n)u_{0,K}^n + (1-u_{0,L}^n)u_{0,L}^n - 2(1-\sqrt{u_{0,K}^n u_{0,L}^n})\sqrt{u_{0,K}^n u_{0,L}^n} \right) \\ = \sum_{\sigma \in \mathcal{E}_{\text{int}}} a_{\sigma} \left( \sqrt{u_{0,K}^n} - \sqrt{u_{0,L}^n} \right)^2 - \sum_{\sigma \in \mathcal{E}_{\text{int}}} a_{\sigma} (u_{0,K}^n - u_{0,L}^n)^2, $$

whence we also deduce that

$$ \mathcal{D}_T^n \ge \kappa \sum_{\sigma \in \mathcal{E}_{\text{int}}} a_{\sigma} \left( \sqrt{u_{0,K}^n} - \sqrt{u_{0,L}^n} \right)^2 - C. $$

To deduce (47), it eventually remains to remark from (42) and Lemma 3.7 that there exists $C$ depending on none of $h$, $\Delta t$, $N$ nor on the initial data $U^0 = (u_i^0)_{1 \le i \le I}$ (provided it fulfills (7)) such that $\sum_{n=1}^N \tau^n \mathcal{D}_T^n \le C$. Combining this with (51) yields (47).

To recover estimate (48), remark that

$$ \left( \sqrt{u_{0,K}^n u_{i,K}^n} - \sqrt{u_{0,L}^n u_{i,L}^n} \right) = \left[ \sqrt{u_{i,K}^n} \left( \sqrt{u_{0,K}^n u_{i,K}^n} - \sqrt{u_{0,L}^n u_{i,L}^n} \right) \right. \\ \left. + u_{i,L}^n \left( \sqrt{u_{0,K}^n} - \sqrt{u_{0,L}^n} \right) + \sqrt{u_{i,L}^n} \left( \sqrt{u_{0,L}^n u_{i,K}^n} - \sqrt{u_{0,K}^n u_{i,L}^n} \right) \right]. $$

As $(a+b+c)^2 \le 3(a^2 + b^2 + c^2)$, and since $0 \le u_{i,K}^n \le 1$, one gets that

$$ \frac{1}{3} \left( \sqrt{u_{0,K}^n u_{i,K}^n} - \sqrt{u_{0,L}^n u_{i,L}^n} \right)^2 \le \left( \sqrt{u_{0,K}^n u_{i,K}^n} - \sqrt{u_{0,L}^n u_{i,L}^n} \right)^2 $$
$$ + (\sqrt{u_{0,K}^n} - \sqrt{u_{0,L}^n})^2 + (\sqrt{u_{0,L}^n u_{i,K}^n} - \sqrt{u_{0,K}^n u_{i,L}^n})^2. $$

The two first terms in the right-hand side are controlled thanks to (47), while, bearing in mind (50), the third term can be overestimated as follows:

$$ \begin{aligned} \left(\sqrt{u_{0,L}^n u_{i,K}^n} - \sqrt{u_{0,K}^n u_{i,L}^n}\right)^2 &\le \sum_{i=1}^I \left(\sqrt{u_{0,L}^n u_{i,K}^n} - \sqrt{u_{0,K}^n u_{i,L}^n}\right)^2 \\ &\le \sum_{i=1}^I \left(\sqrt{u_{0,K}^n u_{i,K}^n} - \sqrt{u_{0,L}^n u_{i,L}^n}\right)^2 + (u_{0,K}^n - u_{0,L}^n)^2, \end{aligned} $$

which again, after summation on $n = 1, \dots, N$ and $\sigma \in \mathcal{E}_{\text{int}}$, can be controlled thanks to (47).

Finally, to establish (49), remark that

$$ u_{i,K}^n u_{0,K}^n - u_{i,L}^n u_{0,L}^n = \left(\sqrt{u_{i,K}^n u_{0,K}^n} - \sqrt{u_{i,L}^n u_{0,L}^n}\right) \left(\sqrt{u_{i,K}^n u_{0,K}^n} + \sqrt{u_{i,L}^n u_{0,L}^n}\right). $$

Since $0 \le u_{0,K}^n, u_{i,K}^n \le 1$, one gets that

$$ (u_{i,K}^n u_{0,K}^n - u_{i,L}^n u_{0,L}^n)^2 \le 4 \left(\sqrt{u_{i,K}^n u_{0,K}^n} - \sqrt{u_{i,L}^n u_{0,L}^n}\right)^2. $$

Summing over $n = 1, \dots, N$ and $\sigma \in \mathcal{E}_{\text{int}}$ and using (47) gives the desired result. $\square$

One also deduces the following discrete $L_{\text{loc}}^2(L^2)^d$ estimates on the fluxes, which amount to some discrete $L_{\text{loc}}^2(H^1)'$ estimate on time increments of the discrete counterpart to $\partial_t u_i$.

**Lemma 3.9** (Uniform bounds on the fluxes). There exists $C$ depending only on $\Omega$, $\phi^D$, $\lambda$, $f$, $(z_i)_i$, $(D_i)_i$, and $\zeta_{\tau}$ such that

$$ (52) \qquad \sum_{i=1}^{I} \sum_{n=1}^{N} \tau^n \sum_{\sigma \in \mathcal{E}_{\text{int}}} \frac{d_{\sigma}}{m_{\sigma}} |F_{i,K\sigma}^n|^2 \le C(1+t^N). $$

*Proof.* To obtain discrete $L_{\text{loc}}^2(L^2)^d$ estimates on the fluxes, we need to exploit the discrete uniform estimates we have for $(U_K^n, \phi_K^n)_{K \in \mathcal{T}, n \ge 1}$. In this regard, we rewrite the fluxes in a different way before taking advantage of some properties of the Bernoulli function (28). One splits the flux (26) into two parts corresponding to convection and diffusion respectively:

$$ F_{i,K\sigma}^n = F_{i,K\sigma}^{\text{conv},n} + F_{i,K\sigma}^{\text{diff},n}, $$

with

$$ \begin{aligned} F_{i,K\sigma}^{\text{conv},n} &= a_{\sigma} D_i \frac{u_{i,K}^n u_{0,L}^n + u_{i,L}^n u_{0,K}^n}{2} [\mathfrak{B}(z_i(\phi_L^n - \phi_K^n)) - \mathfrak{B}(z_i(\phi_K^n - \phi_L^n))], \\ F_{i,K\sigma}^{\text{diff},n} &= a_{\sigma} D_i \frac{u_{i,K}^n u_{0,L}^n - u_{i,L}^n u_{0,K}^n}{2} [\mathfrak{B}(z_i(\phi_L^n - \phi_K^n)) + \mathfrak{B}(z_i(\phi_K^n - \phi_L^n))]. \end{aligned} $$

The flux $(F_{i,K\sigma}^n)_{\sigma,n}$ is bounded in $L_{\text{loc}}^2(L^2)^d$ in the sense of (52) if both $(F_{i,K\sigma}^{\text{conv},n})_{\sigma,n}$ and $(F_{i,K\sigma}^{\text{diff},n})_{\sigma,n}$ are.

As $\mathfrak{B}(-y) - \mathfrak{B}(y) = y$, one gets that

$$ (53) \qquad F_{i,K\sigma}^{\text{conv},n} = a_{\sigma} D_i \frac{u_{i,K}^n u_{0,L}^n + u_{i,L}^n u_{0,K}^n}{2} z_i(\phi_K^n - \phi_L^n). $$

The $L_{\text{loc}}^2(L^2)^d$ character of the above expression directly follows from the uniform bound on $u_{i,K}^n$, $0 \le i \le I$ and from the discrete $L^\infty(H^1)$ bound on $(\phi_K^n)_{K,n}$ stated in Proposition 3.3. Therefore,

$$ \sum_{i=1}^{I} \sum_{n=1}^{N} \tau^n \sum_{\sigma \in \mathcal{E}_{\text{int}}} \frac{d_{\sigma}}{m_{\sigma}} |F_{i,K\sigma}^{\text{conv},n}|^2 \le Ct^N. $$
Concerning the diffusive term, as $2 \le \mathfrak{B}(y) + \mathfrak{B}(-y) = y \coth \frac{y}{2} \le 2 + \frac{y^2}{6}$, one has that

$$ (54) \qquad F_{i,K\sigma}^{\text{diff},n} = a_\sigma D_i (u_{i,K}^n u_{0,L}^n - u_{i,L}^n u_{0,K}^n) \left(1 + \mathcal{O}\left((\phi_K^n - \phi_L^n)^2\right)\right). $$

From the discrete $L^\infty(H^1)$ bound on $(\phi_K^n)_{K,n}$, one can uniformly estimate the remainder. For the other term, since $u_{i,K}^n u_{0,L}^n - u_{i,L}^n u_{0,K}^n = u_{i,K}^n u_{0,K}^n - u_{i,L}^n u_{0,L}^n + (u_{i,K}^n + u_{i,L}^n)(u_{0,L}^n - u_{0,K}^n)$, then

$$ \begin{aligned} |u_{i,K}^n u_{0,L}^n - u_{i,L}^n u_{0,K}^n|^2 &\le C \left( |u_{i,K}^n u_{0,K}^n - u_{i,L}^n u_{0,L}^n|^2 + |(u_{i,K}^n + u_{i,L}^n)(u_{0,L}^n - u_{0,K}^n)|^2 \right) \\ &\le C \left( \left| \sqrt{u_{i,K}^n u_{0,K}^n} - \sqrt{u_{i,L}^n u_{0,L}^n} \right|^2 + |u_{0,L}^n - u_{0,K}^n|^2 \right), \end{aligned} $$

so that we obtain

$$ \sum_{i=1}^{I} \sum_{n=1}^{N} \tau^n \sum_{\sigma \in \mathcal{E}_{\text{int}}} \frac{d_\sigma}{m_\sigma} |F_{i,K\sigma}^{\text{diff},n}|^2 \le C \sum_{i=1}^{I} \sum_{n=1}^{N} \tau^n \sum_{\sigma \in \mathcal{E}_{\text{int}}} a_\sigma \left| \sqrt{u_{i,K}^n u_{0,K}^n} - \sqrt{u_{i,L}^n u_{0,L}^n} \right|^2 \\ + C \sum_{i=1}^{I} \sum_{n=1}^{N} \tau^n \sum_{\sigma \in \mathcal{E}_{\text{int}}} a_\sigma |u_{0,L}^n - u_{0,K}^n|^2, $$

thanks to the uniform bounds on $(U_K^n)_{K \in \mathcal{T}, n \le 1}$. Lemma 3.8 provides the desired $L_{\text{loc}}^2(L^2)$ bound on $F_{i,K\sigma}^{\text{diff},n}$, concluding the proof of Lemma 3.9. $\square$

From Lemma 3.9, we deduce the following discrete $L_{\text{loc}}^2((H^1)')$ estimate:

**Corollary 3.10.** There exists $C$ depending only on $\Omega$, $\phi^K$, $\lambda$, $f$, $(z_i)_i$, $(D_i)_i$, and $\zeta_T$ such that, for all $i = 0, \dots, I$ and all $\varphi_{T,\tau} = \sum_{K \in \mathcal{T}} \sum_{n=1}^N \varphi_K^n \mathbb{1}_{\{t^{n-1}, t^n\} \times K}$, one has

$$ (55) \qquad \sum_{n=1}^{N} \sum_{K \in \mathcal{T}} m_K (u_{i,K}^n - u_{i,K}^{n-1}) \varphi_K^n \le C (1+t^N)^{1/2} \left( \sum_{n=1}^{N} \tau^n \sum_{\sigma \in \mathcal{E}_{\text{int}}} a_\sigma (\varphi_K^n - \varphi_{K\sigma}^n)^2 \right)^{1/2}. $$

*Proof.* Let us first establish (55) for $i \ge 1$. Multiplying (25) by $\tau^n \varphi_K^n$, summing over $K \in \mathcal{T}$ and $n \in \{1, \dots, N\}$, and performing a discrete integration by parts on the contribution of the fluxes, namely $\sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K} F_{i,K\sigma}^n \varphi_K^n = \sum_{\sigma \in \mathcal{E}_{\text{int}}} F_{i,K\sigma}^n (\varphi_{K\sigma}^n - \varphi_K^n)$ for all $i = 1, \dots, I$, provides

$$ (56) \qquad \sum_{n=1}^{N} \sum_{K \in \mathcal{T}} m_K (u_{i,K}^n - u_{i,K}^{n-1}) \varphi_K^n = \sum_{n=1}^{N} \tau^n \sum_{\sigma \in \mathcal{E}_{\text{int}}} F_{i,K\sigma}^n (\varphi_{K\sigma}^n - \varphi_K^n). $$

Applying Cauchy–Schwarz inequality to the right-hand side then using Lemma 3.9 provides the desired result. The recovery of the estimate for $i = 0$ then directly follows from the definition (27) of $u_{0,K}^n$ and from (55) for $i \ge 1$. $\square$

### 4. CONVERGENCE OF THE SCHEMES

This section is devoted to the proof of Theorem 2.2, which relies on compactness arguments. Given sequences $(\mathcal{T}_\ell, \mathcal{E}_\ell, \{x_K\}_{K \in \mathcal{T}_\ell})_{\ell \ge 1}$ and $(\tau_\ell)_{\ell \ge 1} = ((\tau_\ell^n)_{n \ge 1})_{\ell \ge 1}$ of admissible meshes and a sequence of time steps in the sense of Section 2.1, with

$$ (57) \qquad \lim_{\ell \to \infty} h_{\mathcal{T}_\ell} = \lim_{\ell \to \infty} \Delta t_\ell = 0 \quad \text{and} \quad \zeta_{\mathcal{T}_\ell} \le \zeta_\star < +\infty, $$

we define the sequences $(u_{i,\mathcal{T}_\ell,\tau_\ell})_{\ell \ge 1}$ and $(\phi_{\mathcal{T}_\ell,\tau_\ell})_{\ell \ge 1}$ as in (33)&(34). For readability, the $\ell$ index is removed when not essential for understanding.

As usual in the analysis of finite volume schemes, we also need to handle quantities attached to the faces $\sigma \in \mathcal{E}$. To this end, we introduce the so-called diamond mesh of $\Omega$ by associating a diamond cell $\omega_\sigma$ to every $\sigma \in \mathcal{E}$. More precisely, $\omega_\sigma$ is the convex hull of $\{x_K, x_L, \sigma\}$ if $\sigma = K|L \in \mathcal{E}_{\text{int}}$, and $\omega_\sigma = \text{conv}\{x_K, \sigma\}$ if $\sigma \in \mathcal{E}_{\text{ext}} \cap \mathcal{E}_K$, see Figure 1 for an illustration. The Lebesgue measure $m_{\omega_\sigma}$ of $\omega_\sigma$ is then given by

$$ (58) \qquad m_{\omega_\sigma} = \frac{m_\sigma d_\sigma}{d}, \quad \sigma \in \mathcal{E}, $$
where d is the dimension of the space domain $\Omega$.

Among other quantities attached to faces, one defines the inflated fluxes $(F_{i,\varepsilon,\tau})_{1\le i\le I}$ as the piecewise constant in space and time vector fields defined by

$$F_{i,\varepsilon,\tau}(t,x) = \frac{d}{m_{\sigma}} F_{i,K\sigma}^{n} n_{K\sigma} \quad \text{if } (t,x) \in (t^{n-1}, t^n] \times \omega_{\sigma},$$

with $n_{K\sigma}$ the normal to $\sigma$ outward with respect to $K$. A major part of the analysis carried out in this section consists in showing that, up to the extraction of a subsequence,

$$F_{i,\varepsilon,\tau_{\ell}}(t,x) \xrightarrow{\ell \to +\infty} F_i \quad \text{weakly in } L_{\text{loc}}^2(\mathbb{R}_{\ge 0} \times \bar{\Omega})^d.$$

where we make use of the definition (15) of the continuous fluxes.

<figure><img src="image_2.png" /><figcaption>FIGURE 1. Examples of diamond cells $\omega_{\sigma}, \omega_{\sigma'}$ for inner and external faces $\sigma$ and $\sigma'$.</figcaption></figure>

Following [17, 25], to a piecewise constant in space function $v_{\mathcal{T}}(x) = \sum_{K \in \mathcal{T}} v_K 1_K(x)$, we associate the inflated gradient $\nabla_{\mathcal{T}} v_{\mathcal{T}} : \Omega \to \mathbb{R}^d$ defined by

$$(59a) \qquad \nabla_{\mathcal{T}} v_{\mathcal{T}}(x) = d \frac{v_{K\sigma} - v_K}{d_{\sigma}} n_{K\sigma} \quad \text{if } x \in \omega_{\sigma}.$$

This operator straightforwardly extends to piecewise constant functions in space and time:

$$(59b) \qquad \nabla_{\mathcal{T}} v_{\mathcal{T},\tau}(t,x) = d \frac{v_{K\sigma}^n - v_K^n}{d_{\sigma}} n_{K\sigma} \quad \text{if } (t,x) \in (t^{n-1}, t^n] \times \omega_{\sigma}.$$

This section is organized as follows. The needed compactness properties are established in Section 4.1, allowing to claim for the existence of (at least) a limit point $(U, \phi)$ to the sequence $(U_{\mathcal{T}_{\ell},\tau_{\ell}}, \phi_{\mathcal{T}_m,\Delta t_{\ell}})_{\ell \ge 1}$. Section 4.2 concludes the proof of Theorem 2.2 by showing that the limit value $(U, \phi)$ is a weak solution of the continuous problem (1)-(7) in the sense of Definition 1.1.

**4.1. Compactness of the approximate solutions.** We establish here enough compactness properties to pass to the limit in the scheme.

**Proposition 4.1 (Compactness results).** There exist functions $u_0 \in L^\infty(\mathbb{R}_{\ge 0} \times \Omega)$ with $\sqrt{u_0} \in L_{\text{loc}}^2(\mathbb{R}_{\ge 0}; H^1(\Omega))$, $u_1, \dots, u_n \in L^\infty(\mathbb{R}_{\ge 0} \times \Omega)$ with $u_i \sqrt{u_0} \in L_{\text{loc}}^2(\mathbb{R}_{\ge 0}; H^1(\Omega))$, and $\phi \in L^\infty(\mathbb{R}_{\ge 0} \times \Omega)$ with $\phi - \phi^D \in L^\infty(\mathbb{R}_{\ge 0}; V)$ such that, up to a subsequence, as $\ell \to +\infty$, there holds

$$(60) \qquad u_{i,\mathcal{T},\tau} \to u_i \quad \text{in the } L^\infty(\mathbb{R}_{\ge 0} \times \Omega)\text{-weak-\star sense, } 1 \le i \le I,$$

$$(61) \qquad u_{0,\mathcal{T},\tau} \to u_0 \quad \text{in the } L^\infty(\mathbb{R}_{\ge 0} \times \Omega)\text{-weak-\star sense and a.e. in } \mathbb{R}_{\ge 0} \times \Omega,$$

$$(62) \qquad \phi_{\mathcal{T},\tau} \to \phi \quad \text{in the } L^\infty(\mathbb{R}_{\ge 0} \times \Omega)\text{-weak-\star sense and a.e. in } \mathbb{R}_{\ge 0} \times \Omega.$$

Moreover, for $1 \le i \le I$, one has

$$(63) \qquad u_{i,\mathcal{T},\tau} \sqrt{u_{0,\mathcal{T},\tau}} \to u_i \sqrt{u_0} \quad \text{a.e. in } \mathbb{R}_{\ge 0} \times \Omega,$$

$$(64) \qquad u_{0,\mathcal{T},\tau} u_{i,\mathcal{T},\tau} \to u_0 u_i \quad \text{a.e. in } \mathbb{R}_{\ge 0} \times \Omega,$$

and also

$$(65) \qquad \nabla_{\mathcal{T},\tau} \sqrt{u_{0,\mathcal{T},\tau}} \to \nabla \sqrt{u_0} \quad \text{weakly in } L_{\text{loc}}^2(\mathbb{R}_{\ge 0}; L^2(\Omega)^d),$$

$$(66) \qquad \nabla_{\mathcal{T},\tau} [u_{0,\mathcal{T},\tau} u_{i,\mathcal{T},\tau}] \to \nabla [u_0 u_i] \quad \text{weakly in } L_{\text{loc}}^2(\mathbb{R}_{\ge 0}; L^2(\Omega)^d),$$
$$ (67) \qquad \nabla_{\mathcal{T},\tau} \phi_{\mathcal{T},\tau} \to \nabla \phi \quad \text{weakly in } L_{\text{loc}}^2(\mathbb{R}_{\ge 0}; L^2(\Omega)), $$

**Proof.** The $L^\infty(\mathbb{R}_{\ge 0} \times \Omega)$-weak-⋆ convergences (60)-(62) are the consequences of the uniform bounds on $u_{i,\mathcal{T},\tau}$ and $\phi_{\mathcal{T},\tau}$ in $L^\infty(\mathbb{R}_{\ge 0} \times \Omega)$, thanks to the Banach–Alaoglu Theorem. The almost everywhere convergence (61) follows from strong convergence given by the Aubin–Lions lemma, see [29] for a general presentation of the lemma, and [7, Lemma 9] for its application in our context. It can be applied thanks to the uniform estimates of Lemma 3.8 and Corollary 3.10. These estimates also yield the weak convergence in $L_{\text{loc}}^2(\mathbb{R}_{\ge 0}; L^2(\Omega))$ of $\nabla_{\mathcal{T},\tau} u_{0,\mathcal{T},\tau}$ towards $\nabla u_0$. Concerning the point-wise convergence of $\phi_{\mathcal{T},\tau}$ towards $\phi$, it can be proven by using the discrete $L_{\text{loc}}^2((H^1)^\prime)$ estimate on time increments of the right-hand side of the discrete Poisson equation (24) that follows from Corollary 3.10. As the proof is fully similar to the one of [14, Proposition 4.2], we do not provide details here. The weak compactness property (67) is also established in [14, Proposition 4.2].

In order to prove (65), one can either make use of the non-linear Aubin–Simon theorem of [1], or directly remark that, thanks to Lemma 3.8, the vector field $\nabla_{\mathcal{T}}\sqrt{u_{0,\mathcal{T},\tau}}$ is uniformly bounded in $L_{\text{loc}}^2(\mathbb{R}_{\ge 0}; L^2(\Omega)^d)$. We deduce from this boundedness that there exists $G_0 \in L_{\text{loc}}^2(\mathbb{R}_{\ge 0}; L^2(\Omega)^d)$ such that, up to a subsequence, $\nabla_{\mathcal{T}}\sqrt{u_{0,\mathcal{T},\tau}}$ tends to $G_0$ weakly. As $\sqrt{u_{0,\mathcal{T},\tau}}$ converges point-wise (and thus in $L_{\text{loc}}^1(\mathbb{R}_{\ge 0} \times \Omega))$ towards $\sqrt{u_0}$ owing to (61), the weak consistency of the inflated gradient $\nabla_{\mathcal{T}}$, see for instance [17, 22], allows to show that $G_0 = \nabla\sqrt{u_0}$.

Regarding the point-wise convergences (63) and (64), of $u_{i,\mathcal{T},\tau}\sqrt{u_{0,\mathcal{T},\tau}}$, respectively, we can conclude by applying [7, Lemma 10. (Discrete Aubin–Lions of “degenerate” type)]. To prove (63), we use it with $y_\ell = \sqrt{u_{0,\mathcal{T},\tau}}$ and $z_\ell = u_{i,\mathcal{T},\tau}$, taking advantage of the uniform estimate (55), (47), and (48). Whereas, with $y_\ell = u_{0,\mathcal{T},\tau}$ and $z_\ell = u_{i,\mathcal{T},\tau}$, we conclude (64) similarly, replacing (48) by (49).

It finally remains to prove (66) for $i = 1, \dots, I$. As a consequence of (49), $\nabla_{\mathcal{T}}[u_{i,\mathcal{T},\tau}u_{0,\mathcal{T},\tau}]$ is uniformly bounded in $L_{\text{loc}}^2(\mathbb{R}_{\ge 0} \times \bar{\Omega})^d$. Hence there exists $G_i$ such that, up to a subsequence,

$$ \nabla_{\mathcal{T}}[u_{i,\mathcal{T},\tau}u_{0,\mathcal{T},\tau}] \to G_i \quad \text{weakly in } L_{\text{loc}}^2(\mathbb{R}_{\ge 0} \times \bar{\Omega})^d. $$

As $u_{i,\mathcal{T},\tau}u_{0,\mathcal{T},\tau}$ converges towards $u_i u_0$ thanks to (64), we can invoke again the weak consistency of the inflated gradient to conclude that $G_i = \nabla[u_i u_0]$. $\square$

**4.2. Identification of the limit.** In this section, we conclude the proof of Theorem 2.2 by showing the following proposition.

**Proposition 4.2.** Let $(U, \phi)$ be as in Proposition 4.1, then it is a weak solution of the problem (1)-(7) in the sense of Definition 1.1.

**Proof.** The regularity requirements on $U$ and $\phi$ have already been checked in Proposition 4.1. Therefore it only remains to verify that the weak formulations (16) and (17) hold true. The case of the Poisson equation is classical. It will not be detailed here, and we refer to [14, Proposition 4.2] for a synthetic proof. We rather focus our attention on the derivation of (17).

Let $\varphi \in C_c^\infty(\mathbb{R}_{\ge 0} \times \bar{\Omega})$, then, for some admissible mesh $(\mathcal{T}_\ell, \mathcal{E}_\ell, (x_K)_{K \in \mathcal{T}_\ell})$ (we remove the subscript $\ell$ when possible for legibility), denote by $\varphi_K^n = \varphi(t^n, x_K)$. Multiplying (25) by $\tau^n \varphi_K^{n-1}$ and summing over $n \ge 1$ and $K \in \mathcal{T}$ yields, as for (56):

$$ (68) \qquad \sum_{n=1}^{N} \sum_{K \in \mathcal{T}} m_K(u_{i,K}^n - u_{i,K}^{n-1}) \varphi_K^{n-1} = \sum_{n=1}^{N} \tau^n \sum_{\sigma \in \mathcal{E}_{\text{int}}} F_{i,K\sigma}^n (\varphi_{K\sigma}^{n-1} - \varphi_K^{n-1}). $$

Since $\varphi$ is compactly supported in time, $\varphi_K^n = 0$ for $n$ large enough. Therefore, the left-hand side of (68) is rewritten as

$$ \sum_{n=1}^{N} \sum_{K \in \mathcal{T}} m_K(u_{i,K}^n - u_{i,K}^{n-1}) \varphi_K^{n-1} = \sum_{n=1}^{N} \tau^n \sum_{K \in \mathcal{T}} m_K u_{i,K}^n \frac{\varphi_K^{n-1} - \varphi_K^n}{\tau^n} - \sum_{K \in \mathcal{T}} m_K u_{i,K}^0 \varphi_K^n. $$

Since $\varphi$ is smooth, the approximate time derivative $\delta_\tau \varphi_{\mathcal{T},\tau}$ of $\varphi$ defined by

$$ \delta_\tau \varphi_{\mathcal{T},\tau}(t,x) = \frac{\varphi_K^n - \varphi_K^{n-1}}{\tau^n} \quad \text{if } (t,x) \in (t^{n-1}, t^n] \times K $$
converges in $L^1(\mathbb{R}_{\ge 0} \times \Omega)$ towards $\partial_t \varphi$. Therefore, using (60), we get that

$$ (69) \qquad \sum_{n=1}^{N} \tau^n \sum_{K \in \mathcal{T}} m_K u_{i,K}^n \frac{\varphi_K^{n-1} - \varphi_K^n}{\tau^n} \underset{\ell \to +\infty}{\longrightarrow} - \iint_{\mathbb{R}_{\ge 0} \times \Omega} u_i \partial_t \varphi. $$

The function $u_{i,\mathcal{T}}^0 = \sum_{K \in \mathcal{T}} u_{i,K}^0 1_K$ converges strongly in $L^1(\Omega)$ towards $u_i^0$, while $\varphi_{\mathcal{T}}^0 = \sum_{K \in \mathcal{T}} \varphi_K^0 1_K$ converges uniformly towards $\varphi(0, \cdot)$ thanks to the regularity of $\varphi$. Therefore,

$$ (70) \qquad - \sum_{K \in \mathcal{T}} m_K u_{i,K}^0 \varphi_K^0 \underset{\ell \to +\infty}{\longrightarrow} - \int_{\Omega} u_i^0 \varphi(0, \cdot), $$

so we can pass to the limit in the left-hand side of (68).

Let us now focus on its right-hand side. Define $\tilde{\varphi}_{\mathcal{T},\tau}$ by

$$ \tilde{\varphi}_{\mathcal{T},\tau}(t,x) = \varphi(t^{n-1}, x_K) \quad \text{if} \ (t,x) \in [t^{n-1}, t^n) \times K. $$

Following [23] (see [20] for a practical example), one can reconstruct a second approximate gradient operator $\hat{\nabla}_{\mathcal{T}}$ mapping the set of piecewise constant functions in space and time to $\mathbb{R}^d$ such that

$$ \hat{\nabla}_{\mathcal{T}} \tilde{\varphi}_{\mathcal{T},\tau}(t,x) \cdot n_{K\sigma} = \frac{\varphi_L^n - \varphi_K^n}{d_\sigma} \quad \text{if } (t,x) \in [t^{n-1}, t^n) \times \omega_\sigma, \quad \sigma = K|L \in \mathcal{E}_{\text{int}}, $$

and which is strongly consistent, i.e.,

$$ \hat{\nabla}_{\mathcal{T}} \tilde{\varphi}_{\mathcal{T},\tau} \underset{\ell \to +\infty}{\longrightarrow} \nabla \varphi \quad \text{strongly in } L^2(\mathbb{R}_{\ge 0} \times \Omega)^d. $$

The right-hand side of (68) then simply boils down to

$$ \sum_{n=1}^{N} \tau^n \sum_{\sigma \in \mathcal{E}_{\text{int}}} F_{i,K\sigma}^n (\varphi_{K\sigma}^{n-1} - \varphi_K^{n-1}) = \iint_{\mathbb{R}_{\ge 0} \times \Omega} F_{i,\mathcal{E},\tau} \cdot \hat{\nabla}_{\mathcal{T}} \tilde{\varphi}_{\mathcal{T},\tau}. $$

We have shown in Lemma 3.9 that $(F_{i,\mathcal{E},\tau\epsilon})_{\ell \ge 1}$ is uniformly bounded in $L_{\text{loc}}^2(\mathbb{R}_{\ge 0} \times \bar{\Omega})^d$. Therefore, it converges (up to a subsequence) towards some vector field $F_i$ weakly in $L_{\text{loc}}^2(\mathbb{R}_{\ge 0} \times \bar{\Omega})$. We can then pass to the limit in the right-hand side of (68) to obtain

$$ (71) \qquad \sum_{n=1}^{N} \tau^n \sum_{\sigma \in \mathcal{E}_{\text{int}}} F_{i,K\sigma}^n (\varphi_{K\sigma}^{n-1} - \varphi_K^{n-1}) \underset{\ell \to +\infty}{\longrightarrow} \iint_{\mathbb{R}_{\ge 0} \times \Omega} F_i \cdot \nabla \varphi. $$

To conclude the proof of Proposition 4.2, we still have to identify the limiting flux $F_i$ under the form (15). To this end, we split the inflated fluxes into a convective and diffusive part,

$$ F_{i,\mathcal{E},\tau}^{\text{conv}}(t,x) = \frac{d}{m_\sigma} F_{i,K\sigma}^{\text{conv},n} n_{K\sigma} \quad \text{if } (t,x) \in (t^{n-1}, t^n] \times \omega_\sigma, $$

$$ F_{i,\mathcal{E},\tau}^{\text{diff}}(t,x) = \frac{d}{m_\sigma} F_{i,K\sigma}^{\text{diff},n} n_{K\sigma} \quad \text{if } (t,x) \in (t^{n-1}, t^n] \times \omega_\sigma, $$

where $F_{i,K\sigma}^{\text{conv},n}$ and $F_{i,K\sigma}^{\text{diff},n}$ are as in the proof of Lemma 3.9. As both $F_{i,\mathcal{E},\tau}^{\text{conv}}$ and $F_{i,\mathcal{E},\tau}^{\text{diff}}$ have been shown to be bounded in $L_{\text{loc}}^2(\mathbb{R}_{\ge 0} \times \bar{\Omega})^d$, one has, up to a subsequence,

$$ (72) \qquad F_{i,\mathcal{E},\tau}^{\text{diff}} \underset{\ell \to +\infty}{\longrightarrow} F_i^{\text{diff}} \quad \text{and} \quad F_{i,\mathcal{E},\tau}^{\text{conv}} \underset{\ell \to +\infty}{\longrightarrow} F_i^{\text{conv}} \quad \text{weakly in } L_{\text{loc}}^2(\mathbb{R}_{\ge 0} \times \bar{\Omega})^d. $$

Let us first focus on the convective part $F_{i,\mathcal{E},\tau}^{\text{conv}}$. For $i=1, \dots, I$, $\sigma \in \mathcal{E}$ and $n \ge 1$, define

$$ \eta_{i,\sigma}^n = \begin{cases} \frac{u_{i,K}^n u_{0,L}^n + u_{i,L}^n u_{0,K}^n}{2} & \text{if } \sigma = K|L \in \mathcal{E}_{\text{int}}, \\ u_{i,K}^n u_{0,K}^n & \text{if } \sigma \in \mathcal{E}_K \cap \mathcal{E}_{\text{ext}}, \end{cases} $$

and

$$ \tilde{\eta}_{i,\sigma}^n = \begin{cases} \frac{u_{i,K}^n u_{0,K}^n + u_{i,L}^n u_{0,L}^n}{2} & \text{if } \sigma = K|L \in \mathcal{E}_{\text{int}}, \\ u_{i,K}^n u_{0,K}^n & \text{if } \sigma \in \mathcal{E}_K \cap \mathcal{E}_{\text{ext}}, \end{cases} $$
and then $\eta_{i,\mathcal{E},\tau}(t, x) = \eta_{i,\sigma}^n$ and $\tilde{\eta}_{i,\mathcal{E},\tau}(t, x) = \eta_{i,\sigma}^n$ if $(t, x) \in (t^{n-1}, t^n] \times \omega_\sigma$. In view of (53) and of the definition (59) of the inflated gradient, the convective part of the inflated approximate flux is rewritten as

$$ (73) \qquad F_{i, \mathcal{E}, \tau}^{\text{conv}} = -D_i \eta_{i, \mathcal{E}, \tau} z_i \nabla_{\mathcal{T}} \phi_{\mathcal{T}, \tau}. $$

Thanks to (64), proving that

$$ (74) \qquad \eta_{i, \mathcal{E}, \tau} - \tilde{\eta}_{i, \mathcal{E}, \tau} \xrightarrow{\ell \to +\infty} 0 \quad \text{a.e. in } \mathbb{R}_{\ge 0} \times \Omega $$

and

$$ (75) \qquad \tilde{\eta}_{i, \mathcal{E}, \tau} - u_{0, \mathcal{T}, \tau} u_{i, \mathcal{T}, \tau} \xrightarrow{\ell \to +\infty} 0 \quad \text{a.e. in } \mathbb{R}_{\ge 0} \times \Omega $$

is equivalent to showing that

$$ \eta_{i, \mathcal{E}, \tau} \xrightarrow{\ell \to +\infty} u_0 u_i \quad \text{a.e. in } \mathbb{R}_{\ge 0} \times \Omega, $$

and thus strongly in $L_{\text{loc}}^2(\mathbb{R}_{\ge 0} \times \bar{\Omega})$ thanks to the uniform bound $0 \le u_{i,\mathcal{T},\tau} \le 1$. Bearing in mind (67), this allows to pass to the weak limit in (73) and to recover that

$$ (76) \qquad F_i^{\text{conv}} = -D_i u_0 u_i z_i \nabla \phi. $$

To establish (74), remark that

$$ \|\eta_{i, \mathcal{E}, \tau} - \tilde{\eta}_{i, \mathcal{E}, \tau}\|_{L^2((0, t^N) \times \Omega)}^2 = \frac{1}{4} \sum_{n=1}^{N} \tau^n \sum_{\sigma \in \mathcal{E}_{\text{int}}} m_{\omega_\sigma} (u_{i,K}^n - u_{i,L}^n)^2 (u_{0,K}^n - u_{0,L}^n)^2. $$

Since $0 \le u_{i,K}^n \le 1$, owing to (58), and since $d_\sigma \le 2h_\mathcal{T}$ this gives

$$ \|\eta_{i, \mathcal{E}, \tau} - \tilde{\eta}_{i, \mathcal{E}, \tau}\|_{L^2((0, t^N) \times \Omega)}^2 \le \frac{h_{\mathcal{T}}^2}{d} \sum_{n=1}^{N} \tau^n \sum_{\sigma \in \mathcal{E}_{\text{int}}} a_\sigma (u_{0,K}^n - u_{0,L}^n)^2. $$

We can make use of estimate (47) to get that $\eta_{i,\mathcal{E},\tau} - \tilde{\eta}_{i,\mathcal{E},\tau}$ tends to 0 in $L_{\text{loc}}^2(\mathbb{R}_{\ge 0} \times \bar{\Omega})$, hence almost everywhere up to the extraction of yet another subsequence, whence (74).

Concerning (75), we can proceed similarly to obtain

$$ \|\tilde{\eta}_{i, \mathcal{E}, \tau} - u_{0, \mathcal{T}, \tau} u_{i, \mathcal{T}, \tau}\|_{L^2((0, t^N) \times \Omega)}^2 \le \frac{h_{\mathcal{T}}^2}{d} \sum_{n=1}^{N} \tau^n \sum_{\sigma \in \mathcal{E}_{\text{int}}} a_{\sigma} (u_{0,K}^n u_{i,K}^n - u_{0,L}^n u_{i,L}^n)^2 $$

which tends to 0 thanks to (49).

It finally remains to identify the limit of $F_{i,\mathcal{E},\tau}^{\text{diff}}$ as $F_i^{\text{diff}}$. To this end, remark first that, for $\sigma = K|L \in \mathcal{E}_{\text{int}}$ and $n \ge 1$, one has

$$ u_{i,K}^n u_{0,L}^n - u_{i,L}^n u_{0,K}^n = T_{1,\sigma}^n + T_{2,\sigma}^n + T_{3,\sigma}^n $$

with

$$ T_{1,\sigma}^n = u_{i,K}^n u_{0,K}^n - u_{i,L}^n u_{0,L}^n, $$

$$ T_{2,\sigma}^n = -2 \left( u_{i,K}^n \sqrt{u_{0,K}^n} + u_{i,L}^n \sqrt{u_{0,L}^n} \right) \left( \sqrt{u_{0,K}^n} - \sqrt{u_{0,L}^n} \right), $$

$$ T_{3,\sigma}^n = (u_{i,K}^n - u_{i,L}^n) \left( \sqrt{u_{0,K}^n} - \sqrt{u_{0,L}^n} \right)^2. $$

Therefore, in view of (54), we can split

$$ (77) \qquad F_{i, \mathcal{E}, \tau}^{\text{diff}} = -D_i \nabla_{\mathcal{T}}[u_{i, \mathcal{T}, \tau} u_{0, \mathcal{T}, \tau}] + 4D_i \gamma_{i, \mathcal{E}, \tau} \nabla_{\mathcal{T}} \sqrt{u_{0, \mathcal{T}, \tau}} + \mathcal{R}_{i, \mathcal{E}, \tau} + S_{i, \mathcal{E}, \tau}, $$

where we have set

$$ \gamma_{i, \mathcal{E}, \tau}(t, x) = \begin{cases} \frac{1}{2} \left( u_{i,K}^n \sqrt{u_{0,K}^n} + u_{i,L}^n \sqrt{u_{0,L}^n} \right) & \text{if } (t, x) \in (t^{n-1}, t^n] \times \omega_\sigma, \sigma = K|L \in \mathcal{E}_{\text{int}}, \\ u_{i,K}^n \sqrt{u_{0,K}^n} & \text{if } (t, x) \in (t^{n-1}, t^n] \times \omega_\sigma, \sigma \in \mathcal{E}_K \cap \mathcal{E}_{\text{ext}}, \end{cases} $$
and where, for $(t, x) \in (t^{n-1}, t^n] \times \omega_\sigma$, $\sigma = K|L \in \mathcal{E}_{\text{int}}$,

$$ \mathcal{R}_{i, \mathcal{E}, \tau}(t, x) = d D_i (u_{i,K}^n - u_{i,L}^n) \left( \sqrt{u_{0,K}^n} - \sqrt{u_{0,L}^n} \right)^2 n_{K\sigma}, $$

and, in view of (54),

$$ |S_{i, \mathcal{E}, \tau}(t, x)| \le \frac{C}{d_\sigma} (\phi_K^n - \phi_L^n)^2. $$

Both $\mathcal{R}_{i, \mathcal{E}, \tau}$ and $S_{i, \mathcal{E}, \tau}$ are assumed to vanish on $\mathbb{R}_{\ge 0} \times \omega_\sigma$ for all $\sigma \in \mathcal{E}_{\text{ext}}$. Then thanks to Proposition 3.3, one has

$$ \|S_{i, \mathcal{E}, \tau}\|_{L^1((0, t^N) \times \Omega)^d} \le C \sum_{n=1}^{N} \tau^n \sum_{\sigma \in \mathcal{E}_{\text{int}}} m_\sigma (\phi_K^n - \phi_L^n)^2 \le Ch_\tau t^N. $$

On the other hand, since $0 \le u_{i,K}^n \le 1$, we deduce from (47) that

$$ \|\mathcal{R}_{i, \mathcal{E}, \tau}\|_{L^1((0, t^N) \times \Omega)^d} \le D_i \sum_{n=1}^{N} \tau^n \sum_{\sigma \in \mathcal{E}_{\text{int}}} m_\sigma \left( \sqrt{u_{0,K}^n} - \sqrt{u_{0,L}^n} \right)^2 \le Ch_\tau (1 + t^N). $$

The two last terms in (77) then tend to 0, while the first term tends to $-D_i \nabla[u_0 u_i]$ thanks to (66). In view of the weak $L_{\text{loc}}^2(L^2)$ convergence of $\nabla_\tau \sqrt{u_{0,\tau}}$ towards $\nabla \sqrt{u_0}$, cf. (65), then it suffices to show that

$$ \gamma_{i, \mathcal{E}, \tau} \xrightarrow{\ell \to +\infty} u_i \sqrt{u_0} \quad \text{strongly in } L_{\text{loc}}^2(\mathbb{R}_{\ge 0} \times \bar{\Omega}) $$

to pass to the limit in (77). As $u_{i,\tau} \sqrt{u_{0,\tau}}$ converges strongly towards $u_i \sqrt{u_0}$, see (63), then

$$ \begin{aligned} \|\gamma_{i, \mathcal{E}, \tau} - u_{i, \mathcal{T}, \tau} \sqrt{u_{0, \mathcal{T}, \tau}}\|_{L^2((0, t^N) \times \Omega)}^2 &= \frac{1}{2} \sum_{n=1}^{N} \tau^n \sum_{\sigma \in \mathcal{E}_{\text{int}}} m_{\omega_\sigma} \left( u_{i,K}^n \sqrt{u_{0,K}^n} - u_{i,L}^n \sqrt{u_{0,L}^n} \right)^2 \\ &\le \frac{h_\tau^2}{d} \sum_{n=1}^{N} \tau^n \sum_{\sigma \in \mathcal{E}_{\text{int}}} a_\sigma \left( u_{i,K}^n \sqrt{u_{0,K}^n} - u_{i,L}^n \sqrt{u_{0,L}^n} \right)^2 \\ &\le \frac{h_\tau^2}{d} (1 + t^N) \xrightarrow{\ell \to +\infty} 0. \end{aligned} \tag{48} $$

This completes the proof of (72), thus the ones of Proposition 4.2 and Theorem 2.2. $\square$

## 5. LONG-TIME ASYMPTOTIC OF THE SCHEME

Our aim here is to prove our last main theoretical contribution, namely Theorem 2.3. Through-out this section, we will work on a fixed grid $(\mathcal{E}, \mathcal{T}, (x_K)_{K \in \mathcal{T}})$ and for a fixed time discretization. We are interested in steady states $(U_\tau^\infty, \phi_\tau^\infty)$ of the scheme (24)–(27), i.e. solutions to

$$ (78) \qquad \lambda^2 \sum_{\sigma \in \mathcal{E}_K} a_\sigma (\phi_K^\infty - \phi_{K\sigma}^\infty) = m_K \left( f_K + \sum_{i=1}^I z_i u_{i,K}^\infty \right), \quad K \in \mathcal{T}, $$

$$ (79) \qquad \sum_{\sigma \in \mathcal{E}_K} F_{i,K\sigma}^\infty = 0, \quad i = 1, \dots, I, \quad K \in \mathcal{T}, $$

with $F_{i,K\sigma}^\infty = 0$ if $\sigma \in \mathcal{E}_{\text{ext}}$ and

$$ (80) \qquad F_{i,K\sigma}^\infty = a_\sigma D_i \left( u_{i,K}^\infty u_{0,L}^\infty \mathfrak{B} \left( z_i (\phi_L^\infty - \phi_K^\infty) \right) - u_{i,L}^\infty u_{0,K}^\infty \mathfrak{B} \left( z_i (\phi_K^\infty - \phi_L^\infty) \right) \right) $$

for $\sigma = K|L \in \mathcal{E}_{\text{int}}$, where we have set

$$ (81) \qquad u_{0,K}^\infty = 1 - \sum_{i=1}^I u_{i,K}^\infty, \quad K \in \mathcal{T}. $$
The above equations have to be complemented with some discrete counterpart of the constraint on the mass (22), that is

$$ (82) \qquad \sum_{K \in \mathcal{T}} m_K u_{i,K}^{\infty} = \int_{\Omega} u_i^0, \quad i = 1, \dots, I. $$

We are moreover interested in the convergence of $(U_T^n, \phi_T^n)$ towards $(U_T^\infty, \phi_T^\infty)$ as $n$ goes to $+\infty$.

**Proposition 5.1.** There exists a solution to (78)–(82), with constant in space potentials in the sense that there exists $\mu_T^\infty = (\mu_{i,T}^\infty)_{1\le i\le I} \in \mathbb{R}^I$ such that

$$ (83) \qquad \log \frac{u_{i,K}^\infty}{u_{0,K}^\infty} + z_i \phi_K^\infty = \mu_{i,T}^\infty, \quad K \in \mathcal{T}, 1 \le i \le I, $$

such that

$$ (84) \qquad (U_T^n, \phi_T^n) \underset{n \to +\infty}{\longrightarrow} (U_T^\infty, \phi_T^\infty). $$

*Proof.* Owing to (42) and the non-negativity of the discrete dissipations $D_T^n$, the sequence $(H_T^n)_{n \ge 0}$ is decreasing, and moreover bounded thanks to Lemma 3.7. Therefore, it converges towards some finite limit $H_T^\infty$, while $D_T^n$ has to tend to 0 as $n$ goes to $+\infty$. Then so does $D_{i,\sigma}^n$, cf. (43), for all $i=1, \dots, I$ and all $\sigma = K|L \in \mathcal{E}_{\text{int}}$. Since $(U_T^n, \phi_T^n)_{n \ge 0}$ is included in a bounded (thus relatively compact) set, the LaSalle invariance principle yields the existence of (at least) an accumulation point $(U_T^\infty, \phi_T^\infty)$, the energy of which being equal to $H_T^\infty$. Then it follows from the discrete energy / energy dissipation inequality (42) that the dissipation $D_T^\infty$ corresponding to $(U_T^\infty, \phi_T^\infty)$ is equal to 0, then so are all the edge dissipation contributions $D_{i,\sigma}^\infty$ which are defined by (43) but with $(U_T^\infty, \phi_T^\infty)$ instead of $(U_T^n, \phi_T^n)$. Since $D_{i,\sigma}^\infty = 0$ iff the corresponding flux $F_{i,K\sigma}^\infty$ defined by (80) vanishes for all $i=1, \dots, I$ and $\sigma = K|L \in \mathcal{E}_{\text{int}}$, one gets that $(U_T^\infty, \phi_T^\infty)$ is a steady solution, and that the convergence (84) holds true.

Let us now show (83). The condition $F_{i,K\sigma}^\infty = 0$ for $\sigma = K|L$ implies the following alternative:

(i) $u_{0,K}^\infty = 0$ and $u_{0,L}^\infty = 0$,

(ii) $u_{i,K}^\infty = 0$ and $u_{i,L}^\infty = 0$,

(iii) $u_{i,K}^\infty = 0$ and $u_{0,K}^\infty = 0$, or $u_{i,L}^\infty = 0$ and $u_{0,L}^\infty = 0$,

(iv) $\log \frac{u_{i,K}^\infty}{u_{0,K}^\infty} + z_i \phi_K^\infty = \log \frac{u_{i,L}^\infty}{u_{0,L}^\infty} + z_i \phi_L^\infty$.

This can be seen using the two formulations of the fluxes (80) and (31) at the limit $n \to \infty$. The formula (80) clearly yields (i) to (iii) if any of the involved volume fractions vanish. If they are all positive then (31) is well-defined at the limit $n \to \infty$ and yields (iv).

Assume for contradiction that there exists $K \in \mathcal{T}$ such that $u_{0,K}^\infty = 0$, then $F_{i,K\sigma}^\infty = 0$ for all $i$ requires that either $u_{i,K}^\infty = 0$ for all $1 \le i \le I$, which is incompatible with (81), or that $u_{0,L}^\infty = 0$ for all neighboring cells $L$ sharing an edge with $K$. An induction shows relying on the connected character of $\Omega$ shows that $u_{0,T}^\infty = 0$, which contradicts (82). The cases (i) and (iii) are therefore impossible, and $u_{0,T}^\infty > 0$.

Assume now that $u_{i,K}^\infty = 0$ for some $K \in \mathcal{T}$, then as $u_{0,K}^\infty > 0$, a similar reasoning shows that $u_{i,L}^\infty$ also has to be equal to 0 for all neighboring cell of $K$. Then $u_{i,T}^\infty = 0$, contradicting again (82). Therefore, alternative (ii) is not feasible, and (iv) has to hold true. This concludes the proof of Proposition 5.1. $\square$

Given $\phi_T^\infty$ and $\mu_T^\infty$, one can invert the system (81)&(83). As in the continuous case, one then gets that

$$ u_{i,K}^\infty = v_i(\phi_K^\infty, \mu_T^\infty), \quad K \in \mathcal{T}, 0 \le i \le I $$

with $v_i$ defined in (18). Incorporating this relation in (78) provides the following discrete counterpart to the modified Poisson-Boltzmann equation (20):

$$ (85) \qquad \lambda^2 \sum_{\sigma \in \mathcal{E}_K} a_\sigma (\phi_K^\infty - \phi_{K\sigma}^\infty) + m_K r (\phi_K^\infty, \mu_T^\infty) = m_K f_K, \quad K \in \mathcal{T}, $$
to be complemented with the mass constraints (82), with $r$ being the non-decreasing function defined in (21).

**Proposition 5.2.** The solution $(\phi_{\mathcal{T}}^{\infty}, \mu_{\mathcal{T}}^{\infty})$ to (85) &(82) minimizes the strictly convex functional $\Psi_{\mathcal{T}}$ defined by (35). In particular, $(\phi_{\mathcal{T}}^{\infty}, \mu_{\mathcal{T}}^{\infty})$ is unique.

*Proof.* One readily checks that

$$ \frac{\partial \Psi_{\mathcal{T}}}{\partial y_K}(y_{\mathcal{T}}, \boldsymbol{\xi}) = \lambda^2 \sum_{\sigma \in \mathcal{E}_K} a_{\sigma}(y_K - y_{K\sigma}) + m_K r(y_K, \boldsymbol{\xi}) - m_K f_K, \quad K \in \mathcal{T}. $$

and that

$$ \frac{\partial \Psi_{\mathcal{T}}}{\partial \xi_i}(y_{\mathcal{T}}, \boldsymbol{\xi}) = \sum_{K \in \mathcal{T}} m_K [v_i(y_K, \boldsymbol{\xi}) - u_{i,K}^0]. $$

Bearing in mind the definition (29) of $u_{i,K}^0$, one gets that $(\phi_{\mathcal{T}}^{\infty}, \mu_{\mathcal{T}}^{\infty})$ solves (85)&(82) if and only if it is a critical point of $\Psi_{\mathcal{T}}$. In particular, Proposition 5.1 guarantees the existence of such a critical point as it ensures the existence of a solution to (85)&(82).

To conclude the proof of Proposition 5.2, it remains to show the strict convexity of $\Psi_{\mathcal{T}}$. Let $(y_{\mathcal{T}}, \boldsymbol{\xi}_{\mathcal{T}})$ and $(\tilde{y}_{\mathcal{T}}, \tilde{\boldsymbol{\xi}}_{\mathcal{T}})$ be two elements of $\mathbb{R}^{\mathcal{T}} \times \mathbb{R}^{\mathcal{I}}$. Then

$$ (86) \qquad \left( D\Psi_{\mathcal{T}}(y_{\mathcal{T}}, \boldsymbol{\xi}_{\mathcal{T}}) - D\Psi_{\mathcal{T}}(\tilde{y}_{\mathcal{T}}, \tilde{\boldsymbol{\xi}}_{\mathcal{T}}) \right) \cdot (y_{\mathcal{T}} - \tilde{y}_{\mathcal{T}}, \boldsymbol{\xi}_{\mathcal{T}} - \tilde{\boldsymbol{\xi}}_{\mathcal{T}}) = A + B, $$

where

$$ A = \lambda^2 \sum_{K \in \mathcal{T}} (y_K - \tilde{y}_K) \sum_{\sigma \in \mathcal{E}_K} a_{\sigma}((y_K - \tilde{y}_K) - (y_{K\sigma} - \tilde{y}_{K\sigma})), $$

$$ B = \sum_{K \in \mathcal{T}} m_K [r(y_K, \boldsymbol{\xi}) - r(\tilde{y}_K, \tilde{\boldsymbol{\xi}})](y_K - \tilde{y}_K) \\ + \sum_{K \in \mathcal{T}} \sum_{i=1}^{I} m_K [(v_i(y_K, \boldsymbol{\xi}) - v_i(\tilde{y}_K, \tilde{\boldsymbol{\xi}}))](\xi_i - \tilde{\xi}_i). $$

Performing a discrete integration by parts gives

$$ A = \lambda^2 \sum_{\sigma \in \mathcal{E}} a_{\sigma}((y_K - \tilde{y}_K) - (y_{K\sigma} - \tilde{y}_{K\sigma}))^2 \geq 0. $$

As we assumed that $\mathcal{E}_{\text{ext}}^{D}$ is not empty, it follows from the discrete Poincaré inequality (see Lemma 9.1 and Remark 9.4 of [26]) that

$$ (87) \qquad A \text{ vanishes if and only if } y_{\mathcal{T}} = \tilde{y}_{\mathcal{T}}. $$

Using the relation (21) between $r$ and $v_i$, we rewrite the term $B$ as

$$ B = \sum_{K \in \mathcal{T}} m_K \sum_{i=1}^{I} (v_i(y_K, \boldsymbol{\xi}) - v_i(\tilde{y}_K, \tilde{\boldsymbol{\xi}}))(\xi_i - z_i y_K - (\tilde{\xi}_i - z_i \tilde{y}_K)). $$

Denote by

$$ H^*(R) = \sup_{U \in \mathcal{A}} \left[ \sum_{i=1}^{I} u_i r_i - H(U) \right] = \log \left( 1 + \sum_{i=1}^{I} e^{r_i} \right), \quad R = (r_i)_{1 \le i \le I} \in \mathbb{R}^{\mathcal{I}}, $$

the Legendre transform of the mixing entropy functional $H$, cf. (10), then, setting

$$ R_K = (r_{i,K})_{1 \le i \le I} = (\xi_i - z_i y_K)_{1 \le i \le I}, \quad K \in \mathcal{T}, $$

the term $B$ is rewritten as

$$ (88) \qquad B = \sum_{K \in \mathcal{T}} m_K \sum_{i=1}^{I} \left\langle \frac{\partial H^*}{\partial r_i}(R_K) - \frac{\partial H^*}{\partial r_i}(\tilde{R}_K), R_K - \tilde{R}_K \right\rangle. $$
As $H^*$ is a Legendre transform, it is by construction convex and thus $B \ge 0$. By setting $r_{0,K} = \tilde{r}_{0,K} = 0$, the convexity of $H^*$ can also be seen as a consequence of Holder's inequality

$$ \sum_{j=0}^{I} e^{\theta r_{j,K} + (1-\theta)\tilde{r}_{j,K}} \le \left( \sum_{j=0}^{I} e^{r_{j,K}} \right)^{\theta} \left( \sum_{j=0}^{I} e^{\tilde{r}_{j,K}} \right)^{1-\theta} $$

for $\theta \in [0, 1]$. The inequality is an equality only if $(e^{r_{j,K}})_j$ and $(e^{\tilde{r}_{j,K}})_j$ are linearly dependent. As $r_{0,K} = \tilde{r}_{0,K}$ it can only happen if $R_K = \tilde{R}_K$. In other words $H^*$ is actually strictly convex.

It follows from the above discussion that the right-hand side of (86) vanishes if and only if both A and B vanish. Thanks to (87), this requires that $y_T = \tilde{y}_T$. On the other hand because of the strict convexity of $H^*$, the term B (88) vanishes if and only if $R_T = \tilde{R}_T$, thus if $\xi = \tilde{\xi}$. Then $\Psi_T$ is strictly convex and its minimizer $(\phi_T^\infty, \mu_T^\infty)$ is unique. $\square$

## 6. NUMERICAL RESULTS

The non-linear system corresponding to the scheme is solved thanks to a Newton-Raphson method with stopping criterion $||F_T^n((U_K^n)_{K \in T}, (\phi_K^n)_{K \in T})||_\infty \le 10^{-10}$, the components of $F_T^n$ being given by the left-hand side of (25).

### 6.1. Convergence under grid refinement.
The goal of our first numerical test is to show our scheme is second order accurate w.r.t. the mesh size. To this end, we consider the one-dimensional domain $\Omega = (0, 1)$, in which $I = 2$ different ions evolve, both with the same diffusion coefficient $D_1 = D_2 = 1$. Their (normalized) charge is set to $z_1 = 2$ and $z_2 = 1$, yielding repulsive interaction. No background charge is considered, i.e. $f = 0$, whereas Dirichlet boundary conditions are imposed for the electric potential on both sides of the interval, that are $\phi^D(t, 0) = 10$ and $\phi^D(t, 1) = 0$. We consider a moderately small Debye length $\lambda^2 = 10^{-2}$. We start at initial time $t = 0$ with the following configurations: $u_1^0(x) = 0.2 + 0.1(x-1)$ and $u_2^0 = 0.4$.

<figure><img src="image_3.png" /><figcaption>FIGURE 2. Left: Concentration profiles $u_1(T, x)$, $u_2(T, x)$ and $u_0(T, x)$ at time $T = 1$ depicted with solid lines. The corresponding long-time limit $(u_1^\infty, u_2^\infty, u_3^\infty)$ is depicted with dashed lines. Right: Electric potential profile $\phi(T, x)$ at time $T = 1$ (solid line) and in the long-time limit (dashed line).</figcaption></figure>

A reference solution is computed on a grid made of 1638400 cells and with a constant time step $\tau = 10^{-3}$, to which are compared solutions computed on successively refined grids but with the same constant time step. The profile of the solution at the final time $T = 1$ is depicted on Figures 2. The relative space-time $L^1$ error is plotted as a function of the number of cells on Figure 3, showing some second order accuracy in space, as specified in the introductory discussion. The number of Newton iterations required to solve a time step remains between 6 for the very first iterations and 2 for larger times and is mainly insensitive to the mesh size.
<figure><img src="image_4.png" /><figcaption>FIGURE 3. Convergence of the schemes under space grid refinement.</figcaption></figure>

**Remark 6.1.** An optimal error estimate showing order 2 in space and order 1 in time convergence can be established on cartesian grids if one assumes that there exists a smooth strong solution to the continuous problem, see for instance [12] in a scalar context, or [21] for classical Poisson-Nernst-Planck system. The extension to unstructured grids is then much more delicate, as usual for finite volumes, see for instance [24] for a contribution in this direction. However, the existence of such a strong solution to the problem is an open question.

6.2. **Long-time behavior of the scheme.** We now focus on the long-time behavior of the scheme, and more specifically on the convergence towards its unique steady state. Our study relies on a two-dimensional test case, with $I = 3$ charged species (plus the solvent $u_0$). The corresponding charges are $z_1 = 2$, $z_2 = 1$, $z_3 = -1$ (and $z_0 = 0$ as the solvent is electrically neutral in our model), while the diffusion coefficients are set to $D_1 = 1$, $D_2 = 2$, and $D_3 = 2$.

The geometry is the following: $\Omega$ is the unit square $(0,1)^2$, the boundary of which being split into $\Gamma^D = \{y=1\} \times \{0 \le x \le 1/2\}$ and its complement $\Gamma^N$. The electric potential is set to 0 on $\Gamma^D$, i.e. $\phi^D = 0$.

<figure><img src="image_5.png" /><figcaption>FIGURE 4. Geometry of the two-dimensional domain</figcaption></figure>

We run our scheme on a Delaunay mesh made of 7374 triangles and 11167 edges generated with Gmsh for 3 different initial concentration profiles. Two initial profiles $U^{0,(1)}$ and $U^{0,(2)}$ are chosen globally neutral:

$$ (89) \qquad \begin{cases} u_1^{0,(1)}(x) = 0.3 \times 1_{(0,1/2)^2}(x), \\ u_2^{0,(1)}(x) = 0.3 \times 1_{(1/2,1) \times (0,1/2)}(x), \\ u_3^{0,(1)}(x) = 0.9 \times 1_{(1/2,1)^2}(x). \end{cases} $$
and

$$ (90) \qquad \left\{ \begin{array}{l} u_1^{0,(2)}(x) = 0.1 \times u_1^{0,(1)}(x), \\ u_2^{0,(2)}(x) = 0.1 \times u_2^{0,(1)}(x) + 0.9 \times \mathbb{1}_{(0,1)\times(1/2,1)}(x), \\ u_3^{0,(2)}(x) = 0.1 \times u_3^{0,(1)}(x) + 0.9 \times \mathbb{1}_{(0,1)\times(0,1/2)}(x). \end{array} \right. $$

A third initial configuration is chosen globally charged and constant in space:

$$ (91) \qquad u_1^{0,(3)}(x) = 0.2, \quad u_2^{0,(3)}(x) = 0.2, \quad u_3^{0,(3)}(x) = 0.3. $$

The steady state corresponding to (91) is not constant w.r.t. space, as shown on Figure 5.

<figure><img src="image_6.png" /><figcaption>FIGURE 5. Stationary concentration profiles corresponding to the initialization (91) for $\lambda^2 = 0.01$</figcaption></figure>

We run our scheme with a constant time step $\tau = 10^{-4}$ and for two different Debye lengths until a final time $T = 3$, and look for the evolution of the relative energy $H_T^n - H_T^\infty$ along-time. The relative energy is decaying for all the curves plotted on Figure 6, but the velocity at which the decay occurs varies strongly depending on the Debye length and on the initial profile. In the electroneutral configurations, we mainly observe that higher quantities of $u_0$ yield a faster convergence. This is expected as the dissipation (14) contains a $u_0^2$ prefactor, in opposition to a mere $u_0$ for the classical Poisson-Nernst-Planck system. Then, still in the electroneutral regime, the smaller is $\lambda^2$ then the faster is the convergence towards the constant in space equilibrium. The trend is different for globally charged profiles. For the choice (91) of the initial profile, the convergence becomes extremely slow when $\lambda^2$ becomes small. This is due to the fact that, in some areas, the concentration of $u_0$ becomes extremely small. The equilibration between the other phases is then extremely slow in these regions. To illustrate this difficulty, we depict on Figure 7 the concentration profiles at $T = 3$. The solvent concentration $u_0$ is already in good agreement with the equilibrium profile depicted on Figure 5 (see also Figure 8 for a representation of the difference between $U_{T,\tau}(T)$ and $U_T^\infty$), but this is not the case of the other concentrations, the
<figure><img src="image_7.png" /><figcaption>FIGURE 6. Convergence towards the steady long-time behavior in terms of relative energy $H_T^n - H_T^\infty$.</figcaption></figure>

evolution of which being stuck in the zones where the concentration of $u_0$ is small. Obtaining a fine quantitative understanding of the long-time behavior of the model should be investigated in future works.

<figure><img src="image_8.png" /><figcaption>FIGURE 7. Concentration profiles at final time $T = 3$ corresponding to the initialization (91) for $\lambda^2 = 0.01$.</figcaption></figure>

**Acknowledgement.** CC and MH acknowledge support from the LabEx CEMPI (ANR-11-LABX-0007) and the ministries of Europe and Foreign Affairs (MEAE) and Higher Education, Research and Innovation (MESRI) through PHC Amadeus 46397PA. AM acknowledges support from the
<figure><img src="image_9.png" /><figcaption>FIGURE 8. Difference between the final concentration profiles of Figure 7 and the stationary ones, cf. Figure 5(91) for $\lambda^2 = 0.01$.</figcaption></figure>

multilateral project of the Austrian Agency for International Co-operation in Education and Research (OeAD), grant FR 01/2021, as well as partial support from the Austrian Science Fund (FWF), grant DOI 10.55776/P33010 and 10.55776/F65. This work has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme, ERC Advanced Grant NEUROMORPH, no. 101018153. CC also acknowledges support from the European Union's Horizon 2020 research and innovation programme under grant agreement No 847593 (EURAD program, WP DONUT) and from the COMODO (ANR-19-CE46-0002) project.

## REFERENCES

[1] B. Andreianov, C. Cancès, and A. Moussa. A non-linear time compactness result and applications to discretization of degenerate parabolic-elliptic PDEs. *J. Funct. Anal.*, 273:3633-3670, 2017.

[2] R. Bailo, J. A. Carrillo, and J. Hu. Bound-preserving finite-volume schemes for systems of continuity equations with saturation. *SIAM J. Appl. Math.*, 83(3):1315-1339, 2023.

[3] M. Bessemoulin-Chatard and C. Chainais-Hillairet. Exponential decay of a finite volume scheme to the thermal equilibrium for drift-diffusion systems. *J. Numer. Math.*, 25(3):147-168, September 2017.

[4] M. Bessemoulin-Chatard, C. Chainais-Hillairet, and M.-H. Vignal. Study of a finite volume scheme for the drift-diffusion system. Asymptotic behavior in the quasi-neutral limit. *SIAM J. Numer. Anal.*, 52(4):1666-1691, 2014.

[5] M. Burger, M. Di Francesco, J.-F. Pietschmann, and B. Schlake. Nonlinear cross-diffusion with size-exclusion. *SIAM J. Math. Anal.*, 46(6):2842-2871, 2010.

[6] M. Burger, B. Schlake, and M.-T. Wolfram. Nonlinear Poisson-Nernst-Planck equations for ion flux through confined geometries. *Nonlinearity*, 25(4):961, 2012.

[7] C. Cancès, C. Chainais-Hillairet, A. Gerstenmayer, and A. Jüngel. Finite-volume scheme for a degenerate cross-diffusion model motivated from ion transport. *Numer. Methods Partial Differential Equations*, 35(2):545-575, 2019.
[8] C. Cancès, C. Chainais-Hillairet, B. Merlet, F. Raimondi, and J. Venel. Mathematical analysis of a thermo-dynamically consistent reduced model for iron corrosion. *Z. Angew. Math. Phys.*, 74, 2023. Article number 96.

[9] C. Cancès, V. Ehrlacher, and L. Monasse. Finite Volumes for the Stefan-Maxwell cross-diffusion system. *IMA J. Numer. Anal.*, 44(2):1029-1060, 2024.

[10] C. Cancès and B. Gaudeul. A convergent entropy diminishing finite volume scheme for a cross-diffusion system. *SIAM J. Numer. Anal.*, 58(5):2684-2710, 2020.

[11] C. Cancès, M. Herda, and A. Massimini. Finite volumes for a generalized Poisson-Nernst-Planck system with cross-diffusion and size exclusion. In E. Franck, J. Fuhrmann, V. Michel-Dansac, and L. Navoret, editors, *Finite Volumes for Complex Applications X—Volume 1, Elliptic and Parabolic Problems. FVCA 2023*, volume 432 of Springer Proceedings in Mathematics & Statistics, pages 57-73, Cham, 2023. Springer.

[12] C. Cancès and J. Venel. On the square-root approximation finite volume scheme for nonlinear drift-diffusion equations. *Comptes Rendus. Mathématique*, 361:525-558, 2023.

[13] C. Cancès and A. Zurek. A convergent finite volume scheme for dissipation driven models with volume filling constraint. *Numer. Math.*, 151:279-328, 2022.

[14] C. Cancès, C. Chainais-Hillairet, J. Fuhrmann, and B. Gaudeul. A numerical analysis focused comparison of several finite volume schemes for a unipolar degenerated drift-diffusion model. *IMA J. Numer. Anal.*, 41(1):271-314, 2021.

[15] J. A. Carrillo, F. Filbet, and M. Schmidtchen. Convergence of a finite volume scheme for a system of interacting species with cross-diffusion. *Numer. Math.*, 145:473-511, 2020.

[16] C. Chainais-Hillairet. Entropy method and asymptotic behaviours of finite volume schemes. In *Finite volumes for complex applications. VII. Methods and theoretical aspects*, volume 77 of Springer Proc. Math. Stat., pages 17-35. Springer, Cham, 2014.

[17] C. Chainais-Hillairet, J.-G. Liu, and Y.-J. Peng. Finite volume scheme for multi-dimensional drift-diffusion equations and convergence analysis. *ESAIM: M2AN*, 37(2):319-338, 2003.

[18] M. Chatard. Asymptotic behavior of the Scharfetter-Gummel scheme for the drift-diffusion model. In *Finite volumes for complex applications VI. Problems & perspectives. Volume 1, 2*, volume 4 of Springer Proc. Math., pages 235-243. Springer, Heidelberg, 2011.

[19] X. Chen, A. Jüngel, X. Lin, and L. Liu. Large-time asymptotics for degenerate cross-diffusion population models with volume filling. *J. Differential Equations*, 386:1-15, 2024.

[20] Y. Coudière, J.-P. Vila, and P. Villedieu. Convergence rate of a finite volume scheme for a two dimensional convection-diffusion problem. *ESAIM Mathematical Modelling and Numerical Analysis*, 33:493-516, 1999.

[21] Ding, J., Wang, C., and Zhou, S. Convergence analysis of structure-preserving numerical methods based on Slotboom transformation for the Poisson-Nernst-Planck equations. *Comm. Math. Sci.*, 21(2):459-484, 2023.

[22] J. Droniou and R. Eymard. Study of the mixed finite volume method for Stokes and Navier-Stokes equations. *Numerical Methods for Partial Differential Equations*, 25:137-171, 2009.

[23] J. Droniou and R. Eymard. The asymmetric gradient discretisation method. In C. Cancès and P. Omnes, editors, *Finite volumes for complex applications VIII - methods and theoretical aspects*, volume 199 of Springer Proc. Math. Stat., pages 311-319, Cham, 2017. Springer.

[24] J. Droniou and N. Nataraj. Improved $L^2$ estimate for gradient schemes and super-convergence of the TPFA finite volume scheme. *IMA J. Numer. Anal.*, 38(3):1254-1293, 2018.

[25] R. Eymard and T. Gallouët. H-convergence and numerical schemes for elliptic problems. *SIAM J. Numer. Anal.*, 41(2):539-562, 2003.

[26] R. Eymard, T. Gallouët, and R. Herbin. Finite volume methods. In *Handbook of numerical analysis, Vol. VII*, pages 713-1020. North-Holland, Amsterdam, 2000.

[27] H. Gajewski. On a variant of monotonicity and its application to differential equations. *Nonlinear Anal.*, 22(1):73-80, 1994.

[28] H. Gajewski and I. V. Skrypnik. To the uniqueness problem for nonlinear parabolic equations. *Discr. Cont. Dyn. Syst.*, 10(1&2):315-336, 2003.

[29] T. Gallouët and J.-C. Latché. Compactness of discrete approximate solutions to parabolic PDEs—application to a turbulence model. *Commun. Pure Appl. Anal.*, 11(6):2371-2391, 2012.

[30] B. Gaudeul and J. Fuhrmann. Entropy and convergence analysis for two finite volume schemes for a Nernst-Planck-Poisson system with ion volume constraints. *Numer. Math.*, 151(1):99-149, 2022.

[31] A. Gerstenmayer and A. Jüngel. Analysis of a degenerate parabolic cross-diffusion system for ion transport. *J. Math. Anal. Appl.*, 461(1):523-543, 2018.

[32] A. Gerstenmayer and A. Jüngel. Comparison of a finite-element and finite-volume scheme for a degenerate cross-diffusion system for ion transport. *Comput. Appl. Math.*, 38(3):Art. 108, 23, 2019.

[33] M. Heida. Convergences of the squareroot approximation scheme to the Fokker-Planck operator. *Math. Models Methods Appl. Sci.*, 28(13):2599-2635, 2018.

[34] M. Heida, M. Kantner, and A. Stephan. Consistency and convergence of a family of finite volume discretizations of the Fokker-Planck operator. *ESAIM: Math. Model Numer. Anal.*, 55(6):3017-3042, 2021.

[35] A. Jüngel. The boundedness-by-entropy method for cross-diffusion systems. *Nonlinearity*, 28(6):1963, 2015.

[36] A. Jüngel and A. Zurek. A convergent structure-preserving finite-volume scheme for the Shigesada-Kawasaki-Teramoto population system. *SIAM J. Numer. Anal.*, 59(4):2286-2309, 2021.
[37] A. Jüngel and A. Zurek. A discrete boundedness-by-entropy method for finite-volume approximations of cross-diffusion systems. *IMA J. of Numer. Anal.*, 43(1):560–589, 2023.

[38] H. C. Lie, K. Fackeldey, and M. Weber. A square root approximation of transition rates for a Markov state model. *SIAM J. Matrix Anal. Appl.*, 34:738–756, 2013.

[39] D. L. Scharfetter and H. K. Gummel. Large-signal analysis of a silicon read diode oscillator. *Electron Devices, IEEE Transactions on*, 16(1):64–77, 1969.

[40] A. Schlichting and C. Seis. The Scharfetter–Gummel scheme for aggregation–diffusion equations. *IMA J. Numer. Anal.*, 42(3):2361–2402, 2022.

UNIV. LILLE, CNRS, INRIA, UMR 8524 - LABORATOIRE PAUL PAINLEVÉ, F-59000 LILLE, FRANCE
*Email address: clement.cances@inria.fr, maxime.herda@inria.fr*

CERMICS, ECOLE DES PONTS PARIS-TECH AND INRIA MATHERIALS TEAM-PROJECT, 6 ET 8 AV. BLAISE PASCAL CITÉ DESCARTES 77455 MARNE-LA-VALLÉE, FRANCE
*Email address: annamaria.massimini@enpc.fr*