Clément Cancès
email: clement.cances@inria.fr

Benoît Gaudeul
email: benoit.gaudeul@univ-lille.fr

A CONVERGENT ENTROPY DIMINISHING FINITE VOLUME SCHEME FOR A CROSS-DIFFUSION SYSTEM

Keywords: Cross-diffusion system, Finite Volumes, discrete entropy method, convergence AMS subject classifications. 35K51, 65M08, 65M12

We study a two-point flux approximation finite volume scheme for a cross-diffusion system. The scheme is shown to preserve the key properties of the continuous systems, among which the decay of the entropy. The convergence of the scheme is established thanks to compactness properties based on the discrete entropy -entropy dissipation estimate. Numerical results illustrate the behavior of our scheme.

Introduction.

1.1. The system under study. The system studied in this paper has been originally introduced by [START_REF] Bakhta | Cross-diffusion systems with non-zero flux and moving boundary conditions[END_REF] to model the production of solar panels using vapor deposition. In this system, we study the diffusion of N species whose respective concentrations are U = (u 1 , . . . , u N ) in a (nonempty) connected bounded open domain Ω of R d for a fixed time T . We denote by Q T = (0, T )×Ω. The diffusion occurs through exchanges between different species which are quantified by the matrix A = (a i,j ) of cross-diffusion coefficients. It leads to the following system of partial differential equations:

(1.1)

∂ t u i -div   N j=1
a i,j (u j ∇u i -u i ∇u j )

  = 0 in Q T for i ∈ [[1, N ]].
The matrix A is assumed to be symmetric with nonnegative coefficients, i.e. a i,j = a j,i ≥ 0. A does not depend on U and thus differs from the diffusion matrix D(U ) = (d i,j (U )) defined by

d i,j (U ) = δ i,j k =i a i,k u k -a i,j u i ,
where δ i,j stands for Kronecker symbol, such that the problem (1.1) rewrites

(1.2) ∂ t U -div (D(U )∇U ) = 0.
System (1.2) enters the family of the nonlinear cross-diffusion systems since D depends on U and has nonzero off-diagonal entries. Challenges both from the analytical and numerical points of view come from the presence of off-diagonal zeros in A. In the previous contributions [START_REF] Burger | Nonlinear Poisson-Nernst-Planck equations for ion flux through confined geometries[END_REF][START_REF] Gerstenmayer | Analysis of a degenerate parabolic cross-diffusion system for ion transport[END_REF][START_REF] Cancès | Finite-volume scheme for a degenerate cross-diffusion model motivated from ion transport[END_REF], the zeros are integrated through the assumption that the cross-diffusion occurs with and only with a solvent specie. Until Section 4 we will not make any assumption about the zeros of A. A non-degeneracy assumption will be further assumed in Section 4, but our convergence result could extend to the particular cross-diffusion matrices considered in [START_REF] Gerstenmayer | Analysis of a degenerate parabolic cross-diffusion system for ion transport[END_REF][START_REF] Cancès | Finite-volume scheme for a degenerate cross-diffusion model motivated from ion transport[END_REF][START_REF] Gerstenmayer | Comparison of a finite-element and finite-volume scheme for a degenerate cross-diffusion system for ion transport[END_REF]].

We supplement system (1.1) with no-flux boundary conditions

(1.3) N j=1
a i,j (u j ∇u i -u i ∇u j ) • n = 0 on (0, T ) × ∂Ω, i ∈ [ [1, N ]].

The initial concentration U 0 = (u 0 1 , . . . , u 0 N ) is supposed to be measurable and to map Ω into

A = U = (u 1 , . . . , u N ) ∈ R N + N i=1 u i = 1 ,
so we write in the condensed form U 0 ∈ L ∞ (Ω; A), which means that U 0 is measurable and takes its values in the bounded subset A of R N . Finally, we assume that all the chemical species under consideration are present:

(1.4)

Ω u 0 i dx > 0, ∀i ∈ [[1, N ]].
1.2. Formal structure. This system has several structural properties, the goal of this subsection is to exhibit them. The calculations presented in this section are formal: we assume that the solutions to (1.1) enjoy enough regularity to justify the calculations below. Rigorous proofs at the continuous level for the system under consideration here can be found in [START_REF] Bakhta | Cross-diffusion systems with non-zero flux and moving boundary conditions[END_REF][START_REF] Berendsen | Strong solutions and weak-strong stability in a system of cross-diffusion equations[END_REF] (see also [START_REF] Gerstenmayer | Analysis of a degenerate parabolic cross-diffusion system for ion transport[END_REF]). The properties listed here can also be obtained by passing to the limit in the numerical scheme. The first property we point out is the conservation of mass for all the species involved in System (1.1). Proof. Let U be a solution of (1.1), t ∈ [0, T ], i ∈ [[1, N ]], and let ϕ(x, s) = 1 [0,t] (s). With this particular choice of ϕ, we have for all s that

Ω div   N j=1 a i,j (u j ∇u i -u i ∇u j )   ϕ(x, s)dx = - Ω N j=1
a i,j (u j ∇u i -u i ∇u j ) ∇ϕ(x, s)dx = 0.

Hence, using ϕ as a test function in (1.1), we have:

t 0 d ds Ω u i (x, s)dx ds = 0.
The fundamental theorem of calculus yields the desired lemma.

The symmetry of the matrix A = (a i,j ) yields: N i=1 N j=1 a i,j (u j ∇u i -u i ∇u j ) = 0.

Therefore, a solution U to (1.1) satisfies ∂ t N i=1 u i = 0. Admit that u i (t, x) ≥ 0 for all t > 0 (this will be proved in the discrete setting and is proved in [START_REF] Berendsen | Strong solutions and weak-strong stability in a system of cross-diffusion equations[END_REF]Proposition 2.2] in the continuous setting), then the admissibility condition encoded in A is preserved along time.

Lemma 1.2. Let U be a solution to (1.1) and (1.3) corresponding to an initial data U 0 ∈ L ∞ (Ω; A), then U (t, x) ∈ A for all (t, x) ∈ A, i.e., U ∈ L ∞ (Q T ; A).

The system can be derived by passing to the macroscopic limit from a random jump process in the spirit of [START_REF] Giacomin | Phase segregation dynamics in particle systems with long range interactions. I. Macroscopic limits[END_REF][START_REF] Burger | Nonlinear cross-diffusion with size exclusion[END_REF]. As expected because of this derivation from statistical physic considerations, the system fulfills Onsager's reciprocal relation [START_REF] Onsager | Reciprocal relations in irreversible processes[END_REF][START_REF] Onsager | Reciprocal relations in irreversible processes[END_REF] and has a formal gradient flow structure. The driving functional is the mixing entropy

(1.5) E : L ∞ (Ω; A) → R, U → Ω N i=1 u i ln(u i )dx.
The next property we want to highlight at the continuous level is the decay of this entropy. Using the chain rule ∇c = c∇ ln(c), the system (1.1) is formally equivalent to

(1.6) ∂ t u i -div   N j=1 a i,j u i u j (∇ ln(u i ) -∇ ln(u j ))   = 0, i ∈ [[1, N ]].
Proposition 1.3. E is a Lyapunov functional for the system (1.3)-(1.6). More precisely, the following entropy -entropy dissipation estimate holds:

(1.7) d dt E(U ) + Ω   1≤i<j≤N a i,j u i u j |∇ ln(u i ) -∇ ln u j | 2   dx = 0.
Proof. First, we notice that thanks to the conservation of mass:

d dt E(U ) = d dt Ω N i=1 u i (ln(u i ) -1) = Ω N i=1 ln(u i )∂ t u i .
Then multiply Equation (1.6) by ln(u i ) and integrate by part in order to get:

Ω ln(u i )∂ t u i + Ω   N j=1 a i,j u i u j ∇ ln(u i ) • (∇ ln(u i ) -∇ ln(u j ))   = 0. Summing over i ∈ [[1, N ]
] yields the announced result thanks to the symmetry of A.

The entropy -entropy dissipation relation (1.7) is key in the analysis of many cross-diffusion systems, as exposed in [START_REF] Jüngel | The boundedness-by-entropy method for cross-diffusion systems[END_REF][START_REF] Jüngel | Entropy methods for diffusive partial differential equations[END_REF]. It will also play a central role in this paper. Assume that (1.8) min i =j a i,j > 0 as it will be done in Section 4. As a consequence of the inequality

N i=1 Ω |∇u i | 2 ≤ 4 N i=1 Ω |∇ √ u i | 2 ≤ 1 min i =j a i,j Ω 1≤i<j≤N a i,j u i u j |∇ ln(u i ) -∇ ln(u j )| 2 ,
we deduce from (1.7) a L 2 (0, T ; H 1 (Ω)) N estimate on U . This motivates the following notion of weak solution.

Definition 1.4. A weak solution U to (1.1) and (1.3) corresponding to the initial profile

U 0 ∈ L ∞ (Ω; A) is a function of L ∞ (Q T ; A) ∩ L 2 [0, T ]; H 1 (Ω) N satisfying, ∀i ∈ [[1, N ]], ∀ϕ ∈ C ∞ c ([0, T ) × Ω):
(1.9) -

Q T u i ∂ t ϕdxdt - Ω u 0 i ϕ(0, •)dx + Q T N j=1 a i,j (u j ∇u i -u i ∇u j ) ∇ϕ = 0.
The regularity requirement on a weak solution U is natural in the setting where Assumption (1.8) holds. In this case, the solution even enjoys a stronger regularity, as established in the recent contribution [START_REF] Berendsen | Strong solutions and weak-strong stability in a system of cross-diffusion equations[END_REF]. In the case where (1.8) is not fulfilled (but under a structural assumption on the matrix A), a more involved notion of weak solution has to be introduced, cf. [START_REF] Gerstenmayer | Analysis of a degenerate parabolic cross-diffusion system for ion transport[END_REF].

There is an important property that relates the model (1.1) to classical Fickian diffusion. As a consequence of Lemma 1.2, one can rewrite

(1.10) div   N j=1 (u j ∇u i -u i ∇u j )   = ∆ u i , i ∈ [[1, N ]].
As a consequence, if all the a i,j are equal to some a ∈ R, then the system (1.1) reduces to N uncoupled heat equations ∂ t u i = a∆u i . Based on the identity (1.10), we can rewrite the system (1.1) under the form

(1.11) ∂ t u i -a ∆ u i -div   N j=1 (a i,j -a ) (u j ∇u i -u i ∇u j )   = 0, i ∈ [[1, N ]],
where a ∈ R is arbitrary for the moment. The formulation (1.11) is at the basis of our discretization.

1.3. Objectives. The goal of this paper is to build and analyze a numerical scheme preserving the properties discussed in the previous section, namely:

• the non-negativity of the concentrations;

• the conservation of mass (Lemma 1.1);

• the preservation of the volume filling constraint (Lemma 1.2);

• the entropy-entropy dissipation relation (Proposition 1.3). The construction of our scheme is the purpose of Section 2. In Section 3, we will show the existence of solutions to this scheme and the preservation of discrete counterparts to the previously listed physical properties. Section 4 is devoted to the convergence of the numerical scheme toward weak solutions provided Assumption (1.8) is satisfied. Finally, in Section 5, we show the outcomes of some numerical experiments.

Before entering the core of the paper, let us mention that the development of numerical analysis for cross-diffusion systems is quite recent. To our knowledge, the first convergence study of a finite volume approximation for a non-degenerate cross-diffusion problem was carried out in [START_REF] Andreianov | Analysis of a finite volume method for a cross-diffusion model in population dynamics[END_REF]. This contribution is based on classical quadratic energy estimate, similarly to what is proposed in [START_REF] Zhou | A second order numerical scheme for the annealing of metal-intermetallic laminate composite: A ternary reaction system[END_REF]. The implementation of the discrete entropy method [START_REF] Chainais-Hillairet | Entropy method and asymptotic behaviours of finite volume schemes[END_REF] for cross-diffusion systems is more recent. Let us cite [START_REF] Ait | Numerical analysis of a finite volume scheme for a seawater intrusion model with cross-diffusion in an unconfined aquifer[END_REF][START_REF] Hammou Oulhaj | Convergence of a positive nonlinear control volume finite element scheme for an anisotropic seawater intrusion model with sharp interfaces[END_REF] where upstream mobility finite volume and control volume finite element schemes for a multiphase extension of the porous medium equation are studied. Upwinding is also used in [START_REF] Cancès | Finite-volume scheme for a degenerate cross-diffusion model motivated from ion transport[END_REF] to approximate the solution of a system which is very close to the problem (1.1) under study, or in [START_REF] Carrillo | Convergence of a finite volume scheme for a system of interacting species with cross-diffusion[END_REF] for a problem in which nonlocal interactions are also considered. As a consequence of the upwind choice for the mobility, the schemes presented in [START_REF] Ait | Numerical analysis of a finite volume scheme for a seawater intrusion model with cross-diffusion in an unconfined aquifer[END_REF][START_REF] Hammou Oulhaj | Convergence of a positive nonlinear control volume finite element scheme for an anisotropic seawater intrusion model with sharp interfaces[END_REF][START_REF] Cancès | Finite-volume scheme for a degenerate cross-diffusion model motivated from ion transport[END_REF] and [START_REF] Carrillo | Convergence of a finite volume scheme for a system of interacting species with cross-diffusion[END_REF] are first-order accurate in space. A natural solution to pass to order two is to rather consider mobilities given by arithmetic means [START_REF] Daus | Convergence of a finite-volume scheme for a degenerate-singular cross-diffusion system for biofilms[END_REF]. The motivation of the finite element scheme proposed in [START_REF] Jüngel | Convergence of an implicit Euler Galerkin scheme for Poisson-Maxwell-Stefan systems[END_REF] is also the same. However, the scheme proposed in [START_REF] Jüngel | Convergence of an implicit Euler Galerkin scheme for Poisson-Maxwell-Stefan systems[END_REF] is expressed in entropy (or dual) variables (in our context log(u i )) leading to computational difficulties when the concentrations are close to 0. Other entropy stable numerical schemes have been proposed for cross-diffusion systems, as for instance discontinuous Galerkin schemes in [START_REF] Sun | An entropy stable high-order discontinuous Galerkin method for cross-diffusion gradient flow systems[END_REF], or finite volumes on staggered cartesian grids for Maxwell-Stefan cross-diffusion in [START_REF] Huo | An energy stable and positivity-preserving scheme for the Maxwell-Stefan diffusion system[END_REF]. Finally, let us point out that the design of entropy (or energy) stable numerical schemes for dissipative systems with formal gradient flow structure in a Riemannian geometry have been the purpose of intense research in the recent years, as shows the extensive (but not exhaustive) recent literature [START_REF] Cancès | Numerical analysis of a robust free energy diminishing finite volume scheme for parabolic equations with gradient structure[END_REF][START_REF] Kou | Linearly decoupled energy-stable numerical methods for multicomponent two-phase compressible flow[END_REF][START_REF] Cancès | Convergence and a posteriori error analysis for energystable finite element approximations of degenerate parabolic equations[END_REF][START_REF] Almeida | Energy and implicit discretization of the Fokker-Planck and Keller-Segel type equations[END_REF][START_REF] Murphy | Control volume approximation of degenerate two-phase porous flows[END_REF][START_REF] Cancès | A numerical analysis focused comparison of several finite volume schemes for an unipolar degenerated drift-diffusion model[END_REF][START_REF] Cancès | Finite volume approximation of a two-phase two fluxes degenerate cahn-hilliard model[END_REF][START_REF] Gu | Bound preserving and energy dissipative schemes for porous medium equation[END_REF][START_REF] Bailo | Convergence of a fully discrete and energy-dissipating finite-volume scheme for aggregation-diffusion equations[END_REF][START_REF] Carrillo | Lagrangian schemes for Wasserstein gradient flows[END_REF] on this topic. Let us also refer to [START_REF] Chen | Positivity-preserving, energy stable numerical schemes for the Cahn-Hilliard equation with logarithmic potential[END_REF][START_REF] Dong | A positivity-preserving, energy stable and convergent numerical scheme for the Cahn-Hilliard equation with a Flory-Huggins-deGennes energy[END_REF][START_REF] Shen | A new class of efficient and robust energy stable schemes for gradient flows[END_REF][START_REF] Both | The gradient flow structures of thermo-poro-visco-elastic processes in porous media[END_REF][START_REF] Shen | Unconditionally bound preserving and energy dissipative schemes for a class of Keller-Segel equations[END_REF] for the simpler situation of gradient flows in Hilbert spaces.

Remark 1.1. Our study can be extended to the case where reaction terms are incorporated in the system. More precisely, one can consider a system of the form

(1.12) ∂ t U -div(D(U )∇U ) = R(U ),
where D(U ) is as in (1.2), and where the reaction term function

R : R N → R N U → R(U ) = (r i (U )) 1≤i≤N
is continuous and satisfies the following structural properties which are classically satisfied for reactive systems: (i) Isochore process:

N i=1 r i (U ) = 0 for all U ∈ R N ; (ii) Positivity preservation: r i (U ) ≥ 0 for U ∈ R N with u i ≤ 0; (iii) Entropy dissipation: there exists U = (u i ) 1≤i≤N in A with u i > 0 for all i ∈ {1, ,N } such that (1.13) R(U ) • ln(U/U ) = N i=1 r i (U ) (ln(u i ) -ln(u i )) ≤ 0, ∀U ∈ A.
Because of reaction terms, the volume of each specie is no longer conserved, so that Lemma 1.1 does no longer hold true. However, because of Assumption (i) above, the total volume is conserved, hence the condition N i=1 u i (t, x) = 1 remains true for all time. Since Assumption (ii) guarantees the positivity of the solution, one gets that U (t, x) belongs to L ∞ (Q T , A) with the reaction term as well. Finally, Assumption (iii) on the reaction terms ensures that the relative entropy

(1.14) E(U |U ) = N i=1 Ω u i log u i u i ≥ 0
is a Lyapunov functional for the system. This stability property allows to extend our purpose in presence of reaction. Note that in absence of reaction R ≡ 0, this relative entropy (one can for instance set u i = Ω u 0 i ) coincides with the mixing entropy (1.5) up to an additive constant thanks to the conservation of the volume of each specie, cf. Lemma 1.1. Finally, let us note that if U ∈ A is such that u i = 0 for some i ∈ {1, . . . , N }, the relative entropy E(U |U ) is no longer well-defined. Our analysis can still be extended by showing that the mixing entropy E(U ) grows at most linearly with time, which is sufficient for establishing the convergence of the straightforward extension to the case R = 0 of the finite volume scheme to be presented in the next section.

2. Finite Volume approximation. This section is organized as follows. First, in Section 2.1, we state the requirements on the mesh and fix some notations. Then in Section 2.2, we describe the numerical scheme to be studied in this paper. It is based on Formulation (1.11) of the problem. Then in Section 2.3, we state our two main results. The first one, namely Theorem 2.2, focuses on the case of a fixed mesh. We are interested in the existence of a solution to the nonlinear system corresponding to the scheme, and the dissipation of the entropy at the discrete level. More precisely, one establishes that the studied scheme satisfies a discrete entropy -entropy dissipation inequality that can be thought of as a counterpart to Proposition 1.3. Our second main result, namely Theorem 2.3, is devoted to the convergence of the scheme towards a weak solution as the time step and the mesh size tend to 0.

2.1. Discretization of (0, T ) × Ω. The scheme we propose relies on two-point flux approximation (TPFA) finite volumes. As explained in [START_REF] Droniou | Finite volume schemes for diffusion equations: introduction to and review of modern methods[END_REF][START_REF] Eymard | TP or not TP, that is the question[END_REF][START_REF] Gärtner | Why Do We Need Voronoi Cells and Delaunay Meshes[END_REF], this approach appears to be very efficient as soon as the continuous problem to be solved numerically is isotropic and one has the freedom to choose a suitable mesh fulfilling the so-called orthogonality condition [START_REF] Herbin | An error estimate for a finite volume scheme for a diffusion-convection problem on a triangular mesh[END_REF][START_REF] Eymard | Finite volume methods[END_REF]. We recall here the definition of such a mesh, which is illustrated in Figure 1.

Definition 2.1. An admissible mesh of Ω is a triplet T , E, (x K ) K∈T such that the following conditions are fulfilled.

(i) Each control volume (or cell) K ∈ T is non-empty, open, polyhedral and convex. We assume that

K ∩ L = ∅ if K, L ∈ T with K = L, while K∈T K = Ω. (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 satisfy 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

x σ ∈ σ such that x σ -x K is orthogonal to σ.
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 we also define E K,int = E K ∩ E int the set of its faces that belong to E int . For such a face σ ∈ E K,int , we may write σ = K|L, meaning that σ = K ∩ L, where L ∈ T .

E ext = {σ ∈ E | σ ⊂ ∂Ω}. For a given control volume K ∈ T , σ = K | L K x σ x K x L
Given σ ∈ E, we let

d σ = |x K -x L | if σ = K|L ∈ E int , |x K -x σ | if σ ∈ E ext , and
τ σ = m σ d σ .
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

h T = max K∈T diam(K), ζ T = min K∈T min σ∈E K d(x K , σ) d σ .
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 n ∈ {1, . . . , N T }, by ∆t = (∆t n ) 1≤n≤N T , and by h T = max 1≤n≤N T ∆t n . In what follows, we will use boldface notations for mesh-indexed families, typically for elements

of R T , (R T ) N , (R T ) N T , or even (R T ) N ×N T . 2.2. Numerical scheme. The initial data U 0 ∈ L ∞ (Ω; A) is discretized into U 0 = u 0 i i∈[[1,N ]] ∈ (R T ) N = u 0 i,K K∈T ,i∈[[1,N ]]
by setting

(2.1)

u 0 i,K = 1 m K K u 0 i (x)dx, ∀K ∈ T , i ∈ [[1, N ]]. Assume that U n-1 = u n-1 i,K K∈T ,i∈[[1,N ]]
is given for some n ≥ 1, then we have to define how to compute

U n = u n i,K K∈T ,i∈[[1,N ]]
. First, we introduce some notations. Given any discrete scalar field c = (c K ) K∈T ∈ R T , we define for all cell K ∈ T and interface σ ∈ E K the mirror value c Kσ of c K across σ by setting:

(2.2) c Kσ = c L if σ = K|L ∈ E int , c K if σ ∈ E ext .
We also define the oriented and absolute jumps of c across any edge by

D Kσ c = c Kσ -c K , D σ c = |D Kσ c|, ∀K ∈ T , ∀σ ∈ E K .
The scheme is based on the formulation (1.11). It requires the introduction of a parameter a on which we only have the following requirements:

(2.3)
a > 0 and a ≥ min

(i,j)
a i,j .

The conservation laws are discretized in a conservative way with a time discretization relying on the backward Euler scheme:

(2.4a) m K u n i,K -u n-1 i,K ∆t n + σ∈E K F n i,Kσ = 0, ∀K ∈ T , ∀i ∈ [[1, N ]].
The discrete fluxes are computed thanks to a formula based on (1.11) and on TPFA finite volumes:

(2.4b)

F n i,Kσ = -a τ σ D Kσ u n i -τ σ   N j=1 (a i,j -a ) u n j,σ D Kσ u n i -u n i,σ D Kσ u n j   , for all K ∈ T , σ ∈ E K and i ∈ [[1, N ]].
Edge values u n j,σ j of the concentrations u j appears in Formula (2.4b). It is deduced from u n j,K and u n j,Kσ thanks to a logarithmic mean, i.e., (2.4c)

u n j,σ =          0 if min(u n j,K , u n j,Kσ ) ≤ 0, u n j,K if 0 ≤ u n j,K = u n j,Kσ , u n j,K -u n j,Kσ ln(u n j,K ) -ln(u n j,Kσ )
otherwise.

This choice for the edge concentration is crucial for the preservation at the discrete level of a discrete entropy -entropy dissipation inequality similar to the one highlighted in Proposition 1.3. Equations (2.4b) and (2.2) implies that for all σ ∈ E ext : F n i,Kσ = 0, so that the no-flux boundary condition (1.3) is taken into account.

Remark 2.1. Let us highlight why the choice of a strictly positive a is important. Consider a mesh with two cells K, L, and one edge. We consider two species and let u 0 K = (0, 1) and u 0 L = (1, 0). We have: u 0 1,K|L = 0 and u 0 2,K|L = 0, hence, if a = 0, the initial condition is a stationary solution even though this is not expected for a discretization of the heat equation. Setting a > 0 eliminates these spurious solutions. The choice of a has a strong influence on the numerical outcomes, as it will be shown in Section 5, but we don't have a clear understanding yet on the methodology to choose an optimal a . What seems clear is that a has to be chosen in the interval [min i =j a i,j , max i =j a i,j ]. A tentative non-optimal formula is proposed in Section 5.

Remark 2.2. The time discretization in scheme (2.4) is only first-order accurate since it relies on the backward Euler approximation. Going to second-order time discretizations is tempting, but no theoretical guarantees concerning the entropy stability of the scheme can be granted then. This is due to the fact that the entropy (1.5) is not quadratic, hence neither the Crank-Nicolson scheme nor the BDF2 scheme can be shown to be unconditionally stable here. This lack of theoretical foundation for the entropy stability has for instance also been reported in [START_REF] Gu | Bound preserving and energy dissipative schemes for porous medium equation[END_REF].

Main results and organization.

The first theorem proven is this paper concerns the existence of discrete solutions for a given mesh, and the preservation of the structural properties listed in Section 1.3:

• the mass of each specie is conserved along the time steps;

• the concentrations are (strictly) positive and sum to 1 in all the cells, i.e., U n K ∈ A for all K ∈ T and n ≥ 1; • the discrete counterpart of the entropy decays along time. For this last property, we need to introduce the discrete entropy functional E T , which is defined by:

(2.5) E T (U ) = K∈T N i=1 m K u i,K ln u i,K , ∀U = (u i,K ) K∈T ,i∈[[1,N ]] ∈ A T .
As stated in Theorem 2.2 below, the nonlinear system corresponding to our scheme (2.4) admits solutions that preserve the physical bounds on the concentrations and the decay of the entropy.

Theorem 2.2. Let (T , E, (x K ) K∈T ) be an admissible mesh and let U 0 be defined by (2.1). Then, for all 1 ≤ n ≤ N T , the nonlinear system of equations (2.2) -(2.4), has a positive solution U n ∈ A T . Moreover, such a solution satisfies

E T (U n ) ≤ E T (U n-1 ) for all n ∈ [[1, N T ]], K∈T m K u n i,K = Ω u 0 i for all i ∈ [[1, N ]] and n ∈ [[0, N T ]].
The proof of Theorem 2.2 will be the purpose of Section 3. With a discrete solution (U n ) 1≤n≤N T to the scheme (2.4) at hand, we can define the piecewise constant approximate solution

U T ,∆t = (u i,T ,∆t ) i∈[[1,N ]] : Q T → A defined almost everywhere by U T ,∆t (t, x) = U n K if (t, x) ∈ (t n-1 , t n ] × K.
This definition will be developed in Section 4 and supplemented by other reconstruction operators. Let T m , E m , (x K ) K∈Tm m≥1 be a sequence of admissible discretizations with h Tm , h T,m tending to 0 as m tends to +∞, while the regularity ζ Tm remains uniformly bounded from below by a positive constant ζ . Thanks to Theorem 2.2, we dispose of a family U m of solutions to our scheme. The convergence of U m is the purpose of Theorem 2.3 whose proof is detailed in Section 4.

Theorem 2.3. Assume that the nondegeneracy assumption (1.8) holds. Given any sequence of solutions

U m = u n i,K i∈[[1,N ]],K∈Tm,1≤n≤N T ,m
, there exists at least one

U ∈ L ∞ (Q T ; A) ∩ L 2 (0, T ); H 1 (Ω) such that, up to a subsequence, (2.6) U Tm,∆tm -→ m→∞ U strongly in L p (Q T ), for any 1 ≤ p < ∞,
Moreover, U is a weak solution in the sense of Definition 1.4.

Remark 2.3. Theorem 2.3 above establishes the convergence of the scheme, but no error estimate can be deduced from its proof. Indeed, its proof can be thought of as an adaptation to the discrete setting of an existence proof based on compactness arguments, as the for instance proposed in [START_REF] Jüngel | The boundedness-by-entropy method for cross-diffusion systems[END_REF]. The derivation of error estimates is different since it relies on the perturbation of uniqueness proofs. As far as we know, uniqueness for the problem under consideration is an open question excepted in the one-dimensional setting [START_REF] Berendsen | Strong solutions and weak-strong stability in a system of cross-diffusion equations[END_REF], but a natural approach would be to derive error estimates based on the use of the relative entropy, as done for instance in [START_REF] Cancès | Error estimate for time-explicit finite volume approximation of strong solutions to systems of conservation laws[END_REF] for hyperbolic systems or in [START_REF] Gallouët | Error estimates for a numerical approximation to the compressible barotropic Navier-Stokes equations[END_REF][START_REF] Gallouët | Error estimates for the implicit MAC scheme for the compressible Navier-Stokes equations[END_REF] for the compressible Navier-Stokes problem. However, the aforementioned strategy generally leads to under-optimal error estimates. The recovery of optimal error estimates for finite volume approximation of diffusion equations on unstructured grids has only been achieved recently [START_REF] Droniou | Improved L 2 estimate for gradient schemes and super-convergence of the TPFA finite volume scheme[END_REF]. The extension of this optimal result to our much more complex cross-diffusion system appears to be an interesting and challenging issue.

3. Numerical analysis on a fixed mesh. This section is devoted to the proof of Theorem 2.2. In Section 3.1, we establish a priori estimates on a slightly modified scheme that will be shown to reduce to the original scheme (2.4). Then in Section 3.2, we apply a topological degree argument to prove the existence of solutions to our scheme. Section 3.3 is devoted to the proof of the entropy dissipation property.

To prove the existence of solutions to the system of equations (2.4), we need the inequality i u i,σ ≤ 1. We then slightly modify (2.4) by adding the following equation:

u n i,σ = u n i,σ max(1, N j=1 u n j,σ )
, and replacing u n i,σ by u n i,σ in (2.4b). We will denote this new system (S) and see in Proposition 3.2 that its solutions satisfy i u i,σ ≤ 1, so that u n i,σ = u n i,σ . Whence they also satisfy the original system of equations.

A priori estimates. The first lemma shows the nonnegativity of the solutions to (S).

Lemma 3.1. Given a nonnegative U n-1 , any solution U n to (S) is also nonnegative.

Proof. Let U n be a solution of (S) and let i ∈ [ [1, N ]]. We consider a cell K ∈ T where u n i reaches its minimum, i.e., u n i,K ≤ u n i,L for all L ∈ T , and assume for contradiction that u n i,K is (strictly) negative. Equation (2.4b) then gives:

m K u n i,K -u n-1 i,K ∆t n = - σ∈E K F n Kσ .
The term on the left hand side is negative since u n-1 i,K ≥ 0 > u n i,K , whereas the right-hand side may be simplified noticing that u n i,σ = 0:

σ∈E K a τ σ D Kσ u n i + σ∈E K τ σ N j=1 (a i,j -a ) u n j,σ D Kσ u n i = - σ∈E K F n Kσ < 0.
Noticing that D Kσ u n i ≥ 0, u n j,σ ≥ 0, and

N j=1 u n j,σ ≤ 1 we obtain that 0 ≤ σ∈E K a (1 - N j=1 u n j,σ )τ σ D Kσ u n i < 0,
which is absurd, hence the desired result.

Let us now show that the concentrations sum to 1 in all the cells. Lemma 3.2. Given U n-1 in A T , any solution U n to (S) is also in A T .

Proof. Thanks to Lemma 3.1, it suffices to show that N i=1 u n i,K = 1 for all K ∈ T . Let U n be a solution to (S). Using (2.4b) in (2.4a) and summing over the species leads to:

N i=1 u n i,K - N i=1 u n-1 i,K ∆t n m K -a σ∈E K τ σ D Kσ i u i - σ∈E K τ σ i   N j=1 (a i,j -a ) u n j,σ D Kσ u i -u n i,σ D Kσ u j   = 0, ∀K ∈ T .
The third term of the left-hand side vanishes thanks to the symmetry of A, so that

N i=1 u n i,K - N i=1 u n-1 i,K ∆t n m K -a σ∈E K τ σ D Kσ i u i = 0, ∀K ∈ T .
The discrete quantity i u i is solution to the classical backward Euler TPFA scheme for the heat equation, which is well posed. So i u n i = i u n-1 i = 1 is its unique solution, hence the desired result.

Existence of solutions.

Using the tools exposed in the previous subsection, we may derive the existence of a solution to (S): Proposition 3.1. Given U n-1 in A T , there exists at least one solution to (S) in A T .

Proof. The proof relies on a topological degree argument [START_REF] Leray | Topologie et équations fonctionnelles[END_REF][START_REF] Deimling | Nonlinear functional analysis[END_REF]. The idea is to transform continuously our complex nonlinear system into a linear system while guaranteeing that the a priori estimates controlling the solution remain valid all along the homotopy. We sketch the main ideas of the proof, making the homotopy explicit. We are interested in the existence of zeros for a functional

H : [0, 1] × (R N ) T → (R N ) T (λ, U ) → H(λ, U )
that boils down to the scheme (S) when λ = 1. In our case, we set:

(3.1) H(λ, U ) i,K = u i,K -u n-1 i,K ∆t n m K -a σ∈E K τ σ D Kσ u i -λ σ∈E K τ σ   N j=1 (a i,j -a ) ( u j,σ D Kσ u i -u i,σ D Kσ u j )   , ∀K ∈ T , ∀i ∈ [[1, N ]].
One notices that H(0, U ) = 0 is the classical heat equation, the solution of which belongs to A T . Therefore, fixing η > 0, the relatively compact open set

A T η = U ∈ R T inf V ∈A T U -V < η
has a topological degree equal to 1. Note that the choice of the norm in the definition of A T η is not important since the dimension is finite. Moreover, thanks to Lemma 3.2, the solutions u (λ) of H(λ, U ) = 0 remains in A T , thus in the interior of A T η . Thus the topological degree of A T η for λ = 1 is still equal to 1, hence the existence of (at least) one solutions to (S). Since η > 0 is arbitrary, then there is a solution in

A T = η>0 A T η .
To prove the Theorem 2.2, we need to transfer this existence result on the original system. Proposition 3.2. A solution U n of (S) is a solution of (2.4). Reciprocally, a solution of (2.4) in A T is a solution of (S).

Proof. Let U n be a solution of (S). A simple convexity argument shows that the logarithmic mean of two nonnegative number is smaller than the arithmetic mean, so that u n i,σ ≤

u n i,K +u n i,Kσ 2
. Summing w.r.t. i ∈ [[1, N ]] and using that the solution U of (S) belongs to A T , one gets that i u n i,σ ≤ 1 for all σ ∈ E. Therefore u n i,σ = u n i,σ and U n is a also solution to (2.4). The proof of the reverse implication follows the same lines.

Entropy dissipation.

We intend here to prove a discrete counterpart to Proposition 1.3. The proof will be very similar and requires a discrete counterpart of the conservation of mass (Lemma 1.1). Lemma 3.3. Given any U n-1 ∈ A T , any solution U n to (2.4) satisfies:

K∈T m K u n i,K = K∈T m K u n-1 i,K = Ω u 0 i dx, ∀i ∈ [[1, N ]].
The proof of this lemma is a straightforward calculation based on equation (2.4a), the conservativity of the fluxes, and the definition (2.1) of the discrete initial condition. With this lemma and Proposition 3.2, we can refine the result Lemma 3.1 to get the strict positivity of any solution to (2.4) belonging to A T .

Lemma 3.4. Let U n-1 ∈ A T be such that K m K u n-1 i,K > 0 for all i ∈ [[1, N ]], then any solution to (2.4) in A T is positive: u n i,K > 0 for all i ∈ [[1, N ]
] and all K ∈ T . Proof. Let U n ∈ A T be a solution to the scheme (2.4), and let i ∈ [ [1, N ]]. We know from Lemma 3.1 that u n i ≥ 0. Assume for contradiction that there exists one cell K such that u 0 i,K vanishes. Using Lemma 3.3 and the connectivity of Ω, there exists σ = K|L ∈ E int such that u n i,K = 0 and u n i,L > 0. Then u n i,σ = 0 and as in the proof of Lemma 3.1:

a (1 - N j=1 u n j,σ )τ σ D Kσ u n i ≤ 0. Using u n j,σ ≤ u n j,K +u n j,L 2
and u i,σ = 0 we deduce that

N j=1 u n j,σ ≤ N j =i u n j,K + u n j,L 2 ≤ 1 - u n i,L 2 < 1.
Therefore a (1 -N j=1 u n j,σ )τ σ > 0, and since D Kσ u n i > 0, we deduce that:

0 < a (1 - N j=1 u n j,σ )τ σ D Kσ u n i ≤ 0.
As this statement is absurd, our assumption was false, hence the desired result.

As in the continuous case, we will use the conservation of mass (Lemma 3.3) and a discrete equivalent of the chain rule ∇c = c∇ ln c. This equivalently writes

(3.2) D Kσ u n i = u n i,σ D Kσ ln(u n i ), ∀i ∈ [[1, N ]], ∀K ∈ T .
The above discrete chain rule follows from the definition (2.4c) of u n i,σ and the positivity of solutions to (2.4) which gives a sense to ln(u n i ). Using (3.2) in (2.4b), U n satisfies (3.3)

u n i,K -u n-1 i,K ∆t n m K - σ∈E K τ σ   N j=1 (a i,j -a )u n i,σ u n j,σ (D Kσ ln(u i ) -D Kσ ln(u j ))   -a σ∈E K τ σ D Kσ u n i = 0, ∀K ∈ T , ∀i ∈ [[1, N ]].
This reformulation is suitable for proving a discrete entropy -entropy dissipation inequality, which should be seen as a discrete counterpart of Proposition 1.3.

Proposition 3.3. Given U n-1 in A T , any solution U n ∈ A T to (2.4) satisfies (3.4) E T (U n ) -E T (U n-1 ) + ∆t n min 1≤i,j≤N a i,j σ∈E N i=1 τ σ u n i,σ (D Kσ ln(u n i )) 2 ≤ 0. In particular, E T (U n ) ≤ E T (U n-1 ).
Proof. Multiplying equation (3.3) by ∆t n ln(u n i,K ) and summing over the cells and species leads to:

(3.5) K∈T N i=1 (u n i,K ln(u n i,K ) -u n-1 i,K ln(u n i,K ))m K + ∆t n a σ∈E N i=1 τ σ u n i,σ (D Kσ ln(u n i )) 2 -∆t n K∈T N i=1 σ∈E K τ σ   N j=1 (a i,j -a )u n j,σ u n i,σ ln(u n i,K )D Kσ ln(u n i ) -ln(u n j )   = 0.
Using the symmetry of the matrix A and discrete integration by part, both in space and with respect to the species, we have:

(3.6) K∈T N i=1 σ∈E K τ σ   N j=1 (a i,j -a )u n j,σ u n i,σ ln(u n i,K )D Kσ ln(u n i ) -ln(u n j )   = - σ∈E τ σ   1≤i<j≤N ) (a i,j -a )u n j,σ u n i,σ D Kσ ln(u n i ) -ln(u n j ) 2   .
On the other hand, the convexity of c ln(c) yields:

u n i,K -u n-1 i,K + u n i,K ln(u n i,K ) -u n-1 i,K ln(u n i,K ) ≥ u n i,K ln(u n i,K ) -u n-1 i,K ln(u n-1 i,K ).
Combining this inequality with Equation (3.6) and Lemma 3.3 in (3.5) provides:

E T (U n ) -E T (U n-1 ) + ∆t n a σ∈E N i=1 τ σ u n i,σ (D Kσ ln(u n i )) 2 + ∆t n σ∈E τ σ   1≤i<j≤N (a i,j -a )u n j,σ u n i,σ D Kσ ln(u n i ) -ln(u n j ) 2   ≤ 0.
Using the hypothesis 0 ≤ min a i,j ≤ a together with

(3.7) N i=1 u n i,σ (D Kσ ln(u n i )) 2 -   1≤i<j≤N u n j,σ u n i,σ D Kσ ln(u n i ) -ln(u n j ) 2   = N i=1 u n i,σ   1 - N j=1 u n j,σ   (D Kσ ln(u n i )) 2 ≥ 0,
we deduce that (3.4) holds.

The proof of Theorem 2.2 is now complete.

4. Convergence analysis. The goal of this Section is to prove Theorem 2.3, which states the convergence of the approximate solution towards a weak solution to the continuous problem in the sense of Definition 1.4 under the nondegeneracy condition (1.8). We could extend this result on several other special cases including the one treated in [START_REF] Cancès | Finite-volume scheme for a degenerate cross-diffusion model motivated from ion transport[END_REF]. We hint that the optimal assumption would be that the zeros of the diffusion matrix form a cluster-graph. However, we stick to the study of the non-degenerate case for the sake of simplicity.

We consider here a sequence T m , E m , (x K ) K∈Tm m≥1 of admissible discretizations with h Tm , h T,m tending to 0 as m tends to +∞, while the regularity ζ Tm remains uniformly bounded from below by a positive constant ζ . Theorem 2.2 provides the existence of a family of discrete solutions U m = u n i,K i∈[[1,N ]],K∈Tm,1≤n≤Nm . To prove Theorem 2.3, we first establish in Section 4.2 some compactness properties on the family of piecewise constant approximate solutions U Tm,∆tm . Then we identify the limit as a weak solution in Section 4.3. In order to enlighten the notations, we remove the subscript m as soon as it is not necessary for understanding.

Reconstruction operators.

To carry out the convergence analysis, we introduce some reconstruction operators following the methodology proposed in [START_REF] Droniou | The Gradient Discretisation Method[END_REF]. The operators π T : R T → L ∞ (Ω) and π T ,∆t : R

T N T → L ∞ (Q T ) are defined respectively by π T f (x) = f K if x ∈ K, ∀f = (f K ) K∈T , and
π T ,∆t f (t, x) = f n K if (t, x) ∈ (t n-1 , t n ] × K, ∀f = (f n K ) K∈T ,1≤n≤N T .
These operators allow to pass from the discrete solution (U n ) 1≤n≤N T to the approximate solution since

u i,T ,∆t = π T ,∆t (u n i ) n , ∀i ∈ [[1, N ]].
In order to carry out the analysis, we further need to introduce approximate gradient reconstruction. For σ = K|L ∈ E int , we denote by ∆ σ the diamond cell corresponding to σ, which 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 → L 2 (Ω) d we use in the analysis is merely weakly consistent (unless d = 1) and takes its source in [START_REF] Chainais-Hillairet | Finite volume scheme for multi-dimensional drift-diffusion equations and convergence analysis[END_REF][START_REF] Eymard | H-convergence and numerical schemes for elliptic problems[END_REF]. It is piecewise constant on the diamond cells ∆ σ , and it is defined as follows:

∇ T f (x) = d D Kσ f d σ n Kσ if x ∈ ∆ σ , ∀f ∈ R T ,
where n Kσ is the outer-pointing normal of K at σ. We also define ∇ T ,∆t : R

T ×N T → L 2 (Q T ) d by setting ∇ T ,∆t f (t, •) = ∇ T f n if t ∈ (t n-1 , t n ], ∀f = (f n ) 1≤n≤N T ∈ R T ×N T .
It follows from the definition of the approximate gradient that (4.1)

σ∈E τ σ D Kσ f D Kσ g = 1 d Ω ∇ T f • ∇ T gdx, ∀f , g ∈ R T .
This implies in particular that (

σ∈E τ σ |D σ f | 2 = 1 d Ω |∇ T f | 2 dx, ∀f ∈ R T . 4.2)

Compactness properties.

In this subsection, we take advantage of Proposition 3.3 and of the non-degeneracy assumption (1.8) to get enough compactness for the convergence. Lemma 4.1. There exists C depending only on Ω and min i =j a i,j such that

N i=1 Q T |∇ Tm,∆tm √ u i,m | 2 + π Tm,∆tm √ u i,m 2 dxdt ≤ C, ∀m ≥ 1.
Proof. We get rid of the subscript m for the ease of reading. The L ∞ bound on U yields immediately the L 2 estimate on π T ,∆t √ u i . The proof thus consists in proving the bound on the discrete gradient. Le us focus on the proof of

Q T |∇ T ,∆t √ u i | 2 dxdt ≤ C for some fixed i ∈ [[1, N ]]. Thanks to (4.
2), we have

Q T |∇ T ,∆t √ u i | 2 = d N T n=1 ∆t n σ∈Eint τ σ |D σ u n i | 2 , = d N T n=1 ∆t n σ∈Eint τ σ ǔn iσ |D σ ln(u n i )| 2 ,
where ǔn iσ = 4 (Dσ √

u n i ) 2 (Dσ ln(u n i ))
2 . It results from Cauchy-Schwarz inequality that

4 √ a - √ b 2 ≤ (a -b)(ln(a) -ln(b)), ∀(a, b) ∈ (0, +∞),
so that ǔn iσ ≤ u n iσ . Therefore, Proposition 3.3 provides:

min i =j a i,j N i=1 Q T |∇ T ,∆t √ u i | 2 ≤ d 4 E T (U 0 ) -E T (U N T ) .
As E T is bounded between -m Ω and 0 and as, by hypothesis, min a i,j > 0, we obtain the desired bound.

The inequality 2D σ u n i ≥ D σ u n i and Lemma 4.1 yield the following discrete L 2 (0, T ; H 1 (Ω)) estimate on u i . Corollary 4.2. There exists C depending only on Ω and min i =j a i,j such that

N i=1 Q T |∇ Tm,∆tm u i,m | 2 + (π Tm,∆tm u i,m ) 2 dxdt ≤ C, ∀m ≥ 1.
The following proposition is about the relative compactness of the approximate solution and of the weakly consistent approximate gradient. Proposition 4.1. Let (U m ) be the family of discrete solutions. There exists at least one

U ∈ L ∞ (Q T ; A) ∩ L 2 ((0, T ); H 1 (Ω)) such that, up to a subsequence, for all i ∈ [[1, N ]]: π Tm,∆tm u i,m -→ m→∞ u i strongly in L 2 (Q T ), (4.3)
∇ Tm,∆tm u i,m -→ m→∞ ∇u i weakly in L 2 (Q T ) d . (4.4)
Proof. We drop the subscript m for clarity. The proof of this result relies on a discrete Aubin-Lions lemma [START_REF] Gallouët | Compactness of discrete approximate solutions to parabolic PDEs-application to a turbulence model[END_REF]Lemma 3.4] on the particular setting of [START_REF] Cancès | Finite-volume scheme for a degenerate cross-diffusion model motivated from ion transport[END_REF]Lemma 9]. Define the discrete L 2 (0, T ; (H 1 (Ω)) ) norm by duality as follows:

v -1 = sup Ω π T vπ T ϕ, π T ϕ 2 L 2 + ∇ T ϕ 2 L 2 = 1 , ∀v ∈ R T . Therefore if ∇ T ,∆t u i L 2 (Q T ) ≤ C and n u n i -u n-1 i -1 ≤ C,
then, up to a subsequence, π T ,∆t u i tends towards some u i in L 2 (Q T ), while ∇ T ,∆t u i converges weakly towards ∇u i . In particular, U ∈ L 2 (0, T ; H 1 (Ω)) N . Corollary 4.2 provides the L 2 bound on ∇ T ,∆t u i . For the other inequality, we let

ϕ ∈ R T , n ∈ [[1, N T ]] and i ∈ [[1, N ]]. It follows from (2.4a) that Ω π T u n i -u n-1 i π T ϕ = -∆t n K∈T ϕ K σ∈E K F n i,Kσ .
Using (2.4b), this yields

1 ∆t n Ω π T u n i -u n-1 i π T ϕ = σ∈E a τ σ D Kσ u n i D Kσ ϕ + σ∈E τ σ   N j=1 (a i,j -a ) u n j,σ D Kσ u n i -u n i,σ D Kσ u n j   D Kσ ϕ.
Using the Cauchy-Schwarz inequality, the L ∞ bound on u n i,σ σ∈E,i∈[ [1,N ]] and (4.1) then leads to

1 ∆t n Ω π T u n i -u n-1 i π T ϕ ≤ a ∇ T u n i L 2 (Ω) ∇ T ϕ L 2 (Ω) + ∇ T ϕ L 2 (Ω) N j=1 |a i,j -a | ∇ T u n i L 2 (Ω) + ∇ T u n j L 2 (Ω) .
By definition of the discrete (H 1 (Ω)) norm, we have

u n i -u n-1 i ∆t n -1 ≤ a ∇ T u n i L 2 (Ω) + N j=1 |a i,j -a | ∇ T u n i L 2 (Ω) + ∇ T u n j L 2 (Ω) .
Using Corollary 4.2 again provides that

n u n i -u n-1 i -1 ≤ C.
The relative compactness properties on π T ,∆t u i and ∇ T ,∆t u i follow.

We still have to prove that U is in

L ∞ (Q T ; A). Let i ∈ [[1, N ]
] and let ϕ i ∈ L 2 (Q T ) be zero where the limit u i is nonnegative and 1 where the limit is negative, then

Q T ϕ i π T ,∆t u i ≥ 0 and Q T ϕ i π T ,∆t u i -→ m→+∞ Q T u i ϕ i ≤ 0.
Therefore, Q T u i ϕ i = 0, so that u i is nonnegative. Finally, the linearity of the limit yields N i=1 u i = 1. Remark 4.1. The uniform L ∞ (Q T ) bound on π Tm,∆tm U m together with the strong convergence in L 2 (Q T ) yield (2.6) thanks to Hölder's inequality:

π Tm,∆tm U m -→ m→∞ U strongly in L p (Q T ) N , for any 1 ≤ p < ∞.
We also need convergence properties for the face values u i,σ . We can reconstruct an approximate solution u i,E,∆t which is piecewise constant on the diamond cells by setting, for all i ∈ [[1, N ]]:

u i,E,∆t (t, x) = u n i,σ if (t, x) ∈ (t n-1 , t n ] × ∆ σ , σ ∈ E.
u i,Em,∆tm -→ m→∞ u i in L p (Q T ), for any 1 ≤ p < ∞,
where U is as in Proposition 4.1.

Proof. Here again, we get rid of m for clarity, and show the convergence for a specific value of p. The convergence for any finite p follows from the L ∞ (Q T ) bound on u i,Em,∆tm and Hölder's inequality. Since u i,T ,∆t converges towards u i in L 1 (Q T ), and since u i,E,∆t is uniformly bounded, it suffices to show that u i,E,∆t -u i,T ,∆t L 1 (Q T ) tends to 0. Denote by ∆ Kσ the half-diamond cell which is defined as the interior of the convex hull of {x K , σ} for K ∈ T and σ ∈ E K , then the following geometrical relation holds:

m(∆ Kσ ) = 1 d m σ dist(x K , σ) ≤ h T d m σ .
As a consequence,

u i,E,∆t -u i,T ,∆t L 1 (Q T ) = N T n=1 ∆t n K∈T σ∈E K m ∆ Kσ |u n i,K -u n i,σ | ≤ h T d N T n=1 ∆t n K∈T σ∈E K m σ |u n i,K -u n i,σ |.
As we have u n i,K = u n i,σ , the contributions corresponding to the boundary edges vanish. For σ ∈ E int , u i,σ is an average of u i,K and u i,Kσ , hence

|u n i,K -u n i,σ | ≤ |u n i,K -u n i,Kσ |. Therefore, we obtain that u i,E,∆t -u i,T ,∆t L 1 (Q T ) ≤ h T d N T n=1 ∆t n σ∈Eint 2m σ |D σ u n i | ≤ 2 h T d N T n=1 ∆t n σ∈Eint m σ d σ 1 2 N T n=1 ∆t n σ∈Eint τ σ |D σ u n i | 2 1 2
.

We deduce from Corollary 4.2 that u i,E,∆t -u i,T ,∆t L 1 (Q T ) ≤ Ch T , hence u i,E,∆t and u i,T ,∆t share the same limit in L 1 (Q T ).

4.3.

Convergence towards a weak solution. The last step to conclude the proof of Theorem 2.3 is to identify the limit value U exhibited in Proposition 4.1 as a weak solution to (1.1), (1.3) corresponding to the initial profile U ∈ L ∞ (Ω; A). This is the purpose of our last statement. Proof. We drop again the subscript m for the sake of readability, and let

i ∈ [[1, N ]], ϕ ∈ C ∞ c ([0, T ) × Ω)
, then define ϕ = (ϕ n K ) by ϕ n K = ϕ(x K , t n ) for all n ∈ {0, . . . , N T } and K ∈ T . Multiplying (2.4a) by ∆t n ϕ n-1 K , then summing over K ∈ T and n ∈ {1, . . . , N T } leads to (4.5)

T 1 + T 2 + T 3 = 0,
where we have set

T 1 = N T n=1 K∈T m K (u n i,K -u n-1 i,K )ϕ n-1 K , T 2 = N T n=1 ∆t n σ∈E τ σ a D Kσ u n i D Kσ ϕ n-1 , T 3 = N T n=1 ∆t n σ∈E τ σ N j=1 (a i,j -a ) u n jσ D Kσ u n i -u n iσ D Kσ u n j D Kσ ϕ n-1 .
The term T 1 can be rewritten as

T 1 = N T n=1 ∆t n K∈T m K u n i,K ϕ n-1 K -ϕ n K ∆t n - K∈T m K u 0 i,K ϕ 0 K ,
so that it follows from the convergence of π T ,∆t U towards U and of π T U 0 towards U 0 together with the regularity of ϕ that (4.6)

T 1 -→ m→∞ - Q T u i ∂ t ϕdxdt - Ω u 0 i ϕ(0, •)dx.
To treat the term T 2 , we introduce a strongly consistent reconstruction of the gradient. Following [START_REF] Droniou | The asymmetric gradient discretisation method[END_REF] (see [START_REF] Coudière | Convergence rate of a finite volume scheme for a two dimensional convection-diffusion problem[END_REF] for a practical example), one can reconstruct a second approximate gradient operator ∇

T : R T → L ∞ (Ω) d such that ∆σ ∇ T u • ∇ T vdx = τ σ D Kσ uD Kσ v, ∀u, v ∈ R T , ∀σ ∈ E,
and which is strongly consistent, i.e.,

∇ T ϕ n -→ h T →0 ∇ϕ(•, t n ) uniformly in Ω, ∀n ∈ {1, . . . , N T },
thanks to the smoothness of ϕ. Using this tool, the terms T 2 and T 3 , are easy to treat. The first one can be rewritten as:

T 2 = a Q T ∇ T ,∆t u i • ∇ T ,∆t ϕdxdt, so that (4.7) T 2 -→ m→∞ a Q T ∇u i • ∇ϕdxdt.
On the other hand, the term T 3 rewrites

T 3 = Q T N j=1 (a i,j -a ) (u j,E,∆t ∇ T ,∆t u i -u i,E,∆t ∇ T ,∆t u j ) ∇ T ,∆t ϕ, so that (4.8) T 3 -→ m→∞ Q T N j=1 (a i,j -a ) (u j ∇u i -u i ∇u j ) ∇ϕ.
Combining (4.5), (4.6), (4.7), and (4.8), we obtain that

- Q T u i ∂ t ϕdxdt - Ω u 0 i ϕ(0, •)dx + a Q T ∇u i • ∇ϕdxdt + Q T N j=1 (a i,j -a ) (u j ∇u i -u i ∇u j ) ∇ϕdxdt = 0, ∀ϕ ∈ C ∞ c ([0, T ) × Ω).
Using U ∈ A and the relation (1.10), we recover the weak formulation (1.9).

5. Numerical results. The numerical scheme has been implemented using MATLAB. The nonlinear system corresponding to the scheme is solved thanks to a variation of the Newton method with stopping criterion U n,k+1 -U n,k ∞ ≤ 10 -12 for successful convergence or 20 iteration for failed convergence. The solution of the Newton iteration, U n,k+1/3 , is then "projected" on A by setting U n,k+2/3 = 10 2 some cells, so that the choice a = 0 would allow for spurious solutions as highlighted in Remark 2.1. In this situation, the choice of a strongly affects the quality of the results, especially for the first time steps where some concentrations are still close to 0. The numerical experiment and homogeneity considerations suggest the following suboptimal rule for choosing a :

A = A reg A = A sing U 0 = U 0 s U 0 = U 0 r U 0 = U 0 s U 0 = U 0 r nb.
(5.1)

a = min max i =j a i,j ; max min i =j a i,j , h 2 T τ ,
where h T is the mesh size, τ the current time step and a small parameter to be tuned by the user. Another interesting feature of Figure 3 is the behavior of the curve corresponding to A = A Lap . The classical TPFA scheme for the heat equation corresponding to a = 1 is outperformed in terms of accuracy by the scheme corresponding to higher value of a 9. The introduction of enhanced diffusion by picking high values of a is not covered by formula (5.1).

A 2D test case with reaction.

Our second test is two-dimensional. We choose A sing as the diffusion matrix and a = 0.1. The domain Ω = (0, 22) × (0, 16) is discretized into a cartesian grid made of 110 × 80 cells which is presented along with initial condition in Figure 4. We use a uniform time stepping with τ = 2 -3 . To The reaction term R obviously satisfies Assumptions (i) and (ii) of Remark 1.1. Let us now discuss Assumption (iii). We have for all U, U ∈ A: R(U ). ln(U/U ) = r 1 (U ) ln(u 1 u 3 ) -ln(u 2 2 ) -ln(u 1 u 3 ) -ln(u 2 2 ) .

The above expression is then nonpositive for any U such that ln(u 1 u 3 ) -ln(u 2 2 ) = -ln(1000), or equivalently R(U ) = 0. A particular choice for such a U is the steady state U ∞ , which is constant w.r.t. to space and determined as follows. Denote by α the average advancement of the reaction, then

u ∞ 1 (x) = 9 44 -α, u ∞ 2 (x) = 2 11 + 2α, u ∞ 3 (x) = 27 44 -α.
For α = 0, the fraction of each specie is the ratio of corresponding occupied area in Figure 4, i.e. 1 mΩ Ω U 0 = ( 9 44 , 2 11 , 27 44 ) T . The value of α is determined by imposing that R(U ∞ ) = 0, which amounts to find a root of a polynomial of degree two. Among the two roots, only the choice α = -5 √ 206530+4504 10956 yields a non-negative U ∞ . The time evolution of the relative energy E T (U |U ∞ ) is plotted on Figure 5, showing exponential decay to the steady-state even though the diffusion matrix is singular. Snapshots showing the evolution of the concentration profiles are presented in Figure 6.

To compute the solutions to our numerical scheme, we have to adapt the continuation procedure sketched at the beginning of Section 5 to include source terms. Roughly speaking, we solve the discrete counterpart to

∂ t u i -a ∆ u i -λ div   N j=1
(a i,j -a ) u j ∇u i -u i ∇u j   = µr i (U ), where the source terms are discretized in a fully implicit way. Due to the stiffness of the reaction terms, we have to treat the reaction first then the cross-diffusion effects. More precisely, given U n-1 ∈ A T , we initialise the iterative method for the computation of U n with U n,0 defined as the unique solution to the N uncoupled heat equations corresponding to λ = µ = 0. Then one tries to solve the system corresponding to λ = µ = 1. If the modified Newton's method with truncation and rescaling sketched at the beginning of Section 5 fails to converge, then one sets λ = 0 and µ = 1 2 . Then one use a similar continuation method to the one described at the beginning of Section 5 to increase µ until it reaches the value 1. Then the continuation method is used again to increase the value of λ until λ = 1 is reached.

u 0 1 = 1 u 0 2 = 1 u 0 3 = 1

Conclusion.

We proposed a finite volume scheme based on two-point flux approximation for a degenerate cross-diffusion system. The scheme was designed to preserve the key properties of the continuous system, namely the positivity of the solutions, the constraint on the composition and the decay of the entropy. The scheme requires the introduction of a positive parameter a to avoid unphysical solutions. This parameter plays an important role in the convergence proof, which is carried out under a non-degeneracy assumption. Its importance is also confirmed in the numerical experiments, in particular in the presence of initial profiles with concentrations vanishing in some parts of the computational domain.

Lemma 1 . 1 (

 11 conservation of mass). (1.1) and (1.3) corresponding to an initial data U 0 ∈ L ∞ (Ω; A), thenΩ u i (t, x)dx = Ω u 0 i (x)dx, ∀t ∈ [0, T ], ∀i ∈ [[1, N ]].

Figure 1 :

 1 Figure 1: Illustration of an admissible mesh as in Definition 2.1.

Lemma 4 . 3 .

 43 We have, for any i ∈ [[1, N ]]:

Proposition 4 . 2 .

 42 Let U be as in Proposition 4.1, then U is a weak solution in the sense of Definition 1.4.

L 2 normA

 2 = A reg and U 0 = U 0 s A = A reg and U 0 = U 0 r A = A sing and U 0 = U 0 s A = A sing and U 0 = U 0 r A = A Lap and U 0 = U 0 s idem with ∆t n = 2 -25

Figure 2 :

 2 Figure 2: Error with respect to reference solution.

10 - 5 AsFigure 3 : 1 2e 2 ,+ 3 r 2 (

 531232 Figure 3: Evolution of the ratio U a -U ref 2 U a opt -U ref 2 , where U a is computed with 2 7 cells and U ref is as in Section 5.1.

Figure 4 :Figure 5 :

 45 Figure 4: Initial configuration U 0 for the concentrations

  Profile of u1 at time t = 20Profile of u2 at time t = 20Profile of u1 at time t = 50Profile of u2 at time t = 50

Figure 6 :

 6 Figure 6: Concentration configurations for various times. The concentration of the third specie can be deduced thanks to u 1 + u 2 + u 3 = 1

Table 1 :

 1 Values of a opt for different parameters. a opt is computed with respect to the reference solution of Section 5.1 for the L 2 norm.

	of cells	32 128	T = 0.125 T = 0.25 T = 0.125 T = 0.25	0.86 0.67 0.86 0.67	0.21 0.13 0.17 0.11	0.79 0.49 0.79 0.49	0.0023 0.00082 0.00050 0.00049

Acknowledgements. The authors acknowledge support from the Labex CEMPI (ANR-11-LABX-0007-01). Clément Cancès also acknowledges support from the COMODO project (ANR-19-CE46-0002), and he warmly thanks Virginie Ehrlacher and Laurent Monasse for stimulating discussions that were at the origin of this work.

max(U n,k+1/3 , 10 -10 τ ), and then for all K ∈ T : U n,k+1

). The seed of the algorithm is the solution to N uncoupled heat equations with diffusion coefficients all equal to a .

For the first time step, we also make use of a continuation method based on the intermediate diffusion coefficients

. The parameter λ is originally set to 1. If Newton's method does not converge, we let λ = (λ + λ prev )/2 where λ prev is originally set to 0. If Newton's method converges, we let λ prev = λ and λ = 1.

5.1.

Convergence under grid refinement. Our first test case is devoted to the convergence analysis of the scheme in a one-dimensional setting Ω = (0, 1). Two different initial conditions are considered: U 0 s is smooth and vanished point-wise at the boundary of Ω, whereas U 0 r is discontinuous and vanishes on intervals of Ω:

. We also consider three cross-diffusion coefficients matrices, a first one A Lap corresponding to 3 uncoupled heat equation, a second one A reg called regular with positive off-diagonal coefficients, and a third one A sing called singular with a few null off-diagonal coefficients:

For the convergence tests, we have let a = 0.1 and the meshes are uniform discretisations of [0, 1] from 2 5 cells to 2 15 cells. The approximate solutions are compared to a reference solution which is analytical when A = A Lap and computed on the finest grid (2 15 cells) when A = A reg or A = A sing . The final time is 0.25, and the time discretisation is fixed with a time step of ∆t = 2 -18 . In the case A = A Lap , we also report in Figure 2 results with 128 times finer time step, i.e. ∆t = 2 -25 . One notices that our scheme is second-order accurate in space in the setting presented in this paper (A = A reg ), but only first-order accurate when confronted to non-diffusive discontinuities. We call non-diffusive discontinuities a spatial discontinuity of u 0 1 and u 0 2 (recall that a 1,2 = 0 in A sing ) for which u 0 3 is equal to 0 on both sides of the discontinuity, so that the contributions corresponding to a 1,3 and a 2,3 vanish at t = 0. The origin of this lower order might lie in the difficulty to compute accurately the near-zero concentrations in the neighborhood of such discontinuities. We also notice in Figure 2 the prevalence of the error in time when comparing with respect to an analytical solution, which will motivate the development of higher-order in time methods already discussed in Remark 2.2. 5.2. On the influence of the parameter a . The choice of a is a natural question concerning our scheme. The equation (2.3) gives a lower bound: a > 0. The existence of an upper bound is not as clear. Equation (3.7) shows that for large a , we over-estimate the diffusion. The optimal value of a depends on many variables such as the initial condition, the final time, and the mesh. Optimal choices of a are reported in Table 1. Notice that the optimal value is test cases dependent, since it is affected by the initial condition and by the final time.

One notices on Fig. 3 that the dependency of the quality of the results is strong for the initial data U 0 r . This is due to the presence of vanishing concentrations in