# On multi-species diffusion with size exclusion

Katharina Hopf  Katharina Hopf, Weierstrass Institute for Applied Analysis and
Stochastics (WIAS), Mohrenstrasse 39, 10117 Berlin, Germany, E-Mail:
<EMAIL_ADDRESS>and  Martin Burger  Martin Burger, Department of
Mathematics, Friedrich-Alexander Universität Erlangen-Nürnberg, Cauerstrasse
11, 91058 Erlangen, Germany, E-mail<EMAIL_ADDRESS>
###### Abstract.

We revisit a classical continuum model for the diffusion of multiple species
with size-exclusion constraint, which leads to a degenerate nonlinear cross-
diffusion system. The purpose of this article is twofold: first, it aims at a
systematic study of the question of existence of weak solutions and their
long-time asymptotic behaviour. Second, it provides a weak-strong stability
estimate for a wide range of coefficients, which had been missing so far.

In order to achieve the results mentioned above, we exploit the formal
gradient-flow structure of the model with respect to a logarithmic entropy,
which leads to best estimates in the full-interaction case, where all cross-
diffusion coefficients are non-zero. Those are crucial to obtain the minimal
Sobolev regularity needed for a weak-strong stability result. For meaningful
cases when some of the coefficients vanish, we provide a novel existence
result based on approximation by the full-interaction case.

###### Key words and phrases:

Cross diffusion, size exclusion, volume filling, gradient-flow structure,
degenerate nonlinear mobility, weak-strong stability, large-data weak
solutions, long-time asymptotics, relative entropy, convexity.

###### 2020 Mathematics Subject Classification:

35K59, 35K65, 35A02, 35B35

##  1\. Introduction

Nonlinear cross-diffusion models and their mathematical analysis have received
considerable attention in the last years, including various degenerate cases
as arising in models with size exclusion or local repulsion (cf. e.g. [4, 3,
6, 9, 8, 15, 16, 14, 24]). While the local-in-time existence and uniqueness of
classical solutions has been established a while ago in a more general setting
(cf. [1]), the questions of global-in-time regularity and qualitative
properties of weak solutions continue to pose a significant challenge due to
non-uniform parabolicity and missing maximum principles for systems. From a
modelling point of view such systems can naturally be derived from lattice
systems for multiple species with size exclusion, which found applications in
particular in cell biology (cf. [7, 18, 19, 22]).

In this paper we will revisit the canonical macroscopic cross-diffusion system
arising from these models and derive new results both on the global existence
of weak solutions as well as their stability and uniqueness (in the sense of a
weak-strong stability result). For $T^{*}\in(0,\infty]$, $n\in\mathbb{N}$, a
smooth bounded domain $\Omega\subset\mathbb{R}^{d}$, $d\in\mathbb{N}$, and
constant coefficients $K_{ij}$, $i,j\in\\{0,\dots,n\\}$, we consider the
cross-diffusion system

$\displaystyle\partial_{t}u_{i}-\nabla\cdot\left[\sum_{j=0}^{n}K_{ij}(u_{j}\nabla
u_{i}-u_{i}\nabla u_{j})\right]=0\quad$ $\displaystyle\text{ in
}(0,T^{*})\times\Omega,\quad i\in\\{0,\ldots,n\\},\qquad$ (1.1a) supplemented
by no-flux boundary conditions
$\displaystyle\nu\cdot\sum_{j=0}^{n}K_{ij}(u_{j}\nabla u_{i}-u_{i}\nabla
u_{j})=0\quad$ $\displaystyle\text{ on }(0,T^{*})\times\partial\Omega,\quad
i\in\\{0,\ldots,n\\},\quad$ (1.1b)

with $\nu$ denoting the outward unit normal to $\partial\Omega$. Here,
$u_{i}=u_{i}(t,x)\geq 0$ denotes the density of the $i$-th species $X_{i}$.
Since effectively only off-diagonal entries of $(K_{ij})_{ij}$ enter into
(1.1a), (1.1b), we may assume, without loss of generality, that $K_{ii}=0$ for
all $i\in\\{0,\dots,n\\}$, a convention to be adopted throughout this paper.

Supposing the natural symmetry hypothesis $K_{ij}=K_{ji}$, the pointwise sum
$\sum_{i=0}^{n}u_{i}$ of all densities is formally preserved under the
evolution. Assuming, for simplicity, that $\sum_{i=0}^{n}{u_{i}}_{|t=0}\equiv
1$, we are thus interested in solutions $u=(u_{0},\dots,u_{n})$ of the above
PDE system taking values in the planar hypersurface

$\displaystyle\mathcal{S}:=\left\\{u^{\prime}:=(u_{0}^{\prime},\ldots,u_{n}^{\prime})\in(0,\infty)^{n+1}:\;\sum_{i=0}^{n}u_{i}^{\prime}=1\right\\}.$
(1.2)

Our study is based on the well-known observation (cf. [6, 16, 24]) that system
(1.1) can be written as a gradient flow of the convex integral functional

$\displaystyle\mathcal{H}(u)=\int_{\Omega}h(u)\,\differential x,\qquad
h(u)=\sum_{i=0}^{n}\lambda(u_{i})$ (1.3)

with $\lambda(s):=s\log s-s+1$, that takes the form

$\displaystyle\begin{aligned}
\partial_{t}u_{i}-\nabla\cdot\left(\sum_{l=0}^{n}\mathbb{M}_{il}(u)\nabla
D_{l}h(u)\right)=0\quad&\text{ in }(0,T^{*})\times\Omega,\quad
i\in\\{0,\ldots,n\\}\end{aligned}$ (1.4)

for a suitable symmetric and positive semi-definite mobility
$\mathbb{M}:\mathcal{\overline{S}}\to\mathbb{R}^{(n+1)\times(n+1)}$.

To see this, we first rewrite equation (1.1a) more concisely as

$\partial_{t}u_{i}-\nabla\cdot(\sum_{j}A_{ij}(u)\nabla u_{j})=0$

with the matrix $A=(A_{ij})_{i,j=0,\dots,n}$ given by

$\displaystyle A_{ij}(u)=\delta_{ij}\sum_{k=0}^{n}K_{ik}u_{k}-K_{ij}u_{i}.$
(1.5)

Here, $\delta_{ij}$ denotes the Kronecker symbol. Then, defining the mobility
matrix $\mathbb{M}(u):=A(u)(D^{2}h(u))^{-1}$, system (1.1a) turns into (1.4).
For later reference, we note that the components of $\mathbb{M}(u)$ take the
form

$\displaystyle\mathbb{M}_{il}(u)=u_{i}(\delta_{il}\sum_{k=0}^{n}K_{ik}u_{k}-K_{il}u_{l}),$
(1.6)

from which we see in particular that $\mathbb{M}(u)$ is symmetric for all
$u\in\mathcal{\overline{S}}$ if and only if $(K_{ij})$ is. Furthermore, as we
will see in Section 2.2, for $(K_{ij})$ symmetric the positive semi-
definiteness of $\mathbb{M}(u)$ for all $u\in\mathcal{\overline{S}}$ is
equivalent to the non-negativity of $K_{ij}$ for all $i\not=j$. This motivates
imposing the following mild hypotheses on the coefficients $K_{ij}$.

##### Hypotheses

  1. (H1)

Symmetry: $K_{ij}=K_{ji}$ for all $i,j\in\\{0,\dots,n\\}$, $i\not=j$.

  2. (H2)

Non-negativity: $K_{ij}\geq 0$ for all $i,j\in\\{0,\dots,n\\}$, $i\not=j$.

Note that from a microscopic point of view the $K_{ij}$ are proportional to
the rate at which particles of type $i$ and $j$ change sites and hence are
non-negative and symmetric. In particular, hypotheses (H1) & (H2) are
fulfilled in all relevant applications and will therefore be assumed
throughout this manuscript without necessarily being referred to. The symmetry
of the $K_{ij}$ implies that, upon adding up the equations (1.1a),

$\frac{\differential}{\differential
t}\sum_{i=0}^{n}u_{i}=\sum_{i=0}^{n}\partial_{t}u_{i}=0,$

and the non-negativity of $K_{ij}$ can be used to conclude preservation of
non-negativity of each $u_{i}$, hence they are directly related to obtaining a
solution in $\mathcal{\overline{S}}$. More specifically, under hypotheses
(H1), (H2) and an additional non-degeneracy condition on $\\{K_{ij}\\}$ (such
as (H2enumi) or (H3) below), a strong maximum-type principle (cf. e.g. [10,
Chapter III], [20]) ensures that smooth flows starting in the surface
$\mathcal{S}$ will not leave it as long as they maintain their regularity.

Our analysis further requires some non-degeneracy condition on the
coefficients $K_{ij}$. While existence and equilibration will be established
under a rather mild condition (to be specified in Section 4, see (H3)), our
stability analysis in Section 3 requires in addition to (H1) and (H2) the
full-interaction hypothesis

  1. (H2enumi)

Positivity: $K_{ij}>0$ for all $i,j\in\\{0,\dots,n\\}$, $i\not=j$.

Note that from a microscopic point of view, the strengthened condition
(H2enumi) means that all combinations of different particles can change sites
and hence there is no situation where particles can get stuck. In this
context, let us recall that the diagonal values $K_{ii}$ are irrelevant for
the dynamics and have been set to zero to simplify notations.

##### Approach

Exploiting the formal gradient-flow structure (1.4) is a key feature of our
analysis. At a formal level, (1.4) implies the following basic entropy
dissipation law

$\displaystyle\frac{\differential}{\differential
t}\mathcal{H}(u)+\int_{\Omega}\frac{1}{2}\sum_{i,j=0}^{n}K_{ij}u_{i}u_{j}\left|\frac{\nabla
u_{i}}{u_{i}}-\frac{\nabla u_{j}}{u_{j}}\right|^{2}\,\differential x=0$ (1.7)

along solutions $u=u(t,x)\in\mathcal{\overline{S}}$, which is the main
building block in our global existence and long-time asymptotic analysis. The
latter will be carried out under a hypothesis far weaker than (H2enumi) (cf.
Section 4, (H3)). To explain the new difficulties arising when relaxing
(H2enumi), let us note that hypothesis (H2enumi) ensures that the dissipated
quantity $\int_{\Omega}D$, where

$D:=\frac{1}{2}\sum_{i,j=0}^{n}K_{ij}u_{i}u_{j}\left|\frac{\nabla
u_{i}}{u_{i}}-\frac{\nabla u_{j}}{u_{j}}\right|^{2},$

provides $L^{2}$ control of the Sobolev gradient of $\sqrt{u_{i}}$, since
then111In fact, by inequality (2.28) in Lemma 2.3 one has $D\geq
4\min_{\\{\alpha,\ell:\alpha\not=\ell\\}}K_{\alpha\ell}\cdot\sum_{i=0}^{n}|\nabla\sqrt{u_{i}}|^{2}$.
$D\gtrsim\sum_{i=0}^{n}|\nabla\sqrt{u_{i}}|^{2}$. Without hypothesis (H2enumi)
this bound is no longer valid and the dissipation law (1.7) provides only poor
gradient control of solutions. Hence, we cannot rely on standard Sobolev space
theory and compactness arguments to pass to the limit in the approximation
scheme. Indeed, taking for instance $K_{i0}>0$ and $K_{ij}=0$ for all
$i,j\in\\{1,\dots,n\\}$, one even has

$\displaystyle\sum_{i,j=0}^{n}K_{ij}u_{i}u_{j}\left|\frac{\nabla
u_{i}}{u_{i}}-\frac{\nabla u_{j}}{u_{j}}\right|^{2}\leq
K\left(\frac{1}{u_{0}}|\nabla
u_{0}|^{2}+u_{0}\sum_{j=0}^{n}\frac{1}{u_{j}}|\nabla u_{j}|^{2}\right),$ (1.8)

where $K:=\max_{1\leq i\leq n}K_{0i}$, so that any gradient control of
$u_{j}$, $j\geq 1$, obtained from (1.7) is lost when $u_{0}$ vanishes. (This
bound follows from an estimate analogous to that in the proof of Lemma 2.3.)

For proving our weak-strong uniqueness and stability result (with the relative
entropy as a generalised distance), the dissipation property (1.7) is not
sufficient. Here, we follow [12] and rely on the structure (1.4) for
controlling the relative entropy of a weak solution with respect to a strong
solution. Developing further the approach in [12], we will see that weak-
strong stability is mostly a consequence of the general structural properties
of (1.4): apart from the smoothness properties of $h$ and $\mathbb{M}$, the
strict convexity222In this article, strict convexity of the function $h\in
C^{2}(\mathcal{S})$ is to be understood in the sense of the Hessian
$D^{2}h(u)$ being positive definite at every point $u\in\mathcal{S}$. of the
entropy, and the parabolicity of the gradient-flow evolution system (1.4), we
only rely on the a priori estimate (1.7) providing sufficient gradient
control. Precisely because of this last prerequisite our stability analysis is
confined to models $\\{K_{ij}\\}$ obeying the full-interaction hypothesis
(H2enumi). However, it significantly improves the perturbative results of [3],
which are based on the same assumption, but also need all values
$\\{K_{ij}\\}_{i\not=j}$ to be close to some positive constant.

A primary purpose of the present manuscript is to expose key aspects of the
entropy method for (degenerate) cross-diffusion systems such as (1.1) rather
than aim for optimal estimates, and we would like to mention that, in the case
of full interactions, a version of the weak-strong stability estimate (1.12)
may be obtained under relaxed hypotheses on the strong solution and the
initial data. In particular, using a suitably adapted relative entropy as
recently introduced in [5], it might be possible to admit general non-negative
initial data that are not necessarily strictly positive. Let us further note
that our existence theory and the weak-strong stability estimate are equally
valid in the presence of external potentials, that is, for models (1.4) where
the Boltzmann-type entropy $\sum_{i}\lambda(u_{i})$ is enhanced by a linear
potential part $\sum_{i}u_{i}\Upsilon_{i}$ for given smooth potentials
$\Upsilon_{i}:\overline{\Omega}\to\mathbb{R}$. This leads to evolution systems
(1.1a) of the form

$\displaystyle\partial_{t}u_{i}-\nabla\cdot\big{[}\sum_{j=0}^{n}K_{ij}(u_{j}\nabla
u_{i}-u_{i}\nabla
u_{j}+u_{i}u_{j}\nabla(\Upsilon_{i}-\Upsilon_{j}))\big{]}=0.$ (1.9)

An extension to mean-field interaction (cf. e.g. [4]) also seems possible in
certain regimes. Let us finally mention a parallel development [13] covering
Maxwell–Stefan diffusions in a non-degenerate regime (corresponding to the
case of full interactions).

###  1.1. Main results

Our main results can be divided into two parts, contained in Sections 3 and 4,
respectively.

  1. (i)

The first part concerns fully interacting systems (H1), (H2enumi): Here,
gradient bounds are available globally and existence of weak solutions is well
established [14]. In this setting, we are able to perform a fairly general
weak-strong stability and asymptotic analysis.

  2. (ii)

The second part considers more general partially interacting systems. These
are quite relevant for applications, but due to the degeneracy of the
parabolic structure even the issue of existence in a weak sense has not been
fully understood so far. We establish the existence of weak solutions (in a
slightly generalised sense) as well as their convergence to the steady state
in the long-time limit.

####  1.1.1. Results on full interactions in part (i)

Let $T^{*}\in(0,\infty]$.

######  Definition 1.1 (Weak solution).

We call $u=(u_{0},\dots,u_{n})$ a weak solution of (1.1) with if $u_{i}\in
L^{2}(0,T;H^{1}(\Omega))$, $\partial_{t}u_{i}\in
L^{2}(0,T;(H^{1}(\Omega))^{*})$ for all $T<T^{*}$ and $i\in\\{0,\dots,n\\}$,
if $u(t,x)\in\mathcal{\overline{S}}$ for a.e. $(t,x)\in(0,T^{*})\times\Omega$,
and if $u$ satisfies (1.1) in the weak sense, i.e. if for all $T<T^{*}$, all
$\psi\in L^{2}(0,T;H^{1}(\Omega))$ and all $i\in\\{0,\dots,n\\}$

$\displaystyle\int_{0}^{T}\langle\partial_{t}u_{i},\psi\rangle\,\differential
t+\int_{0}^{T}\\!\int_{\Omega}\sum_{j=0}^{n}A_{ij}(u)\nabla
u_{j}\cdot\nabla\psi\,\differential x\differential t=0.$ (1.10)

Observe that under the hypotheses of Definition 1.1 one has $0\leq u_{i}\leq
1$ for all $i$.

For systems (1.1) obeying hypothesis (H2enumi), we perform a weak-strong
stability analysis choosing the relative entropy

$\displaystyle\mathcal{H}_{\mathrm{rel}}(u,\tilde{u})$
$\displaystyle:=\int_{\Omega}h_{\mathrm{rel}}(u,\tilde{u})\,\differential
x:=\int_{\Omega}\big{[}h(u)-\sum_{i=0}^{n}D_{i}h(\tilde{u})(u_{i}-\tilde{u}_{i})-h(\tilde{u})\big{]}\,\differential
x\quad$ (1.11)

as a generalised distance between a weak solution $u$ and some strictly
positive reference solution $\tilde{u}$ that enjoys certain extra regularity
properties. Our most interesting result under the full-interaction hypothesis
consists in a weak-strong stability estimate, where ‘strong solutions’ are
defined in the following way.

######  Definition 1.2 (Strong solution).

We call $\tilde{u}=(\tilde{u}_{0},\dots,\tilde{u}_{n})$ a strong solution of
(1.1) if $\tilde{u}\in C^{0,1}([0,T^{*})\times\overline{\Omega})^{1+n}$, if
$\tilde{u}(t,x)\in\mathcal{S}$ for all $(t,x)\in[0,T]\times\overline{\Omega}$
and all $T<T^{*}$, and if $\tilde{u}$ satisfies (1.1) in the weak sense (as
specified in Definition 1.1).

Observe that our notion of a strong solution not only requires Lipschitz
regularity, but also strict positivity of all components $\tilde{u}_{i}$ on
$[0,T^{*})\times\overline{\Omega}$ and as a consequence
$\min_{[0,T]\times\overline{\Omega}}\tilde{u}_{i}>0$ for all $T<T^{*}$ and all
$i\in\\{0,\dots,n\\}$. Let us further remark that under hypotheses (H1),
(H2enumi) (and more generally under (H1), (H2) and the condition (H3) to be
introduced in Section 4) for given Lipschitz continuous initial data
$u^{\mathrm{in}}$ with $u^{\mathrm{in}}(\overline{\Omega})\subset\mathcal{S}$,
local existence of a unique strong solution follows from Amann [1] and the
parabolic structure of the system (see Section 2.1 for details).

Our weak-strong uniqueness and stability result states as follows.

######  Theorem 1.3 (Weak-strong stability).

Assume hypotheses (H1), (H2enumi) and suppose that $\tilde{u}\in
C^{0,1}([0,T^{*})\times\overline{\Omega})^{1+n}$ is a strong solution in the
sense of Definition 1.2. For all $T\in(0,T^{*})$ there exists a constant
$C_{T,\tilde{u}}<\infty$ only depending on
$\|\tilde{u}\|_{C^{0,1}([0,T]\times\overline{\Omega})}$ and
$(\inf_{[0,T]\times\overline{\Omega}}\tilde{u}_{i})_{i=0}^{n}$ such that the
following holds true:

any weak solution $u$ in the sense of Definition 1.1 that has the additional
regularity $\sqrt{u_{i}}\in L^{2}_{\mathrm{loc}}([0,T^{*});H^{1}(\Omega))$ for
all $i\in\\{0,\ldots,n\\}$ obeys the stability estimate

$\displaystyle\mathcal{H}_{\mathrm{rel}}(u(t),\tilde{u}(t))\leq\mathcal{H}_{\mathrm{rel}}(u(0),\tilde{u}(0))\exp(C_{T,\tilde{u}}t)\quad\text{
for al }t\in[0,T],$ (1.12)

where $\mathcal{H}_{\mathrm{rel}}$ denotes the relative entropy functional
(1.11).

If $u_{|t=0}=\tilde{u}_{|t=0}$, the above theorem implies that $u=\tilde{u}$
a.e. in $(0,T^{*})\times\Omega$.

See Section 3.2 for the proof of Theorem 1.3. The regularity hypothesis
$\sqrt{u_{i}}\in L^{2}_{\mathrm{loc}}([0,T^{*});H^{1}(\Omega))$ in Theorem 1.3
is natural under hypothesis (H2enumi), and the existence of weak solutions
obeying this hypothesis is well established, see Theorem 2.5 below.

For completeness, we further note the following asymptotic stability result,
which follows from classical arguments in the full-interaction regime.

######  Theorem 1.4 (Exponential stability of steady states).

Assume hypotheses (H1), (H2enumi) and let
$\overline{u}\in\mathcal{\overline{S}}$. There exist finite positive constants
$\epsilon=\epsilon(\Omega)$ and $C=C(\Omega)$ such that any weak solution
$u=(u_{0},\dots,u_{n})$ of (1.1) in the sense of Definition 1.1 with the extra
regularity $\sqrt{u_{i}}\in L^{2}_{\mathrm{loc}}([0,T^{*});H^{1}(\Omega))$ and
average values $\frac{1}{|\Omega|}\int_{\Omega}u_{i}(0)\,\differential
x=\overline{u}_{i}$ for all $i\in\\{0,\dots,n\\}$ satisfies the decay estimate

$\displaystyle\|u(t)-\overline{u}\|_{L^{1}(\Omega)}^{2}\leq
C\mathcal{H}_{\mathrm{rel}}(u(0),\overline{u})\mathrm{e}^{-\epsilon\kappa
t}\qquad\text{ for all }t>0,$ (1.13)

where $\kappa:=\min_{i\not=j}K_{ij}>0$.

The proof of Theorem 1.4 is given in Section 3.3.

####  1.1.2. Results in part (ii) on partial interactions

Since the precise statement of our results concerning partial interactions
requires several additional concepts and notations that will be introduced in
detail in Section 4, we here prefer to only summarise the main results on part
(ii):

  * •

Large-data global existence of weak solutions (see Theorem 4.2),

  * •

Convergence to equilibrium of the constructed solutions (see Theorem 4.4).

While for specific subclasses of models existence has previously been obtained
in the literature (see [6, 14, 24]), our result concerning the long-time
asymptotics appears to be the first rigorous one for partially interacting
systems.

We should mention that for certain degenerate classes of partially interacting
models the regularity obtained from entropy estimates does not suffice to
uniquely determine all flux terms $K_{ij}(u_{j}\nabla u_{i}-u_{i}\nabla
u_{j})$ appearing in (1.1a). In fact, in contrast to the models covered by
existing literature [6, 14, 24], the present existence analysis faces
situations where, due to poor gradient control, expressions of the form
$K_{ij}u_{j}\nabla u_{i}$ may not even be defined in the distributional sense
(despite the fact that $0\leq u_{l}\leq 1$ for all $l$). Hence, our existence
result, Theorem 4.2, is formulated for a somewhat generalised notion of weak
solutions that, however, reduces to the standard concept of weak or
distributional solutions when restricting to the models considered in [6, 14,
24]. Besides, the generalised notion of weak solutions in Theorem 4.2 respects
mass conservation and allows to uniquely determine the long-time asymptotic
behaviour of solutions (cf. Theorem 4.4).

###  1.2. Previous literature

Below, we provide an overview of the analytical results available for system
(1.1). We should point out that [6] was the first work to analyse such a
system (for $n=2$).
Literature on global existence:

  * •

Full interactions:

    * –

the case $K_{ij}>0$ for all $i\not=j$ has been fully covered by [14, Theorem
2], see [2] for details.

    * –

[2]: global existence of weak solutions for non-zero flux boundary conditions
and time-dependent domains in one space dimension $d=1$. Equilibration for
suitable data.

  * •

Partial interactions:

    * –

[6]: global existence of weak solutions for $n=2$ and $K_{i0}>0,K_{ij}=0$ for
$j,i\geq 1$. (Only the case $d\in\\{1,2,3\\}$ was considered, but the proof
does not rely on this restriction.)

    * –

[14, Theorem 3]: global existence of weak solutions again for $n=2$ and
$K_{i0}>0,K_{ij}=0$ for $j,i\geq 1$, but more general transition rates (not
considered in the present manuscript).

    * –

[24]: global existence of weak solutions for general $n\geq 2$ and
$K_{i0}>0,K_{ij}=0$ for $j,i\geq 1$.

Literature on convergence to equilibrium (partial interactions):

  * •

[6, Section 5]: non-rigorous sketch proof that assumes monotonicity of the
entropy.333The argument can be made rigorous by following the ideas in Section
4, and in particular the proof of Theorem 4.4.

  * •

[24, Theorem 4]: attempt to rely on decoupling to deduce an exponential
convergence rate. (However, the crucial Gronwall type argument therein appears
to be false and fixing it would require a new idea.)

Literature on uniqueness (partial interactions):

  * •

[24, Theorem 5]: uniqueness of weak solutions in a situation where the problem
essentially decouples ($K_{0i}=a_{0}>0$ for all $i\in\\{1,\dots,n\\}$ and
$K_{ij}=0$ for $i,j\geq 1$).

Let us finally note that in perturbative or small data settings classical
techniques apply. Such regimes are mathematically better behaved and easier to
handle:

  * •

[6, Section 3]: wellposedness of strong solutions near equilibrium.

  * •

[3]: existence and uniqueness of strong solutions when the interaction
coefficients $K_{ij},i\not=j$ are positive and close to each other. In this
case, the system is close to decoupled linear diffusion, and the authors apply
classical $L^{2}$ methods to study the wellposedness in this perturbative
setting.

###  1.3. New aspects of the present results

As discussed in the previous paragraph, so far, uniqueness results have only
been available for small data settings or under rather restrictive assumptions
on the model. In the present work we remove such smallness restrictions on the
model and establish a weak-strong stability estimate for system (1.1) with
general coefficients satisfying (H1), (H2enumi), see Theorem 1.3. Our method
relies on an adaptation of the technique in [12] devised for proving weak-
strong uniqueness in energy-reaction-(cross-)diffusion systems. This approach
in turn was motivated by an adjusted relative entropy method due to Fischer
[11] for reaction-diffusion systems. A crucial ingredient allowing us to
handle the new difficulties induced by the size-exclusion effect is a local
coercivity-type estimate (cf. inequality (3.21)), which we deduce from the
strict convexity of the entropy density.

We further establish existence of weak solutions and equilibration in the
long-time limit under a fairly mild hypothesis on $K_{ij}$ arguably covering
all models relevant for applications. Here, our idea is to take advantage of
the existence result for full interactions [14], which provides us with a
family of approximate solutions. Upon an analysis of the compactness and
convergence properties of these approximate solutions, we then take a
‘vanishing-interaction limit’, where the artificially introduced interaction
parameters are being sent to zero, to eventually obtain weak solutions to the
partially interacting system. Owing to poor compactness properties linked to
the degeneracy of such systems, we are not able to establish a strong entropy
dissipation inequality nor monotonicity of the entropy in time. However, we
will show that a version of the strong entropy dissipation inequality holds
true up to an error term (cf. inequality (4.4d)) that is initially only known
to tend to zero along a some sequence of times $t_{k}\to\infty$, along which
stronger convergence properties hold. Inequality (4.4d), however, allows us to
upgrade this result and prove the strong convergence to equilibrium along any
sequence $t\to\infty$.

###  1.4. Outline

In Section 2, we use a change of variables to remove the degeneracy of
equation (1.1) that is linked to the volume constraint. We further establish
basic algebraic identities and coercivity estimates involving the mobility
matrix that will be needed in subsequent parts. The stability results for
fully interacting systems are proved in Section 3. In Section 4, the full-
interaction hypothesis is dropped, and general, partially interacting systems
relevant for applications are introduced. For these problems, we establish
existence of weak solutions and convergence to equilibrium.

This paper has two appendices. Appendix A contains a technical lemma used in
the proof of the weak-strong stability estimate. For details concerning our
notations we refer to Appendix B.

##  2\. Preliminaries

###  2.1. The gradient system for $n$ species

Consistent with the size constraint, the mobility $\mathbb{M}$ underlying
system (1.1) degenerates in the direction
$\textbf{n}_{\mathcal{S}}:=(1,\dots,1)^{T}$ normal to the hyperplane
containing $\mathcal{S}$. This degeneracy can be removed by means of an affine
coordinate transformation mapping the $n$-dimensional domain
$\mathcal{D}\subset\mathbb{R}^{n}$, given by

$\displaystyle\mathcal{D}:=\left\\{U^{\prime}=(U_{1}^{\prime},\ldots,U_{n}^{\prime})\in(0,\infty)^{n}:\;\sum_{i=1}^{n}U_{i}^{\prime}<1\right\\},$
(2.1)

diffeomorphically onto $\mathcal{S}$. Such a transformation can be realised by
the map

$\displaystyle\Psi:\mathcal{D}\to\mathcal{S},\qquad
U=(U_{1},\dots,U_{n})\mapsto(1-\sum_{i=1}^{n}U_{i},U_{1},\dots,U_{n}).$ (2.2)

By means of this change of variables, the gradient structure with respect to
$\mathcal{S}$ gives rise to a gradient flow for the vector of n densities
$U(t,x)\in\mathcal{D}$, as will be detailed below. Of course, we could have
immediately written down a PDE system for $U$ equivalent to (1.1) by replacing
$u_{0}$ by $1-\sum_{i=1}^{n}U_{i}$ and rearranging terms. However, we hope
that the following few expository paragraphs may help better understand some
of the underlying structures.

##### Entropy density

For $U\in\mathcal{D}$ we define the pulled-back entropy density

$\displaystyle\hat{h}(U):=h(\Psi(U))=\sum_{i=1}^{n}\lambda(U_{i})+\lambda\big{(}1-\sum_{i=1}^{n}U_{i}\big{)},$
(2.3)

where as before $\lambda(s)=s\log s-s+1$. Observe that $\hat{h}\in
C(\mathcal{\overline{D}})\cap C^{\infty}(\mathcal{D})$ is strictly convex.
Abbreviating $v:=1-\sum_{i=1}^{n}U_{i}$, we have

$\displaystyle
D_{i}\hat{h}(U)=\log(U_{i})-\log(v)=\log\big(\tfrac{U_{i}}{v}\big{missing})$
(2.4)

and

$\displaystyle
D^{2}\hat{h}(U)=\operatorname{diag}(\tfrac{1}{U_{1}},\dots,\tfrac{1}{U_{n}})+\tfrac{1}{v}E,$
(2.5)

where $E_{ij}=1$ for all $i,j\in\\{1,\dots,n\\}$.

##### Mobility

Since the Onsager operator $\operatorname{div}(\mathbb{M}(u)\nabla\,\cdot\,)$
is acting on the dual variables, computing the appropriate mobility
$\mathbb{\hat{M}}$ on $\mathcal{D}$ requires using the inverse
$\Phi:=\Psi^{-1}$ of $\Psi$,

$\displaystyle\Phi:\mathcal{S}\to\mathcal{D},\quad
u=(u_{0},\dots,u_{n})\mapsto U=(u_{1},\dots,u_{n}).$ (2.6)

We may now define the $(n\times n)$ mobility matrix

$\displaystyle\mathbb{\hat{M}}(U):=D\Phi_{|\Psi(U)}\mathbb{M}(\Psi(U))D\Phi^{T}_{|\Psi(U)},\quad
U\in\mathcal{D},$ (2.7)

which is again symmetric. (Since $\Phi$ is affine, the argument in
$D\Phi_{|\Psi(U)}$ may be dropped.)

##### Gradient-flow equation

The above ingredients determine the gradient-flow equation (1.4) in the new
variables as

$\left\\{\begin{array}[]{rll}\partial_{t}U-\nabla\cdot(\mathbb{\hat{M}}(U)\nabla
D\hat{h}(U))&=0&\text{ in }\quad(0,T^{*})\times\Omega,\\\
\mathbb{\hat{M}}(U)\nabla D\hat{h}(U)\cdot\nu&=0&\text{ on
}\quad(0,T^{*})\times\partial\Omega.\\\ \end{array}\right.$ (2.8)

Introducing the ‘diffusion matrix’

$\displaystyle\hat{A}(U):=\mathbb{\hat{M}}(U)D^{2}\hat{h}(U),$ (2.9)

equation (2.8) may further be written as the PDE system

$\left\\{\begin{array}[]{rll}\partial_{t}U-\nabla\cdot(\hat{A}(U)\nabla
U)&=0&\text{ in }\quad(0,T^{*})\times\Omega,\\\ (\hat{A}(U)\nabla
U)\cdot\nu&=0&\text{ on }\quad(0,T^{*})\times\partial\Omega.\\\
\end{array}\right.$ (2.10)

A direct computation yields the following explicit form for the matrix
$\hat{A}(U)$

$\hat{A}_{ij}(U)=\delta_{ij}\left[\sum_{l=1}^{n}K_{il}u_{l}+K_{i0}(1-\sum_{l=1}^{n}u_{l})\right]-(K_{ij}-K_{i0})u_{i},\quad
i,j\in\\{1,\ldots,n\\}.$ (2.11)

Observe that system (2.10) can also directly be obtained from (1.1) by
inserting $u_{0}=1-\sum_{i=1}^{n}U_{i}$ and rearranging terms.

##### Entropy dissipation

We emphasize that the information contained in (2.8) is the same as in the
original system. In particular, we assert that

$\displaystyle
D^{2}\hat{h}\mathbb{\hat{M}}D^{2}\hat{h}=D\Psi^{T}(D^{2}h\circ\Psi)(\mathbb{M}\circ\Psi)(D^{2}h\circ\Psi)D\Psi.$
(2.12)

This follows from the identities
$D^{2}\hat{h}=D\Psi^{T}(D^{2}h\circ\Psi)D\Psi$ and

$\displaystyle(D\Psi D\Phi)\mathbb{M}(u)(D\Psi D\Phi)^{T}=\mathbb{M}(u),$
(2.13)

which can be verified explicitly by an elementary calculation using the
structure of $\mathbb{M}$. Hence, as a consequence of (2.12), the dissipated
quantities of the flows in transformed and original coordinates are indeed the
same:

$\displaystyle\begin{multlined}\hat{\mathscr{P}}(U):=\sum_{i,j=1}^{n}\nabla
U_{i}\cdot(D^{2}\hat{h}(U)\mathbb{\hat{M}}(U)D^{2}\hat{h}(U))_{ij}\nabla
U_{j}\\\ \hskip 86.72267pt=\sum_{k,l=0}^{n}\nabla
u_{k}\cdot(D^{2}h(u)\mathbb{M}(u)D^{2}h(u))_{kl}\nabla
u_{l}=:\mathscr{P}(u).\end{multlined}\hat{\mathscr{P}}(U):=\sum_{i,j=1}^{n}\nabla
U_{i}\cdot(D^{2}\hat{h}(U)\mathbb{\hat{M}}(U)D^{2}\hat{h}(U))_{ij}\nabla
U_{j}\\\ \hskip 86.72267pt=\sum_{k,l=0}^{n}\nabla
u_{k}\cdot(D^{2}h(u)\mathbb{M}(u)D^{2}h(u))_{kl}\nabla u_{l}=:\mathscr{P}(u).$
(2.16)

##### Parabolic structure

Assuming a mild non-degeneracy condition on the coefficients $K_{ij}$, the
problem formulated in terms of $n$ species enjoys full parabolicity in a
classical sense, as will be explained in the following. For simplicity, we
only consider the case of a fully interacting system obeying (H2enumi), which
guarantees that $\kappa:=\min_{i,j\in\\{0,\dots,n\\},i\not=j}K_{ij}$ is
strictly positive. It has been shown in the proof of [2, Lemma 2.3] that

$\displaystyle\zeta^{T}D^{2}\hat{h}(U)\hat{A}(U)\zeta\geq\kappa\sum_{i=1}^{n}U_{i}^{-1}\zeta_{i}^{2}\qquad\text{
for }U\in\mathcal{D}.$ (2.17)

See also Corollary 2.4 below, which actually implies this bound (although our
reasoning is different and perhaps easier to grasp). In terms of
$\mathbb{\hat{M}}$, inequality (2.17) takes the form

$\displaystyle\zeta^{T}D^{2}\hat{h}(U)\mathbb{\hat{M}}(U)D^{2}\hat{h}(U)\zeta\geq\kappa\sum_{i=1}^{n}U_{i}^{-1}\zeta_{i}^{2}\qquad\text{
for }U\in\mathcal{D}.$ (2.18)

Since for each $U\in\mathcal{D}$ the matrix $D^{2}\hat{h}(U)$ is invertible,
the last inequality implies that $\mathbb{\hat{M}}(U)>0$ for all
$U\in\mathcal{D}$, where we recall that $\mathbb{\hat{M}}(U)$ is symmetric.
Hence, by continuity, for every compact subset
$\mathcal{K}\subset\subset\mathcal{D}$ there exists a positive constant
$c_{*}(\mathcal{K})>0$ such that

$\displaystyle\mathbb{\hat{M}}(U)\geq
c_{*}(\mathcal{K})\operatorname{diag}(1,\dots,1)\qquad\text{ for all
}U\in\mathcal{K}.$ (2.19)

We further note that, since for every $U\in\mathcal{D}$ the diffusion matrix
$\hat{A}(U)$ is the product of the positive definite symmetric matrices
$\mathbb{\hat{M}}(U)$ and $D^{2}\hat{h}(U)$, it has its spectrum
$\sigma(\hat{A}(U))$ contained in $\mathbb{R}_{>0}$. Therefore, thanks to
Amann’s work [1], which requires the real parts of the eigenvectors of
$\hat{A}(U)$ to be strictly positive, the Cauchy problem associated with
system (2.10) has a unique local-in-time smooth solution for Lipschitz regular
initial data $U^{\mathrm{in}}$ satisfying
$U^{\mathrm{in}}(\overline{\Omega})\subset\mathcal{D}$.

######  Remark.

In fact, under hypothesis (H2enumi), the inclusion
$\sigma(\hat{A}(U))\subset\mathbb{R}_{>0}$ sufficient for applying [1] even
holds for all $U\in\mathcal{\overline{D}}$ (and thus by continuity also in an
open neighbourhood of $\mathcal{\overline{D}}$). This is a consequence of [15,
Lemma 2.3]. Hence, local wellposedness is even valid under the more general
condition $U^{\mathrm{in}}(\overline{\Omega})\subset\mathcal{\overline{D}}$.
Preservation of positivity and size constraint are, however, not fully trivial
if $U^{\mathrm{in}}(\overline{\Omega})\cap\partial\mathcal{D}\not=\emptyset$
and will have to be verified a posteriori.

##### Equivalence of the two PDE systems

The $n$-species diffusion system (2.10) precisely agrees with the last $n$
components of the original system (1.1). This may be verified explicitly by
using formula (2.11) and the relation $\sum_{i=0}^{n}u_{i}=1$. Since the
equation for the zeroth component $u_{0}=1-\sum_{i=1}^{n}U_{i}$ in (1.1) is
obtained by summing up the equations for the components $i=1,\dots,n$, system
(2.10) is equivalent to the original problem. For later reference, we
formulate this observation for the concept of weak solutions in the following
lemma.

######  Lemma 2.1.

System (1.10) for $i\in\\{0,\dots,n\\}$ is equivalent to the weak form of
(2.10), i.e. to

$\displaystyle\int_{0}^{T}\langle\partial_{t}U_{i},\psi\rangle\,\differential
t+\int_{0}^{T}\\!\int_{\Omega}\sum_{j=1}^{n}\hat{A}_{ij}(U)\nabla
U_{j}\cdot\nabla\psi\,\differential x\differential t=0$ (2.20)

for $i\in\\{1,\dots,n\\}$. More precisely, $u=(u_{0},\dots,u_{n})$ being a
weak solution of (1.1) in the sense of Definition 1.1 is equivalent to
$U=(u_{1},\dots,u_{n})$ satisfying $u_{i}\in L^{2}(0,T;H^{1}(\Omega))$,
$\partial_{t}u_{i}\in L^{2}(0,T;(H^{1}(\Omega))^{*})$ for all $T<T^{*}$ and
$i\in\\{1,\dots,n\\}$, fulfilling $U(t,x)\in\mathcal{\overline{D}}$ for a.e.
$(t,x)\in(0,T^{*})\times\Omega$, and obeying (2.20) for all $T<T^{*}$ and all
$\psi\in L^{2}(0,T;H^{1}(\Omega))$.

We should note that the dimension reduction leading to full parabolicity comes
at a price: the entropy function in the new coordinates induces a strong
coupling between the $n$ unknowns, potentially introducing new difficulties in
the uniqueness and stability analysis.

###  2.2. Algebraic estimates

Here, we establish coercivity bounds for the quadratic form induced by the
symmetric and positive semi-definite matrix

$\displaystyle\mathbb{P}(u):=D^{2}h(u)\mathbb{M}(u)D^{2}h(u),\quad
u\in\mathcal{S},$ (2.21)

on the tangent space $T_{u}\mathcal{S}$ to $\mathcal{S}$. Since the vector
$\textbf{n}_{\mathcal{S}}:=(1,\dots,1)$ points in the direction perpendicular
to the hyperplane containing $\mathcal{S}$, we can identify

$\displaystyle
T_{u}\mathcal{S}=\\{\xi\in\mathbb{R}^{1+n}:\xi\cdot\textbf{n}_{\mathcal{S}}=0\\}.$
(2.22)

######  Lemma 2.2.

Let $K_{ij}$ satisfy Hypotheses (H1) & (H2). Further recall that the mobility
matrix $\mathbb{M}$ is given by (1.6), and let $\mathbb{P}$ be given by (2.21)
(with $h$ as in (1.3)). For all $u\in\mathcal{S}$ and all $\xi\in
T_{u}\mathcal{S}$

$\displaystyle\xi^{T}\mathbb{P}(u)\xi=\frac{1}{2}\sum_{i,j}K_{ij}u_{i}u_{j}\left|\frac{\xi_{i}}{u_{i}}-\frac{\xi_{j}}{u_{j}}\right|^{2}.$
(2.23)

While identity (2.23) holds for general $\xi\in\mathbb{R}^{1+n}$, only tangent
vectors $\xi$ are relevant in our analysis.

###### Proof of Lemma 2.2.

Observe that $D^{2}h(u)_{ij}=\frac{1}{u_{i}}\delta_{ij}$, and hence

$\displaystyle\mathbb{P}_{ij}(u)=\frac{1}{u_{j}}\delta_{ij}\sum_{k=0}^{n}K_{ik}u_{k}-K_{ij}.$
(2.24)

We then compute

$\displaystyle\sum_{i,j}\mathbb{P}_{ij}(u)\xi_{i}\xi_{j}$
$\displaystyle=\sum_{i,k}K_{ik}\frac{u_{k}}{u_{i}}\xi_{i}^{2}-\sum_{i,j}K_{ij}\xi_{i}\xi_{j}$
(2.25)
$\displaystyle=\frac{1}{2}\sum_{i,j}K_{ij}u_{i}u_{j}\big{(}\left|\tfrac{\xi_{i}}{u_{i}}\right|^{2}+\left|\tfrac{\xi_{j}}{u_{j}}\right|^{2}\big{)}-\frac{1}{2}\sum_{i,j}K_{ij}u_{i}u_{j}2\tfrac{\xi_{i}}{u_{i}}\tfrac{\xi_{j}}{u_{j}},$
(2.26)

where second step uses the symmetry of $K_{ij}$. The last line equals the
right-hand side of (2.23) and the assertion follows. ∎

######  Lemma 2.3 (Matrix coercivity estimates).

Assume the hypotheses and use the notations of Lemma 2.2. Then for all
$u\in\mathcal{S}$ and $\xi\in T_{u}\mathcal{S}$

$\displaystyle\xi^{T}\mathbb{P}(u)\xi\geq\sum_{\alpha=0}^{n}\kappa^{(\alpha)}\bigg{(}\frac{1}{u_{\alpha}}|\xi_{\alpha}|^{2}+u_{\alpha}\sum_{j=0}^{n}\frac{1}{u_{j}}|\xi_{j}|^{2}\bigg{)},$
(2.27)

where
$\kappa^{(\alpha)}:=\frac{1}{2}\min_{\ell:\ell\not=\alpha}K_{\alpha\ell}\geq
0$ for $\alpha\in\\{0,\dots,n\\}$.

In particular,

$\displaystyle\xi^{T}\mathbb{P}(u)\xi\geq\min_{\alpha,\ell:\alpha\not=\ell}K_{\alpha\ell}\cdot\sum_{j=0}^{n}\frac{1}{u_{j}}|\xi_{j}|^{2}.$
(2.28)

###### Proof.

Identity (2.23) allows us to estimate

$\displaystyle\xi^{T}\mathbb{P}(u)\xi$
$\displaystyle=\frac{1}{2}\sum_{\alpha=0}^{n}\sum_{j:j\not=\alpha}K_{\alpha
j}u_{\alpha}u_{j}\left|\frac{\xi_{\alpha}}{u_{\alpha}}-\frac{\xi_{j}}{u_{j}}\right|^{2}$
(2.29)
$\displaystyle\geq\sum_{\alpha=0}^{n}\kappa^{(\alpha)}u_{\alpha}\sum_{j:j\not=\alpha}u_{j}\bigg{(}\frac{1}{u_{\alpha}^{2}}\xi_{\alpha}^{2}+\frac{1}{u_{j}^{2}}\xi_{j}^{2}-2\frac{1}{u_{\alpha}u_{j}}\xi_{\alpha}\xi_{j}\bigg{)}$
(2.30)
$\displaystyle=\sum_{\alpha=0}^{n}\kappa^{(\alpha)}u_{\alpha}\bigg{(}\frac{1-u_{\alpha}}{u_{\alpha}^{2}}\xi_{\alpha}^{2}+\sum_{j:j\not=\alpha}\frac{1}{u_{j}}\xi_{j}^{2}+2\frac{1}{u_{\alpha}}\xi_{\alpha}^{2}\bigg{)}$
(2.31)
$\displaystyle=\sum_{\alpha=0}^{n}\kappa^{(\alpha)}\bigg{(}\frac{1}{u_{\alpha}}|\xi_{\alpha}|^{2}+u_{\alpha}\sum_{j=0}^{n}\frac{1}{u_{j}}|\xi_{j}|^{2}\bigg{)},$
(2.32)

where the penultimate equality uses the fact that
$\xi\cdot\textbf{n}_{\mathcal{S}}=0$, that is, $\sum_{j=0}^{n}\xi_{j}=0$. This
shows (2.27).

To infer inequality (2.28) it suffices to note that the right-hand side of
(2.27) is bounded below by

$\displaystyle\frac{1}{2}\min_{\alpha,\ell:\alpha\not=\ell}K_{\alpha\ell}\cdot\sum_{\beta=0}^{n}\bigg{(}\frac{1}{u_{\beta}}|\xi_{\beta}|^{2}+u_{\beta}\sum_{j=0}^{n}\frac{1}{u_{j}}|\xi_{j}|^{2}\bigg{)}=\frac{1}{2}\min_{\alpha,\ell:\alpha\not=\ell}K_{\alpha\ell}\cdot\bigg{(}2\sum_{j=0}^{n}\frac{1}{u_{j}}|\xi_{j}|^{2}\bigg{)}.$
(2.33)

∎

######  Corollary 2.4.

Let $K_{ij}$ satisfy Hypotheses (H1), (H2) and define

$\displaystyle\mathbb{\hat{P}}(U):=D^{2}\hat{h}(U)\mathbb{\hat{M}}(U)D^{2}\hat{h}(U),\qquad
U\in\mathcal{D}.$ (2.34)

Then, for all $U\in\mathcal{D}$ and $\zeta\in\mathbb{R}^{n}$,

$\displaystyle\zeta^{T}\mathbb{\hat{P}}(U)\zeta\geq\kappa\bigg{(}|\sum_{l=1}^{n}\zeta_{l}|^{2}+\sum_{i=1}^{n}\frac{1}{U_{i}}|\zeta_{i}|^{2}\bigg{)},$
(2.35)

where $\kappa:=\min_{j,k\in\\{0,\dots,n\\},j\not=k}K_{jk}$.

###### Proof.

This follows from Lemma 2.3 by inserting the tangent vector $\xi:=D\Psi\zeta$
and performing some elementary manipulations. Note that here we have also used
the identity (2.12). ∎

###  2.3. Existence of weak solutions under hypothesis (H2enumi)

Observe that under hypothesis (H2enumi), i.e. for
$\kappa:=\min_{i\not=j}K_{ij}>0$, Corollary 2.4 implies that
$\zeta^{T}\mathbb{\hat{P}}(U)\zeta\gtrsim\sum_{i=1}^{n}\frac{1}{U_{i}}|\zeta_{i}|^{2}$.
This estimate, combined with the convexity properties of $\hat{h}$, allows us
to apply Theorem 2 of [14] to infer global existence of weak solutions to
system (1.1) under the condition (H2enumi), where we also rely on Lemma 2.1 on
the equivalence between the original and the transformed system:

######  Theorem 2.5 (See [14] & L. 2.1).

Let $\min_{i\not=j}K_{ij}>0$. Let $u^{\mathrm{in}}\in L^{1}(\Omega)^{1+n}$
satisfy $u^{\mathrm{in}}(x)\in\mathcal{\overline{S}}$ for a.e. $x\in\Omega$.
Then there exists a function $u\in
L^{2}_{\mathrm{loc}}([0,\infty);H^{1}(\Omega))^{1+n}$ with $\partial_{t}u\in
L^{2}_{\mathrm{loc}}([0,\infty);H^{1}(\Omega)^{*})^{1+n}$ and
$u(t,x)\in\mathcal{\overline{S}}$ for a.e. $(t,x)\in(0,\infty)\times\Omega$
that satisfies system (1.1) in the weak sense and takes the initial datum
$u(0)=u^{\mathrm{in}}$. Furthermore, the regularity

$\displaystyle\sqrt{u_{i}}\in
L^{2}_{\mathrm{loc}}([0,\infty);H^{1}(\Omega))\qquad\textit{ for all
}i\in\\{0,\dots,n\\}$ (2.36)

holds true.

The main feature of the proof of Theorem 2.5 consists in a transformation to
the Legendre conjugate variables $w=D\hat{h}(U)$ upon which the system can be
regularised by a standard higher-order elliptic term while preserving the key
entropy estimate.

We should note that the statement of Theorem 2 by Jüngel [14] assumes that
$u^{\mathrm{in}}(x)\in\mathcal{S}$ for a.e. $x\in\Omega$. However, as observed
by the author in the paragraph following the statement of this theorem, an
approximation argument allows to relax this hypothesis to
$u^{\mathrm{in}}(x)\in\mathcal{\overline{S}}$ for a.e. $x\in\Omega$. It is
further necessary to point out that the regularity (2.36) does not appear in
the statement of [14, Theorem 2], but it directly follows from the proof of
this result, see [14, proof of Theorem 2; in particular, equation (23) and
Step 3].

Let us remark that in [2, Theorem 2.2] a result similar to Theorem 2.5 appears
(formulated for the $n$-species system), again deduced by verifying the
hypotheses of Theorem 2 of [14]. But since it is stated in a slightly weaker
form, we will not use it here.

##  3\. Stability results for full interactions

###  3.1. Entropy dissipation balance

Theorems 1.3 and 1.4 both rely on the following entropy dissipation identity.

######  Proposition 3.1.

Let $u=(u_{0},\dots,u_{n})$ be a weak solution according to Definition 1.1
enjoying the regularity $\sqrt{u_{i}}\in
L^{2}_{\mathrm{loc}}([0,T^{*});H^{1}(\Omega))$ for all $i$. Then for all
$T\in(0,T^{*})$

$\displaystyle\int_{\Omega}h(u(T))\,\differential
x=\int_{\Omega}h(u(0))\,\differential
x-\int_{0}^{T}\\!\int_{\Omega}\mathscr{P}(u)\,\differential x\differential t,$
(3.1)

where $\mathscr{P}(u)=\sum_{i,j}\nabla u_{i}\cdot\mathbb{P}_{ij}(u)\nabla
u_{j}=\sum_{i,j}\nabla u_{i}\cdot(D^{2}h(u)\mathbb{M}(u)D^{2}h(u))_{ij}\nabla
u_{j}$ (cf. (2.16)).

As a consequence, for all $0\leq s<T<T^{*}$

$\displaystyle\int_{\Omega}h(u(T))\,\differential
x=\int_{\Omega}h(u(s))\,\differential
x-\int_{s}^{T}\\!\int_{\Omega}\mathscr{P}(u)\,\differential x\differential t.$
(3.2)

###### Proof of Prop. 3.1.

We begin by observing that $u\in C([0,T];L^{2}(\Omega))^{1+n}$. This is a
consequence of the regularity hypotheses in Definition 1.1 combined with
standard embeddings for Bochner spaces. Formally, identity (3.1) follows upon
testing (1.10) with $D_{i}h(u)$, taking the sum over $i$ and applying a chain
rule. However, the Boltzmann-type entropy density $h$ is not smooth up to
$u_{i}=0$, and hence, we need to regularise first. The procedure is very
standard, but since it has not been carried out for this particular system,
some details are required.

For $\varepsilon>0$ let
$h_{\varepsilon}(u):=h(u_{0}+\varepsilon,\dots,u_{n}+\varepsilon)=\sum_{i=0}^{n}\lambda(u_{i}+\varepsilon).$
Then

$\displaystyle
D_{i}h_{\varepsilon}(u)=\log(u_{i}+\varepsilon)\quad\in\,L^{2}_{\mathrm{loc}}([0,T^{*});H^{1}(\Omega))$
(3.3)

and $D_{ik}h_{\varepsilon}(u)=\tfrac{1}{u_{i}+\varepsilon}\delta_{ik}$ for all
$i,k\in\\{0,\dots,n\\}$, where $\delta_{ik}$ denotes the Kronecker symbol.

Choosing $\psi=D_{i}h_{\varepsilon}(u)$ in the weak equation (1.10) for
$u_{i}$ and taking the sum $\sum_{i=0}^{n}$ gives

$\displaystyle\int_{\Omega}h_{\varepsilon}(u(T))\,\differential
x-\int_{\Omega}h_{\varepsilon}(u(0))\,\differential
x=-\int_{0}^{T}\\!\int_{\Omega}\sum_{i,j,k}D_{ik}h_{\varepsilon}(u)A_{ij}(u)\nabla
u_{j}\cdot\nabla u_{k}\,\differential x\differential t.\qquad$ (3.4)

To pass to the limit $\varepsilon\to 0$, we wish to apply the dominated
convergence theorem. In the terms on the left-hand side this is immediate
thanks to the bound $0\leq u_{i}\leq 1$. Concerning the integral on the right-
hand side, we need appropriate information about its integrand. For this
purpose, we compute for $i,j\in\\{0,\dots,n\\}$, using the formula (1.6) for
$\mathbb{M}$,

$\displaystyle\sum_{k}D_{ik}h_{\varepsilon}(u)A_{ij}(u)\nabla u_{j}\cdot\nabla
u_{k}$
$\displaystyle=\frac{1}{u_{i}+\varepsilon}\mathbb{M}_{ij}(u)\frac{1}{u_{j}}\nabla
u_{j}\cdot\nabla u_{i}$ (3.5)
$\displaystyle=\frac{1}{u_{i}+\varepsilon}u_{i}(\delta_{ij}\sum_{l=0}^{n}K_{il}u_{l}-K_{ij}u_{j})\frac{1}{u_{j}}\nabla
u_{j}\cdot\nabla u_{i}$ (3.6)
$\displaystyle=\delta_{ij}\frac{1}{u_{i}+\varepsilon}|\nabla
u_{i}|^{2}\sum_{l=0}^{n}K_{il}u_{l}-\frac{u_{i}}{u_{i}+\varepsilon}K_{ij}\nabla
u_{j}\cdot\nabla u_{i}.$ (3.7)

The terms in the last line are bounded above (in absolute values, pointwise
a.e.) by the integrable function $C\sum_{i}|\nabla\sqrt{u_{i}}|^{2}$ for some
$\varepsilon$-independent constant $C\in(0,\infty)$. Hence, the dominated
convergence theorem allows to infer identity (3.1) upon taking the limit
$\varepsilon\downarrow 0$ in equation (3.4). ∎

###  3.2. Weak-strong stability estimate

In the derivation of the stability estimate asserted in Theorem 1.3, we
roughly follow the technique in [12]. The size-exclusion effect, however,
induces some new difficulties requiring the new estimate (3.21) to handle the
present problem.

###### Proof of Theorem 1.3.

Let $\tilde{u}=(\tilde{u}_{0},\dots,\tilde{u}_{n})$ be a strong solution and
$u=(u_{0},\dots,u_{n})$ a weak solution with the properties as stated in
Theorem 1.3. We recall the relation $U=(u_{1},\dots,u_{n})$,
$\tilde{U}=(\tilde{u}_{1},\dots,\tilde{u}_{n})$ and note that for
$i\in\\{1,\dots,n\\}$

$\displaystyle
D_{i}\hat{h}(U)=\log\bigg(\frac{U_{i}}{u_{0}}\bigg{missing})=\log\big(\tfrac{u_{i}}{u_{0}}\big{missing}).$
(3.8)

In the proof below, we freely switch between the original and the transformed
variables $u,\tilde{u}$ and $U,\tilde{U}$ and choose the form that appears to
be more convenient for the specific argument in question. As explained in
Section 2.1, both forms are equivalent and the specific choice of variables
has no mathematical significance.

We then define

$\displaystyle\hat{h}_{\mathrm{rel}}(U,\tilde{U})$
$\displaystyle:=\hat{h}(U)-\sum_{i=1}^{n}D_{i}\hat{h}(\tilde{U})(U_{i}-\tilde{U}_{i})-\hat{h}(\tilde{U})$
(3.9)
$\displaystyle=h(u)-\sum_{k=0}^{n}D_{k}h(\tilde{u})(u_{k}-\tilde{u}_{k})-h(\tilde{u})=h_{\mathrm{rel}}(u,\tilde{u}).$
(3.10)

Since $\hat{h}$ is strictly convex uniformly in $\mathcal{D}$ and
$u_{0}=1-\sum_{i=1}^{n}u_{i}$, we have the lower bound

$\displaystyle\hat{h}_{\mathrm{rel}}(U,\tilde{U})\gtrsim|U-\tilde{U}|^{2}\sim|u-\tilde{u}|^{2}.$
(3.11)

Let now $T\in(0,T^{*})$. We will establish a weak-strong stability estimate of
the form

$\displaystyle\evaluated{\int_{\Omega}\hat{h}_{\mathrm{rel}}(U,\tilde{U})\,\differential
x}_{\tau=0}^{\tau=t}$ $\displaystyle\leq
C_{T,\tilde{u}}\int_{0}^{t}\\!\\!\int_{\Omega}\;\hat{h}_{\mathrm{rel}}(U,\tilde{U})\;\differential
x\differential\tau,\quad t\in(0,T),$ (3.12)

with
$C_{T,\tilde{u}}=C(\|\tilde{u}\|_{C^{0,1}([0,T]\times\overline{\Omega})},(\inf_{[0,T]\times\overline{\Omega}}\tilde{u}_{i})_{i=0}^{n})$,
upon which the asserted inequality (1.12) follows from Gronwall’s lemma.

To begin with, we compute the time evolution of
$\hat{h}_{\mathrm{rel}}(U,\tilde{U})$ using the entropy dissipation balance
(3.2) and the regularity of $U,\tilde{U}$ as well as integration by parts for
Sobolev–Bochner functions (cf. [21, Lemma 7.3]),

$\displaystyle\evaluated{\int_{\Omega}\hat{h}_{\mathrm{rel}}(U,\tilde{U})\,\differential
x}_{\tau=0}^{\tau=t}$
$\displaystyle=-\int_{0}^{t}\\!\\!\int_{\Omega}\hat{\mathscr{P}}(U)\,\differential
x\differential\tau-\int_{0}^{t}\frac{\differential}{\differential\tau}\\!\int_{\Omega}\sum_{i=1}^{n}D_{i}\hat{h}(\tilde{U})(U_{i}-\tilde{U}_{i})+\hat{h}(\tilde{U})\,\differential
x\,\differential\tau$ (3.13)
$\displaystyle=-\int_{0}^{t}\\!\\!\int_{\Omega}\hat{\mathscr{P}}(U)\,\differential
x\differential\tau-\int_{0}^{t}\\!\\!\int_{\Omega}\big{(}\sum_{i,k=1}^{n}D_{ik}\hat{h}(\tilde{U})\partial_{t}\tilde{U}_{k}(U_{i}-\tilde{U}_{i})\big{)}\,\differential
x\differential\tau$ (3.14) $\displaystyle\hskip
108.405pt-\int_{0}^{t}\sum_{i=1}^{n}\langle\partial_{t}U_{i},D_{i}\hat{h}(\tilde{U})\rangle\,\differential\tau.$
(3.15)

For the penultimate integral on the right-hand side we next use the fact that
$\tilde{u}$ is a strong solution, while the last integral on the right-hand
side is further rewritten appealing to the weak solution property (2.20) of
$U$. Upon insertion of the identity $\hat{A}=\mathbb{\hat{M}}D^{2}\hat{h}$,
some rearrangement of terms and renaming of indices, we arrive at

$\begin{split}\quad\evaluated{\int_{\Omega}\hat{h}_{\mathrm{rel}}(U,\tilde{U})\,\differential
x}_{\tau=0}^{\tau=t}&=\int_{0}^{t}\\!\\!\int_{\Omega}\;\rho\;\differential
x\differential\tau,\quad t\in(0,T),\end{split}$ (3.16)

where444A caveat is in place here: Since the function $D_{l}\hat{h}(U)$ may
not be locally integrable (not to mention weakly differentiable in the Sobolev
sense), the gradient $\nabla D_{l}\hat{h}(U)$ in (3.17) needs to be read in a
symbolic way with the rigorous meaning $\nabla
D_{l}\hat{h}(U):=\sum_{j=1}^{n}D_{j}D_{l}\hat{h}(U)\nabla U_{j}$.

$\begin{split}\rho:=-\hat{\mathscr{P}}(U)&+\sum_{i,l=1}^{n}\nabla
D_{i}\hat{h}(\tilde{U})\cdot\mathbb{\hat{M}}_{il}(U)\nabla D_{l}\hat{h}(U)\\\
&+\sum_{i,l,j=1}^{n}\nabla\big{(}D_{ij}\hat{h}(\tilde{U})(U_{j}-\tilde{U}_{j})\big{)}\cdot\mathbb{\hat{M}}_{il}(\tilde{U})\nabla
D_{l}\hat{h}(\tilde{U})\end{split}$ (3.17)

with the understanding that $\nabla
D_{l}\hat{h}(U):=\sum_{k=1}^{n}D_{lk}\hat{h}(U)\nabla U_{k}$.

Hence, for proving the stability estimate (3.12) it suffices to show the
pointwise bound

$\displaystyle\rho\leq C_{T,\tilde{u}}\,\hat{h}_{\mathrm{rel}}(U,\tilde{U}).$
(3.18)

For this purpose, we fix $\iota:=\iota_{T}\in(0,1)$ such that
$\tilde{u}_{i}\geq 2\iota$ in $[0,T]\times\Omega$ for all
$i\in\\{0,\dots,n\\}$. For a.e. $(t,x)\in[0,T]\times\Omega$ we then
distinguish the following cases:

  * •

Case 1: $\min\\{u_{0}(t,x),\dots,u_{n}(t,x)\\}>\iota$.

  * •

Case 2: $\min\\{u_{0}(t,x),\dots,u_{n}(t,x)\\}\leq\iota$.

Case 1. In this case, $u_{i}\geq\iota$ for all $i\in\\{0,\dots,n\\}$ and there
exists a convex set $\mathcal{K}_{\iota}\subset\subset\mathcal{D}$ with
$\operatorname{dist}(\mathcal{K}_{\iota},\partial\mathcal{D})\gtrsim\iota$
such that $U(t,x)\in\mathcal{K}_{\iota}$. We will show that $\rho\leq
C_{T,\tilde{u}}|u-\tilde{u}|^{2}$ using the nondegeneracy property (2.19), the
strict convexity of $\hat{h}$, as well as the smoothness of $\mathbb{\hat{M}}$
and $\hat{h}$ in $\mathcal{D}$. Dependencies of constants on fixed quantities
such as $\|\tilde{u}\|_{C^{0,1}([0,T]\times\overline{\Omega})}$ and $\iota$
will frequently be omitted.

Recalling the identity $\hat{\mathscr{P}}(U)=\sum_{i,l=1}^{n}\nabla
D_{i}\hat{h}(U)\cdot\mathbb{\hat{M}}_{il}(U)\nabla D_{l}\hat{h}(U)$ and taking
advantage of the fact that, since $U(t,x)$ is away from $\partial\mathcal{D}$,
the Hessian $D^{2}\hat{h}(U(t,x))\in\mathbb{R}^{n\times n}$ is well-defined,
we may rewrite (3.17) in the form

$\displaystyle\begin{aligned}
\rho&=-\sum_{i,l=1}^{n}\nabla(D_{i}\hat{h}(U)-D_{i}\hat{h}(\tilde{U}))\cdot\mathbb{\hat{M}}_{il}(U)\nabla(D_{l}\hat{h}(U)-D_{l}\hat{h}(\tilde{U}))\\\
&\qquad-\sum_{i,l=1}^{n}\nabla(D_{i}\hat{h}(U)-D_{i}\hat{h}(\tilde{U}))\cdot(\mathbb{\hat{M}}_{il}(U)-\mathbb{\hat{M}}_{il}(\tilde{U}))\nabla
D_{l}\hat{h}(\tilde{U})\\\
&\qquad-\sum_{i,j,l=1}^{n}\nabla\big{(}D_{i}\hat{h}(U)-D_{i}\hat{h}(\tilde{U})-D_{ij}\hat{h}(\tilde{U})(U_{j}-\tilde{U}_{j})\big{)}\cdot\mathbb{\hat{M}}_{il}(\tilde{U})\nabla
D_{l}\hat{h}(\tilde{U}),\end{aligned}$ (3.19)

where we continue to employ the short-hand notation $\nabla
D_{i}\hat{h}(U):=\sum_{k=1}^{n}D_{ik}\hat{h}(U)\nabla U_{k}$.

Thanks to inequality (2.19) with $\mathcal{K}=\mathcal{K}_{\iota}$ and the
fact that $h\in C^{4}(\mathcal{K}_{\iota})$, there exists
$\epsilon_{*}=\epsilon_{*}(\mathcal{K}_{\iota})>0$ such that for any
$\delta\in(0,1)$

$\displaystyle\rho$ $\displaystyle\leq-\epsilon_{*}\sum_{i=1}^{n}|\nabla
D_{i}\hat{h}(U)-\nabla D_{i}\hat{h}(\tilde{U})|^{2}+\delta|\nabla
u-\nabla\tilde{u}|^{2}+C_{\delta}|u-\tilde{u}|^{2}.$ (3.20)

Here, the last term in (3.19) has been estimated invoking Taylor’s theorem,
see Appendix A for details.

We assert that there exists $\epsilon(\iota_{T})>0$ and
$C_{T,\tilde{u}}=C(\|\tilde{u}\|_{C^{0,1}([0,T]\times\overline{\Omega})},\iota_{T})<\infty$
such that

$G:=\sum_{i=1}^{n}|\nabla D_{i}\hat{h}(U)-\nabla
D_{i}\hat{h}(\tilde{U})|^{2}\geq\epsilon(\iota_{T})|\nabla
u-\nabla\tilde{u}|^{2}-C_{T,\tilde{u}}|u-\tilde{u}|^{2}.$ (3.21)

First suppose that (3.21) holds true. Then, choosing
$\delta=\epsilon_{*}\epsilon(\iota_{T})$ allows us to deduce (3.18) in Case 1.

Proof of inequality (3.21) (in Case 1). We compute for $i\in\\{1,\dots,n\\}$

$\displaystyle\nabla D_{i}\hat{h}(U)-\nabla D_{i}\hat{h}(\tilde{U})$
$\displaystyle=\sum_{j=1}^{n}\big{(}D_{ij}\hat{h}(U)\nabla
U_{j}-D_{ij}\hat{h}(\tilde{U})\nabla\tilde{U}_{j}\big{)}$ (3.22)
$\displaystyle=\sum_{j=1}^{n}\big{[}D_{ij}\hat{h}(U)(\nabla
U_{j}-\nabla\tilde{U}_{j})+(D_{ij}\hat{h}(U)-D_{ij}\hat{h}(\tilde{U}))\nabla\tilde{U}_{j}\big{]}.$
(3.23)

Since $\hat{h}$ is smooth and strictly convex on $\mathcal{D}$, the inverse
Hessian $B(U):=(D^{2}\hat{h}(U))^{-1}$ is well-defined and for
$U\in\mathcal{K}_{\iota}$ its spectral norm $|||B(U)|||$ is bounded above by a
finite constant $C(\iota)$ only depending on $\iota>0$. Hence,

$\displaystyle|\nabla U-\nabla\tilde{U}|$ $\displaystyle\leq|B(U)(\nabla
D\hat{h}(U)-\nabla
D\hat{h}(\tilde{U}))|+|B(U)(D^{2}\hat{h}(U)-D^{2}\hat{h}(\tilde{U}))\nabla\tilde{U}|$
(3.24) $\displaystyle\leq C(\iota)|\nabla D\hat{h}(U)-\nabla
D\hat{h}(\tilde{U})|+C_{T,\tilde{u}}|U-\tilde{U}|,$ (3.25)

which implies (3.21).

This concludes the proof of inequality (3.21).

Case 2. In this case, we recall the definition of $\rho$ in (3.17) and
Corollary 2.4 (see also inequality (2.18)) to estimate

$\displaystyle\rho$
$\displaystyle\leq-\epsilon_{1}\sum_{i=1}^{n}|\nabla\sqrt{u_{i}}|^{2}+C_{T,\tilde{u}}|\nabla
u|+C_{T,\tilde{u}}$ (3.26)
$\displaystyle\leq-\epsilon_{2}\sum_{i=1}^{n}|\nabla\sqrt{u_{i}}|^{2}+C_{T,\tilde{u}},$
(3.27)

where we have used the boundedness of $u$. Since
$|u-\tilde{u}|\gtrsim_{\iota}1$, we deduce (3.18).

The proof of Theorem 1.3 is now complete. ∎

######  Remark 3.2 (Second proof of inequality (3.21)).

For a somewhat more quantitative estimate in Case 1, the following alternative
proof of inequality (3.21) may be of interest. Let $T\in(0,\infty)$ be fixed
but arbitrary and abbreviate $\iota:=\iota_{T}$. In view of (3.8), we have

$\displaystyle
G=\sum_{i=1}^{n}|\nabla\log(\tfrac{u_{i}}{u_{0}})-\nabla\log(\tfrac{\tilde{u}_{i}}{\tilde{u}_{0}})|^{2},$
(3.28)

which can be rewritten as

$\displaystyle G=\sum_{i=1}^{n}|\nabla\varphi_{i}-\nabla\varphi_{0}|^{2}$
(3.29)

for $\varphi_{i}:=\log\big(\tfrac{u_{i}}{\tilde{u}_{i}}\big{missing})$. We
compute using the identity $u_{0}=1-\sum_{i=1}^{n}u_{i}$

$\displaystyle\nabla\varphi_{0}=\nabla\log(u_{0})-\nabla\log(\tilde{u}_{0})$
$\displaystyle=-\sum_{i=1}^{n}\tfrac{1}{u_{0}}\nabla
u_{i}+\sum_{i=1}^{n}\tfrac{1}{\tilde{u}_{0}}\nabla\tilde{u}_{i}$ (3.30)
$\displaystyle=-\sum_{i=1}^{n}\tfrac{u_{i}}{u_{0}}\nabla\log(u_{i})+\sum_{i=1}^{n}\tfrac{\tilde{u}_{i}}{\tilde{u}_{0}}\nabla\log(\tilde{u}_{i})$
(3.31)
$\displaystyle=-\sum_{i=1}^{n}\tfrac{u_{i}}{u_{0}}\nabla\varphi_{i}+R(u,\tilde{u}),$
(3.32)

where here and below $R(u,\tilde{u})$ denotes a harmless term that satisfies
$|R(u,\tilde{u})|\leq C_{T,\tilde{u}}|u-\tilde{u}|$ and may change from line
to line. Hence,

$\displaystyle\nabla\varphi_{1}-\nabla\varphi_{0}$
$\displaystyle=(1+\tfrac{u_{1}}{u_{0}})\nabla\varphi_{1}+\sum_{i=2}^{n}\tfrac{u_{i}}{u_{0}}\nabla\varphi_{i}+R(u,\tilde{u})$
(3.33)
$\displaystyle=(1+\sum_{i=1}^{n}\tfrac{u_{i}}{u_{0}})\nabla\varphi_{1}+\sum_{i=2}^{n}\tfrac{u_{i}}{u_{0}}(\nabla\varphi_{i}-\nabla\varphi_{1})+R(u,\tilde{u})$
(3.34)
$\displaystyle=(1+\sum_{i=1}^{n}\tfrac{u_{i}}{u_{0}})\nabla\varphi_{1}+\sum_{i=2}^{n}\tfrac{u_{i}}{u_{0}}(\nabla\varphi_{i}-\nabla\varphi_{0}+\nabla\varphi_{0}-\nabla\varphi_{1})+R(u,\tilde{u}).$
(3.35)

Using the inequality $2|a-b|^{2}\geq a^{2}-2b^{2}$, we infer for any
$\delta\in(0,1)$

$|\nabla\varphi_{1}-\nabla\varphi_{0}|^{2}\geq(1-2\delta)|\nabla\varphi_{1}-\nabla\varphi_{0}|^{2}+\delta(1+\tfrac{1-u_{0}}{u_{0}})|\nabla\varphi_{1}|^{2}\\\
-C_{1}(\iota)\delta\sum_{i=1}^{n}|\nabla\varphi_{i}-\nabla\varphi_{0}|^{2}-C_{T,\tilde{u}}|u-\tilde{u}|^{2}.$
(3.36)

For $\delta=\delta(\iota)>0$ small enough, we deduce

$\displaystyle G=\sum_{i=1}^{n}|\nabla\varphi_{i}-\nabla\varphi_{0}|^{2}$
$\displaystyle\geq\epsilon_{1}(\iota)|\nabla\varphi_{1}|^{2}-C_{T,\tilde{u}}|u-\tilde{u}|^{2}$
(3.37) $\displaystyle\geq\epsilon_{2}(\iota)|\nabla
u_{1}-\nabla\tilde{u}_{1}|^{2}-C_{T,\tilde{u}}|u-\tilde{u}|^{2},$ (3.38)

where the last estimate follows from the definition of $\varphi_{1}$ and
elementary manipulations. To infer the full bound (3.21) one argues by
symmetry.

###  3.3. Exponential stability of steady states

###### Proof of Theorem 1.4.

By Proposition 3.1, the weak solution $u$ satisfies identity (3.2). Moreover,
equation (1.10) implies that $\int u_{i}(t)\,\differential
x=\int\overline{u}_{i}\,\differential x$ for all $t\geq 0$ and
$i\in\\{0,\dots,n\\}$. Hence, we have for all $0\leq s<t<\infty$ and a
positive constant $\epsilon=\epsilon(\Omega)$

$\displaystyle\int_{\Omega}h_{\mathrm{rel}}(u(t),\overline{u})\,\differential
x-\int_{\Omega}h_{\mathrm{rel}}(u(s),\overline{u})\,\differential x$
$\displaystyle=-\int_{s}^{t}\\!\int_{\Omega}\mathscr{P}(u)\,\differential
x\differential\tau$ (3.39)
$\displaystyle\leq-4\kappa\int_{s}^{t}\\!\int_{\Omega}\sum_{i=0}^{n}|\nabla\sqrt{u_{i}}|^{2}\,\differential
x\differential\tau$ (3.40)
$\displaystyle\leq-\epsilon\kappa\int_{s}^{t}\\!\int_{\Omega}h_{\mathrm{rel}}(u(\tau),\overline{u})\,\differential
x\differential\tau,$ (3.41)

where the second step uses inequality (2.28) in Lemma 2.3 and the last
estimate follows from the logarithmic Sobolev inequality applied pointwise in
time. A suitable version of Gronwall’s inequality now implies that

$\displaystyle\int_{\Omega}h_{\mathrm{rel}}(u(t),\overline{u})\,\differential
x\leq\int_{\Omega}h_{\mathrm{rel}}(u(0),\overline{u})\,\differential
x\,\exp(-\epsilon\kappa t)$ (3.42)

for all $t\geq 0$. By a classical Csiszár–Kullback type inequality (see e.g.
[23]), the left-hand side is bounded below by
$\frac{1}{C}\|u(t)-\overline{u}\|_{L^{1}(\Omega)}^{2}$ for some
$C=C(\Omega)\in(0,\infty)$. This proves the asserted inequality (1.13). ∎

##  4\. Existence and time asymptotics for partial interactions

In this section, we perform a global weak existence analysis for systems where
not necessarily all species are interacting. Except for the case of full
interactions, which is covered by [14, Theorem 2] (cf. Theorem 2.5), no
systematic existence analysis is available in the literature for system (1.1).
Here, we take advantage of this result for full interactions, which easily
gives rise to a family of approximate solutions to more general, partially
interacting systems, and then perform a vanishing-interaction limit. The weak
solutions obtained upon this construction further enjoy certain entropy
dissipation inequalities that allow us to prove convergence to equilibrium as
$t\to\infty$. Since the gradient bounds obtained from entropy dissipation may
be quite poor, the asymptotic analysis involves some non-standard arguments
when it comes to specifying a weak formulation of the equations and
determining the long-time behaviour of solutions.

###  4.1. Definitions and results

In this section, we only assume in addition to (H1), (H2) the following
hypothesis:

  1. (H3)

There exists $i_{0}\in\\{0,\dots,n\\}$ such that $K_{i_{0}j}>0$ for all
$j\not=i_{0}$.

Hypothesis (H3) is quite natural from a modelling point of view, and one may
think of the species $X_{i_{0}}$ as representing the vacancies.

Depending on their interaction properties $\\{K_{ij}\\}$, the species
$X_{0},\dots,X_{n}$ differ in their (mathematical) character. Our analysis
suggests classifying them as follows.

######  Definition 4.1 (Types $\mathfrak{A}$, $\mathfrak{B}$ and
$\mathfrak{C}$).

Let $\\{K_{ij}\\}_{i,j=0}^{n}$ satisfy (H3) and let $i\in\\{0,\dots,n\\}$.

  * •

We call the species $X_{i}$ of type $\mathfrak{A}$ if $K_{ij}>0$ for all
$j\not=i$.

  * •

Let the species $X_{i}$ not be of type $\mathfrak{A}$. We say that $X_{i}$ is
of type $\mathfrak{B}$ if $K_{i\ell}=0$ for all species $X_{\ell}$,
$\ell\not=i$, that are not of type $\mathfrak{A}$.

  * •

A species that is neither of type $\mathfrak{A}$ nor of type $\mathfrak{B}$
will be referred to as type $\mathfrak{C}$.

We informally write $i\in\mathfrak{A}$ to express the situation that $X_{i}$
is of type $\mathfrak{A}$ and use analogous notations for types $\mathfrak{B}$
and $\mathfrak{C}$. Moreover, the symbol $\alpha$ will be reserved as an index
for species of type $\mathfrak{A}$, and we frequently abbreviate

$\sum_{\alpha\in\mathfrak{A}}u_{\alpha}:=\sum_{\\{\alpha:\,X_{\alpha}\,\text{is
of type }\mathfrak{A}\\}}u_{\alpha}.$

With this convention we have
$\\{0,\dots,n\\}=\mathfrak{A}\,\dot{\cup}\,\mathfrak{B}\,\dot{\cup}\,\mathfrak{C}$.

In terms of regularity, species of type $\mathfrak{A}$ play a distinguished
role. In view of Lemma 2.3 we expect their densities to have gradients in
$L^{2}$.

We should further remark that hypothesis (H3) does not really make sense if
one admits initial data $u^{\mathrm{in}}$ satisfying
$u^{\mathrm{in}}_{i_{0}}\equiv 0$. Indeed, in this case mass conservation and
non-negativity enforce that $u_{i_{0}}\equiv 0$ and problem (1.1) effectively
reduces to an $n$-species system that does not necessarily conform to (H3) any
longer. More generally, the present approach to cross-diffusion system (1.1)
in the (more degenerate) partially interacting regime primarily aims at
understanding flows emanating from suitably regular initial data satisfying in
particular
$\inf_{\Omega}\sum_{\alpha\in\mathfrak{A}}u^{\mathrm{in}}_{\alpha}>0$ or at
least $\sum_{\alpha\in\mathfrak{A}}\overline{u^{\mathrm{in}}_{\alpha}}>0$; the
latter condition turns out to be sufficient for the generalised solutions to
be constructed below to relax to the steady state in the long-time limit (cf.
Theorem 4.4).

Let us now turn to the statement of our main result regarding existence for
general interactions. Loosely speaking, we obtain global existence of weak
solutions extending previous results in [6, 24], where only species of type
$\mathfrak{A}$ and $\mathfrak{B}$ were considered and where
$\mathfrak{A}=\\{0\\}$. (The method of proof in the present paper is, however,
different.) For the precise sense, in which gradients appearing in Theorem 4.2
below are generally to be understood, we refer to Remark 4.3 following the
statement of this theorem. We further recall that
$C([0,\infty);L^{2}_{\mathrm{weak}}(\Omega))$ denotes the space of functions
$v:[0,\infty)\to L^{2}(\Omega)$ such that $t\mapsto v(t)$ is continuous with
respect to the weak topology on $L^{2}(\Omega)$.

######  Theorem 4.2 (Existence for general interactions).

Let $\\{K_{ij}\\}_{i,j=0}^{n}$ satisfy (H1), (H2) and (H3). Let further
$u^{\mathrm{in}}\in L^{1}(\Omega)^{1+n}$ satisfy
$u^{\mathrm{in}}(x)\in\mathcal{\overline{S}}$ for a.e. $x\in\Omega$. There
exists $u=(u_{i})_{i=0}^{n}\in L^{\infty}((0,\infty)\times\Omega)^{1+n}\cap
C([0,\infty);L^{2}_{\mathrm{weak}}(\Omega))^{1+n}$ with
$u(t,x)\in\mathcal{\overline{S}}$ for a.e. $(t,x)\in(0,\infty)\times\Omega$
enjoying the regularity

  * $\partial_{t}u_{i}\in L^{2}_{\mathrm{loc}}([0,\infty);H^{1}(\Omega)^{*})$ for all $i\in\\{0,\dots,n\\}$,

  * $\sqrt{u_{\alpha}}\in L^{2}_{\mathrm{loc}}([0,\infty);H^{1}(\Omega))$ for all $\alpha\in\mathfrak{A}$,

  * $\sqrt{u_{\alpha}u_{i}}\in L^{2}_{\mathrm{loc}}([0,\infty);H^{1}(\Omega))$ for all $\alpha\in\mathfrak{A}$ and all $i\in\\{0,\dots,n\\}$,

taking the initial data $u^{\mathrm{in}}$ and obeying for every
$i\in\\{0,\dots,n\\}$ the conservation law

$\displaystyle\int_{\Omega}u_{i}(t,x)\,\differential
x=\int_{\Omega}u^{\mathrm{in}}_{i}(x)\,\differential x\text{ for all }t>0,$
(4.1)

and satisfying the cross-diffusion system (1.1) in the following weak sense:

for all $i,k\in\\{0,\dots,n\\}$ such that $K_{ik}>0$, there exists
$\mathbb{Y}_{ik}\in L^{2}(0,\infty;L^{2}(\Omega))^{d}$ with
$\mathbb{Y}_{ik}=-\mathbb{Y}_{ki}$ and the properties that555See Remark 4.3
for the meaning of the terms $u_{k}\nabla u_{i}$.

$\displaystyle\mathbb{Y}_{ik}$ $\displaystyle=u_{k}\nabla u_{i}-u_{i}\nabla
u_{k}\qquad\text{ a.e.\ in
}\big{\\{}\sum_{\alpha\in\mathfrak{A}}u_{\alpha}>0\big{\\}},$ (4.2a) and if
$\\{i,k\\}\cap\mathfrak{A}\not=\emptyset$, $\displaystyle\mathbb{Y}_{ik}$
$\displaystyle=u_{k}\nabla u_{i}-u_{i}\nabla u_{k}\qquad\text{ a.e.\ in
}(0,\infty)\times\Omega,$ (4.2b)

such that for all $T<\infty$ and all $\psi\in
C^{\infty}([0,T]\times\overline{\Omega})$

$\displaystyle\int_{0}^{T}\\!\langle\partial_{t}u_{i},\psi\rangle\,\differential
t=-\int_{0}^{T}\\!\int_{\Omega}\sum_{k}K_{ik}\mathbb{Y}_{ik}\cdot\nabla\psi\,\differential
x\differential t.$ (4.3)

Furthermore, there exists a measurable set $J\subset(0,\infty)$ with
$\mathcal{L}^{1}(J)=0$ such that the following entropy dissipation
inequalities hold for all $T\in(0,\infty):$

$\displaystyle\mathcal{H}(u(T))+\int_{0}^{T}\\!\\!\int_{\Omega}\\!\mathscr{P}(u)\chi_{\\{\sum_{\alpha\in\mathfrak{A}}u_{\alpha}>0\\}}\,\differential
x\differential t\leq\mathcal{H}(u^{\mathrm{in}}),$ (4.4a)
$\displaystyle\begin{multlined}\mathcal{H}(u(T))+\int_{s}^{T}\\!\\!\int_{\Omega}\\!\mathscr{P}(u)\chi_{\\{\sum_{\alpha\in\mathfrak{A}}u_{\alpha}>0\\}}\,\differential
x\differential
t\leq\mathcal{H}(u(s))+C\mathcal{L}^{d}\big{(}\big{\\{}\sum_{\alpha\in\mathfrak{A}}u_{\alpha}(s)=0\big{\\}}\big{)}\\\
\quad\text{ for all }s\in(0,T)\setminus
J.\end{multlined}\mathcal{H}(u(T))+\int_{s}^{T}\\!\\!\int_{\Omega}\\!\mathscr{P}(u)\chi_{\\{\sum_{\alpha\in\mathfrak{A}}u_{\alpha}>0\\}}\,\differential
x\differential
t\leq\mathcal{H}(u(s))+C\mathcal{L}^{d}\big{(}\big{\\{}\sum_{\alpha\in\mathfrak{A}}u_{\alpha}(s)=0\big{\\}}\big{)}\\\
\quad\text{ for all }s\in(0,T)\setminus J.$ (4.4d)

Here, $\mathscr{P}(u)=\frac{1}{2}\sum_{i,j}K_{ij}u_{i}u_{j}|\frac{\nabla
u_{i}}{u_{i}}-\frac{\nabla u_{j}}{u_{j}}|^{2}$ in
${\\{\sum_{\alpha\in\mathfrak{A}}u_{\alpha}>0\\}}$ and $C<\infty$ is a fixed
constant.

######  Remark 4.3 (Generalised Sobolev gradients).

It is necessary to comment on the meaning of the gradients appearing in
identity (4.2a) (resp. (4.4)), since for $i,k\not\in\mathfrak{A}$ terms of the
form $u_{k}\nabla u_{i}$ (resp. $\sqrt{u_{k}}\nabla\sqrt{u_{i}}$) may not
define a distribution in the classical sense. Let us first note that the
regularity $\sqrt{u_{\alpha}}\in
L^{2}_{\mathrm{loc}}([0,\infty);H^{1}(\Omega))$ and $\sqrt{u_{\alpha}u_{i}}\in
L^{2}_{\mathrm{loc}}([0,\infty);H^{1}(\Omega))$ for $\alpha\in\mathfrak{A}$
ensures that the products $\sqrt{u_{\alpha}}\nabla\sqrt{u_{i}}$ and
$u_{\alpha}\nabla u_{i}$ are well-defined in the distributional sense and
agree with the $L^{2}_{t,x}$ functions
$\nabla(\sqrt{u_{\alpha}u_{i}})-\sqrt{u_{i}}\nabla\sqrt{u_{\alpha}}$ and
$\nabla(u_{\alpha}u_{i})-u_{i}\nabla u_{\alpha}$, respectively. This in turn
determines measurable functions ‘ $\nabla\sqrt{u_{i}}$ ’ and ‘ $\nabla u_{i}$
’ almost everywhere in the set $\\{u_{\alpha}>0\\}$ and uniquely in the
weighted Lebesgue space $L^{2}(u_{\alpha}\differential x\differential t)$. The
‘maximal’ set where for all $i\in\\{0,\dots,n\\}$ gradients
$\nabla\sqrt{u_{i}},$ $\nabla u_{i}$ can be defined in this sense is obtained
by repeating the above reasoning for the sums
$\sum_{\alpha\in\mathfrak{A}}\sqrt{u_{\alpha}}$ or
$\sum_{\alpha\in\mathfrak{A}}u_{\alpha}$, which allows us to define
$\nabla\sqrt{u_{i}}$ and $\nabla u_{i}$ in a pointwise sense
$\mathcal{L}^{1+d}$-a.e. in the set
$\\{\sum_{\alpha\in\mathfrak{A}}u_{\alpha}>0\\}$. Simple examples show that
under the general hypotheses of Theorem 4.2, we cannot expect gradients to
exist on a set larger than this. In fact, letting
$\mathfrak{B}\not=\emptyset$, $\mathfrak{C}=\emptyset$ and
$u^{\mathrm{in}}_{\alpha}\equiv 0$ for all $\alpha\in\mathfrak{A}$, the
function $u(t)\equiv u^{\mathrm{in}}\in L^{1}(\Omega)$ would be the natural
candidate for a solution. Yet, under the hypotheses of Theorem 4.2, the
distributional gradient $\nabla u^{\mathrm{in}}$ does not need to have a
pointwise meaning. If one focuses on the question of giving a meaning, on a
set (possibly) larger than $\\{\sum_{\alpha\in\mathfrak{A}}u_{\alpha}>0\\}$,
to _products_ of the form $u_{k}\nabla u_{i}$ for $i,k\in\mathfrak{C}$, the
problem becomes more delicate and is out of the scope of this manuscript.
Here, we content ourselves with the fact that, as we will see in Theorem 4.4
below, the information contained in Theorem 4.2 suffices to uniquely determine
the long-time asymptotic behaviour of such weak solutions, provided they
emanate from data $u^{\mathrm{in}}$ for which
$\sum_{\alpha\in\mathfrak{A}}u^{\mathrm{in}}_{\alpha}$ is non-trivial.

Let us note that the vector fields $\mathbb{Y}_{ik}$ in Theorem 4.2 will be
obtained from a compactness argument (cf. (4.31)) and that, as discussed in
Remark 4.3, when $\\{i,k\\}\cap\mathfrak{A}=\emptyset$, a direct description
of $\mathbb{Y}_{ik}$ in terms of $u$ is only available in the set
$\\{\sum_{\alpha\in\mathfrak{A}}u_{\alpha}>0\\}$. Observe, however, that for
$i\in\mathfrak{A}\cup\mathfrak{B}$ identity (4.2b) implies that

$\displaystyle\sum_{k=0}^{n}K_{ik}\mathbb{Y}_{ik}=\sum_{j=0}^{n}A_{ij}(u)\nabla
u_{j}\quad\text{ a.e.\ in }(0,\infty)\times\Omega$ (4.5)

with $A_{ij}(u)$ as in (1.5), so that in this case eq. (4.3) agrees with the
classical weak form of the $i$-th component of (1.1).

######  Theorem 4.4 (Convergence to equilibrium).

Assume the hypotheses of Theorem 4.2 and let $u\in
L^{\infty}((0,\infty)\times\Omega)^{1+n}\cap
C([0,\infty);L^{2}_{\mathrm{weak}}(\Omega))^{1+n}$ be an entropy-dissipating
weak solution of (1.1) in the sense that it enjoys the properties listed in
Theorem 4.2. Suppose further that
$\sum_{\alpha\in\mathfrak{A}}\overline{u^{\mathrm{in}}_{\alpha}}>0$ with
$\overline{\,\cdot\,}$ denoting the average
$\frac{1}{|\Omega|}\int_{\Omega}\,\cdot\,$. Then

$\displaystyle\lim_{t\to\infty}\|u(t)-\overline{u}\|_{L^{1}(\Omega)}=0,$ (4.6)

where $\overline{u}\equiv\overline{u^{\mathrm{in}}}$ is the steady state
associated with the initial data $u^{\mathrm{in}}$.

See Section 4.3 for the proof of Theorem 4.4. Note that in models without
$\mathfrak{C}$-type species the vanishing of
$\sum_{\alpha\in\mathfrak{A}}\overline{u^{\mathrm{in}}_{\alpha}}$ means that a
trivial ‘solution’ of (1.1) can be obtained by letting $u\equiv
u^{\mathrm{in}}$.

###  4.2. Construction via vanishing-interaction limit

The purpose of this section is to prove Theorem 4.2.

Observe that $\mathfrak{A}$ is non-empty by hypothesis (H3), and abbreviate

$\displaystyle K:=\min_{\alpha\in\mathfrak{A}}\min_{i:i\not=\alpha}K_{\alpha
i}>0.$ (4.7)

Our construction is based on taking a vanishing-interaction limit in the
family of models of full interaction type defined by
$K_{ij}^{\varepsilon}:=\max\\{K_{ij},\varepsilon\\}$ for all $0\leq
i\not=j\leq n$ (and as before $K_{ii}^{\varepsilon}=0$ for all $i$), where
$0<\varepsilon\ll 1$ is a small parameter which we assume to satisfy
$\varepsilon\leq\min\\{K_{ij}:K_{ij}>0\\}$. Then, by Theorem 2.5, there exists
a global weak solution $u^{\varepsilon}$ of system (1.1) with parameters
$\\{K_{ij}\\}:=\\{K_{ij}^{\varepsilon}\\}$ satisfying
$u^{\varepsilon}(t,x)\in\mathcal{\overline{S}}$ for a.e. $(t,x)$ and enjoying
the extra regularity (cf. (2.36))

$\displaystyle\sqrt{u^{\varepsilon}_{i}}\in
L^{2}_{\mathrm{loc}}([0,\infty);H^{1}(\Omega))\quad\text{ for all
}i\in\\{0,\dots,n\\}.$ (4.8)

Thanks to Proposition 3.1 and Lemma 2.2, $u^{\varepsilon}$ obeys the entropy
dissipation balance

$\displaystyle\int_{\Omega}h(u^{\varepsilon}(T))\,\differential
x+\int_{s}^{T}\\!\int_{\Omega}\mathscr{P}^{\varepsilon}(u^{\varepsilon})\,\differential
x\differential t=\int_{\Omega}h(u^{\varepsilon}(s))\,\differential x$ (4.9)

for all $0\leq s<T<\infty$, where

$\displaystyle\mathscr{P}^{\varepsilon}(u^{\varepsilon})=\frac{1}{2}\sum_{i,k}K_{ik}^{\varepsilon}u^{\varepsilon}_{i}u^{\varepsilon}_{k}\big{|}\tfrac{\nabla
u^{\varepsilon}_{i}}{u^{\varepsilon}_{i}}-\tfrac{\nabla
u^{\varepsilon}_{k}}{u^{\varepsilon}_{k}}\big{|}^{2}=2\sum_{i,k}K_{ik}^{\varepsilon}\big{|}\sqrt{u^{\varepsilon}_{k}}\nabla\sqrt{u^{\varepsilon}_{i}}-\sqrt{u^{\varepsilon}_{i}}\nabla\sqrt{u^{\varepsilon}_{k}}\big{|}^{2}.$
(4.10)

Observe that, as a consequence of Lemma 2.3,

$\displaystyle\mathscr{P}^{\varepsilon}(u^{\varepsilon})$
$\displaystyle\gtrsim\sum_{\alpha\in\mathfrak{A}}|\nabla\sqrt{u^{\varepsilon}_{\alpha}}|^{2}+\sum_{\alpha\in\mathfrak{A}}\sum_{i}u^{\varepsilon}_{\alpha}|\nabla\sqrt{u^{\varepsilon}_{i}}|^{2}+\varepsilon\sum_{i}|\nabla\sqrt{u^{\varepsilon}_{i}}|^{2}$
(4.11)
$\displaystyle\qquad+\sum_{i,k}K_{ik}^{\varepsilon}\big{|}\sqrt{u^{\varepsilon}_{k}}\nabla\sqrt{u^{\varepsilon}_{i}}-\sqrt{u^{\varepsilon}_{i}}\nabla\sqrt{u^{\varepsilon}_{k}}\big{|}^{2}+\sum_{\alpha\in\mathfrak{A}}\sum_{i}|\nabla\sqrt{u^{\varepsilon}_{\alpha}u^{\varepsilon}_{i}}|^{2},$
(4.12)

where the bound for
$\sum_{\alpha\in\mathfrak{A}}\sum_{i}|\nabla\sqrt{u^{\varepsilon}_{\alpha}u^{\varepsilon}_{i}}|^{2}$
follows a posteriori by estimating
$|\nabla\sqrt{u^{\varepsilon}_{\alpha}u^{\varepsilon}_{i}}|^{2}\lesssim|\nabla\sqrt{u^{\varepsilon}_{\alpha}}|^{2}+u^{\varepsilon}_{\alpha}|\nabla\sqrt{u^{\varepsilon}_{i}}|^{2}$.
Since $h\geq 0$ and $\int_{\Omega}h(u^{\varepsilon}(0))\,\differential
x=\int_{\Omega}h(u^{\mathrm{in}})\,\differential x<\infty$, we infer the
following $(\varepsilon,T)$-uniform bounds:

$\displaystyle\begin{aligned}
\int_{0}^{T}\\!\int_{\Omega}\bigg{(}\;&\sum_{\alpha\in\mathfrak{A}}|\nabla\sqrt{u^{\varepsilon}_{\alpha}}|^{2}+\sum_{\alpha\in\mathfrak{A}}\sum_{i}u^{\varepsilon}_{\alpha}|\nabla\sqrt{u^{\varepsilon}_{i}}|^{2}+\sum_{\alpha\in\mathfrak{A}}\sum_{i}|\nabla\sqrt{u^{\varepsilon}_{\alpha}u^{\varepsilon}_{i}}|^{2}\\\
&+\sum_{i,k}K_{ik}^{\varepsilon}\big{|}\sqrt{u^{\varepsilon}_{k}}\nabla\sqrt{u^{\varepsilon}_{i}}-\sqrt{u^{\varepsilon}_{i}}\nabla\sqrt{u^{\varepsilon}_{k}}\big{|}^{2}+\varepsilon\sum_{i}|\nabla\sqrt{u^{\varepsilon}_{i}}|^{2}\bigg{)}\,\differential
x\differential t\leq C<\infty.\end{aligned}$ (4.13)

Using the $\varepsilon$-uniform control of the term
$\sqrt{K_{ik}^{\varepsilon}}(\sqrt{u^{\varepsilon}_{k}}\nabla\sqrt{u^{\varepsilon}_{i}}-\sqrt{u^{\varepsilon}_{i}}\nabla\sqrt{u^{\varepsilon}_{k}})$
in $L^{2}_{t,x}$ in (4.13) and the weak equation satisfied by
$u^{\varepsilon}_{i}$, we further deduce the $\varepsilon$-uniform bound

$\displaystyle\|\partial_{t}u^{\varepsilon}_{i}\|_{L^{2}(0,T;H^{1}(\Omega)^{*})}\leq
C<\infty$ (4.14)

for all $i\in\\{0,\dots,n\\}$.

Thanks to the above estimates and the control
$u^{\varepsilon}(t,x)\in\mathcal{\overline{S}}$, there exists a subsequence
$\varepsilon\downarrow 0$ (not relabelled) such that for all
$i\in\\{0,\dots,n\\}$ and suitable $u_{i}\in
L^{\infty}((0,\infty)\times\Omega)$, $0\leq u_{i}\leq 1$, satisfying
$\sqrt{u_{\alpha}}\in L^{2}_{\mathrm{loc}}([0,\infty);H^{1}(\Omega))$ whenever
$\alpha\in\mathfrak{A}$ the following convergence properties hold:

$\displaystyle u^{\varepsilon}_{i}\rightharpoonup u_{i}\quad\text{ in
}L^{2}_{\mathrm{loc}}([0,\infty);L^{2}(\Omega)),$ (4.15)
$\displaystyle\partial_{t}u^{\varepsilon}_{i}\overset{*}{\rightharpoonup}\partial_{t}u_{i}\quad\text{
in }L^{2}_{\mathrm{loc}}([0,\infty);H^{1}(\Omega)^{*})$ (4.16) $\displaystyle
u^{\varepsilon}_{i}(t)\rightharpoonup u_{i}(t)\quad\text{ in
}L^{2}(\Omega),\quad\text{ locally uniformly in }t\in[0,\infty),$ (4.17)

where the last line exploits the fact that for every $T<\infty$, the family
$\\{u^{\varepsilon}_{i}\\}$ is $\varepsilon$-uniformly bounded in
$L^{\infty}(0,T;L^{2}(\Omega))\cap H^{1}(0,T;H^{1}(\Omega)^{*})$, giving
compactness in $C([0,T];L^{2}_{\mathrm{weak}}(\Omega))$ (cf. [17, Lemma C.1]).

Additionally, the subsequence $\varepsilon\downarrow 0$ will be chosen such
that for all $\alpha\in\mathfrak{A}$

$\displaystyle\begin{aligned}
&\sqrt{u^{\varepsilon}_{\alpha}}\rightharpoonup\sqrt{u_{\alpha}}\quad\text{ in
}L^{2}_{\mathrm{loc}}([0,\infty);H^{1}(\Omega)),\\\
&u^{\varepsilon}_{\alpha}\to u_{\alpha}\quad\text{ a.e.\ in
}(0,\infty)\times\Omega.\end{aligned}$ (4.18)

Since the constraint $\sum_{i=0}^{n}u^{\varepsilon}_{i}=1$ is linear, the weak
convergence (4.15) implies that $\sum_{i=0}^{n}u_{i}(t,x)=1$ for a.e.
$(t,x)\in(0,\infty)\times\Omega$. Line (4.15) further yields
$\int_{0}^{T}\\!\\!\int_{\Omega}u_{i}\varphi\,\differential x\differential
t\geq 0$ for all $\varphi\in C_{c}((0,T)\times\Omega)$ with $\varphi\geq 0$;
hence $u_{i}\geq 0$ a.e., allowing us to conclude that
$u(t,x)\in\mathcal{\overline{S}}$ for a.e. $(t,x)\in(0,\infty)\times\Omega$.

Next, we assert that for all $\alpha\in\mathfrak{A}$ and $i\in\\{0,\dots,n\\}$

$\displaystyle{u^{\varepsilon}_{\alpha}}u^{\varepsilon}_{i}\to
u_{\alpha}u_{i}\quad\text{ in
}L^{2}_{\mathrm{loc}}([0,\infty);L^{2}(\Omega)).$ (4.19)

The convergence (4.19) follows from the Aubin–Lions lemma applied to
${u^{\varepsilon}_{\alpha}}u^{\varepsilon}_{i}$. To verify the required
control on the time derivative, we compute in the sense of distributions

$\displaystyle\partial_{t}\big{(}{u^{\varepsilon}_{\alpha}}u^{\varepsilon}_{i}\big{)}=\partial_{t}u^{\varepsilon}_{\alpha}u^{\varepsilon}_{i}+\partial_{t}u^{\varepsilon}_{i}u^{\varepsilon}_{\alpha}.$
(4.20)

Let now $T<\infty$. For the first term on the right-hand side we compute for
$\psi\in C^{\infty}([0,T]\times\overline{\Omega})$ using the equation
satisfied by $u^{\varepsilon}_{\alpha}$

$\displaystyle\int_{0}^{T}\\!\langle\partial_{t}u^{\varepsilon}_{\alpha}u^{\varepsilon}_{i},\psi\rangle\,\differential
t=\int_{0}^{T}\\!\langle\partial_{t}u^{\varepsilon}_{\alpha},u^{\varepsilon}_{i}\psi\rangle\,\differential
t$ $\displaystyle=-\int_{0}^{T}\\!\int_{\Omega}\sum_{j}A_{\alpha
j}(u^{\varepsilon})\nabla
u^{\varepsilon}_{j}\cdot\nabla(u^{\varepsilon}_{i}\psi)\,\differential
x\differential t$ (4.21) $\displaystyle=-\int_{0}^{T}\\!\langle
F^{\varepsilon},\sqrt{u^{\varepsilon}_{\alpha}}\nabla(u^{\varepsilon}_{i}\psi)\rangle\,\differential
t,$ (4.22)

where
$F^{\varepsilon}:=\frac{1}{\sqrt{u^{\varepsilon}_{\alpha}}}\sum_{j}A_{\alpha
j}(u^{\varepsilon})\nabla u^{\varepsilon}_{j}\in L^{2}_{t,x}$. By (4.13) the
term $F^{\varepsilon}$ is $\varepsilon$-uniformly bounded in
$L^{2}(0,T;L^{2}(\Omega))$. Hence, recalling also the uniform bounds

$\displaystyle\|u^{\varepsilon}_{\alpha}\|_{L^{2}([0,T];H^{1}(\Omega))}+\|\sqrt{u^{\varepsilon}_{\alpha}}\nabla
u^{\varepsilon}_{i}\|_{L^{2}(0,T;L^{2}(\Omega))}+\|\partial_{t}u^{\varepsilon}_{i}\|_{L^{2}(0,T;H^{1}(\Omega)^{*})}\leq
C_{T}<\infty,$ (4.23)

we find that
$\partial_{t}\big{(}{u^{\varepsilon}_{\alpha}}u^{\varepsilon}_{i}\big{)}$ is
$\varepsilon$-uniformly bounded in

$L^{1}(0,T;(H^{1}\cap L^{\infty})(\Omega)^{*}).$

Observing further the $\varepsilon$-uniform bound of
$\nabla\big{(}{u^{\varepsilon}_{\alpha}}u^{\varepsilon}_{i}\big{)}$ in
$L^{2}_{\mathrm{loc}}([0,\infty);L^{2}(\Omega))$ (cf. (4.13)), we infer (4.19)
by appealing to the Aubin–Lions lemma (see e.g. [21, Lemma 7.7]).

Thus, upon extraction of a subsequence, we may henceforth assume that

$\displaystyle{u^{\varepsilon}_{\alpha}}u^{\varepsilon}_{i}\to
u_{\alpha}u_{i}\quad\text{a.e.\ in }(0,\infty)\times\Omega.$ (4.24)

We next assert that this convergence implies pointwise convergence a.e. of
$u^{\varepsilon}_{i}$ to $u_{i}$ in the set $V_{\alpha}:=\\{u_{\alpha}>0\\}$,
where for definiteness we let $u_{\alpha}(t),u_{i}(t)$ denote the precise
representatives of the corresponding Lebesgue classes. To show this assertion,
let us fix $N\subset(0,\infty)\times\Omega$ measurable with
$\mathcal{L}^{1+d}(N)=0$ in such a way that $u^{\varepsilon}_{\alpha}\to
u_{\alpha}$ and ${u^{\varepsilon}_{\alpha}}u^{\varepsilon}_{i}\to
u_{\alpha}u_{i}$ pointwise in $((0,\infty)\times\Omega)\setminus N$. In
particular, for every $(t,x)\in V_{\alpha}$, we have
$u^{\varepsilon}_{\alpha}(t,x)>0$ for $\varepsilon>0$ small enough, allowing
us to infer that
$u_{\alpha}u^{\varepsilon}_{i}=\frac{u_{\alpha}}{u^{\varepsilon}_{\alpha}}\cdot{u^{\varepsilon}_{\alpha}}u^{\varepsilon}_{i}\to
u_{\alpha}u_{i}$ pointwise in $V_{\alpha}\setminus N$, since for each
$(t,x)\in V_{\alpha}\setminus N$ the expression
$u^{\varepsilon}_{\alpha}u^{\varepsilon}_{i}\cdot\frac{u_{\alpha}}{u^{\varepsilon}_{\alpha}}$
is well-defined for small enough $\varepsilon$. We have thus shown that for
all $i\in\\{0,\dots,n\\}$

$\displaystyle u^{\varepsilon}_{i}\to u_{i}\text{ pointwise a.e.\ in
}\bigcup_{\alpha\in\mathfrak{A}}\\{u_{\alpha}>0\\}=\Big{\\{}\sum_{\alpha\in\mathfrak{A}}u_{\alpha}>0\Big{\\}}.$
(4.25)

######  Lemma 4.5.

For all $i,k\in\\{0,\dots,n\\}$ with $K_{ik}>0$ there exist vector fields
$\mathbb{W}_{ik},\mathbb{Y}_{ik}\in L^{2}(0,\infty;L^{2}(\Omega))^{d}$,
$\mathbb{W}_{ik}=-\mathbb{W}_{ki}$, $\mathbb{Y}_{ik}=-\mathbb{Y}_{ki}$ that
satisfy

$\displaystyle\mathbb{W}_{ik}$
$\displaystyle=\sqrt{u_{k}}\nabla\sqrt{u_{i}}-\sqrt{u_{i}}\nabla\sqrt{u_{k}}$
$\displaystyle\text{ a.e.\ in
}\big{\\{}\sum_{\alpha\in\mathfrak{A}}u_{\alpha}>0\big{\\}},$ (4.26)
$\displaystyle\mathbb{Y}_{ik}$ $\displaystyle=u_{k}\nabla u_{i}-u_{i}\nabla
u_{k}$ $\displaystyle\text{ a.e.\ in
}\big{\\{}\sum_{\alpha\in\mathfrak{A}}u_{\alpha}>0\big{\\}},$ (4.27)

and if $\\{i,k\\}\cap\mathfrak{A}\not=\emptyset$,

$\displaystyle\mathbb{W}_{ik}$
$\displaystyle=\sqrt{u_{k}}\nabla\sqrt{u_{i}}-\sqrt{u_{i}}\nabla\sqrt{u_{k}}$
$\displaystyle\text{ a.e.\ in }(0,\infty)\times\Omega,$ (4.28)
$\displaystyle\mathbb{Y}_{ik}$ $\displaystyle=u_{k}\nabla u_{i}-u_{i}\nabla
u_{k}$ $\displaystyle\text{ a.e.\ in }(0,\infty)\times\Omega,$ (4.29)

and are such that

$\displaystyle\sqrt{u^{\varepsilon}_{k}}\nabla\sqrt{u^{\varepsilon}_{i}}-\sqrt{u^{\varepsilon}_{i}}\nabla\sqrt{u^{\varepsilon}_{k}}\;\rightharpoonup\;$
$\displaystyle\mathbb{W}_{ik}$ $\displaystyle\text{ in
}L^{2}((0,\infty)\times\Omega),\qquad$ (4.30) $\displaystyle
u^{\varepsilon}_{k}\nabla u^{\varepsilon}_{i}-u^{\varepsilon}_{i}\nabla
u^{\varepsilon}_{k}\;\rightharpoonup\;$ $\displaystyle\mathbb{Y}_{ik}$
$\displaystyle\text{ in }L^{2}((0,\infty)\times\Omega).\qquad$ (4.31)

Moreover, for all $i\in\\{0,\dots,n\\}$ one has

$\displaystyle\varepsilon\nabla\sqrt{u^{\varepsilon}_{i}}\;\to\;0\quad\text{
in }L^{2}((0,\infty)\times\Omega).$ (4.32)

###### Proof.

Properties (4.25), (4.18) and (4.13) imply that for all
$\alpha\in\mathfrak{A}$

$\displaystyle\sqrt{u^{\varepsilon}_{i}}\nabla\sqrt{u^{\varepsilon}_{\alpha}}\;\rightharpoonup\;\sqrt{u_{i}}\nabla\sqrt{u_{\alpha}}\;\;\;\text{
in }L^{2}((0,\infty)\times\Omega)\text{ for all }i\in\\{0,\dots,n\\}.$ (4.33)

Furthermore, using (4.25) and (4.18), it is easy to see for all
$\alpha\in\mathfrak{A}$ that
$\sqrt{u^{\varepsilon}_{\alpha}}\nabla\sqrt{u^{\varepsilon}_{i}}$ converges to
$\sqrt{u_{\alpha}}\nabla\sqrt{u_{i}}$ in the sense of distributions and hence,
by (4.13),

$\displaystyle\sqrt{u^{\varepsilon}_{\alpha}}\nabla\sqrt{u^{\varepsilon}_{i}}\;\rightharpoonup\;\sqrt{u_{\alpha}}\nabla\sqrt{u_{i}}\;\;\;\text{
in }L^{2}((0,\infty)\times\Omega)\text{ for all }i\in\\{0,\dots,n\\}.$ (4.34)

Let us now prove (4.30), (4.26). Since $K_{ik}>0$, the term
$\mathbb{W}^{\varepsilon}_{ik}:=(\sqrt{u^{\varepsilon}_{k}}\nabla\sqrt{u^{\varepsilon}_{i}}-\sqrt{u^{\varepsilon}_{i}}\nabla\sqrt{u^{\varepsilon}_{k}})$
is $\varepsilon$-uniformly bounded in $L^{2}((0,\infty)\times\Omega)$. Hence,
by weak compactness in $L^{2}((0,\infty)\times\Omega)$, there exists
$\mathbb{W}_{ik}\in L^{2}((0,\infty)\times\Omega)$ such that along a
subsequence $\varepsilon\downarrow 0$,
$\mathbb{W}^{\varepsilon}_{ik}\rightharpoonup\mathbb{W}_{ik}$ in
$L^{2}((0,\infty)\times\Omega)$. Thus,

$\displaystyle\eta\big{(}\sum_{\alpha\in\mathfrak{A}}u^{\varepsilon}_{\alpha}\big{)}\mathbb{W}^{\varepsilon}_{ik}\rightharpoonup\eta\big{(}\sum_{\alpha\in\mathfrak{A}}u_{\alpha}\big{)}\mathbb{W}_{ik}$
(4.35)

for all $\eta\in C^{\infty}([0,1])$ with
$\operatorname{supp}\eta\subset(0,1]$. At the same time, for any such $\eta$
we have weakly in $L^{2}_{t,x}$ the convergence

$\eta\big{(}\sum_{\alpha\in\mathfrak{A}}u^{\varepsilon}_{\alpha}\big{)}\mathbb{W}^{\varepsilon}_{ik}=\tfrac{\eta\big{(}\sum_{\alpha\in\mathfrak{A}}u^{\varepsilon}_{\alpha}\big{)}}{\sum_{\alpha\in\mathfrak{A}}u^{\varepsilon}_{\alpha}}\sum_{\alpha\in\mathfrak{A}}(\sqrt{u^{\varepsilon}_{k}}u^{\varepsilon}_{\alpha}\nabla\sqrt{u^{\varepsilon}_{i}}-\sqrt{u^{\varepsilon}_{i}}u^{\varepsilon}_{\alpha}\nabla\sqrt{u^{\varepsilon}_{k}})\\\
\rightharpoonup\tfrac{\eta\big{(}\sum_{\alpha\in\mathfrak{A}}u_{\alpha}\big{)}}{\sum_{\alpha\in\mathfrak{A}}u_{\alpha}}\sum_{\alpha\in\mathfrak{A}}(\sqrt{u_{k}}u_{\alpha}\nabla\sqrt{u_{i}}-\sqrt{u_{i}}u_{\alpha}\nabla\sqrt{u_{k}})\\\
=\eta\big{(}\sum_{\alpha\in\mathfrak{A}}u_{\alpha}\big{)}(\sqrt{u_{k}}\nabla\sqrt{u_{i}}-\sqrt{u_{i}}\nabla\sqrt{u_{k}}),$
(4.36)

where we also used (4.34). Uniqueness of the weak limit in $L^{2}$ implies
that

$\displaystyle\eta\big{(}\sum_{\alpha\in\mathfrak{A}}u_{\alpha}\big{)}\mathbb{W}_{ik}=\eta\big{(}\sum_{\alpha\in\mathfrak{A}}u_{\alpha}\big{)}(\sqrt{u_{k}}\nabla\sqrt{u_{i}}-\sqrt{u_{i}}\nabla\sqrt{u_{k}}).$
(4.37)

Since $\eta\in C^{\infty}_{c}((0,1])$ was arbitrary, we infer that
$\mathbb{W}_{ik}=\sqrt{u_{k}}\nabla\sqrt{u_{i}}-\sqrt{u_{i}}\nabla\sqrt{u_{k}}$
a.e. in $\\{\sum_{\alpha\in\mathfrak{A}}u_{\alpha}>0\\}$. Note that the
antisymmetry of $\mathbb{W}$ follows from that of $\mathbb{W}^{\varepsilon}$.
If $\\{i,k\\}\cap\mathfrak{A}\not=\emptyset$, then (4.33), (4.34) and
uniqueness of the weak limit in $L^{2}$ immediately yield
$\mathbb{W}_{ik}=\sqrt{u_{k}}\nabla\sqrt{u_{i}}-\sqrt{u_{i}}\nabla\sqrt{u_{k}}$
a.e. in $(0,\infty)\times\Omega.$ This proves the first assertions (4.30),
(4.26), (4.28). The assertions (4.31), (4.27), (4.29) concerning
$\mathbb{Y}_{ik}$ can be shown along similar lines.

The convergence (4.32) follows from the fact that
$\sqrt{\varepsilon}\nabla\sqrt{u^{\varepsilon}_{i}}$ is
$\varepsilon$-uniformly bounded in $L^{2}_{t,x}$ (cf. (4.13)). ∎

We are now in a position to complete the proof of Theorem 4.2.

###### Proof of Theorem 4.2.

Let $u$ denote the limiting candidate constructed on page 4.15 (see
(4.15)–(4.18)). We have already seen that it enjoys the asserted regularity
properties and satisfies $u(t,x)\in\mathcal{\overline{S}}$ for a.e. $(t,x)$.

The convergence properties (4.16), (4.31) and (4.32) allow to pass to the
limit $\varepsilon\downarrow 0$ in the weak form satisfied by
$u^{\varepsilon}_{i}$, $i\in\\{0,\dots,n\\}$, for any test function $\psi$,
giving the asserted equation (4.3) with $\mathbb{Y}_{ik}=-\mathbb{Y}_{ki}$
satisfying (4.2a) and, if $\\{i,k\\}\cap\mathfrak{A}\not=\emptyset$, (4.2b)
(cf. (4.27), (4.29)).

We next show the conservation law (4.1). It is a consequence of the
convergence

$\displaystyle\partial_{t}u^{\varepsilon}_{i}\overset{*}{\rightharpoonup}\partial_{t}u_{i}\quad\text{
in }L^{2}_{\mathrm{loc}}([0,\infty);H^{1}(\Omega)^{*}),$ (4.38)

the fact that
$\int_{0}^{t}\langle\partial_{t}u^{\varepsilon}_{i},1\rangle\differential\tau\equiv
0$, which follows from the equation for $u^{\varepsilon}_{i}$, and the
identity666Identity (4.39) can be shown by approximation, using the density of
the set $C^{1}([0,T];L^{2}(\Omega))$ in the Sobolev-Bochner space $\\{v\in
L^{2}(0,T;L^{2}(\Omega)):\frac{\differential v}{\differential t}\in
L^{2}(0,T;H^{1}(\Omega)^{*})\\}$ equipped with the standard norm.

$\displaystyle\int_{0}^{t}\langle\partial_{t}u_{i},1\rangle\differential\tau=\int_{\Omega}u_{i}(t,x)\,\differential
x-\int_{\Omega}u_{i}(0,x)\,\differential x.$ (4.39)

It remains to prove the entropy dissipation inequalities (4.4). The proof of
inequality (4.4a) being somewhat simpler, we only provide details on how to
derive (4.4d), where we rely on the entropy balance law (4.9) for
$u^{\varepsilon}$. For passing to the limit in the entropy at the final time
$T$ it is convenient to use (4.17), which gives, thanks to the weak lower
semicontinuity in $L^{2}(\Omega)$ of the entropy functional,
$\mathcal{H}(u(T))\leq\liminf_{\varepsilon\to
0}\mathcal{H}(u^{\varepsilon}(T))$. Note that here we have specified $u=u(T)$
on a set of times of measure zero in such a way that (4.17) holds for all
$T\in[0,\infty)$. To deal with the entropy at times $s\in(0,T)$, we observe
that thanks to (4.25) there exists a measurable set $J\subset[0,\infty)$ with
$\mathcal{L}^{1}(J)=0$ such that for all $t\in[0,\infty)\setminus J$

$\displaystyle u^{\varepsilon}(t)\to u(t)\qquad\text{ a.e.\ in
}\big{\\{}\sum_{\alpha\in\mathfrak{A}}u_{\alpha}(t)>0\big{\\}}.$ (4.40)

Since $|u^{\varepsilon}|\lesssim 1$, Lebesgue’s dominated convergence theorem
thus gives for $s\not\in J$

$\displaystyle\limsup_{\varepsilon\to
0}\int_{\Omega}h(u^{\varepsilon}(s))\,\differential
x\leq\int_{\Omega}h(u(s))\,\differential
x+C\mathcal{L}^{d}\big{(}\big{\\{}\sum_{\alpha\in\mathfrak{A}}u_{\alpha}(s)=0\big{\\}}\big{)},$
(4.41)

where $C<\infty$ is a universal constant.

We are left with estimating below the liminf of the dissipation term. First
recall that

$\displaystyle\mathscr{P}^{\varepsilon}(u^{\varepsilon})\geq\mathscr{P}(u^{\varepsilon})=2\sum_{i,k}K_{ik}|\sqrt{u^{\varepsilon}_{k}}\nabla\sqrt{u^{\varepsilon}_{i}}-\sqrt{u^{\varepsilon}_{i}}\nabla\sqrt{u^{\varepsilon}_{k}}|^{2}=2\sum_{i,k}K_{ik}|\mathbb{W}_{ik}^{\varepsilon}|^{2},$
(4.42)

where for $K_{ik}>0$ we have
$\mathbb{W}_{ik}^{\varepsilon}:=\sqrt{u^{\varepsilon}_{k}}\nabla\sqrt{u^{\varepsilon}_{i}}-\sqrt{u^{\varepsilon}_{i}}\nabla\sqrt{u^{\varepsilon}_{k}}\rightharpoonup\mathbb{W}_{ik}$
in $L^{2}_{t,x}$ with $\mathbb{W}_{ik}$ as in Lemma 4.5. Thus, weak lower
semicontinuity of the norm in $L^{2}_{t,x}$ yields

$\displaystyle\int_{0}^{T}\\!\\!\int_{\Omega}\\!\mathscr{P}(u)\chi_{\\{\sum_{\alpha\in\mathfrak{A}}u_{\alpha}>0\\}}\,\differential
x\differential
t\leq\int_{0}^{T}\\!\\!\int_{\Omega}\\!2\sum_{i,k}K_{ik}|\mathbb{W}_{ik}|^{2}\,\differential
x\differential t\leq\liminf_{\varepsilon\to
0}\int_{0}^{T}\\!\\!\int_{\Omega}\\!\mathscr{P}^{\varepsilon}(u^{\varepsilon})\,\differential
x\differential t.$ (4.43)

The proof of Theorem 4.2 is now complete. ∎

###  4.3. Convergence to equilibrium

###### Proof of Theorem 4.4.

Below, the mass conservation property (4.1) will, in places, be used without
explicit reference. Let $J$ be such that (4.4) holds. Thanks to (4.4a) there
exists a sequence $\\{t_{k}\\}\subset[0,\infty)\setminus J$ with
$t_{k}\to\infty$ such that

$\displaystyle\lim_{k\to\infty}\int_{\Omega}\mathscr{P}(u(t_{k}))\chi_{\\{\sum_{\alpha\in\mathfrak{A}}u_{\alpha}(t_{k})>0\\}}\,\differential
x=0.$ (4.44)

Similarly to (4.11), we infer from Lemma 2.3 the uniform bound

$\displaystyle\mathscr{P}(u)\chi_{\\{\sum_{\alpha\in\mathfrak{A}}u_{\alpha}>0\\}}$
$\displaystyle\gtrsim\sum_{\alpha\in\mathfrak{A}}|\nabla\sqrt{u_{\alpha}}|^{2}+\sum_{\alpha\in\mathfrak{A}}\sum_{i}|\nabla\sqrt{u_{\alpha}u_{i}}|^{2}$
(4.45)
$\displaystyle\qquad+\sum_{i,l}K_{il}\big{|}\sqrt{u_{l}}\nabla\sqrt{u_{i}}-\sqrt{u_{i}}\nabla\sqrt{u_{l}}\big{|}^{2}\chi_{\\{\sum_{\alpha\in\mathfrak{A}}u_{\alpha}>0\\}},$
(4.46)

where we recall that the gradients in the last line are to be understood as
explained in Remark 4.3.

Hence, (4.44) implies that
$\lim_{k\to\infty}\|\nabla\sqrt{u_{\alpha}}(t_{k})\|_{L^{2}(\Omega)}=0$ for
all $\alpha\in\mathfrak{A}$, which further yields, by classical properties for
Sobolev weak derivatives,

$\displaystyle\sqrt{u_{\alpha}}(t_{k})\to\sqrt{\overline{u}_{\alpha}}\qquad\text{
in }H^{1}(\Omega).$ (4.47)

We assert that this further gives

$\displaystyle u_{i}(t_{k})\to\overline{u_{i}}\qquad\text{ in
}L^{2}(\Omega)\qquad\text{for all }i\in\\{0,\dots,n\\}.$ (4.48)

To show (4.48), pick $\hat{\alpha}\in\mathfrak{A}$ such that
$\overline{u}_{\hat{\alpha}}>0$ (such $\hat{\alpha}$ exists by hypothesis) and
observe that (4.44) combined with (4.45) implies the convergence

$\displaystyle\chi_{\\{u_{\hat{\alpha}}(t_{k})>0\\}}\big{(}\sqrt{u_{\hat{\alpha}}}\nabla\sqrt{u_{i}}-\sqrt{u_{i}}\nabla\sqrt{u_{\hat{\alpha}}}\big{)}|_{t_{k}}\;\overset{k\to\infty}{\to}\;0\qquad\text{in
}L^{2}(\Omega),$ (4.49)

and hence
$\chi_{\\{u_{\hat{\alpha}}(t_{k})>0\\}}\sqrt{u_{\hat{\alpha}}(t_{k})}\nabla\sqrt{u_{i}(t_{k})}=\sqrt{u_{\hat{\alpha}}(t_{k})}\nabla\sqrt{u_{i}(t_{k})}\to
0$ in $L^{2}(\Omega)$. Recalling that
$\sqrt{u_{\hat{\alpha}}(t_{k})}\nabla\sqrt{u_{i}(t_{k})}:=\nabla(\sqrt{u_{\hat{\alpha}}(t_{k})}\sqrt{u_{i}(t_{k})})-(\nabla\sqrt{u_{\hat{\alpha}}(t_{k})})\sqrt{u_{i}(t_{k})}$
(cf. Remark 4.3) and using again classical arguments for Sobolev weak
derivatives, we find that

$\displaystyle\sqrt{u_{\hat{\alpha}}(t_{k})u_{i}(t_{k})}\to
c\equiv\text{const}\qquad\text{ in }H^{1}(\Omega).$ (4.50)

An elementary argument similar to that leading to (4.25) finally gives, upon
possibly passing to a subsequence, $u_{i}(t_{k})\to\overline{u}_{i}$ a.e. in
$\Omega$, which implies (4.48).

Let now $\\{\tau_{j}\\}\subset(0,\infty)$ with $\tau_{j}\to\infty$ be
arbitrary. Then there exists a subsequence $\\{t_{k_{j}}\\}$ of $\\{t_{k}\\}$
(with $k_{j}$ not necessarily strictly increasing, but
$\lim_{j\to\infty}k_{j}=\infty$) and $j_{0}\in\mathbb{N}$ such that
$t_{k_{j}}<\tau_{j}$ for all $j\geq j_{0}$. From the strong entropy
dissipation inequality (4.4d) (and mass conservation) we infer for all
$\tau>0$ and all $s\in(0,\tau)\setminus J$

$\displaystyle\mathcal{H}_{\mathrm{rel}}(u(\tau),\overline{u})\leq\mathcal{H}_{\mathrm{rel}}(u(s),\overline{u})+C\mathcal{L}^{d}\big{(}\big{\\{}\sum_{\alpha\in\mathfrak{A}}u_{\alpha}(s)=0\big{\\}}\big{)}.$
(4.51)

Inserting $\tau:=\tau_{j}$ and $s:=t_{k_{j}}$ in inequality (4.51) and sending
$j\to\infty$ yields

$\displaystyle\lim_{j\to\infty}\mathcal{H}_{\mathrm{rel}}(u(\tau_{j}),\overline{u})=0.$
(4.52)

Here we used the strong convergence (4.48) and the hypothesis
$\sum_{\alpha\in\mathfrak{A}}\overline{u}_{\alpha}>0$ to ensure that the
right-hand side of (4.51) vanishes as $s=t_{k_{j}}\to\infty$. ∎

##  Appendix A An auxiliary estimate

The estimate asserted in the lemma below has previously been employed in [12,
Theorem 2.8, Case $\mathcal{A}_{+}$] in a similar context. For the purpose of
being self-contained, we here also recall its elementary proof.

######  Lemma A.1 (Cf. [12]).

Let $G\subset\mathbb{R}^{N}$ be convex, let $f\in C^{3}_{b}(G)$ and let $v\in
H^{1}(\Omega)^{N},\tilde{v}\in C^{0,1}(\overline{\Omega},G)$. Then, for a.e.
$x\in\Omega$ such that $v(x)\in G$

$\big{|}\nabla\big{(}f(v)-f(\tilde{v})-\sum_{j=1}^{N}D_{j}f(\tilde{v})(v_{j}-\tilde{v}_{j})\big{)}\big{|}\\\
\lesssim_{\|f\|_{C^{3}(G)}}|\nabla
v-\nabla\tilde{v}||v-\tilde{v}|+\|\nabla\tilde{v}\|_{L^{\infty}}|v-\tilde{v}|^{2}\quad\text{at
the point }x,$ (A.1)

where we abbreviated $\nabla f(v):=\sum_{k=1}^{N}D_{k}f(v)\nabla v_{k}$.

###### Proof.

Using the definition of $\nabla f(v)$ and rearranging terms, we may rewrite

$\displaystyle\nabla\big{(}f(v)-f(\tilde{v})-\sum_{j=1}^{N}D_{j}f(\tilde{v})(v_{j}-\tilde{v}_{j})\big{)}$
(A.2) $\displaystyle=\sum_{k=1}^{N}\big{(}D_{k}f(v)\nabla
v_{k}-D_{k}f(\tilde{v})\nabla\tilde{v}_{k}\big{)}-\sum_{j=1}^{N}D_{j}f(\tilde{v})(\nabla
v_{j}-\nabla\tilde{v}_{j})-\sum_{j,l=1}^{N}D_{lj}f(\tilde{v})(v_{j}-\tilde{v}_{j})\nabla\tilde{v}_{l}$
(A.3) $\displaystyle=\sum_{k=1}^{N}(D_{k}f(v)-D_{k}f(\tilde{v}))(\nabla
v_{k}-\nabla\tilde{v}_{k})$ (A.4) $\displaystyle\hskip
108.405pt+\sum_{l=1}^{N}\big{(}D_{l}f(v)-D_{l}f(\tilde{v})-\sum_{j=1}^{N}D_{j}D_{l}f(\tilde{v})(v_{j}-\tilde{v}_{j})\big{)}\nabla\tilde{v}_{l}.$
(A.5)

Notice that, in the last step, some indices have been relabelled. The
assertion now follows by means of the Cauchy–Schwarz inequality and an
application of Taylor’s theorem to $D_{l}f\in C^{2}_{b}(G)$,
$l\in\\{1,\dots,N\\}$. ∎

##  Appendix B Notations

For definiteness, we here provide some details concerning our notations, most
of which are quite classical.

  * •

The convention $K_{ii}:=0$ for all $i\in\\{0,\dots,n\\}$ is adopted throughout
the paper.

  * •

We will further adopt the convention that $K_{ij}\cdot Q\equiv 0$ whenever
$K_{ij}=0$, even if the quantity ‘$Q$’ is not well-defined mathematically.
This will simplify notations as it avoids having to restrict sums involving
$K_{ij}$ to pairs $\\{i,j\\}$ satisfying $K_{ij}>0$.

  * •

Summation: since we are dealing not only with $(n+1)$ but also with $n$
dimensional systems, we avoid the usual summation convention. Unspecified sums
involving a Latin summation index (such as $i,j,k$) always indicate summing
over all indices, that is, we let $\sum_{i}:=\sum_{i=0}^{n}$ etc. Likewise
$\sum_{i\not=j}:=\sum_{i\in\\{0,\dots,n\\}}\sum_{j\in\\{0,\dots,n\\}\setminus\\{i\\}}$
and $\sum_{j:j\not=i}:=\sum_{j\in\\{0,\dots,n\\}\setminus\\{i\\}}$.

  * •

We further let
$\min_{i\not=j}=\min_{i\in\\{0,\dots,n\\}}\min_{j\in\\{0,\dots,n\\}\setminus\\{i\\}}$
and $\min_{i:i\not=j}=\min_{i\in\\{0,\dots,n\\}\setminus\\{j\\}}$.

  * •

The notation $A(u)\nabla u$ for $u=u(x)\in\mathbb{R}^{1+n}$ and
$A(u)\in\mathbb{R}^{(1+n)\times(1+n)}$ will be used as a short hand for the
$d\times(1+n)$ matrix $(A(u)\nabla u)_{i}=\sum_{j}A_{ij}(u)\nabla u_{j}$.
Analogous short-hand notations will be adopted for the matrices
$\mathbb{\hat{M}}(U),\hat{A}(U)\in\mathbb{R}^{n\times n}$.

  * •

For an interval $I\subseteq[0,\infty)$, $C(I;L^{2}_{\mathrm{weak}}(\Omega))$
denotes the space of functions $v:I\to L^{2}(\Omega)$ such that $t\mapsto
v(t)$ is weakly continuous in $L^{2}(\Omega)$.

  * •

For non-negative quantities $A,B$, we write $A\lesssim B$ to indicate that
there exists a constant $C\in(0,\infty)$ that only depends on fixed parameters
such that $A\leq CB$. We define $A\gtrsim B$ by $B\lesssim A$, and $A\sim B$
by [$A\lesssim B$ and $A\gtrsim B$].

  * •

For a smooth manifold $M$ and $p\in M$, we denote by $T_{p}M$ the tangent
space of $M$ at $p$.

  * •

For a sufficiently smooth function $f=f(u),u\in\mathbb{R}^{N},$ and
$i\in\\{1,\dots,n\\}$, we write $D_{i}f:=\frac{\partial}{\partial u_{i}}f$.
Analogous notations are used if $f=f(u_{0},\dots,u_{N})$ etc. We further
abbreviate $D_{ik}f:=D_{k}D_{i}f.$

  * •

Unless specified otherwise, we abbreviate
$\overline{u}_{i}:=\frac{1}{|\Omega|}\int_{\Omega}u_{i}\,\differential x$ for
all $i$.

## References

  * [1] H. Amann.  Dynamic theory of quasilinear parabolic equations. II. Reaction-diffusion systems.  Differential Integral Equations, 3(1):13–75, 1990.
  * [2] A. Bakhta and V. Ehrlacher.  Cross-diffusion systems with non-zero flux and moving boundary conditions.  ESAIM Math. Model. Numer. Anal., 52(4):1385–1415, 2018.
  * [3] J. Berendsen, M. Burger, V. Ehrlacher, and J.-F. Pietschmann.  Uniqueness of strong solutions and weak-strong stability in a system of cross-diffusion equations.  J. Evol. Equ., 20(2):459–483, 2020.
  * [4] J. Berendsen, M. Burger, and J.-F. Pietschmann.  On a cross-diffusion model for multiple species with nonlocal interaction and size exclusion.  Nonlinear Analysis, 159:10–39, 2017.
  * [5] M. Braukhoff, C. Raithel, and N. Zamponi.  Partial Hölder regularity for solutions of a class of cross-diffusion systems with entropy structure.  ArXiv preprint 2204.06080, 2022.
  * [6] M. Burger, M. Di Francesco, J.-F. Pietschmann, and B. Schlake.  Nonlinear cross-diffusion with size exclusion.  SIAM J. Math. Anal., 42(6):2842–2871, 2010.
  * [7] M. Burger, B. Schlake, and M.-T. Wolfram.  Nonlinear Poisson-Nernst-Planck equations for ion flux through confined geometries.  Nonlinearity, 25(4):961–990, 2012.
  * [8] L. Chen, E. S. Daus, A. Holzinger, and A. Jüngel.  Rigorous derivation of population cross-diffusion systems from moderately interacting particle systems.  Journal of Nonlinear Science, 31(6):1–38, 2021.
  * [9] X. Chen and A. Jüngel.  When do cross-diffusion systems have an entropy structure?  Journal of Differential Equations, 278:60–72, 2021.
  * [10] D. Daners and P. Koch Medina.  Abstract evolution equations, periodic problems and applications, volume 279 of Pitman Research Notes in Mathematics Series.  Longman Scientific & Technical, Harlow; copublished in the United States with John Wiley & Sons, Inc., New York, 1992.
  * [11] J. Fischer.  Weak-strong uniqueness of solutions to entropy-dissipating reaction-diffusion equations.  Nonlinear Anal., 159:181–207, 2017.
  * [12] K. Hopf.  Weak–strong uniqueness for energy-reaction-diffusion systems.  Math. Models Methods Appl. Sci., 32(5):1015–1069, 2022.
  * [13] X. Huo, A. Jüngel, and A. E. Tzavaras.  Weak-strong uniqueness for Maxwell–Stefan systems.  ArXiv preprint 2110.05331, 2021.
  * [14] A. Jüngel.  The boundedness-by-entropy method for cross-diffusion systems.  Nonlinearity, 28(6):1963–2001, 2015.
  * [15] A. Jüngel and I. V. Stelzer.  Existence analysis of Maxwell–Stefan systems for multicomponent mixtures.  SIAM J. Math. Anal., 45(4):2421–2440, 2013.
  * [16] M. Liero and A. Mielke.  Gradient structures and geodesic convexity for reaction–diffusion systems.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 371(2005):20120346, 2013.
  * [17] P.-L. Lions.  Mathematical topics in fluid mechanics. Vol. 1, volume 3 of Oxford Lecture Series in Mathematics and its Applications.  The Clarendon Press, Oxford University Press, New York, 1996.  Incompressible models, Oxford Science Publications.
  * [18] S. R. Noureen, J. P. Owen, R. L. Mort, and C. A. Yates.  Swapping in lattice-based cell migration models.  ArXiv preprint 2109.14093, 2021.
  * [19] K. J. Painter.  Continuous models for cell migration in tissues and applications to cell sorting via differential chemotaxis.  Bulletin of Mathematical Biology, 71(5):1117, 2009.
  * [20] P. Quittner and P. Souplet.  Superlinear parabolic problems.  Birkhäuser Advanced Texts: Basler Lehrbücher. Birkhäuser/Springer, Cham, 2019.  Blow-up, global existence and steady states.
  * [21] T. Roubíček.  Nonlinear partial differential equations with applications, volume 153 of International Series of Numerical Mathematics.  Birkhäuser/Springer Basel AG, Basel, second edition, 2013.
  * [22] M. J. Simpson, K. A. Landman, and B. D. Hughes.  Multi-species simple exclusion processes.  Physica A: Statistical Mechanics and its Applications, 388(4):399–406, 2009.
  * [23] A. Unterreiter, A. Arnold, P. Markowich, and G. Toscani.  On generalized Csiszár-Kullback inequalities.  Monatsh. Math., 131(3):235–253, 2000.
  * [24] 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\.