[2]⟨#1,#2⟩ ()^+ ()^- [2]∥#2∥_L^#1 [3]∥#3∥_L^#1(#2) [2]∥#1∥_L^*#2^#1
[3]∥#2∥_L^*(#1)#3^#1 ref eqref Cref

# Bound-Preserving Finite-Volume Schemes for Systems of Continuity Equations
with Saturation

Rafael Bailo  Mathematical Institute, University of Oxford  José A. Carrillo
Mathematical Institute, University of Oxford  Jingwei Hu  Department of
Applied Mathematics, University of Washington

###### Abstract

We propose finite-volume schemes for general continuity equations which
preserve positivity and global bounds that arise from saturation effects in
the mobility function. In the particular case of gradient flows, the schemes
dissipate the free energy at the fully discrete level. Moreover, these schemes
are generalised to coupled systems of non-linear continuity equations, such as
multispecies models in mathematical physics or biology, preserving the bounds
and the dissipation of the energy whenever applicable. These results are
illustrated through extensive numerical simulations which explore known
behaviours in biology and showcase new phenomena not yet described by the
literature.

AMS Subject Classification —  65M08; 35Q70; 35Q92; 45K05.

Keywords —  Gradient flows; systems of conservation laws; finite-volume
methods; bound preservation; energy dissipation; integro-differential
equations.

##  1 Introduction

Systems of aggregation-diffusion equations are ubiquitous in science and
engineering, particularly in mathematical biology, where population dynamics
models have a distinctive importance. These systems faithfully reproduce
attraction effects, usually modelled by non-local terms; repulsion effects,
typically encoded by non-linear diffusion and cross-diffusion; and global
constraints, usually due to volume exclusion effects within the populations.
The balance and interplay between these terms, which leads interesting
phenomena in single population models [10, 9], can produce even richer
behaviours in systems of multiple species [22, 14]. For instance, cell-sorting
phenomena, triggered by differential adhesion, has been reported in
aggregation-reaction-diffusion systems [6, 28]; this biological behaviour is
sharply captured by these models, matching experimental data well [14].

In their most basic form, these models take the form of a system of
conservation laws:

$\begin{cases}\partial_{t}\boldsymbol{\rho}+\nabla\cdot\left(M(\boldsymbol{\rho})\psi(\sigma)\boldsymbol{u}(\boldsymbol{x})\right)\coloneqq\partial_{t}\boldsymbol{\rho}+\sum_{l=1}^{d}\partial_{x_{l}}\left(M(\boldsymbol{\rho})\psi\boldsymbol{u}_{x_{l}}\right)=0,&t>0,\,\boldsymbol{x}\in\Omega\subseteq{\mathbb{R}^{d}},\\\
\boldsymbol{\rho}(0,\boldsymbol{x})=\boldsymbol{\rho}_{0}(\boldsymbol{x}),\end{cases}$
(1.1)

where $\boldsymbol{\rho}$ is a vector of the densities of the $P$ species,
$\boldsymbol{\rho}=(\rho_{1},\rho_{2},\cdots,\rho_{P})^{\top}$ and
$\rho_{p}(t,\boldsymbol{x}):{\mathbb{R}^{+}}\times\Omega\rightarrow{\mathbb{R}^{+}}$.
Here, $M(\boldsymbol{\rho})$ is a $P\times P$, positive-semi-definite
(possibly non-symmetric) matrix. The velocity $\boldsymbol{u}(\boldsymbol{x})$
comprises spatial components $\boldsymbol{u}_{x_{l}}$ along the
$x_{l}$-direction; each is a vector-valued function
$\boldsymbol{u}_{x_{l}}(\boldsymbol{x}):\Omega\rightarrow{\mathbb{R}}^{P}$.
Additionally, $\psi(\sigma):{\mathbb{R}}^{+}\rightarrow{\mathbb{R}}$ is a
saturation, viz. Definition 2.1: a decreasing function of the total density of
the system,
$\sigma(t,\boldsymbol{x})\coloneqq\sum_{p=1}^{P}\rho_{p}(t,\boldsymbol{x})$,
such that $\psi(\alpha)=0$ at a given saturation level $\alpha$.

More complex models require velocity fields which depend (non-linearly) on the
density vector. The theory of systems of conservation laws is not yet ready
for a general, density-dependent $\boldsymbol{u}$. However, there is a family
of models which are better understood: the gradient flow generalisation of
(1.1), given by

$\displaystyle\partial_{t}\boldsymbol{\rho}=\nabla\cdot\left(M(\boldsymbol{\rho})\psi(\sigma)\nabla\left(\frac{\delta
E}{\delta\boldsymbol{\rho}}\right)\right)=\sum_{l=1}^{d}\partial_{x_{l}}\left(M(\boldsymbol{\rho})\psi(\sigma)\partial_{x_{l}}\left(\frac{\delta
E}{\delta\boldsymbol{\rho}}\right)\right),$ (1.2)

for a given energy functional $E$. This functional is key to the analysis of
this problem, since it provides a formal dissipation estimate:

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}E[\boldsymbol{\rho}]=\int_{\Omega}\frac{\delta
E}{\delta\boldsymbol{\rho}}\cdot\partial_{t}\boldsymbol{\rho}=-\int_{\Omega}\psi(\sigma)\nabla\left(\frac{\delta
E}{\delta\boldsymbol{\rho}}\right)\cdot
M(\boldsymbol{\rho})\nabla\left(\frac{\delta
E}{\delta\boldsymbol{\rho}}\right)\mathop{}\\!\mathrm{d}\boldsymbol{x}\leq 0,$
(1.3)

since $M(\boldsymbol{\rho})$ is positive semi-definite. The form of the energy
functional can lead to a large range of phenomena, including linear/non-linear
diffusion and local/non-local drifts [12].

The analysis of cross-diffusion systems which include potential or non-local
terms is nevertheless quite challenging. We highlight a series of recent works
[18, 19, 30] which deal with a family of cross-diffusion systems without drift
terms; they exploit the gradient-flow structure of the systems to define
global solutions with global pointwise bounds. The entropy dissipation
techniques used there are not only essential to define the solutions, but also
to determine their long-time asymptotics. Similar results have only been
achieved for systems with drift terms in a handful of specific cases [6, 5].

The design of numerical schemes that preserve the structural properties of the
solutions to (1.1) and (1.2), and thus, that are suitable to explore their
rich qualitative behaviours, is desirable. To that end, there are two crucial
aspects that beg consideration. The first one is the saturation effect, due to
$\psi(\sigma)$, which arises in these models due to volume exclusion effects,
see for instance [1]. The total density of the system has a natural global
bound, which must be preserved by the numerical scheme, along with the non-
negativity of the species’ densities. The second one is the gradient-flow
structure, present in many models. In these cases, a scheme which captures the
dissipation at the fully discrete level (a structure-preserving scheme) is
desirable. We emphasise that the saturation terms affect the rate of energy
dissipation of the system, and therefore their proper discretisation is
required for this goal.

While numerical works for continuity equations are abundant ([8, 23, 26, 2]
among others), few deal with the saturation effects as well as systems [27,
21], and none do so in full generality. The work [1] proposes an upwind scheme
for scalar Fokker-Planck-like equations in chemotaxis, where the saturation
effects only act on the drift term. A series of works by Jüngel and
collaborators ([20] and the references therein) propose finite-volume schemes
for cross-diffusion systems without drifts; there, a hard volume-filling
constraint is imposed by expressing one solution variable in terms of the
rest.

In this work we will propose finite-volume schemes for the general
formulations (1.1) and (1.2) which preserve both the positivity and the global
bounds induced by the saturation effects. In the case of gradient flows,
(1.2), they also guarantee the dissipation of the free energy at the fully
discrete level. Based on our previous work [2], we consider explicit and
implicit schemes where the numerical fluxes are carefully chosen to preserve
the lower and upper bounds, a property which later plays a crucial role in
preserving the gradient-flow structure. Our numerical methods combine
techniques from systems of hyperbolic conservation laws [16, 15], such as
upwinding, with a careful discretization of the velocity fields [8, 2] in
order to capture the desired dissipative structure.

We will demonstrate the prowess of our methods through a series of numerical
examples. We remark that no general theoretical results for equations or
systems with saturation are available in the literature, except for a few
cases [11, 16, 15], making our numerical explorations very interesting in
terms of new phenomena.

Our work is organised in a constructive way. We begin by discussing the case
of scalar continuity equations with saturation and given velocity fields in
Section 2, where we design explicit and implicit numerical schemes, and show
they preserve the lower and upper bounds on the solution. This property is
generalised in Section 3 to the case of coupled systems with saturation on the
total population. Section 4 deals with the special case of gradient flows with
saturation, both for the scalar equation and the system, showing that a
careful discretisation of the velocity field leads to the dissipation of the
free energy at the fully discrete level. To conclude, Section 5 is devoted to
numerical experiments which validate the numerical schemes and showcase their
use. These experiments illustrate novel phenomena which arises from the
saturation effects, not yet described by the available literature.

##  2 Non-Linear Conservation Laws with Saturation

In this section we study the scalar non-linear conservation law

$\begin{cases}\partial_{t}\rho+\nabla\cdot(\rho\psi(\rho)\boldsymbol{u}(\boldsymbol{x}))=0,&t>0,\,\boldsymbol{x}\in\Omega\subseteq{\mathbb{R}^{d}},\\\
\rho(0,\boldsymbol{x})=\rho_{0}(\boldsymbol{x}),\end{cases}$ (2.1)

where the density
$\rho(t,\boldsymbol{x}):{\mathbb{R}^{+}}\times\Omega\rightarrow{\mathbb{R}^{+}}$
is an unknown non-negative function, the initial datum
$\rho_{0}(\boldsymbol{x}):\Omega\rightarrow{\mathbb{R}^{+}}$ is prescribed,
and the velocity
$\boldsymbol{u}(\boldsymbol{x}):\Omega\rightarrow{\mathbb{R}^{d}}$ is a given
vector-valued function. Equation 2.1 can also be interpreted as a non-linear
continuity equation with saturation. The term $\psi(\rho)$ is a saturation of
the density $\rho$ in the following sense:

######  Definition 2.1 (Saturation).

A saturation is a continuous function
$\psi(s):{\mathbb{R}^{+}}\rightarrow{\mathbb{R}}$ with the properties: $\psi$
is non-increasing ($s_{1}\leq s_{2}$ implies $\psi(s_{1})\geq\psi(s_{2})$);
and there exists $\alpha>0$, called the saturation level, such that
$\psi(\alpha)=0$ and $(\alpha-s)\psi(s)>0$ for $s\neq\alpha$.

Typical examples are $\psi(s)=\alpha-s$ and
$\psi(s)=(\alpha-s)|\alpha-s|^{m-1}$, $m>1$, commonly used in population
models in mathematical biology [17, 14], as well as the social sciences [7,
4].

Entropy solutions to the non-linear conservation law (2.1) satisfy a
comparison principle whereby the density remains bounded between $0$ and
$\alpha$, see [16, 15] for instance. We recall:

######  Theorem 2.2 (Entropy solutions with saturation).

Let $\rho$ be an entropy solution to Eq. 2.1 with initial datum
$0\leq\rho_{0}(\boldsymbol{x})\leq\alpha$ for all $\boldsymbol{x}\in\Omega$
and no-flux boundary condition $\rho\psi\boldsymbol{u}\cdot{\bf n}=0$ on
$\partial\Omega$, where $\bf n$ is the unit outward normal vector to
$\partial\Omega$. Assume the flux
$\rho\psi(\rho)\boldsymbol{u}(\boldsymbol{x})$ is locally Lipschitz with
respect to $\rho$ and continuous in $\boldsymbol{x}$. Then the unique entropy
solution satisfies

$0\leq\rho(t,\boldsymbol{x})\leq\alpha,\quad\forall
t>0,\,\boldsymbol{x}\in\Omega,$ (2.2)

and

$\displaystyle\|\rho(t,\cdot)\|_{{L^{1}}}=\|\rho_{0}(\cdot)\|_{{L^{1}}},\quad\forall
t>0.$ (2.3)

Furthermore, if the initial datum satisfies either of the bounds strictly, the
solution satisfies the same bound strictly too.

######  Remark 2.3 (Reduction to one dimension).

The bounds of Theorem 2.2 hold for scalar conservation laws in any dimension.
In the following, we will introduce numerical schemes in one dimension which
preserve said bounds. Said methods easily generalise to higher-dimensional
scalar conservation laws in quadrangular domains through dimensional
splitting; therefore, we safely reduce the discussion to the one-dimensional
case henceforth.

The aim of this section is to develop numerical schemes which preserve the
saturation effect in the family of one-dimensional scalar conservation laws
(2.1), which reads

$\begin{cases}\partial_{t}\rho+\partial_{x}(\rho\psi(\rho)u(x))=0,&t\in(0,T),\,x\in(0,L),\\\
\rho(0,x)=\rho_{0}(x).\end{cases}$ (2.4)

We will prescribe a full discretisation in the finite-volume setting. The
spatial domain is divided into $M$ uniform volumes of size $\Delta x=L/M$. We
shall define the points $x_{i}=\Delta x(i-1/2)$, $i=1,\cdots,M$ to construct
the $i$th cell $C_{i}=(x_{i-1/2},x_{i+1/2})$ with centre $x_{i}$. The time
domain $(0,T)$ is discretised through equispaced points with separation
$\Delta t=T/N$, with the $n$th point given by $t^{n}=n\Delta t$,
$n=0,\cdots,N$.

The solution $\rho(t,x)$ is to be approximated by a function defined at each
time step and constant over each cell; we denote by $\rho_{i}^{n}$ its value
at time $t^{n}$ over the cell $C_{i}$.

###  2.1 Explicit Scheme

We propose an explicit scheme for (2.4) as

$\displaystyle\begin{split}&\frac{\rho_{i}^{n+1}-\rho_{i}^{n}}{\Delta
t}+\frac{F_{i+1/2}^{n}-F_{i-1/2}^{n}}{\Delta x}=0,\\\
&F_{i+1/2}^{n}=\rho_{i}^{E}\psi_{i+1}^{W}u_{i+1/2}^{+}+\rho_{i+1}^{W}\psi_{i}^{E}u_{i+1/2}^{-},\\\
&u_{i+1/2}=u(x_{i+1/2}),\quad(\cdot)^{+}=\max\\{\cdot,\,0\\},\
(\cdot)^{-}=\min\\{\cdot,\,0\\},\\\ &\rho_{i}^{E}=\rho_{i}^{n}+\frac{\Delta
x}{2}(\rho_{x})_{i}^{n},\quad\rho_{i}^{W}=\rho_{i}^{n}-\frac{\Delta
x}{2}(\rho_{x})_{i}^{n},\\\
&(\rho_{x})_{i}^{n}=\operatorname*{minmod}\left(\theta\frac{\rho_{i+1}^{n}-\rho_{i}^{n}}{\Delta
x},\frac{\rho_{i+1}^{n}-\rho_{i-1}^{n}}{2\Delta
x},\theta\frac{\rho_{i}^{n}-\rho_{i-1}^{n}}{\Delta x}\right),\\\
&\psi_{i}^{E}=\psi(\rho_{i}^{E}),\quad\psi_{i}^{W}=\psi(\rho_{i}^{W}).\end{split}$
(2.5)

We shall prescribe no-flux boundary conditions: $F_{1/2}^{n}=F_{M+1/2}^{n}=0$.
The minmod limiter is defined as

$\displaystyle\operatorname*{minmod}(a,b,c)\coloneqq\begin{cases}\min(a,b,c)&\text{if
}a,b,c>0,\\\ \max(a,b,c)&\text{if }a,b,c<0,\\\ 0&\text{otherwise},\end{cases}$
(2.6)

and we choose $\theta=2$. This scheme is first-order accurate in time and
second-order accurate in space.

######  Proposition 2.4 (Boundedness and non-negativity).

For a given $n$, suppose $0\leq\rho_{i}^{n}\leq\alpha$ for all $i$. Then
scheme (2.5) satisfies $0\leq\rho_{i}^{n+1}\leq\alpha$ for all $i$ if

$\Delta t\leq\Gamma\frac{\Delta x}{2\max_{i}|u_{i+1/2}|},$ (2.7)

where

$\Gamma=\min\left\\{\frac{1}{\psi(0)},\gamma\right\\},\quad\gamma=\inf_{s\in[0,\alpha]}\frac{\alpha-s}{\alpha\psi(s)}.$
(2.8)

###### Proof.

We first show that $0\leq\rho_{i}^{n}\leq\alpha$ for all $i$ implies
$0\leq\rho_{i}^{E}\leq\alpha$ and $0\leq\rho_{i}^{W}\leq\alpha$ for all $i$.
To prove so, we will study $(\rho_{x})_{i}^{n}$, which may only be non-zero if
$\rho_{i-1}^{n}<\rho_{i}^{n}<\rho_{i+1}^{n}$ or
$\rho_{i-1}^{n}>\rho_{i}^{n}>\rho_{i+1}^{n}$ by construction. Assuming that
$\rho_{i-1}^{n}<\rho_{i}^{n}<\rho_{i+1}^{n}$, we find

$\displaystyle
0\leq\rho_{i-1}^{n}\leq\rho_{i}^{n}-\frac{\theta}{2}(\rho_{i}^{n}-\rho_{i-1}^{n})\leq\rho_{i}^{W}\leq\rho_{i}^{n}\leq\rho_{i}^{E}\leq\rho_{i}^{n}+\frac{\theta}{2}(\rho_{i+1}^{n}-\rho_{i}^{n})\leq\rho_{i+1}^{n}\leq\alpha$
(2.9)

since $\theta\in[0,2]$, establishing
$0\leq\rho_{i}^{E},\rho_{i}^{W}\leq\alpha$. Similarly, if
$\rho_{i-1}^{n}>\rho_{i}^{n}>\rho_{i+1}^{n}$, all the inequalities of (2.9)
are reversed, and again we find $0\leq\rho_{i}^{E},\rho_{i}^{W}\leq\alpha$. In
any other case, $0\leq\rho_{i}^{W}=\rho_{i}^{n}=\rho_{i}^{E}\leq\alpha$, which
finally proves the claim.

To address non-negativity, we note that the numerical fluxes admit the bounds

$\displaystyle\rho_{i+1}^{W}\psi_{i}^{E}u_{i+1/2}^{-}\leq
F_{i+1/2}^{n}\leq\rho_{i}^{E}\psi_{i+1}^{W}u_{i+1/2}^{+},$ (2.10)

since $\rho_{i+1}^{W}\psi_{i}^{E}u_{i+1/2}^{-}$ is non-positive and
$\rho_{i}^{E}\psi_{i+1}^{W}u_{i+1/2}^{+}$ is non-negative. Here, we have used
the properties of the saturation $\psi$ listed in Definition 2.1 to conclude
that $0\leq\rho_{i}^{E},\rho_{i}^{W}\leq\alpha$ implies
$0\leq\psi_{i}^{E},\psi_{i}^{W}\leq\psi(0)$. Thus, we write scheme (2.5) in
full and produce a lower bound:

$\displaystyle\rho_{i}^{n+1}$ $\displaystyle=\rho_{i}^{n}-\frac{\Delta
t}{\Delta x}[F_{i+1/2}^{n}-F_{i-1/2}^{n}],$ (2.11)
$\displaystyle\geq\rho_{i}^{n}-\frac{\Delta t}{\Delta
x}[\rho_{i}^{E}\psi_{i+1}^{W}u_{i+1/2}^{+}-\rho_{i}^{W}\psi_{i-1}^{E}u_{i-1/2}^{-}].$
(2.12)

Using the identity $\rho_{i}^{n}=(\rho_{i}^{E}+\rho_{i}^{W})/2$, the
inequality becomes

$\displaystyle\rho_{i}^{n+1}$
$\displaystyle\geq\rho_{i}^{E}\left(\frac{1}{2}-\frac{\Delta t}{\Delta
x}\psi_{i+1}^{W}u_{i+1/2}^{+}\right)+\rho_{i}^{W}\left(\frac{1}{2}+\frac{\Delta
t}{\Delta x}\psi_{i-1}^{E}u_{i-1/2}^{-}\right),$ (2.13)

and the right hand side must be non-negative whenever

$\displaystyle\Delta t\leq\frac{\Delta
x}{2\max\left\\{\psi_{i+1}^{W}u_{i+1/2}^{+},-\psi_{i-1}^{E}u_{i-1/2}^{-}\right\\}}.$
(2.14)

Taking the minimum over $i$ and using $\psi_{i}^{E},\psi_{i}^{W}\leq\psi(0)$
yields

$\Delta t\leq\frac{\Delta x}{2\psi(0)\max_{i}|u_{i+1/2}|}.$ (2.15)

as a sufficient condition for non-negativity.

Regarding boundedness: from the scheme written in full and the flux bounds,
Eqs. 2.11 and 2.10, we can deduce the upper bound

$\displaystyle\rho_{i}^{n+1}$ $\displaystyle\leq\rho_{i}^{n}-\frac{\Delta
t}{\Delta
x}[\rho_{i+1}^{W}\psi_{i}^{E}u_{i+1/2}^{-}-\rho_{i-1}^{E}\psi_{i}^{W}u_{i-1/2}^{+}].$
(2.16)

Again using the identity $\rho_{i}^{n}=(\rho_{i}^{E}+\rho_{i}^{W})/2$, we
obtain

$\displaystyle\rho_{i}^{n+1}$
$\displaystyle\leq\left(\frac{\rho_{i}^{E}}{2}-\frac{\Delta t}{\Delta
x}\rho_{i+1}^{W}\psi_{i}^{E}u_{i+1/2}^{-}\right)+\left(\frac{\rho_{i}^{W}}{2}+\frac{\Delta
t}{\Delta x}\rho_{i-1}^{E}\psi_{i}^{W}u_{i-1/2}^{+}\right);$ (2.17)

we proceed to show that each parenthesis is bounded by $\alpha/2$. In the
leftmost term, either $\rho_{i}^{E}=\alpha$ or $\rho_{i}^{E}<\alpha$. If
$\rho_{i}^{E}=\alpha$, then $\psi_{i}^{E}=0$, and the parenthesis is exactly
equal to $\alpha/2$. If $\rho_{i}^{E}<\alpha$, we require

$\displaystyle\Delta t\leq(\alpha-\rho_{i}^{E})\frac{\Delta
x}{-2\rho_{i+1}^{W}\psi_{i}^{E}u_{i+1/2}^{-}}.$ (2.18)

Using $\rho_{i+1}^{W}\leq\alpha$, we arrive at the stricter condition

$\displaystyle\Delta t\leq\gamma\frac{\Delta x}{2\max_{i}|u_{i+1/2}|},$ (2.19)

where

$\displaystyle\gamma=\inf_{s\in(0,\alpha)}\frac{\alpha-s}{\alpha\psi(s)}.$
(2.20)

The same analysis can be performed for the rightmost parenthesis, thus
establishing Eq. 2.19 as a sufficient condition to guarantee
$\rho_{i}^{n+1}\leq\alpha$.

The combination of Eqs. 2.15 and 2.19 yields Eq. 2.7 as a sufficient condition
for non-negativity and boundedness of the scheme. ∎

######  Remark 2.5 (The character of the saturation, CFL condition).

It is crucial to verify that the quantity $(\alpha-\rho)/{\alpha\psi(\rho)}$
is bounded away from zero in the $\rho\rightarrow\alpha$ limit in order to
ensure that scheme (2.5) can be used in practice. Assuming $\psi$ has a Taylor
expansion about $\alpha$, we consider $\rho=\alpha-\epsilon$ for $\epsilon>0$;
the quantity becomes

$\displaystyle\frac{\alpha-\rho}{\alpha\psi(\rho)}=\frac{\epsilon}{\alpha\psi(\alpha-\epsilon)}=\frac{1}{\alpha\left(-\psi^{\prime}(\alpha)+\frac{\epsilon}{2}\psi^{\prime\prime}(\alpha)-\frac{\epsilon^{2}}{6}\psi^{\prime\prime\prime}(\alpha)+\cdots\right)},$
(2.21)

where we have used $\psi(\alpha)=0$ (viz. Definition 2.1). Therefore, it is
sufficient to require

$\displaystyle(-1)^{k}\psi^{(k)}(\alpha)>0,$ (2.22)

where $\psi^{(k)}$ is the first of the derivatives of $\psi$ to be non-zero at
$\rho=\alpha$. An example which verifies condition (2.22) is $\psi(s)=\cos(\pi
s/2\alpha)$; indeed, one readily finds $\gamma=2/\pi$. Another important
example is $\psi(s)=(\alpha-s)^{m}/\alpha^{m}$ for $m\geq 1$, where we find
$\gamma=1$; conveniently, the CFL condition of Proposition 2.4 reduces to that
of a typical hyperbolic scheme.

###  2.2 Implicit Scheme

In order to bypass the CFL condition altogether, we propose an implicit scheme
for (2.4) as

$\displaystyle\begin{split}&\frac{\rho_{i}^{n+1}-\rho_{i}^{n}}{\Delta
t}+\frac{F_{i+1/2}^{n+1}-F_{i-1/2}^{n+1}}{\Delta x}=0,\\\
&F_{i+1/2}^{n+1}={\rho_{i}^{n+1}}(\psi_{i+1}^{n+1})^{+}u_{i+1/2}^{+}+{\rho_{i+1}^{n+1}}(\psi_{i}^{n+1})^{+}u_{i+1/2}^{-},\\\
&\psi_{i}^{n}=\psi(\rho_{i}^{n}),\quad u_{i+1/2}=u(x_{i+1/2}),\end{split}$
(2.23)

where $(\cdot)^{+}$, $(\cdot)^{-}$, and the boundary conditions are defined as
in Section 2.1. This scheme is first-order accurate in both time and space.

######  Proposition 2.6 (Boundedness and non-negativity).

For a given $n$, suppose $0\leq\rho_{i}^{n}\leq\alpha$ for all $i$. Then
scheme (2.23) satisfies $0\leq\rho_{i}^{n+1}\leq\alpha$ unconditionally for
all $i$.

###### Proof.

We will prove non-negativity followed by boundedness. Seeking a contradiction,
we will assume $\rho_{i}^{n+1}<0$ or $\rho_{i}^{n+1}>\alpha$ for some values
of $i$. Without loss of generality we may assume that these “pathological”
values lie on contiguous cells; if there are two or more separate pathological
clusters, the proof we present below may be applied several times to deal with
each one individually. Thus we may simply assume the values in question are
$\rho_{j}^{n+1},\rho_{j+1}^{n+1},\cdots,\rho_{k}^{n+1}$. Summing scheme (2.23)
over the corresponding cells yields

$\displaystyle\sum_{i=j}^{k}(\rho_{i}^{n+1}-\rho_{i}^{n})\frac{\Delta
x}{\Delta t}$ $\displaystyle=-F_{k+1/2}^{n+1}+F_{j-1/2}^{n+1}.$ (2.24)

To prove non-negativity, we assume $\rho_{i}^{n+1}$ is strictly negative for
$j\leq i\leq k$. Since $\rho_{i}^{n}\geq 0$, the left hand side of Eq. 2.24 is
strictly negative also. The right hand side is comprised of

$\displaystyle-
F_{k+1/2}^{n+1}=-{\rho_{k}^{n+1}}(\psi_{k+1}^{n+1})^{+}u_{k+1/2}^{+}-{\rho_{k+1}^{n+1}}(\psi_{k}^{n+1})^{+}u_{k+1/2}^{-}$
(2.25)

and

$\displaystyle
F_{j-1/2}^{n+1}={\rho_{j-1}^{n+1}}(\psi_{j}^{n+1})^{+}u_{j-1/2}^{+}+{\rho_{j}^{n+1}}(\psi_{j-1}^{n+1})^{+}u_{j-1/2}^{-}.$
(2.26)

The first and fourth terms,
$-{\rho_{k}^{n+1}}(\psi_{k+1}^{n+1})^{+}u_{k+1/2}^{+}$ and
${\rho_{j}^{n+1}}(\psi_{j-1}^{n+1})^{+}u_{j-1/2}^{-}$, are guaranteed to be
non-negative, since $\rho_{k}^{n+1},\rho_{j}^{n+1}<0$; the second and third
terms, $-{\rho_{k+1}^{n+1}}(\psi_{k}^{n+1})^{+}u_{k+1/2}^{-}$ and
${\rho_{j-1}^{n+1}}(\psi_{j}^{n+1})^{+}u_{j-1/2}^{+}$ are also non-negative.
This follows because, outside the range $j\leq i\leq k$, the densities
$\rho_{i}^{n+1}$ are non-negative (whether the upper bound is respected here
is of no consequence). Therefore, the right hand side of Eq. 2.24 is also
shown to be non-negative, producing a contradiction.

Boundedness is proven in a similar fashion. We assume $\rho_{i}^{n+1}$ is
strictly larger than $\alpha$ for $j\leq i\leq k$, resulting on
$\psi_{i}^{n+1}\leq 0$, and rendering the left hand side of Eq. 2.24 strictly
positive. On the right hand side, the second and third terms,
$-{\rho_{k+1}^{n+1}}(\psi_{k}^{n+1})^{+}u_{k+1/2}^{-}$ and
${\rho_{j-1}^{n+1}}(\psi_{j}^{n+1})^{+}u_{j-1/2}^{+}$, are identically zero
because the saturation is non-positive. Furthermore, the first and fourth
terms, $-{\rho_{k}^{n+1}}(\psi_{k+1}^{n+1})^{+}u_{k+1/2}^{+}$ and
${\rho_{j}^{n+1}}(\psi_{j-1}^{n+1})^{+}u_{j-1/2}^{-}$, are non-positive. Thus,
the right hand side is also non-positive, yielding a contradiction. ∎

######  Remark 2.7 (Strict bounds).

It is possible to show a strict inequality version of Proposition 2.6: for a
given $n$, suppose $0<\rho_{i}^{n}<\alpha$ for all $i$; then scheme (2.23)
satisfies $0<\rho_{i}^{n+1}<\alpha$ unconditionally for all $i$. This can be
accomplished through an M-matrix argument, as was done for the positivity
results in [2]; for positivity, one proves $A\rho^{n+1}=\rho^{n}$ where $A$ is
an inverse-positive matrix, whereas for boundedness one shows
$B[\alpha\boldsymbol{1}-\rho^{n+1}]=\alpha\boldsymbol{1}-\rho^{n}$, where $B$
is also inverse positive and $\boldsymbol{1}=(1,\cdots,1)^{\top}$. For the
latter, one needs the additional assumption discussed in Remark 2.5.

######  Remark 2.8 (Reduction to conservation laws without saturation).

By taking $\psi\equiv 1$, the explicit and implicit schemes introduced above
reduce to schemes for the usual advection equation in one dimension,

$\begin{cases}\partial_{t}\rho+\partial_{x}(\rho
u(x))=0,&t\in(0,T),\,x\in(0,L),\\\ \rho(0,x)=\rho_{0}(x).\end{cases}$ (2.27)

Indeed, they can be viewed as the flux splitting finite-difference scheme (see
[25], for instance).

##  3 Non-Linear Systems with Saturation

In this section we turn to the study of systems of conservation laws with $P$
species,

$\begin{cases}\partial_{t}\boldsymbol{\rho}+\nabla\cdot\left(\operatorname{diag}(\boldsymbol{\rho})\psi(\sigma)\boldsymbol{u}(\boldsymbol{x})\right)\coloneqq\partial_{t}\boldsymbol{\rho}+\sum_{l=1}^{d}\partial_{x_{l}}\left(\operatorname{diag}(\boldsymbol{\rho})\psi\boldsymbol{u}_{x_{l}}\right)=0,&t>0,\,\boldsymbol{x}\in\Omega\subseteq{\mathbb{R}^{d}},\\\
\boldsymbol{\rho}(0,\boldsymbol{x})=\boldsymbol{\rho}_{0}(\boldsymbol{x}),\end{cases}$
(3.1)

where the density
$\boldsymbol{\rho}=(\rho_{1},\rho_{2},\cdots,\rho_{P})^{\top}$ is a vector
with entries
$\rho_{p}(t,\boldsymbol{x}):{\mathbb{R}^{+}}\times\Omega\rightarrow{\mathbb{R}^{+}}$
corresponding to the density of the $p$th species.
$\operatorname{diag}(\boldsymbol{\rho})$ is a $P\times P$ diagonal matrix with
diagonal entries given by $\rho_{1},\rho_{2},\cdots,\rho_{P}$. The velocity
$\boldsymbol{u}(\boldsymbol{x})$ now comprises spatial components
$\boldsymbol{u}_{x_{l}}$ along the $x_{l}$-direction; each is a vector-valued
function
$\boldsymbol{u}_{x_{l}}(\boldsymbol{x}):\Omega\rightarrow{\mathbb{R}}^{P}$.
Additionally, $\psi(\sigma):{\mathbb{R}}^{+}\rightarrow{\mathbb{R}}$ is once
again a saturation, viz. Definition 2.1; in this case, the saturation depends
on the total density of the system,
$\sigma(t,\boldsymbol{x})\coloneqq\sum_{p=1}^{P}\rho_{p}(t,\boldsymbol{x})$.

Though the theory of Eq. 3.1 is less developed than that of the scalar case
(2.1), analogous bounds on the solution can be established. The formulation of
the system would suggest that if the initial datum consists of positive
densities whose sum is bounded above by the saturation level $\alpha$, these
bounds will persist for all time. This can be proven, for instance, through
the comparison principle of [16], which we generalise to include the direct
dependency of the flux on the spatial coordinate. Given a system of
conservation laws,

$\displaystyle\partial_{t}\boldsymbol{\rho}+\nabla\cdot\boldsymbol{F}(\boldsymbol{\rho},\boldsymbol{x})=0,$
(3.2)

suppose there exists a value $\rho_{p}^{*}$ such that $\boldsymbol{F}$
satisfies

$\displaystyle\begin{split}\boldsymbol{F}(\rho_{1},\cdots,\rho_{p-1},\rho_{p}^{*},\rho_{p+1},\cdots,\rho_{P},\boldsymbol{x})\equiv
C\end{split}$ (3.3)

for some constant $C$, for any values of
$\rho_{1},\cdots,\rho_{p-1},\rho_{p+1},\cdots,\rho_{P},\boldsymbol{x}$.
Furthermore, assume the necessary smoothness so that the first partial
derivatives of $\boldsymbol{F}$ with respect to $\rho_{p}$ and
$\boldsymbol{x}$ are defined. If $\boldsymbol{\rho}$ is a solution to the
system, and if $\rho_{p}>\rho_{p}^{*}$ holds initially, then it holds for all
times (similarly for $\rho_{p}<\rho_{p}^{*}$).

The result of [16], which immediately yields the positivity of solutions to
Eq. 3.1, can also be used to recover the upper bound on the sum; we illustrate
here the argument in a system of two species in one spatial dimension for the
sake of readability, though the proof can be similarly performed in full
generality. The simplified system reads

$\begin{cases}\partial_{t}\rho+\partial_{x}(\rho\psi(\sigma)u)=0,\\\
\partial_{t}\eta+\partial_{x}(\eta\psi(\sigma)v)=0,\end{cases}$ (3.4)

where $\rho$ and $\eta$ are respectively the first and second species’
densities, $u(x)$ and $v(x)$ are their corresponding velocities, and $\psi$
depends on $\sigma=\rho+\eta$. The key to the upper bound is to replace
$\sigma$ by $\alpha-\theta$, where $\alpha$ is the saturation level and
$\theta$ is a new, independent variable. Letting
$\tilde{\psi}(\theta)=\psi(\alpha-\theta)$, we define a new system

$\begin{cases}\partial_{t}\rho+\partial_{x}[\rho\tilde{\psi}(\theta)u]=0,\\\
\partial_{t}\eta+\partial_{x}[\eta\tilde{\psi}(\theta)v]=0,\\\
\partial_{t}\theta-\partial_{x}[\tilde{\psi}(\theta)\left(\rho u+\eta
v\right)]=0.\end{cases}$ (3.5)

Crucially, any weak solution $(\rho,\eta,\theta)$ with datum
$\rho_{0}+\eta_{0}+\theta_{0}\equiv\alpha$ satisfies
$\rho+\eta+\theta\equiv\alpha$, establishing a one-to-one correspondence
between solutions of Eq. 3.4 and solutions of Eq. 3.5 via
$\sigma\equiv\alpha-\theta$.

The new system is in the correct form to apply the comparison principle of
[16]: the flux,

$\displaystyle\boldsymbol{F}(\rho,\eta,\theta)=\begin{pmatrix}F_{\rho}(\rho,\eta,\theta)\\\
F_{\eta}(\rho,\eta,\theta)\\\
F_{\theta}(\rho,\eta,\theta)\end{pmatrix}=\begin{pmatrix}\rho\tilde{\psi}(\theta)u\\\
\eta\tilde{\psi}(\theta)v\\\ -\tilde{\psi}(\theta)\left(\rho u+\eta
v\right)\end{pmatrix},$ (3.6)

satisfies the condition (3.3). We therefore recover, for any datum
$\rho_{0},\eta_{0},\theta_{0}>0$, a viscosity solution $\rho,\eta,\theta>0$,
understood here as the $\varepsilon\rightarrow 0$ limit of the solutions to
the family of parabolic systems

$\begin{cases}\partial_{t}\rho+\partial_{x}[\rho\tilde{\psi}(\theta)u]=\varepsilon\partial_{xx}\rho,\\\
\partial_{t}\eta+\partial_{x}[\eta\tilde{\psi}(\theta)v]=\varepsilon\partial_{xx}\eta,\\\
\partial_{t}\theta-\partial_{x}[\tilde{\psi}(\theta)\left(\rho u+\eta
v\right)]=\varepsilon\partial_{xx}\theta.\end{cases}$ (3.7)

In particular, solutions to Eq. 3.7 with datum
$\rho_{0}+\eta_{0}+\theta_{0}\equiv\alpha$ correspond to solutions
$(\rho,\eta)$ to Eq. 3.4 which satisfy $0<\rho,\eta$ and $\rho+\eta<\alpha$.
Moreover, these are also viscosity solutions of (3.4), which can be verified
by establishing a similar correspondence between Eq. 3.7 and the parabolic
version of Eq. 3.4,

$\begin{cases}\partial_{t}\rho+\partial_{x}(\rho\psi(\sigma)u)=\varepsilon\partial_{xx}\rho,\\\
\partial_{t}\eta+\partial_{x}(\eta\psi(\sigma)v)=\varepsilon\partial_{xx}\eta.\end{cases}$
(3.8)

In summary:

######  Theorem 3.1 (Viscosity solutions with saturation).

Let $\boldsymbol{\rho}$ be a viscosity solution to Eq. 3.1 with initial datum
$0\leq\rho_{p}(0,\boldsymbol{x})$ for all $p$ and such that
$\sigma(0,\boldsymbol{x})\leq\alpha$ for all $\boldsymbol{x}\in\Omega$, and
with no-flux boundary conditions
$\rho_{p}\psi\boldsymbol{u}_{p}\cdot\boldsymbol{n}=0$ on $\partial\Omega$,
where $\boldsymbol{n}$ is the unit outward normal vector to $\partial\Omega$.
Assume the flux $\operatorname{diag}(\boldsymbol{\rho})\psi\boldsymbol{u}$ is
locally Lipschitz with respect to $\boldsymbol{\rho}$ and continuous in
$\boldsymbol{x}$. Then the viscosity solution satisfies

$0\leq\rho_{p}(t,\boldsymbol{x}),\quad\sigma(t,\boldsymbol{x})\leq\alpha,\quad\forall
p,\,\forall t>0,\,\boldsymbol{x}\in\Omega,$ (3.9)

and

$\displaystyle\|\rho_{p}(t,\cdot)\|_{{L^{1}}}=\|\rho_{p}(0,\cdot)\|_{{L^{1}}}\quad\forall
p,\,\forall t>0.$ (3.10)

Furthermore, if the datum satisfies any of the bounds strictly, then the
solution satisfies the same bound strictly too.

######  Remark 3.2 (General motility matrices).

System (3.1) is the restriction of the system of conservation laws discussed
in the introduction,

$\partial_{t}\boldsymbol{\rho}+\nabla\cdot\left(M(\boldsymbol{\rho})\psi(\sigma)\boldsymbol{u}(\boldsymbol{x})\right)=0,$
(1.1 revisited)

by choosing $M(\boldsymbol{\rho})=\operatorname{diag}(\boldsymbol{\rho})$.
Nevertheless, the schemes that we develop in this section may be used for the
general case by performing a change of variables:

$\displaystyle
M(\boldsymbol{\rho})\psi(\sigma)\boldsymbol{u}(\boldsymbol{x})=\operatorname{diag}(\boldsymbol{\rho})\psi(\sigma)\boldsymbol{v}(\boldsymbol{\rho},\boldsymbol{x}),\quad\text{where
}\boldsymbol{v}(\boldsymbol{\rho},\boldsymbol{x})=-\operatorname{diag}(\boldsymbol{\rho})^{-1}M(\boldsymbol{\rho})\boldsymbol{u}(\boldsymbol{x}).$
(3.11)

The dependency of $\boldsymbol{v}$ on $\boldsymbol{\rho}$ does not present a
problem at the level of the numerical scheme. This technique is detailed at
length in Eq. 4.16, where we do consider a system with a general motility
matrix.

Once again we are concerned with the development of numerical schemes for
(3.1) which preserve the saturation effect. As in the scalar case, we safely
restrict the discussion to one-dimensional problems, and generalise the
schemes through dimensional splitting, viz. Remark 2.3.

In one dimension, system (3.1) considered over a finite domain reads

$\displaystyle\begin{cases}\partial_{t}\boldsymbol{\rho}+\partial_{x}\left(\operatorname{diag}(\boldsymbol{\rho})\psi\boldsymbol{u}\right)=0,&t\in(0,T),\,x\in(0,L),\\\
\boldsymbol{\rho}(0,x)=\boldsymbol{\rho}_{0}(x).\end{cases}$ (3.12)

The finite-volume discretisation will be done similarly as in Section 2.

###  3.1 Explicit Scheme

We propose an explicit scheme for (3.12) as

$\displaystyle\begin{split}&\frac{\boldsymbol{\rho}_{i}^{n+1}-\boldsymbol{\rho}_{i}^{n}}{\Delta
t}+\frac{\boldsymbol{F}_{i+1/2}^{n}-\boldsymbol{F}_{i-1/2}^{n}}{\Delta
x}=0,\\\
&\boldsymbol{F}_{i+1/2}^{n}=\operatorname{diag}(\boldsymbol{\rho}_{i}^{E})\psi_{i+1}^{W}\boldsymbol{u}_{i+1/2}^{+}+\operatorname{diag}(\boldsymbol{\rho}_{i+1}^{W})\psi_{i}^{E}\boldsymbol{u}_{i+1/2}^{-},\\\
&\boldsymbol{u}_{i+1/2}=\boldsymbol{u}(x_{i+1/2}),\quad(\cdot)^{+}=\max\\{\cdot,\,0\\},\
(\cdot)^{-}=\min\\{\cdot,\,0\\},\\\
&\boldsymbol{\rho}_{i}^{E}=\boldsymbol{\rho}_{i}^{n}+\frac{\Delta
x}{2}(\boldsymbol{\rho}_{x})_{i}^{n},\quad\boldsymbol{\rho}_{i}^{W}=\boldsymbol{\rho}_{i}^{n}-\frac{\Delta
x}{2}(\boldsymbol{\rho}_{x})_{i}^{n},\\\
&(\boldsymbol{\rho}_{x})_{i}^{n}=\operatorname*{minmod}\left(\theta\frac{\boldsymbol{\rho}_{i+1}^{n}-\boldsymbol{\rho}_{i}^{n}}{\Delta
x},\frac{\boldsymbol{\rho}_{i+1}^{n}-\boldsymbol{\rho}_{i-1}^{n}}{2\Delta
x},\theta\frac{\boldsymbol{\rho}_{i}^{n}-\boldsymbol{\rho}_{i-1}^{n}}{\Delta
x}\right),\\\
&\psi_{i}^{E}=\psi(\sigma_{i}^{E}),\quad\psi_{i}^{W}=\psi(\sigma_{i}^{W}),\\\
&\sigma_{i}^{E}=\sigma_{i}^{n}+\frac{\Delta
x}{2}(\sigma_{x})_{i}^{n},\quad\sigma_{i}^{W}=\sigma_{i}^{n}-\frac{\Delta
x}{2}(\sigma_{x})_{i}^{n},\quad\sigma_{i}^{n}=\sum_{p=1}^{P}(\boldsymbol{\rho}_{i}^{n})_{p},\\\
&(\sigma_{x})_{i}^{n}=\operatorname*{minmod}\left(\theta\frac{\sigma_{i+1}^{n}-\sigma_{i}^{n}}{\Delta
x},\frac{\sigma_{i+1}^{n}-\sigma_{i-1}^{n}}{2\Delta
x},\theta\frac{\sigma_{i}^{n}-\sigma_{i-1}^{n}}{\Delta x}\right),\\\
\end{split}$ (3.13)

where $\theta=2$ as before. In the above, $(\cdot)^{+}$, $(\cdot)^{-}$, and
the minmod limiter are applied entry-wise to vectors. The no-flux boundary
conditions are prescribed by
$\boldsymbol{F}_{1/2}^{n}=\boldsymbol{F}_{M+1/2}^{n}=0$. This scheme is first-
order accurate in time and second-order accurate in space.

######  Proposition 3.3 (Boundedness and non-negativity).

For a given $n$, suppose $0\leq\boldsymbol{\rho}_{i}^{n}$ (entry-wise) as well
as $\sigma_{i}^{n}\leq\alpha$, for all $i$. Then scheme (3.13) satisfies
$0\leq\boldsymbol{\rho}_{i}^{n+1}$ as well as $\sigma_{i}^{n+1}\leq\alpha$,
for all $i$, if

$\Delta t\leq\frac{\Gamma}{4\|\boldsymbol{u}\|}\Delta x,$ (3.14)

where

$\Gamma=\min\left\\{\frac{2}{\psi(0)},\gamma\right\\},\quad\gamma=\inf_{s\in[0,\alpha]}\frac{\alpha-s}{\alpha\psi(s)},\quad\text{and
}\|\boldsymbol{u}\|=\max_{i,p}|(\boldsymbol{u}_{i+1/2})_{p}|.$ (3.15)

###### Proof.

We first show that $0\leq\boldsymbol{\rho}_{i}^{n}$ and
$\sigma_{i}^{n}\leq\alpha$ imply
$0\leq\sigma_{i}^{E},\sigma_{i}^{W}\leq\alpha$ as well as

$\displaystyle 0\leq\sum_{p=1}^{P}(\boldsymbol{\rho}_{i}^{E})_{p}\leq
2\alpha,\quad 0\leq\sum_{p=1}^{P}(\boldsymbol{\rho}_{i}^{W})_{p}\leq 2\alpha,$
(3.16)

for all $i$. The statement about the total density $\sigma$ is proven using
the technique already used in the proof of Proposition 2.4; namely, using the
fact that the minmod limiter guarantees the monotonicity of the
reconstruction. Equation 3.16 is also proven by monotonicity: we note
$\min\left\\{(\boldsymbol{\rho}_{i}^{n})_{p},(\boldsymbol{\rho}_{i+1}^{n})_{p}\right\\}\leq(\boldsymbol{\rho}_{i}^{E})_{p}\leq\max\left\\{(\boldsymbol{\rho}_{i}^{n})_{p},(\boldsymbol{\rho}_{i+1}^{n})_{p}\right\\}$,
which leads to

$\displaystyle
0\leq\sum_{p=1}^{P}(\boldsymbol{\rho}_{i}^{E})_{p}\leq\sum_{p=1}^{P}\max\left\\{(\boldsymbol{\rho}_{i}^{n})_{p},(\boldsymbol{\rho}_{i+1}^{n})_{p}\right\\}\leq\sum_{p=1}^{P}\left[(\boldsymbol{\rho}_{i}^{n})_{p}+(\boldsymbol{\rho}_{i+1}^{n})_{p}\right]=\sigma_{i}^{n}+\sigma_{i+1}^{n}\leq
2\alpha;$ (3.17)

the estimate for $\boldsymbol{\rho}_{i}^{W}$ is identical. As a remark, note
that $\sum_{p}(\boldsymbol{\rho}_{i}^{E})_{p}\leq\alpha$ is not true in
general.

To address non-negativity, just as in the one-species case, we note that the
numerical fluxes admit the entry-wise bounds

$\displaystyle\operatorname{diag}(\boldsymbol{\rho}_{i+1}^{W})\psi_{i}^{E}\boldsymbol{u}_{i+1/2}^{-}\leq\boldsymbol{F}_{i+1/2}^{n}\leq\operatorname{diag}(\boldsymbol{\rho}_{i}^{E})\psi_{i+1}^{W}\boldsymbol{u}_{i+1/2}^{+},$
(3.18)

since
$\operatorname{diag}(\boldsymbol{\rho}_{i+1}^{W})\psi_{i}^{E}\boldsymbol{u}_{i+1/2}^{-}$
and
$\operatorname{diag}(\boldsymbol{\rho}_{i}^{E})\psi_{i+1}^{W}\boldsymbol{u}_{i+1/2}^{+}$
contain respectively only non-positive and non-negative entries, having
deduced $0\leq\psi_{i}^{E},\psi_{i}^{W}\leq\psi(0)$ from
$0\leq\sigma_{i}^{E},\sigma_{i}^{W}\leq\alpha$. Thus, we write scheme (3.13)
in full and produce a lower bound:

$\displaystyle\boldsymbol{\rho}_{i}^{n+1}$
$\displaystyle=\boldsymbol{\rho}_{i}^{n}-\frac{\Delta t}{\Delta
x}[\boldsymbol{F}_{i+1/2}^{n}-\boldsymbol{F}_{i-1/2}^{n}],$ (3.19)
$\displaystyle\geq\boldsymbol{\rho}_{i}^{n}-\frac{\Delta t}{\Delta
x}[\operatorname{diag}(\boldsymbol{\rho}_{i}^{E})\psi_{i+1}^{W}\boldsymbol{u}_{i+1/2}^{+}-\operatorname{diag}(\boldsymbol{\rho}_{i}^{W})\psi_{i-1}^{E}\boldsymbol{u}_{i-1/2}^{-}].$
(3.20)

Using the identity
$\boldsymbol{\rho}_{i}^{n}=(\boldsymbol{\rho}_{i}^{E}+\boldsymbol{\rho}_{i}^{W})/2$,
the entry-wise inequality becomes

$\displaystyle\boldsymbol{\rho}_{i}^{n+1}$
$\displaystyle\geq\operatorname{diag}(\boldsymbol{\rho}_{i}^{E})\left(\frac{\boldsymbol{1}}{2}-\frac{\Delta
t}{\Delta
x}\psi_{i+1}^{W}\boldsymbol{u}_{i+1/2}^{+}\right)+\operatorname{diag}(\boldsymbol{\rho}_{i}^{W})\left(\frac{\boldsymbol{1}}{2}+\frac{\Delta
t}{\Delta x}\psi_{i-1}^{E}\boldsymbol{u}_{i-1/2}^{-}\right),$ (3.21)

where $\boldsymbol{1}=(1,\cdots,1)^{\top}$ is a $P$-vector of ones. The right
hand side must be non-negative whenever

$\displaystyle\Delta t\leq\frac{\Delta
x}{2\max\left\\{\psi_{i+1}^{W}{(\boldsymbol{u}_{i+1/2})^{+}_{p}},-\psi_{i-1}^{E}{(\boldsymbol{u}_{i+1/2})^{-}_{p}}\right\\}},$
(3.22)

for every $p$. Taking the minimum over $i$ and $p$, and using
$\psi_{i}^{E},\psi_{i}^{W}\leq\psi(0)$ yields

$\Delta t\leq\frac{\Delta
x}{2\psi(0)\max_{i,p}|(\boldsymbol{u}_{i+1/2})_{p}|},$ (3.23)

which is a sufficient condition for non-negativity.

Regarding boundedness: from the scheme written in full and the flux bounds,
Eqs. 3.19 and 3.18, we can deduce the upper bound

$\displaystyle\boldsymbol{\rho}_{i}^{n+1}$
$\displaystyle\leq\boldsymbol{\rho}_{i}^{n}-\frac{\Delta t}{\Delta
x}[\operatorname{diag}(\boldsymbol{\rho}_{i+1}^{W})\psi_{i}^{E}\boldsymbol{u}_{i+1/2}^{-}-\operatorname{diag}(\rho_{i-1}^{E})\psi_{i}^{W}(\boldsymbol{u}_{i-1/2})^{+}].$
(3.24)

Summing over the vector entries, we recover

$\displaystyle\sigma_{i}^{n+1}$ $\displaystyle\leq\sigma_{i}^{n}-\frac{\Delta
t}{\Delta
x}\sum_{p=1}^{P}\left[(\boldsymbol{\rho}_{i+1}^{W})_{p}\psi_{i}^{E}(\boldsymbol{u}_{i+1/2})^{-}_{p}-(\boldsymbol{\rho}_{i-1}^{E})_{p}\psi_{i}^{W}(\boldsymbol{u}_{i-1/2})^{+}_{p}\right],$
(3.25)

which, upon rewriting $\sigma_{i}^{n}=(\sigma_{i}^{E}+\sigma_{i}^{W})/2$,
becomes

$\displaystyle\sigma_{i}^{n+1}$
$\displaystyle\leq\left(\frac{\sigma_{i}^{E}}{2}-\frac{\Delta t}{\Delta
x}\sum_{p=1}^{P}(\boldsymbol{\rho}_{i+1}^{W})_{p}\psi_{i}^{E}(\boldsymbol{u}_{i+1/2})^{-}_{p}\right)+\left(\frac{\sigma_{i}^{W}}{2}+\frac{\Delta
t}{\Delta
x}\sum_{p=1}^{P}(\boldsymbol{\rho}_{i-1}^{E})_{p}\psi_{i}^{W}(\boldsymbol{u}_{i-1/2})^{+}_{p}\right);$
(3.26)

we proceed to show that each parenthesis is bounded by $\alpha/2$. In the
leftmost term, either $\sigma_{i}^{E}=\alpha$ or $\sigma_{i}^{E}<\alpha$. If
$\sigma_{i}^{E}=\alpha$, then $\psi_{i}^{E}=0$, and the parenthesis is exactly
equal to $\alpha/2$. If $\sigma_{i}^{E}<\alpha$, we require

$\displaystyle\Delta t\leq(\alpha-\sigma_{i}^{E})\frac{\Delta
x}{-2\sum_{p}(\boldsymbol{\rho}_{i+1}^{W})_{p}\psi_{i}^{E}(\boldsymbol{u}_{i+1/2})^{-}_{p}}.$
(3.27)

Using $\sum_{p}(\boldsymbol{\rho}_{i}^{W})_{p}\leq 2\alpha$, we arrive at the
stricter condition

$\displaystyle\Delta t\leq\gamma\frac{\Delta
x}{4\max_{i,p}|(\boldsymbol{u}_{i+1/2})_{p}|},$ (3.28)

where

$\displaystyle\gamma=\inf_{s\in(0,\alpha)}\frac{\alpha-s}{\alpha\psi(s)},$
(3.29)

just as in the proof of Proposition 2.4. The same analysis can be performed
for the rightmost parenthesis of Eq. 3.26, thus establishing (3.28) as a
sufficient condition to guarantee $\sigma_{i}^{n+1}\leq\alpha$.

Combining Eqs. 3.23 and 3.28 yields Eq. 3.14 as a sufficient condition for
non-negativity and boundedness of the scheme. ∎

###  3.2 Implicit Scheme

As in the scalar setting, we bypass the CFL condition completely by proposing
an implicit scheme for (3.12) as

$\displaystyle\begin{split}&\frac{\boldsymbol{\rho}_{i}^{n+1}-\boldsymbol{\rho}_{i}^{n}}{\Delta
t}+\frac{\boldsymbol{F}_{i+1/2}^{n+1}-\boldsymbol{F}_{i-1/2}^{n+1}}{\Delta
x}=0,\\\
&\boldsymbol{F}_{i+1/2}^{n+1}=\operatorname{diag}({\boldsymbol{\rho}_{i}^{n+1}})(\psi_{i+1}^{n+1})^{+}\boldsymbol{u}_{i+1/2}^{+}+\operatorname{diag}({\boldsymbol{\rho}_{i+1}^{n+1}})(\psi_{i}^{n+1})^{+}\boldsymbol{u}_{i+1/2}^{-},\\\
&\psi_{i}^{n}=\psi(\sigma_{i}^{n}),\quad\sigma_{i}^{n}=\sum_{p=1}^{P}(\boldsymbol{\rho}_{i}^{n})_{p},\quad\boldsymbol{u}_{i+1/2}=\boldsymbol{u}(x_{i+1/2}),\end{split}$
(3.30)

where $(\cdot)^{+}$, $(\cdot)^{-}$, and the boundary conditions are defined as
in Section 3.1. This scheme has first-order accuracy in time and space.

######  Proposition 3.4 (Boundedness and non-negativity).

For a given $n$, suppose $0\leq\boldsymbol{\rho}_{i}^{n}$ (entry-wise) as well
as $\sigma_{i}^{n}\leq\alpha$, for all $i$. Then scheme (3.30) satisfies
$0\leq\boldsymbol{\rho}_{i}^{n+1}$ as well as $\sigma_{i}^{n+1}\leq\alpha$,
unconditionally for all $i$.

###### Proof.

We will argue by contradiction, proving non-negativity followed by
boundedness. As in the proof of Proposition 2.6, we will only deal with one
cluster of contiguous pathological values, $j\leq i\leq k$. Summing scheme
(3.30) over the corresponding cells yields

$\displaystyle\sum_{i=j}^{k}(\boldsymbol{\rho}_{i}^{n+1}-\boldsymbol{\rho}_{i}^{n})\frac{\Delta
x}{\Delta t}$
$\displaystyle=-\boldsymbol{F}_{k+1/2}^{n+1}+\boldsymbol{F}_{j-1/2}^{n+1}.$
(3.31)

We prove non-negativity one species at a time. Given $p$, suppose
$(\boldsymbol{\rho}_{i}^{n+1})_{p}$ is strictly negative for $j\leq i\leq k$,
and non-negative otherwise. Whether the entries outside the pathological range
satisfy the boundedness property is, for this part, irrelevant, as are the
values of the other species. The $p$th entry of Eq. 3.31 is given by

$\displaystyle\sum_{i=j}^{k}\left[(\boldsymbol{\rho}_{i}^{n+1})_{p}-(\boldsymbol{\rho}_{i}^{n})_{p}\right]\frac{\Delta
x}{\Delta t}$
$\displaystyle=-(\boldsymbol{F}_{k+1/2}^{n+1})_{p}+(\boldsymbol{F}_{j-1/2}^{n+1})_{p}.$
(3.32)

Since $(\boldsymbol{\rho}_{i}^{n})_{p}\geq 0$, the left hand side of this
equation is strictly negative. The right hand side is comprised of

$\displaystyle-(\boldsymbol{F}_{k+1/2}^{n+1})_{p}=-(\boldsymbol{\rho}_{k}^{n+1})_{p}(\psi_{k+1}^{n+1})^{+}(\boldsymbol{u}_{k+1/2})^{+}_{p}-(\boldsymbol{\rho}_{k+1}^{n+1})_{p}(\psi_{k}^{n+1})^{+}(\boldsymbol{u}_{k+1/2})^{-}_{p}$
(3.33)

and

$\displaystyle(\boldsymbol{F}_{j-1/2}^{n+1})_{p}=(\boldsymbol{\rho}_{j-1}^{n+1})_{p}(\psi_{j}^{n+1})^{+}(\boldsymbol{u}_{j-1/2})^{+}_{p}+(\boldsymbol{\rho}_{j}^{n+1})_{p}(\psi_{j-1}^{n+1})^{+}(\boldsymbol{u}_{j-1/2})^{-}_{p}.$
(3.34)

Just as in the proof of Proposition 2.6, the first and fourth terms,
$-(\boldsymbol{\rho}_{k}^{n+1})_{p}(\psi_{k+1}^{n+1})^{+}(\boldsymbol{u}_{k+1/2})^{+}_{p}$
and
$(\boldsymbol{\rho}_{j}^{n+1})_{p}(\psi_{j-1}^{n+1})^{+}(\boldsymbol{u}_{j-1/2})^{-}_{p}$,
are non-negative, since
$(\boldsymbol{\rho}_{k}^{n+1})_{p},(\boldsymbol{\rho}_{j}^{n+1})_{p}<0$; the
second and third terms,
$-(\boldsymbol{\rho}_{k+1}^{n+1})_{p}(\psi_{k}^{n+1})^{+}(\boldsymbol{u}_{k+1/2})^{-}_{p}$
and
$(\boldsymbol{\rho}_{j-1}^{n+1})_{p}(\psi_{j}^{n+1})^{+}(\boldsymbol{u}_{j-1/2})^{+}_{p}$
are guaranteed to be non-negative also. Therefore, the right hand side of Eq.
3.32 is shown to be non-negative, producing a contradiction.

Boundedness is now proven in a similar fashion, by considering instead the sum
of Eq. 3.31 over all species:

$\displaystyle\sum_{i=j}^{k}(\sigma_{i}^{n+1}-\sigma_{i}^{n})\frac{\Delta
x}{\Delta t}$
$\displaystyle=-\sum_{p=1}^{P}(\boldsymbol{F}_{k+1/2}^{n+1})_{p}+\sum_{p=1}^{P}(\boldsymbol{F}_{j-1/2}^{n+1})_{p}.$
(3.35)

Here, we assume that $\sigma_{i}^{n+1}$ is strictly larger than $\alpha$ for
$j\leq i\leq k$, resulting on $\psi_{i}^{n+1}\leq 0$, and rendering the left
hand side of (3.35) strictly positive. The right hand side now comprises four
terms, corresponding to the sums over $p$ of Eqs. 3.33 and 3.34. Two of the
terms,
$-\sum_{p}(\boldsymbol{\rho}_{k+1}^{n+1})_{p}(\psi_{k}^{n+1})^{+}(\boldsymbol{u}_{k+1/2})^{-}_{p}$
and
$\sum_{p}(\boldsymbol{\rho}_{j-1}^{n+1})_{p}(\psi_{j}^{n+1})^{+}(\boldsymbol{u}_{j-1/2})^{+}_{p}$,
are identically zero; the other two,
$-\sum_{p}(\boldsymbol{\rho}_{k}^{n+1})_{p}(\psi_{k+1}^{n+1})^{+}(\boldsymbol{u}_{k+1/2})^{+}_{p}$
and
$\sum_{p}(\boldsymbol{\rho}_{j}^{n+1})_{p}(\psi_{j-1}^{n+1})^{+}(\boldsymbol{u}_{j-1/2})^{-}_{p}$,
are guaranteed to be non-positive. Thus, the right hand side is also non-
positive, yielding a contradiction. ∎

######  Remark 3.5 (Strict bounds).

It is possible to show a strict inequality version of Proposition 3.4, as it
was for Proposition 2.6: for a given $n$, suppose
$0<\boldsymbol{\rho}_{i}^{n}$ (entry-wise) as well as $\sigma_{i}^{n}<\alpha$,
for all $i$; then scheme (3.30) satisfies $0<\boldsymbol{\rho}_{i}^{n+1}$ as
well as $\sigma_{i}^{n+1}<\alpha$, unconditionally for all $i$. Again this is
accomplished through an M-matrix argument, in the style of [2]; see Remark
2.7.

##  4 Saturation in Gradient Flows

The proofs of positivity and boundedness for the schemes introduced in
Sections 2 and 3 rely only on the discretisation of the advecting densities
$\boldsymbol{\rho}$ and the saturation $\psi$; no properties of the velocities
$\boldsymbol{u}$ are used. This suggests that the schemes might admit a
generalisation to systems of continuity equations where the velocity depends
on the solution.

The question of whether solutions to Eqs. 2.1 and 3.1 retain the positivity
and boundedness properties under more general velocities is a subtle one,
particularly for the case of systems, and a general characterisation of
admissible velocity functions is not within the purview of this work. Instead,
we will focus on a class of systems of continuity equations with saturation
and a gradient-flow structure.

In this section, we combine the virtues of the scalar implicit scheme, (2.23),
with the schemes developed in [2], to arrive at a scheme that captures the
positivity and boundedness of solutions while preserving the dissipation of
the free energy at the fully discrete level. Furthermore, we generalise [2] to
consider also systems of equations with a gradient-flow structure, and we
adapt (3.30) to similarly preserve the positivity, boundedness, and
dissipation properties.

###  4.1 Scalar Gradient Flows

We begin the discussion with the scalar case. The recent work [2] studied the
family of equations

$\begin{cases}\partial_{t}\rho=\nabla\cdot[\rho\nabla(H^{\prime}(\rho)+V+W\ast\rho)],&t>0,\,\boldsymbol{x}\in\Omega\subseteq\mathbb{R}^{d},\\\
\rho(0,\boldsymbol{x})=\rho_{0}(\boldsymbol{x}),\end{cases}$ (4.1)

which describe the evolution in time of a density of particles
$\rho(t,\boldsymbol{x})$. The internal energy density $H(\rho)$, a convex
function, models linear or non-linear diffusion on $\rho$; the confining
potential $V(\boldsymbol{x})$ models external forces acting on the density;
the symmetric interaction potential $W(\boldsymbol{x})$, typically a radial
function, models the attraction or repulsion between particles.

This problem, which can be viewed as the continuity equation with velocity
field $\boldsymbol{u}=-\nabla\xi$ for $\xi=H^{\prime}(\rho)+V+W\ast\rho$,
possesses an interesting structural property: it is the gradient flow of the
energy functional

$E[\rho]=\int_{\Omega}[H(\rho)+V\rho+\frac{1}{2}(W\ast\rho)\rho]\mathop{}\\!\mathrm{d}\boldsymbol{x},$
(4.2)

see [12, 13]. In the Wasserstein sense, this means that Eq. 4.1 can be written
as

$\displaystyle\partial_{t}\rho=\nabla\cdot\left(\rho\nabla\left(\frac{\delta
E}{\delta\rho}\right)\right),$ (4.3)

where $\frac{\delta E}{\delta\rho}$ is the first variation of the energy
functional (under zero-mass perturbations). This structure leads to the
dissipation of the energy along solution of the equation,

$\frac{\mathrm{d}E}{\mathrm{d}t}=-\int_{\Omega}\rho|\nabla\xi|^{2}\,\mathop{}\\!\mathrm{d}{\boldsymbol{x}}\leq
0,$ (4.4)

and thus $E$ is a Lyapunov functional for the problem. Due to the relevance of
Eq. 4.1 in applications, numerical solutions which preserve this dissipation
property are desirable; indeed [2] introduces schemes to that effect.

The formulation of Eq. 4.1 can be extended to include the applications
discussed in the introduction which exhibit saturation effects. The new family
of models reads

$\begin{cases}\partial_{t}\rho=\nabla\cdot[\rho\psi(\rho)\nabla(H^{\prime}(\rho)+V+W\ast\rho)],&t>0,\,\boldsymbol{x}\in\Omega\subseteq\mathbb{R}^{d},\\\
\rho(\boldsymbol{x},0)=\rho_{0}(\boldsymbol{x}).\end{cases}$ (4.5)

This equation can be written as

$\displaystyle\partial_{t}\rho=\nabla\cdot\left(\rho\psi(\rho)\nabla\left(\frac{\delta
E}{\delta\rho}\right)\right),$ (4.6)

and thus the energy (4.2) is still formally dissipated along solutions at the
rate

$\frac{\mathrm{d}E}{\mathrm{d}t}=-\int_{\Omega}\rho\psi(\rho)|\nabla\xi|^{2}\,\mathop{}\\!\mathrm{d}{\boldsymbol{x}}\leq
0.$ (4.7)

Specific cases of gradient flows have been dealt with in [11] due to their
convexity properties. We can address the question of positivity and
boundedness of solutions to this equation by verifying the properties
discussed at the beginning of Section 3, in order to apply the comparison
principle of [16]. Once again we limit the discussion to one spatial dimension
for readability. Rewriting the equation as

$\displaystyle\left\\{\begin{array}[]{{>{\displaystyle\vphantom{\frac{1}{1}}}l}}\partial_{t}\rho+\partial_{x}F(\rho,x)=0,\\\
F(\rho,x)=\rho\psi(\rho)u(\rho,x),\\\
u(\rho,x)=-\partial_{x}(H^{\prime}(\rho)+V(x)+W\ast\rho),\end{array}\right.$
(4.11)

we must assume enough regularity on $H$, $V$, and $W$ to ensure continuity of
$F(\rho,x)$ at $\rho=0$ and $\rho=\alpha$, and that $\partial_{x}u$ is
defined.

In applications, the diffusion term is typically $H(\rho)=\rho\log\rho-\rho$
(the Boltzmann entropy, yielding linear diffusion) or $H(\rho)=\rho^{m}/(m-1)$
for some $m>1$ (non-linear diffusion). These forms present no difficulty at
$\rho=\alpha$, but require attention at $\rho=0$. In the linear case, it is
preferable to rewrite the equation as
$\partial_{t}\rho=\partial_{x}\left(\psi(\rho)\partial_{x}\rho\right)$ and to
apply parabolic theory in order to obtain lower bounds. In the non-linear
case, the question of regularity at the origin is quite involved, see [29] for
an exhaustive survey. This complexity permeates also the analysis of numerical
schemes, as explored in [3], which establishes the convergence of the scheme
of [2].

Concerning the potentials, it is sufficient to assume $V,W\in C^{2}(\Omega)$.
This is stringent and could be relaxed in some cases; for instance, if the
form of $H$ leads to establishing some regularity on the solution, it might be
possible to interpret $\partial_{x}(-\partial_{x}(W\ast\rho))$ as
$-\partial_{x}W\ast\partial_{x}\rho$, requiring milder assumptions on $W$.

Once the saturation bound has been established on the solutions of Eq. 4.5, we
can adapt the numerical schemes from [2] to the problem with saturation by
including the discretisation of $\psi$ as above. In the interest of brevity,
we describe the adaptation of one scheme only, the so-called scheme S2; the
remaining schemes can be found analogously. The scheme reads:

$\displaystyle\begin{split}&\frac{\rho_{i}^{n+1}-\rho_{i}^{n}}{\Delta
t}+\frac{F_{i+1/2}^{n+1}-F_{i-1/2}^{n+1}}{\Delta x}=0,\\\
&F_{i+1/2}^{n+1}=\rho_{i}^{n+1}(\psi_{i+1}^{n+1})^{+}(u_{i+1/2}^{n+1})^{+}+\rho_{i+1}^{n+1}(\psi_{i}^{n+1})^{+}(u_{i+1/2}^{n+1})^{-},\\\
&u_{i+1/2}^{n+1}=-\frac{\xi_{i+1}^{n+1}-\xi_{i}^{n+1}}{\Delta x},\quad\\\
&\xi_{i}^{n+1}=H^{\prime}(\rho_{i}^{n+1})+V_{i}+(W\ast\rho^{**})_{i},\end{split}$
(4.12)

for a symmetric potential $W$. The confining potential terms are defined as
$V_{i}=V(x_{i})$ and the convolution is given by
$(W\ast\rho^{**})_{i}=\sum_{k}W_{i-k}\rho_{k}^{**}\Delta x$, where
$W_{i-k}=W(x_{i}-x_{k})$ and $\rho^{**}=(\rho^{n+1}+\rho^{n})/{2}$. We have
discretised space and time as in the previous sections, with one difference:
we assume the domain is of the form $\Omega=(-L,L)$, and discretise it into
$2M$ uniform volumes of size $\Delta x=L/M$.

We present the energy dissipation result as a straightforward corollary to [2,
Theorem 3.9]:

######  Corollary 4.1 (Energy dissipation).

Scheme (4.12) satisfies the energy dissipation property unconditionally, i.e.,
the discrete energy defined by

$E_{\Delta}[\rho^{n}]=\sum_{i=1}^{2M}H(\rho_{i}^{n})\Delta
x+\sum_{i=1}^{2M}V_{i}\rho_{i}^{n}\Delta
x+\frac{1}{2}\sum_{i=1}^{2M}\sum_{k=1}^{2M}W_{i-k}\rho_{i}^{n}\rho_{k}^{n}\Delta
x^{2},$ (4.13)

satisfies $E_{\Delta}[\rho^{n+1}]\leq E_{\Delta}[\rho^{n}]$.

###### Proof.

The proof of the original theorem remains valid until the last two lines,
where the inequality becomes

$\displaystyle\frac{E_{\Delta}[\rho^{n+1}]-E_{\Delta}[\rho^{n}]}{\Delta
t\Delta x}\leq$
$\displaystyle-\sum_{i=1}^{2M-1}\left(\rho_{i}^{n+1}(\psi_{i+1}^{n+1})^{+}(u_{i+1/2}^{n+1})^{+}+\rho_{i+1}^{n+1}(\psi_{i}^{n+1})^{+}(u_{i+1/2}^{n+1})^{-}\right)u_{i+1/2}^{n+1},$
(4.14) $\displaystyle\leq$
$\displaystyle-\sum_{i=1}^{2M-1}\min\left\\{\rho_{i}^{n+1}(\psi_{i+1}^{n+1})^{+},\rho_{i+1}^{n+1}(\psi_{i}^{n+1})^{+}\right\\}|u_{i+1/2}^{n+1}|^{2}\leq
0,$ (4.15)

by virtue of Proposition 2.6. ∎

###  4.2 Systems with a Gradient Flow Structure

The gradient flow (4.5) can be generalised to a system of multiple species:

$\displaystyle\partial_{t}\boldsymbol{\rho}=\nabla\cdot\left(M(\boldsymbol{\rho})\psi(\sigma)\nabla\left(\frac{\delta
E}{\delta\boldsymbol{\rho}}\right)\right)=\sum_{l=1}^{d}\partial_{x_{l}}\left(M(\boldsymbol{\rho})\psi(\sigma)\partial_{x_{l}}\left(\frac{\delta
E}{\delta\boldsymbol{\rho}}\right)\right),$ (4.16)

where, as in Section 3, $\boldsymbol{\rho}=(\rho_{1},\dots,\rho_{P})^{T}$, is
a vector of densities, $\psi(\sigma):{\mathbb{R}}^{+}\rightarrow{\mathbb{R}}$
is a saturation which depends on
$\sigma(t,\boldsymbol{x})\coloneqq\sum_{p=1}^{P}\rho_{p}(t,\boldsymbol{x})$,
and $M(\boldsymbol{\rho})$ is a $P\times P$, positive-semi-definite (possibly
non-symmetric) matrix. This equation will formally verify a dissipation
property:

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}E[\boldsymbol{\rho}]=\int_{\Omega}\frac{\delta
E}{\delta\boldsymbol{\rho}}\cdot\partial_{t}\boldsymbol{\rho}=-\int_{\Omega}\psi(\sigma)\nabla\left(\frac{\delta
E}{\delta\boldsymbol{\rho}}\right)\cdot
M(\boldsymbol{\rho})\nabla\left(\frac{\delta
E}{\delta\boldsymbol{\rho}}\right)\mathop{}\\!\mathrm{d}\boldsymbol{x}\leq 0,$
(4.17)

since $M(\boldsymbol{\rho})$ is positive semi-definite.

For brevity, we consider here the case of two species, letting
$\boldsymbol{\rho}=(\rho,\eta)^{\top}$, and $\sigma=\rho+\eta$. We define the
energy functional as

$E[\boldsymbol{\rho}]=\mathcal{H}[\boldsymbol{\rho}]+\mathcal{V}[\boldsymbol{\rho}]+\mathcal{W}[\boldsymbol{\rho}],$
(4.18)

where

$\displaystyle\mathcal{H}[\boldsymbol{\rho}]$
$\displaystyle=\int_{\Omega}\left[H_{\rho}(\rho)+H_{\eta}(\eta)+H_{\sigma}(\sigma)\right]\mathop{}\\!\mathrm{d}x,$
(4.19) $\displaystyle\mathcal{V}[\boldsymbol{\rho}]$
$\displaystyle=\int_{\Omega}\left[V_{\rho}\rho+V_{\eta}\eta\right]\mathop{}\\!\mathrm{d}x,$
(4.20) $\displaystyle\mathcal{W}[\boldsymbol{\rho}]$
$\displaystyle=\int_{\Omega}\left[\frac{1}{2}\rho(W_{\rho}\ast\rho)+\frac{1}{2}\eta(W_{\eta}\ast\eta)+\rho(W_{\sigma}\ast\eta)\right]\mathop{}\\!\mathrm{d}x.$
(4.21)

The internal energy densities $H_{\rho}(\rho)$, $H_{\eta}(\eta)$, and
$H_{\sigma}(\sigma)$ are convex functions which model diffusion, respectively,
on the first species, on the second, and on the total density of the system.
The confining potentials $V_{\rho}$ and $V_{\eta}$ model external forces
acting on each of the species. The interaction potentials $W_{\rho}$,
$W_{\eta}$, and $W_{\sigma}$, model attraction or repulsion between particles;
$W_{\rho}$ and $W_{\eta}$ represent the forces within each of the species,
whereas $W_{\sigma}$ models the interaction between the species, assumed here
to be symmetric.

The inclusion of the saturation term, $\psi(\sigma)$, ensures that the
gradient flow can be cast in the form of Eq. 3.12; in one dimension,

$\displaystyle\left\\{\begin{array}[]{{>{\displaystyle\vphantom{\frac{1}{1}}}l}}\partial_{t}\boldsymbol{\rho}+\partial_{x}\boldsymbol{F}(\boldsymbol{\rho},x)=0,\\\
\boldsymbol{F}(\boldsymbol{\rho},x)=\operatorname{diag}(\boldsymbol{\rho})\psi\boldsymbol{u}(\boldsymbol{\rho},x),\\\\[5.0pt]
\boldsymbol{u}(\boldsymbol{\rho},x)=-\operatorname{diag}(\boldsymbol{\rho})^{-1}M(\boldsymbol{\rho})\partial_{x}\left(\frac{\delta
E}{\delta\boldsymbol{\rho}}\right).\end{array}\right.$ (4.25)

Once again, in order to apply the comparison principle of [16], we must verify
that $\boldsymbol{F}(\boldsymbol{\rho},x)$ is continuous at $\rho=0$, at
$\eta=0$, and at $\sigma=\alpha$, and that
$\partial_{x}\boldsymbol{u}(\boldsymbol{\rho},x)$ is defined; see the
discussion in the previous section.

If the positivity and boundedness of solutions to Eq. 4.16 are established, we
can extend scheme (3.30) to capture these properties while preserving the
variational structure. The new scheme reads:

$\displaystyle\begin{split}&\frac{\boldsymbol{\rho}_{i}^{n+1}-\boldsymbol{\rho}_{i}^{n}}{\Delta
t}+\frac{\boldsymbol{F}_{i+1/2}^{n+1}-\boldsymbol{F}_{i-1/2}^{n+1}}{\Delta
x}=0,\\\
&\boldsymbol{F}_{i+1/2}^{n+1}=\operatorname{diag}({\boldsymbol{\rho}_{i}^{n+1}})(\psi_{i+1}^{n+1})^{+}(\boldsymbol{u}_{i+1/2}^{n+1})^{+}+\operatorname{diag}({\boldsymbol{\rho}_{i+1}^{n+1}})(\psi_{i}^{n+1})^{+}(\boldsymbol{u}_{i+1/2}^{n+1})^{-},\\\
&\boldsymbol{u}_{i+1/2}^{n+1}=-S_{i+1/2}^{n+1}\frac{\boldsymbol{\xi}_{i+1}^{n+1}-\boldsymbol{\xi}_{i}^{n+1}}{\Delta
x},\\\
&S_{i+1/2}^{n+1}=\frac{1}{2}\left(\operatorname{diag}(\boldsymbol{\rho}_{i}^{n+1})^{-1}M(\boldsymbol{\rho}_{i}^{n+1})+\operatorname{diag}(\boldsymbol{\rho}_{i+1}^{n+1})^{-1}M(\boldsymbol{\rho}_{i+1}^{n+1})\right),\\\
&\boldsymbol{\xi}_{i}^{n+1}=\left(\xi_{\rho,i}^{n+1},\xi_{\eta,i}^{n+1}\right)^{\top},\\\
&\xi_{\rho,i}^{n+1}=H_{\rho}^{\prime}(\rho_{i}^{n+1})+H_{\sigma}^{\prime}(\sigma_{i}^{n+1})+V_{\rho,i}+(W_{\rho}\ast\rho^{**})_{i}+(W_{\sigma}\ast\eta^{**})_{i},\\\
&\xi_{\eta,i}^{n+1}=H_{\eta}^{\prime}(\eta_{i}^{n+1})+H_{\sigma}^{\prime}(\sigma_{i}^{n+1})+V_{\eta,i}+(W_{\eta}\ast\eta^{**})_{i}+(W_{\sigma}\ast\rho^{**})_{i},\\\
\end{split}$ (4.26)

for symmetric potentials $W_{\rho}$, $W_{\eta}$, and $W_{\sigma}$. The
definitions of $V_{\rho,i}$, $V_{\eta,i}$, $(W_{\rho}\ast\rho^{**})_{i}$,
$(W_{\eta}\ast\eta^{**})_{i}$, $(W_{\sigma}\ast\rho^{**})_{i}$,
$(W_{\sigma}\ast\eta^{**})_{i}$, $\rho^{**}$, and $\eta^{**}$ are the
immediate generalisations of the terms of scheme (4.12), and the spatial
discretisation is identical.

This scheme preserves the variational structure, though the argument is more
involved than the scalar case. The discrete analogue of the energy (4.18) is
given by

$E_{\Delta}[\boldsymbol{\rho}^{n}]=\mathcal{H}_{\Delta}[\boldsymbol{\rho}^{n}]+\mathcal{V}_{\Delta}[\boldsymbol{\rho}^{n}]+\mathcal{W}_{\Delta}[\boldsymbol{\rho}^{n}],$
(4.27)

where

$\displaystyle\mathcal{H}_{\Delta}[\boldsymbol{\rho}^{n}]$
$\displaystyle=\sum_{i=1}^{2M}\left[H_{\rho}(\rho_{i}^{n})+H_{\eta}(\eta_{i}^{n})+H_{\sigma}(\sigma_{i}^{n})\right]\Delta
x,$ (4.28) $\displaystyle\mathcal{V}_{\Delta}[\boldsymbol{\rho}^{n}]$
$\displaystyle=\sum_{i=1}^{2M}\left[V_{\rho,i}\rho_{i}^{n}+V_{\eta,i}\eta_{i}^{n}\right]\Delta
x,$ (4.29) $\displaystyle\mathcal{W}_{\Delta}[\boldsymbol{\rho}^{n}]$
$\displaystyle=\frac{1}{2}\sum_{i=1}^{2M}\sum_{k=1}^{2M}\left[\rho_{i}^{n}\rho_{k}^{n}W_{\rho,i-k}+\eta_{i}^{n}\eta_{k}^{n}W_{\eta,i-k}+2\rho_{i}^{n}\eta_{k}^{n}W_{\sigma,i-k}\right]\Delta
x^{2}.$ (4.30)

######  Theorem 4.2 (Energy Dissipation).

Scheme (4.26) satisfies the energy dissipation property unconditionally, i.e.,
the discrete energy (4.27) satisfies $E_{\Delta}[\boldsymbol{\rho}^{n+1}]\leq
E_{\Delta}[\boldsymbol{\rho}^{n}]$.

###### Proof.

Multiply the scheme by the entropy variable and integrate in space to find

$\displaystyle\sum_{i=1}^{2M}\boldsymbol{\xi}_{i}^{n+1}\cdot(\boldsymbol{\rho}_{i}^{n+1}-\boldsymbol{\rho}_{i}^{n})=-\frac{\Delta
t}{\Delta
x}\sum_{i=1}^{2M}\boldsymbol{\xi}_{i}^{n+1}\cdot(\boldsymbol{F}_{i+1/2}^{n+1}-\boldsymbol{F}_{i-1/2}^{n+1}),$
(4.31)

which can be rewritten as

$\displaystyle\sum_{i=1}^{2M}\boldsymbol{V}_{i}\cdot(\boldsymbol{\rho}_{i}^{n+1}-\boldsymbol{\rho}_{i}^{n})$
(4.32) $\displaystyle=-\frac{\Delta t}{\Delta
x}\sum_{i=1}^{2M}\boldsymbol{\xi}_{i}^{n+1}\cdot(\boldsymbol{F}_{i+1/2}^{n+1}-\boldsymbol{F}_{i-1/2}^{n+1})-\sum_{i=1}^{2M}H_{\rho}^{\prime}(\rho_{i}^{n+1})(\rho_{i}^{n+1}-\rho_{i}^{n})$
(4.33)
$\displaystyle\quad-\sum_{i=1}^{2M}H_{\eta}^{\prime}(\eta_{i}^{n+1})(\eta_{i}^{n+1}-\eta_{i}^{n})-\sum_{i=1}^{2M}H_{\sigma}^{\prime}(\sigma_{i}^{n+1})(\sigma_{i}^{n+1}-\sigma_{i}^{n})$
(4.34)
$\displaystyle\quad-\sum_{i=1}^{2M}(\rho_{i}^{n+1}-\rho_{i}^{n})(W_{\rho}\ast\rho^{**})_{i}-\sum_{i=1}^{2M}(\rho_{i}^{n+1}-\rho_{i}^{n})(W_{\sigma}\ast\eta^{**})_{i}$
(4.35)
$\displaystyle\quad-\sum_{i=1}^{2M}(\eta_{i}^{n+1}-\eta_{i}^{n})(W_{\eta}\ast\eta^{**})_{i}-\sum_{i=1}^{2M}(\eta_{i}^{n+1}-\eta_{i}^{n})(W_{\sigma}\ast\rho^{**})_{i},$
(4.36)

where $\boldsymbol{V}_{i}=(V_{\rho,i},V_{\eta,i})^{\top}$. On the other hand,
the evolution of the discrete energy,
$E_{\Delta}(\boldsymbol{\rho}^{n+1})-E_{\Delta}(\boldsymbol{\rho}^{n})$, is
equal to the sum of
$\mathcal{H}_{\Delta}(\boldsymbol{\rho}^{n+1})-\mathcal{H}_{\Delta}(\boldsymbol{\rho}^{n})$,
$\mathcal{V}_{\Delta}(\boldsymbol{\rho}^{n+1})-\mathcal{V}_{\Delta}(\boldsymbol{\rho}^{n})$,
and
$\mathcal{W}_{\Delta}(\boldsymbol{\rho}^{n+1})-\mathcal{W}_{\Delta}(\boldsymbol{\rho}^{n})$.
The contribution corresponding to $\mathcal{H}$ is

$\displaystyle\mathcal{H}_{\Delta}(\boldsymbol{\rho}^{n+1})-\mathcal{H}_{\Delta}(\boldsymbol{\rho}^{n})$
(4.37)
$\displaystyle=\sum_{i=1}^{2M}\left[H_{\rho}(\rho_{i}^{n+1})-H_{\rho}(\rho_{i}^{n})+H_{\eta}(\eta_{i}^{n+1})-H_{\eta}(\eta_{i}^{n})+H_{\sigma}(\sigma_{i}^{n+1})-H_{\sigma}(\sigma_{i}^{n})\right]\Delta
x;$ (4.38)

the contribution of $\mathcal{V}$ is

$\displaystyle\mathcal{V}_{\Delta}(\boldsymbol{\rho}^{n+1})-\mathcal{V}_{\Delta}(\boldsymbol{\rho}^{n})$
$\displaystyle=\sum_{i=1}^{2M}\left[V_{\rho,i}(\rho_{i}^{n+1}-\rho_{i}^{n})+V_{\eta,i}(\eta_{i}^{n+1}-\eta_{i}^{n})\right]\Delta
x$ (4.39)
$\displaystyle=\sum_{i=1}^{2M}\boldsymbol{V}_{i}\cdot(\boldsymbol{\rho}_{i}^{n+1}-\boldsymbol{\rho}_{i}^{n});$
(4.40)

finally, the contribution of $\mathcal{W}$ is

$\displaystyle\mathcal{W}_{\Delta}(\boldsymbol{\rho}^{n+1})-\mathcal{W}_{\Delta}(\boldsymbol{\rho}^{n})$
$\displaystyle=\frac{1}{2}\sum_{i=1}^{2M}\sum_{k=1}^{2M}\left(\rho_{i}^{n+1}\rho_{k}^{n+1}-\rho_{i}^{n}\rho_{k}^{n}\right)W_{\rho,i-k}\Delta
x^{2}$ (4.41)
$\displaystyle\quad+\frac{1}{2}\sum_{i=1}^{2M}\sum_{k=1}^{2M}\left(\eta_{i}^{n+1}\eta_{k}^{n+1}-\eta_{i}^{n}\eta_{k}^{n}\right)W_{\eta,i-k}\Delta
x^{2}$ (4.42)
$\displaystyle\quad+\sum_{i=1}^{2M}\sum_{k=1}^{2M}\left(\rho_{i}^{n+1}\eta_{k}^{n+1}-\rho_{i}^{n}\eta_{k}^{n}\right)W_{\sigma,i-k}\Delta
x^{2}.$ (4.43)

The second contribution, that pertaining to $\mathcal{V}$, corresponds to the
identity Eq. 4.32. Through its substitution, we may rewrite the evolution of
the discrete energy as the sum of three terms:

$\displaystyle
E_{\Delta}[\boldsymbol{\rho}^{n+1}]-E_{\Delta}[\boldsymbol{\rho}^{n}]=I+II+III.$
(4.44)

The first energy term involves only entropies:

$\displaystyle I$
$\displaystyle=\sum_{i=1}^{2M}\left[H_{\rho}(\rho_{i}^{n+1})-H_{\rho}(\rho_{i}^{n})-H_{\rho}^{\prime}(\rho_{i}^{n+1})(\rho_{i}^{n+1}-\rho_{i}^{n})\right]\Delta
x$ (4.45)
$\displaystyle\quad+\sum_{i=1}^{2M}\left[H_{\eta}(\eta_{i}^{n+1})-H_{\eta}(\eta_{i}^{n})-H_{\eta}^{\prime}(\eta_{i}^{n+1})(\eta_{i}^{n+1}-\eta_{i}^{n})\right]\Delta
x$ (4.46)
$\displaystyle\quad+\sum_{i=1}^{2M}\left[H_{\sigma}(\sigma_{i}^{n+1})-H_{\sigma}(\sigma_{i}^{n})-H_{\sigma}^{\prime}(\sigma_{i}^{n+1})(\sigma_{i}^{n+1}-\sigma_{i}^{n})\right]\Delta
x.$ (4.47)

This quantity is immediately controlled, using the convexity of $H_{\rho}$,
$H_{\eta}$, and $H_{\sigma}$; for instance,
$H_{\rho}(\rho_{i}^{n+1})-H_{\rho}(\rho_{i}^{n})-H_{\rho}^{\prime}(\rho_{i}^{n+1})(\rho_{i}^{n+1}-\rho_{i}^{n})\leq
0$. Applying this to every term, we find $I\leq 0$.

The second energy term involves only potentials:

$\displaystyle II$
$\displaystyle=\frac{1}{2}\sum_{i=1}^{2M}\sum_{k=1}^{2M}\left(\rho_{i}^{n+1}\rho_{k}^{n+1}-\rho_{i}^{n}\rho_{k}^{n}-(\rho_{i}^{n+1}-\rho_{i}^{n})(\rho_{k}^{n+1}+\rho_{k}^{n})\right)W_{\rho,i-k}\Delta
x^{2}$ (4.48)
$\displaystyle\quad+\frac{1}{2}\sum_{i=1}^{2M}\sum_{k=1}^{2M}\left(\eta_{i}^{n+1}\eta_{k}^{n+1}-\eta_{i}^{n}\eta_{k}^{n}-(\eta_{i}^{n+1}-\eta_{i}^{n})(\eta_{k}^{n+1}+\eta_{k}^{n})\right)W_{\eta,i-k}\Delta
x^{2}$ (4.49)
$\displaystyle\quad+\frac{1}{2}\sum_{i=1}^{2M}\sum_{k=1}^{2M}\left(\rho_{i}^{n+1}\eta_{k}^{n+1}-\rho_{i}^{n}\eta_{k}^{n}-(\rho_{i}^{n+1}-\rho_{i}^{n})(\eta_{k}^{n+1}+\eta_{k}^{n})\right)W_{\sigma,i-k}\Delta
x^{2}$ (4.50)
$\displaystyle\quad+\frac{1}{2}\sum_{i=1}^{2M}\sum_{k=1}^{2M}\left(\rho_{i}^{n+1}\eta_{k}^{n+1}-\rho_{i}^{n}\eta_{k}^{n}-(\eta_{i}^{n+1}-\eta_{i}^{n})(\rho_{k}^{n+1}+\rho_{k}^{n})\right)W_{\sigma,i-k}\Delta
x^{2}.$ (4.51)

Exploiting the symmetric character of the potentials, e.g.
$W_{\sigma,i-k}=W_{\sigma,k-i}$, this term is shown to be identically zero.

The last energy term involves only fluxes:

$\displaystyle III$ $\displaystyle=-\Delta
t\sum_{i=1}^{2M}\boldsymbol{\xi}_{i}^{n+1}\cdot(\boldsymbol{F}_{i+1/2}^{n+1}-\boldsymbol{F}_{i-1/2}^{n+1})$
(4.52) $\displaystyle=\Delta
t\sum_{i=1}^{2M-1}(\boldsymbol{\xi}_{i+1}^{n+1}-\boldsymbol{\xi}_{i}^{n+1})\cdot\boldsymbol{F}_{i+1/2}^{n+1},$
(4.53)

performing summation by parts and invoking the zero-flux boundary conditions.
Applying the definition of the numerical flux, each summand is seen to equal

$\displaystyle(\boldsymbol{\xi}_{i+1}^{n+1}-\boldsymbol{\xi}_{i}^{n+1})\cdot\left(\operatorname{diag}({\boldsymbol{\rho}_{i}^{n+1}})(\psi_{i+1}^{n+1})^{+}(\boldsymbol{u}_{i+1/2}^{n+1})^{+}+\operatorname{diag}({\boldsymbol{\rho}_{i+1}^{n+1}})(\psi_{i}^{n+1})^{+}(\boldsymbol{u}_{i+1/2}^{n+1})^{-}\right).$
(4.54)

Assuming this quantity is not zero, the term $\boldsymbol{u}_{i+1/2}^{n+1}$
will either be positive or negative. If it is positive, the summand equals

$\displaystyle-(\boldsymbol{\xi}_{i+1}^{n+1}-\boldsymbol{\xi}_{i}^{n+1})\cdot\operatorname{diag}({\boldsymbol{\rho}_{i}^{n+1}})(\psi_{i+1}^{n+1})^{+}S_{i+1/2}^{n+1}\frac{\boldsymbol{\xi}_{i+1}^{n+1}-\boldsymbol{\xi}_{i}^{n+1}}{\Delta
x},$ (4.55)

which is found to be non-positive, since the matrix
$\operatorname{diag}({\boldsymbol{\rho}_{i}^{n+1}})(\psi_{i+1}^{n+1})^{+}S_{i+1/2}^{n+1}$
is positive semi-definite. If, on the other hand,
$\boldsymbol{u}_{i+1/2}^{n+1}$ is negative, the summand is again non-positive,
this time using the definiteness of the matrix
$\operatorname{diag}({\boldsymbol{\rho}_{i+1}^{n+1}})(\psi_{i}^{n+1})^{+}S_{i+1/2}^{n+1}$.
All in all, the sum is proven to be non-positive, finally yielding

$\displaystyle
E_{\Delta}[\boldsymbol{\rho}^{n+1}]-E_{\Delta}[\boldsymbol{\rho}^{n}]=I+II+III\leq
III\leq 0.$ (4.56)

∎

######  Remark 4.3 (Schemes for gradient flows in higher dimensions).

The generalisation of scheme (4.26) to higher dimensions is non-trivial. One
option is to prescribe a full scheme in higher dimensions; the proofs of
positivity, boundedness, and energy dissipation all generalise nicely, but the
computational cost of an implicit numerical scheme with non-local terms in
multiple dimensions is impractically high. A direct dimensional splitting
generalisation will capture the positivity and boundedness properties, but
will not preserve the energy dissipation whenever the non-local terms are
present.

In the non-local case, we will trivially adapt the sweeping dimensional
splitting technique of [2], which preserves the dissipative structure on top
of the other properties while remaining computationally efficient. This
technique is employed in the simulations of Section 5.4.

##  5 Numerical Experiments

This section is devoted to showcasing the properties of our numerical schemes
in several challenging problems. We will use the implicit scheme (4.26) for
all examples, except for Section 5.2, where we compare the performance of both
the implicit scheme (4.12) and its explicit version. This choice is based on
computational performance; the implicit schemes consistently outperform their
explicit counterparts, since they admit much larger time steps while
preserving all the properties of interest. We shall prioritise numerical
experiments involving systems of equations, as these are both more challenging
and more interesting.

We first verify the numerical accuracy of our implicit scheme on systems with
and without saturation in Section 5.1. Section 5.2 describes the effect of
saturation in a gradient flow with an external potential, a case not directly
covered by the theory of gradient flows with non-linear mobility [11]. We
numerically explore the qualitative properties of the solutions, showing that
saturation leads to the convergence towards stationary states with kinks at
the free boundary of the fully saturated region. Section 5.3 is devoted to
another novel qualitative effect for systems due to saturation, what we call
the freeze-in-place phenomenon. Because of the upper bound on the density,
certain populations, whose fate was to segregate in the absence of saturation
effects, do not achieve complete segregation, as the total population reaches
the saturation level, leading to a “frozen” configuration. Finally, Section
5.4 showcases cell-cell adhesion sorting phenomena, based on the different
strengths of attraction between two subpopulations. The different attractions
between the species lead to segregation, partial engulfment, or complete
engulfment, as was reported in [14] for compactly supported attractive
kernels.

###  5.1 Cross-Diffusion Experiment

Figure 1:  SKT model. ${L^{1}}$, ${L^{2}}$, and ${L^{\infty}}$ errors of the
numerical solution of scheme (4.26) with respect to the manufactured solution
(5.5). Left: regular model (5.2). Right: model with saturation (5.10). $\Delta
x=2^{-k}\pi$ and $\Delta t=2^{-k}/10$ for $k=1,\cdots,8$.

The first experiment aims to verify the numerical accuracy of the scheme
(4.26) on a problem with cross diffusion, a feature not present in the scalar
tests found in [2]. To that end, we shall employ the SKT population model
[24], given by

$\left\\{\begin{array}[]{{>{\displaystyle\vphantom{\frac{1}{1}}}l}}\partial_{t}\rho=\nabla\cdot\left[\left(2\rho+\eta\right)\nabla\rho+\rho\nabla\eta\vphantom{{x^{2}}^{2}}\right]+s_{\rho},\\\
\partial_{t}\eta=\nabla\cdot\left[\eta\nabla\rho+\left(\rho+2\eta\right)\nabla\eta\vphantom{{x^{2}}^{2}}\right]+s_{\eta},\end{array}\right.$
(5.1)

where $\boldsymbol{s}=(s_{\rho},s_{\eta})^{\top}$ is a given source term. We
will consider the problem on $\Omega=(-\pi,\pi)^{2}$, posed with periodic
boundary conditions. Equation 5.1 can be cast in an analogous form to (4.16)
with right-hand sides, written as

$\partial_{t}\begin{pmatrix}\rho\\\
\eta\end{pmatrix}=\nabla\cdot\left[M(\rho,\eta)\nabla\begin{pmatrix}\log\rho\\\
\log\eta\end{pmatrix}\right]+\begin{pmatrix}s_{\rho}\\\
s_{\eta}\end{pmatrix},$ (5.2)

where

$\displaystyle M(\rho,\eta)=\begin{pmatrix}\rho(2\rho+\eta)&\rho\eta\\\
\rho\eta&(\rho+2\eta)\eta\end{pmatrix}=\begin{pmatrix}\rho&0\\\
0&\eta\end{pmatrix}\begin{pmatrix}2\rho+\eta&\eta\\\
\rho&\rho+2\eta\end{pmatrix}.$ (5.3)

If $\boldsymbol{s}=0$, then this is a gradient flow with respect to the energy

$\displaystyle
E=\int_{\Omega}\rho\left(\log\rho-1\right)+\eta\left(\log\eta-1\right)\mathop{}\\!\mathrm{d}\boldsymbol{x}.$
(5.4)

In order to perform the validation, we will use manufactured solutions. We
assume the solution is of the form

$\displaystyle\rho(t,x,y)=\frac{1}{4}\left[1+\sin(x+t)\right],\quad\eta(t,x,y)=\frac{1}{4}\left[1+\cos(y+t)\right],$
(5.5)

and compute the corresponding source terms,

$\displaystyle s_{\rho}$
$\displaystyle=\frac{1}{4}\cos\left(t+x\right)+\frac{3}{16}\sin\left(t+x\right)+\frac{1}{16}\cos\left(t+y\right)$
(5.6)
$\displaystyle\quad-\frac{1}{8}\cos^{2}\left(t+x\right)+\frac{1}{8}\sin^{2}\left(t+x\right)+\frac{1}{8}\sin\left(t+x\right)\cos\left(t+y\right),$
(5.7)

and

$\displaystyle s_{\eta}$
$\displaystyle=\frac{1}{16}\sin\left(t+x\right)+\frac{3}{16}\cos\left(t+y\right)-\frac{1}{4}\sin\left(t+y\right)$
(5.8)
$\displaystyle\quad+\frac{1}{8}\cos^{2}\left(t+y\right)-\frac{1}{8}\sin^{2}\left(t+y\right)+\frac{1}{8}\sin\left(t+x\right)\cos\left(t+y\right).$
(5.9)

We will also perform the validation on a version of Equation 5.2 with a
saturation term:

$\partial_{t}\begin{pmatrix}\rho\\\
\eta\end{pmatrix}=\nabla\cdot\left[\left(1-\rho-\eta\right)M(\rho,\eta)\nabla\begin{pmatrix}\log\rho\\\
\log\eta\end{pmatrix}\right]+\begin{pmatrix}s_{\rho}\\\
s_{\eta}\end{pmatrix}.$ (5.10)

Assuming the same solution, (5.5), yields the sources

$\displaystyle s_{\rho}$
$\displaystyle=\frac{1}{4}\cos\left(t+x\right)+\frac{3}{32}\sin\left(t+x\right)+\frac{1}{32}\cos\left(t+y\right)-\frac{1}{64}\cos^{2}\left(t+x\right)$
(5.11)
$\displaystyle\quad+\frac{1}{64}\sin^{2}\left(t+x\right)-\frac{1}{32}\sin^{3}\left(t+x\right)-\frac{1}{64}\cos^{2}\left(t+y\right)+\frac{1}{64}\sin^{2}\left(t+y\right)$
(5.12)
$\displaystyle\quad+\frac{1}{16}\cos^{2}\left(t+x\right)\sin\left(t+x\right)-\frac{1}{32}\cos^{2}\left(t+y\right)\sin\left(t+x\right)+\frac{1}{64}\sin^{2}\left(t+y\right)\sin\left(t+x\right)$
(5.13)
$\displaystyle\quad+\frac{3}{64}\cos^{2}\left(t+x\right)\cos\left(t+y\right)-\frac{1}{16}\sin^{2}\left(t+x\right)\cos\left(t+y\right),$
(5.14)

and

$\displaystyle s_{\eta}$
$\displaystyle=\frac{1}{32}\sin\left(t+x\right)+\frac{3}{32}\cos\left(t+y\right)-\frac{1}{4}\sin\left(t+y\right)+\frac{1}{64}\cos^{2}\left(t+x\right)$
(5.15)
$\displaystyle\quad-\frac{1}{64}\sin^{2}\left(t+x\right)+\frac{1}{64}\cos^{2}\left(t+y\right)-\frac{1}{32}\cos^{3}\left(t+y\right)-\frac{1}{64}\sin^{2}\left(t+y\right)$
(5.16)
$\displaystyle\quad-\frac{1}{16}\cos^{2}\left(t+y\right)\sin\left(t+x\right)+\frac{3}{64}\sin^{2}\left(t+y\right)\sin\left(t+x\right)+\frac{1}{64}\cos^{2}\left(t+x\right)\cos\left(t+y\right)$
(5.17)
$\displaystyle\quad-\frac{1}{32}\sin^{2}\left(t+x\right)\cos\left(t+y\right)+\frac{1}{16}\sin^{2}\left(t+y\right)\cos\left(t+y\right).$
(5.18)

We solve Eqs. 5.2 and 5.10 over the interval $t\in(0,0.1)$, choosing $\Delta
x=\Delta y=2^{-k}\pi$ and $\Delta t=2^{-k}/10$ for $k=1,\cdots,8$, and compute
the ${L^{p}}$ error ($p=1,2,\infty$) at the final time with respect to the
solution (5.5). Figure 1 shows the results of the validation; scheme (4.26)
clearly approximates the manufactured solution (5.5) with first-order
accuracy.

###  5.2 Saturation Experiment

Figure 2:  Stationary solutions to Eq. 5.19 in one dimension. In the cases
$M<M_{c}$ and $M=M_{c}$, the solution is a half-Gaussian. For $M>M_{c}$, the
solution is made of a straight saturated segment and a Gaussian tail.  Figure
3:  Stationary solution of Eq. 5.19 in one dimension with mass $M=1.66$
corresponding to $l=1$. $\Omega=(0,4)$, $\alpha=1$, $D=1$, $C=1$. Explicit:
$\Delta t=\Delta x^{2}/4$. Implicit: $\Delta t=\Delta x$.  Figure 4:
Stationary solution of Eq. 5.19 in two dimensions with mass $M=4.71$
corresponding to $l=2$. $\Omega=(0,4)^{2}$, $\alpha=1$, $D=1$, $C=1$, $\Delta
t=\Delta x=2^{-7}$.

The second experiment concerns the saturation effects. In order to demonstrate
the boundedness property of the numerical schemes, we consider the equation

$\displaystyle\partial_{t}\rho=\nabla\cdot\left(\rho(\alpha-\rho)\nabla\left(D\log(\rho)+\frac{C}{2}|\boldsymbol{x}|^{2}\right)\right).$
(5.19)

This problem exhibits different steady states which depend on
$M\coloneqq\|\rho_{0}\|_{{L^{1}}}$, the (conserved) mass of the solution. If
$M$ is small ($M<M_{c}$ for some critical mass $M_{c}$), the steady state is
simply

$\displaystyle\rho_{\infty}=A\exp\left(-\frac{C}{2D}|\boldsymbol{x}|^{2}\right),$
(5.20)

where $A$ is a positive constant such that
$\|\rho_{\infty}\|_{{L^{1}}}=M=\|\rho_{0}\|_{{L^{1}}}$. However, if the mass
is beyond the critical value, the steady state consists of two segments:

$\displaystyle\rho_{\infty}=\begin{cases}\alpha&\text{if }|\boldsymbol{x}|\leq
l,\\\ \displaystyle
B\exp\left(-\frac{C}{2D}|\boldsymbol{x}|^{2}\right)&\text{if
}|\boldsymbol{x}|\geq l,\end{cases}$ (5.21)

where $B$ is again a normalisation constant, and $l$ is to be determined, see
Fig. 2 for a diagram. Observing that $\rho_{\infty}=\alpha$ whenever
$|\boldsymbol{x}|=l$, we determine
$B=\alpha\exp\left(\frac{Cl^{2}}{2D}\right)$; thus the steady state is
rewritten as

$\displaystyle\rho_{\infty}(\boldsymbol{x})=\alpha\exp\left(-\frac{C}{2D}(|\boldsymbol{x}|^{2}-l^{2})^{+}\right),$
(5.22)

where $(s)^{+}=\max\\{s,0\\}$ for any $s\in{\mathbb{R}}$. The value of $l$ is
an increasing function of $M$, and can be determined form the initial datum.
It is left to the reader to check (5.22) are indeed weak stationary solutions
to (5.19).

We may further explore the normalization condition. We first consider Eq. 5.19
in one dimension, posed on $\Omega=(0,\infty)$; by virtue of the conservation
of mass, we arrive at the self-consistency equation

$\displaystyle M$
$\displaystyle=\int_{0}^{\infty}\alpha\exp\left(-\frac{C}{2D}(x^{2}-l^{2})^{+}\right)\mathop{}\\!\mathrm{d}x=\alpha\left[l+\int_{l}^{\infty}\exp\left(\frac{C}{2D}(l^{2}-x^{2})\right)\mathop{}\\!\mathrm{d}x\right]$
(5.23) $\displaystyle=\alpha\left[l+\sqrt{\frac{\pi
D}{2C}}\exp\left(\frac{C}{2D}l^{2}\right)\left[1-\operatorname{erf}\left(l\sqrt{\frac{C}{2D}}\right)\right]\right].$
(5.24)

In particular, this expression yields the value of the critical mass,
$M_{c}=\alpha\sqrt{\frac{\pi D}{2C}}$, which corresponds to $l=0$. In the two
dimensional case, we pose Eq. 5.19 on the domain $\Omega=(0,\infty)^{2}$. We
find

$\displaystyle M$
$\displaystyle=\int_{0}^{\infty}\\!\\!\\!\int_{0}^{\infty}\\!\\!\alpha\exp\left(-\frac{C}{2D}(x^{2}+y^{2}-l^{2})^{+}\right)\mathop{}\\!\mathrm{d}x\mathop{}\\!\mathrm{d}y=\int_{0}^{\frac{\pi}{2}}\\!\\!\\!\int_{0}^{\infty}\\!\\!\alpha
r\exp\left(-\frac{C}{2D}(r^{2}-l^{2})^{+}\right)\mathop{}\\!\mathrm{d}r\mathop{}\\!\mathrm{d}\theta$
(5.25)
$\displaystyle=\frac{\alpha\pi}{2}\left[\int_{0}^{l}r\mathop{}\\!\mathrm{d}r+\int_{l}^{\infty}r\exp\left(-\frac{C}{2D}(r^{2}-l^{2})\right)\mathop{}\\!\mathrm{d}r\right]=\frac{\alpha\pi}{2}\left[\frac{l^{2}}{2}+\frac{D}{C}\right],$
(5.26)

yielding $M_{c}=\frac{\alpha\pi D}{2C}.$

We will solve Eq. 5.19 in one and two dimensions, over sufficiently large but
finite domains, with supercritical masses, in order to validate the behaviour
of the schemes in the presence of saturation. The problems are initialised
with a constant density $\rho=M/|\Omega|$, where $|\Omega|$ is the Lebesgue
measure of $\Omega$; the solution is computed for $t\in(0,15)$, and the final
state is taken as an approximation of the steady state.

In one dimension, we let $\Omega=(0,4)$ and choose the supercritical mass
$M=1.66$, which corresponds to $l=1$. Figure 3 shows the numerical steady
state, computed with both the explicit and implicit schemes, each on a coarse
($\Delta x=2^{-3}$) and a fine ($\Delta x=2^{-7}$) mesh. Both schemes capture
the upper bound of the solution, but the implicit scheme permits the use of a
much larger time step. In two dimensions, we set $\Omega=(0,4)^{2}$ and choose
the supercritical mass $M=4.71$, found by letting $l=2$. Figure 4 shows the
numerical steady state computed with the implicit scheme.

We would like to emphasise that there is no general weak solution theory yet
that establishes the kink solutions (5.22) as the global asymptotic profile
for generic initial data to (5.19) with supercritical mass, though this
behaviour was observed numerically for diverse initial data. The Cauchy theory
for these general gradient flows (5.19) with non-linear mobility and external
interaction potentials, not included in [11], is challenging and will be
studied elsewhere.

###  5.3 Freeze-In-Place Experiment

Figure 5:  Stationary solutions of Eq. 5.29 in one dimension. $\Omega=(-1,1)$,
$\alpha=1$, $D=0.1$, $C_{1}=4$, $C_{2}=2$, $\Delta x=2^{-8}$, and $\Delta
t=0.1$. Top: datum (5.32). Middle: “frozen-in-place” stationary state with
saturation effects. Bottom: segregated stationary state without saturation.
Figure 6:  Stationary solutions of Eq. 5.29 in two dimensions.
$\Omega=(-1,1)^{2}$, $\alpha=1$, $D=0.1$, $C_{1}=4$, $C_{2}=2$, $\Delta
x=2^{-7}$, and $\Delta t=0.1$. Top: datum (5.32). Middle: “frozen-in-place”
stationary state with saturation effects. Bottom: segregated stationary state
without saturation.

This test will explore the appearance of infinitely many steady states as a
result of the saturation effects. To that end, we will consider the two
species system

$\displaystyle\left\\{\begin{array}[]{{>{\displaystyle\vphantom{\frac{1}{1}}}l}}\partial_{t}\rho=\nabla\cdot\left(\rho(\alpha-\sigma)\nabla\left(D\sigma+\frac{C_{1}}{2}|\boldsymbol{x}|^{2}\right)\right),\\\\[10.00002pt]
\partial_{t}\eta=\nabla\cdot\left(\eta(\alpha-\sigma)\nabla\left(D\sigma+\frac{C_{2}}{2}|\boldsymbol{x}|^{2}\right)\right).\end{array}\right.$
(5.29)

There are two effects present here: diffusion and confinement. The non-linear
diffusion acts on the sum, $\sigma=\rho+\eta$. This can be seen by ignoring
the $|\boldsymbol{x}|^{2}$ terms; then, the two equations sum to
$\partial_{t}\sigma=D\nabla\cdot\left(\sigma(\alpha-\sigma)\nabla\sigma\right)$.
Hence, the diffusion terms act to smooth irregularities on $\sigma$, but not
necessarily on the individual profiles. Meanwhile, the quadratic potentials
serve to confine each species, driving their densities towards the origin.

Crucially, we will choose $C_{1}=2C_{2}$; the confinement effects will
therefore act more strongly on $\rho$ than on $\eta$. This difference in
strengths, combined with the diffusion, will cause $\rho$ to displace $\eta$
near the origin. Below the saturation level, the species are expected to
segregate completely; however, this might not occur if the species sum reaches
the saturation level. As $\sigma$ approaches $\alpha$, the mobility reduces
drastically, which will cause the solution to “freeze” in whatever
configuration it may find itself.

We will let $D=0.1$, $C_{1}=4$, $C_{2}=2$, and $\alpha=1$. In one dimension,
we solve the problem on $\Omega=(-1,1)$, with the initial datum

$\displaystyle\left\\{\begin{array}[]{{>{\displaystyle\vphantom{\frac{1}{1}}}l}}\rho_{1}(0,x)=\left[f(x)\left(1-\frac{\cos(\omega
x)}{2}\right)\right]^{+},\\\\[10.00002pt]
\rho_{2}(0,x)=\left[f(x)\left(1+\frac{\cos(\omega
x)}{2}\right)\right]^{+},\end{array}\right.\qquad\text{where
}f(x)=\frac{4}{5}\left(1-\left(\frac{4x}{3}\right)^{2}\right),$ (5.32)

and $\omega=16\pi$. We will solve the system for $t\in(0,30)$ and take the
final state as an approximation to the asymptotic steady state which
corresponds to this datum. Figure 5 shows the datum, the “frozen-in-place”
steady state, as well as the segregated steady state of the analogous system
without saturation, all computed with the implicit scheme (4.26).

In two dimensions, we solve the problem on $\Omega=(-1,1)^{2}$. The datum now
reads

$\displaystyle\left\\{\begin{array}[]{{>{\displaystyle\vphantom{\frac{1}{1}}}l}}\rho_{1}(0,x,y)=\left[f\left(\sqrt{x^{2}+y^{2}}\right)\left(1-\frac{\cos(\omega
x)\cos(\omega y)}{2}\right)\right]^{+},\\\\[10.00002pt]
\rho_{2}(0,x,y)=\left[f\left(\sqrt{x^{2}+y^{2}}\right)\left(1+\frac{\cos(\omega
x)\cos(\omega y)}{2}\right)\right]^{+},\end{array}\right.$ (5.35)

where $\omega=4\pi$ and $f$ is defined as above. Once again, we take the
solution at $t=30$ as an approximation to the asymptotic steady state. Figure
6 shows the datum as well as the steady states with and without saturation. As
in the one-dimensional case, the saturation induces a “freezing-in-place”
effect, preventing the complete segregation of the species.

###  5.4 Cell-Cell Adhesion Experiment

Figure 7:  Steady state of Eqs. 5.36 and 5.37 in two dimensions.
$\Omega=(-2,2)^{2}$, $\alpha=1$, $\varepsilon=0.1$,
$c_{\rho\rho}=2c_{\rho\eta}=2c_{\eta\rho}=c_{\eta\eta}=1$, $\Delta x=\Delta
y=\Delta t=2^{-6}$. Left: Eq. 5.37 (saturation). Right: Eq. 5.36 (no
saturation).  Figure 8:  Steady state of Eqs. 5.36 and 5.37 in two dimensions.
$\Omega=(-2,2)^{2}$, $\alpha=1$, $\varepsilon=0.1$,
$4c_{\rho\rho}=2c_{\rho\eta}=2c_{\eta\rho}=c_{\eta\eta}=1$, $\Delta x=\Delta
y=\Delta t=2^{-6}$. Left: Eq. 5.37 (saturation). Right: Eq. 5.36 (no
saturation).  Figure 9:  Dissipation of the free energy in the solutions to
Eqs. 5.36 and 5.37 shown in Figs. 7 and 8. Left: partial engulfment,
corresponding to Fig. 7. Right: complete engulfment, corresponding to Fig. 8.

To conclude the experiments section, we explore a version of the two-species
model for cell-cell adhesion introduced in [14]. The model for the evolution
of the densities of two different species, $\rho$ and $\eta$, is given by

$\left\\{\begin{array}[]{{>{\displaystyle\vphantom{\frac{1}{1}}}l}}\partial_{t}\rho=\nabla\cdot\left(\rho\nabla(\varepsilon\sigma+c_{\rho\rho}K\ast\rho+c_{\rho\eta}K\ast\eta)\right),\\\
\partial_{t}\eta=\nabla\cdot\left(\eta\nabla(\varepsilon\sigma+c_{\eta\rho}K\ast\rho+c_{\eta\eta}K\ast\eta)\right),\end{array}\right.$
(5.36)

in two spatial dimensions, for a positive value of $\varepsilon$, non-negative
values of $c_{\rho\rho}$, $c_{\rho\eta}$, $c_{\eta\rho}$, and $c_{\eta\eta}$,
and where $K(\boldsymbol{x})=|\boldsymbol{x}|^{2}/2$. The parameter
$\varepsilon$ regulates the strength of non-linear diffusion on the total
density of the system, which in this model accounts for localised repulsion
effects. Each of the coefficients $c_{\alpha\beta}$ correspond to the strength
of attraction of the species $\alpha$ to species $\beta$.

This model exhibits different asymptotic behaviours, regulated by the values
of $c_{\alpha\beta}$. As discussed in [14], different values of the parameters
may lead to the mixing of the two species, a complete segregation, or the
engulfment of one of the species (either partial or complete). We will
concentrate here on the latter case, and will observe how the inclusion of
saturation effects affect these behaviours, through the modified model

$\left\\{\begin{array}[]{{>{\displaystyle\vphantom{\frac{1}{1}}}l}}\partial_{t}\rho=\nabla\cdot\left(\rho(\alpha-\sigma)\nabla(\varepsilon\sigma+c_{\rho\rho}K\ast\rho+c_{\rho\eta}K\ast\eta)\right),\\\
\partial_{t}\eta=\nabla\cdot\left(\eta(\alpha-\sigma)\nabla(\varepsilon\sigma+c_{\eta\rho}K\ast\rho+c_{\eta\eta}K\ast\eta)\right).\end{array}\right.$
(5.37)

We will prescribe an initial datum consisting of two touching disks:

$\displaystyle\left\\{\begin{array}[]{{>{\displaystyle\vphantom{\frac{1}{1}}}l}}\rho_{0}=0.95\,\mathds{1}_{\mathcal{D}_{\rho}},\\\
\eta_{0}=0.95\,\mathds{1}_{\mathcal{D}_{\eta}},\end{array}\right.$ (5.40)

where $\mathcal{D}_{\rho}$ is a disk centered at $(-0.5,0)$ with radius $0.5$,
and $\mathcal{D}_{\rho}$ is a disk centred at $(-0.4,0)$ with radius $0.4$. We
will let $\varepsilon=0.1$, and $\alpha=1$ when applicable. We shall solve
Eqs. 5.36 and 5.37 on $\Omega=(-2,2)^{2}$, over the interval $t\in(0,45)$, and
take their final state as an approximation of the asymptotic steady state. We
choose $\Delta x=\Delta y=\Delta t=2^{-6}$ to ensure high accuracy.

We first explore the partial engulfment case, given by
$c_{\rho\rho}=2c_{\rho\eta}=2c_{\eta\rho}=c_{\eta\eta}=1$, where the two
species are expected to partially surround each other in a lunar shape. Figure
7 shows the numerical steady states. In the system with saturation, the steady
state exhibits a saturated, flat region. Nevertheless, both solutions are
relatively similar.

We turn to the complete engulfment case, given by
$4c_{\rho\rho}=2c_{\rho\eta}=2c_{\eta\rho}=c_{\eta\eta}=1$. Here, one species
is expected to completely surround the other; in the long-time limit, one
species will be supported on a disk, and the other, on an annulus. Figure 8
shows the numerical steady states. This time, the effects of the saturation
are very visible: the initial concentration of the density creates a region of
saturation where the two species remain mixed for all time, in the same vein
as the example from Section 5.3. This prevents the full segregation seen in
the unsaturated case, and yields a cardioid-like free boundary between the
species.

We conclude by showing the evolution of the (relative) free energy
corresponding to the both the partial engulfment and the complete engulfment
scenarios in Fig. 9. The energy is clearly dissipated by the numerical scheme,
and the dissipation rates are affected by the presence of the saturation
terms. In the partial engulfment case, the energies with and without
saturation evolve very differently, but eventually reach a comparable value;
this is consistent with the similar (though not equal) steady states reached
in either setting. Meanwhile, the energies of the complete engulfment case
demonstrate similar evolutions, but the dissipation rate is higher in the
unsaturated setting.

## Acknowledgements

RB was funded by Labex CEMPI (ANR-11-LABX-0007-01). RB and JAC were supported
by the Advanced Grant Nonlocal-CPD (Nonlocal PDEs for Complex Particle
Dynamics: Phase Transitions, Patterns and Synchronization) of the European
Research Council Executive Agency (ERC) under the European Union’s Horizon
2020 research and innovation programme (grant agreement No. 883363). JH was
partially supported under the NSF CAREER grant DMS-1654152 and AFOSR grant
FA9550-21-1-0358.

## References

  * [1] L. Almeida, F. Bubba, B. Perthame, and C. Pouchol.  Energy and implicit discretization of the fokker-planck and keller-segel type equations.  Networks Heterog. Media, 14(1):23–41, 2019.
  * [2] R. Bailo, J. A. Carrillo, and J. Hu.  Fully Discrete Positivity-Preserving and Energy-Dissipating Schemes for Aggregation-Diffusion Equations with a Gradient Flow Structure.  Commun. Math. Sci., 18(5):1259–1303, Sept. 2020.
  * [3] R. Bailo, J. A. Carrillo, H. Murakawa, and M. Schmidtchen.  Convergence of a fully discrete and energy-dissipating finite-volume scheme for aggregation-diffusion equations.  Math. Models Methods Appl. Sci., 30(13):2487–2522, Dec. 2020.
  * [4] M. Bruna, M. Burger, H. Ranetbauer, and M.-T. Wolfram.  Cross-Diffusion Systems with Excluded-Volume Effects and Asymptotic Gradient Flow Structures.  J Nonlinear Sci, 27(2):687–719, Apr. 2017.
  * [5] 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.
  * [6] M. Burger, M. Di Francesco, S. Fagioli, and A. Stevens.  Sorting phenomena in a mathematical model for two mutually attracting/repelling species.  SIAM J. Math. Anal., 50(3):3210–3250, 2018.
  * [7] M. Burger and J.-F. Pietschmann.  Flow characteristics in a crowded transport model.  Nonlinearity, 29(11):3528–3550, Oct. 2016.
  * [8] 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.
  * [9] J. A. Carrillo, K. Craig, and Y. Yao.  Aggregation-diffusion equations: dynamics, asymptotics, and singular limits.  In Active particles, Vol. 2, Model. Simul. Sci. Eng. Technol., pages 65–108. Birkhäuser/Springer, Cham, 2019.
  * [10] 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.
  * [11] J. A. Carrillo, S. Lisini, G. Savaré, and D. Slepčev.  Nonlinear mobility continuity equations and generalized displacement convexity.  J. Funct. Anal., 258(4):1273–1309, 2010.
  * [12] J. A. Carrillo, R. J. McCann, and C. Villani.  Kinetic equilibration rates for granular media and related equations: Entropy dissipation and mass transportation estimates.  Rev. Mat. Iberoamericana, 19:971–1018, 2003.
  * [13] J. A. Carrillo, R. J. McCann, and C. Villani.  Contractions in the 2-Wasserstein Length Space and Thermalization of Granular Media.  Arch. Ration. Mech. Anal., 179(2):217–263, Feb. 2006.
  * [14] J. A. Carrillo, H. Murakawa, M. Sato, H. Togashi, and O. Trush.  A population dynamics model of cell-cell adhesion incorporating population pressure and density saturation.  Journal of Theoretical Biology, 474:14–24, Aug. 2019.
  * [15] H. Holden and N. H. Risebro.  Front Tracking for Hyperbolic Conservation Laws, volume 152 of Applied Mathematical Sciences.  Springer Berlin Heidelberg, Berlin, Heidelberg, 2015.
  * [16] H. Holden, N. H. Risebro, and A. Tveito.  Maximum Principles for a Class of Conservation Laws.  SIAM Journal on Applied Mathematics, 55(3):651–661, 1995.
  * [17] T. L. Jackson and H. M. Byrne.  A mechanical model of tumor encapsulation and transcapsular spread.  Mathematical Biosciences, 180(1):307–328, Nov. 2002.
  * [18] A. Jüngel.  The boundedness-by-entropy method for cross-diffusion systems.  Nonlinearity, 28(6):1963–2001, 2015.
  * [19] A. Jüngel and N. Zamponi.  Qualitative behavior of solutions to cross-diffusion systems from population dynamics.  J. Math. Anal. Appl., 440(2):794–809, 2016.
  * [20] A. Jüngel and A. Zurek.  A discrete boundedness-by-entropy method for finite-volume approximations of cross-diffusion systems, 2021.
  * [21] N. Loy and M. Zanella.  Structure preserving schemes for Fokker-Planck equations with nonconstant diffusion matrices.  Math. Comput. Simulation, 188:342–362, 2021.
  * [22] H. Murakawa and H. Togashi.  Continuous models for cell–cell adhesion.  Journal of Theoretical Biology, 374:1–12, 2015.
  * [23] L. Pareschi and M. Zanella.  Structure preserving schemes for nonlinear Fokker-Planck equations and applications.  J. Sci. Comput., 74(3):1575–1600, 2018.
  * [24] N. Shigesada, K. Kawasaki, and E. Teramoto.  Spatial segregation of interacting species.  Journal of Theoretical Biology, 79(1):83–99, July 1979.
  * [25] C.-W. Shu.  Essentially Non-Oscillatory and Weighted Essentially Non-Oscillatory Schemes for Hyperbolic Conservation Laws, volume 1697, pages 325–432.  Springer Berlin Heidelberg, Berlin, Heidelberg, 1998.
  * [26] Z. Sun, J. A. Carrillo, and C.-W. Shu.  A discontinuous Galerkin method for nonlinear parabolic equations and gradient flow problems with interaction potentials.  J. Comput. Phys., 352:76–104, 2018.
  * [27] Z. Sun, J. A. Carrillo, and C.-W. Shu.  An entropy stable high-order discontinuous Galerkin method for cross-diffusion gradient flow systems.  Kinetic & Related Models, 12(4):885–908, 2019.
  * [28] O. Trush, C. Liu, X. Han, Y. Nakai, R. Takayama, H. Murakawa, J. A. Carrillo, H. Takechi, S. Hakeda-Suzuki, T. Suzuki, and M. Sato.  N-cadherin orchestrates self-organization of neurons within a columnar unit in the drosophila medulla.  Journal of Neuroscience, 39(30):5861–5880, 2019.
  * [29] J. L. Vazquez.  The Porous Medium Equation.  Oxford University Press, Oct. 2006.
  * [30] N. Zamponi and A. Jüngel.  Analysis of degenerate cross-diffusion population models with volume filling.  Ann. Inst. H. Poincaré Anal. Non Linéaire, 34(1):1–29, 2017\.