On the Existence of Global Smooth Solutions to the Parabolic–Elliptic Keller–Segel System with Irregular Initial Data

We consider the parabolic–elliptic Keller–Segel system *ut=Δu-χ∇·(u∇v),0=Δv-v+u\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \left\{ \begin{aligned} u_t&= \Delta u - \chi \nabla \cdot (u \nabla v), \\ 0&= \Delta v - v + u \end{aligned} \right. \end{aligned}$$\end{document}in a smooth bounded domain Ω⊆Rn\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega \subseteq {\mathbb {R}}^n$$\end{document}, n∈N\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n\in {\mathbb {N}}$$\end{document}, with Neumann boundary conditions. We look at both chemotactic attraction (χ>0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\chi > 0$$\end{document}) and repulsion (χ<0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\chi < 0$$\end{document}) scenarios in two and three dimensions. The key feature of interest for the purposes of this paper is under which conditions said system still admits global classical solutions due to the smoothing properties of the Laplacian even if the initial data is very irregular. Regarding this, we show for initial data μ∈M+(Ω¯)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu \in {\mathcal {M}}_+({\overline{\Omega }})$$\end{document} that, if either n=2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n = 2$$\end{document}, χ<0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\chi < 0$$\end{document} or n=2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n = 2$$\end{document}, χ>0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\chi > 0$$\end{document} and the initial mass is small or n=3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n = 3$$\end{document}, χ<0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\chi < 0$$\end{document} and μ=f∈Lp(Ω)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu = f \in L^p(\Omega )$$\end{document}, p>1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p > 1$$\end{document} holds, it is still possible to construct global classical solutions to (⋆\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\star $$\end{document}), which are continuous in t=0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t = 0$$\end{document} in the vague topology on M+(Ω¯)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {M}}_+({\overline{\Omega }})$$\end{document}. n=2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n = 2$$\end{document}, χ<0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\chi < 0$$\end{document} or n=2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n = 2$$\end{document}, χ>0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\chi > 0$$\end{document} and the initial mass is small or n=3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n = 3$$\end{document}, χ<0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\chi < 0$$\end{document} and μ=f∈Lp(Ω)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu = f \in L^p(\Omega )$$\end{document}, p>1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p > 1$$\end{document}


Introduction
In this paper, we are primarily concerned with systems of partial differential equations used in the study of biological systems.More specifically, we are interested in systems modeling chemotaxis, the directed movement of cells along a chemical gradient.This use of partial differential equations in the biological study of chemotactic processes was mostly initiated by the seminal work of Keller and Segel in 1970 (cf.[12]), in which Keller and Segel used them to model certain slime molds in an effort to understand their aggregation behavior.The popularity of this approach was further bolstered by the subsequent successful mathematical analysis of said model, which confirmed the presence of aggregation in the sense that under appropriate initial conditions the solutions to the system blow up in finite time while retaining their initial mass (cf.[24], [39]).This success in modeling and mathematical analysis has led to many more biological processes (and sometimes even comparable processes from other fields, such as criminology, cf.[35]) to be modeled in a similar fashion.For a survey, we refer the reader to [2].
As already alluded to, there exist many scenarios (generally dependent on the initial data and dimension of the domain), in which solutions to these kinds of chemotaxis models blow up in finite time (cf.[9], [10], [21], [24], [33], [39]).It is further known that in some of these scenarios this blowup takes the form of the solution converging to approximately a Dirac mass or some function of potentially very little regularity (cf.[25], [36]).What we are now interested in is in a sense the opposite scenario, which has to our knowledge been less often discussed thus far.Namely, we consider the case of starting with initial data of very little regularity, such as a Dirac measure, and are then concerned with deriving whether or under which conditions there still exist sensible global classical solutions for such a model, which are still connected to the initial data in a reasonable fashion.This can be interpreted as essentially starting our analysis at the point in time when aggregation occurred and then investigating under which circumstances the model still yields sensible results from that point onward.
The model we want to analyze in this regard will be a variation on the original (here somewhat simplified) Keller-Segel model (cf.[12]): In this model, the function u represents the density of the organism under consideration while the function v represents the density of the attractant substance.Both are under the influence of diffusion modeled by the terms ∆u and ∆v.The central term representing the chemotatic interaction is ∇ • (u∇v).The remaining linear terms in the second equation then model the degradation of the attractant over time and the production of the attractant by the organisms, respectively.
We then consider the similar system (with the same roles for u and v) in which the second equation is only of elliptic type.This modification of the original Keller-Segel system is in fact not uncommon (cf.[10], [21]) and can be understood as reflecting that the chemical responds immediately everywhere to changes in the population density of the modeled organism as opposed to the organism density only influencing its time evolution and therefore having a delayed effect on it.Additionally, we also add a potentially negative coefficient χ to the chemotaxis term ∇ • (u∇v), which models that the chemical substance cannot only be an attractant but also possibly a repellent (cf.[17], [18], [29] for some biological processes involving chemorepulsion).
Main result.We consider the system (1.1) with Neumann boundary conditions ∇u(x, t) • ν = 0 for all x ∈ ∂Ω, t > 0, ∇v(x, t) • ν = 0 for all x ∈ ∂Ω, t > 0 (1.2) in a bounded domain Ω ⊆ R n , n ∈ {2, 3}, with smooth boundary.Concerning the initial data, we assume µ to be an element of M + (Ω), the set of all positive Radon measures with the vague topology.The vague topology on M + (Ω) is characterized as follows (cf.[1, Definition 30.1]):A sequence (µ n ) n∈N ⊆ M + (Ω) converges to µ ∈ M + (Ω) in the vague topology if and only if For the purposes of this paper, whenever necessary we identify nonnegative functions ϕ ∈ L 1 (Ω) with the measure ϕ(x) dx, where dx represents the standard Lebesgue measure.
Under these conditions, we then investigate some scenarios under which the smoothing properties of the Laplacian are sufficient to counteract the irregularity of the initial data and the destabilizing effects of the taxis term and therefore make it possible to still construct smooth solutions that attain the initial data in a sensible way.These scenarios are repulsive chemotaxis in two and three dimensions (with slightly higher regularity needed for the initial data in the three dimensional case) and attractive chemotaxis in two dimensions (with an additional initial mass condition).More precisely, we prove the following or n = 3 and χ < 0 and µ = f for some f ∈ L p (Ω) with p > 1, (S3) that solve (1.1) on Ω × (0, ∞) with boundary conditions (1.2) classically and attain the initial data µ in the following way: Prior work.To give some context for our existence result, we will now give a brief overview over some notable prior work in this area.
We begin by reviewing some results concerning smooth initial data as opposed to the irregular initial data considered here.This of course makes the construction of solutions easier and thus, especially in the repulsive case, there are quite strong existence results available.Namely, it can be shown that problems of type (1.1) with χ < 0 and smooth initial data have global classical solutions in domains of arbitrary dimension, which converge to their steady states at an exponential rate (cf.[19], [20], [37]).The attractive case is somewhat more complex regarding existence theory as existence here generally depends on properties of the initial data.
In two dimensions existence centrally depends on the initial mass, while in higher dimensions existence can only be ensured for much stronger initial data smallness conditions (cf.[11], [22], [23], [26]).This already suggests that the mass condition for the two-dimensional attractive case (S2) is certainly necessary.
While we are not as concerned with the parabolic-parabolic case, we still want to mention that similar, but maybe not always quite as strong, results are available in this case as well (cf.[5] for existence results in the repulsive case and [9], [27], [28], [39] for discussions of the attractive case to only list a few).We again refer to the survey [2] for a broader overview.
There have also been efforts to analyze chemotaxis systems with two associated elliptic or parabolic equations with one modeling an attractant and the other a repellent with the key result being that existence of solutions is ensured as long are the repellent forces as stronger than their attractive counterparts (cf.[37]).We mention this result as it already illustrates how repulsive chemotaxis generally poses less of problem when constructing solutions as opposed to its attractive counterpart, which is in a sense mirrored in our result.
We now transition to some prior work concerning results about systems similar to (1.1) with irregular initial data.For the two-dimensional whole space case, there are in fact existence results available for a system similar to (1.1) with measure valued initial data with results e.g. based on methods from harmonic analysis (cf.[3], [31]).Moreover, a weak solution construction on the torus is presented in [34].In [40], a system similar to (1.1) with an added logistic source is discussed under the assumption that the initial data is radially symmetric and has a singularity in x = 0 but is otherwise fairly regular.As all of these results restrict themselves to very specific settings in an effort to make use of these restrictions to construct solutions, we can generally not translate the methods employed in them to our setting, which is concerned with general bounded domains and fairly general initial data.
As for the parabolic-parabolic case for the classic Keller-Segel model in bounded domains, it is know that smooth solutions with irregular initial data exist in either the one dimensional case (cf.[41]) or in the two dimensional case with an added logistic source term acting as an additional regularizing factor (cf. [15]).

Approach.
As is often the case when constructing solutions, our basic approach will be looking at approximate solutions (u ε , v ε ) ε∈(0,1) that solve a certain regularized version of (1.1), which can easily be seen to admit global classical solutions, and then gaining our desired solutions (u, v) as limits of (u ε , v ε ) as ε ց 0.
The key regularizations employed by us for this approach are approximating the initial data by smooth functions and replacing the linear growth term u in the second equation by a term that is bounded independent of u, but approaches the original linear term as ε ց 0. For the exact system, see (2.1).
Our next step then is deriving a priori estimates for the approximate solutions that do not depend on ε.This is made particularly challenging due to the fact that we cannot rely on much initial data regularity, which is normally a key part of most testing or semigroup based approaches.When using testing based methods, this is due to the fact that the a priori information gained using them is often based on deriving ordinary differential equations with at best a linear decay term for some of the norms of the solution components and then using comparison arguments, which still take the norm of the initial data into account.
As such, our testing approaches are focused on deriving ordinary differential equations for terms of the form Ω u p ε , which have a superlinear decay term, see Lemma 4.3 and Lemma 4.5.It is the arguments used in the proofs of both of these lemmas where most of the restrictions on the allowed dimension n and values of χ, as well as the initial mass restriction in Theorem 1.1 originate from (Only the need for higher regularity of the initial data in the three dimensional case stems from a later argument in Lemma 7.1, which is apparently necessary to ensure that the constructed solutions still attain the initial data in a sensible fashion).As proven in Lemma 3.1, these ordinary differential equations then allow us to gain uniform (in regards to ε) u ε L p (Ω) bounds for all p ∈ [1, ∞) on (t 0 , ∞), t 0 > 0, for the approximate solutions, which by standard bootstrap arguments combined with some compact embedding properties of Hölder spaces and standard regularity theory yield sufficiently regular classical solutions (u, v) of (1.1) and (1.2) along a suitable sequence ε j ց 0. Additionally, the same argument also gives us certain uniform time integrability properties for Ω u p ε on (0, 1), which are then used in Section 7 to conclude that Ω u εt ϕ has similar uniform time integrability properties.This implies that the approximate solutions are continuous in t = 0 regarding the vague topology in a uniform sense.In Lemma 7.3, we then use this to argue that this continuity therefore survives the limit process and is thus still present in the actual solutions (u, v).

A regularized version of (1.1) with approximated initial data
The key to our construction of solutions to the system (1.1) with irregular initial data will lie in framing said solutions as the limits of approximate solutions to a similar system, which is regularized in two key ways to make the existence of classical solutions much more obvious.The first regularization we employ will be to approximate the measure-valued initial data by smooth functions while the second is replacing the linear growth term u in the second equation of (1.1) by a term that is always uniformly bounded independent of the value of u but results in the original term after a limit processes.
More precisely, we will use the following approximate system with smooth initial data u 0,ε : For this system, it is fairly straightforward to construct unique global classical solutions by first constructing local solutions by standard contraction mapping methods (cf.[8]) and then arguing that finite-time blowup is in fact impossible.The latter step is made easier by the fact that uε 1+εuε ≤ 1 ε , which immediately gives us quite strong bounds for the second solution component.
As both this approach and the employed regularization are quite standard, we will only give the following existence argument in brief.
Lemma 2.1.Let n ∈ N, χ ∈ R and let Ω ⊆ R n be a bounded domain with smooth boundary.Then for ε ∈ (0, 1) and nonnegative u 0,ε ∈ C ∞ (Ω), there exist nonnegative functions that are a global classical solution of (2.1) with the following additional mass conservation property: Proof.By an adaption of standard contraction mapping and maximum principle arguments used in similar settings as seen e.g. in [8, Proposition 3.1], we gain a maximal T max ∈ (0, ∞] and nonnegative functions that are a classical solution to (2.1) on [0, T max ) and adhere to (2.2) on [0, T max ).As a consequence of this standard construction, we further know that, if We will now briefly sketch why this blowup of the L ∞ (Ω) norm of u ε in finite time is in fact impossible, which is of course sufficient to complete this proof.First and foremost due to the fact that uε 1+εuε ≤ 1 ε , standard elliptic regularity theory (cf.[7,Theorem 19.1]) applied to the second equation in (2.1) immediately yields a W 2,p (Ω) bound for v ε on [0, T max ) for all p ∈ (1, ∞), which by the well-known embedding properties of Sobolev spaces in turn translates to a W 1,∞ (Ω) bound for v ε on [0, T max ).With this, we can now use the variation-of-constants representation of u ε relative to the semigroup (e t∆ ) t>0 in combination with wellknown smoothness estimates (cf.[38,Lemma 1.3]) of said semigroup, the maximum principle and the Hölder inequality to gain constants C, λ > 0, α ∈ (0, 1) and q ∈ (n, ∞) such that As the remaining integral term above is bounded independent of t by prior arguments, taking the supremum over t ∈ [0, T ] then immediately yields an L ∞ (Ω) bound for u ε on [0, T ] for all T ∈ [0, T max ), which is in fact independent of T .This implies that u ε L ∞ (Ω×(0,Tmax)) < ∞ and therefore completes the proof.
Given that we have now established the necessary existence theory for our approximate solutions, let us fix some functions, measures and parameters for the remainder of this paper in a effort to not unnecessarily clutter later results and arguments.
First, let n ∈ N, χ ∈ R be fixed and let Ω ⊆ R n always be a bounded domain with smooth boundary.We then fix some initial data µ ∈ M + (Ω) and a family (u 0,ε ) ε∈(0,1) ⊆ C ∞ (Ω) of nonnegative functions with that approximate µ in the following way: Let us give a brief argument as to how such an approximation of Radon measures can be achieved: It is fairly easy to see that approximating Dirac measures δ x by smooth functions in this way is indeed possible given sufficient boundary regularity (e.g. by using f ε (y) := C(ε)e − 1 ε |x−y| 2 with C(ε) > 0 some normalization constant).This implies that the Dirac measures are contained in the closure (relative to the vague topology) of the set Further, one can show that the Dirac measures are the extreme points of the convex set of probability measures M 1 (Ω) ⊆ M + (Ω), which is compact in the vague topology.This makes it accessible to the Krein-Millman theorem (cf.[32,Theorem 3.23]) implying that and therefore that M 1 (Ω) = F (see also [1,Corollary 30.5]).As M + (Ω) is metrizable (cf.[1, Theorem 31.5]), this is sufficient to gain our desired approximation after a straightforward scaling argument.
If µ = f for some f ∈ L p (Ω), p > 1, we further assume that This additional approximation property can be achieved by standard methods for approximating L p (Ω) functions by smooth functions combined with a straightforward normalization argument to ensure (2.3).

An initial data independent estimate for an ordinary differential equation
Deriving ordinary differential equations for key norms by testing partial differential equations with carefully chosen functions is often one of the first steps in the process of gaining sufficient a priori information about said partial differential equations.But, if the decay terms in these ordinary differential equations are not sufficiently strong, the estimates gained from them are often still dependent on the initial data.In general, this poses not much of a problem as the initial data is in many cases assumed to be fairly regular, but in our case this means estimates of this type are mostly useless because our set of approximate initial data (u ε,0 ) ε∈(0,1) is in general not even bounded in any L p (Ω) with p > 1.As such, the ordinary differential equations we will derive in this paper for norms of our approximate solutions will need to have decay terms strong enough to allow for initial data (and therefore ε) independent estimates.Note that, by their very nature, these estimates will always break down for t ց 0, but may still yield certain integrability properties up to zero instead.
For this purpose, we will in this section consider the initial value problem with y 0 ≥ 0, A > 0, B ≥ 0, α > 1 and prove that the above superlinear decay term is in fact enough to grant us initial data independent estimates for its solution.We will even quantify somewhat how severely these estimates deteriorate for t ց 0.
While the proof presented here is fairly straightforward, we will still present the argument leading to the following estimate for (3.1) in full due to its important role for the central results of this paper.
That the solution y exists locally and is unique is ensured by the Picard-Lindelöf theorem.Global existence and nonnegativity then follow by comparing with a sufficiently large constant or zero. If As this is already sufficient for our desired result, we will now focus on the remaining case y 0 > (B/A) 1 α .In this case, we immediately gain for all t > 0 by essentially the same comparison argument.Given this, we define z(t) := y(t) − (B/A) 1 α > 0 for all t ≥ 0 and then note that By now comparing z with the explicit solution to the initial value problem for all t > 0 and therefore that y(t) ≤ (A(α − 1)) for all t > 0. (3.3) As we have now covered all necessary cases, combining (3.2) and (3.3) completes the proof with .

Deriving our central differential inequality for Ω u p ε
This next section is now devoted to deriving exactly the type of differential inequalities discussed in the previous one for terms of the form Ω u p ε , p ∈ (1, ∞), with constants independent of ε.The basic approach for this is testing the first equation in (2.1) with u p−1 ε , employing partial integration and applying the second equation in (2.1) to the resulting ∆v ε terms.Due to the change of sign of the terms originating from the second equation depending on the sign of χ (resulting in different problematic terms), we will need to treat the repulsive case (χ < 0, see Lemma 4.3) and attractive case (χ > 0, see Lemma 4.5) somewhat separately.Nonetheless, we still manage to achieve essentially the same result for both barring some additional restrictions on the initial data or dimension.
Before we dive into the actual derivation of the central lemmas of this section, we will first establish some preliminary results about the approximate solutions.The first such result is the derivation of some bounds for the second solution component v ε .For this, we employ a classic result from elliptic regularity theory dealing with L 1 (Ω) source terms (cf.[4]) to the second equation in (2.1) to gain baseline bounds for later interpolation arguments.Further, we also test the same equation with v r−1 ε to derive some additional estimates, which will prove helpful when applying some other regularity results later on.
Proof.Given the mass conservation property (2.2), we can apply elliptic regularity theory that can be used with L 1 (Ω) source terms (cf.[4, Lemma 23]) to the operator −∆ + 1 to gain the desired uniform W 1,p (Ω) bounds for the second solution component for all p ∈ [1, n n−1 ).The remaining L q (Ω) bounds then follow from the Sobolev embedding theorem (cf.[6, Theorem 2.72]).
For r ∈ [1, ∞), the inequality (4.1) is a consequence of testing the second equation in (2.1) with v r−1 ε and integrating by parts to first gain ) r due to Young's inequality and therefore (4.1) for all t > 0 and ε ∈ (0, 1).The case r = ∞ then follows by taking the limit r ր ∞ in (4.1).
As our second preliminary result of this section, we now further derive two interpolation inequalities for the first solution component u ε based on a slightly extended variant of the Gagliardo-Nirenberg inequality found in [16], which are used in both the repulsive and attractive case.and for all t > 0 and ε ∈ (0, 1) with m as defined in (2.3).
we can use the Gagliardo-Nirenberg inequality, or rather a variant of it from [16, Lemma 2.3], which allows for some of the parameters to be from the interval (0, 1) in a way not covered by the original inequality, to gain for all t > 0 and ε ∈ (0, 1) with with K 2 := (2K 1 ) 2+ 4 np due to the mass conservation property seen in Lemma 2.1.
Because moreover the same Gagliardo-Nirenberg type inequality from reference [16] is again applicable and gives us for all t > 0 and ε ∈ (0, 1) with . Having established all the necessary preliminaries, we will now begin deriving the core results of this section by first considering the repulsive case (χ < 0).

Remark 4.4.
Let us now briefly discuss why the above argument breaks down for dimensions greater than three.To do this, we start by identifying its linchpin, namely the question whether or not the v ε integral term in (4.9) can be absorbed by the only term with negative sign in said same inequality.The argument employed by us to answer this question positively in two and three dimensions centrally relies on interpolation and embedding properties of certain Sobolev spaces combined with some elliptic regularity theory to gain an estimate of the form for all p ∈ (1, ∞), ε ∈ (0, 1) and t > 0 with some α < 1 + 2 n(p−1) and C > 0 (see the inequality (4.11) above).But as for higher dimensions the exponent of v ε increases and the type of interpolation properties used diminish in effectiveness, using the same argument in dimension four and higher only yields the above estimate with α > 1 + 2 n(p−1) , which is insufficient for our approach.
While generally the more difficult case to handle and thus needing some additional restrictions on the initial data and dimension to make work, the attractive case (χ > 0) can be handled in a very similar manner to the repulsive case discussed above.The main difference is that the key problematic term changes from due to χ's change of sign.While this reduces proof length as we do not first need to split up the new problematic term to separate u ε and v ε , the Ω u p+1 ε term itself is already much more difficult to handle and can therefore only be compensated for by the dissipative term in two dimensions and for small initial mass, at least when using our approach based on the Gagliardo-Nirenberg inequality.
Proof.Fix p ∈ (1, ∞).Similar to the argument in Lemma 4.3, we start again by testing the first equation in (2.1) with u p−1 ε to gain for all t > 0 and ε ∈ (0, 1).Due to the interpolation result (4.2) from Lemma 4.2, we further know that there exists for all t > 0 and ε ∈ (0, 1) with for all t > 0 and ε ∈ (0, 1) because m ≤ C 1 .Applying the interpolation inequality (4.3) from Lemma 4.2 to the above then immediately gives us our desired result.
Given that we will in two dimensions actually only need the inequality (4.12) to be true for the values p = 5 2 and p = 8 (cf.Lemma 5.1, Lemma 7.1), we can now fix the relevant minimal constant C m > 0 for these two cases.Note that this is exactly the constant C m mentioned in Theorem 1.1.

Corollary 4.6.
There exists a minimal C m > 0 such that, if χ > 0, n = 2 and m ≤ C m , the inequality (4.12) holds for p ∈ { 5 2 , 8} and all t > 0, ε ∈ (0, 1).As the differential inequalities (4.6) and (4.12) are essential for all further arguments (either to ensure sufficient regularity from a time t 0 > 0 onward or to ensure uniform continuity in t = 0), we will from now on always assume one of (S1)-(S3) to be true (with the constant C m in (S2) chosen exactly as in the corollary above).We do this, as this means that either Lemma 4.3 or Corollary 4.6 will be always usable from now on.

Initial data independent a priori estimates for u ε and v ε on
Ω × (t 0 , ∞) for all t 0 > 0 Given that we have now established the central differential inequalities (4.6) and (4.12) in all of the cases (S1)-( S3) and for all relevant values of p used in this section, we will now prove certain u 0,ε -independent smoothing properties of the approximate solutions.That is, we will essentially show that, from each t 0 > 0 onward, both solution components are bounded in sufficiently good function spaces independent of the initial data u 0,ε .
It will be these bounds combined with the compact embedding properties of said function spaces that will then allow us in the following section to gain a null sequence (ε j ) j∈N , along of which the approximate solutions converge to a tuple of functions (u, v), which serves as a candidate for our actual solution.
As the first step in the bootstrap argument leading us toward this goal, we begin by extracting the following L ∞ (Ω) boundedness property from Lemma 4.3 or Corollary 4.6.
We can now use this result to gain a uniform, global C 1+α (Ω)-type bound for v ε due to standard elliptic regularity theory and embedding properties of Sobolev spaces into Hölder spaces.This in turn then allows us to apply parabolic regularity theory from [30] to achieve C α, α 2 (Ω × [t 0 , t 1 ])-type bounds for u ε .
The inequality (5.2) is a straightforward consequence of standard elliptic regularity theory (cf.After this straightforward application of elliptic regularity theory, we will now transition to its parabolic counterpart found in [30,Theorem 1.3] to gain (5.1).For this, we first fix K 1 > 0 such that for all t ≥ 1 2 t 0 and ε ∈ (0, 1) according to Lemma 5.1 and (5.2).In the notation of reference [30], the first equation of (2.1) considered in isolation can be then be written as 3), it is easy to see that a and b have all necessary structure conditions with constants and parameters only depending on t 0 and K 1 .An application of Theorem 1.3 from [30] on each interval [t − 1 2 t 0 , t + 1], t ≥ t 0 , therefore yields constants α ∈ (0, 1) and K 2 > 0, only dependent on t 0 and K 1 , such that for all ε ∈ (0, 1), x, y ∈ Ω and t, s for all ε ∈ (0, 1), x, y ∈ Ω and t, s ∈ [t 0 , ∞) with |t − s| ≥ 1.Combined with (5.3), the inequalities (5.4) and (5.5) yield our desired result.
Given that the critical source term in the second equation in (2.1) is in fact not u ε but uε 1+εuε , we now derive the following corollary translating the above results to said source term.
Proof.This follows directly from (5.1) in Lemma 5.2 due to the fact that While the Hölder bounds derived above are already quite strong and something similar to the following lemma could likely be derived from them by more abstract means, we will now give a very short argument, which provides us with a crucial bound for the gradients of the functions u ε .The argument is based on a straightforward testing procedure for the first equation in (2.1) with u ε .We do this to later be able to translate certain weak solution properties from the approximate solutions to the solution candidates constructed in the following section.
Proof.Fix t 0 > 0. According to Lemma 5.2, there exists K > 0 such that for all t ≥ t 0 and ε ∈ (0, 1).Testing the first equation in (2.1) with u ε then directly yields for all t > t 0 and ε ∈ (0, 1).A straightforward time integration of the above combined with (5.6) completes the proof.
As the last important a priori bound of this section, we again use standard elliptic regularity theory to provide a stronger uniform Hölder bound for the approximate solution components v ε .This will allow us to later use the fact that all approximate solution components v ε are classical solutions of the second equation in (2.1) to argue that their limits are in fact classical solutions of the second equation in (1.1).One key idea here is to apply said elliptic regularity theory not only to the functions v ε themselves but also to the difference functions v ε (•, t) − v ε (•, s) and use the already established parabolic Hölder bounds for u ε to gain time Hölder bounds from regularity theory that is originally only interested in the space variable.

Construction of our solution candidates (u, v) as limits of the approximate solutions
The target of this section will be to construct solution candidates (u, v), which solve (1.1) and adhere to the boundary conditions (1.2), by using the a priori estimates from the previous section.As said estimates are all gained in a fashion that was initial data independent by necessity, the construction in this section will by itself not provide us with any type of continuity in t = 0.This issue will instead be addressed in the section directly following.
We begin by using the compact embedding properties of various function spaces combined with the bounds derived in the previous section to construct our solution candidates as limits of the approximate solutions by multiple subsequence extraction and diagonal sequence arguments.
Proof.Due to Lemma 5.2, there exist β k ∈ (0, 1) and Due to the compact embedding properties of Hölder spaces this implies that, for each k ∈ N, there exist α k ∈ (0, 1), functions ) j∈N is a subsequence of (ε These sequences are gained by multiple subsequence extraction arguments.Due to the uniqueness of pointwise limits, this implies that Therefore, we can define u as u(x, t) := u k (x, t) for all (x, t) ∈ Ω × (0, ∞) and some k ∈ N such that t ∈ [ 1 k , k] (6.4) in well-defined fashion.If we now set ε j := ε (j) j , j ∈ N, we gain our first desired convergence property (6.1) using a standard diagonal sequence argument.
By essentially the same argument as above, but this time derived from the bound established in Lemma 5.5, we can now construct a function v : Ω × (0, ∞) → [0, ∞) and extract another subsequence from the previous one, along of which the convergence property (6.3) holds.
Finally, Lemma 5.1 and Lemma 5.4 allow us to extract a final subsequence with the property (6.2) by using compactness properties (relative to the weak topology) of bounded sets in the spaces L 2 (( 1 k , k); W 1,2 (Ω)), k ∈ N, combined with a similar diagonal sequence argument.Note hereby that the convergence property (6.1) derived above ensures that all these weak limits coincide with the limit function u from (6.4).
These convergence properties are now immediately good enough to ensure that v is a C 2,0 (Ω × (0, ∞)) classical solution of the second equation in (1.1) with Neumann boundary conditions, but are insufficient to derive the same solution property for u and the first equation in (1.1).Instead, we will first show that u is a weak solution compatible with well-known parabolic theory from [13] and then use said theory to argue that it was already classical.This is done by using the fact that such weak solutions are indeed unique and the fact that due to the high regularity of v, the associated classical problem always has a sufficiently smooth solution.
We chose this approach instead of first providing stronger a priori estimates with similar parabolic regularity theory for the approximate solutions, which would yield stronger convergence properties and therefore that u is immediately a classical solution, because it is rather more involved to gain the necessary uniform bounds in an ε independent fashion from said parabolic theory as opposed to the mere fact that the weak solutions considered here are indeed classical.Lemma 6.2.The functions u, v constructed in Lemma 6.1 have the properties and satisfy (1.1) on Ω × (0, ∞) and the boundary conditions (1.2).
Proof.Let u, v be the functions and (ε j ) j∈N be the sequence constructed in Lemma 6.1.
That v ∈ C 2,0 (Ω × (0, ∞)) is immediately obvious from (6.3).Further, all v ε satisfy the second equation in (2.1), which is apart from one term very similar to the second equation in (1.1), and have Neumann boundary conditions.Therefore, the same convergence property directly allows us to conclude that v in fact satisfies the second equation in (1.1) on Ω × (0, ∞) and has Neumann boundary conditions due to the convergence property (6.1) ensuring pointwise convergence of u ε to u on Ω × (0, ∞) and the fact that x 1+εx → x as ε ց 0 for all x ∈ [0, ∞).For the function u, we first argue that it is in fact a weak solution of the first equation in (1.1) on Ω × (t 0 , ∞) for all t 0 > 0 in the sense that u with compact support and t 0 > 0. For this purpose, we now fix t 0 > 0 and then notice that the convergence properties (6.1) and (6.2) already ensure that u Ω)).We observe further that the approximate solutions u ε already are weak solutions of the above type as is easily seen by partial integration.Thus, we now only need to prove that this solution property survives taking the limit ε = ε j ց 0. For the first, second and fourth integral term in (6.5), this is immediately ensured by the convergence properties (6.1) and (6.3).The remaining third integral term converges as desired due to (6.2).
As our final step of this proof, we will now rely on some well-known parabolic regularity results from [13] to show that the weak solution formulation above combined with the already known regularity properties for u and v already imply that u is a classical solution to the corresponding partial differential equation.As (6.3) already implies that the coefficients of the operator L(u, ∇u, ∆u) = −∆u + χ∇u • ∇v + χu∆v are in fact elements of C α, α 2 (Ω × [t 0 , t 1 ]) for all t 1 > t 0 and some α ≡ α(t 0 , t 1 ) ∈ (0, 1), a combination of Theorem 5.1 from [13, p.170], which ensures uniqueness of such weak solutions, Theorem 5.3 from [13, p.320], which ensures the existence of higher regularity classical solutions to this problem, and a standard cut-off argument, which helps deal with the missing initial data regularity, then grants us the remaining desired results.

M + (Ω)-valued continuity of u in t = 0
After having now constructed our solution candidate (u, v) in Lemma 6.1 and shown that it is already a classical solution to (1.1) on Ω × (0, ∞) with Neumann boundary conditions, we now only need to argue that said solution candidate is related to the initial data in a sensible way (as it would be very easy to construct solutions that are not, e.g.certain constant ones).
Because our initial data are in many cases only measure-valued, the way we want to do this is to show that u converges to the initial data µ in the vague topology on M + (Ω) as t ց 0, meaning that it is continuous in t = 0 in said topology with a value of µ at t = 0. We will do this by showing that the first components of the approximate solutions in a sense already are uniformly continuous in t = 0 in the vague topology, which can then be easily translated to a similar continuity property for the limit function u.We do this by proving that terms of the form Ω u ε (•, t)ϕ − Ω u ε (•, 0)ϕ = t 0 Ω u εt ϕ tend uniformly to 0 as t ց 0. For this, we first prepare a lemma, which proves a similar uniform convergence result for the critical term occurring in our later arguments treating the t 0 Ω u εt ϕ term.Note further that it is here, where the last remaining unused assumption in (S1)-(S3), namely the higher initial data regularity in the three dimensional case, comes into play.In three dimensions, it is necessary to have a better than L 3 2 (Ω)-type uniform bound for ∇v ε , which is just about not provided by Lemma 4.1, but available to us for slightly better initial data regularity.Everything before this point was indeed possible for measure-valued initial data.Lemma 7.1.For each δ > 0, there exists t 0 ≡ t 0 (δ) > 0 such that t 0 u ε (•, s)∇v ε (•, s) L 1 (Ω) ds ≤ δ for all t ∈ (0, t 0 ) and ε ∈ (0, 1).
Proof.Let us first treat the simpler case of n = 2. Then a combination of the Hölder inequality and Lemma 4.1 gives us a constant K 1 > 0 such that , for all t > 0 and ε ∈ (0, 1) because 5 3 ≤ 2 = n n−1 .By use of Lemma 3.1, Lemma 4.3 (for the case χ < 0) and Corollary 4.6 (for the case χ > 0), this can be improved to for all t > 0, ε ∈ (0, 1) and K 2 > 0 given by the prior results.As time integration then yields for all t > 0 and ε ∈ (0, 1), this already implies our desired result by setting t 0 := min( δ ).We now treat the more challenging case of n = 3.Given that we are then in scenario (S3), we can assume that µ = f with f ∈ L p (Ω) for some p > 1 and therefore there exists K 3 > 0 such that u 0,ε L p (Ω) ≤ K 3 for all ε ∈ (0, 1) according to (2.5).Due to standard embedding properties of Lebesgue spaces we can assume p < 3 without loss of generality.Integration of (4.6) from Lemma 4.3 then immediately yields K 4 > 0 such that u ε (•, t) L p (Ω) ≤ K 4 for all t ∈ (0, 1) and ε > 0, which due to standard elliptic regularity theory (cf.[7,Lemma 19.1]) combined with (4.1) and Sobolev embedding theorems (cf.[6]) further gives us K 5 > 0 and r ∈ ( 3 2 , 3p 3−p ) such that ∇v ε (•, t) L r (Ω) ≤ K 5 for all t ∈ (0, 1) and ε > 0. By a similar application of the Hölder inequality as in the two-dimensional case, we then gain that with q = r r−1 < 3 for all t ∈ (0, 1) and ε ∈ (0, 1).Again by application of Lemma 4.3 combined with Lemma 3.1, we gain K 6 > 0 such that for all t ∈ (0, 1) and ε ∈ (0, 1).Because this implies that 1 q the estimate (7.1) is sufficient to complete the proof by a similar time integration argument as used in the two-dimensional case.
Given this, we can transition to the analysis of the critical t 0 Ω u εt ϕ term, which is now pretty straightforward.
Proof.Let ϕ ∈ C 2 (Ω) with ∇ϕ • ν = 0 on ∂Ω.We then test the first equation in (2.1) with ϕ and use partial integration to see that Due to Lemma 7.1, this already implies our desired result.
Finally, due to already proving sufficiently strong convergence properties for the approximate solutions in Lemma 6.1, we now only need to show that the uniform continuity in t = 0 hinted at in the previous lemma translates to our solution candidates as proper continuity in t = 0 in the vague topology.We do this as follows: for all ε ∈ (0, 1) and t ∈ (0, t 0 ).Further by (2.4), there exists ε ′ ∈ (0, 1) such that for all ε ∈ (0, ε ′ ).Due to the convergence property (6.1), there moreover exists ε(t) ∈ (0, ε ′ ) for each t ∈ (0, t 0 ) such that This then implies that Combining the final result of the two previous sections now immediately gives us the central result of this paper: Proof of Theorem 1.1.Let (u, v) be as in Lemma 6.1.Then Lemma 6.2 and Lemma 7.3 imply all the desired solution properties of (u, v) when applied in combination.

1 α 1 α
, then by comparison with a constant function of value (B/A) , we immediately know that