Maxime Herda
email: maxime.herda@inria.fr

Antoine Zurek
email: antoine.zurek@utc.fr

STUDY OF A STRUCTURE PRESERVING FINITE VOLUME SCHEME FOR A NONLOCAL CROSS-DIFFUSION SYSTEM

Keywords: Nonlocal cross-diffusion, finite volume schemes, entropy method, convergence. Mathematics Subject Classification: 65M08, 65M12, 35K51, 35Q92, 92D25

In this paper we analyse a finite volume scheme for a nonlocal version of the Shigesada-Kawazaki-Teramoto (SKT) cross-diffusion system. We prove the existence of solutions to the scheme, derive qualitative properties of the solutions and prove its convergence. The proofs rely on a discrete entropy-dissipation inequality, discrete compactness arguments, and on the novel adaptation of the so-called duality method at the discrete level. Finally, thanks to numerical experiments, we investigate the influence of the nonlocality in the system: on convergence properties of the scheme, as an approximation of the local system and on the development of diffusive instabilities.

Introduction

We are interested in the numerical discretization of the following nonlocal cross-diffusion system

∂ t u 1 -∆((d 1 + d 11 σ 1 * u 1 + d 12 ρ 1 * u 2 )u 1 ) = R 1 (u 1 , u 2 ), (1)
∂ t u 2 -∆((d 2 + d 21 ρ 2 * u 1 + d 22 σ 2 * u 2 )u 2 ) = R 2 (u 1 , u 2 ), (2)
on a periodic domain Ω = T d (d ≤ 3). For a given final time T we denote the space time domain by Q T = Ω × (0, T ). The parameters d 1 , d 2 , d 11 , d 12 , d 21 and d 22 are some positive constants and ρ 1 , ρ 2 , σ 1 and σ 2 are non-negative convolution kernels. System (1)-( 2) is supplemented with initial conditions

u 1 (•, 0) = u 0 1 (•), u 2 (•, 0) = u 0 2 (•). (3)
In the case where the convolution kernels are given by the Dirac measure ρ 1 = ρ 2 = σ 1 = σ 2 = δ 0 , the system coincides with the celebrated Shigesada, Kawasaki, and Teramoto (SKT) population model [START_REF] Shigesada | Spatial segregation of interacting species[END_REF] which can describe segregation phenomena between competing species. It writes

∂ t u 1 -∆((d 1 + d 11 u 1 + d 12 u 2 )u 1 ) = R 1 (u 1 , u 2 ), (4)
∂ t u 2 -∆((d 2 + d 21 u 1 + d 22 u 2 )u 2 ) = R 2 (u 1 , u 2 ). ( 5
)
Nonlocal cross-diffusion systems appear naturally as a mean field type of limit of interacting many-particle systems. For instance, the model (1)- [START_REF] Anaya | A convergent finite volume method for a model of indirectly transmitted diseases with nonlocal cross-diffusion[END_REF] was introduced in [START_REF] Fontbona | Non local Lotka-Volterra system with cross-diffusion in an heterogeneous medium[END_REF] as the large population limit of a stochastic individual model. If these particle systems allow a precise description of the interactions between individuals, their numerical approximations are very time-consuming. Then, it is reasonable to investigate simpler macroscopic models. In this context we see nonlocal cross-diffusion models as intermediate models between individual based models and local crossdiffusion models. This interpretation has been mathematically justified in the literature, see [START_REF] Chen | Rigorous derivation of population cross-diffusion systems from moderately interacting particle systems[END_REF][START_REF] Dietert | Persisting entropy structure for nonlocal cross-diffusion systems[END_REF][START_REF] Jüngel | Nonlocal cross-diffusion systems for multi-species populations and networks[END_REF][START_REF] Moussa | From nonlocal to classical Shigesada-Kawasaki-Teramoto systems: triangular case with bounded coefficients[END_REF], where the derivation of some local cross-diffusion models from nonlocal models (some of them derived from microscopic models) are shown.

Besides, nonlocal cross-diffusion models can be more than a mathematical intermediate between two scales. Indeed, in population dynamics, they can model nonlocal sensing, as diffusion of a species is impacted by the population located (respectively to their position) on the support of the convolution kernels, see [START_REF] Giunta | Local and global existence for non-local multi-species advection-diffusion models[END_REF][START_REF] Potts | Spatial memory and taxis-driven pattern formation in model ecosystems[END_REF]. In the model ( 1)-( 2) assume for instance that ρ 1 is supported away from 0. Then the resulting effect of the nonlocal cross-diffusion term is to enhance the diffusion of species 1 when species 2 is away, modeling for instance a hunting behavior in a predator-prey model. This could hardly be reproduced by local cross-diffusion terms.

The ability of the nonlocal cross-diffusion terms to model the dynamics of some natural phenomena explain the use of such models in other contexts. They are for instance applied to describe cell sorting [START_REF] Murakawa | Continuous models for cell-cell adhesion[END_REF][START_REF] Painter | A nonlocal model for contact attraction and repulsion in heterogeneous cell populations[END_REF], tumour growth [START_REF] Domschke | Mathematical modelling of cancer invasion: implications of cell adhesion variability for tumour infiltrative growth patterns[END_REF], opinion formation [START_REF] Düring | Boltzmann and Fokker-Planck equations modelling opinion formation in the presence of strong leaders[END_REF] or interactions between spiking neurons [START_REF] Baladron | Mean-field description and propagation of chaos in networks of Hodgkin-Huxley and FitzHugh-Nagumo neurons[END_REF] (just to name a few). In particular, the development of reliable numerical methods to approximate the solutions of nonlocal cross-diffusion systems can enhance our understanding of the "physical" mechanisms described by them. As a by-product this could also help the development of efficient models describing complex phenomena.

Motivated by these reasons, this manuscript deals with the design and analysis of a robust numerical scheme for (1)- [START_REF] Anaya | Numerical analysis for a three interacting species model with nonlocal and cross diffusion[END_REF]. Our approach is inspired by the analysis performed at the continuous level in [START_REF] Dietert | Persisting entropy structure for nonlocal cross-diffusion systems[END_REF][START_REF] Moussa | From nonlocal to classical Shigesada-Kawasaki-Teramoto systems: triangular case with bounded coefficients[END_REF]. In particular, in [START_REF] Dietert | Persisting entropy structure for nonlocal cross-diffusion systems[END_REF] the authors show that there is a persisting entropy structure in the nonlocal case which yields a crucial a priori estimate for the analysis of the model. This extends for instance the approach developed 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] in the local case. Indeed, it was shown that for the system (1)-(3) without reaction terms R 1 = R 2 = 0 and under the following symmetry hypotheses on the convolution kernels ( 6)

     ρ 1 (x) = ρ 2 (-x) = ρ(x), σ 1 (x) = σ 1 (-x), σ 2 (x) = σ 2 (-x),
that the following entropy functional

H(u 1 , u 2 ) = Ω 1 d 12 [u 1 (log(u 1 ) -1) + 1] dx + Ω 1 d 21 [u 2 (log(u 2 ) -1) + 1] dx ,
is dissipated along solutions of ( 1)-( 2). More precisely one has Observe that if the convolution kernels are given by the Dirac measure ρ 1 = ρ 2 = σ 1 = σ 2 = δ 0 , then

(7) d dt H(u 1 , u 2 ) + 2d 11 d 12 Ω Ω σ 1 (y) u 1 (x -y) ∇ u 1 (x) + u 1 (x)∇ u 1 (x -y) 2 dxdy + 2d 22 d 21 Ω Ω σ 2 (y) u 2 (x -y) ∇ u 2 (x) + u 2 (x)∇ u 2 (x -y)
d dt H(u 1 , u 2 ) + 2 d 11 d 12 Ω |∇u 1 | 2 dx + 2 d 22 d 21 Ω |∇u 2 | 2 dx + 4 d 1 d 12 Ω |∇ √ u 1 | 2 dx + 4 d 2 d 21 Ω |∇ √ u 2 | 2 + 4 Ω |∇ √ u 1 u 2 | 2 dx = 0 ,
which was already known for the local SKT system (see [START_REF] Chen | Analysis of a multidimensional parabolic population model with strong cross-diffusion[END_REF][START_REF] Chen | Analysis of a parabolic cross-diffusion population model without self-diffusion[END_REF][START_REF] Galiano | Semi-discretization in time and numerical convergence of solutions of a nonlinear cross-diffusion population model[END_REF]). The fact that (1)-( 2) admits a Lyapunov functional is crucial for the study of the system. Indeed, in [START_REF] Dietert | Persisting entropy structure for nonlocal cross-diffusion systems[END_REF], the authors used [START_REF] Braukhoff | An entropy structure preserving space-time formulation for crossdiffusion systems: Analysis and galerkin discretization[END_REF] together with the so-called duality method, see [START_REF] Desvillettes | On the entropic structure of reaction-cross diffusion systems[END_REF][START_REF] Lepoutre | Entropic structure and duality for multiple species cross-diffusion systems[END_REF][START_REF] Moussa | From nonlocal to classical Shigesada-Kawasaki-Teramoto systems: triangular case with bounded coefficients[END_REF], in order to prove (assuming [START_REF] Bessemoulin-Chatard | Hypocoercivity and diffusion limit of a finite volume scheme for linear kinetic equations[END_REF] and without reaction terms) the existence of distributional solutions to (1)-(3). Definition 1. Given T > 0, let ρ 1 , ρ 2 , σ 1 and σ 2 be some functions in L ∞ (Ω) and u 0 1 and u 0 2 be some initial functions in L 1 (Ω). Then, we say that the measurable functions u 1 , u 2 :

Q T → R + are distributional solutions to (1)-(3) if for every φ ∈ C ∞ 0 (Ω × [0, T )) it holds Q T u 1 ∂ t φ + (d 1 u 1 + d 11 σ 1 * u 1 u 1 + d 12 ρ 1 * u 2 u 1 ) ∆φ dxdt = - T d
u 0 1 (x)φ(x, 0)dx, [START_REF] Burger | Segregation effects and gap formation in cross-diffusion models[END_REF] and

Q T u 2 ∂ t φ + (d 2 u 2 + d 21 ρ 2 * u 1 u 2 + d 22 σ 2 * u 2 u 2 ) ∆φ dxdt = - T d u 0
2 (x)φ(x, 0)dx. [START_REF] Burger | Nonlinear cross-diffusion with size exclusion[END_REF] In this paper we propose and analyze a finite volume scheme for (1)- [START_REF] Anaya | Numerical analysis for a three interacting species model with nonlocal and cross diffusion[END_REF]. A particular focus is put on (i) the preservation of the entropy dissipation property at the discrete level; (ii) the non-negativity of the solution; (iii) the possibility to use the scheme in both the nonlocal and local regimes.

In order to achieve these goals we will design a fully implicit two point flux approximation (TPFA) finite volume scheme. As in the study of some numerical schemes for local cross-diffusion systems, see for instance the following (non-exhaustive) list of contributions [START_REF] Andreianov | Analysis of a finite volume method for a cross-diffusion model in population dynamics[END_REF][START_REF] Braukhoff | An entropy structure preserving space-time formulation for crossdiffusion systems: Analysis and galerkin discretization[END_REF][START_REF] Cancès | Finite-volume scheme for a degenerate cross-diffusion model motivated from ion transport[END_REF][START_REF] Cancès | A convergent entropy diminishing finite volume scheme for a cross-diffusion system[END_REF][START_REF] Jüngel | A convergent structure-preserving finite-volume scheme for the Shigesada-Kawasaki-Teramoto population system[END_REF][START_REF] Sun | An entropy stable high-order discontinuous Galerkin method for cross-diffusion gradient flow systems[END_REF], the preservation of the entropy dissipation property at the discrete level is crucial. This ensures well-posedness and global stability in time [START_REF] Chainais-Hillairet | Large-time behaviour of a family of finite volume schemes for boundarydriven convection-diffusion equations[END_REF][START_REF] Filbet | A finite volume scheme for boundary-driven convection-diffusion equations with relative entropy structure[END_REF] as well as with respect to the choice of convolution kernels (see Theorem 1). Some of these methods are reminiscent of the second author's work in [START_REF] Jüngel | A convergent structure-preserving finite-volume scheme for the Shigesada-Kawasaki-Teramoto population system[END_REF] concerning the study of a finite volume scheme for the local SKT system. Besides, we are able to obtain additional estimates on the solution (see Theorem 2) by adapting the duality method (see [START_REF] Desvillettes | On the entropic structure of reaction-cross diffusion systems[END_REF][START_REF] Lepoutre | Entropic structure and duality for multiple species cross-diffusion systems[END_REF][START_REF] Moussa | From nonlocal to classical Shigesada-Kawasaki-Teramoto systems: triangular case with bounded coefficients[END_REF]) at the discrete level. This technique relies on the study of a discretized Kolmogorov equation, see Section 4. The convergence of solutions of the numerical scheme towards distributional solutions in the sense of Definition 1 is shown in Theorem 3. Let us mention that only L ∞ regularity of the convolution kernels is required to obtain convergence of the scheme. With smoother kernels one additionally obtains local in time L ∞ bounds on the discrete solution that are uniform in the mesh size (see Theorem 2).

Let us notice that there already exists some works dealing with the design and the analysis of finite volume numerical schemes for nonlocal cross-diffusion systems. Indeed, in [START_REF] Anaya | A convergent finite volume method for a model of indirectly transmitted diseases with nonlocal cross-diffusion[END_REF][START_REF] Anaya | Numerical analysis for a three interacting species model with nonlocal and cross diffusion[END_REF] the convergence of some semi-implicit TPFA finite volume schemes are proved. The convergence proofs are based (as in this work) on the adaptation at the discrete level of a Kruzhkov's compactness result [START_REF] Kruzhkov | Results concerning the nature of the continuity of solutions of parabolic equations and some of their applications[END_REF] obtained in [START_REF] Andreianov | Analysis of a finite volume method for a cross-diffusion model in population dynamics[END_REF]. We mention [START_REF] Burger | Segregation effects and gap formation in cross-diffusion models[END_REF] where numerical experiments are shown to illustrate the formation of gaps for a class of nonlocal cross-diffusion systems. In this paper the authors applied an explicit in time finite volume scheme first introduced in [START_REF] Carrillo | A finite-volume method for nonlinear nonlocal equations with a gradient flow structure[END_REF] and then extended for the multi-species case in [START_REF] Carrillo | Zoology of a nonlocal cross-diffusion model for two species[END_REF]. This scheme is a positivity and entropy preserving method as shown in [START_REF] Carrillo | Zoology of a nonlocal cross-diffusion model for two species[END_REF]. Finally, we also refer to [START_REF] Carrillo | Convergence of a finite volume scheme for a system of interacting species with cross-diffusion[END_REF]. In this contribution the convergence of a semi-discrete finite volume scheme is proved. This scheme is also positivity preserving which allows the authors to establish a discrete energy estimate.

In this paper the system (1)-( 3) is set on a periodic domain. A natural question concerns the extension to different boundary conditions. In the case of Neumann boundary conditions in a bounded domain Ω ⊂ R d the nonlocal entropy preserving version of the cross-diffusion model differs from (4)-( 5) because of boundary effects (see [START_REF] Dietert | Persisting entropy structure for nonlocal cross-diffusion systems[END_REF]Section 1.3]). Since the model is different, it is not straightforward to adapt the definition of the scheme introduced in Section 2.2. We postpone this discussion in Appendix B where a finite volume scheme for the nonlocal SKT system in a bounded domain is designed and a priori entropy estimates are obtained.

In order to illustrate and complement the theoretical results, we present several numerical experiments in the last section of this paper. We compute the experimental order of convergence of the numerical method when the mesh size goes to 0 for various initial data and convolution kernels. Then, for a fixed mesh, we investigate the rate of convergence for different metrics of the discrete solutions of the nonlocal system towards solution of the local system when convolution kernels tends to Dirac measures. Finally, we perform simulations of the model with nonzero reaction terms with parameters chosen to describe a prey-predator system with either linear diffusion or non-local cross diffusion modelling hunting behavior. For these models we illustrate the persistence and the modification of Turing patterns in the presence of cross-diffusion.

The paper is organized as follows, in Section 2 we introduce the scheme and state our main results. Section 3 is concerned with the proof of existence of positive solutions to the scheme. We introduce the discrete Kolmogorov equation in Section 4. Then we deduce from the study of this problem some qualitative properties satisfied by the solutions to our scheme in Section 4.3. Sections 5 deals with the convergence of the scheme. Finally, in Section 6, we discuss the implementation and show some numerical experiments in one and two space dimensions.

Numerical scheme and main results

The results of this paper apply to a periodic domain Ω = d i=1 R/L i Z, where L 1 , . . . , L d > 0. However, for the sake of readability we will assume from now on that d = 1 and L 1 = 1, namely Ω = T = R/Z. The generalization in higher dimensions on Cartesian grid is immediate by defining the scheme as the tensorization of the one dimensional scheme.

2.1. Notations and definitions. Let us define N ≥ 1 and ∆x = 1/N . A uniform mesh T of T consists in a finite sequence of cells denoted by

K i = (x i-1 2 , x i+ 1 2 ) , i ∈ I = Z/N Z. centered at x i = i∆x and with extremities x i± 1 2 = (i ± 1
2 )∆x. For T > 0 given, we define an integer N T and a time step ∆t = T /N T and we introduce the sequence (t k ) 0≤k≤N T with t k = k∆t. We denote by D a space-time discretization of Q T = T×(0, T ) composed of a space discretization T of T and the values (∆t, N T ).

Let us now introduce some discrete norms on the space of piecewise constant functions in space

H T = w : T → R : w(x) = i∈I w i 1 K i (x) .
For p ∈ [1, ∞), we define the discrete W 1,p seminorm and discrete W 1,p norm on H T by

|w| 2 1,p,T = i∈I ∆x w i+1 -w i ∆x p 1 p , w 1,p,T = |w| 1,p,T + w L p (T) ,
where, for p ∈ [1, ∞), the norm • L p (T) denotes the usual L p (T) norm. In the case p = ∞, we denote

• L ∞ (T) the L ∞ (T) norm given by w L ∞ (T) = max i∈I |w i |, ∀w ∈ H T .
Let us also recall the definition of the space BV (T), see [START_REF] Ambrosio | Functions of bounded variation and free discontinuity problems[END_REF] for more details. A function w ∈ L 1 (T) belongs to the space BV (T) if its total variation T V (w) given by

T V (w) = sup T w(x) ∂ x φ(x) dx, φ ∈ C 1 c (T), |φ(x)| ≤ 1 ∀x ∈ T ,
is finite. We endow the space BV (T) with the norm

w BV (T) = w L 1 (T) + T V (w), ∀w ∈ BV (T).
In particular, we notice that for each function w ∈ H T ∩ BV (T) we have w BV (T) = w 1,1,T . Finally we introduce the space H D of piecewise constant in time functions with values in H T ,

H D = w : T × [0, T ] → R : w(x, t) = N T k=1 w k (x)1 (t k-1 ,t k ] (t) .
This space can be equipped, for (p, q) ∈ [1, ∞) 2 , with the following discrete L q (0, T ; W 1,p (T)) norm

N T k=1 ∆t w q 1,p,T 1 q ∀w ∈ H D ,
or with the L q (0, T ; L p (T)) norm

N T k=1 ∆t w q L p (T) 1 q
, ∀w ∈ H D .

In particular, in the case p = q = 2, the L 2 (Q T ) norm can also be defined by duality as

w 2 L 2 (Q T ) = sup T 0 T wf dxdt : f ∈ H D , N T k=1 ∆t i∈I ∆x|f k i | 2 = 1 , ∀w ∈ H D . (10)
This dual formulation of • L 2 (Q T ) will be needed later on. 2.2. Numerical scheme. We discretize the initial conditions (3) as

u 0 j,i = 1 ∆x K i u 0 j (x) dx, ∀i ∈ I, j = 1, 2. ( 11
)
Now for given (u k-1 1 , u k-1 2 ) ∈ R 2N , the implicit in time numerical scheme writes as

u k j,i -u k-1 j,i ∆t -∆ T (µ k j u k j ) i = R j (u k 1,i , u k 2,j ), ∀i ∈ I, j = 1, 2, (12)
where ∆ T denotes the discrete Laplacian, namely

∆ T (µ k j u k j ) i = µ k j,i+1 u k j,i+1 -2µ k j,i u k j,i + µ k j,i-1 u k j,i-1 ∆x 2 , ∀i ∈ I, j = 1, 2, (13)
and

µ k 1,i = d 1 + d 11 n∈I ∆xσ 1,i-n u k 1,n + d 12 n∈I ∆xρ 1,i-n u k 2,n , ∀i ∈ I, (14)
µ k 2,i = d 2 + d 21 n∈I ∆xρ 2,i-n u k 1,n + d 22 n∈I ∆xσ 2,i-n u k 2,n , ∀i ∈ I, (15)
with

ρ j,i-n = 1 ∆x K i-n ρ j (y) dy, σ j,i-n = 1 ∆x K i-n σ j (y) dy, ∀i, n ∈ I, j = 1, 2. ( 16
)
Let us notice, by construction, that if we consider ρ 1 = ρ 2 = σ 1 = σ 2 = δ 0 then (11)-( 15) yields a finite volume scheme for the local SKT model ( 4)- [START_REF] Baladron | Mean-field description and propagation of chaos in networks of Hodgkin-Huxley and FitzHugh-Nagumo neurons[END_REF]. We also remark that we could equivalently rewrite [START_REF] Carrillo | A finite-volume method for nonlinear nonlocal equations with a gradient flow structure[END_REF] as

∆x u k j,i -u k-1 j,i ∆t + F k j,i+ 1 2 -F k j,i-1 2 = R j (u k 1,i , u k 2,i ), ∀i ∈ I, j = 1, 2, (17)
where for all i ∈ I the numerical fluxes F k j,i+ 1 2 are defined by ( 18)

F k j,i+ 1 2 = u k j,i µ k j,i -u k j,i+1 µ k j,i+1 ∆x = µ k j,i+ 1 2 u k j,i -u k j,i+1 ∆x + u k j,i+ 1 2 µ k j,i -µ k j,i+1 ∆x , j = 1, 2,
with the centered approximation at interfaces [START_REF] Bessemoulin-Chatard | Hypocoercivity and diffusion limit of a finite volume scheme for linear kinetic equations[END_REF]. In particular ρ i = ρ 1,i = ρ 2,-i for all i ∈ I. (H4) The initial data u 0 1 and u 0 2 are non-negative L 1 (T) functions with finite entropy, namely

µ k j,i+ 1 2 = µ k j,i + µ k j,i+1 2 , u k j,i+ 1 2 = u k j,i + u k j,i+1
h 1 (u 0 1 ), h 2 (u 0 2 ) ∈ L 1 (T). (H5) The reaction terms satisfy R 1 = R 2 = 0.
As already mentioned, hypothesis (H1) is only made for the convenience of the reader and one can adapt the design of the scheme and the results to a d-dimensional periodic domain Ω = d i=1 R/L i Z. For different boundary conditions, we refer to Appendix B. Observe that by assuming (H2) we require cross-diffusion on both species. While this is crucial in our proofs, the scheme performs well in practice even with d 12 = 0 or d 21 = 0 (see Section 6). The assumption (H3) on the symmetry of the functions ρ 1 , ρ 2 , σ 1 and σ 2 are needed, as at the continuous level, to show the discrete entropy inequality satisfied by the solutions of the scheme (11)- [START_REF] Chainais-Hillairet | Large-time behaviour of a family of finite volume schemes for boundarydriven convection-diffusion equations[END_REF]. However, in terms of practical use, the scheme performs well even when dropping this hypothesis (see Section 6.4). Following for instance [START_REF] Jüngel | A convergent structure-preserving finite-volume scheme for the Shigesada-Kawasaki-Teramoto population system[END_REF], the assumption (H5) can be relaxed and one can extend the proofs of Theorem 1 and Theorem 3 in the case of the Lotka-Volterra source terms:

R j (u 1 , u 2 ) = u j a j0 - 2 k=1 a jk u k , j = 1, 2,
with a j0 and a jk some nonnegative constants for j, k = 1, 2.

Our first main result deals with the existence of solutions to scheme (11)-( 15) at each time step. But first let us recall the definition of the discrete entropy functional

H(u k 1 , u k 2 ) = i∈I ∆x h 1 (u k 1,i ) + i∈I ∆x h 2 (u k 2,i ),
where the functions h 1 and h 2 are defined by

h 1 (x) = 1 d 12 (x (log(x) -1) + 1) , h 2 (x) = 1 d 21 (x (log(x) -1) + 1)
, ∀x ∈ (0, +∞), [START_REF] Chen | Analysis of a parabolic cross-diffusion population model without self-diffusion[END_REF] with the obvious continuous extension at x = 0. The corresponding entropy dissipation functional is defined by

(20) D(u k 1 , u k 2 ) = 2 d 11 d 12 j∈I σ 1,j i∈I u k 1,i+1 u k 1,i+1-j -u k 1,i u k 1,i-j 2 + 2 d 22 d 21 j∈I σ 2,j i∈I u k 2,i+1 u k 2,i+1-j -u k 2,i u k 2,i-j 2 + 4 d 1 d 12 u k 1 2 1,2,T + 4 d 2 d 21 u k 2 2 1,2,T + 4 j∈I ρ j i∈I u k 1,i+1 u k 2,i+1-j -u k 1,i u k 2,i-j 2 .
Theorem 1 (Existence of solutions). Let the assumptions (H1)-(H5) hold. Then, for every 1 ≤ k ≤ N T there exists (at least) one nonnegative solution (u k 1 , u k 2 ) to scheme (11)- [START_REF] Chainais-Hillairet | Large-time behaviour of a family of finite volume schemes for boundarydriven convection-diffusion equations[END_REF]. Moreover, this solution satisfies the following properties:

(i) Mass conservation:

i∈I ∆x u k j,i = T u 0 j (x) dx, ∀k ≥ 0, j = 1, 2. ( 21
)
(ii) Entropy production estimate: for all k ≥ 1 it holds

H(u k 1 , u k 2 ) + ∆tD(u k 1 , u k 2 ) ≤ H(u k-1 1 , u k-1 2 ). ( 22
)
Finally for all j ∈ {1, 2}, if d j is positive then u k j is positive for all 0 < k ≤ N T . The proof of existence of Theorem 1 is based on a consequence (see [START_REF] Evans | Partial differential equations[END_REF]Section 9.1]) of the Brouwer fixed point Theorem. It can be applied thanks to the a priori entropy-dissipation estimate [START_REF] Domschke | Mathematical modelling of cancer invasion: implications of cell adhesion variability for tumour infiltrative growth patterns[END_REF] and regularization inspired by [START_REF] Burger | Nonlinear cross-diffusion with size exclusion[END_REF][START_REF] Jüngel | Entropy methods for diffusive partial differential equations[END_REF]. It also follows the line of the existence proof of [START_REF] Chainais-Hillairet | Long-time behaviour of hybrid finite volume schemes for advection-diffusion equations: linear and nonlinear approaches[END_REF].

The second main result is concerned by some properties satisfied by the solutions of the scheme (11)- [START_REF] Chainais-Hillairet | Large-time behaviour of a family of finite volume schemes for boundarydriven convection-diffusion equations[END_REF]. These estimates are discrete counterparts of [START_REF] Dietert | Persisting entropy structure for nonlocal cross-diffusion systems[END_REF]Theorem 9].

Theorem 2 (Qualitative properties of the solutions). Let the assumptions of Theorem 1 hold. Moreover, assume that u 0 1 , u 0 2 ∈ L 2 (T) and let γ and Γ be some nonnegative constants such that

γ ≤ u 0 1 (x), u 0 2 (x) ≤ Γ a.e. x ∈ T. Finally let us introduce m 0 j = u 0 j L 1 (T) for j = 1, 2.
Then the following properties hold. (i) Maximum principle: If ρ, σ 1 and σ 2 are twice continuously differentiable functions and that the time step satisfies the condition

∆t < 1/ min{d 11 m 0 1 ∆σ 1 L ∞ (T) , d 22 m 0 2 ∆σ 2 L ∞ (T) } + min{d 12 m 0 2 , d 21 m 0 1 } ∆ρ L ∞ (T)
, then for all i ∈ I, k ≥ 1 and j = 1, 2 we have

e k ≤ u k j,i ≤ E k
, where

e k = γ 1 + ∆t min{d 11 m 0 1 ∆σ 1 L ∞ (T) , d 22 m 0 2 ∆σ 2 L ∞ (T) } + min{d 12 m 0 2 , d 21 m 0 1 } ∆ρ L ∞ (T) -k , E k = Γ 1 -∆t min{d 11 m 0 1 ∆σ 1 L ∞ (T) , d 22 m 0 2 ∆σ 2 L ∞ (T) } + min{d 12 m 0 2 , d 21 m 0 1 } ∆ρ L ∞ (T) -k .
(ii) Duality estimate: If d 1 and d 2 are positives, then there exists a constant C > 0 which is independent of the mesh size such that

N T k=1 ∆t i∈I ∆x µ k 1,i u k 1,i + µ k 2,i u k 2,i u k 1,i + u k 2,i ≤ C(1 + T A) u 0 1 2 L 2 (T) + u 0 2 2 L 2 (T) ,
where

A = d 1 + d 2 + d 11 m 0 1 σ 1 L 1 (T) + d 22 m 0 2 σ 2 L 1 (T) + ρ L 1 (T) d 12 m 0 2 + d 21 m 0 1 .
The proof of Theorem 2 relies on a discrete duality method. In Section 4, we define and study the properties satisfied by the finite volume solutions to the Kolmogorov equation. Then, in Section 4.3, we apply these results on the solutions to the scheme ( 11)-( 15) in order to establish the Theorem. Let us emphasis that the duality estimate holds without any assumptions on the time step or the regularity of the convolution kernels. This implies in particular that this discrete estimate also holds for the solutions to the local SKT system.

Finally, we show the convergence of the solutions to the scheme ( 11)-( 15) towards a distributional solution to (1)-( 3) in the sense of Definition 1. However, in order to state precisely our convergence result, we need some notations.

We introduce a family (D m ) m∈N of space-time discretizations of Q T indexed by the size η m = max{∆x m , ∆t m } of the mesh, satisfying η m → 0 as m → ∞. We denote by T m the corresponding mesh of T and by ∆t m the corresponding time step. Finally, for every m ∈ N we set (u 1,m , u 2,m ) ∈ H Dm the picewise constant in space and time reconstruction of the solutions to the scheme ( 11)-( 15) corresponding to the mesh D m .

Theorem 3 (Convergence of the scheme). Let the assumptions of Theorem 1 hold, assume that the coefficients d 1 and d 2 are positives and let (D m ) m∈N be a family of space-time discretizations of Q T with η m → 0 as m → ∞. Then, if we denote by (u 1,m , u 2,m ) a family of finite volume solutions to (11)-( 15) obtained in Theorem 1, there exists (u 1 , u 2 ) ∈ (L p (Q T )) 2 for p ∈ [1, 3) a distributional solutions to (1)-(3) in the sense of Definition 1 such that, up to a subsequence, for j = 1, 2 it holds

u j,m → u j strongly in L p (Q T ) for 1 ≤ p < 3 as m → ∞.
The proof of Theorem 3 is based on uniform estimates w.r.t. ∆x and ∆t, established in Section 5.1. These estimates allow us to apply in Section 5.2 a compactness result obtained in [START_REF] Andreianov | Analysis of a finite volume method for a cross-diffusion model in population dynamics[END_REF] which yields, up to a subsequence, the strong convergence in L p (Q T ) of the sequence (u 1,m , u 2,m ) towards the functions u 1 and u 2 stated in Theorem 3. Then, we identify in Section 5.3 the functions u 1 and u 2 as distributional solutions in the sense of Definition 1 of the nonlocal cross-diffusion system (1)-(3). Remark 4. Let us notice that if d > 1, the convergence of the scheme can also be established. However in this case we obtain, up to a subsequence, for j = 1, 2,

u j,m → u j strongly in L p (Q T ) for 1 ≤ p < d d -1
, as m → ∞, see Remark 13 for more details. If d 1 = d 2 = 0, then it is still possible to conclude if the convolution kernels are smooth enough (C 2 for instance). Indeed in this case from weak compactness on (u j,m ) m , strong compactness can be obtained on its convolution with the smooth kernel. Finally, if the convolution kernels are C 2 , then one can prove a stability estimate (in L 2 -norm) for the solutions to (1)-( 3) which provides uniqueness and continuous dependence on the initial data at the continuous and discrete level. As a by-product we deduce that in this case the whole sequence (u 1,m , u 2,m ) converges as m → ∞.

Existence of solution and entropy dissipation estimate

The problem of existence of solution reduces to the resolution of a nonlinear system of equations. The natural unknowns for which a fixed point theorem will be easily applied are linked to the entropy. In our case, given (u 1 , u 2 ) ∈ ((0, +∞) N ) 2 we define the new unknown X = Φ(u 1 , u 2 ) where Φ :

((0, +∞) N ) 2 → R 2N is the smooth diffeomorphism defined by Φ(u 1 , u 2 ) = (d -1 12 log(u 1,1 ), . . . , d -1 12 log(u 1,N ), d -1 21 log(u 2,1 ), . . . , d -1 21 log(u 2,N )) ∈ R 2N .
From there finding a positive solution to the scheme ( 11)-( 15) amounts to finding a zero

X k = Φ(u k 1 , u k 2 )
of the continuous map P k : R 2N → R 2N defined for any (u 1 , u 2 ) ∈ ((0, +∞) N ) 2 by its components

P k i+N (j-1) (Φ(u 1 , u 2 )) = ∆x(u j,i -u k-1 j,i ) -∆t∆x (∆ T (µ j u j )) i , ∀i ∈ {1, . . . , N }, j ∈ {1, 2}, where (u k-1 1 , u k-1
2 ) are given and µ 1 , µ 2 are related to u 1 , u 2 through the relation ( 14) and ( 15) dropping the exponent k.

5. Let (u k-1 1 , u k-1
2 ) be componentwise non-negative. Then for any

X ∈ R 2N , ( 23
) P k (X), 1 j = i∈I (u j,i -u k-1 j,i )∆x , ∀j ∈ {1, 2}, and (24)
P k (X), X ≥ H(u 1 , u 2 ) -H(u k-1 1 , u k-1 2 ) + ∆tD(u 1 , u 2 ),
where

1 1 = (1, . . . , 1, 0, . . . , 0), 1 2 = (0, . . . , 0, 1, . . . , 1), (u 1 , u 2 ) = Φ -1 (X) and D(u 1 , u 2 )
denotes the entropy dissipation functional given by [START_REF] Desvillettes | On the entropic structure of reaction-cross diffusion systems[END_REF].

Proof. In order to prove [START_REF] Düring | Boltzmann and Fokker-Planck equations modelling opinion formation in the presence of strong leaders[END_REF], it suffices to sum the components of P k (X) and observe that i∈I

(∆ T (µ j u j )) i = 0,
since it is a telescopic sum. Concerning the inequality, first observe that

H(u k-1 1 , u k-1 2 ) is well- defined since (u k-1 1 , u k-1
2 ) is non-negative. Then, using the definition of Φ(u 1 , u 2 ) and P k one obtains

P k (X), X = ∆t(I 1 + I 2 + J 1 + J 2 ) with I 1 = 1 d 12 i∈I ∆x(u 1,i -u k-1 1,i ) log(u 1,i ), I 2 = 1 d 21 i∈I ∆x(u 2,i -u k-1 2,i ) log(u 2,i ), J 1 = 1 d 12 ∆x i∈I (-µ 1,i+1 u 1,i+1 + 2µ 1,i u 1,i -µ 1,i-1 u 1,i-1 ) log(u 1,i ) J 2 = 1 d 21 ∆x i∈I (-µ 2,i+1 u 2,i+1 + 2µ 2,i u 2,i -µ 2,i-1 u 2,i-1 ) log(u 2,i )
Using the convexity of x → (x log(x) -x + 1) to bound both I 1 and I 2 from below, one obtains

I 1 + I 2 ≥ H(u 1 , u 2 ) -H(u k-1 1 , u k-1 2 ).
Then for J 1 , a discrete integration by parts (or summation by parts) yields

J 1 = 1 d 12 ∆x i∈I (u 1,i+1 µ 1,i+1 -u 1,i µ 1,i ) (log(u 1,i+1 ) -log(u 1,i )) ,
and a similar formula holds for J 2 . Using the definitions of µ 1 and µ 2 (see ( 14) and ( 15) without the exponents), one has

J 1 = J diff 1 + J σ 1 1 + J ρ 1 1 and J 2 = J diff 2 + J σ 2 2 + J ρ 2 2 with J diff 1 = d 1 d 12 ∆x i∈I (u 1,i+1 -u 1,i ) (log(u 1,i+1 ) -log(u 1,i )) ≥ 4d 1 d 12 ∆x i∈I √ u 1,i+1 - √ u 1,i 2 ,
and a similar estimate for J diff 2 . For the second term one has

J σ 1 1 = d 11 d 12 j∈I σ 1,j i∈I (u 1,i+1 u 1,i+1-j -u 1,i u 1,i-j ) (log(u 1,i+1 ) -log(u 1,i )) = 1 2 J σ 1 1 + d 11 2d 12 j∈I σ 1,j i∈I (u 1,i+1 u 1,i+1+j -u 1,i u 1,i+j ) (log(u 1,i+1 ) -log(u 1,i )) = 1 2 J σ 1 1 + d 11 2d 12 j∈I σ 1,j i∈I (u 1,i+1-j u 1,i+1 -u 1,i-j u 1,i ) (log(u 1,i+1-j ) -log(u 1,i-j )) = d 11 2d 12 j∈I σ 1,j i∈I (u 1,i+1 u 1,i+1-j -u 1,i u 1,i-j ) (log(u 1,i+1 u 1,i+1-j ) -log(u 1,i u 1,i-j )) ≥ 2d 11 d 12 j∈I σ 1,j i∈I √ u 1,i+1 u 1,i+1-j - √ u 1,i u 1,i-j 2 .
In the previous estimate, the second inequality is obtained by changing j into -j and using the symmetry of σ 1 . For the third equality, one changes i into i-j. The fouth one is the combination of the first and third equalities. Once again a similar estimate holds for J σ 2 2 . Finally with the same changes of indices one can estimate the sum

J ρ 1 1 + J ρ 2 2 = j∈I ρ j i∈I (u 1,i+1 u 2,i+1-j -u 1,i u 2,i-j ) (log(u 1,i+1 ) -log(u 1,i )) + j∈I ρ -j i∈I (u 2,i+1 u 1,i+1-j -u 2,i u 1,i-j ) (log(u 2,i+1 ) -log(u 2,i )) = j∈I ρ j i∈I (u 1,i+1 u 2,i+1-j -u 1,i u 2,i-j ) (log(u 1,i+1 ) -log(u 1,i )) + j∈I ρ j i∈I (u 2,i-j+1 u 1,i+1 -u 2,i-j u 1,i ) (log(u 2,i-j+1 ) -log(u 2,i-j )) = j∈I ρ j i∈I (u 1,i+1 u 2,i+1-j -u 1,i u 2,i-j ) (log(u 1,i+1 u 2,i-j+1 ) -log(u 1,i u 2,i-j )) ≥ j∈I ρ j i∈I √ u 1,i+1 u 2,i+1-j - √ u 1,i u 2,i-j 2 .
By summing all the estimates one obtains [START_REF] Evans | Partial differential equations[END_REF]. The last point of the proposition is obtained by induction.

3.2. Proof of Theorem 1. Let us show that P k • Φ has at least one zero. We use an approximation argument by introducing 24) and the non-negativity of the entropy and the entropy dissipation one has

P k ε (X) = P k (X) + εX, ∀X ∈ R 2N . Using (
P k ε (X), X ≥ ε|X| 2 -H(u k-1 1 , u k-1 2 )
. Therefore, as a consequence of Brouwer fixed point theorem (see [START_REF] Evans | Partial differential equations[END_REF]Section 9.1] for details), there is X ε such that

P k ε (X ε ) = 0 and |X ε | 2 ≤ H(u k-1 1 , u k-1 2 )ε -1 . Let us define the associated (u ε 1 , u ε 2 ) = Φ -1 (X ε ), which is componentwise positive by definition. Observe that ∆x h j (u ε j,i ) ≤ H(u ε 1 , u ε 2 ) ≤ H(u k-1 1 , u k-1 2 ), ∀i ∈ I, j = 1, 2.
where the last inequality is again a consequence of (24) for X = X ε . This shows that for any (u ε 1,i , u ε 2,i ) is uniformly bounded in ε. Therefore, there exists a subsequence (not relabeled) such that u ε j,i → u k j,i ≥ 0 as ε → 0, for every i ∈ I and j = 1, 2. Since

|X ε | = O(ε -1/2 ) one has 0 = lim ε→0 P k ε (X ε ) = P k (Φ(u k 1 , u k 2 ))
Therefore (u k 1 , u k 2 ) solves the scheme ( 11)- [START_REF] Chainais-Hillairet | Large-time behaviour of a family of finite volume schemes for boundarydriven convection-diffusion equations[END_REF]. By taking limits in ( 23) and ( 24) evaluated at X ε as ε → 0 one recovers ( 21) and ( 22) respectively.

Let us finally prove that if d j > 0, u k j,i > 0 for all i ∈ I and 0 < k ≤ N T . This is a consequence of the entropy estimate. For a given 0 < k ≤ N T we notice that (thanks to the term J diff j ) the positive solution (u ε 1 , u ε 2 ) satisfies the following estimate

d j ∆t ∆x i∈I u ε j,i+1 -u ε j,i log u ε j,i+1 -log u ε j,i ≤ max(d 12 , d 21 ) H u k-1 1 , u k-1 2 j = 1, 2.
At the limit ε → 0, let us assume by contradiction that there exists i ∈ I such that u k j,i = 0. Then as the r.h.s. of the previous inequality is finite this implies that u k j,i+1 = 0. Thus repeating this argument we deduce that u k j,i = 0 for all i ∈ I. Consequently we have u k j L 1 (T) = 0 which contradicts the mass conservation property [START_REF] Dietert | Persisting entropy structure for nonlocal cross-diffusion systems[END_REF]. This completes the proof of Theorem 1.

Estimates on the discrete Kolmogorov equation

In this section, we focus on estimates concerning the finite volume discretization of the Kolmogorov equation ∂ t z = ∆(µz). In particular we adapt at the discrete level some properties established in [START_REF] Moussa | From nonlocal to classical Shigesada-Kawasaki-Teramoto systems: triangular case with bounded coefficients[END_REF][START_REF] Dietert | Persisting entropy structure for nonlocal cross-diffusion systems[END_REF].

In the rest of this section, we assume that (µ k i ) i∈I , k = 1, . . . , N T is given and componentwise non-negative. From there, the scheme is given for all k ≥ 1 by

z k i -z k-1 i ∆t -∆ T (µ k z k ) i = 0, ∀i ∈ I, (25)
where ∆ T denotes the discrete Laplacian operator defined by (13).

4.1. Well-posedness of the scheme and L ∞ estimates. Let us first prove that the scheme (25) admits a unique solution at each time step. Lemma 6. For any

(z k-1 i ) i∈I there is a unique (z k i ) i∈I satisfying (25). Moreover, if z k-1 is componentwise nonnegative then so is z k . Proof. Let us write Z k-1 = (z k-1 0 , . . . , z k-1 N -1
) for all k ≥ 1. Observe that the scheme writes

M k Z k = Z k-1 where M k is a N × N tridiagonal matrix defined by M k i,i-1 = - ∆t ∆x 2 µ k i-1 , M k i,i = 1 + 2 ∆t ∆x 2 µ k i , M k i,i+1 = - ∆t ∆x 2 µ k i+1 , ∀i ∈ I. ( 26
)
We notice that M k has positive diagonal terms and non-positive off-diagonal terms. Furthermore the matrix M k is strictly diagonally dominant with respect to its columns. Therefore M k is a nonsingular M-matrix and is thus monotone and invertible. This finishes the proof of Lemma 6.

We prove in the following result some L ∞ estimates for the solution to scheme [START_REF] Filbet | A finite volume scheme for boundary-driven convection-diffusion equations with relative entropy structure[END_REF]. Lemma 7. Let us assume that there exists γ, Γ ≥ 0 such that γ ≤ z 0 i ≤ Γ, ∀i ∈ I. Then for every k ≥ 1 and every ∆t > 0 such that

∆t < 1/ max 1≤k≤N T [∆ T µ k ] + L ∞ (T) , the solution Z k to (25) satisfies γ Π k n=1 1 + ∆t [∆ T µ n ] -L ∞ (T) -1 ≤ z k i ≤ Γ Π k n=1 1 -∆t [∆ T µ n ] + L ∞ (T) -1 , ∀i ∈ I, (27)
where [x] + = max(x, 0) and [x] -= min(x, 0).

Proof. We will only deal with the upper bound in [START_REF] Galiano | Semi-discretization in time and numerical convergence of solutions of a nonlinear cross-diffusion population model[END_REF] and the lower bound is obtained in the same way. Let M k denotes the tridiagonal matrix defined by [START_REF] Fontbona | Non local Lotka-Volterra system with cross-diffusion in an heterogeneous medium[END_REF] and define

Γk = Γ Π k n=1 1 -∆t [∆ T µ n ] + L ∞ (T) -1 .
We proceed by induction. Since Γ0 = Γ the bound holds by hypothesis at k = 0. Then observe that for every i ∈ I

(M k (Z k -Γk )) i = Z k-1 -Γk + ∆t Γk (∆ T µ k ) i ≤ Γk-1 -Γk + ∆t Γk (∆ T µ k ) i .

Now we notice that by construction

Γk-1 -Γk = -∆t Γk [∆ T µ k ] + L ∞ (T) .
Then we easily deduce that for every i ∈ I it holds

(M k (Z k -Γk )) i = -∆t Γk [∆ T µ k ] + L ∞ (T) -(∆ T µ k ) i ≤ 0.
Therefore, since M k is a M-matrix we conclude that z k i ≤ Γk for all i ∈ I which concludes the proof of Lemma 7.

The bounds of Lemma 7 are exactly the discrete equivalent of the L ∞ estimates established at the continuous level in [21, Corollary 18]. Proposition 8. Let us assume that it holds

∆t < 1/ max 1≤k≤N T [∆ T µ k ] + L ∞ (T) .
Then the solution to (25) satisfies the following estimate

z k 2 L 2 (T) + k n=1 ∆t i∈T (µ n i + µ n i+1 ) (z n i+1 -z n i ) 2 ∆x ≤ Π k n=1 1 -∆t [∆ T µ n ] + L ∞ (T) -1 z 0 i 2 L 2 (T) , ∀1 ≤ k ≤ N T . Proof.
Let k ≥ 1 be fixed and let us first notice that we can rewrite for every i ∈ I equation [START_REF] Filbet | A finite volume scheme for boundary-driven convection-diffusion equations with relative entropy structure[END_REF] as

∆x z k i -z k-1 i ∆t + µ k i+ 1 2 (z k i -z k i+1 ) ∆x -µ k i-1 2 (z k i-1 -z k i ) ∆x + z k i+ 1 2 (µ k i -µ k i+1 ) ∆x -z k i-1 2 (µ k i-1 -µ k i ) ∆x = 0,
where

µ k i+ 1 2 = µ k i + µ k i+1 2 , z k i+ 1 2 = z k i + z k i+1 2 , ∀i ∈ I.
Now we multiply the above equation by ∆tz k i and we sum over i ∈ I, we obtain I 3 + I 4 + I 5 = 0, where

I 3 = i∈I ∆x(z k i -z k-1 i )z k i , I 4 = ∆t i∈I µ k i+ 1 2 (z k i -z k i+1 ) ∆x -µ k i-1 2 (z k i-1 -z k i ) ∆x z k i , I 5 = ∆t i∈I z k i+ 1 2 (µ k i -µ k i+1 ) ∆x -z k i-1 2 (µ k i-1 -µ k i ) ∆x z k i .
For I 3 using the inequality (a -b)a ≥ (a 2 -b 2 )/2 we obtain

I 3 ≥ 1 2 i∈I ∆x |z k i | 2 -|z k-1 i | 2 . ( 28
)
For I 4 applying a discrete integration by parts yields

I 4 = ∆t i∈I µ k i+ 1 2 (z k i+1 -z k i ) 2 ∆x . ( 29
)
Now we rewrite I 5 as

I 5 = - ∆t 2 i∈I ∆x|z k i | 2 (∆ T µ k ) i + ∆t 2∆x i∈I z k i+1 z k i (µ k i -µ k i+1 ) -z k i-1 z k i (µ k i-1 -µ k i ) ,
and reordering the terms in the r.h.s. the second sum vanishes and we have

I 5 = - ∆t 2 i∈I ∆x|z k i | 2 (∆ T µ k ) i . (30)
Gathering ( 28)-( 30) we end up with

1 2 i∈I ∆x|z k i | 2 + ∆t i∈I µ k i+ 1 2 (z k i+1 -z k i ) 2 ∆x ≤ 1 2 i∈I ∆x|z k-1 i | 2 + ∆t 2 [∆ T µ k ] + L ∞ (T) i∈I ∆x|z k i | 2 .
We deduce that

1 2 i∈I ∆x|z k i | 2 + k n=1 ∆t i∈I µ n i+ 1 2 (z n i+1 -z n i ) 2 ∆x ≤ 1 2 i∈I ∆x|z 0 i | 2 + k n=1 ∆t 2 [∆ T µ n ] + L ∞ (T) i∈I ∆x|z n i | 2 .
It remains to apply the discrete Grönwall inequality, see Lemma 15 in Appendix A, in order to complete the proof of Proposition 8.

Study of the dual problem.

The main objective of this section is to establish a discrete counterpart of the so-called duality inequality for the solution to [START_REF] Filbet | A finite volume scheme for boundary-driven convection-diffusion equations with relative entropy structure[END_REF], see for instance [START_REF] Moussa | From nonlocal to classical Shigesada-Kawasaki-Teramoto systems: triangular case with bounded coefficients[END_REF]Theorem 3]. In this aim, following [START_REF] Moussa | From nonlocal to classical Shigesada-Kawasaki-Teramoto systems: triangular case with bounded coefficients[END_REF], we introduce a "dual" scheme associated to [START_REF] Filbet | A finite volume scheme for boundary-driven convection-diffusion equations with relative entropy structure[END_REF]. Let v N T +1 i be given for every i ∈ I, then for 1 ≤ k ≤ N T we want to determine the solution to the following implicit backward in time scheme

v k i -v k+1 i ∆t -µ k i ∆ T v k i = S k i ∀i ∈ I, (31)
where µ k i is given and non-negative and S k = (S k 0 , . . . , S k N -1 ) is some given vector in R N for all 1 ≤ k ≤ N T . Let us notice that (31) define a set of linear equation which can be rewritten as

(M k ) V k = V k+1 + ∆tS k , ∀1 ≤ k ≤ N T , (32)
where M k is the tridiagonal matrix given by [START_REF] Fontbona | Non local Lotka-Volterra system with cross-diffusion in an heterogeneous medium[END_REF]. Therefore, it follows directly from the proof of Lemma 6 that the problem (32) admits a unique solution for every 1 ≤ k ≤ N T .

Prior to the proof of the discrete duality estimate, see Theorem 10 below, we establish some uniform estimates satisfied by the solution of (32). Proposition 9. Assume that min i∈I µ k i > 0 for every 0 ≤ k ≤ N T and that v N T +1 i = 0 for every i ∈ I. Then the solution to (32) satisfies for every 1 ≤ k ≤ N T the following estimate

|v k | 2 1,2,T + N T n=k ∆t i∈I µ n i (∆ T v n ) 2 i ∆x ≤ µ -1/2 S 2 L 2 (Q T ) , (33)
and there exists a constant C > 0 independent of ∆x such that

v k 2 L 2 (T) ≤ C(1 + µ L 1 (Q T ) ) µ -1/2 S 2 L 2 (Q T ) , ∀1 ≤ k ≤ N T , (34)
where µ and S denote the piecewise reconstruction functions in H D associated to the vectors (µ k ) 1≤k≤N T and (S k ) 1≤k≤N T .

Proof. Let us first establish estimate [START_REF] Kruzhkov | Results concerning the nature of the continuity of solutions of parabolic equations and some of their applications[END_REF]. In this purpose let 1 ≤ k ≤ N T be fixed. We multiply equation ( 31) by ∆t(-

v k i+1 + 2v k i -v k i-1
)/∆x, we sum over i ∈ I and we apply definition (13) of the operator ∆ T and we obtain 3 ,

I 6 + I 7 = I 8 , with I 6 = i∈I (v k i -v k+1 i ) (-v k i+1 + 2v k i -v k i-1 ) ∆x , I 7 = ∆t i∈I µ k i (v k i+1 -2v k i + v k i-1 ) 2 ∆x
I 8 = ∆t i∈I S k i (-v k i+1 + 2v k i -v k i-1 ) ∆x .
For I 6 reordering the terms leads to

I 6 = 1 ∆x i∈I (v k i+1 -v k i ) (v k i+1 -v k i ) -(v k+1 i+1 -v k+1 i ) ,
and using the inequality a(a -b) ≥ (a 2 -b 2 )/2 we get

I 6 ≥ 1 2 i∈I (v k i+1 -v k i ) 2 ∆x - (v k+1 i+1 -v k+1 i ) 2 ∆x . ( 35
)
For I 8 applying the Cauchy-Schwarz and Young inequality yield

|I 8 | ≤ ∆t 2 i∈I ∆x(µ k i ) -1 |S k i | 2 + ∆t 2 i∈I µ k i (-v k i+1 + 2v k i -v k i-1 ) 2 ∆x 3 . ( 36
)
Collecting ( 35)-( 36) we obtain

i∈I (v k i+1 -v k i ) 2 ∆x + ∆t i∈I µ k i (∆ T v k ) 2 i ∆x ≤ i∈I (v k+1 i+1 -v k+1 i ) 2 ∆x + ∆t i∈I ∆x(µ k i ) -1 |S k i | 2 .
In order to prove (33) it remains to sum over n ∈ {k, . . . , N T }.

We now prove estimate [START_REF] Lepoutre | Entropic structure and duality for multiple species cross-diffusion systems[END_REF]. In this purpose we multiply [START_REF] Jüngel | Nonlocal cross-diffusion systems for multi-species populations and networks[END_REF] by ∆x∆t, we sum over i ∈ I and n ∈ {k, . . . , N T } and we obtain i∈I

∆xv k i = N T n=k ∆t i∈I µ k i v n i+1 -2v n i + v n i-1 ∆x + N T n=k ∆t i∈I ∆xS n i .
Applying the Cauchy-Schwarz inequality leads to

i∈I ∆xv k i ≤ N T n=k ∆t i∈I ∆xµ n i 1/2 N T n=k ∆t i∈I µ n i (∆ T v n ) 2 i ∆x 1/2 + N T n=k ∆t i∈I ∆x(µ n i ) -1 |S n i | 2 1/2
.

Using estimate (33) we obtain

i∈I ∆xv k i ≤ 2 N T n=k ∆t i∈I ∆xµ n i 1/2 N T n=k ∆t i∈I ∆x(µ n i ) -1 |S n i | 2 1/2
. Now it remains to apply the discrete Poincaré-Wirtinger inequality on the torus obtained in [6, Lemma 6] in order to conclude the proof of Proposition 9.

We are now in position to establish the discrete dual estimate.

Theorem 10. Let us assume that min i∈I µ k i > 0 for every 1 ≤ k ≤ N T . Then there exists a constant C > 0 independent of ∆x such that the solution (Z k ) 1≤k≤N T to (25) satisfies

µ 1/2 z L 2 (Q T ) ≤ C 1 + µ 1/2 L 1 (Q T ) z 0 L 2 (T) .
Proof. Let (v k i ) i∈I be given in R N for every 1 ≤ k ≤ N T + 1 with v N T +1 i = 0 for all i ∈ I. Now for 1 ≤ k ≤ N T , we multiply (25) by ∆t∆xv k i , we sum over i ∈ I and k ∈ {1, . . . , N T }, we obtain

N T k=1 i∈I ∆x(z k i -z k-1 i )v k i - N T k=1 ∆t i∈I ∆x ∆ T (z k µ k ) i v k i = 0.
Reordering the terms we have

N T k=1 ∆t i∈I ∆xz k i (v k i -v k+1 i ) ∆t -µ k i ∆ T v k i = i∈I ∆xz 0 i v 1 i . ( 37
)
We define (S k i ) i∈I by

S k i = (v k i -v k+1 i ) ∆t -µ k i ∆ T v k i , ∀i ∈ I, 1 ≤ k ≤ N T .
We first notice that (S k i ) i∈I is well-defined since we know that equation ( 31) is well-posed. Besides applying the Cauchy-Schwarz inequality in [START_REF] Moussa | From nonlocal to classical Shigesada-Kawasaki-Teramoto systems: triangular case with bounded coefficients[END_REF] we get

N T k=1 ∆t i∈I ∆xz k i S k i ≤ z 0 L 2 (T) v 1 L 2 (T) .
Now, thanks to [START_REF] Lepoutre | Entropic structure and duality for multiple species cross-diffusion systems[END_REF] we deduce that

N T k=1 ∆t i∈I ∆x(µ k i ) 1/2 z k i (µ k i ) -1/2 S k i ≤ C 1 + µ 1/2 L 1 (Q T ) µ -1/2 S L 2 (Q T ) z 0 L 2 (T) .
In the remaining of the proof we want to use the dual definition [START_REF] Cancès | Finite-volume scheme for a degenerate cross-diffusion model motivated from ion transport[END_REF] of the norm

• L 2 (Q T ) .
Observe that for any vector

F k = (f k i ) i∈I , there exists a unique V = (v i ) i∈I such that (µ k i ) -1/2 v i -v k+1 i ∆t -(µ k i ) 1/2 (∆ T v) i = f k i , ∀i ∈ I,
where (v k+1 i ) i∈I is a given vector. Indeed, this system rewrites (M k ) V = V k+1 + ∆tD k F k , for all 1 ≤ k ≤ N T where M k is the invertible tridiagonal matrix given by ( 26) and D k = diag((µ k i )) i∈I . We deduce thanks to formula (10) that it holds

N T k=1 ∆t i∈I ∆xµ k i |z k i | 2 1/2 ≤ C 1 + µ 1/2 L 1 (Q T ) z 0 L 2 (T) .
This concludes the proof of Theorem 10.

Proof of Theorem 2.

We are now able to prove Theorem 2.

Step 1: Maximum principle. Let us first prove the maximum principle satisfies by the solutions to ( 11)- [START_REF] Chainais-Hillairet | Large-time behaviour of a family of finite volume schemes for boundarydriven convection-diffusion equations[END_REF]. Let us notice that for every 1 ≤ k ≤ N T we have

max i∈I ∆ T µ k 1 i = max i∈I d 11 j∈I ∆xu k 1,j (∆ T σ 1 ) i-j + d 12 j∈I ∆xu k 2,j (∆ T ρ) i-j ≤ d 11 ∆ T σ 1 L ∞ (T) j∈I ∆xu k 1,j + d 12 ∆ T ρ L ∞ (T) j∈I ∆x u k 2,j .
Now, let us recall that m 0 j = u 0 j L 1 (T) for j = 1, 2, then thanks to the mass conservation property (21) we obtain

max i∈I ∆ T µ k 1 i ≤ d 11 m 0 1 ∆ T σ 1 L ∞ (T) + d 12 m 0 2 ∆ T ρ L ∞ (T) .
Similarly we establish the following bound

max i∈I ∆ T µ k 2 i ≤ d 22 m 0 2 ∆ T σ 2 L ∞ (T) + d 21 m 0 1 ∆ T ρ L ∞ (T) .
As a direct consequence of the previous estimates and ( 27) (with γ = γ and Γ = Γ) one obtains point (i) of Theorem 2.

Step 2: Duality estimate. Let us now show the discrete duality estimate satisfied by the solutions to ( 11)- [START_REF] Chainais-Hillairet | Large-time behaviour of a family of finite volume schemes for boundarydriven convection-diffusion equations[END_REF]. For every 0 ≤ k ≤ N T we define the element

z k i = u k 1,i + u k 2,i for all i ∈ I. Observe that z k i is solution to z k i -z k-1 i ∆t + ∆ T (µ k z k ) i = 0, where µ k i = µ k 1,i u k 1,i + µ k 2,i u k 2,i u k 1,i + u k 2,i
, ∀i ∈ I.

Thanks to Theorem 1, we have u k 1,i , u k 2,i > 0 for all i ∈ I and the element µ k i is well-defined. Besides, applying the discrete duality estimate established in Theorem 10 we deduce the existence of a constant C > 0 independent of ∆x such that (38)

N T k=1 ∆t i∈I ∆x µ k 1,i u k 1,i + µ k 2,i u k 2,i u k 1,i + u k 2,i ≤ C 1 + N T k=1 ∆t i∈I ∆x|µ k i | i∈I ∆x|u 0 1,i | 2 + i∈I ∆x|u 0 2,i | 2 .

Now we notice that

N T k=1 ∆t i∈I ∆x|µ k i | ≤ N T k=1 ∆t i∈I ∆x|µ k 1,i | + N T k=1 ∆t i∈I ∆x|µ k 2,i | = I 9 + I 10 .
For I 9 we have

I 9 = N T k=1 ∆t i∈I ∆x d 1 + d 11 j∈I ∆xσ 1,i-j u k 1,j + d 12 j∈I ∆xρ i-j u k 2,j = N T k=1 ∆t d 1 + d 11 j∈I ∆xu k 1,j i∈I ∆xσ 1,i-j + d 12 j∈I ∆xu k 2,j i∈I ∆xρ i-j .
Thus, bearing in mind the mass conservation property ( 21) we obtain

I 9 ≤ T (d 1 + d 11 m 0 1 σ 1 L 1 (T) + d 12 m 0 2 ρ L 1 (T)
), [START_REF] Painter | A nonlocal model for contact attraction and repulsion in heterogeneous cell populations[END_REF] and similarly

I 10 ≤ T (d 2 + d 22 m 0 2 σ 2 L 1 (T) + d 21 m 0 1 ρ L 1 (T) ). ( 40
)
Collecting ( 38)- [START_REF] Potts | Spatial memory and taxis-driven pattern formation in model ecosystems[END_REF] we conclude that point (ii) of Theorem 2 holds.

Convergence of the scheme

This section is dedicated to the proof of Theorem 3. In the following the subscript m refer to the size η m = max{∆x m , ∆t m } of the family (D m ) of space-time discretizations of Q T . We derive uniform in m a priori estimates in subsection 5.1 in order to obtain compactness in L p (Q T ) of the sequences of constant by part reconstructions (u j,m ) m for both species j = 1, 2. The compactness results are gathered in Section 5.2. A keypoint is a discrete L 1 compactness result obtained in [START_REF] Andreianov | Analysis of a finite volume method for a cross-diffusion model in population dynamics[END_REF]Lemma 9.2]. This result is the adaptation at the discrete level of a compactness lemma established by Kruzhkov in [START_REF] Kruzhkov | Results concerning the nature of the continuity of solutions of parabolic equations and some of their applications[END_REF]. Finally in Section 5.3, we prove Theorem 3.

Uniform estimates.

In this section we establish some uniform estimates w.r.t. ∆x and ∆t fulfilled by the solutions to the scheme ( 11)- [START_REF] Chainais-Hillairet | Large-time behaviour of a family of finite volume schemes for boundarydriven convection-diffusion equations[END_REF]. They rely on the entropy dissipation inequality [START_REF] Domschke | Mathematical modelling of cancer invasion: implications of cell adhesion variability for tumour infiltrative growth patterns[END_REF] and the conservation of mass ( 21

(u 0 1 , u 0 2 ) such that max k=1,...,N T u k j L 1 (T) + d j N T k=1 ∆t u k j 2 1,1,T 1 2 ≤ C 1 , for j = 1, 2. (41)
Moreover, assuming that d 1 and d 2 are positive constants, there exists a constant C 2 > 0 only depending on T ,

d 1 , d 2 , d 11 , d 12 , d 21 , d 22 , σ 1 L ∞ (T) , σ 2 L ∞ (T) , ρ L ∞ (T) , m 0 1 , m 0 2 and H(u 0 1 , u 0 2 ) such that N T k=1 ∆t i∈I ∆x F k j,i+ 1 2 ≤ C 2 , for j = 1, 2, ( 42
)
where the numerical fluxes are defined by [START_REF] Chen | Analysis of a multidimensional parabolic population model with strong cross-diffusion[END_REF].

Proof. The uniform L ∞ (0, T ; L 1 (T)) estimate of the first term in the right hand side of ( 41) is a direct consequence of the conservation of mass [START_REF] Dietert | Persisting entropy structure for nonlocal cross-diffusion systems[END_REF]. Then, for the uniform discrete L 2 (0, T ; W 1,1 (T)) estimate, we first notice, for j = 1 or 2 and k ∈ {1, . . . , N T }, that it holds

|u k j | 1,1,T = i∈I u k j,i+1 -u k j,i = i∈I u k j,i+1 -u k j,i u k j,i+1 + u k j,i
.

Hence, the Cauchy-Schwarz inequality yields

|u k j | 1,1,T ≤ u k j 1,2,T i∈I ∆x u k j,i+1 + u k j,i 2 1 2
Since (a + b) 2 ≤ 2(a 2 + b 2 ) and u k j L 1 (T) = u 0 j L 1 (T) = m 0 j (conservation of mass), one has

|u k j | 1,1,T ≤ 2 m 0 j 1/2 u k j 1,2,T
.

Therefore, applying the entropy inequality [START_REF] Domschke | Mathematical modelling of cancer invasion: implications of cell adhesion variability for tumour infiltrative growth patterns[END_REF], we get for the first species

d 1 N T k=1 ∆t |u k 1 | 2 1,1,T ≤ 4d 1 m 0 1 N T k=1 ∆t u k 1 2 1,2,T ≤ d 12 m 0 1 H(u 0 1 , u 0 2 ),
and the equivalent estimate holds for the second species. This yields the existence of C 1 such that (41) holds. It remains to establish [START_REF] Sun | An entropy stable high-order discontinuous Galerkin method for cross-diffusion gradient flow systems[END_REF]. In this purpose we will consider the case j = 1. Then, using the definition (18) of the numerical fluxes, we estimate

N T k=1 ∆t i∈I ∆x F k 1,i+ 1 2 ≤ N T k=1 ∆t i∈I µ k 1,i+ 1 2 u k 1,i -u k 1,i+1 + N T k=1 ∆t i∈I u k 1,i+ 1 2 µ k 1,i -µ k 1,i+1 = I 11 + I 12 .
For I 11 , applying the regularity of the functions σ 1 and ρ we have

I 11 ≤ N T k=1 ∆t i∈I d 1 + d 11 σ 1 L ∞ (T) u k 1 L 1 (T) + d 12 ρ L ∞ (T) u k 2 L 1 (T) u k 1,i -u k 1,i+1 .
Hence, using the conservativity of the scheme and ( 41), we get

I 11 ≤ d 1 + d 11 m 0 1 σ 1 L ∞ (T) + d 12 m 0 2 ρ L ∞ (T) T 1/2 C 1 . (43)
For I 12 , using the definition of µ k 1,i for i ∈ I, we notice that it holds

I 12 ≤ 1 2 N T k=1 ∆t i∈I ∆x u k 1,i + u k 1,i+1 × d 11 n∈I σ 1,n u k 1,i-n -u k 1,i+1-n + d 12 n∈I ρ n u k 2,i-n -u k 2,i+1-n .
Then, thanks to the conservativity of the scheme, we obtain

I 12 ≤ 2m 0 1 N T k=1 ∆t d 11 σ 1 L ∞ (T) u k 1 1,1,T + d 12 ρ L ∞ (T) u k 2 1,1,T .
Therefore, applying the Cauchy-Schwarz inequality and ( 41) we end up with

I 12 ≤ 2 m 0 1 T 1/2 C 1 d 11 d 1/2 1 σ 1 L ∞ (T) + d 12 d 1/2 2 ρ L ∞ (T) . (44)
Collecting ( 43) and ( 44) and the corresponding inequalities for the second species lead to the existence of C 2 such that [START_REF] Sun | An entropy stable high-order discontinuous Galerkin method for cross-diffusion gradient flow systems[END_REF] holds. This concludes the proof of Proposition 11.

Compactness properties.

Let (u 1,m , u 2,m ) m∈N be a family, constructed in Theorem 1, of finite volume solutions to ( 11)-( 15) associated to the sequence (D m ). In order to be able to apply [START_REF] Andreianov | Analysis of a finite volume method for a cross-diffusion model in population dynamics[END_REF]Lemma 9.2], the first task is to rewrite the scheme ( 11)- [START_REF] Chainais-Hillairet | Large-time behaviour of a family of finite volume schemes for boundarydriven convection-diffusion equations[END_REF] as the discretization of an evolution equation under divergence form. In this purpose we use the equivalent form of [START_REF] Carrillo | A finite-volume method for nonlinear nonlocal equations with a gradient flow structure[END_REF] given by [START_REF] Chen | Rigorous derivation of population cross-diffusion systems from moderately interacting particle systems[END_REF]. In particular, for j = 1, 2 and k = 1, . . . , N T , we associate to the family of fluxes F k j,i+1/2 i∈I the following piecewise reconstruction

F k j,m = i∈Im F k j,i+ 1 2 1 (x i ,x i+1 ) .
Then, for this discrete field F k j,m we define its L 1 norm as

F k j,m L 1 (T) = i∈Im ∆x F k j,i+ 1 2
, and its discrete divergence by

div T F k j,m i = 1 ∆x F k j,i+ 1 2 -F k j,i-1 2 , i ∈ I m .
This definition allows us to rewrite (18) as

u k j,i -u k-1 j,i ∆t + div T F k j,m i = 0, ∀i ∈ I m , j = 1, 2, (45)
and we obtain the following result: Proposition 12. Let the assumptions of Theorem 3 hold and let (u 1,m , u 2,m ) m∈N be a sequence of discrete solutions to (11)-( 15) constructed in Theorem 1. Then there exists a subsequence of (u 1,m , u 2,m ), which is not relabeled, and

(u 1 , u 2 ) ∈ (L p (Q T )) 2 , with p ∈ [1, 3), such that u j,m → u j strongly in L p (Q T ) for 1 ≤ p < 3, as m → ∞,
and almost everywhere.

Proof. A direct consequence of Proposition 11 is that there is a constant C independent of ∆x m and ∆t m such that

N T k=1 ∆t m u k j,m L 1 (T) + N T k=1 ∆t m F k j,m L 1 (T) + N T k=1 ∆t m |u k j,m | 1,1,Tm ≤ C, j = 1, 2.
By [4, Lemma 9.2], which can be applied thanks to (45), there is function

u j ∈ L 1 (Q T ), j = 1, 2,
such that, up to a subsequence,

u j,m → u j strongly in L 1 (Q T ) as m → ∞.
Moreover, Proposition 11 also implies that the sequence (u i,m ) is uniformly bounded in the space L ∞ (0, T ; L 1 (T)) and in L 2 (0, T ; BV (T)). The continuous embedding of BV (T) in L ∞ (T) (see [START_REF] Ambrosio | Functions of bounded variation and free discontinuity problems[END_REF]) implies that the sequence (u j,m ) is uniformly bounded in L 2 (0, T ; L ∞ (T)). Hence, by interpolation, one has a uniform bound of (u j,m ) in L 3 (Q T ). Thus, Vitali's theorem gives the strong convergence of (u j,m ) towards u j in L p (Q T ) for all p ∈ [1, 3). This concludes the proof of Proposition 12.

Remark 13. In dimension d ≥ 2, we have the compact embedding of the space

BV (T d ) in L d d-1 (T d ).
In particular in this case the sequence (u j,m ) is uniformly bounded in L d d-1 (Q T ). Therefore arguing as in the previous proof we deduce the existence for j = 1 and 2 of u j ∈ L p (Q T ) for p ∈ [1, d/(d -1)), such that, up to a subsequence,

u j,m → u j strongly in L p (Q T ) for 1 ≤ p < d d -1
, as m → ∞.

Corollary 14. Let the assumptions of Proposition 12 hold. Then there exists a subsequence of (u 1,m , u 2,m ), such that for any p ∈ [1, 3) one has

µ 1,m → µ 1 = d 1 + d 12 ρ 1 * u 2 + d 11 σ 1 * u 1 strongly in L p (0, T ; L ∞ (T)), as m → ∞, µ 2,m → µ 2 = d 2 + d 21 ρ 2 * u 1 + d 22 σ 2 * u 2 strongly in L p (0, T ; L ∞ (T)), as m → ∞,
where we recall that ρ 1 (x) = ρ 2 (-x) = ρ(x) for a.e. x ∈ T.

Proof. Observe that by definition ( 14)-( 16) one has for x ∈ K i and t ∈ (0, T )

µ 1 (x, t) -µ 1,m (x, t) = d 12 (ρ 1 * (u 2 -u 2,m ))(x, t) + d 11 (σ 1 * (u 1 -u 1,m ))(x, t) + d 12 T (ρ 1 (x -y) -ρ 1 (x i -y)) u 2,m (y, t) dy + d 11 T (σ 1 (x -y) -σ 1 (x i -y)) u 1,m (y, t) dy.
Therefore by dominated convergence one obtains the strong convergence of µ 1,m towards µ 1 in L 1 (Q T ) and almost everywhere in Q T . Besides, for a.e. t ∈ (0, T ) thanks to Young's inequality and for p ∈ (1, 3), it holds

ρ 1 * (u 2 -u 2,m )(t) L ∞ (T) ≤ ρ 1 L p p-1 (T) (u 2 -u 2,m )(t) L p (T) .
Then, applying Proposition 12, we obtain

T 0 ρ 1 * (u 2 -u 2,m )(t) p L ∞ (T) dt ≤ ρ 1 p L p p-1 (T) T 0 (u 2 -u 2,m )(t) p L p (T) dt → 0 as m → ∞.
Let us now setting ξ(x, y) = ρ 1 (x -y) -ρ 1 (x i -y) for a.e. x ∈ K i and y ∈ T. Hence, for a.e. t ∈ (0, T ), the Hölder inequality yields for p ∈ (1, 3)

T ξ(•, y)u 2,m (t, y)dy L ∞ (T) ≤ sup |z|≤∆xm ρ 1 (z + •) -ρ 1 L p p-1 (T) u 2,m (t) L p (T) .
The first factor in the right hand side tends to 0 (by density of continuous functions in L p/(p-1) (T)) while, bearing in mind Proposition 12, the second factor is uniformly bounded in L p (Q T ). Therefore one can conclude the strong convergence in L p (0, T ; L ∞ (T)) by using Young's inequality and the previous argument. This finishes the proof of Corollary 14.

5.3.

Proof of Theorem 3. It remains to prove that the functions u 1 and u 2 constructed in Section 5.2 are distributional solutions to (1)-( 3) in the sense of [START_REF] Burger | Segregation effects and gap formation in cross-diffusion models[END_REF]. Let φ ∈ C ∞ 0 (T × [0, T )), we multiply equation ( 12) by ∆t m ∆x m φ k-1 i , where φ k-1 i = φ(x i , t k-1 ), and we sum over i ∈ I and k ∈ {1, . . . , N T }. We obtain

F m 1 + F m 2 = 0 with F m 1 = N T k=1 i∈I ∆x m (u k 1,i -u k-1 1,i )φ k-1 i = N T k=1 i∈I ∆x m u k 1,i (φ k-1 i -φ k i ) - i∈I ∆x m u 0 1,i φ 0 i , and
F m 2 = N T k=1 ∆t m i∈I -µ k 1,i+1 u k 1,i+1 + 2µ k 1,i u k 1,i -µ k 1,i-1 u k 1,i-1 ∆x m φ k-1 i = N T k=1 ∆t m i∈I -φ k-1 i+1 + 2φ k-1 i -φ k-1 i-1 ∆x m µ k 1,i u k 1,i . Let ψ m (x, t) = (φ k i -φ k-1 i )/∆t m for all x ∈ K i and t ∈ (t k-1 , t k ] and φ 0 m (x) = φ 0 i for all x ∈ K i . Then, since ψ m → ∂ t φ in L ∞ (Q T ) and φ 0 m → φ 0 in L ∞ (T)
, so using the convergence results of Proposition 12 one obtains

F m 1 + Q T u 1 ∂ t φ dxdt + T u 0 1 (x)φ(x, 0) dx = Q T (u 1 ∂ t φ -u 1,m ψ m ) dxdt + T (u 0 1 φ(•, 0) -u 0 1 φ 0 m ) dx → 0 as m → 0. Similarly, if one defines ζ m (x, t) = (φ k-1 i+1 -2φ k-1 i + φ k-1 i-1 )/(∆x m ) 2 for all x ∈ K i and t ∈ (t k-1 , t k ], then ζ m → ∆φ in L ∞ (Q T )
and therefore using Proposition 12 and Corollary 14 one obtains

F m 2 + Q T u 1 µ 1 ∆φ dxdt = Q T (u 1 µ 1 ∆φ -u 1,m µ 1,m ζ m ) dxdt → 0 as m → 0.
This concludes the proof of Theorem 3.

Numerical experiments

In this section, we perform several numerical experiments to illustrate the behavior of the scheme.

6.1. Implementation. The scheme was implemented in dimension d = 1 and d = 2 using Matlab. The code is available at https://gitlab.inria.fr/herda/nonlocal-skt. In order to optimize the computational cost, a number of matrices can be pre-assembled and stored using a sparse matrix structure. This is the case for the matrix of the Laplacian and those related to the convolution kernels. Moreover, the assembling can be performed efficiently using the discrete Fourier transform. At each time step the nonlinear system is solved using a Newton method. Convergence of the Newton method is reached when the ∞ norm of the residue divided by the norm of the first guess gets less than a given tolerance, which we took to be 10 -10 in our experiments. An adaptive time step procedure is implemented in case the Newton method fails to converge. After maximum number of steps (50 in the experiments), if the target error is not attained, ∆t is divided by 2. If there was refinement on a given time step, ∆t is multiplied by two for the next time step. In the experiments below the Newton method never failed to converge and the time step remained constant along all the simulations. 6.2. Test case 1: Convergence for various convolution kernels and initial data. In this first test case, we investigate the convergence of the scheme in the case for the following nonlocal cross-diffusion system

∂ t u 1 -∂ 2 xx ((ρ * u 2 )u 1 ) = 0, ∂ t u 2 -2∂ 2
xx ((ρ * u 1 )u 2 ) = 0. The convolution kernel is taken to be either the Dirac measure, which we denote by ρ 0 , either by an approximation of a Dirac (46)

ρ δ (x) = δ -1 χ [-δ/2,δ/2] (x),
where χ A indicator function of the set A, or the smooth kernel

ρ smooth (x) = cos(ν L x) + 1.
with ν L = 2π/L. We consider two initial data, either the indicator functions (47)

u 0 1 (x) = χ [ L 9 , L 3 ] (x) , u 0 2 (x) = χ [ L 3 , 3L 4 ] (x) , x ∈ R/LZ. or the smooth functions (48) u 0 1 (x) = cos (ν L x) + 1 , u 0 2 (x) = sin (ν L x) + 1 , x ∈ R/LZ.
The final time of simulation is taken to be T = 5 and the domain has length L = 25. We run the scheme for a sequence of decreasing space and time steps. More precisely the number of points is N k = 32 • 2 k-1 for k = 1, . . . , 6 and the corresponding time step ∆t k = ∆t 0 • 4 -(k-1) , with ∆t 0 = 5.

Observe that the refinement of the time step allows to witness experimental convergence in space up to second order accuracy if it is attained. As we do not know the analytical solution for this system, we take as reference solution the computed solution on the finest mesh (N = 1024). Then the error for the k-th mesh is taken to be the ∞ norm between the k-th solution and the reference solution projected on the k-th mesh. From these errors the experimental order is evaluated by linear regression (in log scale). In Table 1, we report the experimental order of convergence and the error between the N = 512 mesh and N = 1024 mesh for each kernel and initial data. 1. Estimated order of convergence in space and absolute error at final time for the mesh N = 512 in L ∞ norm for various convolution kernels and initial data. Reference solution is for N = 1024. 6.3. Test case 2: From nonlocal to local cross-diffusion. As a second test case, we investigate numerically the rate of convergence for different metrics of the so-called localization limit. Namely we study the rate of convergence of solutions of the nonlocal cross-diffusion system (1)-( 2) towards solutions of its corresponding local version ( 4)-( 5) as the convolution kernel tends to a Dirac measure. Indeed, if theoretically this localization limit has been proved in [START_REF] Dietert | Persisting entropy structure for nonlocal cross-diffusion systems[END_REF][START_REF] Moussa | From nonlocal to classical Shigesada-Kawasaki-Teramoto systems: triangular case with bounded coefficients[END_REF], the proofs rely on some compactness method and no explicit "error" bounds are available (see also for instance [START_REF] Jüngel | Nonlocal cross-diffusion systems for multi-species populations and networks[END_REF]). The establishment of such explicit estimates seems to be a complex task. In order to get a better understanding of this problem we aim to study this question thanks to our finite volume scheme.

More precisely, in this test case we consider the same system as in the first test case with ρ = ρ δ for various values of δ ∈ [0, L]. The domain has length L = 25, the final time is T = 1 and the mesh is such that N = 1024 and ∆t = 10 -2 . We evaluate the error at time T between the solution (u

(δ) 1 , u (δ)
2 ) computed for the kernel ρ δ and (u

(0) 1 , u (0)
2 ) computed for the local cross-diffusion system in Wasserstein-1 norm,

W 1 (u (δ) 1 , u (0) 1 ) + W 1 (u (δ) 2 , u (0)
2 ), and in L p norms u

(δ) 1 -u (0) 1 L p + u (δ) 2 -u (0) 2
L p , with p = 1 or p = ∞. For the computation of the Wasserstein-1 norm we recall that in dimension 1, if f and g are non-negative integrable functions on R with the same mass, one has

W 1 (f, g) = F -G L 1 , with F = f , G = g and F (-∞) = G(-∞).
In practice, f and g are piecewise contant functions, thus the previous norm can be computed exactly numerically. On Figure 1, we plot the error as a function of δ/L = 4W 1 (ρ δ , ρ 0 )/L for the two initial data (47) and (48). For the smooth initial data (48) supported on the whole domain (up to one point), we observe convergence with rate O(W 1 (ρ δ , ρ 0 ) 2 ) for all the norms. For the discontinuous initial data (48) supported on part of the domain, there is no experimental convergence in L ∞ norm, and O(W 1 (ρ δ , ρ 0 ) α ) convergence with α ≈ 0.38 in L 1 norm and α = 0.65 in Wasserstein-1 norm. 6.4. Test case 3: Turing instabilities in prey-predator systems with nonlocal crossdiffusion. In this last test case, we consider the following system with nonlocal cross diffusion and reaction modelling a population of preys with density u 1 and predators with density u 2 . The system reads

∂ t u 1 -d 1 ∆u 1 = R 1 (u 1 , u 2 ), ∂ t u 2 -∆((d 2 + d 21 ρ 2 * u 1 )u 2 ) = R 2 (u 1 , u 2 ).
The precise reaction terms will be specified below. On the one hand, preys are subject to linear diffusion with constant diffusivity coefficient d 1 . However, the predators diffuse depending on the presence or the absence of preys. More precisely, the convolution kernel ρ 2 is chosen such that it is close to 0 near the origin and large away form the origin (up to a given distance). This models the fact that predators need not seek for preys when they are available at their position, while they shall diffuse more rapidly if higher densities of preys are ahead. The reaction terms will be chosen following the phytoplankton-herbivore model of Segel and Levin [START_REF] Simon | Hypothesis for origin of planktonic patchiness[END_REF] and a variation of Mimura-Nishiura-Yamaguti [START_REF] Mimura | Some diffusive prey and predator systems and their bifurcation problems. Bifurcation theory and applications in scientific disciplines[END_REF]. In both cases, the particularities are an autocatalytic effect on the phytoplankton's (preys) growth rate and a density-dependent mortality of herbivore (predators).

In the case of linear diffusion, this model is famous for exhibiting diffusive instabilities [START_REF] Simon | Hypothesis for origin of planktonic patchiness[END_REF] around the homogenenous equilibrium. The corresponding Turing patterns have been invoked to justify the patchiness of phytoplankton's distribution in the oceans [START_REF] Simon | Hypothesis for origin of planktonic patchiness[END_REF]. In [START_REF] Simon | Hypothesis for origin of planktonic patchiness[END_REF] Segel and Levin mention that in these models the assumption of passive diffusion is made for simplicity only; more complicated movement patterns can also lead to diffusive instability. Here we propose a more complex description model of the behavior of predators thanks to non-local cross-diffusion.

In the following, we illustrate numerically the persistence and the modification of Turing patterns in the presence of nonlocal cross-diffusion.

6.4.1. One dimensional case: Segel-Levin reaction term. We consider the one-dimensional case with the following reaction terms

R 1 (u 1 , u 2 ) = au 1 + eu 2 1 -bu 1 u 2 , R 2 (u 1 , u 2 ) = -du 2 2 + cu 1 u 2 .
where the parameters are a = b = c = d = 1 and e = 1 3 . Concerning the diffusion we consider two cases. In the first case, both species are driven by linear diffusion with d 1 = 0.05 for preys and d 2 = 2 and without cross-diffusion d 21 = 0. In the second case the preys are driven by linear diffusion with d 1 = 0.05 and the predators by nonlocal cross-diffusion with d 21 = 1 and the kernel

ρ 2 (x) = C r [x 2 χ (-r,r) (x) + (x -2r) 2 χ [r,2r) (x) + (x + 2r) 2 χ (-2r,-r] (x)],
with C r a normalizing constant such that ρ 2 = 1. This kernel vanishes at 0, has support on [-2r, 2r] and is maximal at x = ±r. It is designed to model the hunting behavior of predators which will diffuse if most of the preys are away from their position, with a detection radius equal to r and a maximal distance of detection of 2r. In both cases the simulation is performed on a domain of length L = 25 with N = 500 cells. The radius is taken to be r = 10L/49. The final time of simulation T = 500 and the time step is ∆t = 0.1. The initial data is taken as a small perturbation of the homogeneous equilibrium

u 0 1 (x) = ad bc -de + εχ [L/3,L/9] (x) , u 0 2 (x) = ac bc -de ,
with ε = 10 -2 . With the chosen parameters, the homogeneous equilibrium is linearly unstable in both the linear diffusion and the nonlocal cross-diffusion cases. Numerically we observe the solution converges in time towards an heterogeneous equilibrium in both cases. On Figure 2, we plot the densities of preys and predators at final time. The difference between the patterns in the two cases is illustrated. 6.4.2. Two dimensional case: Mimura-Nishiura-Yamaguti reaction term. Now we consider the two-dimensional case with the following reaction terms

R 1 (u 1 , u 2 ) = au 1 + eu 2 1 -du 3 1 -bu 1 u 2 , R 2 (u 1 , u 2 ) = -f u 2 -gu 2 2 + cu 1 u 2 ,
ρ sym 2 (x, y) = Cχ (3/8,1/2) (x 2 + y 2 )
, with C a normalizing constant such that ρ sym 2 = 1. The third case is the same has the second case with linear diffusion for preys and nonlocal cross-diffusion for predators but the kernel is not symmetric and given by ρ nonsym 2 (x, y) = Cχ (3/8,1/2) (x 2 + y 2 )χ [0,∞) (x)χ [0,∞) (y), with C a normalizing constant such that ρ nonsym 2 = 1. In terms of modelling, it means that predators only sense preys that are north-east of their position (upper right quadrant). The final time of simulation T = 20 and the time step is ∆t = 0.01. In any cases the simulation is performed on a domain of horizontal length L x = 4 and vertical length L y = 3 with 133 × 100 cells. The initial data is taken as a small perturbation of the homogeneous equilibrium u 0 1 (x, y) = 5 + εχ [Lx/9,4Lx/9]×[7Ly/9,8Ly/9] (x, y) , u 0 2 (x, y) = 10, and ε = 10 -2 . Once again with the chosen parameters, the homogeneous equilibrium is linearly unstable in all cases and the solution converges in time towards an heterogeneous equilibrium. On Figure 3, we plot the colormap density of preys at final time. The difference between the patterns in the three cases is illustrated. In the last case the patterns are consistent with the breaking of symmetry in the kernel ρ nonsym Proof. We define w 0 = u 0 and for n ≥ 1

w n = u 0 + ∆t n k=1 a k u k Π n k=1 (1 -∆ta k ) -1 .
we define an integer N T and a time step ∆t = T /N T and we introduce the sequence (t k ) 0≤k≤N T with t k = k∆t. Now, as in Section 2.2, we first discretize the initial conditions (54) as u 0 j,i = 1 ∆x K i u 0 j (x) dx, ∀1 ≤ i ≤ , j = 1, 2. (56) Now for (u k-1 1 , u k-1 2 ) ∈ R 2 given, the scheme writes as ∆x u k 1,i -u k-1

1,i ∆t + F k 1,i+ 1 2 -F k 1,i-1 2 = 0, ∀1 ≤ i ≤ , (57) ∆x u k 2,n -u k-1 2,n ∆t + F k 2,n+ 1 2 -F k 2,n-1 2 = 0, ∀1 ≤ n ≤ , (58)
where the numerical fluxes are defined by Now let us prove that, for a given 0 < k ≤ N T , the solutions to scheme (54)-( 63) satisfy (a priori) the inequality

F 1,i+ 1 2 = d 1 u k 1,i -u k 1,i+1
H(u k 1 , u k 2 ) -H(u k-1 1 , u k-1 2 ) ∆t + 4 d 1 d 12 -1 i=0 u k 1,i -u k 1,i+1 ∆x + 4 d 2 d 21 -1 n=0 u k 2,n -u k 2,n+1 ∆x ≤ 0. ( 64
)
In this aim let us assume that u k j,i > 0 for all 1 ≤ i ≤ and j = 1, 2. We multiply equation (57) by log(u k 1,i )/d 12 and equation (58) by log(u k 2,n )/d 21 , we sum both equations over i, n ∈ {1, . . . , } respectively, we apply some discrete integration by parts, we use definitions (62) and (63) and we sum the two relations to obtain

Q 1 + Q 2 + Q 3 + Q 4 ≤ 0, with Q 1 = H(u k 1 , u k 2 ) -H(u k-1 1 , u k-1 2 ) ∆t , Q 2 = d 1 d 12 -1 i=0 (u k 1,i -u k 1,i+1 )(log(u k 1,i ) -log(u k 1,i+1 )) ∆x ≥ 4 d 1 d 12 -1 i=0 u k 1,i -u k 1,i+1 ∆x , Q 3 = d 2 d 21 -1 n=0 (u k 2,n -u k 2,n+1 )(log(u k 2,n ) -log(u k 2,n+1 )) ∆x ≥ 4 d 2 d 21 -1 n=0 u k 2,n -u k 2,n+1 ∆x , Q 4 = 1 ∆x -1 i=0 -1 n=0 ∆xG x i+ 1 2 , x n+ 1 2 u k 1,i -u k 1,i+1 u k 2,n -u k 2,n+1 A k i,n u k 1,i -u k 1,i+1 u k 2,n -u k 2,n+1
, where we have used for Q 1 the convexity of h 1 and h 2 , the inequality (a -b)(log(a) -log(b)) ≥ 4( √ a -√ b) 2 for Q 2 and Q 3 and where the matrix A k i,n is defined by log(u k j,i ) -log(u k j,i+1 ) if u k j,i > 0, u k j,i+1 > 0, and u k j,i = u k j,i+1 , u k j,i if u k j,i = u k j,i+1 > 0. Now, thanks to the nonnegativity of the kernel G, it remains to show that A k i,n is a positive semidefinite matrix in order to conclude that Q 4 ≥ 0. In this purpose, we notice that det(A k i,n ) = 0 and trace(A k i,n ) > 0 this imply that A k i,n has two nonnegative eigenvalues. Therefore we conclude that for all 0 ≤ i, n ≤ -1 the matrix A k i,n is a positive semidefinite matrix and Q 4 ≥ 0. This implies that the entropy inequality (64) holds which proves the claim.

A k i,n =        u k 2,n+ 1 2 u k,log

|∇ u 1

 1 (x)u 2 (x -y)| 2 dxdy = 0 .

3. 1 .

 1 Entropy dissipation and mass conservation. In the following , denotes the Euclidean scalar product and | • | the Euclidean norm.

Proposition

) of Theorem 1 . 11 .

 111 Proposition Let the assumptions of Theorem 1 hold. Then there exists a constant C 1 > 0 only depending on d 12 , d 21 , m 0 1 , m 0 2 and H

W 1 Figure 1 .

 11 Figure 1. Distance between solution of the nonlocal and local cross-diffusion system at final time versus δ/L. Left: Initial data is indicator function (47); Right: Initial data is the smooth function (48) .

Figure 2 .

 2 Figure 2. Turing patterns at final time for (left) linear diffusion for predators and preys (d 1 = 0.05, d 2 = 2 and d 21 = 0) and (right) cross-diffusion for predators and linear diffusion for preys (d 1 = 0.05, d 2 = 0 and d 21 = 1).

with a = 35 / 9 ,

 359 b = c = f = 1, d = 1/9, e = 16/9 and g = 2/5. For these reaction terms, the homogeneous equilibrium is given for the preys by[f + (c(b 2 c 2 -2bceg + 4df bg + e 2 g 2 + 4adg 2 ) 1/2 -bc 2 + ceg -2df g)/(2dg)]/c = 5,and by[c(b 2 c 2 -2bceg + 4df bg + e 2 g 2 + 4adg 2 ) 1/2 -bc 2 + ceg -2df g]/(2dg 2 ) =10, for predators. Concerning the diffusion we consider three cases. In the first case, both species are driven by linear diffusion with d 1 = 0.001 for preys and d 2 = 4 and without cross-diffusion d 21 = 0. Similarly to the one dimensional test case, in the second case the preys are driven by linear diffusion with d 1 = 0.001 and the predators by nonlocal cross-diffusion with d 21 = 2/5 and the kernel is the indicator function of an annulus

2 .

 2

Figure 3 .

 3 Figure 3. Turing patterns in prey density u 1 at final time: (top left) linear diffusion for predators and preys (d 1 = 0.001, d 2 = 4, d 21 = 0); (top right) nonlocal cross-diffusion for predators with symmetric kernel and linear diffusion for preys (d 1 = 0.001, d 2 = 0, d 21 = 2/5, ρ 2 = ρ sym 2 ); (bottom) nonlocal cross-diffusion for predators with non-symmetric kernel and linear diffusion for preys (d 1 = 0.001, d 2 = 0, d 21 = 2/5, ρ 2 = ρ nonsym 2

2 = 0 ,

 20 ∀1 ≤ k ≤ N T , j = 1, 2. (61)Finally, for every 1 ≤ i, n ≤ we define

  The diffusion coefficients d 1 , d 2 , d 11 and d 22 are non-negative constants and the crossdiffusion coefficients d 12 and d 21 are positive constants. (H3) The convolution kernels ρ 1 , ρ 2 , σ 1 and σ 2 are L ∞ (T) functions that are non-negative and satisfy the symmetry hypotheses

	2	, j = 1, 2.
	2.3. Main results. Let us collect our assumptions.
	(H1) The domain is taken as Ω = T.
	(H2)

  = (u k j,i + u k j,i+1)/2 and a similar expression for u k 2,n+1/2 and where we impose the following boundary conditions

	(59)			∆x	+ d 12 µ k 2,i+ 1 2	u k 1,i -u k 1,i+1 ∆x	+ d 12 u k 1,i+ 1 2	μk 2,i+ 1 2	, ∀1 ≤ i ≤ ,
	(60)	F 2,n+ 1 2	= d 2	u k 2,n -u k 2,n+1 ∆x	+ d 21 µ k 1,n+ 1 2	u k 2,n -u k 2,n+1 ∆x	+ d 21 u k 2,n+ 1 2	μk 1,n+ 1 2	, ∀1 ≤ n ≤ ,
	with u k 1,i+1/2

Acknowledgment. MH acknowledges support from the LabEx CEMPI (ANR-11-LABX-0007) and the ministries of Europe and Foreign Affairs (MEAE) and Higher Education, Research and Innovation (MESRI) through PHC Amadeus 46397PA. AZ acknowledges support from the multilateral project of the Austrian Agency for International Co-operation in Education and Research (OeAD), grant FR 01/2021.

Appendix A. Proof of a discrete Grönwall inequality 28 Appendix B. Discrete entropy inequality on a bounded domain 30 References Appendix A. Proof of a discrete Grönwall inequality For the sake of completness we prove in this section a discrete Grönwall inequality Lemma 15. Let (u n ) n∈N and (a n ) n∈N be some non-negative sequences with

Then, for every n ≥ 1, a simple computation leads to

Thanks to (49) we deduce that the sequence (w n ) n∈N is decreasing and since w 0 = u 0 we end up with

which finishes the proof of Lemma 15.

Appendix B. Discrete entropy inequality on a bounded domain

Following [START_REF] Dietert | Persisting entropy structure for nonlocal cross-diffusion systems[END_REF], let us assume that Ω = (0, 1) and d 11 = d 22 = 0. then we consider, for all x, y ∈ Ω and t ≥ 0, the system

where G is a C 0 nonnegative kernel such that

We supplement the system (50)-( 52) with homogeneous Neumann boundary conditions

and initial conditions

In this case we define the entropy functional H as

where we recall definition [START_REF] Chen | Analysis of a parabolic cross-diffusion population model without self-diffusion[END_REF] of h 1 and h 2 . Then, as proved in [START_REF] Dietert | Persisting entropy structure for nonlocal cross-diffusion systems[END_REF]Proposition 7], the following entropy inequality holds d dt

In this section we design a finite volume scheme for system (50)-( 54) which preserves at the discrete level (at least formally) the entropy inequality (55).

In this purpose we introduce some notations related to the discretization of Ω = (0, 1). An uniform mesh, consists in a finite sequence of cells denoted by K i = x i-1/2 , x i+1/2 , for 1 ≤ i ≤ , where the points x i+1/2 are uniformly distributed such that

We note ∆x = x i+ 1 2 -x i-1 2 , for some 1 ≤ i ≤ , the length of the mesh. Moreover, for 1 ≤ i ≤ , we define x i as the center of the cell K i and x 0 = x 1 2 and x +1 = x + 1 2 . Then, for T > 0 given,