Entropy and convergence analysis for two finite volume schemes for a Nernst–Planck–Poisson system with ion volume constraints

In this paper, we consider a drift-diffusion system with cross-coupling through the chemical potentials comprising a model for the motion of finite size ions in liquid electrolytes. The drift term is due to the self-consistent electric field maintained by the ions and described by a Poisson equation. We design two finite volume schemes based on different formulations of the fluxes. We also provide a stability analysis of these schemes and an existence result for the corresponding discrete solutions. A convergence proof is proposed for non-degenerate solutions. Numerical experiments show the behavior of these schemes.


Introduction
Proper modeling of the motion of ions in electrolytes-mixtures of a solvent and N ionic species which can be described by their concentrations c i -and associated simulations are crucial for the development of efficient batteries, fuel cells, and many other applications commonly considered as key technologies for the twenty first century as well as for the understanding of ion channels and other features of biological systems. The classical Nernst-Planck equation is a linear system which for given electrostatic potential , charge number z i and diffusion coefficient D i describes the evolution of the ion concentration c i via The self-consistent electrostatic potential is described by the Poisson equation This model assumes that ions are infinitely small and that the ions of a given species i interact neither with the solvent nor with other ionic species. However, in reality, ion sizes are finite, and ion motion is only possible with a simultaneous displacement of solvent molecules. Moreover, the effective size of ions is increased by the fact that in a polar solvent like water, they are surrounded by a solvation shell consisting of a certain number of solvent molecules. The inclusion of these effects into the model is particularly important for concentrated electrolytes and in electrode boundary layers with high ion concentrations.
Historically, there have been many, often independent attempts to fix this situation, see e.g. the review in [3], the discussion in [30] or [40]. A comprehensive model of ideal mixtures of solvated ions has been derived in [21,22]. In [30,31], a two point flux finite volume discretization approach for these problems has been derived. Various variants of ionic flux approximations have been investigated for the unipolar case, where only one ionic species is considered, in [9], with the result that the flux approximation approach introduced in [30] has several more accurate alternatives. For two of them, we have been able to find appropriate generalizations to the case of several ionic species. One of these generalizations has been independently introduced in [41]. They are analyzed in the present paper using a similar analysis framework. Focus will be set on the cross-diffusion effects arising in the multi-species case, relying on [9] for the coupling with the potential and partial treatment the non-linearity.
In the sequel of Sect. 1, the continuous problem is formulated, and several key properties of the continuous system are discussed. Among these is the decay of an entropy functional for positive solutions.

The Nernst-Planck-Poisson system with finite ionic volumes
Consider a bounded connected polytopal domain ⊂ R d , and finite simulation horizon T > 0. We model the evolution of the concentration c 0 of a solvent and N dissolved species: c i , i ∈ [ [1, N ]]. The mixture satisfies a volume filling constraint where v i are the molar volumes of the species. We will use this constraint using ratios of molar volumes The coefficients (k 1 , . . . , k N ) are parameters of the problem and k 0 is by definition equal to 1. As the molar volumes are not the same, the total concentration is not uniform. The set of positive concentrations c i , i ∈ [ [1, N ]] such that c 0 is positive is denoted by A = (c 1 , ..., c N ) ∈ (0, +∞) N |c 0 := We also introduce the topological adherence of A: For the sake of clarity, we will let C = (c 1 , ..., c N ) ∈ A and often consider c 0 and c as functions of C thanks to (1.1) and (1.2) without clearly expressing the dependency. The dissolved species follow a conservation equation: where z i = z i − k i z 0 the reduced charge number and D i > 0 the diffusion coefficient are parameters of the problem, while h i (C) the chemical potential depends on all the concentrations through: . (1.4) This system is supplemented with Poisson equation for the potential: where c dop : → R is a constant in time doping profile. To simplify the computations, we let c dop = z 0 v 0 + c dop . Using z 0 = 0, we have: To avoid unnecessary complications of the notations, we will drop the tildas for the reduced molar charges as the real molar charges do not appear anymore. Moreover, to simplify the proofs, we will assume that the solvent carries no charge, hence z 0 = 0 and c dop = 0. Treatment of nonzero c dop can be found in [9]. As in [9], we consider a Dirichlet boundary condition for the potential on a nonnegligible part of the boundary D ⊂ ∂ and homogeneous Neumann boundary condition on N = ∂ \ D : where D is assumed to be constant in time and in H 1 ( ) ∩ L ∞ ( ).
The system is supplemented with the following no flux boundary conditions for the concentrations: and with an initial condition C 0 satisfying: (1.8)

Key properties of the continuous system
In this section, we attempt to exhibit the properties of a smooth enough solution (C, ) to the system (1.3)-(1.8) so that calculations are justified. The first property is the conservation of mass. In other words, thanks to (1.3) and (1.7), C satisfies for any t ∈ [0, T ], i ∈ [ [1, N ]]: Moreover, we need the concentrations to be positive for (1.4) to have a sense. In the discrete setting, we will show that the concentrations belong to A. In the continuous setting, it will be assumed. We hint that it might be possible to do it using the entropy method [38], the flux formulation proposed in [30] and [35,Annex A]. Indeed, another key property of the system is the dissipation of a free energy. In this case, the chemical free energy density H (C) is defined as follows: This function is convex, however, the addition of the term −c log c makes the proof quite intricate. This point is detailed in "Appendix A" along with the proof of the following equations: (1.10) The total free energy is formed by the integral of the chemical free energy density and electrical terms: We have using chain rules and (1.9): We also have using chain rules and integrating by part: Notice that we have ∇ · n = 0 on N and = D on D . Using Eq. (1.5), we have: Using this equation and (1.12), we have: Using now Eq. (1.3) and integration by parts, we have the desired Eq. (1.11). Due to the non-negativity of D i c i , E is a Lyapunov functional. Its convexity follows from the assumption C ∈ A (see Lemma A.1).
Finally, we introduce a notion of weak solution that relies on a reformulation of the fluxes: and the space of H 1 functions satisfying the Dirichlet boundary conditions for the potential: More precisely: (1.13) • for all ψ ∈ H D and almost all t ∈ (0, T ), (1.14)

Positioning and outline
The structure of cross-diffusion systems challenges the maximum principle-based methods. In this paper we aim to discretize the system (1.3)-(1.8). For N = 1 this system is a nonlinear drift-diffusion problem and several discretizations have been proposed in [9]. We focus on the extension of these schemes to the more general setting with N > 1 while adapting the proofs to tackle the challenges introduced by cross-diffusion. More precisely, in Sect. 2, the two point flux based finite volume discretization with two variants of the flux approximation is introduced. The main theorems about the existence of discrete solutions and the convergence of approximate solutions are stated. Existence, free energy decay, and positivity of concentrations are proven in Sect. 3, whereas the convergence is proven in Sect. 4. Several 1D and 2D numerical examples showcasing the proven properties of the discretization scheme are discussed in Sect. 5.

Discretization and main theorems
In this section, we propose two discretizations of (1.3)-(1.8) and discrete counterparts of the continuous properties. First, in Sect. 2.1, we state the requirements on the mesh and fix some notations. Then in Sect. 2.2, we describe the common setting for the two schemes to be studied in this paper. These schemes, presented in Sect. 2.3, rely on so-called two-point flux approximations of different formulations of N i . Then in Sect. 2.4, we state our two main results. The first one, namely Theorem 2.1, focuses on the existence of a solution to the nonlinear system corresponding to the schemes for a given mesh, and the dissipation of the energy at the discrete level. More precisely, one establishes that all the studied schemes satisfy a discrete counterpart to Proposition 1.1. Our second main result, namely Theorem 2.2, is devoted to the convergence of the schemes as the time step and the mesh size tend to 0.

Discretization of (0, T) × Ä
In this paper, we perform a parallel study of two numerical schemes based on twopoint flux approximation (TPFA) finite volume schemes. As explained in [23,28], this approach appears to be very efficient for isotropic continuous problems when one has the freedom to choose a suitable mesh fulfilling the so-called orthogonality condition [29,37]. We recall here the definition of such a mesh, which is illustrated in Fig. 1.

Definition 2 An admissible mesh of
is a triplet T , E, (x K ) K ∈T such that the following conditions are fulfilled.
(i) The set T is finite and each control volume (or cell) K ∈ T is non-empty, open, polyhedral, and convex. We assume that (ii) Each face σ ∈ E is closed and is contained in a hyperplane of R d , with positive (d − 1)-dimensional Hausdorff (or Lebesgue) measure denoted by m σ = H d−1 (σ ) > 0. We assume that H d−1 (σ ∩ σ ) = 0 for σ, σ ∈ E unless σ = σ . For all K ∈ T , we assume that there exists a subset E K of E such that ∂ K = σ ∈E K σ . Moreover, we suppose that K ∈T E K = E. Given two distinct control volumes K , L ∈ T , the intersection K ∩ L either reduces to a single face σ ∈ E denoted by K |L, or its (d − 1)-dimensional Hausdorff measure is 0. (iii) The cell centers (x K ) K ∈T belong to their cell: x K ∈ K , and are such that, if K , L ∈ T share a face K |L, then the vector x L − x K is orthogonal to K |L. (iv) For the boundary faces σ ⊂ ∂ , we assume that either σ ⊂ D or σ ⊂ N . For σ ⊂ ∂ with σ ∈ E K for some K ∈ T , we assume additionally that there exists We denote by m K the d-dimensional Lebesgue measure of the control volume K . The set of the faces is partitioned into two subsets: the set E int of the interior faces defined by E int = {σ ∈ E | σ = K |L for some K , L ∈ T } , and the set E ext of the exterior faces defined by E ext = {σ ∈ E | σ ⊂ ∂ } , which can also be partitioned Notice that for any σ ∈ E ext , the edge center x σ satisfying (2) is unique. Given σ ∈ E, we let Thanks to (2) and (2), we have d σ > 0, thus τ σ is well defined. We finally introduce the size h T and the regularity ζ T (which is assumed to be positive) of a discretization (T , E, (x K ) K ∈T ) of by setting Concerning the time discretization of (0, T ), we consider an increasing finite family of times 0 = t 0 < t 1 < . . . , < t N T = T . We denote by t n = t n − t n−1 for 1 ≤ n ≤ N T , by t = ( t n ) 1≤n≤N T , and by h t = max 1≤n≤N T t n . We will use boldface notations for vectors whose number of components is dependent on the mesh while keeping the uppercase notation C when we also consider different species.

A common setting for the Finite Volume schemes
The initial data C 0 which belongs to L ∞ ( ,Ā) thanks to (1.8) is discretized into Notice that previous equation also holds for i = 0 and that this discretization satisfies: (2.2) Assume that C n−1 = c n−1 is given for some n > 0, then we have to define how to compute (C n , n ) = C n K , n K K ∈T . First, we introduce some notations. For all K ∈ T and all σ ∈ E K , we define the mirror values C n K σ and n K σ of C n K and n K respectively across σ by setting (2.3) Given u = (u K ) K ∈T ∈ R T , we define the oriented and absolute jumps of u across any edge by We may now use these operators to describe our scheme. The potential is approximated using the classic TPFA scheme for the Poisson equation: The conservation equation is approximated using a backward-Euler scheme in time: where F n K σ,i should be a conservative and consistent approximation of − D i t n t n t n−1 σ N i · n K σ (n K σ denotes the normal to σ outward K ). Finally, the concentration of the solvent is computed using a discrete version of the volume filling constraint: It remains to define the numerical fluxes F n K σ,i . Two possible choices are given in the next section.

Numerical fluxes for the conservation equations
To close the system (2.4), we have to define the numerical fluxes F n K σ,i . As we intend to use two point flux approximations, they should be of the form: For the sake of readability, we have chosen to define the flux functions F i for unitary D i . Thus this constant should rarely appear in the functional inequalities of the following sections. To preserve the conservation of mass, all the flux functions F i defined afterward satisfy an anti-symmetry property: so that the fluxes are locally conservative, i.e.:

The centered flux
The first numerical flux we consider is based on the original expression of the flux (1.3): The gradient and edge concentration are independently discretized: This flux is a straightforward generalization of the eponymous flux presented in [9]. As such it is also similar to the fluxes introduced in [7,11,12,15,16].

The "Sedan" flux
The other flux under study is also a generalization of the Sedan flux presented in [9]. It originates from and is named after the SEDAN III semiconductor device simulation code [47] and is used to handle the case of degenerated semiconductors in semiconductor device simulators, see [45,46]. In [41], this approach was applied to ion transport in electrolytes, resulting in a scheme almost identical with the one presented here. The scheme relies on the introduction of the excess chemical potential This excess potential characterizes the non-ideality of the electrolyte leading to the following equivalent continuous flux formulation: The Scharfetter-Gummel-inspired discretization [44] of this expression of the flux leads to the so-called Sedan flux: where B(x) = x e x −1 for all x = 0 is the Bernoulli function. Notice that B can be extended by B(0) = 1 and is in C ∞ (R).

Remark
In [9] we studied two other schemes. One was based on the diffusion enhancement and discretization ideas originating from [4]. The extension of this so-called Bessemoulin-Chatard scheme to the multi-species case appears to be not feasible due to the intrinsic use of one-dimensional chain rules. The other scheme based on activity variables and the averaging of the inverse activity coefficient was introduced for the multi-species case in [30]. Numerical analysis of such a scheme is more intricate and would likely not be satisfactory as we were not able to prove convergence in [9]. Moreover, unless more sophisticated inverse activity coefficient averaging strategies are available, this scheme is considerably less accurate compared to all the others discussed in [9].

Main theorems
We have proposed two schemes (2.4), (2.5) supplemented with either (C) or (S). Both schemes are nonlinear systems. Solutions to this nonlinear system should satisfy discrete equivalents of the properties listed in Sect. 1.2, namely conservation of mass and energy-dissipation. For the latter, we introduce the discrete energy functional E T as a discrete counterpart of the continuous energy functional E. It is defined by: (2.7) The first theorem proven in this paper focuses on the existence of discrete solutions for a given mesh, and the preservation of the physical bounds: positive concentrations, and the properties of Sect. 1.2.
be an admissible mesh and let C 0 be defined by (2.1). Then, for all 1 ≤ n ≤ N T , the nonlinear system of equations (2.4), (2.5) supplemented with either (C) or (S) has a solution Moreover, the solution to the scheme satisfies, for all 1 ≤ n ≤ N T , The proof of this theorem is the purpose of Sect. 3. Knowing a discrete solution to the scheme, (C n , n ) 1≤n≤N , we can define an approximate solution (C T , t , T , t ). It is the piecewise constant function defined almost everywhere by This definition will be developed in Sect. 4 and supplemented by other reconstruction operators.
Using this existence result, we let be a sequence of admissible meshes in the sense of Definition 2 and associated approximate solution. We assume that h T m , h t m −→ m→∞ 0 while the mesh regularity remains bounded, i.e., ζ T m ≥ ζ for some ζ > 0 not depending on m. A natural question is the convergence of (C T m , t m , T m , t m ) towards a weak solution to the continuous problem. The convergence result is stated in Theorem 2.2 which will be proved in Sect. 4. m≥1 satisfies, up to a subsequence:

Fixed Mesh analysis
In this section, we intend to prove Theorem 2.1. To this end, we will use a topological degree argument in Sect. 3.3. This topological degree relies on properties of the fluxes and a priori estimates detailed respectively in the following section and in Sect. 3.2.
The methodology of this proof is very similar to the one done in [9]. The key changes and improvements are concentrated in Proposition 3.2, Lemmas 3.2 and 3.5.

Analysis of numerical flux based functions
In this section, we introduce several functions derived from F i . As in [9], the first functions of interest models the free energy dissipation for each species i ∈ [[1, N ]]: We also introduce the local free energy dissipation D := N i=1 D i . In addition to this function, we can define a reconstruction of the concentration at the interfaces. This is the purpose of the following lemma:

Lemma 3.1 For a flux F i defined either by (C) or (S), the corresponding face concentration functions defined by
Proof The proof of the extension by continuity and the average property (3.2) is highly similar to [9, Lemma 3.1]. For the centered scheme defined by (C), we have by definition: hence the extension by continuity and Eq. (3.2). For the Sedan scheme, defined by (S), we introduce and notice that: Using the following property of the Bernoulli function: we have: Finally using (3.3) and the differentiability of B, we have the desired extension on We also have Eq. (3.2) thanks to the monotony of B and the relation Thanks to this lemma, D i rewrites: (3.5) This new formulation along with (3.2) grants the non-negativity of D i and D. The following coercivity lemma gives more detailed information on the behavior of D: We have, for all δ, , M > 0: As the proof of this lemma is purely technical it has been relegated to Appendix B.

A priori estimates
In this section, we intend to establish uniform a priori estimates on the concentration and the potential, in order to prove the existence of solutions that satisfies the properties of Theorem 2.1.
We assume that we dispose of (C n , n ) n∈ [[0,N max ]] solution of (2.1), (2.4), (2.5) supplemented with either (C) or (S) in A T × R T . Where A, the adherence of A is the set of non-negative concentrations c 0 , ...c N satisfying the volume filling constraint. The first a priori estimate is the conservation of mass (2.9): The proof is straightforward and classical thanks to the local conservativity of the fluxes, the no flux boundary conditions, and the discretization choice for C 0 .
We can also build a discrete equivalent to Theorem 1.1 using E T defined in (2.7) and the dissipation function D i . This is the purpose of the following proposition: Proof The proof is fairly classical once noticed that thanks to Lemma Notice that the left-hand side is the term of interest, we will then focus on the reformulation of the right-hand side. We multiply Eq. (2.4b) by h i (C K ) + z i K and we sum over the cells and species in order to get the following three-terms formula: The term concerning the chemical energy, t n T chem , appears directly in (3.8), thus we focus on T el . Using Eq. (2.4a), we have: (3.10) which is the second line of Eq. (3.8). For T diss , a discrete integration by parts (summation by parts) yields: Using this equation and equations (3.10), (3.9) in (3.8), we have (2.8): which concludes the proof thanks to the preliminary remark.
In the following lemma, we will show several bounds on the potential and then take advantage of them to get a bound on the free energy dissipation: Proof The proof of (3.11) is a straightforward application of [9, Proposition A.1]. As the proof of (3.12) is detailed in [9, Lemma 3.6], we focus on the proof of (3.13), assuming (3.12). Multiplying Eq. (2.4a) by n K and summing over K ∈ T yields, using (2.3): Using Eq. (3.12), (3.11), and C n ∈ A T , we have the desired result. The last result is based on (3.7). Summing that equation, we have: We have thanks to Eqs. (1.10), (3.12), and (3.13): Hence the desired result, up to the choice of a bigger constant M * .
Finally, we use the free energy dissipation result (3.14), and the estimates on the free energy dissipation functional to improve the assumption C n ∈ A T .

Lemma 3.5
There exist 0 , 1 , ..., N positive, depending on, among other things, C 0 and decreasing with min t and min σ ∈E τ σ such that: Proof The proof follows the idea of [14, Lemma 3.10] (see also [15,Lemma 3.7], [9, Lemma 3.7]). We start with the proof for i = 0 and a fixed time step n using ϒ δ,M , then treat the case of i ∈ [[1, N ]] using δ, 0 ,M ,i and finally notice that no assumptions were made on n. Thanks to assumption (1.8) on the initial concentrations, and Lemma 3.3, we dispose of K ∈ T such that: It is well defined thanks to the monotony of ϒ and Lemma 3.2. Moreover, we have for every cell L sharing an edge with K : thanks to the positivity of D i and Eq. (3.14). Similarly we recursively define: 16) and notice that thanks to the connectivity of there exist l ≤ card(T ) such that, for all L ∈ T : c n L,0 ≥ δ l .
Hence a possible choice for 0 . As explained above, the proof is exactly the same for i ∈ [ [1, N ]], with the use of δ, 0 ,M ,i instead of ϒ δ,M in Eq. (3.16) and again does not depend on the time step n ≥ 1.

Existence of solutions
Using the estimates of the previous section we can establish the existence of a solution to our numerical scheme. Thanks to Proposition 3.1 and Lemmas 3.3 and 3.5, this will conclude the proof of Theorem 2.1.
Proof As in [9, Proposition 3.8], we use induction and a topological degree argument to transform continuously the non-linear system (2.4), (2.5) to a linear one. However, the path presented in [9] is no longer valid as we do not have a monotony property on h i . The homotopy follows 3 steps. The first one is sketched in "Appendix C", the second one changes the discretization while maintaining k i , D i to 1 and the potential to zero. The last step corresponds to the activation of the potential and the remaining nonlinearities.
Following these ideas, we follow the zeros of a homotopy H α for α ∈ [0, 3]: which is a standard finite volume scheme for the heat equation for α = 0, and our scheme for α = 3. Thanks to (2.4c), for every α, c 0 is eliminated, and proving c 0 . . . c N > 0 allows to conclude that uniformly in α, Step 1: implementation of the solvent effects using an ad hocscheme.
where D α is set to zero. As expressed in Lemma C.1 we dispose of 1 such that the zeros of H α have a concentration that is bounded away from zero by 1 .
Step 2: change of scheme without potential and for identical species. We change the discretization of c i ∇ log(c i /c 0 ). For α ∈ [1,2], H α = 0 writes: where k i,α is set to 1 and D α is again set to zero. Here again we dispose of 2 such that the zeros of H α have a concentration that is bounded away from zero by 2 .
Step 3: activation of the potential and the difference between the species. For α ∈ [2,3], H α = 0 means: where D α is set to (α − 2) D and k i,α to 3 − α + (α − 2)k i . Thanks to Lemma 3.5, we dispose of 3 such that the zeros of H α have a concentration that is bounded away from zero by 3 . Conclusion. Using a topological degree argument [20,39], we can derive the existence of a solution for H α = 0 for α = 3. A detailed development of this argument applied to a finite volume scheme for a scalar nonlinear convection-diffusion problem can be found in [26], proof of Lemma 4.1. It is straightforward to state the existence of a solution for H 0 = 0-the classical finite volume scheme for the heat equation. This leads to a nonzero toplogical degree of H 0 . We have shown uniform bounds on the concentrations from below by min( 1 , 2 , 3 ) and above by 1 v i and the potential: ±M . Continuity of H α and the homotopy invariance of the topological degree yield a nonzero topologial degree of H 3 and thus the existence of a solution for H 3 = 0.

Remark 3.2
Concerning uniqueness of solutions, we think that, given the corresponding results for semiconductors found e.g. in [34,36], it is worth to investigate uniqueness for the thermodynamic equilibrium, and for small times/applied voltages around thermodynamic equilibrium.

Convergence
In this section we prove Theorem 2.2, which states the convergence of our schemes towards a weak solution. We consider a sequence T m , E m , (x K ) K ∈T m m≥1 of admissible meshes with h T m , h t m tending to 0 as m tends to +∞, while the regularity ζ T m remains uniformly bounded from below by a positive constant ζ .
Thanks to Theorem 2.1, we have a family of discrete solutions (C m , m ) m . We will first propose different reconstructions of approximate solutions in Sect. 4.1, then we show several compactness properties in Sect. 4.2 in order to obtain the convergence of a subsequence of approximated solutions. Section 4.3 is then devoted to the identification of the limit as a weak solution.
To enlighten the notations, we will remove the subscript m as soon as it is not necessary for understanding.

Reconstruction operators
In order to carry out the analysis of convergence, we introduce some reconstruction operators following the methodology proposed in [25].
The operators π T : R T → L ∞ ( ) and π T , t : These operators allow passing from the discrete solution (C n , n ) 1≤n≤N T to the approximate solution since To carry out the analysis, we further need to introduce an approximate gradient reconstruction. Since the boundary conditions play a crucial role in the definition of the gradient, we need to enrich the discrete solution by face values C n σ σ ∈E ext ,1≤n≤N and n σ σ ∈E ext ,1≤n≤N defined by C n σ = C n K σ and n σ = n K σ for σ ∈ E ext ∩ E K . With a slight abuse of notations, we still denote by C n = (C n K ) K ∈T , (C n σ ) σ ∈E ext and n = ( n K ) K ∈T , ( n σ ) σ ∈E ext the elements of A T ∪E ext and R T ∪E ext containing both the cell values and the exterior faces values of the concentration and the potential respectively.
For σ = K |L ∈ E int , we denote by σ the diamond cell corresponding to σ , that is the interior of the convex hull of σ ∪ {x K , x L }. For σ ∈ E ext , the diamond cell σ is defined as the interior of the convex hull of σ ∪ {x K }. The approximate gradient ∇ T : R T ∪E ext → L 2 ( ) d is piecewise constant on the diamond cells σ , and it is defined as follows: We also define ∇ T , t : This reconstruction is merely weakly consistent (unless d = 1) and takes its source in [17,27]. More consistent reconstruction operators will be introduced in Sect. 4.3. Let us recall now some key properties to be used in the analysis. First, for all u, v ∈ R T ∪E ext , This implies in particular that With slight abuse of notations, we extend the reconstructions operator presented above to the vectors of concentrations in A T , A T ×N T , etc.

Compactness
In this section we intend to prove a discrete H 1 estimate on the concentrations using the bound on the free-energy dissipation (3.14). To that extend we will introduce a chemical dissipation D chem as a discrete equivalent to c i |∇h i (c)| 2 and compare it both with the usual distance and the total dissipation D.
As the identification of the limit is only possible for the results of this section are proved under this assumption and complemented with remarks indicating whether the hypothesis is necessary or not. In order to apply chain rules for the convergence, we need to change the face concentration C from the one defined by the numerical scheme through Lemma 3.1 to the logarithmic average: This choice of edge concentration will also be used in the definition of D chem to avoid a dependency on the potential. The following lemma provides an estimate the numerical-flux based averages using this logarithmic average.
Proof For the centered scheme, this inequality is known with α = 1 without assumption on c 0 [42]. For the Sedan scheme the proof is more intricate and uses the hypothesis on c 0 . Equation (4.3) is equivalent to the boundedness of , as in the proof of lemma 3.1. By symmetry, one can assume x i ≥ 0 and thanks to our assumption on the solvent and the potential, y i is bounded by some K . Moreover, we notice that by definition of x i , (3.4) yields: so that we have: Remark This result does not hold for = 0. In that case, we have to introduce: for which this lemma is trivially true with α = 1. That third edge reconstruction should be used in the following definition.
Then we try to take advantage of Proposition 3.1. As Lemma 3.4 already provides satisfying estimates on , we introduce A first interesting result is that D chem is a semimetric on A. The non-negativity and symmetry properties are trivially satisfied, the last property is the subject of the following lemma.

Lemma 4.2 We have D chem (C K , C L ) = 0 if and only if C K = C L
Proof If C K = C L , we obviously have D chem (C K , C L ) = 0, we will then focus on the other implication. Assume that we dispose of C K , C L in A such that D chem (C K , C L ) = 0. We let for i ∈ [[0, N ]]: thus D chem is the sum of non-negative terms. As we have D chem (C K , C L ) = 0, we have: Assume that a K ,0 = a L,0 , then The other case is absurd: using the symmetry of D chem , one can freely assume that a K ,0 > a L,0 . Using k i > 0, we have a K ,i > a L,i ∀i ∈ [[1, N ]] hence: The function D chem cannot be extended by continuity onto A 2 . Some information for near zero concentrations can be inferred from lemma 3.2. The following sequential result means that the semi-metric property is preserved near the boundary ∂A 2 .

Lemma 4.3 Let
Proof For the sake of simplicity, as this result will only be used with a lower bound on c 0 , we keep the proof to this simpler case and assume that inf(c l K ,0 , c l L,0 ) > 0. To prove the limit, we will show that from any sub-sequence, we can extract a sub-sub-sequence such that C l K − C l L → 0. Considering any sub-sequence, thanks to the boundedness of A we can extract a sub-sub-sequence such that C l K and C l L converge. If we dispose of i ∈ [ [1, N ]] such that c l K ,i → c * > 0 while c l L,i → 0 (or the symmetric situation), then we have: and so that D chem,i (C l K , C l L ) ∼ l→∞ −c * log(c l L,i ) → ∞, which is absurd. Necessary, we have c l L,i → 0 if and only of c l K ,i → 0. If {i| inf l (c l L,i ) > 0} is empty we have the desired result. Else, we use D chem,i (C l K , C l L ) → 0, to get: where a K ,i , a L,i are defined by (4.4). As both c l K ,i and c l L,i are bounded away from zero, a l K ,i and a l L,i are convergent up to a subsequence. Using our assumption on c 0 , we have the convergence of a l K ,0 and a l L,0 up to yet another subsequence. Then we conclude using the proof of previous lemma.

Remark
The main simplification brought by the restriction to the setting c 0 > is the use of C i in (4.5). The end of the proof can be adapted by considering the behaviour of a l K ,0 − a l L,0 . This semi-metric is however not commonly used and the following lemma intends to compare it with the usual distance.  and use the boundedness of A to extract a convergent sub-sequence of (C n K , n ) and denote (C * , * ) its limit. As c i is bounded, we have D chem (C n K , C n L ) → 0. Thanks to Lemma 4.3, we have * = 0 so that we will consider first order development in n . We notice that the blow-up of the ratio implies that: For the sake of readability, we will drop from now on the superscript n . We have to consider three cases: we dispose of species such that j = o(| |) and c * j = 0, but for all of them log j +c j c j remains bounded; (3) we dispose of a specie such that j = o(| |), c * j = 0, and up to a subsection, log j +c j c j blows-up.
Preliminary remark about the solvent We consider j ∈ [ [1, N ]] such that c * j > 0 and let j = N i=0 i . If such a j cannot be found, then c 0 → 1 v 0 and the conclusion (4.8) holds. Else we have, thanks to (4.7): so that: thus: First order development of h j gives: Thanks to (4.7) we have the estimation: To correct the effect of the species with negligible concentrations, we let: By construction, we have c = c + o (1). Using the hypothesis (4.2) we have = + o(| |) and 0 = 0 + o(| |). These three results and Eq. (4.9) yield: We let ξ j = j c j − c for j = 0 and ξ 0 = 0 c 0 − c . Previous equation yields: Considering c * j >0 c j ξ j , we have: so that: We conclude the first part of the proof with the following estimate that follows from (4.7): For the sake of readability, we will drop the over 0 in the second part of the proof, introduce S = { j ∈ [[1, N ]]|c * j > 0}∪{0}, and assume by symmetry that ≥ 0. We have: so that (4.10) becomes: The Cauchy-Schwarz inequality yields We intend to use ideas presented in [2] to improve the estimation of . More precisely, the stability version of the Cauchy-Schwarz presented in [1] gives: We intend to show that X |X | − Y |Y | is bounded away from zero. To show this bound we let: and consider a minimizing sequence of K under the conditions As K is invariant by scaling, we can assume that we have a convergent minimizing sequence C l inf , l inf of limit C * inf , * inf and of norm equal to 1. Note that we do not assume C ∈ A, nor C * inf > 0 thus we consider broader options than necessary for use in (4.11) to ensure existence of the minimum. Finally, we notice that K is non-negative, thus its infimum is either zero or positive. We will prove the positivity by contradiction.
Assume that the limit of K (C l inf , l inf ) is zero, we show that |Y l inf | is convergent up to a subsequence. We consider j such that, up to a subsequence, is bounded away from zero. If c * inf, j = 0, |y l inf, j | is bounded thus |Y l inf | is too, and up to another subsequence, the latter is convergent. If c * inf, j = 0 we notice that x l inf, j → 0 and |X l inf | is bounded away from zero, so that → 0, which is absurd. We let γ be the limit of |Y l inf | 2 . As we have assumed the infimum to be zero, we have: * This would imply that * inf is non-negative, however, we have * inf,0 = − j∈S\{0} k j * inf, j and * inf is of norm 1. This is absurd, hence the infimum cannot be zero. Thus we dispose of 0 < α depending only on k 1 , . . . k N and S such that: As we have assumed (using symmetry) ≥ 0, we also have α ≤ 1. So that we have: Thanks to Eqs.

Conclusion of the proof in case 4.2
We dispose of j such that c j → 0 and j = o(| |), thus have up to a sub-sequence: The assumed boundedness of log Moreover, we also dispose of α = min(1, inf n c n j + n j c n j ) > 0 such that: Necessary, we have log which is bigger than | | 2 and thus contradicts (4.7).

Conclusion of the proof in case 4.2
Let j be such that j = o(| |), c * j = 0, and log j +c j c j blows-up. We have: and: so that: which contradicts (4.7) since j = o(| |).

Global conclusion
As each of the cases lead to a contradiction, we have the desired inequality for i ∈ [ [1, N ]]. For the solvent, we see that: thus the announced result up to the choice of a bigger constant M.
Remark Even though this lemma is proved without the assumption on c 0 , its adaptation to C is the most difficult. The shortest sketch of proof we obtained relied on the splitting of case (4.2) in more than 15 sub-cases.
Using these tools, we may now prove the following necessary compactness inequality: Proof We will show the result for i ∈ [ [1, N ]] and use the definition of A to extend it the solvent. For improved readability, we will drop the subscript m, and for σ = . By definition, we have: Thanks to Proposition 4.1 and Lemma 4.1, we have: It is sufficient to bound N T n=1 t n σ ∈E int τ σ C n σ, j (D σ h j (C n )) 2 , for all j ∈ [ [1, N ]] to get the desired result. We have: Thanks to Eq. (3.14) of Lemma 3.4, we dispose of M such that: Moreover, σ ∈E int τ σ C n σ, j (z j D σ n ) 2 is also bounded thanks to (3.12) and the L ∞ bound on C i and z i . Thus we have: (4.12) which in turn yields the desired result.
For the solvent we notice that: so that the bound on all ∇ T m , t m c i transfers into a bound on ∇ T m , t m c 0 .
Using this discrete L 2 (H 1 ) estimate, we use a discrete Aubin-Lions lemma to get the compactness of the sequence of solutions, as stated in following proposition: Proof For improved readability we drop again the subscripts m. The proof of the first two result relies on a discrete Aubin-Lions lemma [33,Lemma 3.4]. We intend to use it in the setting described in [10, Lemma 9]. Proposition 4.2 provides a first property, but we still have to prove that there exist C independent of the mesh such that n c n i − c n−1 i T ,−1 ≤ C, where · T ,−1 is defined by duality: Let ϕ ∈ R T . Tanks to (2.4b), we have: Using the definition of F n K σ,i along with the definition of C i respectively Eqs. (2.5) and (3.1), we have: Thanks to the Cauchy-Schwarz inequality, we have: Using C i ≤ 1 k i v 0 , the definition of D i and ∇ T ϕ 2 L 2 ≤ 1, we have: Using the Cauchy-Schwarz inequality and the Eq. (3.14) of Lemma 3.4, we have: This concludes the proof of Eqs. (4.13) and (4.14).
We may now focus on the convergence of the potential. The existence of satisfying (4.15) is a straightforward consequence of (3.11). Similarly, (3.12) implies the existence of a vector field u such that We have to identify u with ∇ . We let w ∈ C ∞ c (Q T , R d ) and define: and the associated diamond-cell reconstruction: Thanks to the smoothness of w, we have convergence of w E, t toward w and: Using the geometric relation d σ m σ = dm σ and the definition of w n σ , we have: Thanks to the smoothness of w and the convergence of , we have: This concludes the identification of u and the proof of (4.16).
These convergence topologies are sub-optimal and will be improved in Lemma 4.4. First, we notice that for the concentrations, we also dispose of edge averages C and C defined in Eqs.
Similarly, we introduce c E, t,i . As we expect, these reconstructions are convergent and share their limit with π T , t c i . This is the main purpose of the following lemma.  Finally, we show a weak-convergence property on the gradients of the logarithms: Lemma 4.5 Let C be as in Proposition 4.3. We have: Proof Let us start with the proof on Eq. (4.22). By definition (4.2), we have: so that, using (4.20), (4.14), and the assumed bound on c 0 we have: We notice that the limit c 0 is also bounded away from 0 and conclude using the continuous chain-rule. For (4.21), we proceed similarly. Notice that since c ≥ 1 v 0 max k i > 0 the bound does not need to be assumed. We only need the strong L 2 convergence of the reconstruction using the logarithmic average on the diamond cells. This is an application of Lemma D.2, as in the proof of Lemma 4.4.

Identification
In this section we will identify the limits obtained in Proposition 4.3 as weak solutions in the sense of Definition 1. First we improve the convergence topology on the potential and identify it as a weak solution of the Poisson equation.
∇ψ(·, t n ) uniformly in , ∀n ∈ {1, . . . , N T }, (4.23) thanks to the smoothness of ψ. The operator ∇ is also such that The scheme (2.4a) then reduces to Integrating with respect to time over (0, T ) and passing to the limit h T , h t → 0 thanks to Proposition 4.3 Eqs. (4.13) and (4.16) and Eq. (4.23) we have: and continuity of the linear application, we have: In particular, (1.14) holds for almost every t ∈ (0, T ).
Concerning the boundary conditions for , the fact that = D on (0, T ) × D can be proved for instance following the lines of [6,Section 4].
The following theorem focuses on the identification of C as a weak solution satisfying (1.13). As announced in Theorem 2.2 this can only be done with an assumption on the solvent. Remark 1.1 is a first clue of the validity of this assumption. For positive initial condition, this assumption is valid in all the numerical test. In the 1D setting and under a CFL condition, it might be possible to prove it through improvements of Lemmas 3.2 and 3.5. This could be the topic of further research.
where ϕ n−1 K σ = 0 for σ ∈ E ext and C σ is defined by Lemma 3.1. The treatment of terms T 1 and T 3 follows the lines of [9,Proposition 4.7].
More precisely, for T 1 we use the discrete integration by part and notice that ϕ N T = 0 to pass the time derivative on ϕ. Thanks to the smoothness of ϕ and the convergence of π T , t c i and π T c 0 i , we have: For T 3 , we extend in time the reconstruction introduced in the proof of Proposition 4.4 and notice that: The treatment of the term T 2 is more intricate. First we let T 2 be the same term with a different edge concentration: where C n σ,i = C i (C n K , C n K σ ) is the logarithmic mean introduced in (4.2). We will first prove the convergence of T 2 then identify its limit. To this end, we set: For term T 2,1 we use the chain rule C n σ,i D K σ log(c n i ) = D K σ c n i and get: Thanks to the weak convergence of ∇ T m , t m c i and the strong convergence of ∇ T m , t m ϕ, we have: For the other terms, we need the enhanced convergence of gradients provided by Lemma 4.5. So that the terms T 2,2 and T 2,3 have the following limits: Let us now establish that T 2 and T 2 share the same limit.
Thanks to the triangle and Cauchy-Schwarz inequalities, one has The first term in the right-hand side is uniformly bounded thanks to (4.12). Thus our problem amounts to show that Let us reformulate R as Thanks to Lemma 4.1, the quantity 1 − is uniformly bounded, whereas the regularity of ϕ implies that D σ ϕ n−1 ≤ ∇ϕ ∞ d σ . Putting this in the above expression of R, we obtain that thanks to Lemma 4.4. Thus T 2 and T 2 share the same limit, which gives the announced result.

Remark 4.1
As the proof is based on compactness we cannot extract convergence rates for our solutions. The numerous non-linearities involved would render a proof of a convergence rate difficult. The best rate one could hope for given the reconstruction introduced is 1.
To get the second order "super convergence" shown in Fig. 6, we have used a P 1 reconstruction. To use this technique we need to use Voronoï meshes with cell centers on the boundary of our domain so that the dual Delaunay mesh covers the whole domain. Due to the specific treatment of boundary conditions in such meshes (formally, τ σ = +∞) they are not covered by our analysis.
In meshes simply satisfying Definition 2, the key idea is to project the reference solution to the space of constant by cells function, i.e. to consider only the discretization errors. The projection to this P 0 space could be done either using L 2 projections (as for the discretiation of the initial condition) or by taking the value at the cell-centers (to acknowledge the degree of freedom in their choice for orthogonal meshes or simplify the code).

Numerical examples
The numerical examples have been implemented in the Julia language [5] based on the package VoronoiFVM.jl [32] which realizes the implicit Euler Voronoi finite volume method for nonlinear systems of diffusion-convection-reaction equations on simplicial grids. The two two schemes (2.4), (2.5) supplemented with either (C) or For this purpose we use Newton's method with analytical Jacobians with optional homotopy embedding. An advantage of the implementation in Julia is the availability of ForwardDiff.jl [43], a forward mode automatic differentiation package based on dual number arithmetic. This package allows the assembly of the analytical Jacobians based on the implementation of functions calculating two point fluxes and charges, without the need to write source code for derivatives. For the one-and two dimensional examples in this paper, the resulting linear systems are solved using UMFPACK [19] as Julia's built-in sparse direct solver.

Let
= (0, L) with L = 20. As an initial state, assume a binary electrolyte with two ionic species with opposite charges and a solvent. At moment t = 0, we assume a spatially constant, electroneutral distribution of the ions. We apply a potential difference via Dirichlet boundary conditions | x=0 = −10 and x=L = 10 and solve the Poisson equation with these data as initial value. We set homogeneous Neumann boundary conditions for both ionic species. With starting time step size t = 10 −3 we start the evolution until the species distribution reaches its equilibrium under the applied potential difference. As discussed in [8], the time step sizes are controlled such that the energy dissipation per time step is limited: Figure 2 shows the evolution in the case v 0 = v 1 = v 2 = 1, z 0 = 0, z 1 = 1, z 2 = −1. At the end of the time evolution, most of the ions are accumulated in their respective polarization boundary layers, almost completely displacing the solvent. As predicted, the ion concentration is bounded by 1. The computation used the flux (S). Figure 3 shows the evolution in the case v 0 = v 1 = v 2 = 1 and z 0 = 0, z 1 = 2, z 2 = −1. Once again, at the end of the evolution, anions and cations pile up in the corresponding boundary layers. Ion concentrations are bounded by 1, but due to the larger charge of the cation, the corresponding boundary layer becomes smaller.  Figure 4 shows the evolution in the case v 0 = v 2 = 1, v 1 = 2 and z 0 = 0, z 1 = 1, z 2 = −1. Once again, at the end of the evolution, anions and cations pile up in the corresponding boundary layers but now, the cation concentration is bounded by 1 2 . The corresponding evolution of the relative free energy E(t) − E ∞ is shown in Fig. 5. We observe an exponential decay and almost equal behavior for both variants of the flux approximation (S) and (C). Moreover, the time step control algorithm keeps the dissipation per timestep below the intended limit. On the grid of 100 nodes the whole computation took 700ms on an Intel(R) Core(TM) i7-9850H CPU with 2.60GHz. The Newton method used 14 iterations for the first timestep, and at most 4 iterations for the remaining time steps.

1D stationary convergence test
In the same domain as above, we set v 0 = 1, v 1 = 2, v 2 = 1, and z 1 = 1, z 2 = −1. This time, we look for the stationary solution with homogeneous Dirichlet boundary conditions for , and Dirichlet boundary conditions for the concentrations. These

An electrolytic diode
The second example regards a domain = (0, W ) × (0, L) with W = 2 and L = 10. We assume z 0 = 0, z 1 = 1, z 2 = −1 and v 0 = 1, v 1 = 4, v 2 = 4. At y = 0 and y = L we fix concentrations to a value c 1 = c 2 = 0.01 We set | y=0 = 0 and apply a changing value bias at y = L. At x = 0 we apply symmetry (homogeneous Neumann) boundary conditions for , c 1 , c 2 . Homogeneous Neumann boundary conditions are also applied for c 1 , c 2 at x = W . We set Neumann boundary conditions λ∇ · n = q(y) at x = W , where with σ = 5. Figure 7 shows three different states of the electrolytic diode. Figure 8 (left) shows the corresponding current-voltage curve. We see a well developed rectification effect: At reverse bias, ion concentrations under the charged surface are rather low, resulting in low conductance and low ionic current. Whereas at forward bias, larger ion concentrations lead to a larger ionic current. On the 2D grid with 1681 nodes, the whole computation for the IV curve took 19 seconds on the aforementioned system. (1.10) We still have to establish the bounds (1.10). To that end, we notice that: As c is positive and 0 ≤ c i c ≤ 1, we have H (C) ≤ 0. For the lower bound, we notice that − N i=0 c i c log c i c can be interpreted as the entropy of a random variable over a set of N + 1 elements. It is common knowledge that it is maximal for c i c = 1 N +1 thus: Finally, notice that which is the desired bound.
The proof for the centered scheme relies on expression (3.5): We notice that h i (C K ) + z i K − h i (C L ) − z i L blows up and that C i ≥ c K ,i 2 , hence the blow-up of the limit.
For the Sedan scheme, it is more intricate and try to bound F i (C K , C L , K , L ) away from zero to take advantage of the blow-up of h i (C K ) + z i K − h i (C L ) − z i L . The positivity of the product, ensures that the limit will have the right sign. Let δ, , M > 0, i ∈ [ [1, N ]]. We denote by O c the set: We notice that the hypothesis c 0 > yields a bound on ν i . Moreover, this bound is uniform in c. We intend to use this bound to prove that the flux function defined by (S) is bounded away from zero. We let: We have, for all (C K , C L , K , L ) ∈ O c : hence F i is bounded away from zero for c small enough and the desired result.

Limit of 7 ı,M
In this section, we prove the remaining limit: Notice that for all (C K , C L , K , L ) ∈ O c , we dispose of i ∈ [ [1, N ]] such that c L,i ≥ 1−v 0 c N k i v 0 . Notice also that ϒ δ,M is increasing, it is then sufficient to prove the limit for a given sequence. Let c n be sequence that steadily decreases to zero such that for all n ∈ N , c n ≤ 1 2v 0 and there exist i ∈ [ [1, N ]], (C n K , C n L , n K , n L ) ∈ O c n satisfying: D(C n K , C n L , n K , n L ) ≤ ϒ δ,M (c n ) + 1 n and c n L,i ≥ 1 2N k i v 0 .
We have, using c n L,0 ≤ c n , c n K ,i ≤ 1 k i v 0 , 1 v 0 max k j ≤ c ≤ 1 v 0 min k j , c n K ,0 ≥ δ, the bounds on , and c n L,i ≥ 1 2N k i v 0 : As all the terms are bounded except log(c n ) which goes to −∞, we have blow-up of h i (C n K ) + z i n K − h i (C n L ) − z i n L . For the centered scheme, we use C i (C n K , C n L , n K , n L ) ≥ 1 4N v i , and we have: hence the desired result. For the Sedan scheme, we will also only consider D i , but we need a more precise approach: as in previous section, we bound the flux away from zero. We let: We have: As the right-hand side tends to −∞, the left-hand side is bounded away from zero. Using the previously detailed arguments, we have the desired limit.

Appendix C: Study of a numerical scheme for h i = log(c i ) −˛log(c 0 )
To prove the existence of solutions to the Sedan and centered scheme, we introduce this simplified cross diffusion system where the coupling occurs only through the solvent using the chemical potential defined above. This system is discretized using the ideas of the centered scheme and [13]. In detail, we use equations (2.4b), (2.4c) with k i , D i = 1, z i = 0, and: We may now use the decay of this entropy to prove the desired result for i = 0. To that extent we proceed as in Lemma 3.5 and see that D K L c n 0 D K L log c n 0 is clearly coercive in the sense of lemma 3.2 while the part in α has the right sign. This yields the uniform bound for c 0 .
The bound for c i relies on the entropy H = N j=1 c j log(c j ) + αc 0 log(c 0 ). As for Lemma A.1, this entropy restricted to A is convex and its derivatives as a function of R N are the chemical potentials, Thus multiplying the conservation equation by h yields: This new dissipation is also coercive in the sense of Lemma 3.2, thus we can proceed as in Lemma 3.5 to get the announced bounds.
Let p be as in the lemma and q its Hölder conjugate. We have: This concludes the proof of the lemma. where q is the Hölder conjugate of p. Using the assumed estimation of the gradient, we have the announced result.