# Study of a structure preserving finite volume scheme for a nonlocal cross-
diffusion system

Maxime Herda  Inria, Univ. Lille, CNRS, UMR 8524–Laboratoire Paul Painlevé,
59000 Lille<EMAIL_ADDRESS>and  Antoine Zurek  Université de
Technologie de Compiègne, LMAC, 60200 Compiègne, France<EMAIL_ADDRESS>
###### Abstract.

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.

Keywords: Nonlocal cross-diffusion, finite volume schemes, entropy method,
convergence.

Mathematics Subject Classification: 65M08, 65M12, 35K51, 35Q92, 92D25.

###### Contents

  1. 1 Introduction
  2. 2 Numerical scheme and main results
    1. 2.1 Notations and definitions
    2. 2.2 Numerical scheme
    3. 2.3 Main results
  3. 3 Existence of solution and entropy dissipation estimate
    1. 3.1 Entropy dissipation and mass conservation
    2. 3.2 Proof of Theorem 1
  4. 4 Estimates on the discrete Kolmogorov equation
    1. 4.1 Well-posedness of the scheme and $L^{\infty}$ estimates
    2. 4.2 Study of the dual problem
    3. 4.3 Proof of Theorem 2
  5. 5 Convergence of the scheme
    1. 5.1 Uniform estimates
    2. 5.2 Compactness properties
    3. 5.3 Proof of Theorem 3
  6. 6 Numerical experiments
    1. 6.1 Implementation
    2. 6.2 Test case 1: Convergence for various convolution kernels and initial data
    3. 6.3 Test case 2: From nonlocal to local cross-diffusion
    4. 6.4 Test case 3: Turing instabilities in prey-predator systems with nonlocal cross-diffusion
      1. 6.4.1 One dimensional case: Segel-Levin reaction term
      2. 6.4.2 Two dimensional case: Mimura-Nishiura-Yamaguti reaction term
  7. A Proof of a discrete Grönwall inequality
  8. B Discrete entropy inequality on a bounded domain

##  1\. Introduction

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

(1) $\displaystyle\partial_{t}u_{1}-\Delta((d_{1}+d_{11}\,\sigma_{1}\ast
u_{1}+d_{12}\,\rho_{1}\ast u_{2})u_{1})$ $\displaystyle=R_{1}(u_{1},u_{2}),$
(2) $\displaystyle\partial_{t}u_{2}-\Delta((d_{2}+d_{21}\,\rho_{2}\ast
u_{1}+d_{22}\sigma_{2}\ast u_{2})u_{2})$ $\displaystyle=R_{2}(u_{1},u_{2}),$

on a periodic domain $\Omega={\mathbb{T}}^{d}$ ($d\leq 3$). For a given final
time $T$ we denote the space time domain by $Q_{T}=\Omega\times(0,T)$. The
parameters $d_{1}$, $d_{2}$, $d_{11}$, $d_{12}$, $d_{21}$ and $d_{22}$ are
some positive constants and $\rho_{1}$, $\rho_{2}$, $\sigma_{1}$ and
$\sigma_{2}$ are non-negative convolution kernels. System (1)–(2) is
supplemented with initial conditions

(3) $\displaystyle u_{1}(\cdot,0)=u_{1}^{0}(\cdot),\quad
u_{2}(\cdot,0)=u_{2}^{0}(\cdot).$

In the case where the convolution kernels are given by the Dirac measure
$\rho_{1}=\rho_{2}=\sigma_{1}=\sigma_{2}=\delta_{0}$, the system coincides
with the celebrated Shigesada, Kawasaki, and Teramoto (SKT) population model
[41] which can describe segregation phenomena between competing species. It
writes

(4)
$\displaystyle\partial_{t}u_{1}-\Delta((d_{1}+d_{11}\,u_{1}+d_{12}\,u_{2})u_{1})$
$\displaystyle=R_{1}(u_{1},u_{2}),$ (5)
$\displaystyle\partial_{t}u_{2}-\Delta((d_{2}+d_{21}\,u_{1}+d_{22}\,u_{2})u_{2})$
$\displaystyle=R_{2}(u_{1},u_{2}).$

Nonlocal cross-diffusion systems appear naturally as a mean field type of
limit of interacting many-particle systems. For instance, the model (1)–(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 [17,
21, 31, 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)–(2) assume for instance
that $\rho_{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
$R_{1}=R_{2}=0$ and under the following symmetry hypotheses on the convolution
kernels

(6) $\left\\{\begin{aligned} &\rho_{1}(x)=\rho_{2}(-x)=\rho(x),\\\
&\sigma_{1}(x)=\sigma_{1}(-x),\\\
&\sigma_{2}(x)=\sigma_{2}(-x),\end{aligned}\right.$

that the following entropy functional

$\displaystyle
H(u_{1},u_{2})=\int_{\Omega}\frac{1}{d_{12}}\left[u_{1}(\log(u_{1})-1)+1\right]\,{\mathrm{d}}x+\int_{\Omega}\frac{1}{d_{21}}\left[u_{2}(\log(u_{2})-1)+1\right]\,{\mathrm{d}}x\,,$

is dissipated along solutions of (1)–(2). More precisely one has

(7)
$\frac{\mathrm{d}}{\mathrm{d}t}H(u_{1},u_{2})+\frac{2d_{11}}{d_{12}}\int_{\Omega}\int_{\Omega}\sigma_{1}(y)\left(\sqrt{u_{1}(x-y)}\,\nabla\sqrt{u_{1}(x)}+\sqrt{u_{1}(x)}\nabla\sqrt{u_{1}(x-y)}\right)^{2}\,\mathrm{d}x\mathrm{d}y\\\
+\frac{2d_{22}}{d_{21}}\int_{\Omega}\int_{\Omega}\sigma_{2}(y)\left(\sqrt{u_{2}(x-y)}\,\nabla\sqrt{u_{2}(x)}+\sqrt{u_{2}(x)}\nabla\sqrt{u_{2}(x-y)}\right)^{2}\,\mathrm{d}x\mathrm{d}y\\\
+4\frac{d_{1}}{d_{12}}\int_{\Omega}|\nabla\sqrt{u_{1}}|^{2}\,\mathrm{d}x+4\frac{d_{2}}{d_{21}}\int_{\Omega}|\nabla\sqrt{u_{2}}|^{2}+4\int_{\Omega}\rho(y)\int_{\Omega}|\nabla\sqrt{u_{1}(x)u_{2}(x-y)}|^{2}\,\mathrm{d}x\mathrm{d}y=0\,.$

Observe that if the convolution kernels are given by the Dirac measure
$\rho_{1}=\rho_{2}=\sigma_{1}=\sigma_{2}=\delta_{0}$, then

$\frac{\mathrm{d}}{\mathrm{d}t}H(u_{1},u_{2})+2\frac{d_{11}}{d_{12}}\int_{\Omega}|\nabla
u_{1}|^{2}\,\mathrm{d}x+2\frac{d_{22}}{d_{21}}\int_{\Omega}|\nabla
u_{2}|^{2}\,\mathrm{d}x\\\
+4\frac{d_{1}}{d_{12}}\int_{\Omega}|\nabla\sqrt{u_{1}}|^{2}\,\mathrm{d}x+4\frac{d_{2}}{d_{21}}\int_{\Omega}|\nabla\sqrt{u_{2}}|^{2}+4\int_{\Omega}|\nabla\sqrt{u_{1}\,u_{2}}|^{2}\,\mathrm{d}x=0\,,$

which was already known for the local SKT system (see [18, 19, 27]).

The fact that (1)–(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 $T>0$, let $\rho_{1}$, $\rho_{2}$, $\sigma_{1}$ and $\sigma_{2}$ be some
functions in $L^{\infty}(\Omega)$ and $u^{0}_{1}$ and $u^{0}_{2}$ be some
initial functions in $L^{1}(\Omega)$. Then, we say that the measurable
functions $u_{1},u_{2}:Q_{T}\to{\mathbb{R}}_{+}$ are distributional solutions
to (1)–(3) if for every $\phi\in C^{\infty}_{0}(\Omega\times[0,T))$ it holds

(8)
$\displaystyle\int_{Q_{T}}\Big{(}u_{1}\partial_{t}\phi+\left(d_{1}u_{1}+d_{11}\sigma_{1}\ast
u_{1}u_{1}+d_{12}\rho_{1}\ast
u_{2}u_{1}\right)\Delta\phi\Big{)}\mathrm{d}x\mathrm{d}t=-\int_{{\mathbb{T}}^{d}}u^{0}_{1}(x)\phi(x,0)\mathrm{d}x,$

and

(9)
$\displaystyle\int_{Q_{T}}\Big{(}u_{2}\partial_{t}\phi+\left(d_{2}u_{2}+d_{21}\rho_{2}\ast
u_{1}u_{2}+d_{22}\sigma_{2}\ast
u_{2}u_{2}\right)\Delta\phi\Big{)}\mathrm{d}x\mathrm{d}t=-\int_{{\mathbb{T}}^{d}}u^{0}_{2}(x)\phi(x,0)\mathrm{d}x.$

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, 32, 42], 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 Theorem 1).
Some of these methods are reminiscent of the second author’s work in [32]
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 [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 3. Let us mention that only $L^{\infty}$ regularity of the convolution
kernels is required to obtain convergence of the scheme. With smoother kernels
one additionally obtains local in time $L^{\infty}$ 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 [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] obtained in [4]. We mention [8] 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 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
$\Omega\subset{\mathbb{R}}^{d}$ the nonlocal entropy preserving version of the
cross-diffusion model differs from (4)–(5) because of boundary effects (see
[21, 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.

##  2\. Numerical scheme and main results

The results of this paper apply to a periodic domain
$\Omega=\prod_{i=1}^{d}\mathbb{R}/L_{i}\mathbb{Z}$, where
$L_{1},\dots,L_{d}>0$. However, for the sake of readability we will assume
from now on that $d=1$ and $L_{1}=1$, namely
$\Omega={\mathbb{T}}=\mathbb{R}/\mathbb{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\geq 1$ and $\Delta x=1/N$. A uniform mesh $\mathcal{T}$ of
${\mathbb{T}}$ consists in a finite sequence of cells denoted by

$K_{i}=(x_{i-\frac{1}{2}},x_{i+\frac{1}{2}})\,,\quad
i\in{\mathcal{I}}={\mathbb{Z}}/N{\mathbb{Z}}.$

centered at $x_{i}=i\Delta x$ and with extremities
$x_{i\pm\frac{1}{2}}=(i\pm\frac{1}{2})\Delta x$. For $T>0$ given, we define an
integer $N_{T}$ and a time step $\Delta t=T/N_{T}$ and we introduce the
sequence $(t_{k})_{0\leq k\leq N_{T}}$ with $t_{k}=k\Delta t$. We denote by
${\mathcal{D}}$ a space-time discretization of $Q_{T}={\mathbb{T}}\times(0,T)$
composed of a space discretization ${\mathcal{T}}$ of ${\mathbb{T}}$ and the
values $(\Delta t,N_{T})$.

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

$\displaystyle\mathcal{H}_{{\mathcal{T}}}=\Bigg{\\{}w:{\mathbb{T}}\to{\mathbb{R}}\,\,:\,\,w(x)=\sum_{i\in{\mathcal{I}}}w_{i}\mathbf{1}_{K_{i}}(x)\Bigg{\\}}.$

For $p\in[1,\infty)$, we define the discrete $W^{1,p}$ seminorm and discrete
$W^{1,p}$ norm on $\mathcal{H}_{\mathcal{T}}$ by

$\displaystyle|w|^{2}_{1,p,{\mathcal{T}}}=\left(\sum_{i\in{\mathcal{I}}}\Delta
x\left|\frac{w_{i+1}-w_{i}}{\Delta
x}\right|^{p}\right)^{\frac{1}{p}},\quad\|w\|_{1,p,{\mathcal{T}}}=|w|_{1,p,{\mathcal{T}}}+\|w\|_{L^{p}({\mathbb{T}})},$

where, for $p\in[1,\infty)$, the norm $\|\cdot\|_{L^{p}({\mathbb{T}})}$
denotes the usual $L^{p}({\mathbb{T}})$ norm. In the case $p=\infty$, we
denote $\|\cdot\|_{L^{\infty}({\mathbb{T}})}$ the $L^{\infty}({\mathbb{T}})$
norm given by

$\displaystyle\|w\|_{L^{\infty}({\mathbb{T}})}=\max_{i\in{\mathcal{I}}}|w_{i}|,\quad\forall
w\in\mathcal{H}_{\mathcal{T}}.$

Let us also recall the definition of the space $BV({\mathbb{T}})$, see [1] for
more details. A function $w\in L^{1}({\mathbb{T}})$ belongs to the space
$BV({\mathbb{T}})$ if its total variation $TV(w)$ given by

$\displaystyle
TV(w)=\sup\left\\{\int_{\mathbb{T}}w(x)\,\partial_{x}\phi(x)\,\mathrm{d}x,\quad\phi\in
C^{1}_{c}({\mathbb{T}}),\quad|\phi(x)|\leq 1\quad\forall
x\in{\mathbb{T}}\right\\},$

is finite. We endow the space $BV({\mathbb{T}})$ with the norm

$\|w\|_{BV({\mathbb{T}})}=\|w\|_{L^{1}({\mathbb{T}})}+TV(w),\quad\forall w\in
BV({\mathbb{T}}).$

In particular, we notice that for each function
$w\in\mathcal{H}_{\mathcal{T}}\cap BV({\mathbb{T}})$ we have
$\|w\|_{BV({\mathbb{T}})}=\|w\|_{1,1,{\mathcal{T}}}$.

Finally we introduce the space $\mathcal{H}_{\mathcal{D}}$ of piecewise
constant in time functions with values in $\mathcal{H}_{\mathcal{T}}$,

$\displaystyle\mathcal{H}_{{\mathcal{D}}}=\Bigg{\\{}w:{\mathbb{T}}\times[0,T]\to{\mathbb{R}}\,\,:\,\,w(x,t)=\sum_{k=1}^{N_{T}}w^{k}(x)\mathbf{1}_{(t_{k-1},t_{k}]}(t)\Bigg{\\}}.$

This space can be equipped, for $(p,q)\in[1,\infty)^{2}$, with the following
discrete $L^{q}(0,T;W^{1,p}({\mathbb{T}}))$ norm

$\displaystyle\left(\sum_{k=1}^{N_{T}}\Delta
t\,\|w\|^{q}_{1,p,{\mathcal{T}}}\right)^{\frac{1}{q}}\quad\forall
w\in\mathcal{H}_{\mathcal{D}},$

or with the $L^{q}(0,T;L^{p}({\mathbb{T}}))$ norm

$\displaystyle\left(\sum_{k=1}^{N_{T}}\Delta
t\,\|w\|^{q}_{L^{p}({\mathbb{T}})}\right)^{\frac{1}{q}},\quad\forall
w\in\mathcal{H}_{\mathcal{D}}.$

In particular, in the case $p=q=2$, the $L^{2}(Q_{T})$ norm can also be
defined by duality as

(10)
$\displaystyle\|w\|^{2}_{L^{2}(Q_{T})}=\sup\left\\{\int_{0}^{T}\int_{{\mathbb{T}}}wf\,\mathrm{d}x\mathrm{d}t\,\,:\,\,f\in\mathcal{H}_{{\mathcal{D}}},\,\,\sum_{k=1}^{N_{T}}\Delta
t\sum_{i\in{\mathcal{I}}}\Delta x|f^{k}_{i}|^{2}=1\right\\},\quad\forall
w\in\mathcal{H}_{\mathcal{D}}.$

This dual formulation of $\|\cdot\|_{L^{2}(Q_{T})}$ will be needed later on.

###  2.2. Numerical scheme

We discretize the initial conditions (3) as

(11) $\displaystyle u^{0}_{j,i}=\frac{1}{\Delta
x}\int_{K_{i}}u_{j}^{0}(x)\,\mathrm{d}x,\quad\forall
i\in{\mathcal{I}},\,j=1,2.$

Now for given $(u^{k-1}_{1},u^{k-1}_{2})\in{\mathbb{R}}^{2N}$, the implicit in
time numerical scheme writes as

(12) $\displaystyle\frac{u^{k}_{j,i}-u^{k-1}_{j,i}}{\Delta
t}-\left(\Delta_{\mathcal{T}}(\mu_{j}^{k}u_{j}^{k})\right)_{i}=R_{j}(u^{k}_{1,i},u^{k}_{2,j}),\quad\forall
i\in{\mathcal{I}},\,j=1,2,$

where $\Delta_{\mathcal{T}}$ denotes the discrete Laplacian, namely

(13)
$\displaystyle\left(\Delta_{\mathcal{T}}(\mu_{j}^{k}u_{j}^{k})\right)_{i}=\frac{\mu^{k}_{j,i+1}u^{k}_{j,i+1}-2\mu^{k}_{j,i}u^{k}_{j,i}+\mu^{k}_{j,i-1}u^{k}_{j,i-1}}{\Delta
x^{2}},\quad\forall i\in{\mathcal{I}},\,j=1,2,$

and

(14) $\displaystyle\mu_{1,i}^{k}$
$\displaystyle=d_{1}+d_{11}\sum_{n\in{\mathcal{I}}}\Delta
x\sigma_{1,i-n}u^{k}_{1,n}+d_{12}\sum_{n\in{\mathcal{I}}}\Delta
x\rho_{1,i-n}\,u^{k}_{2,n},\quad\forall i\in{\mathcal{I}},$ (15)
$\displaystyle\mu_{2,i}^{k}$
$\displaystyle=d_{2}+d_{21}\sum_{n\in{\mathcal{I}}}\Delta
x\rho_{2,i-n}\,u^{k}_{1,n}+d_{22}\sum_{n\in{\mathcal{I}}}\Delta
x\sigma_{2,i-n}u^{k}_{2,n},\quad\forall i\in{\mathcal{I}},$

with

(16) $\displaystyle\rho_{j,i-n}=\frac{1}{\Delta
x}\int_{K_{i-n}}\rho_{j}(y)\,\mathrm{d}y,\quad\sigma_{j,i-n}=\frac{1}{\Delta
x}\int_{K_{i-n}}\sigma_{j}(y)\,\mathrm{d}y,\quad\forall
i,\,n\in{\mathcal{I}},\,j=1,2.$

Let us notice, by construction, that if we consider
$\rho_{1}=\rho_{2}=\sigma_{1}=\sigma_{2}=\delta_{0}$ then (11)–(15) yields a
finite volume scheme for the local SKT model (4)–(5). We also remark that we
could equivalently rewrite (12) as

(17) $\displaystyle\Delta x\frac{u^{k}_{j,i}-u^{k-1}_{j,i}}{\Delta
t}+\mathcal{F}^{k}_{j,i+\frac{1}{2}}-\mathcal{F}^{k}_{j,i-\frac{1}{2}}=R_{j}(u^{k}_{1,i},u^{k}_{2,i}),\quad\forall
i\in{\mathcal{I}},\,j=1,2,$

where for all $i\in{\mathcal{I}}$ the numerical fluxes
$\mathcal{F}^{k}_{j,i+\frac{1}{2}}$ are defined by

(18)
$\mathcal{F}^{k}_{j,i+\frac{1}{2}}=\frac{u^{k}_{j,i}\mu^{k}_{j,i}-u^{k}_{j,i+1}\mu^{k}_{j,i+1}}{\Delta
x}=\mu^{k}_{j,i+\frac{1}{2}}\frac{u^{k}_{j,i}-u^{k}_{j,i+1}}{\Delta
x}+u^{k}_{j,i+\frac{1}{2}}\frac{\mu^{k}_{j,i}-\mu^{k}_{j,i+1}}{\Delta x},\quad
j=1,2,$

with the centered approximation at interfaces

$\displaystyle\mu^{k}_{j,i+\frac{1}{2}}=\frac{\mu^{k}_{j,i}+\mu^{k}_{j,i+1}}{2},\quad
u^{k}_{j,i+\frac{1}{2}}=\frac{u^{k}_{j,i}+u^{k}_{j,i+1}}{2},\quad j=1,2.$

###  2.3. Main results

Let us collect our assumptions.

(A44)

The domain is taken as $\Omega={\mathbb{T}}$.

The diffusion coefficients $d_{1}$, $d_{2}$, $d_{11}$ and $d_{22}$ are non-
negative constants and the cross-diffusion coefficients $d_{12}$ and $d_{21}$
are positive constants.

The convolution kernels $\rho_{1}$, $\rho_{2}$, $\sigma_{1}$ and $\sigma_{2}$
are $L^{\infty}({\mathbb{T}})$ functions that are non-negative and satisfy the
symmetry hypotheses (6). In particular $\rho_{i}=\rho_{1,i}=\rho_{2,-i}$ for
all $i\in{\mathcal{I}}$.

The initial data $u^{0}_{1}$ and $u^{0}_{2}$ are non-negative
$L^{1}({\mathbb{T}})$ functions with finite entropy, namely
$h_{1}(u_{1}^{0}),h_{2}(u_{2}^{0})\in L^{1}({\mathbb{T}})$.

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
$\Omega=\prod_{i=1}^{d}\mathbb{R}/L_{i}\mathbb{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 $\rho_{1}$,
$\rho_{2}$, $\sigma_{1}$ and $\sigma_{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 Section 6.4). Following for
instance [32], 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:

$\displaystyle
R_{j}(u_{1},u_{2})=u_{j}\,\left(a_{j0}-\sum_{k=1}^{2}a_{jk}\,u_{k}\right),\quad
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

$\displaystyle H(u^{k}_{1},u^{k}_{2})=\sum_{i\in{\mathcal{I}}}\Delta
x\,h_{1}(u^{k}_{1,i})+\sum_{i\in{\mathcal{I}}}\Delta x\,h_{2}(u^{k}_{2,i}),$

where the functions $h_{1}$ and $h_{2}$ are defined by

(19) $\displaystyle
h_{1}(x)=\dfrac{1}{d_{12}}\,\left(x\,(\log(x)-1)+1\right),\quad
h_{2}(x)=\dfrac{1}{d_{21}}\,\left(x\,(\log(x)-1)+1\right),\quad\forall
x\in(0,+\infty),$

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\frac{d_{11}}{d_{12}}\sum_{j\in{\mathcal{I}}}\sigma_{1,j}\sum_{i\in{\mathcal{I}}}\left(\sqrt{u^{k}_{1,i+1}u^{k}_{1,i+1-j}}-\sqrt{u^{k}_{1,i}u^{k}_{1,i-j}}\right)^{2}\\\
+2\frac{d_{22}}{d_{21}}\sum_{j\in{\mathcal{I}}}\sigma_{2,j}\sum_{i\in{\mathcal{I}}}\left(\sqrt{u^{k}_{2,i+1}u^{k}_{2,i+1-j}}-\sqrt{u^{k}_{2,i}u^{k}_{2,i-j}}\right)^{2}+4\frac{d_{1}}{d_{12}}\Big{|}\sqrt{u^{k}_{1}}\Big{|}^{2}_{1,2,{\mathcal{T}}}\\\
+4\frac{d_{2}}{d_{21}}\Big{|}\sqrt{u^{k}_{2}}\Big{|}^{2}_{1,2,{\mathcal{T}}}+4\sum_{j\in{\mathcal{I}}}\rho_{j}\sum_{i\in{\mathcal{I}}}\left(\sqrt{u^{k}_{1,i+1}u^{k}_{2,i+1-j}}-\sqrt{u^{k}_{1,i}u^{k}_{2,i-j}}\right)^{2}.$

######  Theorem 1 (Existence of solutions).

Let the assumptions (H1)–(H5) hold. Then, for every $1\leq k\leq N_{T}$ there
exists (at least) one nonnegative solution $(u^{k}_{1},u^{k}_{2})$ to scheme
(11)–(15). Moreover, this solution satisfies the following properties:

  * (i)

Mass conservation:

(21) $\displaystyle\sum_{i\in{\mathcal{I}}}\Delta
x\,u^{k}_{j,i}=\int_{{\mathbb{T}}}u^{0}_{j}(x)\,\mathrm{d}x,\quad\forall k\geq
0,\,j=1,2.$

  * (ii)

Entropy production estimate: for all $k\geq 1$ it holds

(22) $\displaystyle H(u^{k}_{1},u^{k}_{2})+\Delta tD(u^{k}_{1},u^{k}_{2})\leq
H(u^{k-1}_{1},u^{k-1}_{2}).$

Finally for all $j\in\\{1,2\\}$, if $d_{j}$ is positive then $u_{j}^{k}$ is
positive for all $0<k\leq N_{T}$.

The proof of existence of Theorem 1 is based on a consequence (see [24,
Section 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 [9, 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 [21, 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}\in L^{2}({\mathbb{T}})$ and let $\gamma$ and $\Gamma$ be
some nonnegative constants such that

$\displaystyle\gamma\leq u^{0}_{1}(x),\,u^{0}_{2}(x)\leq\Gamma\quad\mbox{a.e.
}x\in{\mathbb{T}}.$

Finally let us introduce $m^{0}_{j}=\|u^{0}_{j}\|_{L^{1}({\mathbb{T}})}$ for
$j=1,2$. Then the following properties hold.

  * (i)

Maximum principle: If $\rho$, $\sigma_{1}$ and $\sigma_{2}$ are twice
continuously differentiable functions and that the time step satisfies the
condition

$\displaystyle\Delta
t<1/\left(\min\\{d_{11}m^{0}_{1}\|\Delta\sigma_{1}\|_{L^{\infty}({\mathbb{T}})},d_{22}m^{0}_{2}\|\Delta\sigma_{2}\|_{L^{\infty}({\mathbb{T}})}\\}+\min\\{d_{12}m^{0}_{2},d_{21}m^{0}_{1}\\}\|\Delta\rho\|_{L^{\infty}({\mathbb{T}})}\right),$

then for all $i\in{\mathcal{I}}$, $k\geq 1$ and $j=1,2$ we have

$\displaystyle e_{k}\leq u^{k}_{j,i}\leq E_{k},$

where

$\displaystyle e_{k}$ $\displaystyle=\gamma\left(1+\Delta
t\left(\min\\{d_{11}m^{0}_{1}\|\Delta\sigma_{1}\|_{L^{\infty}({\mathbb{T}})},d_{22}m^{0}_{2}\|\Delta\sigma_{2}\|_{L^{\infty}({\mathbb{T}})}\\}+\min\\{d_{12}m^{0}_{2},d_{21}m^{0}_{1}\\}\|\Delta\rho\|_{L^{\infty}({\mathbb{T}})}\right)\right)^{-k},$
$\displaystyle E_{k}$ $\displaystyle=\Gamma\left(1-\Delta
t\left(\min\\{d_{11}m^{0}_{1}\|\Delta\sigma_{1}\|_{L^{\infty}({\mathbb{T}})},d_{22}m^{0}_{2}\|\Delta\sigma_{2}\|_{L^{\infty}({\mathbb{T}})}\\}+\min\\{d_{12}m^{0}_{2},d_{21}m^{0}_{1}\\}\|\Delta\rho\|_{L^{\infty}({\mathbb{T}})}\right)\right)^{-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

$\displaystyle\sum_{k=1}^{N_{T}}\Delta t\sum_{i\in{\mathcal{I}}}\Delta
x\left(\mu^{k}_{1,i}u^{k}_{1,i}+\mu^{k}_{2,i}u^{k}_{2,i}\right)\left(u^{k}_{1,i}+u^{k}_{2,i}\right)\leq
C(1+TA)\left(\|u^{0}_{1}\|_{L^{2}({\mathbb{T}})}^{2}+\|u^{0}_{2}\|_{L^{2}({\mathbb{T}})}^{2}\right),$

where

$\displaystyle
A=d_{1}+d_{2}+d_{11}m^{0}_{1}\,\|\sigma_{1}\|_{L^{1}({\mathbb{T}})}+d_{22}m^{0}_{2}\,\|\sigma_{2}\|_{L^{1}({\mathbb{T}})}+\|\rho\|_{L^{1}({\mathbb{T}})}\left(d_{12}m^{0}_{2}+d_{21}m^{0}_{1}\right).$

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 $(\mathcal{D}_{m})_{m\in{\mathbb{N}}}$ of space-time
discretizations of $Q_{T}$ indexed by the size $\eta_{m}=\max\\{\Delta
x_{m},\Delta t_{m}\\}$ of the mesh, satisfying $\eta_{m}\to 0$ as
$m\to\infty$. We denote by ${\mathcal{T}}_{m}$ the corresponding mesh of
${\mathbb{T}}$ and by $\Delta t_{m}$ the corresponding time step. Finally, for
every $m\in{\mathbb{N}}$ we set
$(u_{1,m},u_{2,m})\in\mathcal{H}_{\mathcal{D}_{m}}$ the picewise constant in
space and time reconstruction of the solutions to the scheme (11)–(15)
corresponding to the mesh ${\mathcal{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 $({\mathcal{D}}_{m})_{m\in{\mathbb{N}}}$ be
a family of space-time discretizations of $Q_{T}$ with $\eta_{m}\to 0$ as
$m\to\infty$. 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})\in\left(L^{p}(Q_{T})\right)^{2}$ for $p\in[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

$\displaystyle u_{j,m}\to u_{j}\quad\mbox{strongly in }L^{p}(Q_{T})\mbox{ for
}1\leq p<3\quad\mbox{as }m\to\infty.$

The proof of Theorem 3 is based on uniform estimates w.r.t. $\Delta x$ and
$\Delta t$, 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 $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$,

$\displaystyle u_{j,m}\rightarrow u_{j}\quad\mbox{strongly in
}L^{p}(Q_{T})\mbox{ for }1\leq p<{\frac{d}{d-1}},\quad\mbox{as }m\to\infty,$

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

##  3\. 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})\in((0,+\infty)^{N})^{2}$ we define the new unknown
$X=\Phi(u_{1},u_{2})$ where $\Phi:((0,+\infty)^{N})^{2}\to{\mathbb{R}}^{2N}$
is the smooth diffeomorphism defined by

$\Phi(u_{1},u_{2})=(d_{12}^{-1}\log(u_{1,1}),\dots,d_{12}^{-1}\log(u_{1,N}),d_{21}^{-1}\log(u_{2,1}),\dots,d_{21}^{-1}\log(u_{2,N}))^{\top}\in{\mathbb{R}}^{2N}.$

From there finding a positive solution to the scheme (11)–(15) amounts to
finding a zero $X^{k}=\Phi(u_{1}^{k},u_{2}^{k})$ of the continuous map
$P^{k}:{\mathbb{R}}^{2N}\to{\mathbb{R}}^{2N}$ defined for any
$(u_{1},u_{2})\in((0,+\infty)^{N})^{2}$ by its components

$P_{i+N(j-1)}^{k}(\Phi(u_{1},u_{2}))=\Delta x(u_{j,i}-u^{k-1}_{j,i})-\Delta
t\Delta x\left(\Delta_{\mathcal{T}}(\mu_{j}u_{j})\right)_{i},\quad\forall
i\in\\{1,\dots,N\\},j\in\\{1,2\\},$

where $(u_{1}^{k-1},u_{2}^{k-1})$ are given and $\mu_{1},\mu_{2}$ are related
to $u_{1},u_{2}$ through the relation (14) and (15) dropping the exponent $k$.

###  3.1. Entropy dissipation and mass conservation

In the following $\left\langle,\right\rangle$ denotes the Euclidean scalar
product and $|\cdot|$ the Euclidean norm.

######  Proposition 5.

Let $(u_{1}^{k-1},u_{2}^{k-1})$ be componentwise non-negative. Then for any
$X\in{\mathbb{R}}^{2N}$,

(23) $\left\langle
P^{k}(X),\mathds{1}_{j}\right\rangle=\sum_{i\in{\mathcal{I}}}(u_{j,i}-u^{k-1}_{j,i})\Delta
x\,,\quad\forall j\in\\{1,2\\},$

and

(24) $\left\langle P^{k}(X),X\right\rangle\geq
H(u_{1},u_{2})-H(u_{1}^{k-1},u_{2}^{k-1})+\Delta tD(u_{1},u_{2}),$

where $\mathds{1}_{1}=(1,\dots,1,0,\dots,0)$,
$\mathds{1}_{2}=(0,\dots,0,1,\dots,1)$, $(u_{1},u_{2})=\Phi^{-1}(X)$ and
$D(u_{1},u_{2})$ denotes the entropy dissipation functional given by (20).

###### Proof.

In order to prove (23), it suffices to sum the components of $P^{k}(X)$ and
observe that

$\sum_{i\in{\mathcal{I}}}\left(\Delta_{\mathcal{T}}(\mu_{j}u_{j})\right)_{i}=0,$

since it is a telescopic sum. Concerning the inequality, first observe that
$H(u_{1}^{k-1},u_{2}^{k-1})$ is well-defined since $(u_{1}^{k-1},u_{2}^{k-1})$
is non-negative. Then, using the definition of $\Phi(u_{1},u_{2})$ and $P^{k}$
one obtains

$\left\langle P^{k}(X),X\right\rangle=\Delta t(I_{1}+I_{2}+J_{1}+J_{2})$

with

$\displaystyle I_{1}$
$\displaystyle=\frac{1}{d_{12}}\sum_{i\in{\mathcal{I}}}\Delta
x(u_{1,i}-u^{k-1}_{1,i})\,\log(u_{1,i}),$ $\displaystyle I_{2}$
$\displaystyle=\frac{1}{d_{21}}\sum_{i\in{\mathcal{I}}}\Delta
x(u_{2,i}-u^{k-1}_{2,i})\,\log(u_{2,i}),$ $\displaystyle J_{1}$
$\displaystyle=\frac{1}{d_{12}\Delta
x}\sum_{i\in{\mathcal{I}}}(-\mu_{1,i+1}u_{1,i+1}+2\mu_{1,i}u_{1,i}-\mu_{1,i-1}u_{1,i-1})\,\log(u_{1,i})$
$\displaystyle J_{2}$ $\displaystyle=\frac{1}{d_{21}\Delta
x}\sum_{i\in{\mathcal{I}}}(-\mu_{2,i+1}u_{2,i+1}+2\mu_{2,i}u_{2,i}-\mu_{2,i-1}u_{2,i-1})\,\log(u_{2,i})$

Using the convexity of $x\mapsto(x\log(x)-x+1)$ to bound both $I_{1}$ and
$I_{2}$ from below, one obtains

$I_{1}+I_{2}\geq H(u_{1},u_{2})-H(u_{1}^{k-1},u_{2}^{k-1}).$

Then for $J_{1}$, a discrete integration by parts (or summation by parts)
yields

$J_{1}=\frac{1}{d_{12}\,\Delta
x}\sum_{i\in{\mathcal{I}}}\left(u_{1,i+1}\mu_{1,i+1}-u_{1,i}\mu_{1,i}\right)\left(\log(u_{1,i+1})-\log(u_{1,i})\right),$

and a similar formula holds for $J_{2}$. Using the definitions of $\mu_{1}$
and $\mu_{2}$ (see (14) and (15) without the exponents), one has
$J_{1}=J_{1}^{\text{diff}}+J_{1}^{\sigma_{1}}+J_{1}^{\rho_{1}}$ and
$J_{2}=J_{2}^{\text{diff}}+J_{2}^{\sigma_{2}}+J_{2}^{\rho_{2}}$ with

$\displaystyle J_{1}^{\text{diff}}$ $\displaystyle=\frac{d_{1}}{d_{12}\Delta
x}\sum_{i\in{\mathcal{I}}}\left(u_{1,i+1}-u_{1,i}\right)\left(\log(u_{1,i+1})-\log(u_{1,i})\right)\geq\frac{4d_{1}}{d_{12}\Delta
x}\sum_{i\in{\mathcal{I}}}\left(\sqrt{u_{1,i+1}}-\sqrt{u_{1,i}}\right)^{2},$

and a similar estimate for $J_{2}^{\text{diff}}$. For the second term one has

$\displaystyle J_{1}^{\sigma_{1}}$
$\displaystyle=\frac{d_{11}}{d_{12}}\sum_{j\in{\mathcal{I}}}\sigma_{1,j}\sum_{i\in{\mathcal{I}}}\left(u_{1,i+1}u_{1,i+1-j}-u_{1,i}u_{1,i-j}\right)\left(\log(u_{1,i+1})-\log(u_{1,i})\right)$
$\displaystyle=\frac{1}{2}J_{1}^{\sigma_{1}}+\frac{d_{11}}{2d_{12}}\sum_{j\in{\mathcal{I}}}\sigma_{1,j}\sum_{i\in{\mathcal{I}}}\left(u_{1,i+1}u_{1,i+1+j}-u_{1,i}u_{1,i+j}\right)\left(\log(u_{1,i+1})-\log(u_{1,i})\right)$
$\displaystyle=\frac{1}{2}J_{1}^{\sigma_{1}}+\frac{d_{11}}{2d_{12}}\sum_{j\in{\mathcal{I}}}\sigma_{1,j}\sum_{i\in{\mathcal{I}}}\left(u_{1,i+1-j}u_{1,i+1}-u_{1,i-j}u_{1,i}\right)\left(\log(u_{1,i+1-j})-\log(u_{1,i-j})\right)$
$\displaystyle=\frac{d_{11}}{2d_{12}}\sum_{j\in{\mathcal{I}}}\sigma_{1,j}\sum_{i\in{\mathcal{I}}}\left(u_{1,i+1}u_{1,i+1-j}-u_{1,i}u_{1,i-j}\right)\left(\log(u_{1,i+1}u_{1,i+1-j})-\log(u_{1,i}u_{1,i-j})\right)$
$\displaystyle\geq\frac{2d_{11}}{d_{12}}\sum_{j\in{\mathcal{I}}}\sigma_{1,j}\sum_{i\in{\mathcal{I}}}\left(\sqrt{u_{1,i+1}u_{1,i+1-j}}-\sqrt{u_{1,i}u_{1,i-j}}\right)^{2}.$

In the previous estimate, the second inequality is obtained by changing $j$
into $-j$ and using the symmetry of $\sigma_{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}^{\sigma_{2}}$. Finally with the same changes of indices one can
estimate the sum

$\displaystyle J_{1}^{\rho_{1}}+J_{2}^{\rho_{2}}=$
$\displaystyle\sum_{j\in{\mathcal{I}}}\rho_{j}\sum_{i\in{\mathcal{I}}}\left(u_{1,i+1}u_{2,i+1-j}-u_{1,i}u_{2,i-j}\right)\left(\log(u_{1,i+1})-\log(u_{1,i})\right)$
$\displaystyle+\sum_{j\in{\mathcal{I}}}\rho_{-j}\sum_{i\in{\mathcal{I}}}\left(u_{2,i+1}u_{1,i+1-j}-u_{2,i}u_{1,i-j}\right)\left(\log(u_{2,i+1})-\log(u_{2,i})\right)$
$\displaystyle=$
$\displaystyle\sum_{j\in{\mathcal{I}}}\rho_{j}\sum_{i\in{\mathcal{I}}}\left(u_{1,i+1}u_{2,i+1-j}-u_{1,i}u_{2,i-j}\right)\left(\log(u_{1,i+1})-\log(u_{1,i})\right)$
$\displaystyle+\sum_{j\in{\mathcal{I}}}\rho_{j}\sum_{i\in{\mathcal{I}}}\left(u_{2,i-j+1}u_{1,i+1}-u_{2,i-j}u_{1,i}\right)\left(\log(u_{2,i-j+1})-\log(u_{2,i-j})\right)$
$\displaystyle=$
$\displaystyle\sum_{j\in{\mathcal{I}}}\rho_{j}\sum_{i\in{\mathcal{I}}}\left(u_{1,i+1}u_{2,i+1-j}-u_{1,i}u_{2,i-j}\right)\left(\log(u_{1,i+1}u_{2,i-j+1})-\log(u_{1,i}u_{2,i-j})\right)$
$\displaystyle\geq$
$\displaystyle\sum_{j\in{\mathcal{I}}}\rho_{j}\sum_{i\in{\mathcal{I}}}\left(\sqrt{u_{1,i+1}u_{2,i+1-j}}-\sqrt{u_{1,i}u_{2,i-j}}\right)^{2}.$

By summing all the estimates one obtains (24). The last point of the
proposition is obtained by induction. ∎

###  3.2. Proof of Theorem 1

Let us show that $P^{k}\circ\Phi$ has at least one zero. We use an
approximation argument by introducing

$P_{\varepsilon}^{k}(X)=P^{k}(X)+\varepsilon X,\quad\forall
X\in{\mathbb{R}}^{2N}\,.$

Using (24) and the non-negativity of the entropy and the entropy dissipation
one has

$\left\langle
P^{k}_{\varepsilon}(X),X\right\rangle\geq\varepsilon|X|^{2}-H(u_{1}^{k-1},u_{2}^{k-1})\,.$

Therefore, as a consequence of Brouwer fixed point theorem (see [24, Section
9.1] for details), there is $X_{\varepsilon}$ such that

$P^{k}_{\varepsilon}(X_{\varepsilon})=0\quad\text{and}\quad|X_{\varepsilon}|^{2}\leq
H(u_{1}^{k-1},u_{2}^{k-1})\varepsilon^{-1}\,.$

Let us define the associated
$(u_{1}^{\varepsilon},u_{2}^{\varepsilon})=\Phi^{-1}(X_{\varepsilon})$, which
is componentwise positive by definition. Observe that

$\displaystyle\Delta x\,h_{j}(u_{j,i}^{\varepsilon})\leq
H(u_{1}^{\varepsilon},u^{\varepsilon}_{2})\leq
H(u^{k-1}_{1},u^{k-1}_{2}),\quad\forall i\in{\mathcal{I}},\,j=1,2.$

where the last inequality is again a consequence of (24) for
$X=X_{\varepsilon}$. This shows that for any
$(u_{1,i}^{\varepsilon},u_{2,i}^{\varepsilon})$ is uniformly bounded in
$\varepsilon$. Therefore, there exists a subsequence (not relabeled) such that
$u_{j,i}^{\varepsilon}\to u_{j,i}^{k}\geq 0$ as $\varepsilon\to 0$, for every
$i\in{\mathcal{I}}$ and $j=1,2$. Since
$|X_{\varepsilon}|=O(\varepsilon^{-1/2})$ one has

$0=\lim_{\varepsilon\to
0}P^{k}_{\varepsilon}(X_{\varepsilon})=P^{k}(\Phi(u_{1}^{k},u_{2}^{k}))$

Therefore $(u_{1}^{k},u_{2}^{k})$ solves the scheme (11)-(15). By taking
limits in (23) and (24) evaluated at $X_{\varepsilon}$ as $\varepsilon\to 0$
one recovers (21) and (22) respectively.

Let us finally prove that if $d_{j}>0$, $u^{k}_{j,i}>0$ for all
$i\in{\mathcal{I}}$ and $0<k\leq N_{T}$. This is a consequence of the entropy
estimate. For a given $0<k\leq N_{T}$ we notice that (thanks to the term
$J^{\mathrm{diff}}_{j}$) the positive solution
$(u^{\varepsilon}_{1},u^{\varepsilon}_{2})$ satisfies the following estimate

$\displaystyle\frac{d_{j}\,\Delta t}{\Delta
x}\sum_{i\in{\mathcal{I}}}\left(u^{\varepsilon}_{j,i+1}-u^{\varepsilon}_{j,i}\right)\,\left(\log\left(u^{\varepsilon}_{j,i+1}\right)-\log\left(u^{\varepsilon}_{j,i}\right)\right)\leq\max(d_{12},d_{21})\,H\left(u^{k-1}_{1},u^{k-1}_{2}\right)\quad
j=1,2.$

At the limit $\varepsilon\to 0$, let us assume by contradiction that there
exists $i\in{\mathcal{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\in{\mathcal{I}}$. Consequently we have
$\|u^{k}_{j}\|_{L^{1}({\mathbb{T}})}=0$ which contradicts the mass
conservation property (21). This completes the proof of Theorem 1.

##  4\. Estimates on the discrete Kolmogorov equation

In this section, we focus on estimates concerning the finite volume
discretization of the Kolmogorov equation $\partial_{t}z=\Delta(\mu z)$. In
particular we adapt at the discrete level some properties established in [37,
21].

In the rest of this section, we assume that
$(\mu^{k}_{i})_{i\in{\mathcal{I}}}$, $k=1,\dots,N_{T}$ is given and
componentwise non-negative. From there, the scheme is given for all $k\geq 1$
by

(25) $\displaystyle\frac{z^{k}_{i}-z^{k-1}_{i}}{\Delta
t}-\left(\Delta_{\mathcal{T}}(\mu^{k}z^{k})\right)_{i}=0,\quad\forall
i\in{\mathcal{I}},$

where $\Delta_{\mathcal{T}}$ denotes the discrete Laplacian operator defined
by (13).

###  4.1. Well-posedness of the scheme and $L^{\infty}$ 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\in{\mathcal{I}}}$ there is a unique
$(z^{k}_{i})_{i\in{\mathcal{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},\ldots,z^{k-1}_{N-1})^{\top}$ for all
$k\geq 1$. Observe that the scheme writes $\mathbb{M}^{k}Z^{k}=Z^{k-1}$ where
$\mathbb{M}^{k}$ is a $N\times N$ tridiagonal matrix defined by

(26) $\displaystyle{\mathbb{M}}^{k}_{i,i-1}=-\frac{\Delta t}{\Delta
x^{2}}\mu^{k}_{i-1},\quad{\mathbb{M}}^{k}_{i,i}=1+2\frac{\Delta t}{\Delta
x^{2}}\mu^{k}_{i},\quad{\mathbb{M}}^{k}_{i,i+1}=-\frac{\Delta t}{\Delta
x^{2}}\mu^{k}_{i+1},\quad\forall i\in{\mathcal{I}}.$

We notice that ${\mathbb{M}}^{k}$ has positive diagonal terms and non-positive
off-diagonal terms. Furthermore the matrix ${\mathbb{M}}^{k}$ is strictly
diagonally dominant with respect to its columns. Therefore ${\mathbb{M}}^{k}$
is a non-singular M-matrix and is thus monotone and invertible. This finishes
the proof of Lemma 6. ∎

We prove in the following result some $L^{\infty}$ estimates for the solution
to scheme (25).

######  Lemma 7.

Let us assume that there exists $\tilde{\gamma},\tilde{\Gamma}\geq 0$ such
that

$\displaystyle\tilde{\gamma}\leq z^{0}_{i}\leq\tilde{\Gamma},\quad\forall
i\in{\mathcal{I}}.$

Then for every $k\geq 1$ and every $\Delta t>0$ such that

$\displaystyle\Delta t<1/\max_{1\leq k\leq
N_{T}}\|[\Delta_{\mathcal{T}}\mu^{k}]_{+}\|_{L^{\infty}({\mathbb{T}})},$

the solution $Z^{k}$ to (25) satisfies

(27) $\displaystyle\tilde{\gamma}\,\Pi_{n=1}^{k}\left(1+\Delta
t\|[\Delta_{\mathcal{T}}\mu^{n}]_{-}\|_{L^{\infty}({\mathbb{T}})}\right)^{-1}\leq
z^{k}_{i}\leq\tilde{\Gamma}\,\Pi_{n=1}^{k}\left(1-\Delta
t\|[\Delta_{\mathcal{T}}\mu^{n}]_{+}\|_{L^{\infty}({\mathbb{T}})}\right)^{-1},\,\,\,\forall
i\in{\mathcal{I}},$

where $[x]_{+}=\max(x,0)$ and $[x]_{-}=\min(x,0)$.

###### Proof.

We will only deal with the upper bound in (27) and the lower bound is obtained
in the same way. Let ${\mathbb{M}}^{k}$ denotes the tridiagonal matrix defined
by (26) and define

$\tilde{\Gamma}^{k}=\tilde{\Gamma}\,\Pi_{n=1}^{k}\left(1-\Delta
t\|[\Delta_{\mathcal{T}}\mu^{n}]_{+}\|_{L^{\infty}({\mathbb{T}})}\right)^{-1}.$

We proceed by induction. Since $\tilde{\Gamma}^{0}=\tilde{\Gamma}$ the bound
holds by hypothesis at $k=0$. Then observe that for every $i\in{\mathcal{I}}$

$\displaystyle({\mathbb{M}}^{k}(Z^{k}-\tilde{\Gamma}^{k}))_{i}=Z^{k-1}-\tilde{\Gamma}^{k}+\Delta
t\tilde{\Gamma}^{k}\,(\Delta_{\mathcal{T}}\mu^{k})_{i}\leq\tilde{\Gamma}^{k-1}-\tilde{\Gamma}^{k}+\Delta
t\tilde{\Gamma}^{k}\,(\Delta_{\mathcal{T}}\mu^{k})_{i}.$

Now we notice that by construction

$\displaystyle\tilde{\Gamma}^{k-1}-\tilde{\Gamma}^{k}=-\Delta
t\tilde{\Gamma}^{k}\,\|[\Delta_{\mathcal{T}}\mu^{k}]_{+}\|_{L^{\infty}({\mathbb{T}})}.$

Then we easily deduce that for every $i\in{\mathcal{I}}$ it holds

$\displaystyle({\mathbb{M}}^{k}(Z^{k}-\tilde{\Gamma}^{k}))_{i}=-\Delta
t\tilde{\Gamma}^{k}\left(\|[\Delta_{\mathcal{T}}\mu^{k}]_{+}\|_{L^{\infty}({\mathbb{T}})}-(\Delta_{\mathcal{T}}\mu^{k})_{i}\right)\leq
0.$

Therefore, since ${\mathbb{M}}^{k}$ is a M-matrix we conclude that
$z^{k}_{i}\leq\tilde{\Gamma}^{k}$ for all $i\in{\mathcal{I}}$ which concludes
the proof of Lemma 7. ∎

The bounds of Lemma 7 are exactly the discrete equivalent of the $L^{\infty}$
estimates established at the continuous level in [21, Corollary 18].

######  Proposition 8.

Let us assume that it holds

$\displaystyle\Delta t<1/\max_{1\leq k\leq
N_{T}}\|[\Delta_{\mathcal{T}}\mu^{k}]_{+}\|_{L^{\infty}({\mathbb{T}})}.$

Then the solution to (25) satisfies the following estimate

$\|z^{k}\|_{L^{2}({\mathbb{T}})}^{2}+\sum_{n=1}^{k}\Delta
t\sum_{i\in{\mathcal{T}}}(\mu^{n}_{i}+\mu^{n}_{i+1})\,\frac{(z^{n}_{i+1}-z^{n}_{i})^{2}}{\Delta
x}\\\ \leq\Pi_{n=1}^{k}\left(1-\Delta
t\|[\Delta_{\mathcal{T}}\mu^{n}]_{+}\|_{L^{\infty}({\mathbb{T}})}\right)^{-1}\,\|z^{0}_{i}\|_{L^{2}({\mathbb{T}})}^{2},\quad\forall
1\leq k\leq N_{T}.$

###### Proof.

Let $k\geq 1$ be fixed and let us first notice that we can rewrite for every
$i\in{\mathcal{I}}$ equation (25) as

$\displaystyle\Delta x\frac{z^{k}_{i}-z^{k-1}_{i}}{\Delta
t}+\mu^{k}_{i+\frac{1}{2}}\frac{(z^{k}_{i}-z^{k}_{i+1})}{\Delta
x}-\mu^{k}_{i-\frac{1}{2}}\frac{(z^{k}_{i-1}-z^{k}_{i})}{\Delta
x}+z^{k}_{i+\frac{1}{2}}\frac{(\mu^{k}_{i}-\mu^{k}_{i+1})}{\Delta
x}-z^{k}_{i-\frac{1}{2}}\frac{(\mu^{k}_{i-1}-\mu^{k}_{i})}{\Delta x}=0,$

where

$\displaystyle\mu^{k}_{i+\frac{1}{2}}=\frac{\mu^{k}_{i}+\mu^{k}_{i+1}}{2},\quad
z^{k}_{i+\frac{1}{2}}=\frac{z^{k}_{i}+z^{k}_{i+1}}{2},\quad\forall
i\in{\mathcal{I}}.$

Now we multiply the above equation by $\Delta tz^{k}_{i}$ and we sum over
$i\in{\mathcal{I}}$, we obtain

$\displaystyle I_{3}+I_{4}+I_{5}=0,$

where

$\displaystyle I_{3}$ $\displaystyle=\sum_{i\in{\mathcal{I}}}\Delta
x(z^{k}_{i}-z^{k-1}_{i})z^{k}_{i},$ $\displaystyle I_{4}$
$\displaystyle=\Delta
t\sum_{i\in{\mathcal{I}}}\left(\mu^{k}_{i+\frac{1}{2}}\frac{(z^{k}_{i}-z^{k}_{i+1})}{\Delta
x}-\mu^{k}_{i-\frac{1}{2}}\frac{(z^{k}_{i-1}-z^{k}_{i})}{\Delta
x}\right)z^{k}_{i},$ $\displaystyle I_{5}$ $\displaystyle=\Delta
t\sum_{i\in{\mathcal{I}}}\left(z^{k}_{i+\frac{1}{2}}\frac{(\mu^{k}_{i}-\mu^{k}_{i+1})}{\Delta
x}-z^{k}_{i-\frac{1}{2}}\frac{(\mu^{k}_{i-1}-\mu^{k}_{i})}{\Delta
x}\right)z^{k}_{i}.$

For $I_{3}$ using the inequality $(a-b)a\geq(a^{2}-b^{2})/2$ we obtain

(28) $\displaystyle I_{3}\geq\frac{1}{2}\sum_{i\in{\mathcal{I}}}\Delta
x\left(|z^{k}_{i}|^{2}-|z^{k-1}_{i}|^{2}\right).$

For $I_{4}$ applying a discrete integration by parts yields

(29) $\displaystyle I_{4}=\Delta
t\sum_{i\in{\mathcal{I}}}\mu^{k}_{i+\frac{1}{2}}\,\frac{(z^{k}_{i+1}-z^{k}_{i})^{2}}{\Delta
x}.$

Now we rewrite $I_{5}$ as

$\displaystyle I_{5}=-\frac{\Delta t}{2}\sum_{i\in{\mathcal{I}}}\Delta
x|z^{k}_{i}|^{2}(\Delta_{\mathcal{T}}\mu^{k})_{i}+\frac{\Delta t}{2\Delta
x}\sum_{i\in{\mathcal{I}}}\left(z^{k}_{i+1}z^{k}_{i}(\mu^{k}_{i}-\mu^{k}_{i+1})-z^{k}_{i-1}z^{k}_{i}(\mu^{k}_{i-1}-\mu^{k}_{i})\right),$

and reordering the terms in the r.h.s. the second sum vanishes and we have

(30) $\displaystyle I_{5}=-\frac{\Delta t}{2}\sum_{i\in{\mathcal{I}}}\Delta
x|z^{k}_{i}|^{2}(\Delta_{\mathcal{T}}\mu^{k})_{i}.$

Gathering (28)–(30) we end up with

$\displaystyle\frac{1}{2}\sum_{i\in{\mathcal{I}}}\Delta
x|z^{k}_{i}|^{2}+\Delta
t\sum_{i\in{\mathcal{I}}}\mu^{k}_{i+\frac{1}{2}}\frac{(z^{k}_{i+1}-z^{k}_{i})^{2}}{\Delta
x}\leq\frac{1}{2}\sum_{i\in{\mathcal{I}}}\Delta
x|z^{k-1}_{i}|^{2}+\frac{\Delta
t}{2}\|[\Delta_{\mathcal{T}}\mu^{k}]_{+}\|_{L^{\infty}({\mathbb{T}})}\sum_{i\in{\mathcal{I}}}\Delta
x|z^{k}_{i}|^{2}.$

We deduce that

$\frac{1}{2}\sum_{i\in{\mathcal{I}}}\Delta
x|z^{k}_{i}|^{2}+\sum_{n=1}^{k}\Delta
t\sum_{i\in{\mathcal{I}}}\mu^{n}_{i+\frac{1}{2}}\frac{(z^{n}_{i+1}-z^{n}_{i})^{2}}{\Delta
x}\\\ \leq\frac{1}{2}\sum_{i\in{\mathcal{I}}}\Delta
x|z^{0}_{i}|^{2}+\sum_{n=1}^{k}\frac{\Delta
t}{2}\|[\Delta_{\mathcal{T}}\mu^{n}]_{+}\|_{L^{\infty}({\mathbb{T}})}\sum_{i\in{\mathcal{I}}}\Delta
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. ∎

###  4.2. 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, Theorem 3]. In this aim, following [37], we introduce a “dual” scheme
associated to (25). Let $v^{N_{T}+1}_{i}$ be given for every
$i\in{\mathcal{I}}$, then for $1\leq k\leq N_{T}$ we want to determine the
solution to the following implicit backward in time scheme

(31) $\displaystyle\frac{v^{k}_{i}-v^{k+1}_{i}}{\Delta
t}-\mu^{k}_{i}\left(\Delta_{\mathcal{T}}v^{k}\right)_{i}=S^{k}_{i}\quad\forall
i\in{\mathcal{I}},$

where $\mu^{k}_{i}$ is given and non-negative and
$S^{k}=(S^{k}_{0},\ldots,S^{k}_{N-1})$ is some given vector in
${\mathbb{R}}^{N}$ for all $1\leq k\leq N_{T}$. Let us notice that (31) define
a set of linear equation which can be rewritten as

(32) $\displaystyle(\mathbb{M}^{k})^{\top}V^{k}=V^{k+1}+\Delta
tS^{k},\quad\forall 1\leq k\leq N_{T},$

where $\mathbb{M}^{k}$ is the tridiagonal matrix given by (26). Therefore, it
follows directly from the proof of Lemma 6 that the problem (32) admits a
unique solution for every $1\leq k\leq 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\in{\mathcal{I}}}\mu^{k}_{i}>0$ for every $0\leq k\leq
N_{T}$ and that $v^{N_{T}+1}_{i}=0$ for every $i\in{\mathcal{I}}$. Then the
solution to (32) satisfies for every $1\leq k\leq N_{T}$ the following
estimate

(33) $\displaystyle|v^{k}|^{2}_{1,2,{\mathcal{T}}}+\sum_{n=k}^{N_{T}}\Delta
t\sum_{i\in{\mathcal{I}}}\mu^{n}_{i}\frac{(\Delta_{\mathcal{T}}v^{n})^{2}_{i}}{\Delta
x}\leq\|\mu^{-1/2}S\|^{2}_{L^{2}(Q_{T})},$

and there exists a constant $C>0$ independent of $\Delta x$ such that

(34) $\displaystyle\|v^{k}\|^{2}_{L^{2}({\mathbb{T}})}\leq
C(1+\|\mu\|_{L^{1}(Q_{T})})\,\|\mu^{-1/2}S\|^{2}_{L^{2}(Q_{T})},\quad\forall
1\leq k\leq N_{T},$

where $\mu$ and $S$ denote the piecewise reconstruction functions in
$\mathcal{H}_{{\mathcal{D}}}$ associated to the vectors $(\mu^{k})_{1\leq
k\leq N_{T}}$ and $(S^{k})_{1\leq k\leq N_{T}}$.

###### Proof.

Let us first establish estimate (33). In this purpose let $1\leq k\leq N_{T}$
be fixed. We multiply equation (31) by $\Delta
t(-v^{k}_{i+1}+2v^{k}_{i}-v^{k}_{i-1})/\Delta x$, we sum over
$i\in{\mathcal{I}}$ and we apply definition (13) of the operator
$\Delta_{\mathcal{T}}$ and we obtain

$\displaystyle I_{6}+I_{7}=I_{8},$

with

$\displaystyle I_{6}$
$\displaystyle=\sum_{i\in{\mathcal{I}}}(v^{k}_{i}-v^{k+1}_{i})\frac{(-v^{k}_{i+1}+2v^{k}_{i}-v^{k}_{i-1})}{\Delta
x},$ $\displaystyle I_{7}$ $\displaystyle=\Delta
t\sum_{i\in{\mathcal{I}}}\mu^{k}_{i}\frac{(v^{k}_{i+1}-2v^{k}_{i}+v^{k}_{i-1})^{2}}{\Delta
x^{3}},$ $\displaystyle I_{8}$ $\displaystyle=\Delta
t\sum_{i\in{\mathcal{I}}}S^{k}_{i}\frac{(-v^{k}_{i+1}+2v^{k}_{i}-v^{k}_{i-1})}{\Delta
x}.$

For $I_{6}$ reordering the terms leads to

$\displaystyle I_{6}=\frac{1}{\Delta
x}\sum_{i\in{\mathcal{I}}}(v^{k}_{i+1}-v^{k}_{i})\left[(v^{k}_{i+1}-v^{k}_{i})-(v^{k+1}_{i+1}-v^{k+1}_{i})\right],$

and using the inequality $a(a-b)\geq(a^{2}-b^{2})/2$ we get

(35) $\displaystyle
I_{6}\geq\frac{1}{2}\sum_{i\in{\mathcal{I}}}\left[\frac{(v^{k}_{i+1}-v^{k}_{i})^{2}}{\Delta
x}-\frac{(v^{k+1}_{i+1}-v^{k+1}_{i})^{2}}{\Delta x}\right].$

For $I_{8}$ applying the Cauchy-Schwarz and Young inequality yield

(36) $\displaystyle|I_{8}|\leq\frac{\Delta t}{2}\sum_{i\in{\mathcal{I}}}\Delta
x(\mu^{k}_{i})^{-1}|S^{k}_{i}|^{2}+\frac{\Delta
t}{2}\sum_{i\in{\mathcal{I}}}\mu^{k}_{i}\frac{(-v^{k}_{i+1}+2v^{k}_{i}-v^{k}_{i-1})^{2}}{\Delta
x^{3}}.$

Collecting (35)–(36) we obtain

$\displaystyle\sum_{i\in{\mathcal{I}}}\frac{(v^{k}_{i+1}-v^{k}_{i})^{2}}{\Delta
x}+\Delta
t\sum_{i\in{\mathcal{I}}}\mu^{k}_{i}\frac{(\Delta_{\mathcal{T}}v^{k})^{2}_{i}}{\Delta
x}\leq\sum_{i\in{\mathcal{I}}}\frac{(v^{k+1}_{i+1}-v^{k+1}_{i})^{2}}{\Delta
x}+\Delta t\sum_{i\in{\mathcal{I}}}\Delta x(\mu^{k}_{i})^{-1}|S^{k}_{i}|^{2}.$

In order to prove (33) it remains to sum over $n\in\\{k,\ldots,N_{T}\\}$.

We now prove estimate (34). In this purpose we multiply (31) by $\Delta
x\Delta t$, we sum over $i\in{\mathcal{I}}$ and $n\in\\{k,\ldots,N_{T}\\}$ and
we obtain

$\displaystyle\sum_{i\in{\mathcal{I}}}\Delta
xv^{k}_{i}=\sum_{n=k}^{N_{T}}\Delta
t\sum_{i\in{\mathcal{I}}}\mu^{k}_{i}\frac{v^{n}_{i+1}-2v^{n}_{i}+v^{n}_{i-1}}{\Delta
x}+\sum_{n=k}^{N_{T}}\Delta t\sum_{i\in{\mathcal{I}}}\Delta xS^{n}_{i}.$

Applying the Cauchy-Schwarz inequality leads to

$\left|\sum_{i\in{\mathcal{I}}}\Delta
xv^{k}_{i}\right|\leq\left(\sum_{n=k}^{N_{T}}\Delta
t\sum_{i\in{\mathcal{I}}}\Delta
x\mu^{n}_{i}\right)^{1/2}\Bigg{[}\left(\sum_{n=k}^{N_{T}}\Delta
t\sum_{i\in{\mathcal{I}}}\mu^{n}_{i}\frac{(\Delta_{\mathcal{T}}v^{n})^{2}_{i}}{\Delta
x}\right)^{1/2}\\\ +\left(\sum_{n=k}^{N_{T}}\Delta
t\sum_{i\in{\mathcal{I}}}\Delta
x(\mu^{n}_{i})^{-1}|S^{n}_{i}|^{2}\right)^{1/2}\Bigg{]}.$

Using estimate (33) we obtain

$\displaystyle\left|\sum_{i\in{\mathcal{I}}}\Delta xv^{k}_{i}\right|\leq
2\left(\sum_{n=k}^{N_{T}}\Delta t\sum_{i\in{\mathcal{I}}}\Delta
x\mu^{n}_{i}\right)^{1/2}\left(\sum_{n=k}^{N_{T}}\Delta
t\sum_{i\in{\mathcal{I}}}\Delta
x(\mu^{n}_{i})^{-1}|S^{n}_{i}|^{2}\right)^{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\in{\mathcal{I}}}\mu^{k}_{i}>0$ for every $1\leq
k\leq N_{T}$. Then there exists a constant $C>0$ independent of $\Delta x$
such that the solution $(Z^{k})_{1\leq k\leq N_{T}}$ to (25) satisfies

$\displaystyle\|\mu^{1/2}z\|_{L^{2}(Q_{T})}\leq
C\left(1+\|\mu\|_{L^{1}(Q_{T})}^{1/2}\right)\,\|z^{0}\|_{L^{2}({\mathbb{T}})}.$

###### Proof.

Let $(v^{k}_{i})_{i\in{\mathcal{I}}}$ be given in ${\mathbb{R}}^{N}$ for every
$1\leq k\leq N_{T}+1$ with $v^{N_{T}+1}_{i}=0$ for all $i\in{\mathcal{I}}$.
Now for $1\leq k\leq N_{T}$, we multiply (25) by $\Delta t\Delta xv^{k}_{i}$,
we sum over $i\in{\mathcal{I}}$ and $k\in\\{1,\ldots,N_{T}\\}$, we obtain

$\displaystyle\sum_{k=1}^{N_{T}}\sum_{i\in{\mathcal{I}}}\Delta
x(z^{k}_{i}-z^{k-1}_{i})v^{k}_{i}-\sum_{k=1}^{N_{T}}\Delta
t\sum_{i\in{\mathcal{I}}}\Delta
x\left(\Delta_{\mathcal{T}}(z^{k}\mu^{k})\right)_{i}v^{k}_{i}=0.$

Reordering the terms we have

(37) $\displaystyle\sum_{k=1}^{N_{T}}\Delta t\sum_{i\in{\mathcal{I}}}\Delta
xz^{k}_{i}\left(\frac{(v^{k}_{i}-v^{k+1}_{i})}{\Delta
t}-\mu^{k}_{i}\left(\Delta_{\mathcal{T}}v^{k}\right)_{i}\right)=\sum_{i\in{\mathcal{I}}}\Delta
xz^{0}_{i}v^{1}_{i}.$

We define $(S^{k}_{i})_{i\in{\mathcal{I}}}$ by

$\displaystyle S^{k}_{i}=\frac{(v^{k}_{i}-v^{k+1}_{i})}{\Delta
t}-\mu^{k}_{i}\left(\Delta_{\mathcal{T}}v^{k}\right)_{i},\quad\forall
i\in{\mathcal{I}},\,1\leq k\leq N_{T}.$

We first notice that $(S^{k}_{i})_{i\in{\mathcal{I}}}$ is well-defined since
we know that equation (31) is well-posed. Besides applying the Cauchy-Schwarz
inequality in (37) we get

$\displaystyle\left|\sum_{k=1}^{N_{T}}\Delta t\sum_{i\in{\mathcal{I}}}\Delta
xz^{k}_{i}S^{k}_{i}\right|\leq\|z^{0}\|_{L^{2}({\mathbb{T}})}\|v^{1}\|_{L^{2}({\mathbb{T}})}.$

Now, thanks to (34) we deduce that

$\displaystyle\left|\sum_{k=1}^{N_{T}}\Delta t\sum_{i\in{\mathcal{I}}}\Delta
x(\mu^{k}_{i})^{1/2}z^{k}_{i}\,(\mu^{k}_{i})^{-1/2}S^{k}_{i}\right|\leq
C\left(1+\|\mu\|_{L^{1}(Q_{T})}^{1/2}\right)\|\mu^{-1/2}S\|_{L^{2}(Q_{T})}\|z^{0}\|_{L^{2}({\mathbb{T}})}.$

In the remaining of the proof we want to use the dual definition (10) of the
norm $\|\cdot\|_{L^{2}(Q_{T})}$. Observe that for any vector
$F^{k}=(f^{k}_{i})_{i\in{\mathcal{I}}}$, there exists a unique
$V=(v_{i})_{i\in{\mathcal{I}}}$ such that

$\displaystyle(\mu^{k}_{i})^{-1/2}\frac{v_{i}-v^{k+1}_{i}}{\Delta
t}-(\mu^{k}_{i})^{1/2}\,\left(\Delta_{\mathcal{T}}v\right)_{i}=f^{k}_{i},\quad\forall
i\in{\mathcal{I}},$

where $(v^{k+1}_{i})_{i\in{\mathcal{I}}}$ is a given vector. Indeed, this
system rewrites $(\mathbb{M}^{k})^{\top}V=V^{k+1}+\Delta
t\mathbb{D}^{k}F^{k}$, for all $1\leq k\leq N_{T}$ where $\mathbb{M}^{k}$ is
the invertible tridiagonal matrix given by (26) and
$\mathbb{D}^{k}=\mathrm{diag}((\mu^{k}_{i}))_{i\in{\mathcal{I}}}$. We deduce
thanks to formula (10) that it holds

$\displaystyle\left(\sum_{k=1}^{N_{T}}\Delta t\sum_{i\in{\mathcal{I}}}\Delta
x\mu^{k}_{i}|z^{k}_{i}|^{2}\right)^{1/2}\leq
C\left(1+\|\mu\|_{L^{1}(Q_{T})}^{1/2}\right)\|z^{0}\|_{L^{2}({\mathbb{T}})}.$

This concludes the proof of Theorem 10. ∎

###  4.3. 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)–(15). Let us notice that for every $1\leq k\leq N_{T}$ we have

$\displaystyle\max_{i\in{\mathcal{I}}}\left|\left(\Delta_{\mathcal{T}}\mu^{k}_{1}\right)_{i}\right|$
$\displaystyle=\max_{i\in{\mathcal{I}}}\Bigg{|}d_{11}\sum_{j\in{\mathcal{I}}}\Delta
xu^{k}_{1,j}\left(\Delta_{\mathcal{T}}\sigma_{1}\right)_{i-j}+d_{12}\sum_{j\in{\mathcal{I}}}\Delta
xu^{k}_{2,j}\left(\Delta_{\mathcal{T}}\rho\right)_{i-j}\Bigg{|}$
$\displaystyle\leq
d_{11}\|\Delta_{\mathcal{T}}\sigma_{1}\|_{L^{\infty}({\mathbb{T}})}\sum_{j\in{\mathcal{I}}}\Delta
xu^{k}_{1,j}+d_{12}\|\Delta_{\mathcal{T}}\rho\|_{L^{\infty}({\mathbb{T}})}\sum_{j\in{\mathcal{I}}}\Delta
x\,u^{k}_{2,j}.$

Now, let us recall that $m^{0}_{j}=\|u^{0}_{j}\|_{L^{1}({\mathbb{T}})}$ for
$j=1,2$, then thanks to the mass conservation property (21) we obtain

$\displaystyle\max_{i\in{\mathcal{I}}}\left|\left(\Delta_{\mathcal{T}}\mu^{k}_{1}\right)_{i}\right|$
$\displaystyle\leq
d_{11}m^{0}_{1}\,\|\Delta_{\mathcal{T}}\sigma_{1}\|_{L^{\infty}({\mathbb{T}})}+d_{12}m^{0}_{2}\,\|\Delta_{\mathcal{T}}\rho\|_{L^{\infty}({\mathbb{T}})}.$

Similarly we establish the following bound

$\displaystyle\max_{i\in{\mathcal{I}}}\left|\left(\Delta_{\mathcal{T}}\mu^{k}_{2}\right)_{i}\right|\leq
d_{22}m^{0}_{2}\,\|\Delta_{\mathcal{T}}\sigma_{2}\|_{L^{\infty}({\mathbb{T}})}+d_{21}m^{0}_{1}\,\|\Delta_{\mathcal{T}}\rho\|_{L^{\infty}({\mathbb{T}})}.$

As a direct consequence of the previous estimates and (27) (with
$\tilde{\gamma}=\gamma$ and $\tilde{\Gamma}=\Gamma$) 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)–(15). For every $0\leq k\leq N_{T}$ we define the element
$z^{k}_{i}=u^{k}_{1,i}+u^{k}_{2,i}$ for all $i\in{\mathcal{I}}$. Observe that
$z^{k}_{i}$ is solution to

$\displaystyle\frac{z^{k}_{i}-z^{k-1}_{i}}{\Delta
t}+\Delta_{\mathcal{T}}(\mu^{k}z^{k})_{i}=0,\quad\text{where}\quad\mu^{k}_{i}=\frac{\mu_{1,i}^{k}u_{1,i}^{k}+\mu_{2,i}^{k}u_{2,i}^{k}}{u_{1,i}^{k}+u_{2,i}^{k}},\quad\forall
i\in{\mathcal{I}}.$

Thanks to Theorem 1, we have $u^{k}_{1,i}$, $u^{k}_{2,i}>0$ for all
$i\in{\mathcal{I}}$ and the element $\mu^{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 $\Delta x$ such that

(38) $\sum_{k=1}^{N_{T}}\Delta t\sum_{i\in{\mathcal{I}}}\Delta
x\left(\mu^{k}_{1,i}u^{k}_{1,i}+\mu^{k}_{2,i}u^{k}_{2,i}\right)\left(u^{k}_{1,i}+u^{k}_{2,i}\right)\\\
\leq C\left(1+\sum_{k=1}^{N_{T}}\Delta t\sum_{i\in{\mathcal{I}}}\Delta
x|\mu^{k}_{i}|\right)\left(\sum_{i\in{\mathcal{I}}}\Delta
x|u^{0}_{1,i}|^{2}+\sum_{i\in{\mathcal{I}}}\Delta x|u^{0}_{2,i}|^{2}\right).$

Now we notice that

$\displaystyle\sum_{k=1}^{N_{T}}\Delta t\sum_{i\in{\mathcal{I}}}\Delta
x|\mu^{k}_{i}|\leq\sum_{k=1}^{N_{T}}\Delta t\sum_{i\in{\mathcal{I}}}\Delta
x|\mu^{k}_{1,i}|+\sum_{k=1}^{N_{T}}\Delta t\sum_{i\in{\mathcal{I}}}\Delta
x|\mu^{k}_{2,i}|=I_{9}+I_{10}.$

For $I_{9}$ we have

$\displaystyle I_{9}$ $\displaystyle=\sum_{k=1}^{N_{T}}\Delta
t\sum_{i\in{\mathcal{I}}}\Delta
x\left(d_{1}+d_{11}\sum_{j\in{\mathcal{I}}}\Delta
x\sigma_{1,i-j}u^{k}_{1,j}+d_{12}\sum_{j\in{\mathcal{I}}}\Delta
x\rho_{i-j}u^{k}_{2,j}\right)$ $\displaystyle=\sum_{k=1}^{N_{T}}\Delta
t\left(d_{1}+d_{11}\sum_{j\in{\mathcal{I}}}\Delta
xu^{k}_{1,j}\sum_{i\in{\mathcal{I}}}\Delta
x\sigma_{1,i-j}+d_{12}\sum_{j\in{\mathcal{I}}}\Delta
xu^{k}_{2,j}\sum_{i\in{\mathcal{I}}}\Delta x\rho_{i-j}\right).$

Thus, bearing in mind the mass conservation property (21) we obtain

(39) $\displaystyle I_{9}\leq
T(d_{1}+d_{11}m^{0}_{1}\,\|\sigma_{1}\|_{L^{1}({\mathbb{T}})}+d_{12}m^{0}_{2}\,\|\rho\|_{L^{1}({\mathbb{T}})}),$

and similarly

(40) $\displaystyle I_{10}\leq
T(d_{2}+d_{22}m^{0}_{2}\,\|\sigma_{2}\|_{L^{1}({\mathbb{T}})}+d_{21}m^{0}_{1}\,\|\rho\|_{L^{1}({\mathbb{T}})}).$

Collecting (38)–(40) we conclude that point _(ii)_ of Theorem 2 holds.

##  5\. Convergence of the scheme

This section is dedicated to the proof of Theorem 3. In the following the
subscript $m$ refer to the size $\eta_{m}=\max\\{\Delta x_{m},\Delta t_{m}\\}$
of the family $({\mathcal{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 [4, Lemma 9.2]. This result is the adaptation
at the discrete level of a compactness lemma established by Kruzhkov in [33].
Finally in Section 5.3, we prove Theorem 3.

###  5.1. Uniform estimates

In this section we establish some uniform estimates w.r.t. $\Delta x$ and
$\Delta t$ 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 1.

######  Proposition 11.

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(u^{0}_{1},u^{0}_{2})$ such that

(41)
$\displaystyle\max_{k=1,\ldots,N_{T}}\|u^{k}_{j}\|_{L^{1}({\mathbb{T}})}+\left(d_{j}\sum_{k=1}^{N_{T}}\Delta
t\,\|u^{k}_{j}\|^{2}_{1,1,{\mathcal{T}}}\right)^{\frac{1}{2}}\leq
C_{1},\quad\mbox{for }j=1,2.$

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}$, $\|\sigma_{1}\|_{L^{\infty}({\mathbb{T}})}$,
$\|\sigma_{2}\|_{L^{\infty}({\mathbb{T}})}$,
$\|\rho\|_{L^{\infty}({\mathbb{T}})}$, $m^{0}_{1}$, $m^{0}_{2}$ and
$H(u^{0}_{1},u^{0}_{2})$ such that

(42) $\displaystyle\sum_{k=1}^{N_{T}}\Delta t\sum_{i\in{\mathcal{I}}}\Delta
x\,\left|\mathcal{F}^{k}_{j,i+\frac{1}{2}}\right|\leq C_{2},\quad\mbox{for
}j=1,2,$

where the numerical fluxes are defined by (18).

###### Proof.

The uniform $L^{\infty}(0,T;L^{1}({\mathbb{T}}))$ estimate of the first term
in the right hand side of (41) is a direct consequence of the conservation of
mass (21). Then, for the uniform discrete $L^{2}(0,T;W^{1,1}({\mathbb{T}}))$
estimate, we first notice, for $j=1$ or $2$ and $k\in\\{1,\ldots,N_{T}\\}$,
that it holds

$\displaystyle|u_{j}^{k}|_{1,1,{\mathcal{T}}}=\sum_{i\in{\mathcal{I}}}\left|u^{k}_{j,i+1}-u^{k}_{j,i}\right|=\sum_{i\in{\mathcal{I}}}\left|\left(\sqrt{u^{k}_{j,i+1}}-\sqrt{u^{k}_{j,i}}\right)\,\left(\sqrt{u^{k}_{j,i+1}}+\sqrt{u^{k}_{j,i}}\right)\right|.$

Hence, the Cauchy-Schwarz inequality yields

$\displaystyle|u^{k}_{j}|_{1,1,{\mathcal{T}}}\leq\left|\sqrt{u^{k}_{j}}\right|_{1,2,{\mathcal{T}}}\,\left(\sum_{i\in{\mathcal{I}}}\Delta
x\,\left(\sqrt{u^{k}_{j,i+1}}+\sqrt{u^{k}_{j,i}}\right)^{2}\right)^{\frac{1}{2}}$

Since $(a+b)^{2}\leq 2(a^{2}+b^{2})$ and
$\|u^{k}_{j}\|_{L^{1}({\mathbb{T}})}=\|u^{0}_{j}\|_{L^{1}({\mathbb{T}})}=m^{0}_{j}$
(conservation of mass), one has

$\displaystyle|u^{k}_{j}|_{1,1,{\mathcal{T}}}\leq
2\left(m^{0}_{j}\right)^{1/2}\,\left|\sqrt{u^{k}_{j}}\right|_{1,2,{\mathcal{T}}}.$

Therefore, applying the entropy inequality (22), we get for the first species

$\displaystyle d_{1}\sum_{k=1}^{N_{T}}\Delta
t\,|u^{k}_{1}|^{2}_{1,1,{\mathcal{T}}}\leq
4d_{1}\,m^{0}_{1}\,\sum_{k=1}^{N_{T}}\Delta
t\,\left|\sqrt{u^{k}_{1}}\right|^{2}_{1,2,{\mathcal{T}}}\leq
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 (42). In this purpose we will consider the case $j=1$.
Then, using the definition (18) of the numerical fluxes, we estimate

$\displaystyle\sum_{k=1}^{N_{T}}\Delta t\sum_{i\in{\mathcal{I}}}\Delta
x\,\left|\mathcal{F}^{k}_{1,i+\frac{1}{2}}\right|\leq\sum_{k=1}^{N_{T}}\Delta
t\sum_{i\in{\mathcal{I}}}\mu^{k}_{1,i+\frac{1}{2}}\,\left|u^{k}_{1,i}-u^{k}_{1,i+1}\right|$
$\displaystyle+\sum_{k=1}^{N_{T}}\Delta
t\sum_{i\in{\mathcal{I}}}u^{k}_{1,i+\frac{1}{2}}\,\left|\mu^{k}_{1,i}-\mu^{k}_{1,i+1}\right|$
$\displaystyle=I_{11}+I_{12}.$

For $I_{11}$, applying the regularity of the functions $\sigma_{1}$ and $\rho$
we have

$\displaystyle I_{11}$ $\displaystyle\leq\sum_{k=1}^{N_{T}}\Delta
t\sum_{i\in{\mathcal{I}}}\left(d_{1}+d_{11}\,\|\sigma_{1}\|_{L^{\infty}({\mathbb{T}})}\,\|u^{k}_{1}\|_{L^{1}({\mathbb{T}})}+d_{12}\,\|\rho\|_{L^{\infty}({\mathbb{T}})}\,\|u^{k}_{2}\|_{L^{1}({\mathbb{T}})}\right)\,\left|u^{k}_{1,i}-u^{k}_{1,i+1}\right|.$

Hence, using the conservativity of the scheme and (41), we get

(43) $\displaystyle
I_{11}\leq\left(d_{1}+d_{11}m^{0}_{1}\,\|\sigma_{1}\|_{L^{\infty}({\mathbb{T}})}+d_{12}m^{0}_{2}\,\|\rho\|_{L^{\infty}({\mathbb{T}})}\right)\,T^{1/2}\,C_{1}.$

For $I_{12}$, using the definition of $\mu^{k}_{1,i}$ for $i\in{\mathcal{I}}$,
we notice that it holds

$I_{12}\leq\frac{1}{2}\sum_{k=1}^{N_{T}}\Delta t\sum_{i\in{\mathcal{I}}}\Delta
x\,\left(u^{k}_{1,i}+u^{k}_{1,i+1}\right)\\\
\times\left(d_{11}\sum_{n\in{\mathcal{I}}}\sigma_{1,n}\,\left|u^{k}_{1,i-n}-u^{k}_{1,i+1-n}\right|+d_{12}\sum_{n\in{\mathcal{I}}}\rho_{n}\,\left|u^{k}_{2,i-n}-u^{k}_{2,i+1-n}\right|\right).$

Then, thanks to the conservativity of the scheme, we obtain

$\displaystyle I_{12}$ $\displaystyle\leq 2m^{0}_{1}\sum_{k=1}^{N_{T}}\Delta
t\,\left(d_{11}\,\|\sigma_{1}\|_{L^{\infty}({\mathbb{T}})}\,\left|u^{k}_{1}\right|_{1,1,{\mathcal{T}}}+d_{12}\,\|\rho\|_{L^{\infty}({\mathbb{T}})}\,\left|u^{k}_{2}\right|_{1,1,{\mathcal{T}}}\right).$

Therefore, applying the Cauchy-Schwarz inequality and (41) we end up with

(44) $\displaystyle I_{12}\leq
2\,m^{0}_{1}\,T^{1/2}\,C_{1}\,\left(\frac{d_{11}}{d_{1}^{1/2}}\,\|\sigma_{1}\|_{L^{\infty}({\mathbb{T}})}+\frac{d_{12}}{d_{2}^{1/2}}\,\|\rho\|_{L^{\infty}({\mathbb{T}})}\right).$

Collecting (43) and (44) and the corresponding inequalities for the second
species lead to the existence of $C_{2}$ such that (42) holds. This concludes
the proof of Proposition 11. ∎

###  5.2. Compactness properties

Let $(u_{1,m},u_{2,m})_{m\in{\mathbb{N}}}$ be a family, constructed in Theorem
1, of finite volume solutions to (11)–(15) associated to the sequence
$({\mathcal{D}}_{m})$. In order to be able to apply [4, Lemma 9.2], the first
task is to rewrite the scheme (11)–(15) as the discretization of an evolution
equation under divergence form. In this purpose we use the equivalent form of
(12) given by (17). In particular, for $j=1,2$ and $k=1,\ldots,N_{T}$, we
associate to the family of fluxes
$\left(\mathcal{F}^{k}_{j,i+1/2}\right)_{i\in{\mathcal{I}}}$ the following
piecewise reconstruction

$\displaystyle\mathcal{F}^{k}_{j,m}=\sum_{i\in{\mathcal{I}}_{m}}\mathcal{F}^{k}_{j,i+\frac{1}{2}}\,\mathbf{1}_{(x_{i},x_{i+1})}.$

Then, for this discrete field $\mathcal{F}^{k}_{j,m}$ we define its $L^{1}$
norm as

$\displaystyle\|\mathcal{F}^{k}_{j,m}\|_{L^{1}({\mathbb{T}})}=\sum_{i\in{\mathcal{I}}_{m}}\Delta
x\,\left|\mathcal{F}^{k}_{j,i+\frac{1}{2}}\right|,$

and its discrete divergence by

$\displaystyle\mathrm{div}_{\mathcal{T}}\left(\mathcal{F}^{k}_{j,m}\right)_{i}=\frac{1}{\Delta
x}\,\left(\mathcal{F}^{k}_{j,i+\frac{1}{2}}-\mathcal{F}^{k}_{j,i-\frac{1}{2}}\right),\quad
i\in{\mathcal{I}}_{m}.$

This definition allows us to rewrite (18) as

(45) $\displaystyle\frac{u^{k}_{j,i}-u^{k-1}_{j,i}}{\Delta
t}+\mathrm{div}_{\mathcal{T}}\left(\mathcal{F}^{k}_{j,m}\right)_{i}=0,\quad\forall
i\in{\mathcal{I}}_{m},\,j=1,2,$

and we obtain the following result:

######  Proposition 12.

Let the assumptions of Theorem 3 hold and let
$(u_{1,m},u_{2,m})_{m\in{\mathbb{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})\in\left(L^{p}(Q_{T})\right)^{2}$, with $p\in[1,3)$, such that

$\displaystyle u_{j,m}\rightarrow u_{j}\quad\mbox{strongly in
}L^{p}(Q_{T})\mbox{ for }1\leq p<3,\quad\mbox{as }m\to\infty,$

and almost everywhere.

###### Proof.

A direct consequence of Proposition 11 is that there is a constant $C$
independent of $\Delta x_{m}$ and $\Delta t_{m}$ such that

$\displaystyle\sum_{k=1}^{N_{T}}\Delta
t_{m}\,\|u^{k}_{j,m}\|_{L^{1}({\mathbb{T}})}+\sum_{k=1}^{N_{T}}\Delta
t_{m}\,\|\mathcal{F}^{k}_{j,m}\|_{L^{1}({\mathbb{T}})}+\sum_{k=1}^{N_{T}}\Delta
t_{m}\,|u^{k}_{j,m}|_{1,1,{\mathcal{T}}_{m}}\leq C,\quad j=1,2.$

By [4, Lemma 9.2], which can be applied thanks to (45), there is function
$u_{j}\in L^{1}(Q_{T})$, $j=1,2$, such that, up to a subsequence,

$\displaystyle u_{j,m}\rightarrow u_{j}\quad\mbox{strongly in
}L^{1}(Q_{T})\mbox{ as }m\to\infty.$

Moreover, Proposition 11 also implies that the sequence $(u_{i,m})$ is
uniformly bounded in the space $L^{\infty}(0,T;L^{1}({\mathbb{T}}))$ and in
$L^{2}(0,T;BV({\mathbb{T}}))$. The continuous embedding of $BV({\mathbb{T}})$
in $L^{\infty}({\mathbb{T}})$ (see [1]) implies that the sequence $(u_{j,m})$
is uniformly bounded in $L^{2}(0,T;L^{\infty}({\mathbb{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\in[1,3)$. This concludes the proof of Proposition
12. ∎

######  Remark 13.

In dimension $d\geq 2$, we have the compact embedding of the space
$BV({\mathbb{T}}^{d})$ in $L^{\frac{d}{d-1}}({\mathbb{T}}^{d})$. In particular
in this case the sequence $(u_{j,m})$ is uniformly bounded in
$L^{\frac{d}{d-1}}(Q_{T})$. Therefore arguing as in the previous proof we
deduce the existence for $j=1$ and $2$ of $u_{j}\in L^{p}(Q_{T})$ for
$p\in[1,{d}/{(d-1)})$, such that, up to a subsequence,

$\displaystyle u_{j,m}\rightarrow u_{j}\quad\mbox{strongly in
}L^{p}(Q_{T})\mbox{ for }1\leq p<{\frac{d}{d-1}},\quad\mbox{as }m\to\infty.$

######  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\in[1,3)$ one has

$\displaystyle\mu_{1,m}\rightarrow\mu_{1}=d_{1}+d_{12}\rho_{1}\ast
u_{2}+d_{11}\sigma_{1}\ast u_{1}\quad\mbox{strongly in
}L^{p}(0,T;L^{\infty}({\mathbb{T}})),\quad\mbox{as }m\to\infty,$
$\displaystyle\mu_{2,m}\rightarrow\mu_{2}=d_{2}+d_{21}\rho_{2}\ast
u_{1}+d_{22}\sigma_{2}\ast u_{2}\quad\mbox{strongly in
}L^{p}(0,T;L^{\infty}({\mathbb{T}})),\quad\mbox{as }m\to\infty,$

where we recall that $\rho_{1}(x)=\rho_{2}(-x)=\rho(x)$ for a.e.
$x\in{\mathbb{T}}$.

###### Proof.

Observe that by definition (14)-(16) one has for $x\in K_{i}$ and $t\in(0,T)$

$\displaystyle\mu_{1}(x,t)-\mu_{1,m}(x,t)$
$\displaystyle=d_{12}\,(\rho_{1}\ast(u_{2}-u_{2,m}))(x,t)+d_{11}\,(\sigma_{1}\ast(u_{1}-u_{1,m}))(x,t)$
$\displaystyle+d_{12}\int_{{\mathbb{T}}}(\rho_{1}(x-y)-\rho_{1}(x_{i}-y))\,u_{2,m}(y,t)\,\mathrm{d}y$
$\displaystyle+d_{11}\int_{{\mathbb{T}}}(\sigma_{1}(x-y)-\sigma_{1}(x_{i}-y))\,u_{1,m}(y,t)\,\mathrm{d}y.$

Therefore by dominated convergence one obtains the strong convergence of
$\mu_{1,m}$ towards $\mu_{1}$ in $L^{1}(Q_{T})$ and almost everywhere in
$Q_{T}$. Besides, for a.e. $t\in(0,T)$ thanks to Young’s inequality and for
$p\in(1,3)$, it holds

$\displaystyle\|\rho_{1}\ast(u_{2}-u_{2,m})(t)\|_{L^{\infty}({\mathbb{T}})}\leq\|\rho_{1}\|_{L^{\frac{p}{p-1}}({\mathbb{T}})}\,\|(u_{2}-u_{2,m})(t)\|_{L^{p}({\mathbb{T}})}.$

Then, applying Proposition 12, we obtain

$\displaystyle\int_{0}^{T}\|\rho_{1}\ast(u_{2}-u_{2,m})(t)\|^{p}_{L^{\infty}({\mathbb{T}})}\mathrm{d}t\leq\|\rho_{1}\|^{p}_{L^{\frac{p}{p-1}}({\mathbb{T}})}\,\int_{0}^{T}\|(u_{2}-u_{2,m})(t)\|^{p}_{L^{p}({\mathbb{T}})}\,\mathrm{d}t\to
0\quad\mbox{as }m\to\infty.$

Let us now setting $\xi(x,y)=\rho_{1}(x-y)-\rho_{1}(x_{i}-y)$ for a.e. $x\in
K_{i}$ and $y\in{\mathbb{T}}$. Hence, for a.e. $t\in(0,T)$, the Hölder
inequality yields for $p\in(1,3)$

$\left\|\int_{{\mathbb{T}}}\xi(\cdot,y)u_{2,m}(t,y)\mathrm{d}y\right\|_{L^{\infty}({\mathbb{T}})}\leq\sup_{|z|\leq\Delta
x_{m}}\|\rho_{1}(z+\cdot)-\rho_{1}\|_{L^{\frac{p}{p-1}}({\mathbb{T}})}\|u_{2,m}(t)\|_{L^{p}({\mathbb{T}})}.$

The first factor in the right hand side tends to $0$ (by density of continuous
functions in $L^{p/(p-1)}({\mathbb{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^{\infty}({\mathbb{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 (8). Let
$\phi\in C^{\infty}_{0}({\mathbb{T}}\times[0,T))$, we multiply equation (12)
by $\Delta t_{m}\Delta x_{m}\,\phi^{k-1}_{i}$, where
$\phi^{k-1}_{i}=\phi(x_{i},t^{k-1})$, and we sum over $i\in{\mathcal{I}}$ and
$k\in\\{1,\ldots,N_{T}\\}$. We obtain $F^{m}_{1}+F^{m}_{2}=0$ with

$\displaystyle F^{m}_{1}$
$\displaystyle=\sum_{k=1}^{N_{T}}\sum_{i\in{\mathcal{I}}}\Delta
x_{m}\,(u^{k}_{1,i}-u^{k-1}_{1,i})\phi^{k-1}_{i}=\sum_{k=1}^{N_{T}}\sum_{i\in{\mathcal{I}}}\Delta
x_{m}\,u^{k}_{1,i}(\phi^{k-1}_{i}-\phi^{k}_{i})-\sum_{i\in{\mathcal{I}}}\Delta
x_{m}u^{0}_{1,i}\phi^{0}_{i},$

and

$\displaystyle F^{m}_{2}$ $\displaystyle=\sum_{k=1}^{N_{T}}\Delta
t_{m}\sum_{i\in{\mathcal{I}}}\frac{-\mu^{k}_{1,i+1}u^{k}_{1,i+1}+2\mu^{k}_{1,i}u^{k}_{1,i}-\mu^{k}_{1,i-1}u^{k}_{1,i-1}}{\Delta
x_{m}}\,\phi^{k-1}_{i}$ $\displaystyle=\sum_{k=1}^{N_{T}}\Delta
t_{m}\sum_{i\in{\mathcal{I}}}\frac{-\phi^{k-1}_{i+1}+2\phi^{k-1}_{i}-\phi^{k-1}_{i-1}}{\Delta
x_{m}}\,\mu^{k}_{1,i}u^{k}_{1,i}.$

Let $\psi_{m}(x,t)=(\phi^{k}_{i}-\phi^{k-1}_{i})/\Delta t_{m}$ for all $x\in
K_{i}$ and $t\in(t_{k-1},t_{k}]$ and $\phi_{m}^{0}(x)=\phi^{0}_{i}$ for all
$x\in K_{i}$. Then, since $\psi_{m}\to\partial_{t}\phi$ in $L^{\infty}(Q_{T})$
and $\phi^{0}_{m}\to\phi^{0}$ in $L^{\infty}({\mathbb{T}})$, so using the
convergence results of Proposition 12 one obtains

$\displaystyle
F^{m}_{1}+\int_{Q_{T}}u_{1}\partial_{t}\phi\,\mathrm{d}x\mathrm{d}t+\int_{{\mathbb{T}}}u^{0}_{1}(x)\phi(x,0)\,\mathrm{d}x$
$\displaystyle=$ $\displaystyle\int_{Q_{T}}(u_{1}\partial_{t}\phi-
u_{1,m}\psi_{m})\,\mathrm{d}x\mathrm{d}t+\int_{{\mathbb{T}}}(u^{0}_{1}\phi(\cdot,0)-u^{0}_{1}\phi_{m}^{0})\,\mathrm{d}x\to
0\quad\mbox{as }m\to 0.$

Similarly, if one defines
$\zeta_{m}(x,t)=(\phi^{k-1}_{i+1}-2\phi^{k-1}_{i}+\phi^{k-1}_{i-1})/(\Delta
x_{m})^{2}$ for all $x\in K_{i}$ and $t\in(t_{k-1},t_{k}]$, then
$\zeta_{m}\to\Delta\phi$ in $L^{\infty}(Q_{T})$ and therefore using
Proposition 12 and Corollary 14 one obtains

$\displaystyle
F^{m}_{2}+\int_{Q_{T}}u_{1}\mu_{1}\Delta\phi\,\mathrm{d}x\mathrm{d}t=\int_{Q_{T}}(u_{1}\mu_{1}\Delta\phi-
u_{1,m}\mu_{1,m}\zeta_{m})\,\mathrm{d}x\mathrm{d}t\to 0\quad\mbox{as }m\to 0.$

This concludes the proof of Theorem 3.

##  6\. 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 $\ell^{\infty}$ 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, $\Delta t$ is divided by $2$. If there was refinement on a
given time step, $\Delta 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

$\displaystyle\partial_{t}u_{1}-\partial_{xx}^{2}((\rho\ast u_{2})u_{1})$
$\displaystyle=0,$
$\displaystyle\partial_{t}u_{2}-2\partial_{xx}^{2}((\rho\ast u_{1})u_{2})$
$\displaystyle=0.$

The convolution kernel is taken to be either the Dirac measure, which we
denote by $\rho_{0}$, either by an approximation of a Dirac

(46) $\rho_{\delta}(x)=\delta^{-1}\chi_{[-\delta/2,\delta/2]}(x),$

where $\chi_{A}$ indicator function of the set $A$, or the smooth kernel

$\rho_{\text{smooth}}(x)=\cos(\nu_{L}x)+1.$

with $\nu_{L}=2\pi/L$. We consider two initial data, either the indicator
functions

(47) $u^{0}_{1}(x)=\chi_{[\frac{L}{9},\frac{L}{3}]}(x)\,,\quad
u^{0}_{2}(x)=\chi_{[\frac{L}{3},\frac{3L}{4}]}(x)\,,\quad
x\in\mathbb{R}/L\mathbb{Z}.$

or the smooth functions

(48) $u^{0}_{1}(x)=\cos\left(\nu_{L}x\right)+1\,,\quad
u^{0}_{2}(x)=\sin\left(\nu_{L}x\right)+1\,,\quad x\in\mathbb{R}/L\mathbb{Z}.$

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\cdot 2^{k-1}$ for
$k=1,\dots,6$ and the corresponding time step $\Delta t_{k}=\Delta t_{0}\cdot
4^{-(k-1)}$, with $\Delta 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
$\ell^{\infty}$ 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.

Convolution kernel: | Smooth $\rho_{\text{smooth}}$ | Indicator $\rho_{L/4}$ | Dirac $\rho_{0}$
---|---|---|---
Initial condition: | order: $1.97$ | order: $1.53$ | order: $1.04$
indicator func. (47) | error: $5\cdot 10^{-3}$ | error: $7\cdot 10^{-2}$ | error: $2.7\cdot 10^{-3}$
Initial condition: | order: $2.32$ | order: $2.02$ | order: $2.32$
smooth func. (48) | error: $4.9\cdot 10^{-4}$ | error: $9.2\cdot 10^{-4}$ | error: $5\cdot 10^{-4}$
Table 1. Estimated order of convergence in space and absolute error at final
time for the mesh $N=512$ in $L^{\infty}$ 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 [21, 37], the proofs rely on some
compactness method and no explicit “error” bounds are available (see also for
instance [31]). 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 $\rho=\rho_{\delta}$ for various values of $\delta\in[0,L]$.
The domain has length $L=25$, the final time is $T=1$ and the mesh is such
that $N=1024$ and $\Delta t=10^{-2}$. We evaluate the error at time $T$
between the solution $(u_{1}^{(\delta)},u_{2}^{(\delta)})$ computed for the
kernel $\rho_{\delta}$ and $(u_{1}^{(0)},u_{2}^{(0)})$ computed for the local
cross-diffusion system in Wasserstein-1 norm,

$W_{1}(u_{1}^{(\delta)},u_{1}^{(0)})+W_{1}(u_{2}^{(\delta)},u_{2}^{(0)}),$

and in $L^{p}$ norms

$\|u_{1}^{(\delta)}-u_{1}^{(0)}\|_{L^{p}}+\|u_{2}^{(\delta)}-u_{2}^{(0)}\|_{L^{p}},$

with $p=1$ or $p=\infty$. For the computation of the Wasserstein-1 norm we
recall that in dimension $1$, if $f$ and $g$ are non-negative integrable
functions on $\mathbb{R}$ with the same mass, one has
$W_{1}(f,g)=\|F-G\|_{L^{1}}$, with $F^{\prime}=f$, $G^{\prime}=g$ and
$F(-\infty)=G(-\infty)$. 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
$\delta/L=4W_{1}(\rho_{\delta},\rho_{0})/L$ for the two initial data (47) and
(48).

$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$$10^{-2}$$10^{-1}$$10^{0}$$10^{1}$$\delta/L$$W_{1}$ (order: $0.65$)$L^{1}$ (order: $0.38$)$L^{\infty}$ |  |  $10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$$10^{-6}$$10^{-5}$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$$10^{1}$$\delta/L$$W_{1}$ (order: $2.01$)$L^{1}$ (order: $1.98$)$L^{\infty}$ (order: $1.97$)
---|---|---

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

.

For the smooth initial data (48) supported on the whole domain (up to one
point), we observe convergence with rate
$O(W_{1}(\rho_{\delta},\rho_{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^{\infty}$ norm, and
$O(W_{1}(\rho_{\delta},\rho_{0})^{\alpha})$ convergence with $\alpha\approx
0.38$ in $L^{1}$ norm and $\alpha=0.65$ in Wasserstein-1 norm.

###  6.4. Test case 3: Turing instabilities in prey-predator systems with
nonlocal cross-diffusion

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

$\displaystyle\partial_{t}u_{1}-d_{1}\Delta u_{1}$
$\displaystyle=R_{1}(u_{1},u_{2}),$
$\displaystyle\partial_{t}u_{2}-\Delta((d_{2}+d_{21}\,\rho_{2}\ast
u_{1})u_{2})$ $\displaystyle=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 $\rho_{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 [35] 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 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

$\displaystyle R_{1}(u_{1},u_{2})=au_{1}+eu_{1}^{2}-bu_{1}u_{2},\quad
R_{2}(u_{1},u_{2})=-du_{2}^{2}+cu_{1}u_{2}.$

where the parameters are $a=b=c=d=1$ and $e=\frac{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

$\rho_{2}(x)=C_{r}[x^{2}\chi_{(-r,r)}(x)+(x-2r)^{2}\chi_{[r,2r)}(x)+(x+2r)^{2}\chi_{(-2r,-r]}(x)],$

with $C_{r}$ a normalizing constant such that $\int\rho_{2}=1$. This kernel
vanishes at $0$, has support on $[-2r,2r]$ and is maximal at $x=\pm 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 $\Delta t=0.1$. The initial data is taken as a small perturbation of the
homogeneous equilibrium

$u^{0}_{1}(x)=\frac{ad}{bc-de}+\varepsilon\chi_{[L/3,L/9]}(x)\,,\quad
u^{0}_{2}(x)=\frac{ac}{bc-de},$

with $\varepsilon=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.

$0$$2$$4$$6$$8$$10$$12$$14$$16$$18$$20$$22$$24$$0$$0.5$$1$$1.5$$2$$2.5$$3$$3.5$$4$$4.5$$5$$x$Preys $u_{1}$Predators $u_{2}$ |  $0$$2$$4$$6$$8$$10$$12$$14$$16$$18$$20$$22$$24$$0$$0.5$$1$$1.5$$2$$2.5$$3$$3.5$$4$$4.5$$5$$x$
---|---
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$).

####  6.4.2. Two dimensional case: Mimura-Nishiura-Yamaguti reaction term

Now we consider the two-dimensional case with the following reaction terms

$\displaystyle
R_{1}(u_{1},u_{2})=au_{1}+eu_{1}^{2}-du_{1}^{3}-bu_{1}u_{2},\quad
R_{2}(u_{1},u_{2})=-fu_{2}-gu_{2}^{2}+cu_{1}u_{2},$

with $a=35/9$, $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

$\displaystyle[f+(c(b^{2}c^{2}-2bceg+4dfbg+e^{2}g^{2}+4adg^{2})^{1/2}-bc^{2}+ceg-2dfg)/(2dg)]/c=5,$

and by

$\displaystyle[c(b^{2}c^{2}-2bceg+4dfbg+e^{2}g^{2}+4adg^{2})^{1/2}-bc^{2}+ceg-2dfg]/(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

$\rho_{2}^{\text{sym}}(x,y)=C\chi_{(3/8,1/2)}(x^{2}+y^{2}),$

with $C$ a normalizing constant such that $\int\rho_{2}^{\text{sym}}=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

$\rho_{2}^{\text{nonsym}}(x,y)=C\chi_{(3/8,1/2)}(x^{2}+y^{2})\chi_{[0,\infty)}(x)\chi_{[0,\infty)}(y),$

with $C$ a normalizing constant such that $\int\rho_{2}^{\text{nonsym}}=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 $\Delta 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\times 100$ cells. The initial data is taken as a small
perturbation of the homogeneous equilibrium

$u^{0}_{1}(x,y)=5+\varepsilon\chi_{[L_{x}/9,4L_{x}/9]\times[7L_{y}/9,8L_{y}/9]}(x,y)\,,\quad
u^{0}_{2}(x,y)=10,$

and $\varepsilon=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
$\rho_{2}^{\text{nonsym}}$.

|  |
---|---
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$,
$\rho_{2}=\rho_{2}^{\text{sym}}$); (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$, $\rho_{2}=\rho_{2}^{\text{nonsym}}$).

### 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

For the sake of completness we prove in this section a discrete Grönwall
inequality

######  Lemma 15.

Let $(u_{n})_{n\in{\mathbb{N}}}$ and $(a_{n})_{n\in{\mathbb{N}}}$ be some non-
negative sequences with

$\displaystyle\Delta t<1/\sup_{k\in{\mathbb{N}}}a_{k},$

and satisfying the inequality

(49) $\displaystyle u_{n}\leq u_{0}+\Delta
t\sum_{k=1}^{n}a_{k}u_{k},\quad\forall n\geq 1.$

Then, the following relation holds

$\displaystyle u_{n}\leq u_{0}\Pi_{k=1}^{n}(1-\Delta ta_{k})^{-1},\quad\forall
n\geq 1.$

###### Proof.

We define $w_{0}=u_{0}$ and for $n\geq 1$

$\displaystyle w_{n}=\frac{u_{0}+\Delta
t\sum_{k=1}^{n}a_{k}u_{k}}{\Pi_{k=1}^{n}(1-\Delta ta_{k})^{-1}}.$

Then, for every $n\geq 1$, a simple computation leads to

$\displaystyle w_{n}-w_{n-1}=-\Delta ta_{n}\Pi_{k=1}^{n-1}(1-\Delta
ta_{k})\,\left(u_{0}-u_{n}+\Delta t\sum_{k=1}^{n}a_{k}u_{k}\right).$

Thanks to (49) we deduce that the sequence $(w_{n})_{n\in{\mathbb{N}}}$ is
decreasing and since $w_{0}=u_{0}$ we end up with

$\displaystyle u_{0}+\Delta t\sum_{k=1}^{n}a_{k}u_{k}\leq
u_{0}\Pi_{k=1}^{n}(1-\Delta ta_{k})^{-1},$

which finishes the proof of Lemma 15. ∎

##  Appendix B Discrete entropy inequality on a bounded domain

Following [21], let us assume that $\Omega=(0,1)$ and $d_{11}=d_{22}=0$. then
we consider, for all $x,\,y\in\Omega$ and $t\geq 0$, the system

(50)
$\displaystyle\partial_{t}u_{1}(x)-d_{1}\partial_{x}^{2}u_{1}(x)-\partial_{x}\left(\int_{y\in\Omega}G(x,y)\Big{(}d_{12}u_{2}(y)\partial_{x}u_{1}(x)+d_{12}u_{1}(x)\partial_{y}u_{2}(y)\Big{)}\mathrm{d}y\right)=0,$
(51)
$\displaystyle\partial_{t}u_{2}(y)-d_{2}\partial_{y}^{2}u_{2}(y)-\partial_{y}\left(\int_{x\in\Omega}G(x,y)\Big{(}d_{21}u_{2}(y)\partial_{x}u_{1}(x)+d_{21}u_{1}(x)\partial_{y}u_{2}(y)\Big{)}\mathrm{d}x\right)=0,$

where $G$ is a $C^{0}$ nonnegative kernel such that

(52) $\displaystyle G(x,\cdot)=0\mbox{ if }x=0,\,1\mbox{ or
}G(\cdot,y)=0\mbox{ if }y=0,\,1.$

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

(53)
$\displaystyle\partial_{x}u_{j}(0,t)=\partial_{x}u_{j}(1,t)=0\quad\mbox{for
}j=1,2,$

and initial conditions

(54) $\displaystyle u_{1}(\cdot,0)=u_{1}^{0}(\cdot),\quad
u_{2}(\cdot,0)=u_{2}^{0}(\cdot),\quad\mbox{on }\Omega.$

In this case we define the entropy functional $H$ as

$\displaystyle
H(u_{1},u_{2})=\int_{\Omega}h_{1}(u_{1})\,\mathrm{d}x+\int_{\Omega}h_{2}(u_{2})\,\mathrm{d}x,$

where we recall definition (19) of $h_{1}$ and $h_{2}$. Then, as proved in
[21, Proposition 7], the following entropy inequality holds

(55)
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}H(u_{1},u_{2})+4\frac{d_{1}}{d_{12}}\int_{\Omega}\left|\partial_{x}\sqrt{u_{1}}\right|^{2}+4\frac{d_{2}}{d_{21}}\int_{\Omega}\left|\partial_{x}\sqrt{u_{2}}\right|^{2}\leq
0.$

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
$\Omega=(0,1)$. An uniform mesh, consists in a finite sequence of cells
denoted by $K_{i}=\left(x_{i-1/2},x_{i+1/2}\right)$, for $1\leq i\leq\ell$,
where the points $x_{i+1/2}$ are uniformly distributed such that

$\displaystyle
0=x_{\frac{1}{2}}<x_{\frac{3}{2}}<...<x_{\ell-\frac{1}{2}}<x_{\ell+\frac{1}{2}}=1.$

We note $\Delta x=x_{i+\frac{1}{2}}-x_{i-\frac{1}{2}},$ for some $1\leq
i\leq\ell$, the length of the mesh. Moreover, for $1\leq i\leq\ell$, we define
$x_{i}$ as the center of the cell $K_{i}$ and $x_{0}=x_{\frac{1}{2}}$ and
$x_{\ell+1}=x_{\ell+\frac{1}{2}}$. Then, for $T>0$ given, we define an integer
$N_{T}$ and a time step $\Delta t=T/N_{T}$ and we introduce the sequence
$(t_{k})_{0\leq k\leq N_{T}}$ with $t_{k}=k\Delta t$.

Now, as in Section 2.2, we first discretize the initial conditions (54) as

(56) $\displaystyle u^{0}_{j,i}=\frac{1}{\Delta
x}\int_{K_{i}}u_{j}^{0}(x)\,\mathrm{d}x,\quad\forall 1\leq i\leq\ell,\,j=1,2.$

Now for $(u^{k-1}_{1},u^{k-1}_{2})\in{\mathbb{R}}^{2\ell}$ given, the scheme
writes as

(57) $\displaystyle\Delta x\frac{u^{k}_{1,i}-u^{k-1}_{1,i}}{\Delta
t}+\mathcal{F}^{k}_{1,i+\frac{1}{2}}-\mathcal{F}^{k}_{1,i-\frac{1}{2}}=0,\quad\forall
1\leq i\leq\ell,$ (58) $\displaystyle\Delta
x\frac{u^{k}_{2,n}-u^{k-1}_{2,n}}{\Delta
t}+\mathcal{F}^{k}_{2,n+\frac{1}{2}}-\mathcal{F}^{k}_{2,n-\frac{1}{2}}=0,\quad\forall
1\leq n\leq\ell,$

where the numerical fluxes are defined by

(59) $\displaystyle\mathcal{F}_{1,i+\frac{1}{2}}$
$\displaystyle=d_{1}\frac{u^{k}_{1,i}-u^{k}_{1,i+1}}{\Delta
x}+d_{12}\,\mu^{k}_{2,i+\frac{1}{2}}\frac{u^{k}_{1,i}-u^{k}_{1,i+1}}{\Delta
x}+d_{12}u^{k}_{1,i+\frac{1}{2}}\,\tilde{\mu}^{k}_{2,i+\frac{1}{2}},\quad\forall
1\leq i\leq\ell,$ (60) $\displaystyle\mathcal{F}_{2,n+\frac{1}{2}}$
$\displaystyle=d_{2}\frac{u^{k}_{2,n}-u^{k}_{2,n+1}}{\Delta
x}+d_{21}\,\mu^{k}_{1,n+\frac{1}{2}}\frac{u^{k}_{2,n}-u^{k}_{2,n+1}}{\Delta
x}+d_{21}u^{k}_{2,n+\frac{1}{2}}\,\tilde{\mu}^{k}_{1,n+\frac{1}{2}},\quad\forall
1\leq n\leq\ell,$

with $u^{k}_{1,i+1/2}=(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

(61)
$\displaystyle\mathcal{F}^{k}_{j,\frac{1}{2}}=\mathcal{F}^{k}_{j,\ell+\frac{1}{2}}=0,\quad\forall
1\leq k\leq N_{T},\,j=1,2.$

Finally, for every $1\leq i,n\leq\ell$ we define

(62) $\displaystyle\mu^{k}_{2,i+\frac{1}{2}}$
$\displaystyle=\sum_{n=0}^{\ell-1}\Delta
xG\left(x_{i+\frac{1}{2}},x_{n+\frac{1}{2}}\right)\,\,u^{k}_{2,n+\frac{1}{2}},\quad\mu^{k}_{1,n+\frac{1}{2}}=\sum_{i=0}^{\ell-1}\Delta
xG\left(x_{i+\frac{1}{2}},x_{n+\frac{1}{2}}\right)\,\,u^{k}_{1,i+\frac{1}{2}},$
(63) $\displaystyle\tilde{\mu}^{k}_{2,i+\frac{1}{2}}$
$\displaystyle=\sum_{n=0}^{\ell-1}G\left(x_{i+\frac{1}{2}},x_{n+\frac{1}{2}}\right)\,(u^{k}_{2,n}-u^{k}_{2,n+1}),\quad\tilde{\mu}^{k}_{1,n+\frac{1}{2}}=\sum_{i=0}^{\ell-1}G\left(x_{i+\frac{1}{2}},x_{n+\frac{1}{2}}\right)\,(u^{k}_{1,i}-u^{k}_{1,i+1}).$

Now let us prove that, for a given $0<k\leq N_{T}$, the solutions to scheme
(54)–(63) satisfy (a priori) the inequality

(64)
$\displaystyle\frac{H(u^{k}_{1},u^{k}_{2})-H(u^{k-1}_{1},u^{k-1}_{2})}{\Delta
t}+4\frac{d_{1}}{d_{12}}\sum_{i=0}^{\ell-1}\frac{\left(\sqrt{u^{k}_{1,i}}-\sqrt{u^{k}_{1,i+1}}\right)}{\Delta
x}+4\frac{d_{2}}{d_{21}}\sum_{n=0}^{\ell-1}\frac{\left(\sqrt{u^{k}_{2,n}}-\sqrt{u^{k}_{2,n+1}}\right)}{\Delta
x}\leq 0.$

In this aim let us assume that $u^{k}_{j,i}>0$ for all $1\leq i\leq\ell$ 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\in\\{1,\ldots,\ell\\}$ respectively, we apply some discrete integration
by parts, we use definitions (62) and (63) and we sum the two relations to
obtain

$\displaystyle Q_{1}+Q_{2}+Q_{3}+Q_{4}\leq 0,$

with

$\displaystyle Q_{1}$
$\displaystyle=\frac{H(u^{k}_{1},u^{k}_{2})-H(u^{k-1}_{1},u^{k-1}_{2})}{\Delta
t},$ $\displaystyle Q_{2}$
$\displaystyle=\frac{d_{1}}{d_{12}}\sum_{i=0}^{\ell-1}\frac{(u^{k}_{1,i}-u^{k}_{1,i+1})(\log(u^{k}_{1,i})-\log(u^{k}_{1,i+1}))}{\Delta
x}\geq
4\frac{d_{1}}{d_{12}}\sum_{i=0}^{\ell-1}\frac{\left(\sqrt{u^{k}_{1,i}}-\sqrt{u^{k}_{1,i+1}}\right)}{\Delta
x},$ $\displaystyle Q_{3}$
$\displaystyle=\frac{d_{2}}{d_{21}}\sum_{n=0}^{\ell-1}\frac{(u^{k}_{2,n}-u^{k}_{2,n+1})(\log(u^{k}_{2,n})-\log(u^{k}_{2,n+1}))}{\Delta
x}\geq
4\frac{d_{2}}{d_{21}}\sum_{n=0}^{\ell-1}\frac{\left(\sqrt{u^{k}_{2,n}}-\sqrt{u^{k}_{2,n+1}}\right)}{\Delta
x},$ $\displaystyle Q_{4}$ $\displaystyle=\frac{1}{\Delta
x}\sum_{i=0}^{\ell-1}\sum_{n=0}^{\ell-1}\Delta
xG\left(x_{i+\frac{1}{2}},x_{n+\frac{1}{2}}\right)\left(\begin{matrix}u^{k}_{1,i}-u^{k}_{1,i+1}\\\
u^{k}_{2,n}-u^{k}_{2,n+1}\end{matrix}\right)^{\top}\,{\mathbb{A}}^{k}_{i,n}\,\left(\begin{matrix}u^{k}_{1,i}-u^{k}_{1,i+1}\\\
u^{k}_{2,n}-u^{k}_{2,n+1}\end{matrix}\right),$

where we have used for $Q_{1}$ the convexity of $h_{1}$ and $h_{2}$, the
inequality $(a-b)(\log(a)-\log(b))\geq 4(\sqrt{a}-\sqrt{b})^{2}$ for $Q_{2}$
and $Q_{3}$ and where the matrix ${\mathbb{A}}^{k}_{i,n}$ is defined by

$\displaystyle{\mathbb{A}}^{k}_{i,n}=\left(\begin{matrix}\dfrac{u^{k}_{2,n+\frac{1}{2}}}{u^{k,\log}_{1,i+\frac{1}{2}}}&&\dfrac{u^{k}_{1,i+\frac{1}{2}}}{u^{k,\log}_{1,i+\frac{1}{2}}}\\\
\dfrac{u^{k}_{2,n+\frac{1}{2}}}{u^{k,\log}_{2,n+\frac{1}{2}}}&&\dfrac{u^{k}_{1,i+\frac{1}{2}}}{u^{k,\log}_{2,n+\frac{1}{2}}}\end{matrix}\right)\quad\forall
0\leq i,n\leq\ell-1,$

with for $j=1,2$ and every $1\leq i\leq\ell$

$u^{k,\log}_{j,i+\frac{1}{2}}=\begin{cases}\dfrac{u^{k}_{j,i}-u^{k}_{j,i+1}}{\log(u^{k}_{j,i})-\log(u^{k}_{j,i+1})}\quad&\mbox{if
}u^{k}_{j,i}>0,\ u^{k}_{j,i+1}>0,\mbox{ and }u^{k}_{j,i}\neq u^{k}_{j,i+1},\\\
u^{k}_{j,i}&\mbox{if }u^{k}_{j,i}=u^{k}_{j,i+1}>0.\end{cases}$

Now, thanks to the nonnegativity of the kernel $G$, it remains to show that
${\mathbb{A}}^{k}_{i,n}$ is a positive semidefinite matrix in order to
conclude that $Q_{4}\geq 0$. In this purpose, we notice that
$\det({\mathbb{A}}^{k}_{i,n})=0$ and
$\mathrm{trace}({\mathbb{A}}^{k}_{i,n})>0$ this imply that
${\mathbb{A}}^{k}_{i,n}$ has two nonnegative eigenvalues. Therefore we
conclude that for all $0\leq i,n\leq\ell-1$ the matrix
${\mathbb{A}}^{k}_{i,n}$ is a positive semidefinite matrix and $Q_{4}\geq 0$.
This implies that the entropy inequality (64) holds which proves the claim.

## References

  * [1] L. Ambrosio, N. Fusco, and D. Pallara.  Functions of bounded variation and free discontinuity problems.  Oxford: Clarendon Press, 2000.
  * [2] V. Anaya, M. Bendahmane, M. Langlais, and M. Sepúlveda.  A convergent finite volume method for a model of indirectly transmitted diseases with nonlocal cross-diffusion.  Comput. Math. Appl., 70(2):132–157, 2015.
  * [3] V. Anaya, M. Bendahmane, and M. Sepúlveda.  Numerical analysis for a three interacting species model with nonlocal and cross diffusion.  ESAIM, Math. Model. Numer. Anal., 49(1):171–192, 2015.
  * [4] B. Andreianov, M. Bendahmane, and R. Ruiz-Baier.  Analysis of a finite volume method for a cross-diffusion model in population dynamics.  Math. Models Methods Appl. Sci., 21(2):307–344, 2011.
  * [5] J. Baladron, D. Fasoli, O. Faugeras, and J. Touboul.  Mean-field description and propagation of chaos in networks of Hodgkin-Huxley and FitzHugh-Nagumo neurons.  J. Math. Neurosci., 2:50, 2012.  Id/No 10.
  * [6] M. Bessemoulin-Chatard, M. Herda, and T. Rey.  Hypocoercivity and diffusion limit of a finite volume scheme for linear kinetic equations.  Math. Comput., 89(323):1093–1133, 2020.
  * [7] M. Braukhoff, I. Perugia, and P. Stocker.  An entropy structure preserving space-time formulation for cross-diffusion systems: Analysis and galerkin discretization.  SIAM J. Num. Anal., 60(1):364–395, 2022.
  * [8] M. Burger, J.A. Carrillo, J.-F. Pietschmann, and M. Schmidtchen.  Segregation effects and gap formation in cross-diffusion models.  Interfaces Free Bound., 22(2):175–203, 2020.
  * [9] M. Burger, M. Di Francesco, J.-F. Pietschmann, and B. Schlake.  Nonlinear cross-diffusion with size exclusion.  SIAM J. Math. Anal., 42(6):2842–2871, 2010.
  * [10] C. Cancès, C. Chainais-Hillairet, A. Gerstenmayer, and A. Jüngel.  Finite-volume scheme for a degenerate cross-diffusion model motivated from ion transport.  Numer. Methods Partial Differ. Equations, 35(2):545–575, 2019.
  * [11] C. Cancès and B. Gaudeul.  A convergent entropy diminishing finite volume scheme for a cross-diffusion system.  SIAM J. Numer. Anal., 58(5):2684–2710, 2020.
  * [12] J. A. Carrillo, A. Chertock, and Y. Huang.  A finite-volume method for nonlinear nonlocal equations with a gradient flow structure.  Commun. Comput. Phys., 17(1):233–258, 2015.
  * [13] J. A. Carrillo, Y. Huang, and M. Schmidtchen.  Zoology of a nonlocal cross-diffusion model for two species.  SIAM J. Appl. Math., 78(2):1078–1104, 2018.
  * [14] J.A. Carrillo, F. Filbet, and M. Schmidtchen.  Convergence of a finite volume scheme for a system of interacting species with cross-diffusion.  Numer. Math., 145(3):473–511, 2020.
  * [15] C. Chainais-Hillairet and M. Herda.  Large-time behaviour of a family of finite volume schemes for boundary-driven convection-diffusion equations.  IMA J. Numer. Anal., 40(4):2473–2504, 2020.
  * [16] Claire Chainais-Hillairet, Maxime Herda, Simon Lemaire, and Julien Moatti.  Long-time behaviour of hybrid finite volume schemes for advection–diffusion equations: linear and nonlinear approaches.  Numerische Mathematik, pages 1–54, 2022.
  * [17] L. Chen, E.S. Daus, A. Holzinger, and A. Jüngel.  Rigorous derivation of population cross-diffusion systems from moderately interacting particle systems.  J. Nonlinear Sci., 31(6):38, 2021.  Id/No 94.
  * [18] L. Chen and A. Jüngel.  Analysis of a multidimensional parabolic population model with strong cross-diffusion.  SIAM J. Math. Anal., 36(1):301–322, 2004.
  * [19] L. Chen and A. Jüngel.  Analysis of a parabolic cross-diffusion population model without self-diffusion.  J. Differ. Equations, 224(1):39–59, 2006.
  * [20] L. Desvillettes, T. Lepoutre, A. Moussa, and A. Trescases.  On the entropic structure of reaction-cross diffusion systems.  Commun. Partial Differ. Equations, 40(9):1705–1747, 2015.
  * [21] H. Dietert and A. Moussa.  Persisting entropy structure for nonlocal cross-diffusion systems.  arXiv preprint arXiv:2101.02893, 2021.
  * [22] P. Domschke, D. Trucu, A. Gerisch, and M. Chaplain.  Mathematical modelling of cancer invasion: implications of cell adhesion variability for tumour infiltrative growth patterns.  J. Theor. Biol., 361:41–60, 2014.
  * [23] B. Düring, , P. Markowich, J.-F. Pietschmann, and M.-T. Wolfram.  Boltzmann and Fokker-Planck equations modelling opinion formation in the presence of strong leaders.  Proc. R. Soc. Lond., Ser. A, Math. Phys. Eng. Sci., 465(2112):3687–3708, 2009.
  * [24] Lawrence C. Evans.  Partial differential equations, volume 19 of Grad. Stud. Math. Providence, RI: American Mathematical Society (AMS), 2nd ed. edition, 2010\.
  * [25] F. Filbet and M. Herda.  A finite volume scheme for boundary-driven convection-diffusion equations with relative entropy structure.  Numer. Math., 137(3):535–577, 2017.
  * [26] J. Fontbona and S. Méléard.  Non local Lotka-Volterra system with cross-diffusion in an heterogeneous medium.  J. Math. Biol., 70(4):829–854, 2015.
  * [27] G. Galiano, M. Garzón, and A. Jüngel.  Semi-discretization in time and numerical convergence of solutions of a nonlinear cross-diffusion population model.  Numer. Math., 93(4):655–673, 2003.
  * [28] V. Giunta, T. Hillen, M.A. Lewis, and J.R. Potts.  Local and global existence for non-local multi- species advection-diffusion models.  Submitted for publication, arXiv:2106.06383, 2021.
  * [29] A. Jüngel.  The boundedness-by-entropy method for cross-diffusion systems.  Nonlinearity, 28(6):1963–2001, 2015.
  * [30] A. Jüngel.  Entropy methods for diffusive partial differential equations.  Cham: Springer; Bilbao: BCAM – Basque Center for Applied Mathematics, 2016.
  * [31] A. Jüngel, S. Portisch, and A. Zurek.  Nonlocal cross-diffusion systems for multi-species populations and networks.  Nonlin. Anal., 219:112800, 2022.
  * [32] A. Jüngel and A. Zurek.  A convergent structure-preserving finite-volume scheme for the Shigesada-Kawasaki-Teramoto population system.  SIAM J. Numer. Anal., 59(4):2286–2309, 2021.
  * [33] S.N. Kruzhkov.  Results concerning the nature of the continuity of solutions of parabolic equations and some of their applications.  Mathematical Notes of the Academy of Sciences of the USSR, 6(1):517–523, 1969.
  * [34] T. Lepoutre and A. Moussa.  Entropic structure and duality for multiple species cross-diffusion systems.  Nonlinear Anal., Theory Methods Appl., Ser. A, Theory Methods, 159:298–315, 2017.
  * [35] Simon A Levin and Lee A Segel.  Hypothesis for origin of planktonic patchiness.  Nature, 259(5545):659–659, 1976.
  * [36] Masayasu Mimura, Yasumasa Nishiura, and Masaya Yamaguti.  Some diffusive prey and predator systems and their bifurcation problems.  Bifurcation theory and applications in scientific disciplines, Pap. Conf., New York 1977, Ann. New York Acad. Sci., Vol. 316, 490-510 (1979)., 1979.
  * [37] A. Moussa.  From nonlocal to classical Shigesada-Kawasaki-Teramoto systems: triangular case with bounded coefficients.  SIAM J. Math. Anal., 52(1):42–64, 2020.
  * [38] H. Murakawa and H. Togashi.  Continuous models for cell-cell adhesion.  J. Theor. Biol., 374:1–12, 2015.
  * [39] K.J. Painter, J.M. Bloomfield, J.A. Sherratt, and A. Gerisch.  A nonlocal model for contact attraction and repulsion in heterogeneous cell populations.  Bull. Math. Biol., 77(6):1132–1165, 2015.
  * [40] J.R. Potts and M.A. Lewis.  Spatial memory and taxis-driven pattern formation in model ecosystems.  Bull. Math. Biol., 81(7):2725–2747, 2019.
  * [41] N. Shigesada, K. Kawasaki, and E. Teramoto.  Spatial segregation of interacting species.  J. Theor. Biol., 79(1):83–99, 1979.
  * [42] Z. Sun, J.A. Carrillo, and C.-W. Shu.  An entropy stable high-order discontinuous Galerkin method for cross-diffusion gradient flow systems.  Kinet. Relat. Models, 12(4):885–908, 2019.