Vanishing cross-diffusion limit in a Keller-Segel system with additional cross-diffusion

TL;DR

This paper rigorously analyzes the vanishing cross-diffusion limit in Keller-Segel systems, demonstrating global existence of solutions and convergence to classical models, supported by theoretical proofs and numerical experiments.

Contribution

It provides the first rigorous proof of the vanishing cross-diffusion limit in Keller-Segel systems with additional cross-diffusion, including convergence rates and numerical validation.

Findings

Cross-diffusion stabilizes the system and ensures global solutions.

Solutions converge to classical Keller-Segel models as cross-diffusion vanishes.

Numerical experiments illustrate cell aggregation behavior as a function of cross-diffusion.

Abstract

Keller-Segel systems in two and three space dimensions with an additional cross-diffusion term in the equation for the chemical concentration are analyzed. The cross-diffusion term has a stabilizing effect and leads to the global-in-time existence of weak solutions. The limit of vanishing cross-diffusion parameter is proved rigorously in the parabolic-elliptic and parabolic-parabolic cases. When the signal production is sublinear, the existence of global-in-time weak solutions as well as the convergence of the solutions to those of the classical parabolic-elliptic Keller--Segel equations are proved. The proof is based on a reformulation of the equations eliminating the additional cross-diffusion term but making the equation for the cell density quasilinear. For superlinear signal production terms, convergence rates in the cross-diffusion parameter are proved for local-in-time smooth solutions (since finite-time blow up is possible). The proof is based on careful $H^s(\Omega)$ estimates and a variant of the Gronwall lemma. Numerical experiments in two space dimensions illustrate the theoretical results and quantify the shape of the cell aggregation bumps as a function of the cross-diffusion parameter.

Equations

$${\partial }_t{\rho }_{\delta }=\mathrm{div}(\nabla {\rho }_{\delta }-{\rho }_{\delta }\nabla c_{\delta }),\epsilon {\partial }_t{c}_{\delta }=\Delta c_{\delta }+\delta \Delta {\rho }_{\delta }-c_{\delta }+{\rho }_{\delta }^{\alpha }\text{mbox}in\Omega ,t>0,$$

$$\nabla {\rho }_{\delta }\cdot \nu =\nabla c_{\delta }\cdot \nu =0\text{mbox}on\partial \Omega ,t>0,{\rho }_{\delta }(0)={\rho }^0,\epsilon c_{\delta }(0)=\epsilon c^0\text{mbox}in\Omega ,$$

$${\partial }_t\rho =\mathrm{div}(\nabla \rho -\rho \nabla c),\epsilon {\partial }_{t}c=\Delta c-c+{\rho }^{\alpha }\text{mbox}in\Omega ,t>0,$$

$${\partial }_t{\rho }_{\delta }\in L^2(0,T;H^1(\Omega {)}^{\prime }),{\rho }_{\delta }\in L^{\infty }(0,T;L^3(\Omega )),c_{\delta }+\delta {\rho }_{\delta }\in L^{\infty }(0,T;W^{1,p}(\Omega )),$$

$${\rho }_{\delta }\to \rho$$

$$\nabla {\rho }_{\delta }\rightharpoonup \nabla \rho$$

$$c_{\delta }+\delta {\rho }_{\delta }{\rightharpoonup }^{\ast }c$$

$${\partial }_t{\rho }_{\delta }=\mathrm{div}((1+\delta {\rho }_{\delta })\nabla {\rho }_{\delta }-{\rho }_{\delta }\nabla v_{\delta }),-\Delta v_{\delta }+v_{\delta }=\delta {\rho }_{\delta }+{\rho }_{\delta }^{\alpha }\text{mbox}in\Omega ,t>0,$$

$$\nabla {\rho }_{\delta }\cdot \nu =\nabla v_{\delta }\cdot \nu =0\text{mbox}in\partial \Omega ,{\rho }_{\delta }(0)={\rho }^0\text{mbox}in\Omega .$$

$$\frac{d{H}_1}{dt}+4{\int }_{\Omega }|\nabla {\rho }_{\delta }^{1/2}{|}^{2}dx+\delta {\int }_{\Omega }|\nabla {\rho }_{\delta }{|}^{2}dx\le {\int }_{\Omega }(\delta {\rho }_{\delta }^2+{\rho }_{\delta }^{\alpha +1})dx.$$

$$\delta {\int }_{\Omega }{\rho }_{\delta }^{2}dx\le \frac{\delta }{2}\Vert \nabla {\rho }_{\delta }{\Vert }_{L^2(\Omega )}^2+C\Vert {\rho }_{\delta }{\Vert }_{L^1(\Omega )}^2.$$

$$\frac{d{H}_p}{dt}+{\int }_{\Omega }|\nabla {\rho }_{\delta }^{p/2}{|}^{2}dx+\delta {\int }_{\Omega }|\nabla {\rho }_{\delta }^{(p+1)/2}{|}^{2}dx\le C\Vert {\rho }_{\delta }{\Vert }_{L^{(p+1)/2}(\Omega )}^{(p+1)/2}+C\Vert {\rho }_{\delta }^{p/2}{\Vert }_{L^1(\Omega )}^{\beta (p)},$$

$$\Vert ({\rho }_{\delta }-\rho ,c_{\delta }-c){\Vert }_{L^{\infty }(0,T;H^2(\Omega ))}\le C{\delta }^{1-\lambda }.$$

$$\Gamma (t)=\Vert ({\rho }_R,c_R)(t){\Vert }_{H^2(\Omega )}^2,G(t)=\Vert ({\rho }_R,c_R)(t){\Vert }_{H^2(\Omega )}^2+\Vert \nabla \Delta ({\rho }_R,c_R)(t){\Vert }_{L^2(\Omega )}^2.$$

$$\Gamma (t)+C_1{\int }_0^{t}G(s)ds\le C_2{\int }_0^t(\Gamma (s)+\Gamma (s{)}^{\max \{2,\alpha \}})ds+C_2{\int }_0^t\Gamma (s)G(s)ds+C_2{\delta }^2,$$

$$\displaystyle\frac{1}{\tau}\int_{\Omega}(\rho^{k}-\rho^{k-1})\phi dx+\int_{\Omega}\big{(}(1+\delta\rho^{k})\nabla\rho^{k}-\rho^{k}\nabla v^{k}\big{)}\cdot\nabla\phi dx$$

$$\displaystyle\phantom{xxxx}{}=-\eta\int_{\Omega}\big{(}\Delta w^{k}\Delta\phi+\delta^{-2}|\nabla w^{k}|^{2}\nabla w^{k}\cdot\nabla\phi+w^{k}\rho^{k}\phi\big{)}dx,$$

$$\displaystyle\int_{\Omega}(\nabla v^{k}\cdot\nabla\theta+v^{k}\theta)dx=\int_{\Omega}\big{(}\delta\rho^{k}+(\rho^{k})^{\alpha}\big{)}\theta dx$$

$$\frac{1}{\tau }{\int }_{\Omega }$$

$$\displaystyle\phantom{xx}{}+\frac{\eta}{\delta}\int_{\Omega}\big{(}(\Delta w^{k})^{2}+\delta^{-2}|\nabla w^{k}|^{4}+(w^{k})^{2}\rho^{k}\big{)}dx$$

$$\displaystyle\leq\int_{\Omega}\big{(}-v^{k}\rho^{k}+\delta(\rho^{k})^{2}+(\rho^{k})^{\alpha+1}\big{)}dx\leq\int_{\Omega}\big{(}\delta(\rho^{k})^{2}+(\rho^{k})^{\alpha+1}\big{)}dx,$$

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

$$\displaystyle\leq\mu\big{(}\|\nabla u\|_{L^{2}(\Omega)}^{2}+\|u\|_{L^{1}(\Omega)}^{2}\big{)}+C(\mu)\|u\|_{L^{1}(\Omega)}^{2}=\mu\|\nabla u\|_{L^{2}(\Omega)}^{2}+C(\mu)\|u\|_{L^{1}(\Omega)}^{2}.$$

$$\delta {\int }_{\Omega }({\rho }^k{)}^{2}dx\le \frac{\delta }{4}\Vert \nabla {\rho }^k{\Vert }_{L^2(\Omega )}^2+C\Vert {\rho }^k{\Vert }_{L^1(\Omega )}^2,$$

$$\int_{\Omega}(\rho^{k})^{\alpha+1}dx\leq\int_{\Omega}\big{(}(\rho^{k})^{2}+1\big{)}dx\leq\frac{\delta}{4}\|\nabla\rho^{k}\|_{L^{2}(\Omega)}^{2}+C\|\rho^{k}\|_{L^{1}(\Omega)}^{2}+C(\Omega).$$

$$\frac{1}{\tau }{\int }_{\Omega }$$

$$\displaystyle\phantom{xx}{}+\frac{\eta}{\delta}\int_{\Omega}\big{(}(\Delta w^{k})^{2}+\delta^{-2}|\nabla w^{k}|^{4}+(w^{k})^{2}\rho^{k}\big{)}dx\leq C(\Omega).$$

$$\frac{1}{\tau }{\int }_{\Omega }$$

$$\displaystyle\phantom{xx}{}+\eta p\int_{\Omega}\big{(}\Delta w^{k}\Delta(\rho^{k})^{p-1}+\delta^{-2}|\nabla w^{k}|^{2}\nabla w^{k}\cdot\nabla(\rho^{k})^{p-1}+w^{k}(\rho^{k})^{p}\big{)}dx$$

$$\displaystyle\leq(p-1)\int_{\Omega}\big{(}\delta(\rho^{k})^{p+1}+(\rho^{k})^{p+\alpha}\big{)}dx.$$

$$\delta {\int }_{\Omega }({\rho }^k{)}^{p+1}dx=\delta \Vert ({\rho }^k{)}^{(p+1)/2}{\Vert }_{L^2(\Omega )}^2\le \frac{\delta }{2}\Vert \nabla ({\rho }^k{)}^{(p+1)/2}{\Vert }_{L^2(\Omega )}^2+C\Vert ({\rho }^k{)}^{(p+1)/2}{\Vert }_{L^1(\Omega )}^2.$$

$${\int }_{\Omega }({\rho }^k{)}^{\alpha +p}dx$$

$$\displaystyle=C\big{(}\|\nabla(\rho^{k})^{p/2}\|_{L^{2}(\Omega)}^{2}+\|(\rho^{k})^{p/2}\|_{L^{1}(\Omega)}^{2}\big{)}^{(\alpha+p)\sigma/p}\|(\rho^{k})^{p/2}\|_{L^{1}(\Omega)}^{2(\alpha+p)(1-\sigma)/p}$$

$$\displaystyle\leq\frac{1}{2}\big{(}\|\nabla(\rho^{k})^{p/2}\|_{L^{2}(\Omega)}^{2}+\|(\rho^{k})^{p/2}\|_{L^{1}(\Omega)}^{2}\big{)}+C\|(\rho^{k})^{p/2}\|_{L^{1}(\Omega)}^{\beta(p)},$$

$$\eta$$

$$\displaystyle=\frac{\eta}{\delta}p(p-1)\int_{\Omega}(\rho^{k})^{p-1}\bigg{[}\bigg{(}\Delta w^{k}+\frac{p-1}{2\delta}|\nabla w^{k}|^{2}\bigg{)}^{2}+\frac{1}{4\delta^{2}}\big{(}4-(p-1)^{2}\big{)}|\nabla w^{k}|^{4}\bigg{]}dx$$

$$+\eta p{\int }_{\Omega }{w}^k{e}^{p{w}^k/\delta }dx.$$

Full text

The first two authors acknowledge partial support from the Austrian Science Fund (FWF), grants P30000, W1245, and F65. The third author is partially supported by the National Natural Science Foundation of China (NSFC), grants 11831003, 11771031, 11531010, and by NSF of Qinghai Province, grant 2017-ZJ-908.

1. Introduction

Chemotaxis describes the directed movement of cells in response to chemical gradients and may be modeled by the (Patlak–) Keller–Segel equations [24, 29]. The aggregation of cells induced by the chemical concentration is counter-balanced by cell diffusion. If the cell density is sufficiently large, the chemical interaction dominates diffusion and results in a blow-up of the cell density. However, a single point blow-up is not very realistic from a biological view point, due to the finite size of the cells. Therefore, chemotaxis models that avoid blow-up were suggested in the literature. Possible approaches are bounded chemotaxis sensibilities to avoid over-crowding [7, 12], degenerate cell diffusion [8, 25, 26], death-growth terms [6, 32], or additional cross-diffusion [10, 18]. In this paper, we continue our study of the Keller–Segel system with additional cross-diffusion, which allows for global weak solutions [10, 18]. The question how well the solutions approximate those from the original Keller–Segel system remained open. Here, we will answer this question by performing rigorously the vanishing cross-diffusion limit and giving an estimate for the difference of the respective solutions.

The equations for the cell density $\rho_{\delta}(x,t)$ and the chemical concentration $c_{\delta}(x,t)$ are given by

[TABLE]

subject to the no-flux and initial conditions

[TABLE]

where $\Omega\subset{\mathbb{R}}^{d}$ ($d=2,3$) is a bounded domain, $\nu$ is the exterior unit normal vector of $\partial\Omega$, $\delta>0$ describes the strength of the additional cross-diffusion, and the term $\rho_{\delta}^{\alpha}$ with $\alpha>0$ is the nonlinear signal production. The case $\varepsilon=1$ is called the parabolic-parabolic case, while $\varepsilon=0$ refers to the parabolic-elliptic case. The existence of global weak solutions to (1)–(2) was proved in [18] in two space dimensions and in [10] in three space dimensions (with degenerate diffusion).

Setting $\delta=0$, we obtain the Keller–Segel system with nonlinear signal production. In the classical Keller–Segel model, the signal production is assumed to be linear, $\alpha=1$. In order to deal with the three-dimensional case, we need sublinear signal productions, $\alpha<1$. In fact, it is shown in [33] that $\alpha=2/d$ is the critical value for global existence versus finite-time blow-up in a slightly modified Keller–Segel system with $\delta=0$. We show that also for $\delta>0$, the condition $\alpha<2/d$ guarantees the global existence of weak solutions, while numerical results indicate finite-time blow-up if $\alpha>2/d$.

We are interested in the limit $\delta\to 0$ in (2) leading to the Keller–Segel equations

[TABLE]

with the initial and boundary conditions (2). We refer to the reviews [3, 16] for references concerning local and global solvability and results for variants of this model. In the following, we recall only some effects related to the blow-up behavior.

In the parabolic-elliptic case, a dichotomy arises for (3) in two space dimensions: If the initial mass $M:=\|\rho^{0}\|_{L^{1}(\Omega)}$ is smaller than $8\pi$, the solutions are global in time, while they blow up in finite time if $M>8\pi$ and the second moment of the initial datum is finite (see, e.g., [5]). The condition on the second moment implies that the initial density is highly concentrated around some point. It is necessary in the sense that there exists a set of initial data with total mass larger than $8\pi$ such that the corresponding solutions are global [2].

In the parabolic-parabolic case, again in two space dimensions and with finite second moment, the solutions exist globally in time of $M<8\pi$ [9]. However, in contrast to the parabolic-elliptic case, the threshold value for $M$ is less precise, and solutions with large mass can exist globally. In dimensions $d\geq 3$, a related critical phenomenon occurs: The solutions exist globally in time if $\|\rho^{0}\|_{L^{d/2}(\Omega)}$ is sufficiently small, but they blow up in finite time if the total mass is large compared to the second moment [11].

In this paper, we prove two results. The first one is the convergence of the solutions to the parabolic-elliptic model (1) with $\varepsilon=0$ to a solution to the parabolic-elliptic Keller–Segel system (3) with $\varepsilon=0$. Since we impose a restriction on the parameter $\alpha$ as mentioned above, this result holds globally in time. We set $Q_{T}=\Omega\times(0,T)$.

Theorem 1** (Convergence for the parabolic-elliptic model).**

Let $\Omega\subset{\mathbb{R}}^{d}$ $(d=2,3)$ be bounded with $\partial\Omega\in C^{1,1}$, $T>0$, $\delta>0$, $\varepsilon=0$, and $0\leq\rho^{0}\in L^{\infty}(\Omega)$. If $\alpha<2/d$, there exists a weak solution $(\rho_{\delta},c_{\delta})\in L^{2}(0,T;H^{1}(\Omega))^{2}$ to (1)-(2) satisfying

[TABLE]

for any $p<\infty$. Furthermore, as $\delta\to 0$,

[TABLE]

and $(\rho,c)$ solves (2)–(3), and it holds that $\rho$, $c\in L^{\infty}(0,T;L^{\infty}(\Omega))$.

The idea of the proof is to reformulate (1) via $v_{\delta}:=c_{\delta}+\delta\rho_{\delta}$ as the system

[TABLE]

together with the initial and boundary conditions

[TABLE]

This reformulation was already used in [18] to prove the existence of weak solutions in the two-dimension case with $\alpha=1$. It transforms the asymptotically singular problem to a quasilinear parabolic equation, thus simplifying considerably the asymptotic limit problem. Still, we need estimates uniform in $\delta$ to apply compactness arguments. For this, we use the “entropy” functional $H_{1}(\rho_{\delta})=\int_{\Omega}\rho_{\delta}(\log\rho_{\delta}-1)dx$:

[TABLE]

By the Gagliardo–Nirenberg inequality, the first term on the right-hand side can be estimated as (see Section 2 for the proof)

[TABLE]

The first term on the right-hand side is absorbed by the left-hand side of (6) and the second term is bounded since the total mass is conserved. To estimate the second term on the right-hand side of (6), we need another “entropy” functional $H_{p}(\rho)=\int_{\Omega}\rho_{\delta}^{p}dx$ with $p=2$ or $p=3$, leading to

[TABLE]

where $\beta(p)\geq 2$ is some function depending on $p$. If $p=2$, the last term is the total mass which is bounded uniformly in time. Moreover, the estimate for $\rho_{\delta}^{1/2}$ in $H^{1}(\Omega)$ from (6) implies that $\rho_{\delta}$ is bounded in $W^{1,1}(\Omega)\hookrightarrow L^{3/2}(\Omega)$ such that the first term on the right-hand side of (7) is uniformly bounded as well. Higher-order integrability is then obtained for $p=3$.

Clearly, these arguments are formal. In particular, the estimate for $H_{1}$ requires the test function $\log\rho_{\delta}$ in (1) which may be not defined if $\rho_{\delta}=0$. Therefore, we consider an implicit Euler discretization in time with parameter $\tau>0$ and an elliptic regularization in space with parameter $\eta>0$ to prove first the existence of approximate weak solutions with strictly positive $\rho_{\delta}$. This is done by using the entropy method of [18]. The approximate solutions also satisfy the $\delta$-uniform bounds, and they hold true when passing to the limit $(\eta,\tau)\to 0$. Then the limit $\delta\to 0$ can be performed by applying the Aubin–Lions lemma and weak compactness arguments.

The second result is concerned with the derivation of a convergence rate both in the parabolic-parabolic and parabolic-elliptic case. For $\alpha\geq 1$, we cannot generally expect global solutions. In this case, it is natural to consider local solutions. The technique requires smooth solutions so that we need some regularity assumptions on the initial data.

Theorem 2** (Convergence rates).**

Let $\Omega\subset{\mathbb{R}}^{d}$ $(d\leq 3)$ be a bounded domain with smooth boundary and let $(\rho^{0},c^{0})\in(W^{1,p}(\Omega))^{2}$ for $p>d$ if $\varepsilon=1$ and $\rho^{0}\in C^{2+\gamma}(\overline{\Omega})$ for some $\gamma\in(0,1)$ if $\varepsilon=0$. Furthermore, let $\alpha=1$ or $\alpha\geq 2$ and let $(\rho_{\delta},c_{\delta})$ and $(\rho,c)$ be (weak) solutions to (1) and (3), (2), respectively, with the same initial data. Then these solutions are smooth locally in time and there exist constants $C>0$ and $\delta_{0}>0$ such that for all $0<\delta<\delta_{0}$ and $\lambda>0$,

[TABLE]

The theorem is proved as in, e.g., [19] by deriving carefully $H^{s}(\Omega)$ estimates for the difference $(\rho_{R},c_{R}):=(\rho_{\delta}-\rho,c_{\delta}-c)$. The index $s\in{\mathbb{N}}$ is chosen such that we obtain $L^{\infty}(\Omega)$ estimates which are needed to handle the nonlinearities. If $\varepsilon=1$, we introduce the functions

[TABLE]

The aim is to prove the inequality

[TABLE]

where $C_{1}$, $C_{2}>0$ are constants independent of $\delta$. This inequality allows us to apply a variant of Gronwall’s lemma (see Lemma 4 in the Appendix), implying (8). In the parabolic-elliptic case $\varepsilon=0$, the functions $\Gamma(t)$ and $G(t)$ are defined without $c_{R}$, and the final inequality contains the additional integral $\int_{0}^{t}\Gamma(s)^{\alpha+1}G(s)ds$, which is still covered by Lemma 4. The condition $\alpha=1$ or $\alpha\geq 2$ comes from the fact that the derivative of the mapping $s\mapsto s^{\alpha}$ is Hölder continuous exactly for these values. The numerical results in Section 4 indicate that the convergence result may still hold for $\alpha\in(1,2)$.

The optimal convergence rate is expected to be one. Theorem 2 provides an almost optimal rate. The reason for the non-optimality comes from the variant of the nonlinear Gronwall lemma proved in Lemma 4. We conjecture that an optimal rate holds (changing the constants in Lemma 4), but since this issue is of less interest, we did not explore it further.

The theoretical results are illustrated by numerical experiments, using the software tool NGSolve/Netgen. We remark that numerical tests for model (1) were already performed in [4, 18]. For positive values of $\delta$, the (globally existing) cell density forms bumps at places where the solution to the classical Keller–Segel system develops an $L^{\infty}(\Omega)$ blow up. Compared to [4, 18], we investigate the dependence of the shape of the bumps on $\delta$. In a radially symmetric situation, it turns out that the radius of the bump (more precisely the diameter of a level set $\rho_{\delta}\approx 0$) behaves like $\delta^{a}$ with $a\approx 0.43$, and the maximum of the bump behaves like $\delta^{-b}$ with $b\approx 1.00$.

The paper is organized as follows. Theorem 1 is proved in Section 2, while Theorem 2 is shown in Section 3. The numerical experiments for model (1) are performed in Section 4. Some technical tools, including the nonlinear Gronwall inequality, are recalled in the Appendix.

2. Proof of Theorem 1

We prove Theorem 1 by an approximation procedure and by deriving the uniform bounds from discrete versions of the entropy inequalities (6) and (7).

Step 1: Solution of a regularized system and entropy estimates. We show the existence of solutions to a discretized and regularized version of (4). For this, let $N\in{\mathbb{N}}$ and $\tau=T/N$, and set $\rho(w)=\exp(w/\delta)$. This means that we transform $w=\delta\log\rho$. Let $w^{k-1}\in H^{2}(\Omega)$ and $v^{k-1}\in H^{1}(\Omega)$ be given and set $\rho^{j}=\rho(w^{j})$ for $j=k,k-1$. Consider for given $\tau>0$ and $\eta>0$ the regularized system

[TABLE]

for $\phi\in H^{2}(\Omega)$ and $\theta\in H^{1}(\Omega)$. The time discretization is needed to handle issues due to low time regularity, while the elliptic regularization guarantees $H^{2}(\Omega)$ solutions which, by Sobolev embedding (recall that $d\leq 3$), are bounded. The higher-order gradient term $|\nabla w^{k}|^{2}\nabla w^{k}\cdot\nabla\phi$ is necessary to derive $L^{p}(\Omega)$ estimates. The existence of a solution $w^{k}\in H^{2}(\Omega)$, $0\leq v^{k}\in H^{1}(\Omega)$ follows from the techniques used in the proof of Proposition 3.1 in [18] employing the Leray–Schauder fixed-point theorem. Since these techniques are rather standard now, we omit the proof and refer to [18, 21, 22] for similar arguments.

Inequality $w^{k}\rho^{k}=w^{k}e^{w^{k}/\delta}\geq e^{w^{k}/\delta}-1=\rho^{k}-1$ allows us to show as in [18, page 1004] that the total mass $\|\rho^{k}\|_{L^{1}(\Omega)}$ is bounded uniformly in $\delta$.

Entropy estimates are derived from (9)–(10) by choosing the test functions $\phi=w^{k}/\delta=\log\rho^{k}$ and $\theta=\rho^{k}$, respectively, and adding both equations. Then the terms involving $\nabla v^{k}$ cancel and after some elementary computations, we end up with

[TABLE]

where $h(s)=s(\log s-1)$ for $s\geq 0$. Using the Gagliardo–Nirenberg inequality with $\sigma=d/(d+2)$, the Poincaré–Wirtinger inequality, and then the Young inequality with $p=1/\sigma$, $p^{\prime}=1/(1-\sigma)$, it holds for any $u\in H^{1}(\Omega)$ and $\mu>0$ that

[TABLE]

We deduce from this inequality that the first term on the right-hand side of (11) can be estimated as

[TABLE]

where here and in the following, $C>0$ denotes a constant, independent of $\eta$, $\tau$, and $\delta$, with values varying from line to line. The first term on the right-hand side can be absorbed by the left-hand side of (11), while the second term is bounded. Using $s^{\alpha+1}\leq s^{2}+1$ for $s\geq 0$ (since $\alpha<1$) and (12), the second term on the right-hand side of (11) becomes

[TABLE]

Inserting these estimations into (11), we conclude that

[TABLE]

Step 2: Further uniform estimates. The estimates from (13) are not sufficient for the limit $\delta\to 0$, therefore, we derive further uniform bounds. Let $p=2$ or $p=3$. We choose the admissible test functions $p(\rho^{k})^{p-1}$ and $(p-1)(\rho^{k})^{p}$ in (9)–(10), respectively, and add both equations. The convexity of $s\mapsto s^{p}$ implies that $s^{p}-t^{p}\leq p(s-t)s^{p-1}$. Then, observing that the terms involving $\nabla v^{k}$ cancel, we find that

[TABLE]

By (12), we find that

[TABLE]

The estimate for $(\rho^{k})^{\alpha+p}$ requires that $\alpha<2/d$. Indeed, we deduce from the Gagliardo–Nirenberg inequality with $\sigma=d(2\alpha+p)/((d+2)(\alpha+p))$ and similar arguments as in (12) that

[TABLE]

where $\beta(p):=2(\alpha+p)(1-\sigma)/((\alpha+p)(1-\sigma)-\alpha)\geq 2$ (the exact value of $\beta(p)$ is not important in the following). For the last step, we used the crucial inequality $(\alpha+p)\sigma/p<1$ (which is equivalent to $\alpha<p/d$).

Because of $\rho^{k}=\exp(w^{k}/\delta)$, a computation shows that the integral with factor $\eta$ can be written as

[TABLE]

The last integral is bounded from below, independently of $\delta$. Since $p=2$ or $p=3$, the first integral on the right-hand side is nonnegative. (At this point, we need the term $|\nabla w^{k}|^{2}\nabla w^{k}\cdot\nabla\phi$.) Summarizing these estimates, we infer that

[TABLE]

Step 3: Limit $(\eta,\tau)\to 0$. Let $w_{\tau}(x,t)=w^{k}(x)$, $\rho_{\tau}(x,t)=\rho(w^{k}(x))$, $v_{\tau}(x,t)=v^{k}(x)$ for $x\in\Omega$ and $t\in((k-1)\tau,k\tau]$, $k=1,\ldots,N$, be piecewise constant functions in time. At time $t=0$, we set $w_{\tau}(x,0)=\log\rho^{0}(x)$ and $\rho_{\tau}(x,0)=\rho^{0}(x)$ for $x\in\Omega$. (Here, we need $\rho^{0}\geq C>0$ in $\Omega$ and another approximation procedure which we omit; see, e.g., [22, proof of Theorem 4.1].) Furthermore, we introduce the shift operator $\pi_{\tau}\rho_{\tau}(x,t)=\rho_{\tau}(x,t-\tau)$ for $x\in\Omega$, $t\geq\tau$. Then the weak formulation (9)–(10) can be written as

[TABLE]

where $\phi:(0,T)\to H^{2}(\Omega)$ and $\theta:(0,T)\to H^{1}(\Omega)$ are piecewise constant functions.

Multiplying (13) by $\tau$, summing over $k=1,\ldots,N$, and applying the discrete Gronwall inequality [22, Lemma A.2] provides the following uniform estimates:

[TABLE]

The (simultaneous) limit $(\eta,\tau)\to 0$ does not require estimates uniform in $\delta$. Therefore, we can exploit the bound for $\rho_{\tau}$ in $L^{2}(0,T;H^{1}(\Omega))$. Together with the uniform $L^{\infty}(0,T;L^{1}(\Omega))$ bound, we obtain from the Gagliardo–Nirenberg inequality as in [22, page 95] that $(\rho_{\tau})$ is bounded in $L^{2+2/d}(Q_{T})$, recalling that $Q_{T}=\Omega\times(0,T)$. Since $-\Delta v_{\tau}+v_{\tau}=\delta\rho_{\tau}+\rho_{\tau}^{\alpha}$ is bounded in $L^{2+2/d}(Q_{T})$, we deduce from elliptic regularity a uniform bound for $v_{\tau}$ in $L^{2+2/d}(0,T;W^{2,2+2/d}(\Omega))$. Therefore, $\rho_{\tau}\nabla v_{\tau}$ is uniformly bounded in $L^{1+1/d}(Q_{T})$ and $\rho_{\tau}\nabla\rho_{\tau}$ is uniformly bounded in $L^{(2d+2)/(2d+1)}(Q_{T})$. Consequently, $(\rho_{\tau}-\pi_{\tau}\rho_{\tau})/\tau=\operatorname{div}((1+\delta\rho_{\tau})\nabla\rho_{\tau}-\rho_{\tau}\nabla v_{\tau})$ is uniformly bounded in $L^{(2d+2)/(2d+1)}(0,T;W^{-1,(2d+2)/(2d+1)}(\Omega))$.

The Aubin–Lions lemma in the version of [13] shows that there exists a subsequence, which is not relabeled, such that, as $(\eta,\tau)\to 0$,

[TABLE]

for any $p<6$ and in $L^{q}(Q_{T})$ for any $q<2+2/d$. Moreover, because of the bounds (18), again for a subsequence, as $(\eta,\tau)\to 0$,

[TABLE]

We deduce that $\rho_{\tau}\nabla\rho_{\tau}\rightharpoonup\rho\nabla\rho$ and $\rho_{\tau}\nabla v_{\tau}\rightharpoonup\rho\nabla v$ weakly in $L^{1}(Q_{T})$ as well as $\rho_{\tau}^{\alpha}\to\rho^{\alpha}$ strongly in $L^{2}(Q_{T})$.

The limit in the term involving $\eta$ is performed as in [18]: Estimates (19) imply that, for any $\phi\in L^{4}(0,T;H^{2}(\Omega))$,

[TABLE]

Thus, performing the limit $(\eta,\tau)\to 0$ in (16)–(17), it follows that $(\rho,c)$ solve

[TABLE]

where, by density, we can choose test functions $\phi\in L^{\infty}(0,T;H^{1}(\Omega))$ and $\theta\in L^{2}(0,T;$ $H^{1}(\Omega))$. The initial datum $\rho(0)=\rho^{0}$ is satisfied in the sense of $H^{1}(\Omega)^{\prime}$; see, e.g., [21, pp. 1980f.] for a proof.

Step 4: Limit $\delta\to 0$. For this limit, we need further uniform estimates. Let $(\rho_{\tau},v_{\tau})$ be a solution to (16)–(17). We formulate (15) as

[TABLE]

where we recall that $C>0$ is independent of $(\eta,\tau,\delta)$. The $L^{\infty}(0,T;L^{2}(\Omega))$ and $L^{2}(0,T;$ $H^{1}(\Omega))$ bounds for $(\rho_{\tau}^{1/2})$ show that

[TABLE]

Thus, $(\rho_{\tau})$ is bounded in $L^{2}(0,T;W^{1,1}(\Omega))\hookrightarrow L^{2}(0,T;L^{3/2}(\Omega))$, as $d\leq 3$.

Let $p=2$ in (22). As the right-hand side of (22) is uniformly bounded, we infer the bounds

[TABLE]

Choosing $p=3$ in (22), the right-hand side is again bounded, yielding the estimates

[TABLE]

By elliptic regularity, since $-\Delta v_{\tau}+v_{\tau}=\delta\rho_{\tau}+\rho_{\tau}^{\alpha}\in L^{\infty}(0,T;L^{3}(\Omega))$, the family $(v_{\tau})$ is bounded in $L^{\infty}(0,T;W^{2,3}(\Omega))\hookrightarrow L^{\infty}(0,T;W^{1,p}(\Omega))$ for all $p<\infty$.

We know from Step 3 that $(\rho_{\tau},v_{\tau})$ converges in some norms to $(\rho_{\delta},v_{\delta}):=(\rho,v)$ solving (20)–(21). By the weakly lower semicontinuity of the norm and the a.e. convergence $\rho_{\tau}^{2}\to\rho^{2}$ in $Q_{T}$, it follows that, after performing the limit $(\eta,\tau)\to 0$,

[TABLE]

We wish to derive a uniform estimate for the time derivative $\partial_{t}\rho_{\delta}$. Let $\phi\in L^{2}(0,T;H^{1}(\Omega))$. Then

[TABLE]

This shows that $(\partial_{t}\rho_{\delta})$ is bounded in $L^{2}(0,T;H^{1}(\Omega)^{\prime})$. By the Aubin–Lions lemma in the version of [30], there exists a subsequence, which is not relabeled, such that, as $\delta\to 0$,

[TABLE]

Furthermore, we deduce from the bounds (23), again for a subsequence, that

[TABLE]

In particular, $\rho_{\delta}\nabla v_{\delta}\rightharpoonup\rho\nabla v$ weakly in $L^{2}(Q_{T})$. Thus, we can perform the limit $\delta\to 0$ in (20)–(21), which gives

[TABLE]

for all $\phi$, $\theta\in L^{2}(0,T;H^{1}(\Omega))$. Furthermore, we show as in [21, pp. 1980f.] that $\rho(0)=\rho^{0}$ in the sense of $H^{1}(\Omega)^{\prime}$.

Step 5: Convergence of the whole sequence. The whole sequence $(\rho_{\delta},v_{\delta})$ converges if the limit problem has a unique solution. Uniqueness follows by standard estimates if $\rho$, $\nabla c\in L^{\infty}(0,T;L^{\infty}(\Omega))$. Since $-\Delta c+c=\rho^{\alpha}\in L^{\infty}(0,T;L^{3/\alpha}(\Omega))$ and $3/\alpha>3\geq d$, elliptic regularity shows that $c\in L^{\infty}(0,T;W^{2,3/\alpha}(\Omega))\hookrightarrow L^{\infty}(0,T;W^{1,\infty}(\Omega))$. Then [17, Lemma 1] shows that $\rho\in L^{\infty}(0,T;L^{\infty}(\Omega))$, finishing the proof.

3. Proof of Theorem 2

It is shown in [18] that (1)–(2) with $\alpha=1$ has a global weak solution in two space dimensions. Since the solutions to the limiting Keller–Segel system may blow up after finite time, we cannot generally expect estimates that are uniform in $\delta$ globally in time. Moreover, we need higher-order estimates not provided by the results of [18]. Therefore, we show first the local existence of smooth solutions and then uniform $H^{s}(\Omega)$ bounds.

Step 1: Local existence of smooth solutions. Let $\varepsilon=1$. The eigenvalues of the diffusion matrix associated to (1),

[TABLE]

equal $\lambda=1\pm\mathrm{i}\sqrt{\delta\rho}$, and they have a positive real part for all $\rho>0$, i.e., $A(\rho,c)$ is normally elliptic. Therefore, according to [1, Theorem 14.1] (also see [23, Theorem 3.1]), there exists a unique maximal solution to (1)–(2) satisfying $(\rho,c)\in C^{\infty}(\overline{\Omega}\times(0,T^{*});{\mathbb{R}}^{2})$, where $0<T^{*}\leq\infty$.

Next, let $\varepsilon=0$. We use the Schauder fixed-point theorem to prove the regularity of the solutions to (4)–(5). We only sketch the proof, since the arguments are rather standard. We introduce the set

[TABLE]

for some $R>0$ and $M>0$. Let $\widetilde{\rho}\in S$. By elliptic regularity (combining Theorems 2.4.2.7 and 2.5.1.1 in [15]), the unique solution to

[TABLE]

satisfies $v\in C^{0}([0,T];W^{2,p}(\Omega))$ for all $p<\infty$. Hence, by Sobolev embedding, $h:=\widetilde{\rho}\nabla v\in C^{0}(\overline{\Omega}\times[0,T])$. Thus, using [27, Lemma 2.1iv], the unique solution to

[TABLE]

satisfies $\rho\in C^{\gamma,\gamma/2}(\overline{\Omega}\times[0,T])$. By elliptic regularity again, $v\in C^{2,\gamma/2}(\overline{\Omega}\times[0,T])$. Consequently, $h\in C^{\gamma,\gamma/2}(\overline{\Omega}\times[0,T])$ and applying [27, Lemma 2.1iv] again, we infer that $\rho\in C^{2,1}(\overline{\Omega}\times[0,T])$. It is possible to show that $\rho\in S$ for suitable $R>0$ and $M>0$. Hence, the existence of a solution to (4)–(5) follows from the Schauder fixed-point theorem.

Elliptic regularity implies that $v\in C^{4,1}(\overline{\Omega}\times[0,T])$. Then $f:=\operatorname{div}(\rho\nabla v)\in C^{1,1}(\overline{\Omega}\times[0,T])$ and the solution $u=\rho$ to the linear parabolic equation $\partial_{t}u-\Delta u-\operatorname{div}(\rho\nabla u)=f$ in $\Omega$, $t>0$, with no-flux boundary conditions satisfies $u\in C^{2+\gamma,1+\gamma/2}(\overline{\Omega}\times[0,T])$ [28, Corollary 5.1.22] (here we need $\rho^{0}\in C^{2+\gamma}(\overline{\Omega})$). Thus, the regularity of $f$ improves to $f\in C^{1+\gamma,1+\gamma/2}(\overline{\Omega}\times[0,T])$. By parabolic regularity [14, Theorem 9.2, p. 137], we infer that $u\in C^{2+\gamma,2}(\overline{\Omega}\times[0,T])$. Bootstrapping this argument and using [14, Theorems 10.1–10.2, pp. 139f.], we find that $\rho=u\in C^{\infty}(\overline{\Omega}\times(0,T])$ and consequently $v\in C^{\infty}(\overline{\Omega}\times(0,T])$.

Step 2: Preparations. Let $\varepsilon=1$, let $(\rho_{\delta},c_{\delta})$ be a local smooth solution to (1)-(2) with $0<\delta<1$, and let $(\rho,c)$ be a local smooth solution to (2)–(3). Then $\rho_{R}:=\rho_{\delta}-\rho$ and $c_{R}:=c_{\delta}-c$ solve

[TABLE]

$(\rho_{R},c_{R})$ satisfies homogeneous Neumann boundary conditions and vanishing initial conditions:

[TABLE]

The aim is to prove a differential inequality for

[TABLE]

where $\|(\rho_{R},c_{R})\|^{2}_{X}=\|\rho_{R}\|_{X}^{2}+\|c_{R}\|_{X}^{2}$ for suitable norms $\|\cdot\|_{X}$.

Step 3: $H^{1}(\Omega)$ estimates. We use $\rho_{R}$ as a test function in (24):

[TABLE]

By Young’s inequality, for any $\eta>0$, we have

[TABLE]

where here and in the following, $C>0$ and $C(\eta)>0$ denote generic constants independent of $\delta$ but depending on suitable norms of $(\rho,c)$. The embedding $H^{2}(\Omega)\hookrightarrow L^{\infty}(\Omega)$ (for $d\leq 3$) gives

[TABLE]

Combining the estimates for $I_{1}$ and $I_{2}$ and choosing $\eta>0$ sufficiently small, we find that

[TABLE]

Next, we use the test function $c_{R}$ in (25):

[TABLE]

Thus, for any $\eta>0$,

[TABLE]

For the estimate of $I_{4}$, we apply the mean-value theorem to the function $s\mapsto s^{\alpha}$ (recalling that $\alpha\geq 1$):

[TABLE]

Hence, together with the embedding $H^{2}(\Omega)\hookrightarrow L^{\infty}(\Omega)$,

[TABLE]

Collecting these estimates and choosing $\eta>0$ sufficiently small, it follows that

[TABLE]

Thus, summing (26) and (27),

[TABLE]

Step 4: $H^{2}(\Omega)$ estimates. We multiply (24) by $-\Delta\rho_{R}$ and integrate by parts in the expression with the time derivative:

[TABLE]

Then, taking into account inequality (37) in the Appendix,

[TABLE]

and choosing $\eta>0$ sufficiently small, we end up with

[TABLE]

We multiply (25) by $-\Delta c_{R}$ and estimate similarly as in Step 3:

[TABLE]

Hence, for sufficiently small $\eta>0$,

[TABLE]

Adding this inequality and (29), adding $\|c_{R}\|_{L^{2}(\Omega)}^{2}$ on both sides, using (38) in the Appendix, and choosing $\delta>0$ sufficiently small to absorb the term $C\delta^{2}\|\Delta\rho_{R}\|_{L^{2}(\Omega)}^{2}$, we infer that

[TABLE]

Step 5: $H^{3}(\Omega)$ estimates. We apply the Laplacian to (24) and (25) and multiply both equations by $\Delta\rho_{R}$, $\Delta c_{R}$, respectively:

[TABLE]

We estimate

[TABLE]

Taking into account

[TABLE]

and inequalities (37) and (39) in the Appendix as well as the embedding $H^{2}(\Omega)\hookrightarrow L^{\infty}(\Omega)$, we obtain

[TABLE]

which gives

[TABLE]

This shows that, again for sufficiently small $\eta>0$,

[TABLE]

Next, we estimate $I_{9}$ and $I_{10}$:

[TABLE]

We find that

[TABLE]

By the Hölder continuity of $s\mapsto s^{\alpha-1}$, it follows that

[TABLE]

Consequently, for sufficiently small $\eta>0$,

[TABLE]

Adding the previous inequality and (32) and taking $\delta>0$ sufficiently small such that the term $C\delta^{2}\|\nabla\Delta\rho_{R}\|_{L^{2}(\Omega)}^{2}$ is absorbed by the corresponding term on the left-hand side of (32), we infer that

[TABLE]

Step 6: End of the proof for $\varepsilon=1$. We sum inequalities (28), (31), and (34):

[TABLE]

To get rid of the term $\Gamma(t)^{\alpha-1}$, we need the condition $\alpha\geq 2$. Indeed, under this condition,

[TABLE]

We also remove the term $\Gamma(t)^{2}$ by defining $\kappa:=\max\{\alpha,2\}$ and estimating

[TABLE]

We deduce from elliptic regularity that

[TABLE]

Therefore, integrating (35) over $(0,t)$ and observing that $(\rho_{R},c_{R})(0)=0$, (35) becomes

[TABLE]

Lemma 4 proves the result for $\varepsilon=1$.

Step 7: Parabolic-elliptic case $\varepsilon=0$. Since there is no time derivative of $c_{R}$ anymore, we need to change the definition of the functionals $\Gamma(t)$ and $G(t)$:

[TABLE]

The estimates are very similar to the parabolic-parabolic case with two exceptions: In (34), we have estimated the terms $\|\rho_{R}\|_{H^{2}(\Omega)}^{2}\|\nabla\Delta c_{R}\|_{L^{2}(\Omega)}^{2}$ and $\|\rho_{R}\|_{H^{2}(\Omega)}^{2}\|c_{R}\|_{H^{2}(\Omega)}^{2}$ from above by $\Gamma(t)G(t)$. In the present case, we cannot estimate $\|\nabla\Delta c_{R}\|_{L^{2}(\Omega)}^{2}$ by $G_{0}(t)$ and we need to proceed in a different way.

Estimates (30) and (33), adapted to the case $\varepsilon=0$, become

[TABLE]

The term $\delta^{2}\|\nabla\Delta\rho_{R}\|_{L^{2}(\Omega)}^{2}$ can be absorbed by the corresponding term on the left-hand side of (32). The critical term $\|\rho_{R}\|_{H^{2}(\Omega)}^{2(\alpha-1)}\|\nabla\rho_{R}\|_{L^{2}(\Omega)}^{2}$ is bounded from above by $\Gamma(t)^{\alpha}$. Thus, $\|\nabla\Delta c_{R}\|_{L^{2}(\Omega)}^{2}$ is estimated by $\Gamma_{0}(t)^{\alpha}$ and lower-order terms, and consequently, $\|\rho_{R}\|_{H^{2}(\Omega)}^{2}$ $\times\|\nabla\Delta c_{R}\|_{L^{2}(\Omega)}^{2}$ in (32) is estimated by $\Gamma_{0}(t)^{\alpha}G_{0}(t)$, together with lower-order terms. Furthermore, $\|\rho_{R}\|_{H^{2}(\Omega)}^{2}\|c_{R}\|_{H^{2}(\Omega)}^{2}$ is bounded by $\Gamma_{0}(t)^{\alpha+1}$, up to lower-order terms. More precisely, a computation shows that

[TABLE]

Observing that $\Gamma_{0}(t)\leq C(\|\rho_{R}\|_{H^{1}(\Omega)}^{2}+\|\Delta\rho_{R}\|_{L^{2}(\Omega)}^{2})$ and $\Gamma_{0}(t)^{\alpha-1}+\Gamma_{0}(t)^{\alpha}\leq\Gamma_{0}(t)+\Gamma_{0}(t)^{\alpha+1}$, choosing $\delta>0$ sufficiently small, integrating in time, and using $\Gamma_{0}(0)=0$, we arrive at

[TABLE]

An application of Lemma 4 finishes the proof.

4. Numerical experiments

We present some numerical examples for system (1)–(2) in two space dimensions and for various choices of the parameters $\alpha$ and $\delta$. Equations (1) are discretized by the implicit Euler method in time and by cubic finite elements in space. The scheme is implemented by using the finite-element library NGSolve/Netgen (http://ngsolve.org). The mesh is refined in regions where large gradients are expected. The number of vertices is between 2805 and 12,448, and the number of elements is between 5500 and 24,030. The time step is chosen between $10^{-3}$ and $10^{-4}$ when no blow up is expected and is decreased down to $10^{-13}$ close to expected blow-up times. The resulting nonlinear discrete systems are solved by the standard Newton method. The Jacobi matrix is computed by the NGSolve routine AssembleLinearization. The surface plots are generated by the Python package Matplotlib [20]. We do not use any kind of additional regularizations, smoothing tools, or slope-limiters. All experiments are performed for the parabolic-parabolic equations with $\delta>0$.

We choose the same domain and initial conditions as in [10], i.e. $\Omega=\{x\in{\mathbb{R}}^{2}:|x|<1\}$ and

[TABLE]

A computation shows that the total mass $M=\int_{\Omega}\rho^{0}dx=25\pi/3>8\pi$ is supercritical, i.e., the solution to the classical Keller–Segel system can blow up in the interior of the domain. A sufficient condition is that the initial density is suffíciently concentrated in the sense that $\int_{\Omega}|x-x_{0}|^{2}\rho^{0}dx$ is sufficiently small for some $x_{0}\in\Omega$. Blow up at the boundary can occur if $x_{0}\in\partial\Omega$ and $M>4\pi$.

Experiment 1: $\alpha=1$. We choose the initial datum (36) and the values $\alpha=1$, $\delta=10^{-3}$. In this nonsymmetric setting, the solution exists for all time and the density is expected to concentrate at the boundary [18]. Figure 1 shows the surface plots for the cell density at various times. Since the total mass is initially concentrated near the boundary, we observe a boundary peak. Observe that there is no $L^{\infty}$ blow-up. The steady state is reached at approximately $T=2.5$. By Theorem 2, the peak approximates the blow-up solution to the classical Keller–Segel system in the $L^{\infty}$ norm; see Figure 2. We see that the $L^{\infty}$ norm of the density becomes larger with decreasing values of $\delta$.

Experiment 2: $\alpha>1$. First, we choose the value $\alpha=1.5$. The initial datum is still given by (36). Since $\alpha>1$, we cannot exclude finite-time blow-up, which is confirmed by the numerical experiments in Figure 3. Numerically, the solution seems to exist until time $T^{*}\approx 0.079$. The numerical scheme breaks down at slightly smaller times when $\delta$ becomes smaller. This may indicate that the numerical break-down is an upper bound for the blow-up time of the classical Keller–Segel model. The break-down time becomes smaller for larger values of $\alpha$. Indeed, Figure 4 shows a stronger and faster concentration behavior when we take $\alpha=2.5$.

Experiment 3: Multi-bump initial datum. We take $\alpha=1$ and $\delta=5\cdot 10^{-3}$. As initial datum, we choose a linear combination of the bump function

[TABLE]

where $(x_{0},y_{0})\in\Omega$, $M>0$, and $\theta>0$. Setting $\theta=10^{-2}$, we define $c^{0}=0$ and

[TABLE]

The evolution of the density is presented in Figure 5. The density concentrates in the interior of the domain and the peak travels to the boundary. At time $t=1$, the peak is close to the boundary which is reached later at $t=2.5$ (not shown). A similar behavior was already mentioned in [4] for the parabolic-elliptic model using a single-bump initial datum.

Experiment 4: Shape of peaks. The previous experiments show that the shape of the peaks depends on the value of $\delta$. In this experiment, we explore this dependence in more detail. We claim that the diameter and the height of the bump can be controlled by $\delta$. We choose $\alpha=1$ and the initial datum $\rho^{0}=W_{0,0,20\pi}$ with $\theta=1/400$ and $c^{0}=0$. Furthermore, we prescribe homogeneous Dirichlet boundary conditions for $c$ to avoid that the aggregated bump of cells moves to the boundary. Figure 6 (top row) shows the stationary cell densities for two values of $\delta$. As expected, the maximal diameter of the peak (defined at height $10^{-2}$) becomes smaller and the maximum of the peak becomes larger for decreasing values of $\delta$. The level sets show that the solutions are almost radially symmetric and the level set for $\rho=10^{-2}$ is approximately a circle. This behavior is quantified in Figure 6 (bottom row). We observe that the radius depends on $\delta$ approximately as $r\sim\delta^{0.43}$ and the height approximately as $\rho_{\rm max}\sim\delta^{-1.00}$.

We remark that under no-flux boundary conditions for the chemical concentration, the same behavior of the bumps can be observed for intermediate times. However, the bump will eventually move to the boundary (as in Figure 5), since the chemical substance is not absorbed by the boundary as in the Dirichlet case.

Appendix A Some technical tools

For the convenience of the reader, we collect some technical results.

Lemma 3** (Inequalities).**

Let $d\leq 3$, $\Omega\subset{\mathbb{R}}^{d}$ be a bounded domain, and $\partial\Omega\in C^{2,1}$. There exists a constant $C>0$ such that for all $u$, $v\in H^{1}(\Omega)$,

[TABLE]

for all $u\in H^{2}(\Omega)$ with $\nabla u\cdot\nu=0$ on $\partial\Omega$,

[TABLE]

and for all $u\in H^{3}(\Omega)$ with $\nabla u\cdot\nu=0$ on $\partial\Omega$,

[TABLE]

Inequality (37) follows after applying the Cauchy–Schwarz inequality and then the continuous embedding $H^{1}(\Omega)\hookrightarrow L^{4}(\Omega)$; (38) is proved in [15, Theorem 2.3.3.6], while (39) is a consequence of [31, Theorem 2.24].

Lemma 4** (Nonlinear Gronwall inequality).**

Let $\delta>0$ and $\Gamma$, $G\in C^{0}([0,T])$ be nonnegative functions, possibly depending on $\delta$, satisfying

[TABLE]

where $\alpha>1$, $\beta\geq 0$, $\gamma>0$, $\nu>0$, and $C_{0},\ldots,C_{4}>0$ are constants independent of $\delta$. Furthermore, let $\Gamma(0)\leq C_{5}\delta^{\nu}$ for some $C_{5}>0$. Then there exists $\delta_{0}>0$ such that for all $0<\delta<\delta_{0}$, $0\leq t\leq T$, and $0<\varepsilon<\nu$,

[TABLE]

Proof.

A slightly simpler variant of the lemma was proved in [19, Lemma 10]. Assume, by contradiction, that for all $\delta_{0}\in(0,1)$, there exist $\delta\in(0,\delta_{0})$, $t_{0}\in[0,T]$, and $\varepsilon\in(0,\nu)$ such that $\Gamma(t_{0})>C_{5}\delta^{\nu-\varepsilon}$. Since $\Gamma(0)\leq C_{5}\delta^{\nu}$ by assumption and $\Gamma$ is continuous, there exists $t_{1}\in[0,t_{0})$ such that $\Gamma(t_{1})=C_{5}\delta^{\nu-\varepsilon}$ and $\Gamma(t)\leq C_{5}\delta^{\nu-\varepsilon}$ for all $t\in[0,t_{1}]$. This leads for $t\in[0,t_{1}]$ to

[TABLE]

Since $\nu-\varepsilon>0$, the integral over $G(s)$ on the right-hand side can be absorbed for sufficiently small $\delta>0$ by the corresponding term on the left-hand side. This implies that

[TABLE]

Then Gronwall’s lemma gives, for sufficiently small $\delta_{0}>0$ and $0<\delta<\delta_{0}$,

[TABLE]

which contradicts $\Gamma(t_{1})=C_{5}\delta^{\nu-\varepsilon}$. ∎