# Finite volumes for a generalized Poisson-Nernst-Planck system with cross-diffusion and size exclusion

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

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

## To cite this version:

Clément Cancès, Maxime Herda, Annamaria Massimini. Finite volumes for a generalized Poisson-Nernst-Planck system with cross-diffusion and size exclusion. Finite Volumes for Complex Applications X, Oct 2023, Strasbourg, France. pp.57-73, 10.1007/978-3-031-40864-9_4. hal-04022357v2

HAL Id: hal-04022357

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

Submitted on 8 Nov 2023

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Finite volumes for a generalized Poisson-Nernst-Planck
system with cross-diffusion and size exclusion

Clément Cancès<fnref n="1" />, Maxime Herda<fnref n="1" />, and Annamaria Massimini<fnref n="2" />

**Abstract.** We present two finite volume approaches 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. Both methods utilize a two-point flux approximation and are part of the exponentially fitted scheme framework. The only difference between the two is the selection of a Stolarsky mean for the drift term originating from a self-consistent electric potential. The first version of the scheme, referred to as (SQRA), uses a geometric mean and is an extension of the squareroot approximation scheme. The second scheme, (SG), utilizes an inverse logarithmic mean to create a generalized version of the Scharfetter-Gummel scheme. Both approaches ensure the decay of some discrete free energy. Classical numerical analysis results – existence of discrete solution, convergence of the scheme as the grid size and the time step go to 0 – follow. Numerical simulations show that both schemes are effective for moderately small Debye lengths, with the (SG) scheme demonstrating greater robustness in the small Debye length regime.

**Keywords:** Drift-diffusion, cross-diffusion, exponential fitting, free energy decay, convergence

# 1 The continuous generalized Poisson-Nernst-Planck model

Motivated by the transfer of ions in confined geometries, Burger *et al.* introduced in [3] a model accounting for cross-diffusion and size-exclusion effects. 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* writes

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

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

$$ 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). \quad (2) $$

<footnote n="1">¹ Inria, Univ. Lille, CNRS, UMR 8524 - Laboratoire Paul Painlevé, F-59000 Lille, France
clement.cances@inria.fr, maxime.herda@inria.fr</footnote>

<footnote n="2">² Technische Universität Wien, Institute of Analysis and Scientific Computing, Wiedner
Hauptstraße 8–10, 1040 Wien, Austria
annamaria.massimini@asc.tuwien.ac.at</footnote>
In the above expression, $D_i > 0$ denotes the diffusion coefficient of the species $i$. The quantity

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

shall be thought 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 nonnegative, leading to size exclusion for the other species $u_i$, $i = 1, \dots, I$. Denoting by $z_i$ the charge of species $i$ and by $\lambda > 0$ the (scaled) Debye length, then the electrostatic potential solves the Poisson equation

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

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 be split into an insulator part $\Gamma^N$ and its complement $\Gamma^D$ on which Dirichlet boundary condition is imposed:

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

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), like for instance Robin type boundary condition modeling electrochemical reaction thanks to Butler-Volmer type formula, see for instance [5], or boundary conditions of mixed Dirichlet-Neumann type as in [10]. In the presentation of the scheme, we assume for simplicity that the system is isolated, in the sense that

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

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

$$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. \qquad (7)$$

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

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

Multiplying (1) by $\log w_i = \log \frac{u_i}{u_0} + z_i \phi$, integrating over $\Omega$ and summing over $i = 1, \dots, I$ yields

$$\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, \qquad (9)$$

where, denoting the mixing (neg)entropy density function $H: \mathbb{R}_+^{I+1} \to \mathbb{R}_+^{I+1}$ by

$$H(U) = u_0 \log(u_0) + \sum_{i=1}^{I} u_i \log(u_i) + \log(I+1) \ge 0,$$
the free energy $\mathcal{H}$ is given by

$$ \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$.

A finite volume scheme has been studied in [4]. Even though the scheme mainly behaves well in practice, its mathematical study is very partial since requiring strong assumptions such as constant diffusion coefficients $D_i = D$ for all $i$, or no charge $z_i = 0$. Moreover, since the scheme proposed in [4] uses upwinding for the mobilities, numerical experiments exhibit a mere first order convergence in space. An alternative finite element method using the so-called electrochemical potentials $\mu_i = \log(w_i)$ rather than the $u_i$ as primary variables has been analyzed in [11]. This latter scheme is by construction free energy diminishing without further restriction on the physical parameters, and is shown to converge towards a weak solution as the mesh size and the time step tend to 0 (up to quadrature error terms). Second order convergence w.r.t. the mesh size is observed, but the nonlinear system to be solved at each time step is stiffer than for the finite volume scheme because of the use of the electrochemical potentials as variables, so that no clear gain was observed in comparison with the upstream mobility finite volumes. The finite volume scheme proposed in [1], in which the fluxes $F_i$ are approximated thanks to the second expression of (2) also leads to singular numerical fluxes expressions. Our goal here is to propose and to analyze a scheme which shares the best with the aforementioned approaches: decay of the free energy and unconditional convergence are established, second order accuracy in space and well-behaved nonlinear system for moderately small Debye length.

## 2 Two finite volume schemes

First, we introduce the time discretization and the spatial mesh of the domain $\Omega$. The mesh will be assumed to be admissible in the sense of [9], in the sense that it fulfills the so-called orthogonality condition, which is usual for two-point flux approximation finite volumes.

Let $\mathcal{T}$ denote a family of non-empty, disjointed, 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 in the sense that $\bar{\Omega} = \bigcup_{K \in \mathcal{T}} \bar{K}$. Further, we call $\mathcal{E}$ a family of edges/faces, where $\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} \text{ s. t. } \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, it follows that $\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 $\overline{x_K x_L}$ is orthogonal to $\sigma = K|L$. The triplet $(\mathcal{T}, \mathcal{E}, \{x_K\}_{K \in \mathcal{T}})$ is called an admissible mesh.

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 = \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 = \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}_+ := [0, +\infty)$ into a sequence of increasing number of time steps $0 = t^0 < t^1 < \dots$ with a stepsize

$$ \tau^n = t^n - t^{n-1} $$

at time step $n \in \mathbb{N} \setminus \{0\}$. We finally introduce $\Delta t = \sup_{n \in \mathbb{N} \setminus \{0\}} \tau^n$, which we assume to be finite.

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

$$ \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}, \quad (10) $$

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

$$ \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, m, \quad K \in \mathcal{T}. \quad (11) $$

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

$$ 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))), \quad (12) $$
with

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

Formula (12) 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 continuous system (1)–(2) was originally derived in [3] thanks to a hopping process, suggesting the choice

$$\mathfrak{B}(y) = e^{-y/2}, \qquad (\text{SQRA})$$

leading to a scheme referred to as the square-root approximation (SQRA) scheme in what follows, in reference to [14, 12, 6]. Another natural choice for the function $\mathfrak{B}$ is the Bernoulli function

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

the corresponding scheme being referred to as the Scharfetter-Gummel (SG) scheme although its construction is not based on the original idea of [15]. We rather take advantage of the free-energy diminishing character of the SG scheme highlighted in [8].

In order to close the system, it remains to define the discrete counterpart to $u_0$ as follows:

$$u_{i,K}^0 = \frac{1}{m_K} \int_K u_i^0, \quad K \in \mathcal{T}, i = 0, \dots, I. \qquad (14)$$

Then we infer from (7) that

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

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

The consistency of the discrete fluxes (12) with the continuous ones (2) might not look completely obvious. For the particular choice (SQRA) of the function $\mathfrak{B}$, such a consistency proof is given in [6], with second order accuracy in space. Readers can also convince themselves that the reformulation (22) of the discrete fluxes is consistent with the reformulation (8) of the continuous fluxes.

### 3 Stability and convergence properties of the schemes

The goal of this section is to show that the nonlinear system corresponding to the scheme (10)–(13) 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 and the decay of the free energy. The grid $\mathcal{T}$ and the time steps $(\tau^n)_{n \ge 1}$ remain fixed.

Since our scheme is locally conservative, i.e., $F_{K\sigma}^n + F_{L\sigma}^n = 0$ for all $\sigma = K|L \in \mathcal{E}_{\text{int}}$, then summing (11) over $K$ shows by induction and thanks to (14) that

$$\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. \qquad (16)$$
Since we are interested in discrete solution with positive volume fractions $u_{i,K}^n$, we perform an eventually harmless modification of the flux formula (12) into

$$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). \quad (17)$$

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

$$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. \quad (18)$$

Then any solution $(U_K^n, \phi_K^n)_{K \in \mathcal{T}, n \ge 1}$ to the modified scheme with (17) instead of (12) 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 (17) that $F_{i,K\sigma}^n \ge 0$ for all $\sigma \in \mathcal{E}_K$ and all $i = 1, \dots, I$. Because of (13) and (18), 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 (17) 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 (13), 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 (16), 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 (17), 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 (17) that $u_{i,L}^n \le 0$ for all cell $L$ sharing a cell $\sigma = K|L$ with $K$. As above, this implies as $u_{0,K}^n = 0$ for all $K \in \mathcal{T}$, which contradicts (16). Then $u_{i,K}^n > 0$ for all $K \in \mathcal{T}$, concluding the proof of Proposition 1.

A consequence of previous proposition is that a solution to the modified scheme with (17) instead of (12) is also a solution to the original scheme (10)–(13). 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)$. Therefrom, we deduce some uniform discrete $L^\infty(H^1(\Omega))$ estimate on $(\phi_K)_{K \in \mathcal{T}}$ from [9, Lemma 9.4], while [7, Proposition A.1] gives a uniform bound

$$|\phi_K^n| \le C, \quad K \in \mathcal{T}, n \ge 0, \quad (19)$$
since the right-hand side of the discrete Poisson equation (10) is uniformly bounded. These a priori estimates are sufficient to prove the existence of a solution to the scheme thanks to a topological degree argument we do not detail here. We end up with the following proposition. The proof, which will be detailed in a forthcoming contribution, is an easy generalization to the one given in [6].

**Proposition 2** There exists at least one solution to the numerical scheme (10)–(13) such that $u_{i,K}^n > 0$ for all $i = 0, \dots, I$, for all $K \in \mathcal{T}$ and all $n \ge 1$.

Next proposition is about the thermodynamical consistency of our scheme and the decay of a discrete counterpart of the free energy.

**Proposition 3** Let $(U_K^n, \phi_K^n)_{K \in \mathcal{T}, n \ge 1}$ be a solution to the scheme (10)–(13) as in Proposition 2, then define for $n \ge 0$ the discrete free energy at the $n^{th}$ time step

$$ \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_K^n - \phi_{\sigma}^D), \quad (20) $$

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 is nonnegative for both choices (SQRA) and (SG) of function $\mathfrak{B}$. Then there holds

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

*Proof.* With both choices (SQRA) and (SG) for the function $\mathfrak{B}$, the fluxes (12) enter the framework of the exponentially fitted schemes. Indeed, denoting by

$$ w_{i,K}^n = \frac{u_{i,K}^n}{u_{0,K}^n} e^{z_i \phi_K^n} = \exp(\mu_{i,K}^n) \quad \text{for } K \in \mathcal{T} \text{ and } i = 1, \dots, I $$

(which is well defined since $u_{0,K}^n > 0$), then the fluxes (12) can be reformulated as

$$ 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) \quad (22) $$

for some mean function $\mathfrak{M}$ depending on the choice of $\mathfrak{B}$ (see [13]). More precisely,

$$ \mathfrak{M}(a, b) = \sqrt{ab} \quad \text{for (SQRA), and} \quad \mathfrak{M}(a, b) = \frac{\log(1/a) - \log(1/b)}{1/a - 1/b} \quad \text{for (SG),} $$

for $a, b > 0$ with $a \ne b$, and $\mathfrak{M}(a, a) = a$. As a consequence of the positivity of $u_{0,K}^n$ and of the monotonicity of the exponential function, one easily infers that

$$ \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}}, \quad (23) $$
whence the nonnegativity of $\mathcal{D}^n$.

Define by $\mu_{i,K}^n = \log(\frac{u_{i,K}^n}{u_{0,K}^n}) + z_i \phi_K^n = \log(w_{i,K}^n)$ the electrochemical potential of species $i$, then multiplying the discrete conservation law (11) by $\tau^n \mu_{i,K}^n$, and summing over $i = 1, \dots, I$ and $K \in \mathcal{T}$ provides thanks to discrete integration by parts

$$ \mathcal{A}_{\mathcal{T}}^n + \mathcal{B}_{\mathcal{T}}^n + \tau^n \mathcal{D}_{\mathcal{T}}^n = 0, \quad (24) $$

where we have set

$$ \begin{aligned} \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{(13)}{=} \sum_{i=0}^{I} \sum_{K \in \mathcal{T}} (u_{i,K}^n - u_{i,K}^{n-1}) \log (u_{i,K}^n) m_K, \end{aligned} $$

and

$$ \begin{aligned} \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{(10)}{=} \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{aligned} $$

Then we deduce from the convexity of $H$ that

$$ \mathcal{A}_{\mathcal{T}}^n \geq \sum_{K \in \mathcal{T}} (H(U_K^n) - H(U_K^{n-1})) m_K, \quad (25) $$

while reorganizing the term $\mathcal{B}^n$ gives

$$ \begin{aligned} \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) \\ &+ \lambda^2 \sum_{\sigma \in \mathcal{E}^D} a_{\sigma} \phi_{\sigma}^D (\phi_K^n - \phi_{K\sigma}^{n-1}). \end{aligned} $$

Then using the elementary convexity inequality $a(a-b) \geq (a^2-b^2)/2$ in the above term and combining the result with (25) in (24) provides the desired result (21).

Proposition 3 is interesting in itself, but it also contains important information for proving the convergence of the scheme, as in particular the discrete $L^2_{\text{loc}}(H^1)$ estimates on the discrete counterparts of $u_0$ and $\sqrt{u_i u_0}$. We prove these estimates in Lemma 2. As an intermediate result we need a uniform bound on the discrete free energy.

**Lemma 1** There exists $C > 0$ depending only on $\Omega$, $\phi^D$, $\lambda$, $f$, $(z_i)_i$, and $\zeta$ such that, for all $N \geq 1$, there holds $|\mathcal{H}_{\mathcal{T}}^N| \leq C$.

*Proof.* Because of the bound $0 \leq u_{i,K}^n \leq 1$ for all $i$ and $K$, it is clear that the first two contributions of (20) remain uniformly bounded. 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

$$ \sum_{\sigma \in \mathcal{E}^D} a_{\sigma} \phi_{\sigma}^{D} (\phi_{K}^{n} - \phi_{\sigma}^{D}) = \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}) + \sum_{K \in \mathcal{T}} \phi_{K}^{D} \sum_{\sigma \in \mathcal{E}^D} a_{\sigma} (\phi_{K}^{n} - \phi_{K\sigma}^{n}) $$

which yields using Young's inequality for the first term and the Poisson equation for the second one that

$$ \left| \lambda^2 \sum_{\sigma \in \mathcal{E}^D} a_{\sigma} \phi_{\sigma}^{D} (\phi_{K}^{n} - \phi_{\sigma}^{D}) \right| \leq C (\|\nabla \phi^D\|_{L^2}^2 + \|\phi^D\|_{L^1} (1 + \|f\|_{L^{\infty}})) + \frac{\lambda^2}{4} \sum_{\sigma \in \mathcal{E}} a_{\sigma} (\phi_{K}^{n} - \phi_{K\sigma}^{n})^2 $$

for some $C$ depending only on the domain, $\lambda$, $\zeta$ and $(z_i)_i$.

**Lemma 2** There exists $C > 0$ depending only on $\Omega$, $\phi^D$, $\lambda$, $f$, $(z_i)_i$, $(D_i)_i$ and $\zeta$ such that, for all $N \ge 1$, there holds

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

*Proof.* One gets from the elementary inequality $(a-b)(\log(a) - \log(b)) \ge 4(\sqrt{a} - \sqrt{b})^2$ applied to (23) 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_L^n - \phi_K^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)}$ being equal to 1 for the choice (SQRA) of $\mathfrak{B}$ but not for (SG). However, thanks to (19) 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 \ge \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 - \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 (19), one has $\sinh^2(\frac{z_i}{4}(\phi_K^n - \phi_L^n)) \le C(\phi_K^n - \phi_L^n)^2$. Using more- over 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

$$ \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 \in \mathcal{E}_{\text{int}}$ and using (13) leads to

$$ \begin{aligned} \mathcal{D}_{\mathcal{F}}^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^n)^2. \end{aligned} $$

Invoked again the arguments developed in the discussion preceding Proposition 2 to get a uniform discrete $L^\infty(H^1)$ estimate on $(\phi_K^n)_{K,n}$, we obtain that

$$ \mathcal{D}_{\mathcal{F}}^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. \quad (26) $$

Moreover, the inequality $\sum_{i=0}^{I} \sqrt{u_{i,K}^n u_{i,L}^n} \le 1$ gives that

$$ \begin{aligned} & \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} \\ & \geq \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}, \end{aligned} $$

whence we also deduce that

$$ \mathcal{D}_{\mathcal{F}}^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 conclude the proof, it eventually remains to remark from (21) and Lemma 1 that there exists $C$ depending neither on $h$, $\Delta t$, $N$ nor on the initial data $U^0 = (u_i^0)_{0 \le i \le I}$ (provided it fulfills (7)) such that $\sum_{n=1}^N \tau^n \mathcal{D}_{\mathcal{F}}^n \le C$. Combining this with (26) yields the desired result.

One also deduces the following discrete $L^2_{\text{loc}}(L^2)^d$ estimates on the fluxes, which amount to some discrete $L^2_{\text{loc}}(H^1)'$ estimate on time increments of the discrete counterpart to $\partial_t u_i$.
**Lemma 3** There exists C depending only on $\Omega$, $\phi^D$, $\lambda$, $f$, $(z_i)_i$, $(D_i)_i$ and $\zeta$ such that

$$ \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 + \sum_{n=1}^{N} \tau^n). \quad (27) $$

*Proof.* One splits the flux (2) 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

$$ 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))]. $$

The flux $(F_{i,K\sigma}^n)_{\sigma,n}$ is bounded in $L_{\text{loc}}^2(L^2)^d$ in the sense of (27) if both $(F_{i,K\sigma}^{\text{conv},n})_{\sigma,n}$ and $(F_{i,K\sigma}^{\text{diff},n})_{\sigma,n}$ are. For the choice (SG) of the function $\mathfrak{B}$, then $\mathfrak{B}(-y) - \mathfrak{B}(y) = y$, while $\mathfrak{B}(-y) - \mathfrak{B}(y) = y + \mathcal{O}(y^2)$ for (SQRA), so that

$$ 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) + \mathcal{O}(a_{\sigma}(\phi_K^n - \phi_L^n)^2), \quad (28) $$

the remainder term being null for (SG). 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}$ inherited from the control of the energy $\mathcal{H}_{\mathcal{F}}^n$, to be combined with (19) to control the remainder term.

Concerning the diffusive term, one has for both choices (SQRA) and (SG) of the function $\mathfrak{B}$ that

$$ 1 \le \frac{1}{2} [\mathfrak{B}(z_i(\phi_L^n - \phi_K^n)) + \mathfrak{B}(z_i(\phi_K^n - \phi_L^n))] \le 1 + \mathcal{O}((\phi_K^n - \phi_L^n)^2). $$

Therefore, one gets that

$$ 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) (1 + \mathcal{O}((\phi_K^n - \phi_L^n)^2)). \quad (29) $$

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)$, Lemma 2 provides the desired $L_{\text{loc}}^2(L^2)$ bound on $F_{i,K\sigma}^{\text{diff},n}$, hence Lemma 3.

The above estimates are sufficient to establish the convergence of the numerical scheme. For a given mesh $\mathcal{T}$ and a given time discretization $\tau = (\tau^n)_{n \ge 1}$, we denote by $u_{i,\mathcal{T},\tau}$ and $\phi_{\mathcal{T},\tau}$ the piecewise constant reconstructions defined by

$$ 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]. $$
**Theorem 1** Let $(\mathcal{T}_\ell)_{\ell \ge 1}$ be a sequence of admissible discretizations of $\Omega$ (satisfying the orthogonality condition), such that $h_\ell$ goes to 0 as $\ell$ tends to $+\infty$ while the mesh regularity factor $\zeta_\ell$ remains bounded uniformly w.r.t. $\ell$, and let $(\tau_\ell)_{\ell \ge 1} = ((\tau_\ell^n)_{n \ge 1})_{\ell \ge 1}$ be a sequence of sequences of time steps such that $\Delta t_\ell = \sup_n \tau_\ell^n$ goes to 0 as $\ell$ tends to $+\infty$. Then there exists a weak solution $(U, \phi)$ such that, up to a subsequence,

$$ \phi_{\mathcal{T},\tau} \underset{h, \Delta t \to 0}{\longrightarrow} \phi \quad \text{in } L^\infty(\mathbb{R}_+ \times \Omega)\text{-weak-\star sense and a.e. in } \mathbb{R}_+ \times \Omega, \quad (30) $$

$$ u_{i, \mathcal{T}, \tau} \underset{h, \Delta t \to 0}{\longrightarrow} u_i \quad \text{in the } L^\infty(\mathbb{R}_+ \times \Omega)\text{-weak-\star sense,} \quad i = 0, \dots, I, \quad (31) $$

with furthermore $u_{0, \mathcal{T}, \tau}$ and $u_{i, \mathcal{T}, \tau}(u_{0, \mathcal{T}, \tau})^{1/2}$ converging a.e. in $\mathbb{R}_+ \times \Omega$ towards their respective limits $u_0$ and $u_i(u_0)^{1/2}$ which belong to $L_{\text{loc}}^2(H^1)$.

The proof is technical and will be detailed in a forthcoming contribution. It borrows ideas to the proof proposed in [4] and relies on compactness arguments (in particular on the degenerate Aubin-Lions lemma [4, Lemma 10]) as well as on a suitable notion of weak solution. Indeed, yet another reformulation of the fluxes is needed, like for instance

$$ 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). $$

This last formulation is suitable to establish the convergence since it clearly belongs to $L_{\text{loc}}^2(L^2)^d$ as the product of gradient terms the approximation of which being weakly convergent in $L_{\text{loc}}^2(L^2)^d$ with bounded zeroth order term the approximation of which being strongly convergent.

# 4 Numerical results

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

The goal of our first numerical test is to show that both schemes corresponding to (SQRA) and (SG) are 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$.

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 times $T = 1$ and $T = 5000$ is depicted on Figures 1, 2 and 3. The relative space-time $L^1$ error is plotted as a function of the number of cells on Figure 4, showing some second order accuracy in space, as specified in the introductory discussion. For such a moderately
<figure><img src="image_3.png" /><figcaption>Fig. 1: Concentration profiles $u_1(T,x)$, $u_2(T,x)$ and $u_0(T,x)$ at time $T = 1$, $\lambda^2 = 10^{-2}$.</figcaption></figure>

small value of $\lambda^2 = 10^{-2}$, both schemes exhibit a very similar behavior in terms of accuracy, but also in terms of nonlinear resolution. More precisely, the number of Newton iterations required to solve a time step remains between 6 for the very first iterations and 2 for larger times is mainly insensitive to the mesh size.

Nevertheless, there is an important difference in the numerical behavior of the two schemes in the small Debye length regime. Indeed, when $\lambda^2$ become small, then excepted for very particular values of the data, the variations of $\phi_{\mathcal{T},\tau}$ across the interfaces $\mathcal{E}$ become very large because of (10). Therefore, the drift becomes too large to evaluate its exponential, making the computation with the (SQRA) scheme fail. Since $B(y) \sim -y$ as $y$ tends to $-\infty$, the situation is much less problematic with the (SG) scheme, for which computation of the solution corresponding to $\lambda = 10^{-6}$ is feasible without any specific treatment. However, since the drift becomes large, the use of a reduce time step is required to ensure the convergence of Newton's methods.

The long-time limit of the continuous model has been exhibited in [3]. The model reduces to a nonlinear elliptic equation on the electric potential $\phi$, from which one deduces the concentration profiles. However, no quantitative estimate concerning the convergence towards equilibrium. We then perform a numerical study still with the same parameters as previsously (in particular with $\lambda^2 = 10^{-2}$). The steady solution is computed by choosing a very large final time $T_\infty = 5.10^5$ in the simulation. We denote by $\mathcal{H}_{\mathcal{T}}^\infty$ the corresponding discrete free energy. The relative energy a time $t^n$ is the defined as $\mathcal{H}_{\mathcal{T}}^{\text{rel},n} = \mathcal{H}_{\mathcal{T}}^n - \mathcal{H}_{\mathcal{T}}^\infty$. The energy decay stated in Proposition 3 ensures that $\mathcal{H}_{\mathcal{T}}^{\text{rel},n} \ge 0$ up to numerical errors related to the resolution of the nonlinear systems. One observes on Figure 5 that the (SQRA) scheme dissipates faster energy than the (SG) scheme, the latter exhibiting an almost perfect but rather slow exponential convergence towards the steady state as long as the numerical precision has not been reached. The rigorous proof of such an exponential convergence in the continuous setting can
<figure><img src="image_4.png" /><figcaption>Fig. 2: Concentration profiles at time $T = 5000$ for $\lambda^2 = 10^{-2}$.</figcaption></figure>

be deduced from [16]. Its discrete counterpart should be investigated in future works building on the methodology presented in [2].

**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 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), grants P33010 and 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 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. Bailo, R., Carrillo, J.A., Hu, J.: Bound-preserving finite-volume schemes for systems of continuity equations with saturation. SIAM J. Appl. Math. **83**(3), 1315–1339 (2023).

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

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

4. Cancès, C., Chainais-Hillairet, C., Gerstenmayer, A., Jüngel, A.: Finite-volume scheme for a degenerate cross-diffusion model motivated from ion transport. Numer. Methods Partial Differential Equations **35**(2), 545–575 (2019).
<figure><img src="image_5.png" /><figcaption>Fig. 3: Electric potential profile at times $T = 1$ (solid) and $T = 5000$ (dashed) for $\lambda^2 = 10^{-2}$.</figcaption></figure>

5. Cancès, C., Chainais-Hillairet, C., Merlet, B., Raimondi, F., Venel, J.: Mathematical analysis of a thermodynamically consistent reduced model for iron corrosion. Z. Angew. Math. Phys. **74**(96) (2023).

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

7. Cancès, C., Chainais-Hillairet, C., Fuhrmann, J., Gaudeul, B.: 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).

8. Chatard, M.: Asymptotic behavior of the Scharfetter-Gummel scheme for the drift-diffusion model. In: Finite volumes for complex applications VI. Problems & perspectives. Volume 1, 2, Springer Proc. Math., vol. 4, pp. 235–243. Springer, Heidelberg (2011).

9. Eymard, R., Gallouët, T., Herbin, R.: Finite volume methods. In: Handbook of numerical analysis, Vol. VII, pp. 713–1020. North-Holland, Amsterdam (2000).

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

11. Gerstenmayer, A., Jüngel, A.: 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).

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

13. Heida, M., Kantner, M., Stephan, A.: 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).

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

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

16. Zamponi, N., Jüngel, A.: Analysis of degenerate cross-diffusion population models with volume filling. Ann. I.H.Poincaré-AN **34**, 1–29 (2017).
<figure><img src="image_6.png" /><figcaption>Fig. 4: Convergence of the schemes under space grid refinement, ($\lambda^2 = 10^{-2}$).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 5: Convergence towards the steady long-time behavior in terms of relative energy $\mathcal{H}_{\mathcal{G}}^{\text{rel},n}$ for $\lambda^2 = 10^{-2}$.</figcaption></figure>