Vanishing cross-diffusion limit in a Keller-Segel system with additional cross-diffusion
Ansgar Jüngel Oliver Leingang Shu Wang
math.AP 35B40, 35K51, 35Q92, 92C17
  

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

A. Jüngel]Ansgar Jüngel
Institute for Analysis and Scientific Computing, Vienna University of  
	Technology, Wiedner Hauptstraße 8–10, 1040 Wien, Austria
juengel@tuwien.ac.at 

O. Leingang]Oliver Leingang
Institute for Analysis and Scientific Computing, Vienna University of  
	Technology, Wiedner Hauptstraße 8–10, 1040 Wien, Austria
oliver.leingang@tuwien.ac.at 

S. Wang]Shu Wang
College of Applied Science, Beijing University of Technology, 
Beijing, PR China
wangshu@bjut.edu.cn 



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. 


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(Ω) 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.



  
 

[2000]35B40, 35K51, 35Q92, 92C17.

[
    [
    December 30, 2023
=====================







§ INTRODUCTION


Chemotaxis describes the directed movement of cells in response to chemical gradients
and may be modeled by the (Patlak–) Keller–Segel equations <cit.>. 
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 <cit.>, degenerate cell diffusion 
<cit.>, death-growth terms <cit.>,
or additional cross-diffusion <cit.>. In this paper, we continue
our study of the Keller–Segel system with additional cross-diffusion, which allows
for global weak solutions <cit.>.
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 ρ_δ(x,t) and the chemical concentration
c_δ(x,t) are given by

    _tρ_δ = (ρ_δ-ρ_δ c_δ),   _t c_δ = Δ c_δ + δΔρ_δ 
    	- c_δ + ρ_δ^α  Ω, t>0,

subject to the no-flux and initial conditions

    ρ_δ·ν =  c_δ·ν = 0  Ω, t>0,   ρ_δ(0)=ρ^0,  c_δ(0)= c^0  Ω,

where Ω⊂^d (d=2,3) is a bounded domain,
ν is the exterior unit normal vector of Ω,
δ>0 describes the strength of the additional cross-diffusion, and 
the term ρ_δ^α with α>0 is the nonlinear signal production.
The case =1 is called the parabolic-parabolic case, while =0
refers to the parabolic-elliptic case. 
The existence of global weak solutions to (<ref>)–(<ref>)
was proved in <cit.> in two space dimensions and in <cit.> 
in three space dimensions (with degenerate diffusion). 

Setting δ=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, α=1. In order to deal
with the three-dimensional case, we need sublinear signal productions, 
α<1.
In fact, it is shown in <cit.> that α=2/d is the critical
value for global existence versus finite-time blow-up 
in a slightly modified Keller–Segel system with 
δ=0. We show that also for δ>0, the condition α<2/d
guarantees the global existence of weak solutions, while numerical results
indicate finite-time blow-up if α>2/d.

We are interested in the limit δ→ 0 in (<ref>) leading to the 
Keller–Segel equations

    _tρ = (ρ-ρ c),   _t c = Δ c - c + ρ^α  Ω, t>0,

with the initial and boundary conditions (<ref>). 
We refer to the reviews <cit.> 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 (<ref>) 
in two space dimensions: If the initial mass M:=ρ^0_L^1(Ω) is smaller 
than 8π, the solutions are global in time, while they blow up in finite time
if M>8π and the second moment of the initial datum is finite (see, e.g.,
<cit.>). 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π such that
the corresponding solutions are global <cit.>.

In the parabolic-parabolic case, again in two space dimensions
and with finite second moment, the solutions exist globally in time of 
M<8π <cit.>. 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≥ 3, a related critical phenomenon occurs: The solutions
exist globally in time if ρ^0_L^d/2(Ω) is sufficiently small,
but they blow up in finite time if the total mass is large compared to the
second moment <cit.>.

In this paper, we prove two results. The first one is the convergence of the solutions
to the parabolic-elliptic model (<ref>) with =0 to a solution
to the parabolic-elliptic Keller–Segel system (<ref>) with =0.
Since we impose a restriction on the parameter α as mentioned above,
this result holds globally in time. We set Q_T=Ω×(0,T).


Let Ω⊂^d (d=2,3) be bounded with Ω∈ C^1,1,
T>0, δ>0, =0, and 0≤ρ^0∈ L^∞(Ω). 
If α<2/d, there exists
a weak solution (ρ_δ,c_δ)∈ L^2(0,T;H^1(Ω))^2 to 
(<ref>)-(<ref>) satisfying

    _tρ_δ∈ L^2(0,T;H^1(Ω)'),   ρ_δ∈ L^∞(0,T;L^3(Ω)),   
    	c_δ+δρ_δ∈ L^∞(0,T;W^1,p(Ω)),

for any p<∞. Furthermore, as δ→ 0,

    ρ_δ→ρ     L^2(Q_T), 
    ρ_δ⇀ρ     L^2(Q_T), 
    
    	c_δ+δρ_δ⇀^* c      
    	L^∞(0,T;W^1,p(Ω)),  p<∞,

and (ρ,c) solves (<ref>)–(<ref>), and it holds that
ρ, c∈ L^∞(0,T;L^∞(Ω)).


The idea of the proof is to reformulate (<ref>) via 
v_δ:=c_δ+δρ_δ as the system

    _tρ_δ = ((1+δρ_δ)ρ_δ
    	-ρ_δ v_δ),   
    	-Δ v_δ+v_δ = δρ_δ + ρ_δ^α  Ω, t>0,

together with the initial and boundary conditions

    ρ_δ·ν =  v_δ·ν = 0  Ω,   ρ_δ(0)=ρ^0  Ω.

This reformulation was already used in <cit.> to prove the existence of
weak solutions in the two-dimension case with α=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
δ to apply compactness arguments. For this, we use the “entropy”
functional H_1(ρ_δ)=∫_Ωρ_δ(logρ_δ-1)dx:

    dH_1/dt + 4∫_Ω|ρ_δ^1/2|^2 dx
    	+ δ∫_Ω|ρ_δ|^2 dx 
    	≤∫_Ω(δρ_δ^2 + ρ_δ^α+1)dx.

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

    δ∫_Ωρ_δ^2 dx 
    	≤δ/2ρ_δ_L^2(Ω)^2
    	+ Cρ_δ_L^1(Ω)^2.

The first term on the right-hand side is absorbed by the left-hand side of
(<ref>) and the second term is bounded since the total mass is conserved.
To estimate the second term on the right-hand side of (<ref>), we need
another “entropy” functional H_p(ρ)=∫_Ωρ_δ^p dx with
p=2 or p=3, leading to

    dH_p/dt + ∫_Ω|ρ_δ^p/2|^2 dx
    	+ δ∫_Ω|ρ_δ^(p+1)/2|^2 dx
    	≤ Cρ_δ_L^(p+1)/2(Ω)^(p+1)/2
    	+ Cρ_δ^p/2_L^1(Ω)^β(p),

where β(p)≥ 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 
ρ_δ^1/2 in H^1(Ω) from (<ref>)
implies that ρ_δ is bounded
in W^1,1(Ω)↪ L^3/2(Ω) such that the first term on the
right-hand side of (<ref>) 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ρ_δ in (<ref>) which may be
not defined if ρ_δ=0. 
Therefore, we consider an implicit Euler discretization in time
with parameter τ>0 and an elliptic regularization in space with parameter
η>0 to prove first the existence of approximate
weak solutions with strictly positive ρ_δ. This is done by using
the entropy method of <cit.>. The approximate solutions also satisfy
the δ-uniform bounds, and they hold true when passing to the
limit (η,τ)→ 0. Then the limit δ→ 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 α≥ 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.


Let Ω⊂^d (d≤ 3) be a bounded domain with smooth boundary
and let (ρ^0,c^0)∈ (W^1,p(Ω))^2 for p>d if =1 and
ρ^0∈ C^2+γ(Ω) for some γ∈(0,1) if =0. 
Furthermore,
let α=1 or α≥ 2 and let (ρ_δ,c_δ) and (ρ,c)
be (weak) solutions to (<ref>) and (<ref>), (<ref>), respectively,
with the same initial data.
Then these solutions are smooth locally in time and
there exist constants C>0 and δ_0>0
such that for all 0<δ<δ_0 and λ>0,

    (ρ_δ-ρ,c_δ-c)_L^∞(0,T;H^2(Ω))≤ Cδ^1-λ.



The theorem is proved as in, e.g., <cit.> by deriving carefully 
H^s(Ω) estimates for the difference
(ρ_R,c_R):=(ρ_δ-ρ,c_δ-c).
The index s∈ is chosen such that we obtain L^∞(Ω)
estimates which are needed to handle the nonlinearities. 
If =1, we introduce the functions

    Γ(t) = (ρ_R,c_R)(t)_H^2(Ω)^2,   
    	G(t) = (ρ_R,c_R)(t)_H^2(Ω)^2 
    	+ Δ(ρ_R,c_R)(t)_L^2(Ω)^2.

The aim is to prove the inequality

    Γ(t) + C_1∫_0^t G(s)ds 
    	≤ C_2∫_0^t(Γ(s)+Γ(s)^max{2,α})ds
    	+ C_2∫_0^t Γ(s)G(s)ds + C_2δ^2,

where C_1, C_2>0 are constants independent of δ. 
This inequality allows us to apply
a variant of Gronwall's lemma (see Lemma <ref> in the Appendix), 
implying (<ref>). 
In the parabolic-elliptic case =0, the functions Γ(t) and G(t)
are defined without c_R, and the final inequality contains the
additional integral ∫_0^tΓ(s)^α+1G(s)ds, which is still
covered by Lemma <ref>. The condition α=1 or α≥ 2 comes
from the fact that the derivative of the mapping s↦ s^α is 
Hölder continuous
exactly for these values. The numerical results in Section <ref> indicate
that the convergence result may still hold for α∈(1,2). 

The optimal convergence rate is expected to be one. Theorem <ref> provides 
an almost optimal rate. The reason for the non-optimality comes from the variant
of the nonlinear Gronwall lemma proved in Lemma <ref>. We conjecture that
an optimal rate holds (changing the constants in Lemma <ref>), 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
(<ref>) were already performed in <cit.>.
For positive values of δ, the (globally existing)
cell density forms bumps at places
where the solution to the classical Keller–Segel system develops an 
L^∞(Ω) blow up.
Compared to <cit.>, we investigate the dependence of the shape
of the bumps on δ. In a radially symmetric situation,
it turns out that the radius of the bump (more precisely the diameter of
a level set ρ_δ≈ 0) behaves like δ^a with a≈ 0.43,
and the maximum of the bump behaves like δ^-b with 
b≈ 1.00. 

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





§ PROOF OF THEOREM <REF>


We prove Theorem <ref> by an approximation procedure and by deriving the
uniform bounds from discrete versions of the entropy inequalities
(<ref>) and (<ref>).

Step 1: Solution of a regularized system and entropy estimates.
We show the existence of solutions to a discretized and regularized
version of (<ref>). For this, let N∈ and τ=T/N, and set
ρ(w)=exp(w/δ). This means that we transform w=δlogρ.
Let w^k-1∈ H^2(Ω) and v^k-1∈ H^1(Ω) be given
and set ρ^j=ρ(w^j) for j=k,k-1. Consider for given τ>0
and η>0 the regularized system

    1/τ∫_Ω(ρ^k-ρ^k-1)ϕ dx
      + ∫_Ω((1+δρ^k)ρ^k - ρ^k v^k)·ϕ dx 
    	
       xxxx= -η∫_Ω(Δ w^k Δϕ
    	+ δ^-2| w^k|^2 w^k·ϕ + w^kρ^k ϕ)dx, 
       ∫_Ω( v^k·θ + v^kθ)dx 
    	= ∫_Ω(δρ^k + (ρ^k)^α)θ dx

for ϕ∈ H^2(Ω) and θ∈ H^1(Ω). 
The time discretization is needed to handle issues due to low time regularity,
while the elliptic regularization guarantees H^2(Ω) solutions which, by
Sobolev embedding (recall that d≤ 3), are bounded. The higher-order gradient
term | w^k|^2 w^k·ϕ is necessary to derive 
L^p(Ω) estimates. The existence of a solution
w^k∈ H^2(Ω), 0≤ v^k∈ H^1(Ω) follows from the techniques
used in the proof of Proposition 3.1 in <cit.> 
employing the Leray–Schauder fixed-point theorem. Since these techniques are
rather standard now, we omit the proof and refer to <cit.>
for similar arguments.

Inequality w^kρ^k=w^k e^w^k/δ≥ e^w^k/δ-1=ρ^k-1 allows
us to show as in <cit.> that the total mass 
ρ^k_L^1(Ω) is bounded uniformly in δ.

Entropy estimates are derived from (<ref>)–(<ref>) 
by choosing the test functions ϕ=w^k/δ=logρ^k and θ=ρ^k, 
respectively, and adding both equations.
Then the terms involving v^k cancel and after some elementary computations,
we end up with

    1/τ∫_Ω   (h(ρ^k)-h(ρ^k-1))dx
    	+ 4∫_Ω|(ρ^k)^1/2|^2 dx + δ∫_Ω|ρ^k|^2 dx 
    	
       xx+ η/δ∫_Ω((Δ w^k)^2
    	+ δ^-2| w^k|^4 + (w^k)^2ρ^k)dx 
       ≤∫_Ω(-v^kρ^k + δ(ρ^k)^2 + (ρ^k)^α+1) dx
    	≤∫_Ω(δ(ρ^k)^2 + (ρ^k)^α+1) dx,

where h(s)=s(log s-1) for s≥ 0. Using the Gagliardo–Nirenberg inequality
with σ=d/(d+2), the Poincaré–Wirtinger inequality, and then the
Young inequality with p=1/σ, p'=1/(1-σ), it holds for any
u∈ H^1(Ω) and μ>0 that

    u_L^2(Ω)^2
    	   ≤ Cu_H^1(Ω)^2σu_L^1(Ω)^2(1-σ)
    	= C( u_L^2(Ω)^2 + u_L^1(Ω)^2)^σu_L^1(Ω)^2(1-σ)
       ≤μ( u_L^2(Ω)^2 + u_L^1(Ω)^2) 
    	+ C(μ)u_L^1(Ω)^2
    	= μ u_L^2(Ω)^2 + C(μ)u_L^1(Ω)^2.

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

    δ∫_Ω(ρ^k)^2 dx
    	≤δ/4ρ^k_L^2(Ω)^2 + Cρ^k_L^1(Ω)^2,

where here and in the following, C>0 denotes a constant, independent of 
η, τ, and δ, with values varying from line to line.
The first term on the right-hand side can be absorbed by the left-hand side 
of (<ref>), while the second term is bounded.
Using s^α+1≤ s^2 +1 for s≥ 0 (since α<1) and (<ref>),
the second term on the right-hand side of (<ref>) becomes

    ∫_Ω(ρ^k)^α+1dx ≤∫_Ω((ρ^k)^2+1) dx 
    	≤δ/4ρ^k_L^2(Ω)^2 
    	+ Cρ^k_L^1(Ω)^2 + C(Ω).

Inserting these estimations into (<ref>), we conclude that

    1/τ∫_Ω   (h(ρ^k)-h(ρ^k-1))dx
    	+ 4∫_Ω|(ρ^k)^1/2|^2 dx 
    	+ δ/2∫_Ω|ρ^k|^2 dx 
    	
       xx+ η/δ∫_Ω((Δ w^k)^2
    	+ δ^-2| w^k|^4 + (w^k)^2ρ^k)dx
    	≤ C(Ω).


Step 2: Further uniform estimates.
The estimates from (<ref>) are not sufficient for the limit δ→ 0,
therefore, we derive further uniform bounds.
Let p=2 or p=3. We choose the admissible test functions 
p(ρ^k)^p-1 and (p-1)(ρ^k)^p
in (<ref>)–(<ref>), respectively, 
and add both equations. The convexity of s↦ s^p implies that
s^p-t^p≤ p(s-t)s^p-1. Then, observing that the terms involving v^k
cancel, we find that

    1/τ∫_Ω   ((ρ^k)^p-(ρ^k-1)^p) dx
    	+ 4/p(p-1)∫_Ω|(ρ^k)^p/2|^2 dx
    	+ δ4p(p-1)/(p+1)^2∫_Ω|(ρ^k)^(p+1)/2|^2dx 
       xx
    	+ η p∫_Ω(Δ w^kΔ (ρ^k)^p-1
    	+ δ^-2| w^k|^2 w^k·(ρ^k)^p-1 + w^k(ρ^k)^p)dx 
    	
       ≤ (p-1)∫_Ω(δ(ρ^k)^p+1 + (ρ^k)^p+α)dx.

By (<ref>), we find that

    δ∫_Ω(ρ^k)^p+1dx = δ(ρ^k)^(p+1)/2_L^2(Ω)^2
    	≤δ/2(ρ^k)^(p+1)/2_L^2(Ω)^2
    	+ C(ρ^k)^(p+1)/2_L^1(Ω)^2.

The estimate for (ρ^k)^α+p requires that α<2/d. Indeed,
we deduce from the Gagliardo–Nirenberg inequality with 
σ=d(2α+p)/((d+2)(α+p)) and similar arguments as in (<ref>) that

    ∫_Ω(ρ^k)^α+pdx 
    	   = (ρ^k)^p/2_L^2(α+p)/p(Ω)^2(α+p)/p≤ C(ρ^k)^p/2_H^1(Ω)^2(α+p)σ/p(ρ^k)^p/2_L^1(Ω)^2(α+p)(1-σ)/p
       = C((ρ^k)^p/2_L^2(Ω)^2 + (ρ^k)^p/2_L^1(Ω)^2
    	)^(α+p)σ/p(ρ^k)^p/2_L^1(Ω)^2(α+p)(1-σ)/p
       ≤1/2((ρ^k)^p/2_L^2(Ω)^2 
    	+ (ρ^k)^p/2_L^1(Ω)^2
    	) + C(ρ^k)^p/2_L^1(Ω)^β(p),

where β(p):=2(α+p)(1-σ)/((α+p)(1-σ)-α)≥ 2
(the exact value of β(p) is not important in the following).
For the last step, we used
the crucial inequality (α+p)σ/p<1 (which is equivalent to
α<p/d).

Because of ρ^k=exp(w^k/δ), a computation shows that the integral 
with factor η can be written as

    η   p∫_Ω(Δ w^kΔ (ρ^k)^p-1
    	+ δ^-2| w^k|^2 w^k·(ρ^k)^p-1 + w^k(ρ^k)^p)dx 
       = η/δ p(p-1)∫_Ω(ρ^k)^p-1[
    	(Δ w^k + p-1/2δ| w^k|^2)^2
    	+ 1/4δ^2(4-(p-1)^2)| w^k|^4]dx 
       xx+ η p∫_Ω w^k e^pw^k/δdx.

The last integral is bounded from below, independently of δ. 
Since p=2 or p=3, 
the first integral on the right-hand side is nonnegative. (At this point, we need
the term | w^k|^2 w^k·ϕ.) Summarizing these estimates, we infer that

    1/τ∫_Ω   ((ρ^k)^p-(ρ^k-1)^p) dx
    	+ ∫_Ω|(ρ^k)^p/2|^2 dx
    	+ δ/4∫_Ω|(ρ^k)^(p+1)/2|^2dx 
       ≤ C(ρ^k)^(p+1)/2_L^1(Ω)^2 
    	+ C(ρ^k)^p/2_L^1(Ω)^β(p).


Step 3: Limit (η,τ)→ 0. 
Let w_τ(x,t)=w^k(x), ρ_τ(x,t)=ρ(w^k(x)), v_τ(x,t)=v^k(x) 
for x∈Ω and t∈((k-1)τ,kτ], k=1,…,N, 
be piecewise constant functions in time. 
At time t=0, we set w_τ(x,0)=logρ^0(x) and ρ_τ(x,0)=ρ^0(x)
for x∈Ω. (Here, we need ρ^0≥ C>0 in Ω and another approximation
procedure which we omit; see, e.g., <cit.>.)
Furthermore, we introduce the shift operator
π_τρ_τ(x,t)=ρ_τ(x,t-τ) for x∈Ω, t≥τ.
Then the weak formulation (<ref>)–(<ref>) can be written as

    1/τ∫_0^T∫_Ω(ρ_τ-π_τρ_τ)ϕ dxdt
    	+ ∫_0^T∫_Ω((1+δρ_τ)ρ_τ-ρ_τ v_τ)
    	·ϕ dxdt 
       xxx= -η∫_0^T∫_Ω(Δ w_τΔϕ
    	+ δ^-2| w_τ|^2 w_τ·ϕ + w_τρ_τϕ)dxdt, 
    	
       ∫_0^T∫_Ω( v_τ·θ + v_τθ)dxdt
    	= ∫_0^T∫_Ω(δρ_τ + ρ_τ^α)θ dxdt,

where ϕ:(0,T)→ H^2(Ω) and θ:(0,T)→ H^1(Ω)
are piecewise constant functions.

Multiplying (<ref>) by τ, summing over k=1,…,N, and applying the
discrete Gronwall inequality <cit.> provides the following uniform
estimates:

    ρ_τ_L^∞(0,T;L^1(Ω)) 
    	+ ρ_τ^1/2_L^2(0,T;H^1(Ω)) 
    	+ δ^1/2ρ_τ_L^2(0,T;H^1(Ω))   ≤ C, 
    η^1/2Δ w_τ_L^2(Q_T)
    	+ η^1/4δ^-1/2 w_τ_L^4(Q_T)
    	+ η^1/2w_τρ_τ^1/2_L^2(Q_T)   ≤ C.


The (simultaneous) limit (η,τ)→ 0 does not require estimates 
uniform in δ.
Therefore, we can exploit the bound for ρ_τ in L^2(0,T;H^1(Ω)).
Together with the uniform L^∞(0,T;L^1(Ω)) bound, we obtain from
the Gagliardo–Nirenberg inequality as in <cit.> that
(ρ_τ) is bounded in L^2+2/d(Q_T), recalling that Q_T=Ω×(0,T).
Since -Δ v_τ + v_τ = δρ_τ + ρ_τ^α is
bounded in L^2+2/d(Q_T), we deduce from elliptic regularity a uniform bound for 
v_τ in L^2+2/d(0,T;W^2,2+2/d(Ω)). Therefore,
ρ_τ v_τ is uniformly bounded in L^1+1/d(Q_T) and
ρ_τρ_τ is uniformly bounded in L^(2d+2)/(2d+1)(Q_T).
Consequently, (ρ_τ-π_τρ_τ)/τ=((1+δρ_τ)ρ_τ
- ρ_τ v_τ) is uniformly bounded in L^(2d+2)/(2d+1)
(0,T;W^-1,(2d+2)/(2d+1)(Ω)).

The Aubin–Lions lemma in the version of <cit.> shows that there exists
a subsequence, which is not relabeled, such that, as (η,τ)→ 0,

    ρ_τ→ρ  L^2(0,T;L^p(Ω))

for any p<6 and in L^q(Q_T) for any q<2+2/d. 
Moreover, because of the bounds (<ref>), 
again for a subsequence, as (η,τ)→ 0,

    ρ_τ⇀ρ     L^2(Q_T), 
    τ^-1(ρ_τ-π_τρ_τ) ⇀_tρ     L^(2d+2)/(2d+1)(0,T;W^-1,(2d+2)/(2d+1)(Ω)), 
    
      v_τ⇀ v      L^2+2/d(0,T;W^2,2+2/d(Ω)).

We deduce that ρ_τρ_τ⇀ρρ and
ρ_τ v_τ⇀ρ v weakly in L^1(Q_T) as well as
ρ_τ^α→ρ^α strongly in L^2(Q_T). 

The limit in the term involving η is performed as in <cit.>:
Estimates (<ref>) imply that, for any ϕ∈ L^4(0,T;H^2(Ω)),

    |η   ∫_0^T∫_Ω(Δ w_τΔϕ
    	+ δ^-2| w_τ|^2 w_τ·ϕ + w_τρ_τϕ)
    	dxdt| 
       ≤ηΔ w_τ_L^2(Q_T)Δϕ_L^2(Q_T)
    	+ ηδ^-2 w_τ_L^4(Q_T)^3ϕ_L^4(Q_T)
       xx+ ηw_τρ_τ^1/2_L^2(Q_T)ρ_τ^1/2_L^4(Q_T)ϕ_L^4(Q_T)
       ≤ C(δ)(η^1/2+η^1/4)ϕ_L^4(0,T;H^2(Ω))→ 0  η→ 0.


Thus, performing the limit (η,τ)→ 0 in (<ref>)–(<ref>), 
it follows that (ρ,c) solve

    ∫_0^T⟨_tρ,ϕ⟩ dt
    	+ ∫_0^T∫_Ω((1+δρ)ρ - ρ v)·ϕ dxdt = 0, 
    	
       ∫_0^T∫_Ω( v·θ + vθ)dxdt
    	= ∫_0^T∫_Ω(δρ + ρ^α)θ dxdt,

where, by density, we can choose test functions ϕ∈ 
L^∞(0,T;H^1(Ω)) and θ∈ L^2(0,T; H^1(Ω)).
The initial datum ρ(0)=ρ^0 is satisfied in the sense of H^1(Ω)'; 
see, e.g., <cit.> for a proof.

Step 4: Limit δ→ 0. For this limit, we need further uniform estimates.
Let (ρ_τ,v_τ) be a solution to (<ref>)–(<ref>).
We formulate (<ref>) as

    ∫_Ω   (ρ_τ^p-(π_τρ_τ)^p)dx
    	+ ∫_0^T∫_Ω|ρ_τ^p/2|^2 dxdt
    	+ δ/4∫_0^T∫_Ω|ρ_τ^(p+1)/2|^2 dxdt 
       ≤ C∫_0^Tρ_τ^(p+1)/2_L^1(Ω)^2 dt
    	+ C∫_0^Tρ_τ^p/2_L^1(Ω)^β(p)dt,

where we recall that C>0 is independent of (η,τ,δ). 
The L^∞(0,T;L^2(Ω)) and L^2(0,T; H^1(Ω)) bounds 
for (ρ_τ^1/2) show that

    ∫_0^Tρ_τ_L^1(Ω)^2 dt
    	   = 4∫_0^Tρ_τ^1/2ρ_τ^1/2_L^1(Ω)^2 dt
    	≤ 4∫_0^Tρ_τ^1/2_L^2(Ω)^2
    	ρ_τ^1/2_L^2(Ω)^2 dt 
       ≤ 4ρ_τ_L^∞(0,T;L^1(Ω))∫_0^Tρ_τ^1/2_L^2(Ω)^2 dt ≤ C.

Thus, (ρ_τ) is bounded in L^2(0,T;W^1,1(Ω))↪
L^2(0,T;L^3/2(Ω)), as d≤ 3.

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

    ρ_τ_L^∞(0,T;L^2(Ω)) 
      + ρ_τ_L^2(0,T;H^1(Ω)) 
    	+ δ^1/2ρ_τ^3/2_L^2(0,T;H^1(Ω))≤ C.

Choosing p=3 in (<ref>), the right-hand side is again bounded,
yielding the estimates

    ρ_τ_L^∞(0,T;L^3(Ω)) 
    	+ δ^1/2ρ_τ^2_L^2(0,T;H^1(Ω))≤ C.

By elliptic regularity, since -Δ v_τ+v_τ=δρ_τ+ρ_τ^α∈ L^∞(0,T;L^3(Ω)), the family (v_τ) is bounded in 
L^∞(0,T;W^2,3(Ω))↪ L^∞(0,T;W^1,p(Ω)) 
for all p<∞.

We know from Step 3 that (ρ_τ,v_τ) converges in some norms to
(ρ_δ,v_δ):=(ρ,v) solving (<ref>)–(<ref>).
By the weakly lower semicontinuity of the norm and the a.e. convergence
ρ_τ^2→ρ^2 in Q_T, it follows that, after performing
the limit (η,τ)→ 0,

    ρ_δ_L^∞(0,T;L^3(Ω))
    	+ ρ_δ_L^2(0,T;H^1(Ω))
    	+ δ^1/2ρ_δ^2_L^2(0,T;H^1(Ω)) 
    	+ v_δ_L^∞(0,T;W^1,p(Ω))≤ C.


We wish to derive a uniform estimate for the time derivative _tρ_δ.
Let ϕ∈ L^2(0,T;H^1(Ω)). Then

    |∫_0^T⟨_tρ_δ,ϕ⟩ dt|
    	   ≤(ρ_δ_L^2(Q_T) 
    	+ δ/2(ρ_δ^2)_L^2(Q_T))ϕ_L^2(Q_T)
       xx+ ρ_δ_L^2(0,T;L^4(Ω)) v_δ_L^∞(0,T;L^4(Ω))ϕ_L^2(Q_T)≤ C.

This shows that (_tρ_δ) is bounded in L^2(0,T;H^1(Ω)').
By the Aubin–Lions lemma in the version of <cit.>, there exists a subsequence,
which is not relabeled, such that, as δ→ 0,

    ρ_δ→ρ  L^2(0,T;L^p(Ω)),   
    	p<6.

Furthermore, we deduce from the bounds (<ref>), again for a subsequence, that

    ρ_δ⇀ρ     L^2(Q_T), 
    δ(ρ_δ^2)→ 0      L^2(Q_T), 
    _tρ_δ⇀_tρ     L^2(0,T;H^1(Ω)'),
    	
    
    	v_δ⇀^* v      L^∞(0,T;W^1,p(Ω)),
    	   p<∞.

In particular, ρ_δ v_δ⇀ρ v weakly in
L^2(Q_T). Thus, we can perform the limit δ→ 0 in 
(<ref>)–(<ref>), which gives

    ∫_0^T⟨_tρ,ϕ⟩ dt
    	+ ∫_0^T∫_Ω(ρ-ρ v)·ϕ dxdt    = 0, 
    ∫_0^T∫_Ω( v·θ + vθ)dxdt
    	   = ∫_0^T∫_Ωρ^αθ dx

for all ϕ, θ∈ L^2(0,T;H^1(Ω)). Furthermore, we show as in
<cit.> that ρ(0)=ρ^0 in the sense of H^1(Ω)'.

Step 5: Convergence of the whole sequence. The whole sequence
(ρ_δ,v_δ) converges if the limit problem has a unique solution.
Uniqueness follows by standard estimates if ρ, c∈ 
L^∞(0,T;L^∞(Ω)). Since -Δ c+c=ρ^α∈ 
L^∞(0,T;L^3/α(Ω)) and 3/α>3≥ d, elliptic regularity shows that
c∈ L^∞(0,T;W^2,3/α(Ω))↪ 
L^∞(0,T;W^1,∞(Ω)).
Then <cit.> shows that ρ∈ L^∞(0,T;L^∞(Ω)),
finishing the proof.





§ PROOF OF THEOREM <REF>


It is shown in <cit.> that (<ref>)–(<ref>) with α=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 δ globally in time. 
Moreover, we need higher-order estimates not provided by
the results of <cit.>. Therefore, we show first the local existence of
smooth solutions and then uniform H^s(Ω) bounds.

Step 1: Local existence of smooth solutions.
Let =1. The eigenvalues of the diffusion matrix associated to
(<ref>), 

    A(ρ,c) = [  1 -ρ;  δ  1 ],

equal λ=1±i√(δρ), 
and they have a positive real part for all
ρ>0, i.e., A(ρ,c) is normally elliptic.
Therefore, according to <cit.> (also see <cit.>), 
there exists a unique maximal solution to (<ref>)–(<ref>) satisfying
(ρ,c)∈ C^∞(Ω×(0,T^*);^2),
where 0<T^*≤∞.

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

    S = {ρ∈ C^0(Ω×[0,T]):
    	0≤ρ≤ R, ρ_C^γ,γ/2(Ω×[0,T])≤ K}

for some R>0 and M>0. Let ρ∈ S. By elliptic regularity
(combining Theorems 2.4.2.7 and 2.5.1.1 in <cit.>), the unique solution to

    -Δ v+v = δρ+ρ^α  Ω,
    	  ρ·ν=0  Ω,

satisfies v∈ C^0([0,T];W^2,p(Ω)) for all p<∞. Hence, by
Sobolev embedding, h:=ρ v∈ C^0(Ω×[0,T]).
Thus, using <cit.>, the unique solution to

    _tρ = ((1+δρ)ρ - h)  Ω,
    	 t>0,   ρ·ν=0  Ω,

satisfies ρ∈ C^γ,γ/2(Ω×[0,T]).
By elliptic regularity again, v∈ C^2,γ/2(Ω×[0,T]).
Consequently, h∈ C^γ,γ/2(Ω×[0,T])
and applying <cit.> again, we infer that
ρ∈ C^2,1(Ω×[0,T]). It is possible to show that
ρ∈ S for suitable R>0 and M>0. Hence,
the existence of a solution to (<ref>)–(<ref>)
follows from the Schauder fixed-point theorem. 

Elliptic regularity implies that v∈ C^4,1(Ω×[0,T]). 
Then f:=(ρ v)∈ C^1,1(Ω×[0,T]) and
the solution u=ρ to the linear parabolic equation 
_t u-Δ u - (ρ u)=f in Ω, t>0, with no-flux
boundary conditions satisfies u∈ 
C^2+γ,1+γ/2(Ω×[0,T]) <cit.>
(here we need ρ^0∈ C^2+γ(Ω)).
Thus, the regularity of f improves to f∈ C^1+γ,1+γ/2
(Ω×[0,T]). By parabolic regularity 
<cit.>, we infer that 
u∈ C^2+γ,2(Ω×[0,T]). Bootstrapping this
argument and using <cit.>, we find that
ρ=u∈ C^∞(Ω×(0,T]) and consequently
v∈ C^∞(Ω×(0,T]). 

Step 2: Preparations. Let =1, let (ρ_δ,c_δ) be a local
smooth solution to (<ref>)-(<ref>) with 0<δ<1,
and let (ρ,c) be a local smooth solution to (<ref>)–(<ref>).
Then ρ_R:=ρ_δ-ρ and c_R:=c_δ-c solve

    _tρ_R    = (ρ_R - ρ_R(c+c_R) - ρ c_R), 
    	
    _t c_R    = Δ c_R - c_R + δΔ(ρ+ρ_R)
    	+ (ρ+ρ_R)^α - ρ^α  Ω, t>0,

(ρ_R,c_R) satisfies homogeneous Neumann boundary conditions and
vanishing initial conditions:

    ρ_R·ν =  c_R·ν = 0  Ω, t>0,   ρ_R(0) = c_R(0)=0  Ω.

The aim is to prove a differential inequality for

    Γ(t) = (ρ_R,c_R)(t)_H^2(Ω)^2,   
    	G(t) = (ρ_R,c_R)(t)_H^2(Ω)^2 
    	+ Δ(ρ_R,c_R)(t)_L^2(Ω)^2,

where (ρ_R,c_R)^2_X = ρ_R_X^2 + c_R_X^2 
for suitable norms ·_X.

Step 3: H^1(Ω) estimates. We use ρ_R as a test function
in (<ref>):

    1/2d/dtρ_R_L^2(Ω)^2
    	+ ρ_R_L^2(Ω)^2
    	   = ∫_Ωρ_R(c+c_R)·ρ_R dx 
    	+ ∫_Ωρ c_R·ρ_R dx 
       =: I_1 + I_2.

By Young's inequality, for any η>0, we have

    I_2 ≤η/2ρ_R_L^2(Ω)^2
    	+ 1/2ηρ_L^∞(Ω)^2 c_R_L^2(Ω)^2
    	≤η/2ρ_R_L^2(Ω)^2
    	+ C(η) c_R_L^2(Ω)^2,

where here and in the following, C>0 and C(η)>0 denote generic constants
independent of δ but depending on suitable norms of (ρ,c).
The embedding H^2(Ω)↪ L^∞(Ω) (for d≤ 3) gives

    I_1    ≤ηρ_R_L^2(Ω)^2
    	+ 1/2η c_L^∞(Ω)^2ρ_R_L^2(Ω)^2
    	+ 1/2ηρ_R_L^∞(Ω)^2 c_R_L^2(Ω)^2 
       ≤ηρ_R_L^2(Ω)^2
    	+ C(η)ρ_R_L^2(Ω)^2 
    	+ C(η)ρ_R_H^2(Ω)^2 c_R_L^2(Ω)^2.

Combining the estimates for I_1 and I_2 and choosing η>0 sufficiently
small, we find that

    d/dtρ_R_L^2(Ω)^2
    	+ Cρ_R_L^2(Ω)^2
    	≤ C(1+ρ_R_H^2(Ω)^2) c_R_L^2(Ω)^2
    	+ Cρ_R_L^2(Ω)^2.


Next, we use the test function c_R in (<ref>):

    1/2d/dtc_R_L^2(Ω)^2
    	+ c_R_H^1(Ω)^2 
    	   = -δ∫_Ω(ρ+ρ_R)· c_R dx 
    	+ ∫_Ω((ρ+ρ_R)^α-ρ^α) c_Rdx 
       = I_3 + I_4.

Thus, for any η>0,

    I_3    ≤η/2 c_R_L^2(Ω)^2
    	+ δ^2/2η(ρ_L^2(Ω)^2
    	+ρ_R_L^2(Ω)^2) 
       ≤η/2 c_R_L^2(Ω)^2 
    	+ C(η)ρ_R_L^2(Ω)^2 + C(η)δ^2.


For the estimate of I_4, we apply the mean-value theorem to the function
s↦ s^α (recalling that α≥ 1):

    |(ρ+ρ_R)^α-ρ^α| 
    	≤ C(1+ρ_R_L^∞(Ω))^α-1|ρ_R| 
    	≤ C(1+ρ_R_L^∞(Ω)^α-1)|ρ_R|.

Hence, together with the embedding H^2(Ω)↪ L^∞(Ω),

    I_4 ≤ηc_R_L^2(Ω)^2 
    	+ C(η)(1+ρ_R_H^2(Ω)^2(α-1))
    	ρ_R_L^2(Ω)^2.

Collecting these estimates and choosing η>0 sufficiently small, it follows that

    1/2d/dtc_R_L^2(Ω)^2 + Cc_R_H^1(Ω)^2
    	≤ Cρ_R_L^2(Ω)^2
    	+ C(1+ρ_R_H^2(Ω)^2(α-1))ρ_R_L^2(Ω)^2
    	+ Cδ^2.


Thus, summing (<ref>) and (<ref>), 

    d/dt(ρ_R,c_R)_L^2(Ω)^2
    	+ C(ρ_R,c_R)_H^1(Ω)^2
    	≤ C(Γ(t) + Γ(t)^2 + Γ(t)^α) + Cδ^2.


Step 4: H^2(Ω) estimates. We multiply (<ref>) by -Δρ_R
and integrate by parts in the expression with the time derivative:

    1/2d/dtρ_R_L^2(Ω)^2
    	+ Δρ_R_L^2(Ω)^2
    	   = ∫_Ω(ρ_R(c+c_R))Δρ_R dx
    	+ ∫_Ω(ρ c_R)Δρ_R dx 
       = ∫_Ω(ρ_R·(c+c_R)+ρ_RΔ(c+c_R))Δρ_R dx 
       xx+ ∫_Ω(ρ· c_R + ρΔ c_R)Δρ_R dx
    	=: I_5 + I_6.

Then, taking into account inequality (<ref>) in the Appendix,

    I_5    ≤ηΔρ_R_L^2(Ω)^2
    	+ C(η)ρ_R_L^2(Ω)^2 
    	+ C(η)ρ_R_H^1(Ω)^2 c_R_H^1(Ω)^2 
       xx+ C(η)ρ_R_L^2(Ω)^2
    	+ C(η)ρ_R_H^1(Ω)^2Δ c_R_H^1(Ω)^2, 
    
    	I_6    ≤ηΔρ_R_L^2(Ω)^2
    	+ C(η) c_R_L^2(Ω)^2 + C(η)Δ c_R_L^2(Ω)^2,

and choosing η>0 sufficiently small, we end up with

    1/2d/dtρ_R_L^2(Ω)^2 
    	+ CΔρ_R_L^2(Ω)^2
    	   ≤ Cρ_R_H^1(Ω)^2 
    	+ C(1+ρ_R_H^1(Ω)^2)c_R_H^2(Ω)^2 
       xx+ Cρ_R_H^1(Ω)^2Δ c_R_H^1(Ω)^2.


We multiply (<ref>) by -Δ c_R and estimate similarly as in Step 3:

    1/2d/dt    c_R_L^2(Ω)^2
    	+ Δ c_R_L^2(Ω)^2 +  c_R_L^2(Ω)^2 
       = δ∫_ΩΔ(ρ+ρ_R)Δ c_R dx
      + ∫_Ω((ρ+ρ_R)^α-ρ^α)Δ c_R dx 
       ≤ηΔ c_R_L^2(Ω)^2
    	+ C(1+ρ_R_H^2(Ω)^2(α-1))ρ_R_L^2(Ω)^2
    	+ Cδ^2Δρ_R_L^2(Ω)^2 + Cδ^2.

Hence, for sufficiently small η>0,

    1/2d/dt    c_R_L^2(Ω)^2
    	+ Δ c_R_L^2(Ω)^2 +  c_R_L^2(Ω)^2 
       ≤ C(1+ρ_R_H^2(Ω)^2(α-1))ρ_R_L^2(Ω)^2
    	+ Cδ^2Δρ_R_L^2(Ω)^2 + Cδ^2.

Adding this inequality and (<ref>), adding c_R_L^2(Ω)^2 on
both sides, using (<ref>) in the Appendix, and choosing δ>0
sufficiently small to absorb the term Cδ^2Δρ_R_L^2(Ω)^2, 
we infer that

    d/dt(ρ_R,c_R)_L^2(Ω)^2
    	   + C(ρ_R,c_R)_H^2(Ω)^2 
    	≤ C(ρ_R,c_R)_H^2(Ω)^2
    	+ Cρ_R_H^1(Ω)^2 c_R_H^1(Ω)^2 
       + Cρ_R_H^1(Ω)^2Δ c_R_H^1(Ω)^2
    	+ Cρ_R_H^2(Ω)^2(α-1)ρ_R_L^2(Ω)^2 
    	+ Cδ^2 
       ≤ C(Γ(t) + Γ(t)^2 + Γ(t)^α + Γ(t)G(t))
    	+ Cδ^2.


Step 5: H^3(Ω) estimates.
We apply the Laplacian to (<ref>) and (<ref>) and multiply both equations
by Δρ_R, Δ c_R, respectively:

    1/2d/dt   Δρ_R_L^2(Ω)^2
    	+ Δρ_R_L^2(Ω)^2 
    	= ∫_Ω(ρ_R(c+c_R))·Δρ_R dx 
       + ∫_Ω(ρ c_R)·Δρ_R dx =: I_7 + I_8, 
    1/2d/dt   Δ c_R_L^2(Ω)^2
    	+ Δ c_R_L^2(Ω)^2 + Δ c_R_L^2(Ω)^2 
      = δ∫_ΩΔ(ρ+ρ_R)·Δ c_R dx 
       + ∫_Ω((ρ+ρ_R)^α-ρ^α)·Δ c_R dx
    	=: I_9 + I_10.

We estimate

    I_8    ≤ηΔρ_R_L^2(Ω)^2
    	+ C(η)(ρ c_R)_L^2(Ω)^2
    	+ C(η)ρ_R_H^2(Ω)^2 
       ≤ηΔρ_R_L^2(Ω)^2
    	+ C(η)c_R_H^2(Ω)^2
    	+ C(η)Δ c_R_L^2(Ω)^2
    	+ C(η)ρ_R_H^2(Ω)^2.

Taking into account

    (ρ_R c_R)
    	   = (ρ_R· c_R + ρ_RΔ c_R) 
       = ( c_R·)ρ_R + (ρ_R·) c_R
    	+ ρ_RΔ c_R + ρ_RΔ c_R,

and inequalities (<ref>) and (<ref>) in the Appendix
as well as the embedding H^2(Ω)↪ L^∞(Ω), we obtain

    (ρ_R c_R)_L^2(Ω)^2
    	   ≤ c_R_L^∞(Ω)^2^2ρ_R_L^2(Ω)^2
    	+ ρ_R_H^1(Ω)^2^2 c_R_H^1(Ω)^2 
       xx+ ρ_R_H^1(Ω)^2Δ c_R_H^1(Ω)^2
    	+ ρ_R_L^∞(Ω)^2Δ c_R_L^2(Ω)^2 
       ≤ Cρ_R_H^2(Ω)^2(Δ c_R_L^2(Ω)^2
    	+ c_R_H^2(Ω)^2),

which gives

    I_7 ≤ηΔ c_R_L^2(Ω)^2
    	+ C(η)ρ_R_H^2(Ω)^2 
    	+ C(η)ρ_R_H^2(Ω)^2(Δ c_R_L^2(Ω)^2
    	+ c_R_H^2(Ω)^2).

This shows that, again for sufficiently small η>0,

    1/2d/dtΔρ_R_L^2(Ω)^2
    	+ CΔρ_R_L^2(Ω)^2 
         ≤ C(ρ_R,c_R)_H^2(Ω)^2 + C(1+ρ_R_H^2(Ω)^2)
    	Δ c_R_L^2(Ω)^2 
       xx+ Cρ_R_H^2(Ω)^2c_R_H^2(Ω)^2.


Next, we estimate I_9 and I_10:

    I_9    ≤ηΔ c_R_L^2(Ω)^2 
    	+ C(η)δ^2Δρ_R_L^2(Ω)^2 + Cδ^2, 
    
      I_10   = α∫_Ω((ρ+ρ_R)^α-1-ρ^α-1)
    	ρ·Δ c_R dx 
       xx
    	+ α∫_Ω(ρ+ρ_R)^α-1ρ_R·Δ c_R dx 
    	= J_1 + J_2.

We find that

    J_2 ≤ηΔ c_R_l^2(Ω)^2
    	+ C(η)(1+ρ_R_H^2(Ω)^2(α-1))
    	ρ_R_L^2(Ω)^2.

By the Hölder continuity of s↦ s^α-1, it follows that

    J_2    ≤ C∫_Ω|ρ_R|^α-1|Δ c_R|dx
    	≤ηΔ c_R_L^2(Ω)^2
    	+ C(η)ρ_R_L^2(α-1)(Ω)^2(α-1)
       ≤ηΔ c_R_L^2(Ω)^2
    	+ C(η)ρ_R_H^2(Ω)^2(α-1).

Consequently, for sufficiently small η>0,

    1/2d/dt   Δ c_R_L^2(Ω)^2
    	+ CΔ c_R_L^2(Ω)^2 + Δ c_R_L^2(Ω)^2 
       ≤ Cρ_R_H^2(Ω)^2(α-1)
    	+ C(1+ρ_R_H^2(Ω)^2(α-1))ρ_R_L^2(Ω)^2
    	+ Cδ^2Δρ_R_L^2(Ω)^2 + Cδ^2.


Adding the previous inequality and (<ref>) and taking δ>0 sufficiently
small such that the term Cδ^2Δρ_R_L^2(Ω)^2 is absorbed
by the corresponding term on the left-hand side of (<ref>), we infer that

    d/dt   (ρ_R,c_R)_L^2(Ω)^2
    	+ CΔ(ρ_R,c_R)_L^2(Ω)^2 + CΔ c_R_L^2(Ω)^2 
    	
       ≤ C(ρ_R,c_R)_H^2(Ω)^2	+ C(1+ρ_R_H^2(Ω)^2)
    	Δ c_R_L^2(Ω)^2 
    	+ Cρ_R_H^2(Ω)^2c_R_H^2(Ω)^2 
       xx+ Cρ_R_H^2(Ω)^2(α-1)
    	+ C(1+ρ_R_H^2(Ω)^2(α-1))ρ_R_L^2(Ω)^2
    	+ Cδ^2 
       ≤ C(Γ(t) + Γ(t)^α-1 + Γ(t)^2 + Γ(t)^α
    	+ Γ(t)G(t)) + Cδ^2.


Step 6: End of the proof for =1.
We sum inequalities (<ref>), (<ref>), and (<ref>):

    d/dt(   (ρ_R,c_R)_H^1(Ω)^2
    	+ Δ(ρ_R,c_R)_L^2(Ω)^2) 
       xx+ C((ρ_R,c_R)_H^2(Ω)^2
    	+ Δ(ρ_R,c_R)_L^2(Ω)^2) 
       ≤ C(Γ(t) + Γ(t)^α-1 + Γ(t)^2 + Γ(t)^α
    	+ Γ(t)G(t)) + Cδ^2.

To get rid of the term Γ(t)^α-1, we need the condition α≥ 2.
Indeed, under this condition,

    Γ(t)^α-1≤Γ(t) + Γ(t)^α.

We also remove the term Γ(t)^2 by defining 
κ:= max{α,2} and estimating

    Γ(t)^2 ≤Γ(t) + Γ(t)^κ.

We deduce from elliptic regularity that

    Γ(t) ≤ CΔ(ρ_R,c_R)_L^2(Ω)^2
    	+ C(ρ_R,c_R)_H^1(Ω)^2.

Therefore, integrating (<ref>) over (0,t) and observing that
(ρ_R,c_R)(0)=0, (<ref>) becomes

    Γ(t) + C∫_0^t G(s)ds ≤ C∫_0^t(Γ(s)+Γ(s)^κ)ds
    	+ C∫_0^t Γ(s)G(s)ds + Cδ^2.

Lemma <ref> proves the result for =1.

Step 7: Parabolic-elliptic case =0.
Since there is no time derivative of c_R anymore, we need to change the
definition of the functionals Γ(t) and G(t):

    Γ_0(t) = ρ_R_H^2(Ω)^2,   
    	G_0(t) = ρ_R_H^2(Ω)^2 + Δρ_R_L^2(Ω)^2.

The estimates are very similar to the parabolic-parabolic case with two exceptions:
In (<ref>), we have estimated the terms ρ_R_H^2(Ω)^2
Δ c_R_L^2(Ω)^2 and
ρ_R_H^2(Ω)^2c_R_H^2(Ω)^2 from above by Γ(t)G(t).
In the present case, we cannot estimate Δ c_R_L^2(Ω)^2 
by G_0(t) and we need to proceed in a different way.

Estimates (<ref>) and (<ref>), adapted to the case =0, become

    Δ c_R_L^2(Ω)^2 +  c_R_L^2(Ω)^2
    	   ≤ C(1+ρ_R_H^2(Ω)^2(α-1))ρ_R_L^2(Ω)^2
    	+ Cδ^2Δρ_R_L^2(Ω)^2 + Cδ^2 
       ≤ C(Γ_0(t) + Γ_0(t)^α) + Cδ^2, 
    Δ c_R_L^2(Ω)^2 + Δ c_R_L^2(Ω)^2 
     	   ≤ Cρ_R_H^2(Ω)^2(α-1)
    	+ C(1+ρ_R_H^2(Ω)^2(α-1))ρ_R_L^2(Ω)^2 
       xx+ Cδ^2Δρ_R_L^2(Ω)^2 + Cδ^2.

The term δ^2Δρ_R_L^2(Ω)^2 can be absorbed
by the corresponding term on the left-hand side of (<ref>). The
critical term ρ_R_H^2(Ω)^2(α-1)ρ_R_L^2(Ω)^2
is bounded from above by Γ(t)^α. Thus, Δ c_R_L^2(Ω)^2
is estimated by Γ_0(t)^α and lower-order terms, and consequently,
ρ_R_H^2(Ω)^2 ×Δ c_R_L^2(Ω)^2 in (<ref>)
is estimated by Γ_0(t)^α G_0(t), together with lower-order terms. 
Furthermore, ρ_R_H^2(Ω)^2c_R_H^2(Ω)^2 is bounded by
Γ_0(t)^α+1, up to lower-order terms. More precisely,
a computation shows that

    d/dt   (ρ_R_H^1(Ω)^2 + Δρ_R_L^2(Ω)^2)
    	+ C(ρ_R_H^2(Ω)^2 + Δρ_R_L^2(Ω)^2)
    	+ CG_0(t) 
       ≤ C(Γ_0(t) + Γ_0(t)^2 + Γ_0(t)^α-1 + Γ_0(t)^α
    	+ Γ_0(t)^α+1 + Γ_0(t)G_0(t) 
       xx+ Γ_0(t)^α G_0(t)
    	+ δ^2 G_0(t)) + Cδ^2.

Observing that 
Γ_0(t)≤ C(ρ_R_H^1(Ω)^2 + Δρ_R_L^2(Ω)^2)
and Γ_0(t)^α-1+Γ_0(t)^α≤Γ_0(t)+Γ_0(t)^α+1, 
choosing δ>0 sufficiently small, integrating in time, and using
Γ_0(0)=0, we arrive at

    Γ_0(t) + C∫_0^t G_0(s)ds 
    	   ≤ C∫_0^t(Γ_0(s)+Γ_0(t)^α+1)ds 
       xx+ C∫_0^t(Γ_0(s)+Γ_0(s)^α+1)G_0(s)ds + Cδ^2.

An application of Lemma <ref> finishes the proof.





§ NUMERICAL EXPERIMENTS


We present some numerical examples for system (<ref>)–(<ref>) 
in two space dimensions and for various choices of the parameters 
α and δ. Equations (<ref>) 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 <cit.>. We do not use any kind of additional regularizations,
smoothing tools, or slope-limiters. All experiments are performed for the
parabolic-parabolic equations with δ>0.

We choose the same domain and initial conditions as in <cit.>, i.e.
Ω={x∈^2:|x|<1} and 

    ρ^0(x,y) = 80(x^2+y^2-1)^2(x-0.1)^2 + 5,   
    	c^0(x,y) = 0,    (x,y)∈Ω.

A computation shows that the total mass M=∫_Ωρ^0 dx=25π/3>8π
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
∫_Ω|x-x_0|^2ρ^0 dx is sufficiently small for some x_0∈Ω.
Blow up at the boundary can occur if x_0∈Ω and M>4π.

Experiment 1: α=1. We choose the initial datum (<ref>)
and the values α=1, δ=10^-3. In this nonsymmetric setting, 
the solution exists for all time and the density is expected to concentrate 
at the boundary <cit.>.
Figure <ref> 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^∞ blow-up. 
The steady state is reached at approximately T=2.5.
By Theorem <ref>, the peak approximates the blow-up solution to the
classical Keller–Segel system in the L^∞ norm;
see Figure <ref>.
We see that the L^∞ norm of the density becomes larger with decreasing values
of δ.





Experiment 2: α>1. First, we choose
the value α=1.5. The initial datum is still given by (<ref>). 
Since α>1, we cannot exclude finite-time blow-up, which is
confirmed by the numerical experiments in Figure <ref>. Numerically,
the solution seems to exist until time T^*≈ 0.079. The numerical scheme
breaks down at slightly smaller times when δ 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 α. Indeed, Figure <ref> shows a stronger and faster
concentration behavior when we take α=2.5.





Experiment 3: Multi-bump initial datum. We take α=1 and
δ=5· 10^-3. As initial datum, we choose a linear combination of
the bump function

    W_x_0,y_0,M(x,y) = M/2πθexp(
    	-(x-x_0)^2+(y-y_0)^2/2θ),   (x,y)∈Ω,

where (x_0,y_0)∈Ω, M>0, and θ>0. Setting θ=10^-2,
we define c^0=0 and

    ρ^0    = W_0.25,0,10π + W_-0.25,0,4π + W_0,-0.25,4π + W_0,0.25,4π
       xx+ W_0,0.5,4π + W_0,0.35,4π + W_0.5,0,4π 
    	+ W_0.5,0.25,4π.

The evolution of the density is presented in Figure <ref>. 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 <cit.> 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 δ. 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 δ.
We choose α=1 and the initial datum ρ^0=W_0,0,20π with θ=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 <ref> (top row) shows the stationary cell densities for two values of
δ. 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 δ. The level sets show that the solutions
are almost radially symmetric and the level set for ρ=10^-2 is approximately
a circle. This behavior is quantified in Figure
<ref> (bottom row). We observe that the radius depends on δ
approximately as r∼δ^0.43 and the height approximately as
ρ_ max∼δ^-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 <ref>),
since the chemical substance is not absorbed by the boundary as in the
Dirichlet case.








§ SOME TECHNICAL TOOLS


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


Let d≤ 3, Ω⊂^d be a bounded domain, and Ω∈ C^2,1. 
There exists a constant C>0 such that for all u, v∈ H^1(Ω),

    uv_L^2(Ω)≤ Cu_H^1(Ω)v_H^1(Ω),

for all u∈ H^2(Ω) with u·ν=0 on Ω,

    u_H^2(Ω)^2 ≤ C(Δ u_L^2(Ω)^2
    	+u_L^2(Ω)^2),

and for all u∈ H^3(Ω) with u·ν=0 on Ω,

    u_H^3(Ω)^2 ≤ C(Δ u_L^2(Ω)^2
    	+u_H^2(Ω)^2).



Inequality (<ref>) follows after applying the Cauchy–Schwarz inequality
and then the continuous embedding H^1(Ω)↪ L^4(Ω); 
(<ref>) is proved
in <cit.>, while (<ref>) is a consequence of
<cit.>.


Let δ>0 and Γ, G∈ C^0([0,T]) be nonnegative functions,
possibly depending on δ, satisfying

    Γ(t) + C_0∫_0^t G(s)ds
    	   ≤ C_1Γ(0) + C_2∫_0^t(Γ(s)+Γ(s)^α)ds 
       xx
    	+ C_3δ^β∫_0^t(Γ(s)+Γ(s)^γ)G(s)ds + C_4δ^ν,

where α>1, β≥ 0, γ>0, ν>0, and C_0,…,C_4>0 
are constants independent of δ. Furthermore, let 
Γ(0)≤ C_5δ^ν for some C_5>0. Then there exists
δ_0>0 such that for all 0<δ<δ_0, 0≤ t≤ T, and
0<<ν,

    Γ(t) ≤ C_5δ^ν-.




A slightly simpler variant of the lemma was proved in <cit.>.
Assume, by contradiction, that for all δ_0∈(0,1), there exist 
δ∈(0,δ_0), t_0∈[0,T], and ∈(0,ν) such that
Γ(t_0)> C_5δ^ν-. Since Γ(0)≤ C_5δ^ν by
assumption and Γ is continuous, there exists t_1∈[0,t_0) such that
Γ(t_1)=C_5δ^ν- and Γ(t)≤ C_5δ^ν- for
all t∈[0,t_1]. This leads for t∈[0,t_1] to

    Γ(t) + C_0∫_0^t G(s)ds 
    	   ≤ C_1C_5δ^ν
    	+ C_2(1+(C_5δ^ν-)^α-1)∫_0^tΓ(s)ds 
       xx
    	+ C_3δ^β(C_5δ^ν- + (C_5δ^ν-)^γ)
    	∫_0^t G(s)ds + C_4δ^ν.

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

    Γ(t) ≤ (C_1C_5+C_4)δ^ν + 2C_2∫_0^tΓ(s)ds,   
    	0≤ t≤ t_1.

Then Gronwall's lemma gives, for sufficiently small δ_0>0 and
0<δ<δ_0,

    Γ(t) ≤ (C_1C_5+C_4)δ^ν e^2C_2 T≤C_5/2δ^ν-
    	< C_5δ^ν-,    0≤ t≤ t_1.

which contradicts Γ(t_1)=C_5δ^ν-.






11
Ama93 H. Amann. Nonhomogeneous linear and quasilinear elliptic and parabolic 
boundary value problems. In: H. J. Schmeisser and H. Triebel (eds), 
Function Spaces, Differential Operators and Nonlinear Analysis, pp. 9–126. 
Teubner, Stuttgart, 1993.

BaCa15 D. Bartolucci and D. Castorina. A global existence result for a
Keller–Segel type system with supercritical initial data.
J. Ellipt. Parabol. Eqs. 1 (2015), 243–262.

BBTW15 N. Bellomo, A. Bellouquid, Y. Tao, and M. Winkler. Toward a 
mathematical theory of Keller–Segel models of pattern formation in biological tissues. 
Math. Models Meth. Appl. Sci. 25 (2015), 1663–1763.

BeJu13 M. Bessemoulin-Chatard and A. Jüngel. A finite volume scheme for 
a Keller–Segel model with additional cross-diffusion. 
IMA J. Numer. Anal. 34 (2013), 96–122.

Bla12 A. Blanchet. On the parabolic-elliptic Patlak–Keller–Segel system 
in dimension 2 and higher. Séminaire L. Schwartz – EDP Appl. 8 (2011–2012),
26 pages.

BHLM96 A. Bonami, D. Hilhorst, E. Logak, and M. Mimura. A free boundary 
problem arising in a chemotaxis model. In: Free Boundary Problems, Theory and 
Applications (Zakopane, 1995). Pitman Res. Notes Math. Ser. 363, pp. 363–373. 
Longman, Harlow, 1996.

BDD06 M. Burger, M. Di Francesco, and Y. Dolak-Struss. The Keller–Segel
model for chemotaxis with prevention of overcrowding: linear vs. nonlinear
diffusion. SIAM J. Math. Anal. 38 (2006), 1288–1315.

CaCa06 V. Calvez and J. A. Carrillo. Volume effects in the Keller–Segel
model: energy estimates preventing blow-up. J. Math. Pure Appl. 86 (2006),
155–175.

CaCa08 V. Calvez and L. Corrias. The parabolic-parabolic Keller–Segel 
model in ^2. Commun. Math. Sci. 6 (1008), 417–447.

CHJ12 J. A. Carrillo, S. Hittmeir, and A. Jüngel. Cross diffusion and 
nonlinear diffusion preventing blow up in the Keller–Segel model. 
Math. Models Meth. Appl. Sci. 22, (2012), 1250041, 35 pages. 




CPZ04 L. Corrias, B. Perthame, and H. Zaag. Global solutions of some
chemotaxis and angiogenesis systems in high space dimensions.
Milan J. Math. 72 (2004), 1–28.

DiRo08 M. Di Francesco and J. Rosado. Fully parabolic Keller–Segel model
for chemotaxis with prevention of overcrowding. Nonlinearity 21 (2008), 
2715–2730.

DrJu12 M. Dreher and A. Jüngel. Compact families of piecewise constant 
functions in L^p(0,T;B). Nonlin. Anal. 75 (2012), 3072–3077.

Fri08 A. Friedman. Partial Differential Equations. Dover, Mineola,
2008.

Gri85 P. Grisvard. Elliptic Problems in Nonsmooth Domains.
Pitman, Boston, 1985.

HiPa09 T. Hillen and K. Painter. A user's guide to PDE models for chemotaxis. 
J. Math. Biol. 58 (2009), 183–217.

HPS07 T. Hillen, K. Painter, and C. Schmeiser. Global existence for
chemotaxis with finite sampling radius. Discrete Contin. Dyn. Sys.
7 (2007), 125–144.

HiJu11 S. Hittmeir and A. Jüngel. Cross diffusion preventing blow up in the
two-dimensional Keller–Segel model. SIAM J. Math. Anal. 43 (2011), 997–1022.

HsWa06 L. Hsiao and S. Wang. Quasineutral limit of a time-dependent
drift-diffusion-Poisson model for p-n junction semiconductor devices.
J. Diff. Eqs. 225 (2006), 411–439.

Hun07 J. Hunter. Matplotlib: A 2D graphics environment. 
Comput. Sci. Engin. 9 (2007), 90–95.

Jue15 A. Jüngel. The boundedness-by-entropy method for cross-diffusion
systems. Nonlinearity 28 (2015), 1963–2001.

Jue16 A. Jüngel. Entropy Methods for Diffusive Partial Differential
Equations. SpringerBriefs in Mathematics, Springer, 2016.

Jue17 A. Jüngel. Cross-diffusion systems with entropy structure. 
In: K. Mikula, D. Ševčovič, and J. Urbán (eds.), 
Proceedings of Equadiff 2017 Conference, pp. 181–190. 
Spektrum Stu Publishing, Bratislava, 2017.

KeSe70 E. Keller and L. Segel. Initiation of slime mold aggregation viewed 
as an instability. J. Theor. Biol. 26 (1970), 399–415.

Kow05 R. Kowalczyk. Preventing blow-up in a chemotaxis model.
J. Math. Anal. Appl. 305 (2005), 566–588.

KoSz08 R. Kowalczyk and Z. Szymańska. On the global existence of
solutions to an aggregation model. J. Math. Anal. Appl. 343 (2008), 379–398.

Lan17 J. Lankeit. Locally bounded global solutions to a chemotaxis consumption 
model with singular sensitivity and nonlinear diffusion. 
J. Diff. Eqs. 262 (2017), 4052–4084.

Lun95 A. Lunardi. Analytic Semigroups and Optimal Regularity in 
Parabolic Problems. Springer, Basel, 1995.

Pat53 C. Patlak. Random walk with persistence and external bias. 
Bull. Math. Biophys. 15 (1953), 311–338.

Sim87 J. Simon. Compact sets in the space L^p(0,T;B).
Ann. Math. Pura. Appl. 146 (1987), 65–96.

Tro87 G. M. Troianiello. Elliptic Differential Equations and
Obstacle Problems. Plenum Press, New York, 1987.

Win08 M. Winkler. Chemotaxis with logistic source: Very weak solutions
and their boundedness properties. J. Math. Anal. Appl. 348 (2008), 708–729.

Win18 M. Winkler. A critical blow-up exponent in a chemotaxis system with 
nonlinear signal production. Nonlinearity 31 (2018), 2031–2056.




