---
abstract: 'This paper presents new analytical results for a class of nonlinear parabolic systems of partial different equations with small cross-diffusion which describe the macroscopic dynamics of a variety of large systems of interacting particles. Under suitable assumptions, we prove existence of classical solutions and we show exponential convergence in time to the stationary state. Furthermore, we consider the special case of one mobile and one immobile species, for which the system reduces to a nonlinear equation of Fokker-Planck type. In this framework, we improve the convergence result obtained for the general system and we derive sharper $L^{\infty}$-bounds for the solutions in two spatial dimensions. We conclude by illustrating the behaviour of solutions with numerical experiments in one and two spatial dimensions.'
address:
- 'Gran Sasso Science Institute, Viale Francesco Crispi 7, L’Aquila, 67100, Italy'
- 'University of Vienna, Oskar-Morgenstern-Platz 1, 1090 Vienna, Austria'
- 'Laboratoire Jacques-Louis Lions, Sorbonne-Université, 4 place Jussieu, 75252 Paris'
- 'Mathematics Institute, University of Warwick, Gibbet Hill Road, CV47AL Coventry, UK'
- ' Radon Institute for Computational and Applied Mathematics, Altenbergerstr. 69, 4040 Linz, Austria.'
author:
- Luca Alasio$^1$
- Helene Ranetbauer$^2$
- Markus Schmidtchen$^3$
- 'Marie-Therese Wolfram$^4$'
title: 'Trend to equilibrium for systems with small cross-diffusion'
---

Introduction
============

Background and motivation
-------------------------

In this paper we focus on a class of nonlinear cross-diffusion systems with sufficiently small off-diagonal diffusion terms. The smallness assumption ensures that the system we consider is [“close”]{}, in a suitable sense, to a linear, decoupled system. This allows us to adapt techniques developed in the theory of linear, parabolic systems of partial differential equations (PDEs) in order to study long-time behaviour of solutions of the nonlinear system. More specifically, we consider a parabolic system of PDEs of the form

$$\label{main0}
    {\dfrac{\partialu}{\partial t}} 
    -\nabla \cdot \left\{
    D(x) [ (I + \delta \Phi(u)) \nabla u 
    + (\diag(\nabla V) + \delta \Psi(u) ) u ]
    \right\}
    =
    0,$$

which is a compact notation for the system $${\displaystyle}{\dfrac{\partialu_{i}}{\partial t}} - {\displaystyle}\sum_{\alpha, \beta} {\dfrac{\partial}{\partial x_\alpha}} 
    \left\{ D_{i}^{\alpha\beta}(x) 
    \left[ \left(
    {\dfrac{\partialu_{i}}{\partial x_\beta}} + \delta \sum_j \Phi_{ij}(u){\dfrac{\partialu_{j}}{\partial x_\beta}}
    \right)
    +
    \left({\dfrac{\partialV_{i}}{\partial x_\beta}} u_i 
    + \delta \sum_j \Psi_{ij}(u)u_{j} \right)
    \right] \right\}  =  0,$$ for $1\leq i, j \leq m,$ and $1\leq \alpha, \beta \leq d$. Here $m\geq 1$ represents the number of components (or species) and $d\in\{1,2,3\}$ the number of space dimensions. We denote by $u_i$ the $i$-th species. We will specify the assumptions on the diffusion tensor $D$, nonlinear mobilities $\Phi$ and $\Psi$ and the potential $V$ in Section \[sec:model\].

Cross-diffusion systems arise in multiple contexts in Physics, Life Sciences and Social Sciences; in particular, they have been derived as formal macroscopic limits of several microscopic models describing multi-species systems in the presence of finite volume effects, size exclusion or joint population pressures (see, for example, [@Bruna:2012cg; @Burger:2010gb; @perthame2015parabolic; @simpson2009multi]). Models with finite volume effects, which ensure physical bounds on the density, received significant attention in the past years. Recently Gavish et al. derived a mean-field model for a one-dimensional hard rod system rigorously, see [@gavish2019large]. However these techniques can be used in 1D only; in higher space dimensions only approximate models have been developed so far. These mean field models share common features, such as degenerate diffusion and small cross-diffusion terms (see, e.g., Examples 1.1-1.3 in [@alasio2018stability]). They can be derived from a lattice based microscopic description, see, for instance, [@bodnar2005] and a subsequent formal Taylor expansion of the associated master equation, as considered in [@burger2012nonlinear; @burger2016]. Another derivation was performed in [@Bruna:2012cg; @Bruna:2012wu], where the authors derived a cross-diffusion system from an underlying stochastic microscopic representation using the method of matched asymptotic expansions. We highlight that the two approaches yield different continuum models, however, they share the property of having small cross-diffusion terms. The smallness assumption is justified since the cross-diffusion terms are of the same order of magnitude as the microscopic particle size.

Gradient flow techniques
------------------------

In recent years, gradient flow methods have been successfully employed to study certain families of cross-diffusion systems, see, among others, [@jungel2016entropy; @difrancesco2018nonlinear; @desvillettes2015entropic]. Techniques such as the boundedness-by-entropy principle introduced in [@jungel2015boundedness], provide a mathematical framework to ensure existence and uniqueness of solutions to general nonlinear cross-diffusion systems that exhibit a gradient flow structure. In general, entropy methods have been proven to be a very useful tool to analyse the long time behavior of evolution equations. Especially the Bakry-Emery strategy (see [@bakry1985]) provides the necessary convex Sobolev inequalities to quantify the trend to equilibrium. Then entropy-entropy dissipation estimates allow to deduce exponential decay rates for general classes of linear and nonlinear scalar evolution equations, see for example [@arnold2001convex]. The cross-diffusion systems mentioned above lack a full gradient flow structure (in the Wasserstein sense) in certain parameter ranges, even though the underlying microscopic system possesses a natural one. This lack is caused by approximations made due to the finite volume effects. A first connection between the large deviations of stochastic particle systems and macroscopic Wasserstein gradient flows was established in [@adams2011large] without finite volume constraints. Questions related to structural features of cross-diffusion systems to be interpreted as macroscopic gradient flows (as in [@zamponi2017analysis]) or to have a strong solution (see [@Berendsen2019]) were investigated rather recently.\
Unfortunately, the results and techniques mentioned above cannot be applied at this stage to the systems of PDEs we want to study, nevertheless in this work we present an alternative strategy that relies on the smallness of the off-diagonal cross-diffusion terms.

Numerical methods
-----------------

The development of computational methods for non-linear cross-diffusion systems, especially structure preserving schemes, advanced significantly with the recent analytic progress. Structure preserving methods are designed in such a fashion that they preserve important physical and structural features such as positivity, conservation of mass or the dissipation of the associated entropy. Owing to the fact that Wasserstein gradient flows are posed in the set of probability measures, conservation of mass (or probability) of solutions is an important physical feature and finite volume discretisations are a natural framework to guarantee this. In addition, in recent years, several advances have been made in designing flux approximations that are in agreement with the energy dissipation. Bessemoulin-Chatard and Filbet, see [@bessemoulin2012finite], were among the first to present a finite volume method for nonlinear degenerate parabolic equations, which resolved the long-time behaviour correctly. Based on their scheme, different finite volume schemes have been proposed for systems, see for example [@CHS18; @carrillo2018fvconvergence]. Other numerical approximations are based on the underlying variational Wasserstein gradient flow structure. These so-called variational schemes are often restricted to one space dimension, as the computational complexity of computing the Wasserstein distance in higher space dimension is significant, see, for instance, [@carrillo2016diffeo]. Also, convergence results are, to the best of our knowledge, restricted to one spatial dimension; cf. [@matthes2014].

Summary of the main results
---------------------------

We consider a family of systems of PDEs with small cross-diffusion terms (namely system ) for which existence and uniqueness of solutions were investigated in [@alasio2018stability], as recalled in Proposition \[prop:ABC\]. We extend the analysis presented in [@alasio2018stability] by showing in Theorem \[thm:main1\] that (under suitable assumptions) the solutions are global-in-time and classical. Furthermore, we provide insights into the equilibration behaviour of such cross-diffusion systems in Theorem \[lem:energysys\]. A special instance of system , to which we shall return in Section \[sec:agf\], was analysed in [@bruna2017cross] using gradient flow techniques. This reduced PDE (namely problem ) can be interpreted as the diffusion of hard-core interacting particles through a domain with obstacles distributed according to the given immobile function. We recall existing results for the full cross-diffusion system and its stationary states in Proposition \[prop:stationary\] and how they can be improved for the reduced PDE at hand. In addition we present a more general equilibration result as well as sharper $L^{\infty}$-bounds for the reduced non-linear PDE in Theorem \[thm1\] and Theorem \[theor:maximum\] respectively. We conclude by illustrating our analytic results with numerical simulations. In doing so we construct a finite volume scheme in one and two spatial dimensions, investigate the exponential rate of convergence to equilibrium numerically and illustrate the dynamics of the system in physically relevant scenarios. The impact of the immobile species is clearly visible in that the immobile species diminishes the mobility of the mobile species; cf. Section \[sec:numerics\].

The rest of paper is organised as follows: we start by recalling well-posedness results for general systems with small cross-diffusion and prove existence of global classical solutions in Section \[sec:model\]. Section \[sec:expdec\] states the main exponential convergence result for the full system. In Section \[sec:agf\] we consider a specific example of an asymptotic gradient flow system, describing a mobile species diffusing in the presence of an immobile given species. In particular we are able to prove $L^\infty$-bounds for such scalar, nonlinear equations of Fokker-Planck type. Numerical simulations illustrating the analytic results are presented in Section \[sec:numerics\].

Analysis of systems with small cross-diffusion {#sec:model}
==============================================

We start by introducing the cross-diffusion system under consideration and by discussing analytic properties of solutions. We recall existing well-posedness results and establish existence of classic solutions under certain conditions on the initial datum.

A well-posedness result {#ssec:wellposed}
-----------------------

Throughout the paper $m\geq 1$ represents the number of components and $d\in\{1,2,3\}$ corresponds to the number of space dimensions. In addition, we shall denote the $i$-th species by $u_i$ and the related initial distribution by $u_i^0\geq 0$. We call $\nu$ the outward normal of $\Omega$.

For $1\leq i, j \leq m,$ and $1\leq \alpha, \beta \leq d$, we introduce the following notation: $$\begin{aligned}
    \label{eq:1011}
    \left\{
    \begin{aligned}
        {\mathcal{A}}_{ij}^{\alpha\beta}(x,u) 
        &=  D_{i}^{\alpha \beta}(x)
        ( I_{ij}
        +\delta \Phi_{ij}(u)
        ),
        \\[\eqnspace]
        {\mathcal{B}}_{ij}^{\alpha}(x,u) 
        &=  \sum_\beta
        D_{i}^{\alpha \beta}(x)
        \left( I_{ij} {\dfrac{\partialV_j(x)}{\partial x_\beta}}
        +\delta \Psi_{ij}(u)
        \right),
    \end{aligned}
    \right.\end{aligned}$$ where $\delta$ is a small parameter such that $$0 \leq \delta \leq \delta_0,$$ for a suitable constant $\delta_0$ depending on $u_0$, $D$, $V$, $\Phi$, $\Psi$, and $\Omega$ (see [@alasio2018stability] for further details on the smallness assumption). We consider the following initial-boundary value problem for system : $$\begin{aligned}
    \label{main-general}
    \left\{\;
    \begin{aligned}
        {\displaystyle}{\dfrac{\partialu_{i}}{\partial t}} - {\displaystyle}\sum_{\alpha, \beta, j} {\dfrac{\partial}{\partial x_\alpha}} \left[{\mathcal{A}}_{ij}^{\alpha\beta}(x,u) {\dfrac{\partialu_{j}}{\partial x_\beta}} - {\mathcal{B}}_{ij}^{\alpha}(x,u)u_{j}\right] & =  0 &\text{in }& \Omega,\\[\eqnspace]
        \sum_{\alpha, \beta, j} \nu_{\alpha} \left[ {\mathcal{A}}_{ij}^{\alpha\beta}(x,u){\dfrac{\partialu_{j}}{\partial x_\beta}} - {\mathcal{B}}_{ij}^{\alpha}(x,u)u_{j}\right] & =  0 &\text{on }& \partial\Omega,\\[\eqnspace]
        u_i(0,\cdot) & = u_i^{0} &\text{in }& \Omega.
    \end{aligned}
    \right.\end{aligned}$$

Note that, due to the no-flux boundary condition, the total [“mass”]{} (or number of particles in the microscopic case) of each species is conserved.

Global existence and regularity of solutions for system was proved in [@alasio2018stability]. In this framework, it is essential to assume that the system is close to a diagonal, decoupled, linear problem. Such [“reference problem”]{} is given by the weak formulation of $$\begin{aligned}
    \label{eq:ModelZero}
    \left\{
    \begin{aligned}
       {\displaystyle}{\dfrac{\partialu_{i}}{\partial t}} - \sum_{\alpha, \beta} {\dfrac{\partial}{\partial x_\alpha}}
        \left[
        D_{i}^{\alpha\beta}(x)
            \left({\dfrac{\partialu_{i}}{\partial x_\beta}} - {\dfrac{\partialV_i}{\partial x_\beta}} u_{i}\right)
        \right] 
        &= {\displaystyle}0 &\text{in }&\Omega,\\[\eqnspace]
         \sum_{\alpha, \beta} \nu_{\alpha} \cdot \left[
            D_{i}^{\alpha\beta}(x)
            \left({\dfrac{\partialu_{i}}{\partial x_\beta}} - {\dfrac{\partialV_i}{\partial x_\beta}} u_{i}\right)
        \right]
        &=  0 &\text{on }&\partial\Omega, \\[\eqnspace]
        u_i(0,\cdot) &=  u_i^{0} &\text{in }& \Omega,
    \end{aligned}
    \right.\end{aligned}$$ for $1\leq i \leq m,$ and $1\leq \alpha, \beta \leq d$.\
We now state our assumptions in detail (notice that they are strictly related to those in [@alasio2018stability]).

1.  We assume that $D_{i}^{\alpha\beta}\in C^1(\bar{\Omega})$ is symmetric in the space indices $\alpha$ and $\beta$. Furthermore, we suppose that there exist two constants $\Lambda \geq \lambda >0$ such that for every $x\in\Omega$, $1\leq i \leq m$, and $\xi\in\mathbb{R}^{d}$, it holds that $$\label{eq:D-Elliptic}
        \lambda\left|\xi\right|^{2}
        \leq \sum_{\alpha,\beta}
        D_{i}^{\alpha \beta}(x)\xi^{\alpha}\xi^{\beta}
        \leq
        \Lambda\left|\xi\right|^{2},$$ for any $i=1,\ldots,m$.

2.  The dependence on $u$ of the nonlinear terms is (at least) of class $C^2$: $$\Phi,\Psi\in C^{2}\left(\RR^{m}\right)^{m\times m},  \qquad \text{and} \qquad \Phi(0)=\Psi(0)=0.\label{eq:LipschitzPerturbationAssumption}$$

3.  We assume that, for $i=1,\dots,m$, the potentials $V_i:\RR^d\to\RR$ are (at least) of class $C^2$.

4.  The domain $\Omega$ is bounded, connected and sufficiently smooth (e.g. of class $C^2$). The initial datum $u^{0}$ in is assumed to be non-negative and to belong to $H^{2}(\Omega)$ satisfying the following compatibility condition on the boundary, $\partial\Omega$, for $i=1,\ldots,m$: $$\label{eq:compatcondu_0}
        \sum_{\alpha,j}
        \left[
        \sum_{\beta}
        {\mathcal{A}}_{ij}^{\alpha\beta}(x,u^0){\dfrac{\partialu^0_{j}}{\partial x_\beta}}
        - 
        {\mathcal{B}}_{ij}^{\alpha}(x,u^0)u^0_{j} \right] \nu_\alpha=0.$$

Before recalling the well-posedness result from [@alasio2018stability] we introduce the following notation.

\[def:weaksol\] We denote by $W(Q_{T})$ the Banach space of functions with two spatial weak derivatives taking values in $L^{2}(\Omega)$ continuously in time, and one time derivative in $L^2(0,T;H^1(\Omega))$, that is, $$W\left(Q_{T}\right)  = C\left(\left[0,T\right];H^{2}(\Omega)\right)
    \cap H^{1}
    \left( 0,T; H^1(\Omega)\right).$$ Note that, given $T>0$, we denote the parabolic cylinder by $Q_{T}=(0,T)\times\Omega$. Also, we write $H^2(\Omega)$ for $H^2(\Omega; \mathbb R^m)$, and similarly for other spaces.

\[prop:ABC\] Let hypotheses **H1**-**H4** hold, then system admits a unique, global solution $u\in W(Q_{T})$ for all $T>0$. Furthermore, there exist constants $\Gamma_0, \Gamma_1 > 0$ independent of $T$ such that: $${\left\| u \right\|}_{C^0(Q_{T})}
    \leq
    \Gamma_0 {\left\| u \right\|}_{W(Q_{T})}
    \leq
    \Gamma_1 {\left\| u_0 \right\|}_{H^2(\Omega)}.$$

Existence of global classical solutions {#sec:reg}
---------------------------------------

We are now going to show that solutions (in the sense of Definition \[def:weaksol\]) of problem are, indeed, classical solutions if the initial data is sufficiently smooth. The arguments in this section follow closely those presented in [@alasio2019global] for a different class of parabolic systems. We start by showing that classical solutions exist for a short time (see Proposition \[thm:AT\]) and then extend their maximal time interval of definition to the positive half-line using uniform estimates (see Proposition \[thm:amann\]). This yields the main result of this section:

\[thm:main1\] Let hypotheses **H1**-**H4** hold and assume that $u_0$ is (at least) of class $C^2(\bar{\Omega})$. Then problem admits a unique, bounded, global, classical solution in the space $ C^0([0,\infty);C^2(\bar{\Omega})) \cap C^1([0,\infty);C^0(\bar{\Omega}))$.

In order to prove Theorem \[thm:main1\], we will use the following two fundamental building blocks:

\[thm:AT\] Under the hypotheses of Theorem \[thm:main1\], there exists a time $\tau \in(0,T)$ such that, given a sufficiently smooth initial datum $u_0$ that is compatible with the boundary conditions, problem has a unique solution $u$ in the interval $[0,\tau]$ which satisfies $$u\in C^{1+\alpha_0}([0,\tau];L^2(\Omega))\cap C^{\alpha_0}([0,\tau];H^2(\Omega)),$$ where $\alpha_0\in (0,\frac{1}{2})$ depends on $u_0$, ${\mathcal{A}}$ and ${\mathcal{B}}$. Furthermore, for all $\alpha_1\in(0,\alpha_0)$, we have $${\dfrac{\partialu}{\partial t}}, \; \nabla^2 u \; \in C^{\alpha_1}([0,\tau];C^0(\bar{\Omega})).$$

We denote by $J(u^0) \subseteq [0,\infty)$ the maximal interval of definition of classical solutions of given an initial datum $u^0$.

\[thm:amann\] Let all assumptions of Proposition \[thm:main1\] hold and consider a solution $u$ of problem . If an exponent $\eps\in (0,1)$ (not depending on time) exists, such that for any $T>0$, $$u\in C^\eps(J(u_0)\cap[0,T];C^0(\bar{\Omega})),$$ then $u$ is a global solution, i.e., $J(u_0)=[0,\infty)$.

Proposition \[thm:amann\] is a rephrasing of Theorem 2 in [@amann1989dynamic] with $\theta=0$. Notice also that the [“convective term”]{}, indicated by $f$ in the notation of [@amann1989dynamic], is affine in the gradient in our case. Additionally, we do not use the notation $BUC^\eps$ to denote the space of [“bounded, uniformly $\eps$-Hölder continuous functions”]{} but simply write that the Hölder exponent does not depend on time.

The following interpolation results provides sufficient conditions for the criterion in Proposition \[thm:amann\] to be satisfied.

\[lem:fractional\] Let $f: Q_T \to\RR$ be a function in $X = L^2(0,T;H^2(\Omega)) \cap H^1(0,T;L^2(\Omega))$, then $$f\in H^r(0,T;H^s(\Omega)), \;\;\; \text{ for all } \; r,s\geq 0 
    \; \text{ such that } \; r + \frac{s}{2} \leq 1,$$ and, in turn, $$f\in  C^{0,\eta}([0,T];C^{0,\theta}(\bar{\Omega})), 
\; \text{ for all } \; 
\eta,\theta\geq 0 
\; \text{ such that } \; 2\eta + \theta \leq \frac{1}{2}.$$

Thanks to the higher order extensions for Sobolev functions, we can define $f$ on a larger cylindrical domain $R\subseteq \RR^{d+1}$ containing $Q_T$. Introducing a cut-off function, we further extend $f$ to the whole space ensuring sufficiently fast decay at infinity. Let us call $g$ such an extension. We observe that the norm of $g$ in $X'=L^2(\RR;H^2(\RR^d))\cap H^1(\RR;L^2(\RR^d))$ is controlled by the corresponding norms of $f$ on $Q_T$. Let $\langle \kappa \rangle = (1+|\kappa|^2)^{1/2}$. Denoting by $(\omega, \kappa)$ the conjugate variables of $(t,x)$ in Fourier space, we have that $$\langle \omega \rangle \hat{g} \in L^2(\RR^{d+1}), 
    \quad
    \text{as well as}
    \quad 
    \langle \kappa \rangle^2 \hat{g} \in L^2(\RR^{d+1}).$$ This means that $
    \left(
    \langle \omega \rangle + \langle \kappa \rangle^2
    \right) \hat{g} \in L^2(\RR^{d+1})
$ and we obtain $
\langle \omega \rangle^r \langle \kappa \rangle^s |\hat{g}|
\leq
\left(
\langle \omega \rangle + \langle \kappa \rangle^2
\right)^{r+\frac{s}{2}} |\hat{g}|.
$ Thus we obtain the desired fractional Sobolev regularity provided that $r+\frac{s}{2}\leq 1$. We are also using the following inequality, relating the norms of $g$ and $f$, $${\left\| g \right\|}_{X'} \leq 2{\left\| f \right\|}_{X}.$$ Then the Hölder regularity follows from standard embeddings for fractional Sobolev spaces (see, e.g., [@dinezza2012hitchhiker]). In particular, for $r,s>\frac{1}{2}$, we take $\eta = r - \frac{1}{2}$ and $\theta = s - \frac{1}{2}$.

Thanks to Proposition \[thm:AT\], we know that classical solutions exist for short times and that, as explained in [@acquistapace1988quasilinear], they may be extended to a maximal interval of existence denoted by $J(u_0)$ by standard methods. In order to show that such solutions exist for arbitrarily large time we are going to use the criterion provided by Proposition \[thm:amann\]. In particular, we need Hölder continuity of $u$ with respect to time, as well as a uniform $L^\infty$-bound in the space variable. Thanks to Proposition \[prop:ABC\], we know that $u\in W$ with a bound that is uniform in time. This implies that $u\in X$ and hence we apply Lemma \[lem:fractional\] in order to obtain obtain uniform Hölder estimates. Thus Proposition \[thm:amann\] allows us to conclude the proof.

Poincaré’s inequality
---------------------

In the next sections we will need the following simple modification of Poincaré inequality (see e.g. [@troianiello Lemma 1.36]). We recall that, if $\Omega$ is convex, we may take $C_P\leq \diam(\Omega) \pi^{-1}$ (see [@payne1960optimal]); otherwise a detailed treatment of this type of inequalities for more general $\Omega$ can be found in [@ziemer2012weakly Chapter 4].

\[lem:poincare\] Consider a bounded, connected, Lipschitz domain $\Omega$, a function $f\in H^1(\Omega)$ and a subset $S\subset\Omega$ with positive Lebesgue measure. Suppose that $\int_S f(x) \d x=0$, then there exists a constant $C_P>0$ (depending on $\Omega$ and $S$) such that $$\label{poincare}
    {\left\| f \right\|}_{L^2(\Omega)} \leq C_P{\left\| \nabla f \right\|}_{L^2(\Omega)}.$$ Furthermore, given a family of functions $(f_\ell)_{\ell\in\RR} \subset H^1(\Omega)$, such that $\sigma_0 = \mathrm{inf}_{\ell\in\RR} \{
{|\Omega|^{-1}}{|\mathrm{ess\, supp}(f_\ell)|}
\} \in (0,1)$, then inequality holds for any function $f_\ell$ and the constant $C_P$ depends only on $\Omega$ and $\sigma_0$.

We argue by contradiction. If the statement did not hold, then for each $n>0$ we could find a function $f_n\in H^1(\Omega)$ such that $\int_S f_n(x) \d x=0$ and $${\left\| f_n \right\|}_{L^2(\Omega)} > n{\left\| \nabla f_n \right\|}_{L^2(\Omega)}.$$ Without loss of generality, we can assume that ${\left\| f_n \right\|}_{L^2(\Omega)}=1$ (by homogeneity of the norm). This implies that ${\left\| \nabla f_n \right\|}_{L^2(\Omega)}<n^{-1}$. Therefore, as $n\to\infty$, we can extract a subsequence converging strongly in $L^2(\Omega)$ to a limit function $f_\infty\in H^1(\Omega)$ which satisfies the following properties $${\left\| \nabla f_\infty \right\|}_{L^2(\Omega)}= 0, \quad 
    \int_S f_\infty(x) \d x=0,
    \quad \text{and} \quad
    {\left\| f_\infty \right\|}_{L^2(\Omega)}=1.$$ The first two conditions above imply that $f_\infty = 0$, which contradicts the third. The proof of the second statement is very similar and we therefore omit it.

\[rem:P\] Under the same assumptions of Lemma \[lem:poincare\], we have $$\label{eq:estim1}
    {\left\| \nabla f \right\|}_{L^2(\Omega)}
    \leq \frac{1}{1-C_P{\left\| \nabla V \right\|}_\infty}
    {\left\| \nabla f + f\nabla V \right\|}_{L^2(\Omega)},$$ provided that $C_P{\left\| \nabla V \right\|}_\infty < 1$. More explicitly, we have $$\begin{aligned}
    {\left\| \nabla f \right\|}_{L^2(\Omega)} 
    &\leq {\left\| \nabla f + f\nabla V \right\|}_{L^2(\Omega)}
    \,+\, {\left\| f\nabla V \right\|}_{L^2(\Omega)}
    \\[\alnspace]
    &\leq {\left\| \nabla f + f\nabla V \right\|}_{L^2(\Omega)}
    \,+\, {\left\| \nabla V \right\|}_{\infty}{\left\| f \right\|}_{L^2(\Omega)}
    \\[\alnspace]
    &\leq {\left\| \nabla f + f\nabla V \right\|}_{L^2(\Omega)}
    \,+\, C_P{\left\| \nabla V \right\|}_{\infty}{\left\| \nabla f \right\|}_{L^2(\Omega)}.\end{aligned}$$ Notice that the condition $C_P{\left\| \nabla V \right\|}_\infty < 1$ gives a restriction on the size of $\Omega$ (since the Poincaré constant depends on $\diam(\Omega)$) or, alternatively, on the $L^\infty$-norm of the drift term $\nabla V$. We will use this relation in estimate .

Trend to equilibrium {#sec:expdec}
====================

In this section we investigate the equilibration behaviour of solutions to system . We show that they converge exponentially fast to the stationary state in the $L^2$-norm if assumptions **H1**-**H5** are satisfied.\
We choose a more compact notation in the following. Let $u = (u_1, \ldots u_m)$, then system can be written as $$\label{eq:mainsys}
    \left\{
    \begin{array}{rl}
    {\dfrac{\partialu}{\partial t}} 
    -\nabla \cdot \left\{
    D(x) [ (I + \delta \Phi(u)) \nabla u 
    + (\diag(\nabla V) + \delta \Psi(u) ) u ]
    \right\}\remspace
    &=
    0
    \quad \text{ in } \Omega,
    \\[\eqnspace]
    \nu \cdot \left\{  D(x) [ (I + \delta \Phi(u)) \nabla u 
    + (\diag(\nabla V) + \delta \Psi(u) ) u ] \right\} \remspace
    &= 0
    \quad \text{ on } \de\Omega,
    \\[\eqnspace]
    u(0,x)\remspace
    &=
    u^0(x),
    \end{array}
    \right.$$ and the corresponding stationary problem: $$\label{eq:mainsysstat}
    \left\{
    \begin{array}{rl}
    \nabla \cdot \left\{ 
    D(x) [ (I + \delta \Phi(\us)) \nabla \us 
    + (\diag(\nabla V) + \delta \Psi(\us) ) \us ]
    \right\}\remspace
    &=
    0
    \quad \text{ in } \Omega,
    \\[\eqnspace]
    \nu \cdot \left\{D(x) [ (I + \delta \Phi(\us)) \nabla \us 
    + (\diag(\nabla V) + \delta \Psi(\us) ) \us ] \right\} \remspace
    &= 0
    \quad \text{ on } \de\Omega,
    \end{array}
    \right.$$ where we used the compact notation $$D = \left\{ D_{i}^{\alpha \beta} \right\}^{1\leq\alpha,\beta\leq d}_{ 1\leq i\leq m}, \qquad\text{and}\qquad \diag(\nabla V) = \left\{ {\dfrac{\partialV_i}{\partial x_\beta}} \right\}^{1\leq\beta\leq d}_{1\leq i\leq m},$$ as well as $$\Phi = \{ \Phi_{ij} \}_{1\leq i,j\leq m}, \qquad\text{and}\qquad \Psi = \{ \Psi_{ij} \}_{1\leq i,j\leq m}.$$

Recall that $\Omega$ is a bounded domain of class $C^2$ in $\RR^d$ (with $d=1,2,3$), and $0\leq\delta\leq \delta_0\leq 1$, $u^0_i\in C^2(\Omega)$ and $V_i\in C^{2}(\bar{\Omega})$.

Before stating and proving the main theorem of this section, we introduce some notation and simple inequalities that we will use in the proof.

\[rem:L\] Given a function $F:\RR^m\to\RR^{m\times m}$ of class at least $C^k$ ($k\geq 0$) such that $F(0)=0$, we define the quantity $$L_k(F,R) := {\left\| F \right\|}_{C^k(\overline{B_{R}(0)})};$$ in particular, given two functions $w,z:\RR^d\to\RR^m$ such that $w,z\leq M$ for some $M>0$ and for a.e. $x\in\RR^d$, then we have (by Taylor expansion) $$|F(w)|\leq L_0(F,M), \quad \text{ and } \quad 
    |F(w)-F(z)|\leq L_1(F,M) |w-z|.$$

\[rem:V\] We shall use the following notation: $${\left\| V \right\|}_{\infty} = \max_{1\leq i \leq m} {\left\|  V_i \right\|}_{L^\infty(\Omega)},
    \quad
    \text{ and } 
    \quad
    {\left\| \nabla V \right\|}_{\infty} = \max_{1\leq i \leq m} {\left\| \nabla V_i \right\|}_{L^\infty(\Omega)}.$$

We now state the main result of this section:

\[lem:energysys\] Let assumptions **H1**-**H4** hold and suppose that $u\in W(Q_T)$ and $\us \in W^{1,\infty}(\Omega)$ are solutions of problem and respectively. Let us denote by $M>0$ a constant not depending on $T$ such that $$\label{eq:Mbound2}
    |u(t,x)|\leq M,$$ for any $(t,x)\in[0,T]\times\Omega$. Furthermore, suppose that $$\label{eq:deltabound2}
    \delta \leq 
    \min \left\{
    \delta_0,\,\,
    \frac{\lambda (1-C_P{\left\| \nabla V \right\|}_{\infty})}{2\Lambda e^{{\left\| V \right\|}_{\infty}}\Gamma_M(1+C_P{\left\| \nabla V \right\|}_{\infty})}
    \right\},
    \quad
    \text{ and }
    \quad
    C_P{\left\| \nabla V \right\|}_\infty < 1,$$ where $\Gamma_M$ is given by . Then the following inequality holds for any $t_1, t_2 \geq 0$: $$\label{eq:eneryest2}
    {\left\| u(t_2) - \us \right\|}_{L^2(\Omega)}
    \leq 
    \exp\left(3{\left\| V \right\|}_\infty- \frac{\lambda}{2C_P^2} (t_2 - t_1)\right)
    {\left\| u(t_1) - \us \right\|}_{L^2(\Omega)}.$$

Let $h := u - \us$, then $h$ satisfies the following system (with no-flux boundary conditions) in $\Omega$:

$$\label{eq:mainsysh}
    \begin{split}
    {\dfrac{\partialh}{\partial t}} 
    =\nabla \cdot \left\{
    D(x) \left[ (I + \delta \Phi(u)) \nabla h
    +\delta ( \Phi(u)-\Phi(\us) ) \nabla \us
    + \left(\diag(\nabla V) + \delta \Psi(u) \right) h
    + \delta (\Psi(u) - \Psi(\us))\us\right]
    \right\}.
    \end{split}$$

Testing system against $h \exp(V)$ gives $$\begin{aligned}
    \label{eq:testhV}
    {\frac{\mathrm{d}}{\mathrm{d}t}}\int_{\Omega} \frac{1}{2} \sum_{i} e^{V_i} h_i^2 \d x
    \,+\,
    \int_{\Omega} \sum_{ij} e^{V_i}
    \left[ \nabla h_i + h_i\nabla V_i \right]
    \cdot D_{ij} \left[ \nabla h_i + h_i\nabla V_i 
    + \delta \nn_j
    \right] \d x
    = 0,\end{aligned}$$ where $$\begin{aligned}
    \nn_j &= \nn_j^{\mathrm{I}} + \nn_j^{\mathrm{II}} + \nn_j^{\mathrm{III}} + \nn_j^{\mathrm{IV}}
    \\[\alnspace]
    &= \sum_k \left\{
    \Phi_{jk}(u)\nabla h_k
    + [ \Phi_{jk}(u)-\Phi_{jk}(\us) ] \nabla \us_k
    + \Psi_{jk}(u) h_k
    + [\Psi_{jk}(u) - \Psi_{jk}(\us)]\us_{k}
    \right\}.\end{aligned}$$ Next we estimate the order $\delta$ terms in and define $$\int_{\Omega} \sum_{ij} e^{V_i}
    \left[ \nabla h_i + h_i\nabla V_i \right]
    \cdot D_{ij} 
    \left[
    \nn_j^{\mathrm{I}} + \nn_j^{\mathrm{II}} + \nn_j^{\mathrm{III}} + \nn_j^{\mathrm{IV}}
    \right] 
    \d x
    =
    \jj^{\mathrm{I}} + \jj^{\mathrm{II}} + \jj^{\mathrm{III}} + \jj^{\mathrm{IV}}.$$ In particular, making extensive use of Remarks \[rem:L\], \[rem:V\], and \[rem:P\], we derive the following four estimates: $$\begin{aligned}
    \jj^{\mathrm{I}} 
    &= \int_{\Omega} 
    \sum_{ijk} e^{V_i}
    \left[ \nabla h_i + h_i\nabla V_i \right]
    \cdot D_{ij}
     \Phi_{jk}(u)\nabla h_k
    \d x
    \\[\alnspace]
    & \leq
    \Lambda \,L_0(\Phi,M) e^{{\left\| V \right\|}_{\infty}}\, 
    \left[ 
    {\left\| \nabla h \right\|}_{L^2(\Omega)}^2
    + {\left\| \nabla V \right\|}_{\infty}{\left\| \nabla h \right\|}_{L^2(\Omega)}{\left\| h \right\|}_{L^2(\Omega)}
    \right]
    \\[\alnspace]
    & \leq
    \Lambda \,L_0(\Phi,M) e^{{\left\| V \right\|}_{\infty}}\,(1+C_P{\left\| \nabla V \right\|}_{\infty}){\left\| \nabla h \right\|}_{L^2(\Omega)}^2,
    \\
    & \quad
    \\[\alnspace]
    \jj^{\mathrm{II}}
    &=\int_{\Omega} 
    \sum_{ijk} e^{V_i}
    \left[ \nabla h_i + h_i\nabla V_i \right]
    \cdot D_{ij}
     ( \Phi_{jk}(u)-\Phi_{jk}(\us) ) \nabla \us_k
    \d x
    \\[\alnspace]
    & \leq
    \Lambda \,L_1(\Phi,M) \,e^{{\left\| V \right\|}_{\infty}}\, {\left\| \nabla \us \right\|}_{L^\infty(\Omega)}\,
    \left[ 
    {\left\| \nabla h \right\|}_{L^2(\Omega)}{\left\| h \right\|}_{L^2(\Omega)}
    +{\left\| \nabla V \right\|}_{\infty}{\left\| h \right\|}_{L^2(\Omega)}^2
    \right]
    \\[\alnspace]
    & \leq
    \Lambda \,L_1(\Phi,M)\, e^{{\left\| V \right\|}_{\infty}}\,
    {\left\| \nabla \us \right\|}_{L^\infty(\Omega)}\,
    C_P \,(1+C_P{\left\| \nabla V \right\|}_{\infty}){\left\| \nabla h \right\|}_{L^2(\Omega)}^2,\end{aligned}$$ as well as $$\begin{aligned}
    \jj^{\mathrm{III}}
    &=\int_{\Omega} 
    \sum_{ijk} e^{V_i}
    \left[ \nabla h_i + h_i\nabla V_i \right]
    \cdot D_{ij}
     \Psi_{jk}(u) h_k
    \d x
    \\[\alnspace]
    & \leq
    \Lambda\, L_0(\Psi,M) \,e^{{\left\| V \right\|}_{\infty}} \,
    \left[ 
    {\left\| \nabla h \right\|}_{L^2(\Omega)}\,{\left\| h \right\|}_{L^2(\Omega)}
    \,+\,{\left\| \nabla V \right\|}_{\infty}\,{\left\| h \right\|}_{L^2(\Omega)}^2
    \right]
    \\
    & \leq
    \Lambda\, L_0(\Psi,M) \,e^{{\left\| V \right\|}_{\infty}}\,C_P\,(1+C_P{\left\| \nabla V \right\|}_{\infty}){\left\| \nabla h \right\|}_{L^2(\Omega)}^2,
    \\
    & \quad
    \\[\alnspace]
    \jj^{\mathrm{IV}}
    &=\int_{\Omega} 
    \sum_{ijk} e^{V_i}
    \left[ \nabla h_i + h_i\nabla V_i \right]
    \cdot D_{ij}
    (\Psi_{jk}(u) - \Psi_{jk}(\us))\us_{k}
    \d x
    \\[\alnspace]
    & \leq
    \Lambda \,L_1(\Psi,M)\, e^{{\left\| V \right\|}_{\infty}} \,
    {\left\| \us \right\|}_{L^\infty(\Omega)}\,
    \left[ 
    {\left\| \nabla h \right\|}_{L^2(\Omega)}{\left\| h \right\|}_{L^2(\Omega)}
    +{\left\| \nabla V \right\|}_{\infty}{\left\| h \right\|}_{L^2(\Omega)}^2
    \right]
    \\[\alnspace]
    & \leq
    \Lambda\, L_1(\Psi,M) \,e^{{\left\| V \right\|}_{\infty}}\,
    {\left\| \us \right\|}_{L^\infty(\Omega)}
    \,C_P\, (1+C_P{\left\| \nabla V \right\|}_{\infty}){\left\| \nabla h \right\|}_{L^2(\Omega)}^2.\end{aligned}$$ Since $\exp(V_i) \geq 1$, combining equations and yields $$\begin{aligned}
    \label{eq:hLtwoestimate}
    \begin{split}
        {\frac{\mathrm{d}}{\mathrm{d}t}}\int_{\Omega} \frac{1}{2} |h|^2 \d x
        \,+\,
        &\int_{\Omega} \lambda
        |\nabla h + h\,\diag(\nabla V) |^2
        \d x
        \\[\alnspace]
        &\leq
        \delta
        \Lambda e^{{\left\| V \right\|}_{\infty}}
        \Gamma_M \,(1+C_P{\left\| \nabla V \right\|}_{\infty})
        {\left\| \nabla h \right\|}_{L^2(\Omega)}^2
        \\[\alnspace]
        & \leq
         \delta
        \Lambda e^{{\left\| V \right\|}_{\infty}}
        \Gamma_M \,
        \frac{1+C_P{\left\| \nabla V \right\|}_{\infty}}{1-C_P{\left\| \nabla V \right\|}_{\infty}}
        \int_{\Omega}
        |\nabla h + h\,\diag(\nabla V) |^2
        \d x,
    \end{split}\end{aligned}$$ where we used Remark \[rem:P\] in the Appendix and the constant $\Gamma_M$ is given by $$\begin{aligned}
    \label{const:gamma}
    \Gamma_M := 
    \max \big\{ &
    L_0(\Phi,M),\;\;
    L_1(\Phi,M)
    {\left\| \nabla \us \right\|}_{L^\infty(\Omega)}C_P,\;\;
    L_0(\Psi,M) C_P,\;\;
    L_1(\Psi,M) 
    {\left\| \us \right\|}_{L^\infty(\Omega)}
    C_P
    \big\}.\end{aligned}$$ By assumption , we have the bound $$\delta
    \Lambda e^{{\left\| V \right\|}_{\infty}}
    \Gamma_M
    \frac{1+C_P{\left\| \nabla V \right\|}_{\infty}}{1-C_P{\left\| \nabla V \right\|}_{\infty}} \leq \frac{\lambda}{2}.$$ This bound in conjunction with estimate yields $${\frac{\mathrm{d}}{\mathrm{d}t}}\int_{\Omega}  |h|^2 \d x
    +
    \lambda \int_{\Omega} 
    |\nabla h + h\,\diag(\nabla V) |^2
    \d x
    \leq 0.$$ Letting $g=h \exp(V)$, we obtain $${\frac{\mathrm{d}}{\mathrm{d}t}}\int_{\Omega}  e^{-2V}|g|^2 \d x
    +
    \lambda \int_{\Omega} 
    e^{-2V}
    |\nabla g |^2
    \d x
    \leq 0.$$ Using Poincaré’s inequality for $g$ and integrating in time over the interval $(t_1,t_2)\subseteq [0,\infty)$ yields $$\int_{\Omega}  |g(t_2)|^2 \d x
    +
    \frac{\lambda}{C_P^2} \int_{t_1}^{t_2}\int_{\Omega} 
    | g(\tau) |^2
    \d x \d \tau
    \leq e^{4{\left\| V \right\|}_\infty} \int_{\Omega}  |g(t_1)|^2 \d x.$$ Thanks to Gronwall’s inequality we obtain the following continuous-dependence estimate $$\int_{\Omega}  |g(t_2)|^2 \d x
    \leq 
    \exp\left(4{\left\| V \right\|}_\infty- \frac{\lambda}{C_P^2} (t_2 - t_1)\right) \int_{\Omega}  |g(t_1)|^2 \d x.$$ Finally, let us express this continuous-dependence estimate in terms of the original variables. To this end we recall that $$u - u^* =  h = g \exp(-V).$$ Substituting this relation into the preceding inequality, we deduce that $$\begin{aligned}
    {\left\| u(t_2) - \us \right\|}_{L^2(\Omega)}
    &\leq
    {\left\| h(t_2) e^V \right\|}_{L^2(\Omega)}
    \\[\alnspace]
    &\leq \exp\left(2{\left\| V \right\|}_\infty- \frac{\lambda}{2C_P^2} (t_2 - t_1)\right) 
    {\left\| h(t_1) e^V \right\|}_{L^2(\Omega)}
    \\[\alnspace]
    &\leq 
    \exp\left(3{\left\| V \right\|}_\infty- \frac{\lambda}{2C_P^2} (t_2 - t_1)\right)
    {\left\| u(t_1) - \us \right\|}_{L^2(\Omega)},\end{aligned}$$ and this concludes the proof.

Study of an asymptotic gradient flow for one or two species {#sec:agf}
===========================================================

In this section we focus on an asymptotic gradient flow system describing the macroscopic dynamics of diffusing species in the presence of volume-exclusion constraints for the agents or particles. Such systems often have small cross-diffusion terms and therefore fall into the class of problems introduced in Section \[sec:model\]. We will discuss the equilibration behaviour of their solutions and establish sharper $L^{\infty}$-bounds (recall that, thanks to Proposition \[prop:ABC\], we already know that solutions are essentially bounded if $\delta < \delta_0$).

A cross-diffusion model for two interacting diffusing species {#sec:agf:twospecies}
-------------------------------------------------------------

We consider the asymptotic gradient flow system for two species that was introduced by Bruna and Chapman in [@Bruna:2012wu]. It describes the behaviour of the densities of two different types of hard spheres and it was derived from a stochastic system of hardcore-interacting Brownian particles using the method of matched asymptotic expansions. Few analytic results are available, and the little results that exist seem to rely on the assumption that the asymptotic gradient flow system is close to its stationary state (see [@bruna2017cross]) or that it is close to being a diagonal, decoupled linear problem; cf. [@alasio2019global].

The asymptotic gradient flow system presented in [@Bruna:2012wu] can be written in the general form of . In particular:

$$\label{eq:amatrix}
    {\mathcal{A}}(x,u) = \left(
    \begin{array}{cc}
        D_{1}&0\\[\alnspace]
        0 & D_{2}
    \end{array}
    \right)
     \left[I +  
    \left(
    \begin{array}{cc}
        \delta_{1,1} u_1 - \delta_{1,2} u_2 & \delta_{1,3} u_1\\[\alnspace]
        \delta_{2,3} u_2 & \delta_{2,1} u_2 - \delta_{2,2} u_1
    \end{array}
    \right)
    \right],$$

and $$\label{eq:drift}
    {\mathcal{B}}(x, u)=
    \left(
    \begin{array}{cc}
        D_1 & 0\\[\alnspace]
        0 & D_2
    \end{array}
    \right)
    \left(
    \begin{array}{cc}
        \nabla V_1 & (\delta_{2,2}\nabla  V_2-\delta_{1,2} \nabla V_1 )u_1\\[\alnspace]
        (\delta_{1,2}\nabla  V_1-\delta_{2,2} \nabla V_2 )u_2  &\nabla V_2
    \end{array}
    \right),$$ where $\delta_{i,j}$ for $i\in\{1,2\}, \, j\in \{1,2,3\}$ are asymptotically small parameters with $\max_{i,j} \delta_{i,j}=:\delta <\delta_0$. The coefficients $\delta_{i,j}$ depend on the size of the particles and on the corresponding diffusion coefficients $D_1, D_2\geq 0$, for their specific values we refer to [@Bruna:2012wu]. In [@bruna2017cross] it was noted that system with the nonlinearities in and has an *asymptotic gradient flow structure*. In particular, it is well known that certain cross-diffusion systems possess a formal gradient-flow structure, that is, they can be formulated as $$\label{grad_flow_general}
    {\dfrac{\partialu}{\partial t}} - \nabla\cdot \left ( \mm \nabla \frac{\delta E}{\delta u} \right) = 0,$$ where $\mm \in \mathbb R^{m\times m}$ is usually referred to as the *mobility matrix* and $\delta E/\delta u$ is the *variational derivative* of the entropy functional, $E$. In order to highlight the connection between gradient flows and asymptotic gradient flows, let us consider the following entropy functional

\[eq:grad\_flow1\] $$E[u] = \int_{\Omega} \bigg [ u_1\log u_1 + u_2\log u_2 + u_1 V_1 + u_2 V_2 
    +\frac{1}{2} \left(  \delta_{1,1} u_1^2 +2(d-1)(\delta_{1,2}+\delta_{2,2}) u_1 u_2 + \delta_{2,1} u_2^2\right)  \bigg ] \d x,\label{eq:entropy_general}$$ as well as the mobility matrix $$\mm(u)=\begin{pmatrix}D_1 u_1(1-\delta_{1,2} u_2) & D_1 \delta_{2,2} u_1 u_2\\[\alnspace]
    D_2 \delta_{1,2} u_1 u_2 & D_2 u_2(1- \delta_{2,2} u_1)
    \end{pmatrix}.\label{eq:mobility_general}$$

The cross-diffusion system with diffusion and drift matrices and , can be rewritten as $$\label{gradflow_generalasy}
    {\dfrac{\partialu}{\partial t}} = \nabla\cdot\left ( \mm \nabla \frac{\delta E }{\delta u} -G\right),$$

where $G$ is given by $$G= \begin{pmatrix}
 D_1u_1u_2\left(\left(\delta_{1,1}\delta_{1,2}-(d-1)(\delta_{1,2}+\delta_{2,2})\delta_{1,2}\right)\nabla u_1+((d-1)\delta_{1,2}(\delta_{1,2}+\delta_{2,2})-\delta_{2,1}\delta_{2,2} )\nabla u_2 \right) \\[\alnspace]
 D_2u_1u_2\left(\left(\delta_{2,1}\delta_{2,2}-(d-1)(\delta_{1,2}+\delta_{2,2})\delta_{2,2}\right)\nabla u_2+((d-1)\delta_{2,2}(\delta_{1,2}+\delta_{2,2})-\delta_{1,1}\delta_{1,2} )\nabla u_1 \right) 
\end{pmatrix}.$$ For a precise definition of all the parameters we refer to the derivation of the model in [@Bruna:2012wu].

The discrepancy between the system above and the gradient-flow induced by is of order $\delta^{2}$, and one can show that they actually coincide if and only if both species have the same diffusivities $D_1 = D_2 \geq 0$ as well as the same particle sizes, cf. [@bruna2017cross]. In particular, having the same size and same diffusivities implies that $\delta_{i,1}=(d-1)(\delta_{1,2}+\delta_{2,2})$ for $i=1,2$ as well as $\delta_{1,j}=\delta_{2,j}$ for $j=1,2,3$, cf. [@bruna2017cross] and therefore $G \equiv 0$. Moreover, one can show that this is the only possibility for $G$ to vanish. We refer to [@bruna2017cross] for more details as well as the existence proof of unique stationary solutions in both cases. More specifically, we recall the following result:

\[prop:stationary\] Let $\delta =  \max\{\delta_{i,j}\}$, assume $0<\delta<\delta_0$ and suppose that the potentials in satisfy $V_i\in H^3(\Omega)$. Then Theorem \[lem:energysys\] applies to systems and and they admit unique stationary states in $H^3(\Omega)$, which we denote by $u_\infty$ and $u_*$, respectively. Additionally, there exists a constant $C>0$ such that $${\left\| u_* - u_\infty \right\|}_{H^3(\Omega)} \leq C\delta^2.$$

As a consequence, we deduce that the stationary states $u_\infty$ and $u_*$ are in the space $W^{1,\infty}(\Omega)$; this fact will be useful in the next subsections.

The scalar problem with a “frozen” component
--------------------------------------------

In the following we focus on a special case of system , namely $$\begin{aligned}
    \label{eq:main-equation}
    \left\{
    \begin{array}{rl}
    {\dfrac{\partialr}{\partial t}} 
    -\nabla \cdot \left\{(1 + \delta_1 r - \delta_2 b) \nabla r + \delta_3 r \nabla b+r(1-\delta_2 b)\nabla V\right\}\remspace
    &=0
    \quad \text{ in } \Omega,
    \\[\eqnspace]
     \nu \cdot \left\{(1 + \delta_1 r - \delta_2 b) \nabla r + \delta_3 r \nabla b+r(1-\delta_2 b)\nabla V\right\}\remspace
    &= 0
    \quad \text{ on } \de\Omega,
    \\[\eqnspace]
    r(0,x)\remspace{}
    &=
    r_0(x),
    \end{array}
    \right.\end{aligned}$$ where $r=r(t,x)$ describes the density of the mobile species diffusing in the presence of a given immobile species $b(x)$ (the frozen species). Notice that we have simplified the notation by letting $\delta_j = \delta_{i,j}$ and $\delta = \max_j \delta_j$ for $i,j=1,2$. We chose to use a different notation to emphasise that is a scalar equation. It can be obtained from system by setting $r=u_1, b=u_2$ and setting the diffusion coefficient as well as the external potential of the species $b$ to zero. Note that this equation can be derived as a macroscopic limit of a stochastic system with two types of particles, one species diffusing in a domain with fixed obstacles of a certain size interacting via hard-core collisions. Hence, the particles diffuse in a [“perforated domain”]{} with obstructions distributed according to the density $b(x)$; cf. [@bruna2017asymptotic] for more details.

As equation is a special case of , namely with immobile blue particles, Proposition \[prop:stationary\] ensures that equation admits a stationary state $r^*\in H^3(\Omega)$ provided that that $\delta$ is sufficiently small, that $b, V \in H^{3}(\Omega)$ and that for some $M_0>0$ it holds that $$0\leq r_0(x)\leq M_0.$$ From the discussion in Section \[sec:agf:twospecies\] we know that equation does not exhibit a full gradient flow structure, but only an asymptotic one. Hence, taking the entropy functional $$\begin{aligned}
    \label{entropy_fp}
    E[r]&=\int_\Omega r\log r +\delta_1 \frac{r^2}{2}+\delta_3 r b +rV\, \d x,\end{aligned}$$ and mobility matrix $    \mathcal{M}(r)=r(1-\delta_2 b),$ we observe that equation can be cast in the form where $G(r,\delta)=rb (\delta_1 \delta_2\nabla r+\delta_2 \delta_3 \nabla b).$

As a consequence of Proposition \[prop:ABC\], we have that Problem \[eq:main-equation\] admits bounded, regular solutions in the space $W$, however, in the rest of this section we will need less regularity for the solutions and we will work in a weaker [“parabolic”]{} space, $Z$, that we introduce in Definition \[def:parabolic\].

\[def:parabolic\] Consider the space $Z = C^0([0,T]; L^2(\Omega))\cap L^2(0,T;H^1(\Omega))$. We use the so-called [“parabolic norm”]{} $${\left\| f \right\|}_{Z}^2
    = \|f\|_{L^\infty(0,T, L^2(\Omega))}^2 + \|\nabla f\|_{L^2(0,T; L^2(\Omega))}^2.$$

Similarly to Theorem \[lem:energysys\], we obtain a result on exponential convergence to equilibrium. Note that in the scalar case, we can avoid making the smallness assumption on the $\nabla V(x)$ (see ), as well as the use of the $L^\infty$-bound on $r$.

\[thm1\] Let $d\in\{1,2,3\}$ and assume that $b(x), V(x) \in H^3(\Omega)$ are such that $V_l\leq V(x)\leq V_u$, for some $V_l,V_u \in {\mathbb{R}}$ and for all $x\in \Omega$. Moreover, if $d=1,2$, we suppose that $r\in Z$ is a weak solution of problem such that $\|r\|_Z \leq L$, for some constant $L>0$; if $d=3$ we also suppose that $0\leq r\leq M$ for some constant $M>0$ and a.e. $(t,x)\in[0,T]\times\Omega$. Let the stationary state $r^* \in W^{1,\infty}(\Omega)$. Then the following inequality holds for any $t_1, t_2 \geq 0$: $$\begin{aligned}
    \|r(t_2, \cdot) - r^*\|_{L^2(\Omega)} \leq  \|r(t_1,\cdot) - r^*\|_{L^2(\Omega)}
    \exp\left(
    \frac{3}{2}(V_u - V_l) - \frac{t_2-t_1}{2C_P^2}\right),\end{aligned}$$ for $\delta$ sufficiently small, satisfying relation in the proof.

In order to transform the terms in front of the potential $V$ to order $\delta$, we perform the change of variables $w= r\exp(V)$ implying $\nabla r =(\nabla w-w\nabla V) \exp(-V)$ and $$\begin{aligned}
\label{equ_potential}
    e^{-V}{\dfrac{\partialw}{\partial t}}&=\nabla \cdot \left\{e^{-V}((1+\delta_1 e^{-V} w-\delta_2 b)\nabla w-\delta_1 e^{-V}w^2\nabla V+\delta_3 w\nabla b)\right\}.\end{aligned}$$ Now let us consider a perturbation, $h$, of the stationary state $w^*$ of (which exists by Proposition \[prop:stationary\]), namely $$\begin{aligned}
    w(x,t) = w^*(x) + h(x,t),\end{aligned}$$ where $\int_\Omega h(x,t) = 0$ since the mass is conserved. We observe that $h$ satisfies the following equation: $$\begin{aligned}
\label{equ_linstab}
    e^{-V}{\dfrac{\partialh}{\partial t}}&= \nabla \cdot \left\{e^{-V} \left[(1+\delta_1 e^{-V}(h+w^*)-\delta_2 b)\, \nabla h\;+\;(\delta_1 e^{-V}\nabla w^*-\delta_1 e^{-V} (2w^*+h)\nabla V +\delta_3  \nabla b) \,h\right]\right\}.\end{aligned}$$ We test equation against $h$ to obtain $$\begin{aligned}
    \label{equ4}
    \begin{aligned}
        \half {\frac{\mathrm{d}}{\mathrm{d}t}}\int_\Omega e^{-V} h^2 \d x 
        &=-\int_\Omega e^{-V}  (1+\delta_1 e^{-V}(h+w^*)-\delta_2 b)|\nabla h|^2\,\d x\\[\alnspace]
        & \quad-\int_\Omega e^{-V}(\delta_1 e^{-V}\nabla w^*-\delta_1 e^{-V} (2w^*+h)\nabla V +\delta_3  \nabla b) h\cdot \nabla h \,\d x.
    \end{aligned}\end{aligned}$$

Integration in time over $(t_1,t_2)\subset [0,T]$ yields

$$\begin{aligned}
    \label{equ5}
    \begin{aligned}
        \half\int_\Omega e^{-V} h^2(t_2) \d x 
        &=-\int_{t_1}^{t_2}\int_\Omega e^{-V}  (1+\delta_1 e^{-V}(h+w^*)-\delta_2 b)|\nabla h|^2\,\d x\, \d t\\[\alnspace]
        & \quad-\int_{t_1}^{t_2}\int_\Omega e^{-V}(\delta_1 e^{-V}\nabla w^*-\delta_1 e^{-V} (2w^*+h)\nabla V +\delta_3  \nabla b) h\cdot \nabla h \,\d x
        + \half\int_\Omega e^{-V} h^2(t_1) \d x.
    \end{aligned}\end{aligned}$$

Using the fact that $\nabla V, w^*, \nabla w^*, \nabla b \in L^\infty(\Omega),$ we deduce that $$\begin{aligned}
\label{eq:someintermediateestimate}
\begin{split}
-&\int_{t_1}^{t_2}\int_\Omega e^{-V}
\left(\delta_1 e^{-V}\nabla w^*-2\delta_1 e^{-V} w^*\nabla V +\delta_3  \nabla b\right) h\cdot \nabla h   \d x\, \d t
\\[\alnspace]
&\leq \delta\, e^{-V_l}
\left[
e^{-V_l}\|\nabla w^*\|_{L^\infty(\Omega)}+ 
2e^{-V_l}
\|w^*\|_{L^\infty(\Omega)}
\|\nabla V\|_{L^\infty(\Omega)}
+\|\nabla b\|_{L^\infty(\Omega)}
\right]
\int_{t_1}^{t_2} \|h\|_{L^2(\Omega)}
\|\nabla h\|_{L^2(\Omega)} \d t
\\[\alnspace]
&\leq \delta\, K_1
\int_{t_1}^{t_2}\|\nabla h\|_{L^2(\Omega)}^2 \d t,
\end{split}\end{aligned}$$ where $K_1 = C_P\left[
e^{-V_l}\|\nabla w^*\|_{L^\infty(\Omega)}+ 
2e^{-V_l}
\|w^*\|_{L^\infty(\Omega)}
\|\nabla V\|_{L^\infty(\Omega)}
+\|\nabla b\|_{L^\infty(\Omega)}
\right]$. Furthermore, we have that $$\begin{aligned}
\label{eq:someintermediateestimate2}
\begin{split}
\int_{t_1}^{t_2}\int_\Omega 
\delta_1 e^{-2V} h^2 \nabla V \cdot \nabla h   \d x\, \d t 
&\leq 
\delta e^{-2V_l}\|\nabla V\|_{L^\infty(\Omega)}
\int_{t_1}^{t_2}\int_\Omega 
h^2 |\nabla h|   \d x\, \d t
\\[\alnspace]
& \leq
\delta e^{-2V_l}\|\nabla V\|_{L^\infty(\Omega)}
\int_{t_1}^{t_2}{\left\| h \right\|}_{L^4(\Omega)}^2{\left\| \nabla h \right\|}_{L^2(\Omega)}\, \d t
\\[\alnspace]
& \leq
\delta e^{-2V_l}\|\nabla V\|_{L^\infty(\Omega)}
C_{GN}\int_{t_1}^{t_2}{\left\| h \right\|}_{L^2(\Omega)}{\left\| \nabla h \right\|}_{L^2(\Omega)}^2\, \d t
\\[\alnspace]
& \leq
\delta K_2 \int_{t_1}^{t_2}{\left\| \nabla h \right\|}^2_{\Omega} \d t,
\end{split}\end{aligned}$$ where, for $d=2$, we set $K_2 = e^{-2V_l}\|\nabla V\|_{L^\infty(\Omega)}
C_{GN} L$ and we used a special case of Gagliardo-Nirenberg interpolation inequality known as Ladyzhenskaya’s inequality to bound ${\left\| h \right\|}_{L^4(\Omega)}$, as well as $\|h\|_{Z}\leq L$; on the other hand, for $d=3$, we set $K_2= 2e^{-2V_l}\|\nabla V\|_{L^\infty(\Omega)}C_S M$ and we used Sobolev’s embedding in order to estimate ${\left\| h \right\|}_{L^4(\Omega)}$. We define $$K(\delta) =   \delta(K_1 + K_2),$$ and we substitute equations and into equation , which, in turn, becomes $$\begin{aligned}
\half\int_\Omega e^{-V} h^2(t_2) \d x 
+\int_{t_1}^{t_2}\int_\Omega 
e^{-V}  \left[1+\delta_1 e^{-V}(h+w^*)-\delta_2 b - K(\delta)\right]|\nabla h|^2 
\d x \d t
\leq
\half\int_\Omega e^{-V} h^2(t_1) \d x.\end{aligned}$$ Recalling that $h+w^* = w$ is non-negative, we choose $\delta$ sufficiently small so that the following relation is satisfied: $$\label{eq:kdelta}
 K(\delta) < \frac{1}{2}-\delta_2b,$$ which implies $$1+\delta_1 e^{-V}(h+w^*)-\delta_2 b - K(\delta) > \frac{1}{2}.$$ Thus we have obtained the following inequality: $$\int_\Omega e^{-V} h^2(t_2) \d x 
+\int_{t_1}^{t_2}\int_\Omega 
e^{-V}|\nabla h|^2 
\d x \d t
\leq
\int_\Omega e^{-V} h^2(t_1) \d x.$$ Using the Poincare’s inequality and the fact that $V_l<V<V_u$, we have $$\begin{aligned}
\int_\Omega h^2(t_2) \d x 
+C_P^{-2} \int_{t_1}^{t_2}\int_\Omega 
|h|^2 
\d x \d t
\leq
e^{V_u - V_l}\int_\Omega  h^2(t_1) \d x.\end{aligned}$$ Thanks to Gronwall’s inequality we obtain the following estimate: $$\begin{aligned}
    \|h(t_2, \cdot)\|_{L^2(\Omega)} \leq  \|h(t_1, \cdot)\|_{L^2(\Omega)} \exp\left(\frac{1}{2}(V_u-V_l)
    - \frac{t_2-t_1}{2C_P^2}\right).\end{aligned}$$ Finally, let us switch back to the original variables by recalling the fact that $h = w-w^*$. We obtain $$\begin{aligned}
    \|w(t_2, \cdot) - w^*\|_{L^2(\Omega)} \leq  \|w(t_1,\cdot) - w^*\|_{L^2(\Omega)} \exp\left(\frac{1}{2}(V_u-V_l)
- \frac{t_2-t_1}{2C_P^2}\right).\end{aligned}$$ Since $V$ is bounded from above and below by assumption, a simple change of variable $r = w e^{-V}$ implies the exponential convergence to equilibrium for the variable $r$ as claimed in the statement, which concludes the proof.

A maximum principle for the scalar problem ($d=2$)
--------------------------------------------------

In this section we will improve the $L^\infty$-bounds for solutions of problems of type , in dimension $d=2$ (the same result can be obtained in one dimension). We adapt the strategy presented in [@ladyzhenskaia1988linear], Chapter V[^1]. In the context of elliptic equations, this approach is referred to as De Giorgi’s method (see, e.g., [@han2011elliptic]).

\[theor:maximum\] Let $d=2$ and suppose that $\delta$ is sufficiently small and let $r\in Z$ be a solution of problem . Then we have that $$\label{eq:max-princ}
    r(t,x) \leq (1+C\delta^{\sigma})\,
    \exp\left({\left\| V \right\|}_{L^\infty(\Omega)} - V(x)\right)\,
    {\left\| r_0 \right\|}_{L^\infty(\Omega)},$$ for a.e. $t\in [0,\infty)$, $x\in\Omega$ and for some $\sigma\in(\frac{1}{2}, 1)$. The constant $C>0$ does not depend on time.

The proof of Theorem \[theor:maximum\] uses several technical results which we shall recall in the following. Notice that, if $x$ is close to a maximum point of $V$, estimate is [“almost sharp”]{}, in the sense that $r$ is estimated by $ {\left\| r_0 \right\|}_{L^\infty(\Omega)}$ with a multiplicative constant that is very close to $1$.

The following result is essential in order to prove Theorem \[theor:maximum\].

\[lem:maxiimum\] Let $d=2$ and suppose that $z\in Z$ is a weak solution of the problem

$$\begin{aligned}
    \label{eq:generalscalar}
    \left\{
    \begin{array}{rl}
    \omega(x) {\dfrac{\partialz}{\partial t}} -\nabla \cdot \left\{\omega(x) ( (1+\delta A(x,z))\nabla z + \delta F(x,z) ) \right\} \remspace &= 0  \quad \text{ in } \Omega,\\[\eqnspace]
    \nu \cdot \left\{(1+\delta A(x,z) \nabla z + \delta F(x,z))\right\} \remspace &= 0  \quad \text{ on } \de\Omega, \\[\eqnspace]
    z(0,x)\remspace  &= z_0(x).
    \end{array}
    \right.\end{aligned}$$

We assume that $\delta\in(0,1)$ is sufficiently small (see and the assumptions below) and, furthermore, we suppose that

(I) $\omega:\RR^d\to\RR$ satisfies $0 < \mu \leq \omega(x) < \mu^{-1}$, for a.e. $x\in\Omega$,\[en:UpperLowerOnOmega\]\

(II) $A:\RR^d\times\RR\to\RR$ and $\delta$ satisfy $0<\lambda \leq 1 + \delta A(x,z)$, for a.e. $x\in\Omega$,\

(III) $F(x,z) = z F_1(x) + z^2 F_2(x)$, where $F_i\in L^{\infty}(\Omega)$, and $i=1,2$,\

(IV) there exists $\bar{M}>0$ such that $0 \leq z_0(x)\leq \bar{M}$, for a.e. $x\in\Omega$,\

(V) there exists a constant $L>0$ depending only on $z_0, \Omega, F, \alpha, \omega, \delta$ such that ${\left\|  z  \right\|}_{Z} \leq L$.

Then there exists a constant $c>0$, not depending on time and given in , such that $${\left\| z \right\|}_{L^\infty(\Omega)}\leq c \bar{M}.$$

Before proving Lemma \[lem:maxiimum\], we state a preliminary result from Chapter II in [@ladyzhenskaia1988linear]. It can be proved directly by induction.

\[lem:sequence\] Suppose that a sequence $(a_n)_{n\in\NN} \subset {\mathbb{R}}_+$, of non-negative numbers satisfies the recursive relation $$a_{n+1} \leq \kappa \, \zeta^n \, a_n^{1+\eps}
\quad \text{ and } \quad
a_0 \leq \kappa^{-\frac{1}{\eps}} \zeta^{-\frac{1}{\eps^2}},$$ for $\kappa, \eps > 0$ and $\zeta>1$. Then we have that $$a_n\leq \kappa^{-\frac{1}{\eps}} \, \zeta^{-\frac{1}{\eps^2}-\frac{n}{\eps}}.$$

We will use the result above in order to conclude the following proof.

The proof consists of several steps. We begin with an energy estimate involving Stampacchia’s truncation and continue with an estimate for the measure of the superlevel sets $$S_{z>k}(t):=\{x\in \Omega \, : \, z(x,t) > k \},$$ for $k\in \RR_+$ and $t>0$. We will often use the abbreviation $\{z > k \}$ instead of $\{x\in \Omega \, : \, z(x,t) > k \}$. The strategy of the proof is the following: via careful use of a priori estimates, we will construct a sequence of the form $$a_n = \left(
    \int_0^T|S_{z>M(2-2^{-n})}(t)| \d t
    \right)^q,$$ for $n\in\NN$ and suitable $M>0$ as well as $q>0$. Subsequently, we will apply Lemma \[lem:sequence\] to $a_n$ in order to deduce that $|S_{z>M}(t)| = 0$ for a.e. $t\in[0,T]$. This implies that $z$ is essentially bounded.

**Step 0: Stampacchia’s truncation**\
Let $\tilde{M}\geq \bar{M}$ (to be chosen later) and $k\in\RR_+$, $k>\max(\tilde{M}, 1)$. It is easily verified that $(z-k)_+ \in Z$ (see, e.g. Theorem 2.1.11 in [@ziemer2012weakly] or Lemma 4.4 in Chapter II of [@ladyzhenskaia1988linear]). Thus we may test equation against the truncated function $(z-k)_+$ and obtain $${\frac{\mathrm{d}}{\mathrm{d}t}}\int_\Omega
    \frac{1}{2}\omega(x)  (z-k)_+^2 
    \d x
    =
    -\int_\Omega
     [ \omega(x) (1+\delta  A(x,z))|\nabla (z-k)_+|^2 + \delta \omega(x) F(x,z)\cdot\nabla (z-k)_+ ]
     \d x.$$ Rearranging this relation we get $$\label{eq:stamptrunc}
    {\frac{\mathrm{d}}{\mathrm{d}t}}\int_\Omega
    \frac{1}{2}\omega(x)  (z-k)_+^2 
    \d x
    +
    \int_\Omega
    \omega(x) \lambda |\nabla (z-k)_+|^2
       \d x
       \leq
     \frac{\delta}{\mu}  
     \left|  \int_\Omega   
     F(x,z)\cdot\nabla (z-k)_+ 
     \d x  \right|.$$

**Step 1: upper bound for the RHS in terms of $|\{z\geq k\}|$**\
Let us focus on the estimate of the right-hand side: $$\begin{aligned}
     \int_\Omega   
     F(x,z)\cdot\nabla (z-k)_+ 
     \d x
     &=
        \int_\Omega   
        ( z F_1(x) + z^2 F_2(x) )\cdot\nabla (z-k)_+
        ) \d x \\[\alnspace]
     & \leq
         {\left\| F_1 \right\|}_{L^\infty(\Omega)} 
         \int_\Omega
         ((z-k)_+ + k)|\nabla (z-k)_+| 
          \d x\\[\alnspace]
     & \qquad
        +   
        {\left\| F_2 \right\|}_{L^\infty(\Omega)}
        \int_\Omega
        ((z-k)_+ + k)^2 |\nabla (z-k)_+|
        \d x\\[\alnspace]
     & \leq
        C_{V}
       \int_\Omega
       \left\{
       (z-k)_+^2|\nabla (z-k)_+|
       + 3k(z-k)_+|\nabla (z-k)_+|
       \right\}
        \d x\\[\alnspace]
    &\qquad    
        + C_{V}
       \int_\Omega (k^2+k)|\nabla (z-k)_+|   \d x\\[\alnspace]
    & \leq
        \frac{5}{2} C_{V}
        \int_\Omega
        (z-k)_+^2|\nabla (z-k)_+|
        \d x
        + \frac{7}{2} C_{V}
       \int_\Omega k^2|\nabla (z-k)_+|   \d x,\end{aligned}$$ where we set $C_{V}= \max\{ {\left\|  F_1 \right\|}_{L^\infty(\Omega)}, {\left\|  F_2 \right\|}_{L^\infty(\Omega)} \}$ and used the fact that $k\geq 1$.\
Notice that the integral on the right-hand side of equation may be estimated further by $$\begin{aligned}
    \label{eq:13180706}
    \begin{split}
         \int_\Omega   
         F(x,z)\cdot\nabla (z-k)_+ 
         \d x
         & \leq
           \frac{C_{V}}{2}
            \int_\Omega
            \left\{
           5(z-k)_+^2|\nabla (z-k)_+|
           +
        7 k^2|\nabla (z-k)_+| 
        \right\} \d x
        \\[\alnspace]
         & \leq
           \frac{C_{V}}{2}
            \left\{
           5{\left\| (z-k)_+ \right\|}_{L^4(\Omega)}^2{\left\| \nabla (z-k)_+ \right\|}_{L^2(\Omega)}
           +
        7 k^2{\left\| \nabla (z-k)_+ \right\|}_{L^1(\Omega)} 
        \right\}.
    \end{split}\end{aligned}$$ Since $z$ is a Lebesgue-integrable function, it cannot be infinite on a region of full measure, therefore, for a sufficiently large $k$ (depending only on $\int_\Omega z_0 \d z$, thanks to Chebyshev’s inequality), the function $(z-k)_+$ vanishes on a region of positive measure. This implies that the modified version of Poincaré’s inequality given in Lemma \[lem:poincare\] holds, provided we can show that the Poincaré constant $C_P$ is uniformly bounded for all possible solutions $z$ (we do not assume uniqueness). In order to control the Poincaré constant it is essential to obtain a lower bound on the measure of the set $\{x \,\big| \, (z-k)_+ = 0\}$. To do so, let us make the following remark:

\[rem:cheb\] Let $f\in L^1(\Omega)$ be any non-negative function such that $\int_\Omega f \d x = m$, for some positive constant $m$ depending only on $\Omega$. Recall that, by Chebyshev’s inequality, for any $q>0$ it holds that $$|\{ x \,\big|\, f \geq q, \text{ a.e.}\}|
\leq \frac{1}{q} \int_\Omega f \d x.$$ Then we have, for $q=\frac{m}{\alpha |\Omega|}$ and an arbitrary $\alpha\in(0,1)$, $$\left| \left\{ x \,\big|\, f(x) < \frac{m}{\alpha |\Omega|} 
    \right\} \right|
    \,=\, |\Omega| - 
    \left| \left\{ x \,\big|\, f(x) \geq  \frac{m}{\alpha |\Omega|}
    \right\} \right|
    \,\geq\, |\Omega| - \frac{\alpha |\Omega|}{m} \int_\Omega f \d x
    \,=\, |\Omega|(1-\alpha).$$

We recall that integrating over the whole domain gives $\int_\Omega \omega z \d x = \int_\Omega \omega z_0 \d x$, for any $t>0$, as $z$ is continuous with respect to time. We can now argue as in Remark \[rem:cheb\] with $m=\int_\Omega \omega z_0 \d x$, $f=\omega z$, $\alpha \in(0,1)$ and $k \geq \frac{1}{\alpha \mu|\Omega|}\int_\Omega \omega z_0 \d x$, to the set $\{x \,\big|\, \omega z \leq \mu k\}$, which, by assumption (\[en:UpperLowerOnOmega\]), is contained in the set we are interested in, i.e., $\{x \,\big|\, z(t,x) \leq k \}$. Consequently, the Poincaré constant $C_P$ given in the second part of Lemma \[lem:poincare\] (with $\sigma_0 = 1-\alpha$) is always bounded independently of $z$.

Since $d=2$, we derive a bound for ${\left\| (z-k)_+ \right\|}_{L^4(\Omega)}^2$ applying Gagliardo-Nirenberg and Poincaré’s inequalities (see Lemma \[lem:poincare\]); namely we get: $$\begin{aligned}
\label{eq:Lineq}
     {\left\| (z-k)_+ \right\|}_{L^4(\Omega)}^2 
    & \leq 
    C_{GN}
    {\left\| (z-k)_+ \right\|}_{L^2(\Omega)} {\left\| (z-k)_+ \right\|}_{H^1(\Omega)}
    \\[\alnspace]
    & \leq 
    C_{GN}
    {\left\| (z-k)_+ \right\|}_{L^2(\Omega)} \left(
    {\left\| \nabla(z-k)_+ \right\|}_{L^2(\Omega)}
    +
    {\left\| (z-k)_+ \right\|}_{L^2(\Omega)}
    \right)
    \nonumber \\[\alnspace]
    & \leq
      C_{GN} (1+C_P)  {\left\| (z-k)_+ \right\|}_{L^2(\Omega)} {\left\| \nabla(z-k)_+ \right\|}_{L^2(\Omega)}, 
      \nonumber \end{aligned}$$ This may be further estimated upon observing that $$\int_\Omega |z|^2\d x \geq \int_{\{z>k\}} (z-k+k)^2\d x
    =
    \int_{\{z>k\}}\left( 
    (z-k)^2 +2k(z-k) + k^2
    \right)\d x
    \geq {\left\| (z-k)_+ \right\|}_{L^2(\Omega)}^2.$$ Using this in conjunction with assumption (5) in the statement we obtain ${\left\| (z-k)_+ \right\|}_{L^2(\Omega)}\leq {\left\| z \right\|}_{Z}  \leq L$.

Then the previous inequality leads to the following estimate: $$\begin{aligned}
     {\left\| (z-k)_+ \right\|}_{L^4(\Omega)}^2 \leq  C_{GN} (1+C_P) L  {\left\| \nabla(z-k)_+ \right\|}_{L^2(\Omega)}.\end{aligned}$$

Notice that for the remaining term in estimate , we have $$k^2 \int_\Omega|\nabla (z-k)_+| \d x
    \leq 
    \frac{k^2}{2}  \left(|\{z\geq k\}| + \int_\Omega|\nabla (z-k)_+|^2 \d x\right).$$ Hence estimate becomes $$\int_\Omega   
     F(x,z)\cdot\nabla (z-k)_+ 
     \d x
     \leq
       \frac{C_{V}}{2} 
        \left\{
       \left[ 5 C_{GN} (1+C_P)L +\frac72 k^2\right]
       {\left\| \nabla(z-k)_+ \right\|}_{L^2(\Omega)}^2
       +
    \frac72 k^2 |\{z\geq k\}|
    \right\}.$$

Altogether we have obtained the inequality: $${\frac{\mathrm{d}}{\mathrm{d}t}}\int_\Omega
    \frac{1}{2}\omega(x)  (z-k)_+^2 
    \d x
    +
    \int_\Omega
     (\omega(x)\lambda - \delta C_F )|\nabla (z-k)_+|^2
       \d x
       \leq
    \frac{7\delta C_V}{4\mu} k^2 |\{z\geq k\}|,$$ where $C_F =  \frac{C_V}{2\mu} \left[5L  C_{GN} (1+C_P) + \frac{7}{2}M^2 \right]$. Since $\omega \geq \mu$ and $k\geq \bar{M}$, after an integration in time, we obtain $$\frac{\mu}{2}\max_{t\in[0,T]} \int_\Omega
    (z-k)_+^2 
    \d x
    +
    \int_0^T \int_\Omega
    (\mu\lambda - \delta C_F )|\nabla (z-k)_+|^2
       \d x \d t
       \leq
    \frac{7\delta C_V}{4\mu} k^2 \int_0^T|\{z\geq k\}|\d t.$$ Supposing that $\delta$ is such that $$\label{eq:deltabound1}
    \mu \lambda - \delta C_F >0,$$ we have $$\begin{aligned}
    \label{eq:12061640}
    {\left\|  (z-k)_+  \right\|}_Z^2
       \leq
    \delta \Gamma  k^2 \int_0^T|\{z\geq k\}| \d t,\end{aligned}$$ where $\Gamma = 7C_V(4\mu^2\min\{\mu/2,\; \mu\lambda - \delta C_F   \} )^{-1}.$\
**Step 2: lower bound for the LHS in terms of $|\{z\geq k\}|$**\
Notice that, as above, there exists a constant $C_e>0$ such that the following inequality holds for any $1\leq p < \infty$ (it is a direct consequence of Sobolev’s and Poincaré’s inequality): $${\left\| f \right\|}_{L^2(0,T;L^p(\Omega))} \leq C_e {\left\|  f  \right\|}_Z.$$ We also know that, for any $f\in Z$, and thus in particular for $f=(z-k)_+$, there holds ${\left\| f \right\|}_{L^\infty(0,T;L^2(\Omega))} \leq  {\left\|  f  \right\|}_Z$. Thanks to the Riesz–Thorin interpolation theorem, we deduce that, for $\nu=p^{-1} \in (0,1)$, $$\begin{aligned}
    \label{eq:finL4nu}
    {\left\| f \right\|}_{L^{4-\nu}(Q_T)} < \max\{1,C_e\}{\left\|  f  \right\|}_Z,\end{aligned}$$ where $C_e$ does not depend on time. By Chebyshev’s inequality, for any $a\geq 1$, we obtain $$|\{ (z-h)_+ > 0 \}| = |\{ z > h \}| = |\{ z-k > h-k \}| 
    \leq 
    \frac{1}{(h-k)^a} 
    \int_{\{ z > h \}}  (z-k)_+^a \d x,$$ whence, upon integration, we obtain $$\int_0^T|\{ z > h \}| \d t
    \leq
    \frac{1}{(h-k)^a} 
    \int_0^T\int_\Omega  (z-k)_+^a \d x \d t.$$ Choosing $a=4-\nu$, we get $$\begin{aligned}
    \label{eq:1206191631}
    \left(
    \int_0^T|\{ z > h \}| \d t
    \right)^{\frac{2}{4-\nu}}
    \leq
    \frac{1}{(h-k)^2} 
    \left(
    \int_0^T\int_\Omega  (z-k)_+^{4-\nu} \d x \d t\right)^{\frac{2}{4-\nu}}.\end{aligned}$$ **Step 3: combining all previous estimates**\
Combining , , and we obtain $$\begin{aligned}
    \label{eq:12061615}
    \begin{split}
    \left(
\int_0^T|\{ z > h \}| \d t
\right)^{\frac{2}{4-\nu}}
    &\leq
    \frac{1}{(h-k)^2} \max\{1,C_e^2\}
    {\left\|  (z-k)_+  \right\|}_Z^2 \\
    &\leq  \frac{ k^2}{(h-k)^2} {\tilde{\Gamma}_\delta}  \int_0^T|\{z\geq k\}| \d t,
    \end{split}\end{aligned}$$ where ${\tilde{\Gamma}_\delta} =  \max\{1,C_e^2\}   \delta \Gamma$. For $M > \bar{M}$, consider the increasing sequence defined by $$k_n = M(2-2^{-n}),$$ and $$a_n = \left(
    \int_0^T|\{ z > k_n \}| \d t
    \right)^{\frac{2}{4-\nu}}.$$ for any $n\in\NN$. In order to apply Lemma \[lem:sequence\] we need to define a recursive relation for the sequence $(a_n)_{n\in \NN}$. Let us set $$k= k_{n}, \quad h = k_{n+1},$$ then we observe that $$\begin{aligned}
    \frac{k^2}{(h - k)^2} = \frac{k_{n}^2}{(k_{n+1} - k_{n})^2} = \left(\frac{2 -2^{-n}}{2^{-n} - 2^{-n-1}}\right)^2
    =\left( 2^{(n+1)}(2-2^{-n}) \right)^2 \leq 4^{n+2},\end{aligned}$$ we observe that, from , $$\begin{aligned}
    a_{n+1} = \left(
    \int_0^T|\{ z > k_{n+1} \}| \d t
    \right)^{\frac{2}{4-\nu}} \leq 4^{n+2} \tilde{\Gamma}_\delta \int_0^T |\{z\geq k_{n}\}|\d t =  4^{n+2} \tilde{\Gamma}_\delta a_n^{\frac{4-\nu}{2}}.\end{aligned}$$ Thus the sequence $(a_n)_{n\in \NN}$ satisfies the following relation

$$\begin{aligned}
    a_{n+1} \leq   4^{n+2} \tilde{\Gamma}_\delta a_n^{\frac{4-\nu}{2}}.\end{aligned}$$

Let us set $\varepsilon = 1- \nu/2$, $\zeta=4$, and $\kappa=16 \tilde{\Gamma}_\delta$ and apply Lemma \[lem:sequence\]. We deduce that $a_n\to 0$ as $n\to\infty$, and that the following inequality holds: $$a_n\leq {(16\tilde{\Gamma}_\delta)} ^{-\frac{1}{\eps}} 4^{-\frac{1}{\eps^2}-\frac{n}{\eps}},$$ provided that $$a_0 = \left(
    \int_0^T|\{ z > \tilde{M} \}| \d t
    \right)^{\frac{2}{4-\nu}} \leq 
    {(16\tilde{\Gamma}_\delta)}^{-\frac{1}{\eps}} 4^{-\frac{1}{\eps^2}}.$$ This will be satisfied choosing $\tilde{M}$ sufficiently large. Thus, we shall now concentrate on giving an explicit lower bound for $\tilde{M}$. To this end, let $\theta >1$ be such that $\tilde{M} = \theta \bar{M}$. From and and upon choosing $h=\theta \bar{M}$ and $k= \bar{M}$, we also know that, $$\begin{aligned}
    \left(
    \int_0^T|\{ z > \theta \bar{M} \}| \d t
    \right)^{\frac{2}{4-\nu}}
    &\leq
    \frac{\max\{1,C_e^2\}}{\bar{M}^2 (\theta - 1)^2} 
    \| (z-\bar{M})_+ \|_{Z}^2.\end{aligned}$$ Notice that the quantity ${\left\|  (z-\bar{M})_+  \right\|}_Z$ is bounded uniformly in time by $ C_0 = L + \bar{M}$.

Therefore we impose that $\theta$ satisfies $$\frac{ C_0^2  }{\bar{M}^2(\theta-1)^2} \leq{(16\tilde{\Gamma}_\delta)}^{-\frac{1}{\eps}} 4^{-\frac{1}{\eps^2}},$$ which, in turn, gives $$\theta \geq 1 + \frac{ \bar{M}  }{C_0} 
    \sqrt{ 
        {(16\tilde{\Gamma}_\delta)}^{-\frac{1}{\eps}} 4^{-\frac{1}{\eps^2}} 
    },$$ such that the upper bound for $a_0$ is met. We conclude that $$\int_0^T|\{ z > \theta \bar{M} \}| \d t = 0,$$ after passing $n\to \infty$, provided that $$\tilde{M}\geq \bar{M} \left(
    1 + \frac{ \bar{M}  }{C_0} 
    \sqrt{ 
        {(16\tilde{\Gamma}_\delta)}^{-\frac{1}{\eps}} 4^{-\frac{1}{\eps^2}} 
    },  \right),$$ and consequently we have obtained that, for a.e. $t,x$, $$\label{eq:estmax}
    z\leq \left( 1 + \frac{ \bar{M}  }{C_0}  {\tilde{\Gamma}_\delta} ^{\frac{1}{2\eps}} 4^{\frac{1}{\eps^2}}  \right) \bar{M} 
    = c\bar{M}.$$

\[rem:3d\] Lemma \[lem:maxiimum\] does not hold in general for $d > 2$. For example, estimate is not valid in the same form if $d=3$. The main difficulty consists in finding bounds for the term $\delta \nabla \cdot(\omega(x) z(t,x)^2 F_2(x))$ in equation . Sufficient conditions for $L^\infty$ estimates are presented in [@ladyzhenskaia1988linear], Chap. V, Sec. 2, p. 423, equation (2.3), and they impose restrictions on the growth rates of the flux terms. Such conditions could not be verified in our case.

With Lemma \[lem:maxiimum\] at hand, the proof of Theorem \[theor:maximum\] is immediate.

We apply Lemma \[lem:maxiimum\] making the following choices: $$\begin{aligned}
    & z=r \exp(V),
    && \omega = \exp(-V),
    \\
    & A(\cdot,r) = r-b,
    && L=\Gamma_0 \text{ (see Proposition \ref{prop:ABC})},
    \\
    & F_1 = \nabla(b + V),
    && F_2 = \nabla V.\end{aligned}$$ Recall that rewriting equation in terms of the new unknown $z = u \exp(V)$ we obtain $$e^{-V}{\dfrac{\partialz}{\partial t}}
    = \nabla \cdot \left\{e^{-V}((1+\delta e^{-V} z-\delta b)\nabla z-\delta (e^{-V}z^2\nabla V+z\nabla b))\right\},$$ with suitable initial and boundary conditions.

Numerical Simulations for the Scalar Problem {#sec:numerics}
============================================

Let us recall the equation $$\begin{aligned}
    \frac{\partial r}{\partial t} = \nabla \cdot ((1 + \delta_1 r - \delta_2 b) \nabla r + \delta_3 r \nabla b + r (1-\delta_2 b)\nabla V),\end{aligned}$$ with zero flux conditions at the boundary. We note that the equation can be cast into a transport part and a remainder part, i.e., $$\begin{aligned}
    \frac{\partial r}{\partial t}  = \nabla \cdot \left(r [\nabla \left(\log(r) + \delta_1 r + \delta_3 b\right) + (1-\delta_2 b)\nabla V]   - \delta_2 b \nabla r \right).\end{aligned}$$ This formulation is the basis of the finite volume scheme used for the following section. The scheme is based on the schemes studied in [@carrillo2018fvconvergence; @CHS18; @bessemoulin2012finite]. For simplicity we introduce it here in one space dimension but an extension to two (or more) dimensions is immediate. We note that we could have used a different formulation to define the numerical fluxes, for example, based on with entropy . However, the formulation chosen by us, is somewhat more natural, since it clearly highlights the transport part.\
In order to discretise the spatial domain, $\Omega = (-L,L)$, for some $L>0$, we introduce the computational mesh consisting of the control volumes $C_i=[x_{{i-1/2}}, x_{{i+1/2}})$, for all $i \in I:=\{1,\ldots, N\}$. The measure of each control volume is given by $|C_i| = \Delta x_i = x_{i+1/2} - x_{i-1/2}>0$, for all $i\in I$. Note that $x_{1/2}=-L$, and $x_{N+1/2}=L$. We also define $x_i = (x_{i+1/2} + x_{i-1/2})/2$ the centre of cell $C_i$ and set $\Delta x_{i+1/2}=x_{i+1}-x_i$ for $i=1,\ldots,N-1$.

Next, we discretise the initial data by computing the cell averages of the continuous initial data on each cell, i.e., $$\begin{aligned}
    r_i^0 := \frac{1}{\Delta x_i}\int_{C_i} r_{\mathrm{init}}(x)\,\d x,\end{aligned}$$ for all $i\in I$. Upon integrating the equation over $[t^n, t^{n+1}) \times C_i$, we obtain the following finite volume approximation \[eq:scheme\] $$\begin{aligned}
    \label{eq:scheme_evol}
    \frac{r_i^{n+1} - r_i^n}{\Delta t} = - \frac{{\mathbf{F}}_{{{i+1/2}}}^n - {\mathbf{F}}_{i-1/2}^n}{\Delta x_i},\end{aligned}$$ for $i\in I$. Here the numerical fluxes are given by $$\begin{aligned}
    \label{eq:numerical_fluxes}
    {\mathbf{F}}_{{{i+1/2}}}^n &= r_i^n \left[( -\d \xi)_{{{i+1/2}},+} + (1 -\delta_2 \,b_{{{i+1/2}}}\,) \left(-\frac{V_{i+1} - V_i}{\Delta x_{{i+1/2}}}\right)_+\right]\\
    &\quad + r_{i+1}^{n} \left[(-\d \xi)_{{{i+1/2}},-} + (1 -\delta_2 \,b_{{{i+1/2}}}\,) \left(-\frac{V_{i+1} - V_i}{\Delta x_{{i+1/2}}}\right)_-\right]\\
    &\quad + \delta_2 \,b_{{{i+1/2}}}\, \frac{r_{i+1}^n - r_i^n}{\Delta x_{{{i+1/2}}}},\end{aligned}$$ for $i = 1,\ldots, N-1$, and $$\begin{aligned}
    b_{{i+1/2}}= b(x_{{i+1/2}}), \end{aligned}$$ as well as $$\begin{aligned}
    \xi_i^n := \log(r_i^n + \varepsilon) \,+\, \delta_1 \, r_i^n \,+\, \delta_3 \, b(x_i),\end{aligned}$$ where $0< \varepsilon = 10^{-7} \ll 1$ is a small constant that is commonly chosen to regularise the logarithm. In addition, we impose the numerical no-flux boundary conditions ${\mathbf{F}}_{1/2} = {\mathbf{F}}_{N+1/2} =0$. As usual, we use $(z)_\pm$ to denote the positive (resp. negative) part of $z$, i.e., $$\begin{aligned}
    (z)_+ := \mathrm{max}(z,0), \qquad \text{and} \qquad (z)_-:=\mathrm{min}(z,0).\end{aligned}$$

One Dimensional Explorations
----------------------------

Let us fix $\Omega = [-5,5]$ as computational domain and set $\delta_1 = 10^{-1}$ in this subsection. The immobile species is given by $b(x) = 1/\sqrt{2\sigma^2 \pi}  \exp(- x^2/(2\sigma^2))$, with $\sigma = 10^{-1}$. The active species is initially distributed according to $r_0(x) = 1/\sqrt{2 \sigma^2 \pi}  \exp(-(x + 3)^2 / (2\sigma^2))$, with $\sigma = \sqrt{1/2}$. The external potential is chosen as $V(x)= (7/2 - x)^2 \chi_{\{x < 7/2\}}(x) + (x - 7/2)^5 \chi_{\{x\geq 7/2\}}(x)$; cf. Figure \[fig:explanation\_choices\]. This choice forces the mobile species to penetrate the immobile species as it migrates towards the minimum of the external potential. In Figure \[fig:voleffects\_inhibition\] we present the evolution of the mobile species for three different choices of $\delta_2 = \delta_3$. It is apparent that an increase in the parameters $\delta_2= \delta_3$ leads to a more and more inhibited migration. This behaviour is perfectly physical. In fact, the parameters $\delta_i$ are directly linked to the radii of particles of the mobile (resp. immobile) species. Since it is harder for large particles to traverse other large particles the motion is increasingly slowed down. For the second numerical exploration we choose $\tilde V(x)=x^2/2$ and assume that the mobile species is initially distributed according to $\tilde r_0(x) = \frac12 r_0(x-3) + \frac12 r_0(x+3)$, with $r_0$ as above; cf. Figure \[fig:explanation\_choices\]. Figure \[fig:voleffects\_inhibition2\] shows the migration of the two initial bumps of the mobile species towards the minimum of the external potential. Since the immobile species inhabits this region we observe a competition for space leading to an indentation in the density of the mobile species, showcasing the strong effect of finite size exclusion effects. In Figure \[fig:expconvoneD\] we show the exponential convergence of the solution towards equilibrium.

![The $L^2$-distance (blue) between the evolution and the (numerically computed) minimiser $u^\star$ on a semi-log scale. A line of slope $-1.4$ (green) is superimposed to highlight the exponential convergence rate. The parameter choice corresponds to that of Figure \[fig:voleffects\_inhibition2\] (a).[]{data-label="fig:expconvoneD"}](sim_13_log_L2_decay){width="60.00000%"}

Two Dimensional Explorations
----------------------------

In this subsection we set $\Omega = (-1.75,1.75)^2$ and $\delta_2 = \delta_3 = 10^{-1}$. Note that this way both parameters are in the physical range in the sense of [@Bruna:2012wu].

### Immobile species as porous medium

In this example we assume the passive species is spread ‘heterogeneously’ according to the distribution $b(x,y) = N^{-1}(1 + \cos(5 x) + \cos(5y))$, where $N$ is chosen such that $\int_{\Omega} b(x,y)\, dxdy = 1$. Initially the active species is spread around the origin, i.e., $r_0(x,y)=\frac{1}{\sqrt{2 \pi \sigma^2}}\exp\left(- \frac{r^2}{2 \sigma^2}\right)$, with $\sigma^2 = 0.5$. Moreover, we also choose $\delta_1 = 10^{-1}$ and we assume the absence of any external forces. Thus the dynamics are dictated only by the internal dispersal of the mobile species and its interaction with the immobile species. As the time evolution continues the volume exclusion effects imposed on the active species by the passive species become visible. In the final state, complementary regions are occupied by the passive species and the active species, respectively; cf. Figure \[fig:sim\_4\_evolution\].

![Initially the active species (red, left) is spread around the origin while the passive species is spread ‘heterogeneously’ (blue, right). As the time evolution continues the volume exclusion effects imposed on the active species by the passive species become visible. In the final state, complementary regions are occupied by the passive species and the active species, respectively.[]{data-label="fig:sim_4_evolution"}](sim_4_r_init "fig:"){width="40.00000%"} ![Initially the active species (red, left) is spread around the origin while the passive species is spread ‘heterogeneously’ (blue, right). As the time evolution continues the volume exclusion effects imposed on the active species by the passive species become visible. In the final state, complementary regions are occupied by the passive species and the active species, respectively.[]{data-label="fig:sim_4_evolution"}](sim_4_b_immobile "fig:"){width="40.00000%"} ![Initially the active species (red, left) is spread around the origin while the passive species is spread ‘heterogeneously’ (blue, right). As the time evolution continues the volume exclusion effects imposed on the active species by the passive species become visible. In the final state, complementary regions are occupied by the passive species and the active species, respectively.[]{data-label="fig:sim_4_evolution"}](sim_4_r_final "fig:"){width="40.00000%"} ![Initially the active species (red, left) is spread around the origin while the passive species is spread ‘heterogeneously’ (blue, right). As the time evolution continues the volume exclusion effects imposed on the active species by the passive species become visible. In the final state, complementary regions are occupied by the passive species and the active species, respectively.[]{data-label="fig:sim_4_evolution"}](sim_4_b_immobile "fig:"){width="40.00000%"}

After some time the evolution slows down as the mobile species reaches a stationary state. In Figure \[fig:sim\_4\_asymptotics\] we show the numerical long-time asymptotics on a semi-log scale. It can be seen that the dynamic relaxes to the numerically computed stationary state at an exponential rate.

![The $L^2$-distance (blue) between the evolution and the (numerically computed) minimiser $u^\star$ on a semi-log scale. A line of slope $-5/2$ (green) is superimposed to highlight the exponential convergence rate.[]{data-label="fig:sim_4_asymptotics"}](sim_4_log_L2_decay){width="60.00000%"}

### Immobile species as barrier

In this section we present another stunning example of how the immobile species can inhibit the evolution of the mobile species. We set $\delta_1 = 10^{-2}$ and consider the initial datum $$\begin{aligned}
    r_0(x,y)=\frac{1}{\sqrt{2 \pi \sigma^2}}\exp\left(- \frac{(x+3/4)^2 + (y+3/4)^2}{2 \sigma^2}\right),\end{aligned}$$ with $\sigma = 0.3$. The immobile species is fixed and given by $$\begin{aligned}
    b(x,y) = &N^{-1}\left(\sin\left(5 \rho(x,y) \right)\right) \chi_{\{-7/4 +7\pi /5 \leq \rho(x,y) \leq -7/4 + 7.5\pi /5\}}\\[0.3em]
            &\times 
                (1 + \cos(4\pi/3(y+1)))\chi_{\{-1\leq y \leq 1\}}(x,y)\\[0.3em]
            &\times 
                (1 + \cos(4\pi/3(x+1)))\chi_{\{-1\leq x \leq 1\}}(x,y),\end{aligned}$$ where $\rho(x,y) = \sqrt{(x+7/4)^2 + (y+7/4)^2}$ and $N$ normalises the mass to one; cf. Figure \[fig:sim\_8\_b\_acts\_as\_barrier\].

![Initial data of the active species (red, left) and the immobile species (blue, right) acting as a barrier.[]{data-label="fig:sim_8_b_acts_as_barrier"}](sim_7_u_init_3D "fig:"){width="45.00000%"} ![Initial data of the active species (red, left) and the immobile species (blue, right) acting as a barrier.[]{data-label="fig:sim_8_b_acts_as_barrier"}](sim_7_b_init_3D "fig:"){width="45.00000%"}

Below we consider two different choices of external potentials: $V_\mathrm{s}(x) =  \frac1{20} \left((x-7/4)^2 + (y-7/4)^2\right)$ or $V_\mathrm{w}(x) = 1/10 \, V_\mathrm{s}(x)$. In the subsequent simulations the immobile species acts as a barrier as the active species tries to reach the minimum of the external potential centred at the upper right corner of the domain.

In Figure \[fig:barrier\_strong\_potential\] it can be observed that the motion of the active species is slowed down and that its density is lower in the region occupied by the passive species compared to the unoccupied regions. Since the potential is relatively strong the active species moves through the immobile one only slightly changing its shape. In Figure \[fig:barrier\_weak\_potential\] we rescaled the potential by a factor of ten. Again, we see that the motion of the mobile species is slowed down and that its density is reduced in the region occupied by the passive species in comparison with unoccupied regions. Since, now, the potential is relatively weak the motion of the active species is inhibited by the passive species which incentivises circumnavigating the immobile species instead of penetrating it. This behaviour is reflected in the numerical simulation; cf. Figure \[fig:barrier\_weak\_potential\]. Compared to Figure \[fig:barrier\_strong\_potential\] the active species closes in from the side as it approaches the minimum of $V$, rather than directly moving into it.

**Acknowledgements.** The basic ideas of this paper were developed at the Oxford working group meeting [“Asymptotic gradient flows”]{} in 2017. The authors are very grateful to Maria Bruna, José Antonio Carrillo, Jon Chapman and Ende Süli for the inspiring and valuable discussions during this meeting. M.T. Wolfram acknowledges partial support by the Austrian Academy of Sciences via the New Frontier’s Group NST-001 and by the EPSRC vis the First grant EP/P01240X/1. H. Ranetbauer was funded by the Austrian Science Fund (FWF) project F 65.

[10]{}

P. Acquistapace and B. Terreni, *On quasilinear parabolic systems*, Mathematische Annalen **282** (1988), no. 2, 315–335.

S. Adams, N. Dirr, M. A. Peletier, and J. Zimmer, *From a large-deviations principle to the [W]{}asserstein gradient flow: a new micro-macro passage*, Communications in Mathematical Physics **307** (2011), no. 3, 791.

L. Alasio, M. Bruna, and Y. Capdeboscq, *Stability estimates for systems with small cross-diffusion*, ESAIM: Mathematical Modelling and Numerical Analysis **52** (2018), no. 3, 21109 – 1135.

L. Alasio and S. Marchesani, *Global existence for a class of viscous systems of conservation laws*, Nonlinear Differential Equations and Applications **26** (2019), no. 32.

H. Amann, *Dynamic theory of quasilinear parabolic systems*, Mathematische Zeitschrift **202** (1989), no. 2, 219–250.

A. Arnold, P. Markowich, G. Toscani, and A. Unterreiter, *On convex [S]{}obolev inequalities and the rate of convergence to equilibrium for [Fokker–Planck]{} type equations*, Communications in Partial Differential Equations **26** (2001), no. 1-2.

D. Bakry and M. [É]{}mery, *Diffusions hypercontractives*, S[é]{}minaire de Probabilit[é]{}s XIX 1983/84 (Berlin, Heidelberg) (Jacques Az[é]{}ma and Marc Yor, eds.), Springer Berlin Heidelberg, 1985, pp. 177–206.

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*, Journal of Evolution Equations.

M. Bessemoulin-Chatard and F. Filbet, *A finite volume scheme for nonlinear degenerate parabolic equations*, SIAM Journal on Scientific Computing **34** (2012), no. 5, B559–B583.

M. Bodnar and J. J. L. Velazquez, *Derivation of macroscopic equations for individual cell-based models: a formal approach*, Mathematical Methods in the Applied Sciences **28** (2005), no. 15, 1757–1779.

M. Bruna, M. Burger, H. Ranetbauer, and M.-T. Wolfram, *Asymptotic gradient flow structures of a nonlinear [F]{}okker-[P]{}lanck equation*, arXiv preprint arXiv:1708.07304 (2017).

, *Cross-diffusion systems with excluded-volume effects and asymptotic gradient flow structures*, Journal of Nonlinear Science **27** (2017), no. 2, 687–719.

M. Bruna and S. J. Chapman, *[Diffusion of multiple species with excluded-volume effects]{}*, The Journal of Chemical Physics **137** (2012), no. 20, 204116–204116–16.

, *[Excluded-volume effects in the diffusion of hard spheres]{}*, Phys. Rev. E **85** (2012), no. 1.

M. Burger, M. Di Francesco, J.-F. Pietschmann, and B. Schlake, *[Nonlinear [C]{}ross-[D]{}iffusion with [S]{}ize [E]{}xclusion]{}*, SIAM Journal on Mathematical Analysis **42** (2010), no. 6, 2842.

M. Burger, S. Hittmeir, H. Ranetbauer, and M.-T. Wolfram, *Lane formation by [S]{}ide-[S]{}tepping*, SIAM Journal on Mathematical Analysis **48** (2016), no. 2, 981–1005.

M. Burger, B. Schlake, and M.-T. Wolfram, *Nonlinear [P]{}oisson–[N]{}ernst–[P]{}lanck equations for ion flux through confined geometries*, Nonlinearity **25** (2012), no. 4, 961.

J. A. Carrillo, F. Filbet, and M. Schmidtchen, *Convergence of a finite volume scheme for a system of interacting species with cross-diffusion*, arXiv preprint arXiv:1804.04385 (2018).

J. A. Carrillo, Y. Huang, and M. Schmidtchen, *Zoology of a [N]{}onlocal [C]{}ross-[D]{}iffusion [M]{}odel for [T]{}wo [S]{}pecies*, SIAM Journal on Applied Mathematics **78** (2018), no. 2.

J. A Carrillo, H. Ranetbauer, and M.-T. Wolfram, *Numerical simulation of nonlinear continuity equations by evolving diffeomorphisms*, Journal of Computational Physics **327** (2016), 186–202.

L. Desvillettes, T. Lepoutre, A. Moussa, and A. Trescases, *On the entropic structure of reaction-cross diffusion systems*, Communications in Partial Differential Equations **40** (2015), no. 9, 1705–1747.

M. Di Francesco, A. Esposito, and S. Fagioli, *Nonlinear degenerate cross-diffusion systems with nonlocal interaction*, Nonlinear Analysis **169** (2018), 94–117.

E. Di Nezza, G. Palatucci, and E. Valdinoci, *Hitchhiker’s guide to the fractional [S]{}obolev spaces*, Bulletin des Sciences Math[é]{}matiques **136** (2012), no. 5, 521–573.

Nir Gavish, Pierre Nyquist, and Mark Peletier, *Large deviations and gradient flows for the brownian one-dimensional hard-rod system*, 2019.

Q. Han and F. Lin, *Elliptic partial differential equations*, vol. 1, American Mathematical Soc., 2011.

A. J[ü]{}ngel, *Entropy methods for diffusive partial differential equations*, Springer.

, *The boundedness-by-entropy method for cross-diffusion systems*, Nonlinearity **28** (2015), no. 6, 1963.

O. A. Ladyzhenskaia, V. A. Solonnikov, and N. N. Ural’tseva, *Linear and quasi-linear equations of parabolic type*, vol. 23, American Mathematical Soc., 1988.

D. Matthes and H. Osberger, *Convergence of a variational [L]{}agrangian scheme for a nonlinear drift diffusion equation*, ESAIM: Mathematical Modelling and Numerical Analysis **48** (2014), no. 3, 697–726.

L. E. Payne and H. F. Weinberger, *An optimal poincar[é]{} inequality for convex domains*, Archive for Rational Mechanics and Analysis **5** (1960), no. 1, 286–292.

B. Perthame, *Parabolic equations in biology*, Parabolic Equations in Biology, Springer, 2015, pp. 1–21.

M. J. Simpson, K. A. Landman, and B. D. Hughes, *Multi-species simple exclusion processes*, Physica A: Statistical Mechanics and its Applications **388** (2009), no. 4, 399–406.

G. M. Troianiello, *Elliptic differential equations and obstacle problems*, Springer Science & Business Media, 2013.

N. Zamponi and Ansgar J., *Analysis of degenerate cross-diffusion population models with volume filling*, **34** (2017), no. 1, 1–29.

W. P. Ziemer, *Weakly differentiable functions: Sobolev spaces and functions of bounded variation*, vol. 120, Springer Science & Business Media, 1989.

[^1]: Namely, Chapter V, Section 2, p. 425, Theorem 2.1; Chapter II, Section 6, p. 102, Theorem 6.1 (see also Lemma 5.6 on p. 95).
