Study of an entropy dissipating finite volume scheme for a nonlocal cross-diffusion system

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 on a periodic domain Ω = T  ( ≤ 3).For a given final time  we denote the space time domain by   = Ω × (0,  ).The parameters  1 ,  2 ,  11 ,  12 ,  21 and  22 are some positive constants and  1 ,  2 ,  1 and  2 are non-negative convolution kernels.System (1) and ( 2) is supplemented with initial conditions 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 [42] which can describe segregation phenomena between competing species.It writes Nonlocal cross-diffusion systems appear naturally as a mean field type of limit of interacting many-particle systems.For instance, the model ( 1) and (2) was introduced in [26] 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 cross-diffusion models.This interpretation has been mathematically justified in the literature, see [19,21,32,37], 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 [28,40].In the model ( 1) and ( 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 [38,39], tumour growth [22], opinion formation [23] or interactions between spiking neurons [5] (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)- (3).Our approach is inspired by the analysis performed at the continuous level in [21,37].In particular, in [21] 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 [29,30] in the local case.Indeed, it was shown that for the system (1)-(3) without reaction terms  1 =  2 = 0 and under the following symmetry hypotheses on the convolution kernels ⎧ ⎪ ⎨ ⎪ ⎩  1 () =  2 (−) = (),  1 () =  1 (−),  2 () =  2 (−), (6) that the following entropy functional [ 1 (log( which was already known for the local SKT system (see [17,18,27]).The fact that (1) and ( 2) admits a Lyapunov functional is crucial for the study of the system.Indeed, in [21], the authors used (7) together with the so-called duality method, see [20,34,37], in order to prove (assuming (6) and without reaction terms) the existence of distributional solutions to (1)-(3).Definition 1.Given  > 0, let  1 ,  2 ,  1 and  2 be some functions in  ∞ (Ω) and  0 1 and  0 2 be some initial functions in  1 (Ω).Then, we say that the measurable functions  1 ,  2 :   → R + are distributional solutions to (1)-( 3) if for every  ∈  ∞ 0 (Ω × [0,  )) it holds and In this paper we propose and analyze a finite volume scheme for (1)- (3).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 [4,7,10,11,31,43], the preservation of the entropy dissipation property at the discrete level is crucial.This ensures well-posedness and global stability in time [15,25] as well as with respect to the choice of convolution kernels (see Thm. 2).Some of these methods are reminiscent of the second author's work in [31] 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 Thm. 3) by adapting the duality method (see [20,34,37]) 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 4. Let us mention that only  ∞ regularity of the convolution kernels is required to obtain convergence of the scheme.With smoother kernels one additionally obtains local in time  ∞ bounds on the discrete solution that are uniform in the mesh size (see Thm. 3).
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 [2,3] 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 [33] (see also [41]) obtained in [4].We mention [9] 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 [12] and then extended for the multi-species case in [13].This scheme is a positivity and entropy preserving method as shown in [13].Finally, we also refer to [14].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 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 Ω = ∏︀  =1 R/  Z, where  1 , . . .,   > 0. However, for the sake of readability we will assume from now on that  = 1 and  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.

Notations and definitions
Let us define  ≥ 1 and ∆ = 1/ .A uniform mesh  of T consists in a finite sequence of cells denoted by centered at   = ∆ and with extremities  ± 1 2 = ( ± 1 2 )∆.For  > 0 given, we define an integer   and a time step ∆ =  /  and we introduce the sequence (  ) 0≤≤  with   = ∆.We denote by  a space-time discretization of   = T × (0,  ) composed of a space discretization  of T and the values (∆,   ).
Let us now introduce some discrete norms on the space of piecewise constant functions in space }︃ .
For  ∈ [1, ∞), we define the discrete  1, seminorm and discrete  1, norm on ℋ  by where, for the  ∞ (T) norm given by Let us also recall the definition of the space  (T), see [1] for more details.A function  ∈  1 (T) belongs to the space  (T) if its total variation   () given by is finite.We endow the space  (T) with the norm In particular, we notice that for each function  ∈ ℋ  ∩  (T) we have ‖‖  (T) = ‖‖ 1,1, .Finally we introduce the space ℋ  of piecewise constant in time functions with values in ℋ  , }︃ .

2
) ∈ R 2 , the implicit in time numerical scheme writes as where ∆  denotes the discrete Laplacian, namely and with 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) and (5).We also remark that we could equivalently rewrite (12) as where for all  ∈ ℐ the numerical fluxes ℱ  ,+ with the centered approximation at interfaces ,  = 1, 2.

Main results
Let us collect our assumptions.
(H1) The domain is taken as Ω = T.
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 -dimensional periodic domain Ω = ∏︀  =1 R/  Z. 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  12 = 0 or  21 = 0 (see Sect. 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)- (15).However, in terms of practical use, the scheme performs well even when dropping this hypothesis (see Sect. 6.4).Following for instance [31], the assumption (H5) can be relaxed and one can extend the proofs of Theorems 2 and 4 in the case of the Lotka-Volterra source terms: with  0 and   some nonnegative constants for ,  = 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 where the functions ℎ 1 and ℎ 2 are defined by ( (log() − 1) + 1), ∀ ∈ (0, +∞), (19) with the obvious continuous extension at  = 0.The corresponding entropy dissipation functional is defined by Theorem 2 (Existence of solutions).Let the assumptions (H1)-(H5) hold.Then, for every 1 ≤  ≤   there exists (at least) one nonnegative solution (  1 ,   2 ) to scheme (11)- (15).Moreover, this solution satisfies the following properties: (ii) Entropy production estimate: for all  ≥ 1 it holds Finally for all  ∈ {1, 2}, if   is positive then    is positive for all 0 <  ≤   .The proof of existence of Theorem 2 is based on a consequence (see [24], Sect.9.1) of the Brouwer fixed point Theorem.It can be applied thanks to the a priori entropy-dissipation estimate (22) and regularization inspired by Burger et al. [8] and Jüngel [30].It also follows the line of the existence proof of [16].
The second main result is concerned by some properties satisfied by the solutions of the scheme ( 11)-( 15).These estimates are discrete counterparts of Theorem 9 from [21].
Theorem 3 (Qualitative properties of the solutions).Let the assumptions of Theorem 2 hold.Moreover, assume that  0 1 ,  0 2 ∈  2 (T) and let  and Γ be some nonnegative constants such that

Finally let us introduce 𝑚
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 then for all  ∈ ℐ,  ≥ 1 and  = 1, 2 we have where (ii) Duality estimate: If  1 and  2 are positives, then there exists a constant  > 0 which is independent of the mesh size such that )︁ , where The proof of Theorem 3 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 emphasize 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 (  ) ∈N of space-time discretizations of   indexed by the size   = max{∆  , ∆  } of the mesh, satisfying   → 0 as  → ∞.We denote by   the corresponding mesh of T and by ∆  the corresponding time step.Finally, for every  ∈ N we set ( 1, ,  2, ) ∈ ℋ  the picewise constant in space and time reconstruction of the solutions to the scheme ( 11)-( 15) corresponding to the mesh   .
Theorem 4 (Convergence of the scheme).Let the assumptions of Theorem 2 hold, assume that the coefficients  1 and  2 are positives and let (  ) ∈N be a family of space-time discretizations of   with   → 0 as  → ∞.Then, if we denote by ( 1, ,  2, ) a family of finite volume solutions to (11)-( 15) obtained in Theorem 2, there exists ( 1 ,  2 ) ∈ (  (  )) 2 for  ∈ [1, 3) a distributional solutions to (1)-( 3) in the sense of Definition 1 such that, up to a subsequence, for  = 1, 2 it holds The proof of Theorem 4 is based on uniform estimates w.r.t.∆ and ∆, established in Section 5.1.These estimates allow us to apply in Section 5.2 a compactness result obtained in [4] which yields, up to a subsequence, the strong convergence in   (  ) of the sequence ( 1, ,  2, ) towards the functions  1 and  2 stated in Theorem 4.Then, we identify in Section 5.3 the functions  1 and  2 as distributional solutions in the sense of Definition 1 of the nonlocal cross-diffusion system (1)-(3).
Remark 5. Let us make few remarks concerning Theorem 4.
-If  > 1, the convergence of the scheme can also be established.However in this case we obtain, up to a subsequence, for  = 1, 2, see Remark 15 for more details.-If  1 =  2 = 0, then it is still possible to conclude if the convolution kernels are smooth enough ( 2 for instance).Indeed in this case from weak compactness on ( , )  , strong compactness can be obtained on its convolution with the smooth kernel.Moreover, if the convolution kernels are  2 , then one can prove a stability estimate (in  2 -norm on the difference between two solutions) for the solutions to (1)-( 3) which provides uniqueness and continuous dependence on the initial data at the continuous level.As a by-product we deduce that in this case the whole sequence ( 1, ,  2, ) converges as  → ∞ instead of only a subsequence.In the discrete setting, a counterpart of the stability estimate in  2 -norm can also be established uniformly in ∆ at the price of additional  1 regularity on the initial data.This difference with the continuous setting comes from the fact that a Grönwall argument with implicit schemes requires a condition on the time step (see Prop. 10 for an illustration of this fact).Under these assumptions, one also gets uniqueness of solutions to the scheme.Uniqueness may also be obtained without additional regularity assumptions provided that a CFL condition holds (see Rem. 7 for details).-Finally, with enough regularity on the data, one could derive quantitative error estimates between approximate and continuous solutions.

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 ( 1 ,  2 ) ∈ ((0, +∞)  ) 2 we define the new unknown  = Φ( 1 ,  2 ) where Φ : From there finding a positive solution to the scheme ( 11)-( 15) amounts to finding a zero   = Φ(  1 ,   2 ) of the continuous map   : R 2 → R 2 defined for any ( 1 ,  2 ) ∈ ((0, +∞)  ) 2 by its components where ) are given and  1 ,  2 are related to  1 ,  2 through the relation ( 14) and ( 15) dropping the index .) be componentwise non-negative.Then for any  ∈ R 2 ,
Proof.In order to prove (23), it suffices to sum the components of   () and observe that since it is a telescopic sum.Concerning the inequality, first observe that ( −1 ) is non-negative.Then, using the definition of Φ( 1 ,  2 ) and   one obtains with Using the convexity of  ↦ → ( log() −  + 1) to bound both  1 and  2 from below, one obtains Then for  1 , a discrete integration by parts (or summation by parts) yields and a similar formula holds for  2 .Using the definitions of  1 and  2 (see ( 14) and ( 15) without the exponents), one has and a similar estimate for  diff 2 .For the second term one has In the previous estimate, the second inequality is obtained by changing  into − and using the symmetry of  1 .
For the third equality, one changes  into −.The fourth one is the combination of the first and third equalities.
Once again a similar estimate holds for  2 2 .Finally with the same changes of indices one can estimate the sum By summing all the estimates one obtains (24).The last point of the proposition is obtained by induction.
Let us finally prove that if   > 0,   , > 0 for all  ∈ ℐ and 0 <  ≤   .This is a consequence of the entropy estimate.For a given 0 <  ≤   we notice that (thanks to the term  diff  ) the positive solution (  1 ,   2 ) satisfies the following estimate At the limit  → 0, let us assume by contradiction that there exists  ∈ ℐ such that   , = 0. Then as the r.h.s. of the previous inequality is finite this implies that   ,+1 = 0. Thus repeating this argument we deduce that   , = 0 for all  ∈ ℐ. Consequently we have ‖   ‖  1 (T) = 0 which contradicts the mass conservation property (21).This completes the proof of Theorem 2.
Remark 7 (Uniqueness under CFL).Under a parabolic CFL condition, uniqueness of a solution to the scheme can also be proven.Indeed, if one denotes by   a vector of solution given by Theorem 2, then the scheme may be rewritten as M (  )  =  −1 .The matrix M (which has a similar definition as (26) hereafter) depends on   through   1 and   2 and is a perturbation of the identity matrix of size ‖  ∞ is bounded uniformly in terms of the initial masses,  ∞ norms of the convolution kernels and diffusion coefficients, Δ Δ 2 can be made small enough with respect to these quantities only so that a contraction argument yields the uniqueness of solutions.

Estimates on the discrete Kolmogorov equation
In this section, we focus on estimates concerning the finite volume discretization of the Kolmogorov equation    = ∆().In particular we adapt at the discrete level some properties established in [21,37].
In the rest of this section, we assume that (   ) ∈ℐ ,  = 1, . . .,   is given and componentwise non-negative.From there, the scheme is given for all  ≥ 1 by where ∆  denotes the discrete Laplacian operator defined by (13).

Well-posedness of the scheme and 𝐿 ∞ estimates
Let us first prove that the scheme (25) admits a unique solution at each time step.
Proof.Let us write  −1 = ( −1 0 , . . .,  −1  −1 ) ⊤ for all  ≥ 1. Observe that the scheme writes M    =  −1 where M  is a  ×  tridiagonal matrix defined by We notice that M  has positive diagonal terms and non-positive off-diagonal terms.Furthermore the matrix M  is strictly diagonally dominant with respect to its columns.Therefore M  is a non-singular M-matrix and is thus monotone and invertible.This finishes the proof of Lemma 8.
We prove in the following result some  ∞ estimates for the solution to scheme (25).
Proof.We will only deal with the upper bound in (27) and the lower bound is obtained in the same way.Let M  denote the tridiagonal matrix defined by (26) and define We proceed by induction.Since Γ0 = Γ the bound holds by hypothesis at  = 0. Then observe that for every  ∈ ℐ Now we notice that by construction .
Then we easily deduce that for every  ∈ ℐ it holds Therefore, since M  is a M-matrix we conclude that    ≤ Γ for all  ∈ ℐ which concludes the proof of Lemma 9.
The bounds of Lemma 9 are exactly the discrete equivalent of the  ∞ estimates established at the continuous level in Corollary 18 of [21].
Proposition 10.Let us assume that it holds .
Then the solution to (25) satisfies the following estimate Proof.Let  ≥ 1 be fixed and let us first notice that we can rewrite for every  ∈ ℐ equation (25) as where Now we multiply the above equation by ∆   and we sum over  ∈ ℐ, we obtain where For  3 using the inequality ( − ) ≥ ( 2 −  2 )/2 we obtain For  4 applying a discrete integration by parts yields Now we rewrite  5 as and reordering the terms in the r.h.s. the second sum vanishes and we have Gathering ( 28)-( 30) we end up with 1 2 We deduce that 1 2 One ends the proof of Proposition 10 thanks to a discrete Grönwall inequality.

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 (25), see for instance ( [37], Thm. 3).In this aim, following [37], we introduce a "dual" scheme associated to (25).Let    +1  be given for every  ∈ ℐ, then for 1 ≤  ≤   we want to determine the solution to the following implicit backward in time scheme where    is given and non-negative and   = (  0 , . . .,    −1 ) is some given vector in R  for all 1 ≤  ≤   .Let us notice that (31) define a set of linear equation which can be rewritten as where M  is the tridiagonal matrix given by (26).Therefore, it follows directly from the proof of Lemma 8 that the problem (32) admits a unique solution for every 1 ≤  ≤   .
Prior to the proof of the discrete duality estimate, see Theorem 12 below, we establish some uniform estimates satisfied by the solution of (32).Proposition 11.Assume that min ∈ℐ    > 0 for every 0 ≤  ≤   and that    +1  = 0 for every  ∈ ℐ.Then the solution to (32) satisfies for every 1 ≤  ≤   the following estimate and there exists a constant  > 0 independent of ∆ such that where  and  denote the piecewise reconstruction functions in ℋ  associated to the vectors (  ) 1≤≤  and (  ) 1≤≤  .
We now prove estimate (34).In this purpose we multiply (31) by ∆∆, we sum over  ∈ ℐ and  ∈ {, . . .,   } and we obtain Applying the Cauchy-Schwarz inequality leads to Using estimate (33) we obtain Now it remains to apply the discrete Poincaré-Wirtinger inequality on the torus obtained in Lemma 6 of [6] in order to conclude the proof of Proposition 11.
We are now in position to establish the discrete dual estimate.

Proof of Theorem 3
We are now able to prove Theorem 3.
Step 1: Maximum principle.Let us first prove the maximum principle satisfies by the solutions to ( 11)-( 15).Let us notice that for every 1 ≤  ≤   we have max Now, let us recall that  0  = ‖ 0  ‖  1 (T) for  = 1, 2, then thanks to the mass conservation property (21) we obtain max Similarly we establish the following bound As a direct consequence of the previous estimates and ( 27) (with γ =  and Γ = Γ) one obtains point (i) of Theorem 3.
Thanks to Theorem 2, we have   1, ,   2, > 0 for all  ∈ ℐ and the element    is well-defined.Besides, applying the discrete duality estimate established in Theorem 12 we deduce the existence of a constant  > 0 independent of ∆ such that Now we notice that For  9 we have Thus, bearing in mind the mass conservation property ( 21) we obtain and similarly Collecting ( 38)- (40) we conclude that point (ii) of Theorem 3 holds.

Convergence of the scheme
This section is dedicated to the proof of Theorem 4. In the following the subscript  refer to the size   = max{∆  , ∆  } of the family (  ) of space-time discretizations of   .We derive uniform in  a priori estimates in Section 5.1 in order to obtain compactness in   (  ) of the sequences of constant by part reconstructions ( , )  for both species  = 1, 2. The compactness results are gathered in Section 5.2.A keypoint is a discrete  1 compactness result obtained in Lemma 9.2 of [4].This result is the adaptation at the discrete level of a compactness lemma established by Kruzhkov in [33] (see also [41]).Finally in Section 5.3, we prove Theorem 4.

Uniform estimates
In this section we establish some uniform estimates w.r.t.∆ and ∆ fulfilled by the solutions to the scheme (11)- (15).They rely on the entropy dissipation inequality (22) and the conservation of mass (21) of Theorem 2.
Therefore, applying the entropy inequality (22), we get for the first species and the equivalent estimate holds for the second species.This yields the existence of  1 such that (41) holds.It remains to establish (42).In this purpose we will consider the case  = 1.Then, using the definition (18) of the numerical fluxes, we estimate For  11 , applying the regularity of the functions  1 and  we have Hence, using the conservativity of the scheme and (41), we get For  12 , using the definition of   1, for  ∈ ℐ, we notice that it holds Then, thanks to the conservativity of the scheme, we obtain Therefore, applying the Cauchy-Schwarz inequality and (41) we end up with Collecting ( 43) and ( 44) and the corresponding inequalities for the second species lead to the existence of  2 such that (42) holds.This concludes the proof of Proposition 13.
Remark 15.In dimension  ≥ 2, we have the compact embedding of the space  (T  ) in   −1 (T  ).In particular in this case the sequence ( , ) is uniformly bounded in   −1 (  ).Therefore arguing as in the previous proof we deduce the existence for  = 1 and 2 of   ∈   (  ) for  ∈ [1, /( − 1)), such that, up to a subsequence, Corollary 16.Let the assumptions of Proposition 14 hold.Then there exists a subsequence of ( where we recall that  1 () =  2 (−) = () for a.e. ∈ T.
The first factor in the right hand side tends to 0 (by density of continuous functions in  /(−1) (T)) while, bearing in mind Proposition 14, the second factor is uniformly bounded in   (  ).Therefore one can conclude the strong convergence in   (0,  ;  ∞ (T)) by using Young's inequality and the previous argument.This finishes the proof of Corollary 16.

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

Implementation
The scheme was implemented in dimension  = 1 and  = 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, ∆ is divided by 2. If there was refinement on a given time step, ∆ 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.

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 The convolution kernel is taken to be either the Dirac measure, which we denote by  0 , either by an approximation of a Dirac where   indicator function of the set , or the smooth kernel  smooth () = cos(  ) + 1.
with   = 2/.We consider two initial data, either the indicator functions or the smooth functions The final time of simulation is taken to be  = 5 and the domain has length  = 25.We run the scheme for a sequence of decreasing space and time steps.More precisely the number of points is   = 32 • 2 −1 for  = 1, . . ., 6 and the corresponding time step ∆  = ∆ 0 • 4 −(−1) , with ∆ 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 ( = 1024).Then the error for the -th mesh is taken to be the ℓ ∞ norm between the -th solution and the reference solution projected on the -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  = 512 mesh and  = 1024 mesh for each kernel and initial data.

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) and (2) towards solutions of its corresponding local version (4) and ( 5) as the convolution kernel tends to a Table 1.Estimated order of convergence in space and absolute error at final time for the mesh  = 512 in  ∞ norm for various convolution kernels and initial data.Reference solution is for  = 1024.
1. Distance between solution of the nonlocal and local cross-diffusion system at final time versus /.Left: initial data is indicator function (47); Right: initial data is the smooth function (48).and a variation of Mimura-Nishiura-Yamaguti [36].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 [35] around the homogenenous equilibrium.The corresponding Turing patterns have been invoked to justify the patchiness of phytoplankton's distribution in the oceans [35].In [35] 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 crossdiffusion.In the following, we illustrate numerically the persistence and the modification of Turing patterns in the presence of nonlocal cross-diffusion.

One dimensional case: Segel-Levin reaction term
We consider the one-dimensional case with the following reaction terms where the parameters are  =  =  =  = 1 and  = 1 3 .Concerning the diffusion we consider two cases.In the first case, both species are driven by linear diffusion with  1 = 0.05 for preys and  2 = 2 and without cross-diffusion  21 = 0.In the second case the preys are driven by linear diffusion with  1 = 0.  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.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 2 .