Trend to Equilibrium for Systems with Small Cross-Diffusion

TL;DR

This paper analyzes nonlinear parabolic systems with small cross-diffusion, proving existence and exponential convergence to equilibrium, and provides sharper bounds and numerical illustrations for specific cases.

Contribution

It offers new analytical results on existence, convergence, and bounds for systems with small cross-diffusion, including specialized cases and numerical experiments.

Findings

Existence of classical solutions under certain conditions

Exponential convergence to stationary states

Sharper $L^{inity}$-bounds in two dimensions

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.

Equations

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

$$\frac{\partial u_i}{\partial t}-\sum _{\alpha ,\beta }\frac{\partial }{\partial x_{\alpha }}\left\{D_i^{\alpha \beta }(x)\left[\left(\frac{\partial u_i}{\partial x_{\beta }}+\delta \sum _j{\Phi }_{ij}(u)\frac{\partial u_j}{\partial x_{\beta }}\right)+\left(\frac{\partial V_i}{\partial x_{\beta }}{u}_i+\delta \sum _j{\Psi }_{ij}(u)u_j\right)\right]\right\}=0,$$

$$\{\begin{array}{rl}{\mathcal{A}}_{ij}^{\alpha \beta }(x,u)& =D_i^{\alpha \beta }(x)(I_{ij}+\delta {\Phi }_{ij}(u)),\\ {\mathcal{B}}_{ij}^{\alpha }(x,u)& =\sum _{\beta }{D}_i^{\alpha \beta }(x)\left(I_{ij}\frac{\partial V_j(x)}{\partial x_{\beta }}+\delta {\Psi }_{ij}(u)\right),\end{array}$$

$$0\le \delta \le {\delta }_0,$$

$$\{\begin{array}{rlrl}\frac{\partial u_i}{\partial t}-\sum _{\alpha ,\beta ,j}\frac{\partial }{\partial x_{\alpha }}\left[{\mathcal{A}}_{ij}^{\alpha \beta }(x,u)\frac{\partial u_j}{\partial x_{\beta }}-{\mathcal{B}}_{ij}^{\alpha }(x,u)u_j\right]& =0& \text{in }& \Omega ,\\ \sum _{\alpha ,\beta ,j}{\nu }_{\alpha }\left[{\mathcal{A}}_{ij}^{\alpha \beta }(x,u)\frac{\partial u_j}{\partial x_{\beta }}-{\mathcal{B}}_{ij}^{\alpha }(x,u)u_j\right]& =0& \text{on }& \partial \Omega ,\\ u_i(0,\cdot )& =u_i^0& \text{in }& \Omega .\end{array}$$

$$\{\begin{array}{rlrl}\frac{\partial u_i}{\partial t}-\sum _{\alpha ,\beta }\frac{\partial }{\partial x_{\alpha }}\left[D_i^{\alpha \beta }(x)\left(\frac{\partial u_i}{\partial x_{\beta }}-\frac{\partial V_i}{\partial x_{\beta }}{u}_i\right)\right]& =0& \text{in }& \Omega ,\\ \sum _{\alpha ,\beta }{\nu }_{\alpha }\cdot \left[D_i^{\alpha \beta }(x)\left(\frac{\partial u_i}{\partial x_{\beta }}-\frac{\partial V_i}{\partial x_{\beta }}{u}_i\right)\right]& =0& \text{on }& \partial \Omega ,\\ u_i(0,\cdot )& =u_i^0& \text{in }& \Omega ,\end{array}$$

$$\lambda {\left|\xi \right|}^2\le \sum _{\alpha ,\beta }{D}_i^{\alpha \beta }(x){\xi }^{\alpha }{\xi }^{\beta }\le \Lambda {\left|\xi \right|}^2,$$

$$\Phi ,\Psi \in C^2{\left({\mathbb{R}}^m\right)}^{m\times m},\text{and}\Phi (0)=\Psi (0)=0.$$

$$\sum _{\alpha ,j}\left[\sum _{\beta }{\mathcal{A}}_{ij}^{\alpha \beta }(x,u^0)\frac{\partial u_j^0}{\partial x_{\beta }}-{\mathcal{B}}_{ij}^{\alpha }(x,u^0)u_j^0\right]{\nu }_{\alpha }=0.$$

$$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).$$

$${\Vert u\Vert }_{C^0(Q_T)}\le {\Gamma }_0{\Vert u\Vert }_{W(Q_T)}\le {\Gamma }_1{\Vert u_0\Vert }_{H^2(\Omega )}.$$

$$u\in C^{1+{\alpha }_0}([0,\tau ];L^2(\Omega ))\cap C^{{\alpha }_0}([0,\tau ];H^2(\Omega )),$$

$$\frac{\partial u}{\partial t},{\nabla }^{2}u\in C^{{\alpha }_1}([0,\tau ];C^0(\bar{\Omega })).$$

$$u\in C^{\epsilon }(J(u_0)\cap [0,T];C^0(\bar{\Omega })),$$

$$f\in H^r(0,T;H^s(\Omega )),\text{ for all }r,s\ge 0\text{ such that }r+\frac{s}{2}\le 1,$$

$$f\in C^{0,\eta }([0,T];C^{0,\theta }(\bar{\Omega })),\text{ for all }\eta ,\theta \ge 0\text{ such that }2\eta +\theta \le \frac{1}{2}.$$

$$⟨\omega ⟩\hat{g}\in L^2({\mathbb{R}}^{d+1}),\text{as well as}⟨\kappa {⟩}^2\hat{g}\in L^2({\mathbb{R}}^{d+1}).$$

$${\Vert g\Vert }_{X^{\prime }}\le 2{\Vert f\Vert }_X.$$

$${\Vert f\Vert }_{L^2(\Omega )}\le C_P{\Vert \nabla f\Vert }_{L^2(\Omega )}.$$

$${\Vert f_n\Vert }_{L^2(\Omega )}>n{\Vert \nabla f_n\Vert }_{L^2(\Omega )}.$$

$${\Vert \nabla f_{\infty }\Vert }_{L^2(\Omega )}=0,{\int }_S{f}_{\infty }(x)\mathrm{d}x=0,\text{and}{\Vert f_{\infty }\Vert }_{L^2(\Omega )}=1.$$

$${\Vert \nabla f\Vert }_{L^2(\Omega )}\le \frac{1}{1-C_P{\Vert \nabla V\Vert }_{\infty }}{\Vert \nabla f+f\nabla V\Vert }_{L^2(\Omega )},$$

$${\Vert \nabla f\Vert }_{L^2(\Omega )}$$

$$\le {\Vert \nabla f+f\nabla V\Vert }_{L^2(\Omega )}+{\Vert \nabla V\Vert }_{\infty }{\Vert f\Vert }_{L^2(\Omega )}$$

$$\le {\Vert \nabla f+f\nabla V\Vert }_{L^2(\Omega )}+C_P{\Vert \nabla V\Vert }_{\infty }{\Vert \nabla f\Vert }_{L^2(\Omega )}.$$

$$\left\{\begin{array}[]{rl}\dfrac{\partial u}{\partial t}-\nabla\cdot\left\{D(x)[(I+\delta\Phi(u))\nabla u+(\mathrm{diag}(\nabla V)+\delta\Psi(u))u]\right\}\kern-5.0pt&=0\quad\text{ in }\Omega,\\[5.0pt] \nu\cdot\left\{D(x)[(I+\delta\Phi(u))\nabla u+(\mathrm{diag}(\nabla V)+\delta\Psi(u))u]\right\}\kern-5.0pt&=0\quad\text{ on }\partial\Omega,\\[5.0pt] u(0,x)\kern-5.0pt&=u^{0}(x),\end{array}\right.$$

$$\left\{\begin{array}[]{rl}\nabla\cdot\left\{D(x)[(I+\delta\Phi(u^{*}))\nabla u^{*}+(\mathrm{diag}(\nabla V)+\delta\Psi(u^{*}))u^{*}]\right\}\kern-5.0pt&=0\quad\text{ in }\Omega,\\[5.0pt] \nu\cdot\left\{D(x)[(I+\delta\Phi(u^{*}))\nabla u^{*}+(\mathrm{diag}(\nabla V)+\delta\Psi(u^{*}))u^{*}]\right\}\kern-5.0pt&=0\quad\text{ on }\partial\Omega,\end{array}\right.$$

$$D={\left\{D_i^{\alpha \beta }\right\}}_{1\le i\le m}^{1\le \alpha ,\beta \le d},\text{and}\mathrm{diag}(\nabla V)={\left\{\frac{\partial V_i}{\partial x_{\beta }}\right\}}_{1\le i\le m}^{1\le \beta \le d},$$

$$\Phi =\{{\Phi }_{ij}{\}}_{1\le i,j\le m},\text{and}\Psi =\{{\Psi }_{ij}{\}}_{1\le i,j\le m}.$$

$$L_k(F,R):={\Vert F\Vert }_{C^k(\overline{B_R(0)})};$$

$$|F(w)|\le L_0(F,M),\text{ and }|F(w)-F(z)|\le L_1(F,M)|w-z|.$$

Full text

1. Introduction

1.1. 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

[TABLE]

which is a compact notation for the system

[TABLE]

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 2.

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, [14, 15, 31, 32]). 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 [24]. 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 [3]). They can be derived from a lattice based microscopic description, see, for instance, [10] and a subsequent formal Taylor expansion of the associated master equation, as considered in [17, 16]. Another derivation was performed in [14, 13], 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.

1.2. Gradient flow techniques

In recent years, gradient flow methods have been successfully employed to study certain families of cross-diffusion systems, see, among others, [26, 22, 21]. Techniques such as the boundedness-by-entropy principle introduced in [27], 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 [7]) 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 [6]. 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 [2] without finite volume constraints. Questions related to structural features of cross-diffusion systems to be interpreted as macroscopic gradient flows (as in [34]) or to have a strong solution (see [8]) 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.

1.3. 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 [9], 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 [19, 18]. 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, [20]. Also, convergence results are, to the best of our knowledge, restricted to one spatial dimension; cf. [29].

1.4. Summary of the main results

We consider a family of systems of PDEs with small cross-diffusion terms (namely system (1)) for which existence and uniqueness of solutions were investigated in [3], as recalled in Proposition 2.1. We extend the analysis presented in [3] by showing in Theorem 2.2 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 3.1. A special instance of system (1), to which we shall return in Section 4, was analysed in [12] using gradient flow techniques. This reduced PDE (namely problem (30)) 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 4.1 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 4.2 and Theorem 4.3 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 5.

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 2. Section 3 states the main exponential convergence result for the full system. In Section 4 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 5.

2. Analysis of systems with small cross-diffusion

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.

2.1. A well-posedness result

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:

[TABLE]

where $\delta$ is a small parameter such that

[TABLE]

for a suitable constant $\delta_{0}$ depending on $u_{0}$, $D$, $V$, $\Phi$, $\Psi$, and $\Omega$ (see [3] for further details on the smallness assumption). We consider the following initial-boundary value problem for system (1):

[TABLE]

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 (4) was proved in [3]. 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

[TABLE]

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 [3]).

1. H1

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

[TABLE]

for any $i=1,\ldots,m$. 2. H2

The dependence on $u$ of the nonlinear terms is (at least) of class $C^{2}$:

[TABLE] 3. H3

We assume that, for $i=1,\dots,m$, the potentials $V_{i}:\mathbb{R}^{d}\to\mathbb{R}$ are (at least) of class $C^{2}$. 4. H4

The domain $\Omega$ is bounded, connected and sufficiently smooth (e.g. of class $C^{2}$). The initial datum $u^{0}$ in (4) 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$:

[TABLE]

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

Definition 2.1**.**

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,

[TABLE]

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.

Proposition 2.1** (See [3]).**

Let hypotheses H1-H4 hold, then system (4) 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:

[TABLE]

2.2. Existence of global classical solutions

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

Theorem 2.2**.**

Let hypotheses H1-H4 hold and assume that $u_{0}$ is (at least) of class $C^{2}(\bar{\Omega})$. Then problem (4) 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 2.2, we will use the following two fundamental building blocks:

Proposition 2.3** (Existence of classical solutions for short time, [1]).**

Under the hypotheses of Theorem 2.2, 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 (4) has a unique solution $u$ in the interval $[0,\tau]$ which satisfies

[TABLE]

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

[TABLE]

Remark 1**.**

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

Proposition 2.4** (Criterion for global existence, [5]).**

Let all assumptions of Proposition 2.2 hold and consider a solution $u$ of problem (4). If an exponent $\varepsilon\in(0,1)$ (not depending on time) exists, such that for any $T>0$,

[TABLE]

then $u$ is a global solution, i.e., $J(u_{0})=[0,\infty)$.

Remark 2** (Notation).**

Proposition 2.4 is a rephrasing of Theorem 2 in [5] with $\theta=0$. Notice also that the “convective term”, indicated by $f$ in the notation of [5], is affine in the gradient in our case. Additionally, we do not use the notation $BUC^{\varepsilon}$ to denote the space of “bounded, uniformly $\varepsilon$-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 2.4 to be satisfied.

Lemma 2.5**.**

Let $f:Q_{T}\to\mathbb{R}$ be a function in $X=L^{2}(0,T;H^{2}(\Omega))\cap H^{1}(0,T;L^{2}(\Omega))$, then

[TABLE]

and, in turn,

[TABLE]

Proof.

Thanks to the higher order extensions for Sobolev functions, we can define $f$ on a larger cylindrical domain $R\subseteq\mathbb{R}^{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^{\prime}=L^{2}(\mathbb{R};H^{2}(\mathbb{R}^{d}))\cap H^{1}(\mathbb{R};L^{2}(\mathbb{R}^{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

[TABLE]

This means that $\left(\langle\omega\rangle+\langle\kappa\rangle^{2}\right)\hat{g}\in L^{2}(\mathbb{R}^{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$,

[TABLE]

Then the Hölder regularity follows from standard embeddings for fractional Sobolev spaces (see, e.g., [23]). In particular, for $r,s>\frac{1}{2}$, we take $\eta=r-\frac{1}{2}$ and $\theta=s-\frac{1}{2}$. ∎

Proof of Theorem 2.2.

Thanks to Proposition 2.3, we know that classical solutions exist for short times and that, as explained in [1], 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 2.4. 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 2.1, 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 2.5 in order to obtain obtain uniform Hölder estimates. Thus Proposition 2.4 allows us to conclude the proof. ∎

2.3. Poincaré’s inequality

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

Lemma 2.6** (A Poincaré–type inequality).**

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)\mathrm{d}x=0$, then there exists a constant $C_{P}>0$ (depending on $\Omega$ and $S$) such that

[TABLE]

Furthermore, given a family of functions $(f_{\ell})_{\ell\in\mathbb{R}}\subset H^{1}(\Omega)$, such that $\sigma_{0}=\mathrm{inf}_{\ell\in\mathbb{R}}\{{|\Omega|^{-1}}{|\mathrm{ess\,supp}(f_{\ell})|}\}\in(0,1)$, then inequality (10) holds for any function $f_{\ell}$ and the constant $C_{P}$ depends only on $\Omega$ and $\sigma_{0}$.

Proof.

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)\mathrm{d}x=0$ and

[TABLE]

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

[TABLE]

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. ∎

Remark 3**.**

Under the same assumptions of Lemma 2.6, we have

[TABLE]

provided that $C_{P}\left\|\nabla V\right\|_{\infty}<1$. More explicitly, we have

[TABLE]

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 $\mathrm{diam}(\Omega)$) or, alternatively, on the $L^{\infty}$-norm of the drift term $\nabla V$. We will use this relation in estimate (19).

3. Trend to equilibrium

In this section we investigate the equilibration behaviour of solutions to system (4). 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 (5) can be written as

[TABLE]

and the corresponding stationary problem:

[TABLE]

where we used the compact notation

[TABLE]

as well as

[TABLE]

Recall that $\Omega$ is a bounded domain of class $C^{2}$ in $\mathbb{R}^{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.

Remark 4**.**

Given a function $F:\mathbb{R}^{m}\to\mathbb{R}^{m\times m}$ of class at least $C^{k}$ ($k\geq 0$) such that $F(0)=0$, we define the quantity

[TABLE]

in particular, given two functions $w,z:\mathbb{R}^{d}\to\mathbb{R}^{m}$ such that $w,z\leq M$ for some $M>0$ and for a.e. $x\in\mathbb{R}^{d}$, then we have (by Taylor expansion)

[TABLE]

Remark 5**.**

We shall use the following notation:

[TABLE]

We now state the main result of this section:

Theorem 3.1** (Exponential convergence).**

Let assumptions H1-H4 hold and suppose that $u\in W(Q_{T})$ and $u^{*}\in W^{1,\infty}(\Omega)$ are solutions of problem (12) and (13) respectively. Let us denote by $M>0$ a constant not depending on $T$ such that

[TABLE]

for any $(t,x)\in[0,T]\times\Omega$. Furthermore, suppose that

[TABLE]

where $\Gamma_{M}$ is given by (20). Then the following inequality holds for any $t_{1},t_{2}\geq 0$:

[TABLE]

Proof.

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

[TABLE]

Testing system (17) against $h\exp(V)$ gives

[TABLE]

where

[TABLE]

Next we estimate the order $\delta$ terms in (18) and define

[TABLE]

In particular, making extensive use of Remarks 4, 5, and 3, we derive the following four estimates:

[TABLE]

as well as

[TABLE]

Since $\exp(V_{i})\geq 1$, combining equations (18) and (11) yields

[TABLE]

where we used Remark 3 in the Appendix and the constant $\Gamma_{M}$ is given by

[TABLE]

By assumption (15), we have the bound

[TABLE]

This bound in conjunction with estimate (19) yields

[TABLE]

Letting $g=h\exp(V)$, we obtain

[TABLE]

Using Poincaré’s inequality for $g$ and integrating in time over the interval $(t_{1},t_{2})\subseteq[0,\infty)$ yields

[TABLE]

Thanks to Gronwall’s inequality we obtain the following continuous-dependence estimate

[TABLE]

Finally, let us express this continuous-dependence estimate in terms of the original variables. To this end we recall that

[TABLE]

Substituting this relation into the preceding inequality, we deduce that

[TABLE]

and this concludes the proof. ∎

4. Study of an asymptotic gradient flow for one or two species

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 2. We will discuss the equilibration behaviour of their solutions and establish sharper $L^{\infty}$-bounds (recall that, thanks to Proposition 2.1, we already know that solutions are essentially bounded if $\delta<\delta_{0}$).

4.1. A cross-diffusion model for two interacting diffusing species

We consider the asymptotic gradient flow system for two species that was introduced by Bruna and Chapman in [13]. 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 [12]) or that it is close to being a diagonal, decoupled linear problem; cf. [4].

The asymptotic gradient flow system presented in [13] can be written in the general form of (4). In particular:

[TABLE]

and

[TABLE]

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 [13]. In [12] it was noted that system (4) with the nonlinearities in (21) and (22) 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

[TABLE]

where $\mathcal{M}\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

[TABLE]

The cross-diffusion system (4) with diffusion and drift matrices (21) and (22), can be rewritten as

[TABLE]

where $G$ is given by

[TABLE]

For a precise definition of all the parameters we refer to the derivation of the model in [13].

The discrepancy between the system above and the gradient-flow induced by (24) 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. [12]. 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. [12] and therefore $G\equiv 0$. Moreover, one can show that this is the only possibility for $G$ to vanish. We refer to [12] for more details as well as the existence proof of unique stationary solutions in both cases. More specifically, we recall the following result:

Proposition 4.1**.**

Let $\delta=\max\{\delta_{i,j}\}$, assume $0<\delta<\delta_{0}$ and suppose that the potentials in (4) satisfy $V_{i}\in H^{3}(\Omega)$. Then Theorem 3.1 applies to systems (23) and (25) 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

[TABLE]

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.

4.2. The scalar problem with a “frozen” component

In the following we focus on a special case of system (25), namely

[TABLE]

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 (30) is a scalar equation. It can be obtained from system (25) 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. [11] for more details.

Remark 6**.**

As equation (30) is a special case of (25), namely with immobile blue particles, Proposition 4.1 ensures that equation (30) 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

[TABLE]

From the discussion in Section 4.1 we know that equation (30) does not exhibit a full gradient flow structure, but only an asymptotic one. Hence, taking the entropy functional

[TABLE]

and mobility matrix $\mathcal{M}(r)=r(1-\delta_{2}b),$ we observe that equation (30) can be cast in the form (25) where $G(r,\delta)=rb(\delta_{1}\delta_{2}\nabla r+\delta_{2}\delta_{3}\nabla b).$

Remark 7**.**

As a consequence of Proposition 2.1, we have that Problem 30 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 4.1.

Definition 4.1** (Parabolic space).**

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”

[TABLE]

Similarly to Theorem 3.1, 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 (15)), as well as the use of the $L^{\infty}$-bound on $r$.

Theorem 4.2** (Exponential convergence to equilibrium).**

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 (30) 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$:

[TABLE]

for $\delta$ sufficiently small, satisfying relation (40) in the proof.

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

[TABLE]

Now let us consider a perturbation, $h$, of the stationary state $w^{*}$ of (32) (which exists by Proposition 4.1), namely

[TABLE]

where $\int_{\Omega}h(x,t)=0$ since the mass is conserved. We observe that $h$ satisfies the following equation:

[TABLE]

We test equation (33) against $h$ to obtain

[TABLE]

Integration in time over $(t_{1},t_{2})\subset[0,T]$ yields

[TABLE]

Using the fact that $\nabla V,w^{*},\nabla w^{*},\nabla b\in L^{\infty}(\Omega),$ we deduce that

[TABLE]

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

[TABLE]

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

[TABLE]

and we substitute equations (36) and (37) into equation (35), which, in turn, becomes

[TABLE]

Recalling that $h+w^{*}=w$ is non-negative, we choose $\delta$ sufficiently small so that the following relation is satisfied:

[TABLE]

which implies

[TABLE]

Thus we have obtained the following inequality:

[TABLE]

Using the Poincare’s inequality and the fact that $V_{l}<V<V_{u}$, we have

[TABLE]

Thanks to Gronwall’s inequality we obtain the following estimate:

[TABLE]

Finally, let us switch back to the original variables by recalling the fact that $h=w-w^{*}$. We obtain

[TABLE]

Since $V$ is bounded from above and below by assumption, a simple change of variable $r=we^{-V}$ implies the exponential convergence to equilibrium for the variable $r$ as claimed in the statement, which concludes the proof. ∎

4.3. 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 (30), in dimension $d=2$ (the same result can be obtained in one dimension). We adapt the strategy presented in [28], Chapter V111Namely, 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).. In the context of elliptic equations, this approach is referred to as De Giorgi’s method (see, e.g., [25]).

Theorem 4.3** (Maximum principle).**

Let $d=2$ and suppose that $\delta$ is sufficiently small and let $r\in Z$ be a solution of problem (30). Then we have that

[TABLE]

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 4.3 uses several technical results which we shall recall in the following. Notice that, if $x$ is close to a maximum point of $V$, estimate (43) 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 4.3.

Lemma 4.4**.**

Let $d=2$ and suppose that $z\in Z$ is a weak solution of the problem

[TABLE]

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

1. (I)

$\omega:\mathbb{R}^{d}\to\mathbb{R}$* satisfies $0<\mu\leq\omega(x)<\mu^{-1}$, for a.e. $x\in\Omega$,*

** 2. (II)

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

** 3. (III)

$F(x,z)=zF_{1}(x)+z^{2}F_{2}(x)$, where $F_{i}\in L^{\infty}(\Omega)$, and $i=1,2$,

** 4. (IV)

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

** 5. (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 (56), such that

[TABLE]

Before proving Lemma 4.4, we state a preliminary result from Chapter II in [28]. It can be proved directly by induction.

Lemma 4.5**.**

Suppose that a sequence $(a_{n})_{n\in\mathbb{N}}\subset\mathbb{R}_{+}$, of non-negative numbers satisfies the recursive relation

[TABLE]

for $\kappa,\varepsilon>0$ and $\zeta>1$. Then we have that

[TABLE]

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

Proof of Lemma 4.4.

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

[TABLE]

for $k\in\mathbb{R}_{+}$ 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

[TABLE]

for $n\in\mathbb{N}$ and suitable $M>0$ as well as $q>0$. Subsequently, we will apply Lemma 4.5 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\mathbb{R}_{+}$, $k>\max(\tilde{M},1)$. It is easily verified that $(z-k)_{+}\in Z$ (see, e.g. Theorem 2.1.11 in [35] or Lemma 4.4 in Chapter II of [28]). Thus we may test equation (47) against the truncated function $(z-k)_{+}$ and obtain

[TABLE]

Rearranging this relation we get

[TABLE]

**Step 1: upper bound for the RHS in terms of $|\{z\geq k\}|$

**Let us focus on the estimate of the right-hand side:

[TABLE]

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 (48) may be estimated further by

[TABLE]

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}\mathrm{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 2.6 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:

Remark 8**.**

Let $f\in L^{1}(\Omega)$ be any non-negative function such that $\int_{\Omega}f\mathrm{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

[TABLE]

Then we have, for $q=\frac{m}{\alpha|\Omega|}$ and an arbitrary $\alpha\in(0,1)$,

[TABLE]

We recall that integrating (47) over the whole domain gives $\int_{\Omega}\omega z\mathrm{d}x=\int_{\Omega}\omega z_{0}\mathrm{d}x$, for any $t>0$, as $z$ is continuous with respect to time. We can now argue as in Remark 8 with $m=\int_{\Omega}\omega z_{0}\mathrm{d}x$, $f=\omega z$, $\alpha\in(0,1)$ and $k\geq\frac{1}{\alpha\mu|\Omega|}\int_{\Omega}\omega z_{0}\mathrm{d}x$, to the set \{x\,\big{|}\,\omega z\leq\mu k\}, which, by assumption (I), 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 2.6 (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 2.6); namely we get:

[TABLE]

This may be further estimated upon observing that

[TABLE]

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:

[TABLE]

Notice that for the remaining term in estimate (49), we have

[TABLE]

Hence estimate (49) becomes

[TABLE]

Altogether we have obtained the inequality:

[TABLE]

where $C_{F}=\frac{C_{V}}{2\mu}\left[5LC_{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

[TABLE]

Supposing that $\delta$ is such that

[TABLE]

we have

[TABLE]

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):

[TABLE]

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)$,

[TABLE]

where $C_{e}$ does not depend on time. By Chebyshev’s inequality, for any $a\geq 1$, we obtain

[TABLE]

whence, upon integration, we obtain

[TABLE]

Choosing $a=4-\nu$, we get

[TABLE]

**Step 3: combining all previous estimates

**Combining (52), (53), and (54) we obtain

[TABLE]

where ${\tilde{\Gamma}_{\delta}}=\max\{1,C_{e}^{2}\}\delta\Gamma$. For $M>\bar{M}$, consider the increasing sequence defined by

[TABLE]

and

[TABLE]

for any $n\in\mathbb{N}$. In order to apply Lemma 4.5 we need to define a recursive relation for the sequence $(a_{n})_{n\in\mathbb{N}}$. Let us set

[TABLE]

then we observe that

[TABLE]

we observe that, from (55),

[TABLE]

Thus the sequence $(a_{n})_{n\in\mathbb{N}}$ satisfies the following relation

[TABLE]

Let us set $\varepsilon=1-\nu/2$, $\zeta=4$, and $\kappa=16\tilde{\Gamma}_{\delta}$ and apply Lemma 4.5. We deduce that $a_{n}\to 0$ as $n\to\infty$, and that the following inequality holds:

[TABLE]

provided that

[TABLE]

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 (53) and (54) and upon choosing $h=\theta\bar{M}$ and $k=\bar{M}$, we also know that,

[TABLE]

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

[TABLE]

which, in turn, gives

[TABLE]

such that the upper bound for $a_{0}$ is met. We conclude that

[TABLE]

after passing $n\to\infty$, provided that

[TABLE]

and consequently we have obtained that, for a.e. $t,x$,

[TABLE]

∎

Remark 9** (Higher dimensions).**

Lemma 4.4 does not hold in general for $d>2$. For example, estimate (50) 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 (47). Sufficient conditions for $L^{\infty}$ estimates are presented in [28], 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 4.4 at hand, the proof of Theorem 4.3 is immediate.

Proof of Theorem 4.3..

We apply Lemma 4.4 making the following choices:

[TABLE]

Recall that rewriting equation (30) in terms of the new unknown $z=u\exp(V)$ we obtain

[TABLE]

with suitable initial and boundary conditions. ∎

5. Numerical Simulations for the Scalar Problem

Let us recall the equation

[TABLE]

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.,

[TABLE]

This formulation is the basis of the finite volume scheme used for the following section. The scheme is based on the schemes studied in [18, 19, 9]. 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 (23) with entropy (31). 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.,

[TABLE]

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

[TABLE]

for $i\in I$. Here the numerical fluxes are given by

[TABLE]

for $i=1,\ldots,N-1$, and

[TABLE]

as well as

[TABLE]

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.,

[TABLE]

5.1. 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 1. This choice forces the mobile species to penetrate the immobile species as it migrates towards the minimum of the external potential. In Figure 2 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)=\frac{1}{2}r_{0}(x-3)+\frac{1}{2}r_{0}(x+3)$, with $r_{0}$ as above; cf. Figure 1. Figure 3 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 4 we show the exponential convergence of the solution towards equilibrium.

5.2. 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 [13].

5.2.1. 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(5x)+\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 5.

After some time the evolution slows down as the mobile species reaches a stationary state. In Figure 6 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.

5.2.2. 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

[TABLE]

with $\sigma=0.3$. The immobile species is fixed and given by

[TABLE]

where $\rho(x,y)=\sqrt{(x+7/4)^{2}+(y+7/4)^{2}}$ and $N$ normalises the mass to one; cf. Figure 7.

Below we consider two different choices of external potentials: $V_{\mathrm{s}}(x)=\frac{1}{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 8 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 9 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 9. Compared to Figure 8 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.