EVENTUAL SMOOTHNESS AND ASYMPTOTIC BEHAVIOUR OF SOLUTIONS TO A CHEMOTAXIS SYSTEM PERTURBED BY A LOGISTIC GROWTH

. In this paper we study the chemotaxis-system (cid:40) deﬁned in a convex smooth and bounded domain Ω of R n , n ≥ 1, with χ > 0 and endowed with homogeneous Neumann boundary conditions. The source g behaves similarly to the logistic function and satisﬁes g ( s ) ≤ a − bs α , for s ≥ 0, with a ≥ 0, b > 0 and α > 1. Continuing the research initiated in [33], where for appropriate 1 < p < α < 2 and ( u 0 ,v 0 ) ∈ C 0 (¯Ω) × C 2 (¯Ω) the global existence of very weak solutions ( u,v ) to the system (for any n ≥ 1) is shown, we principally study boundedness and regularity of these solutions after some time. More precisely, when n = 3, we establish that - for all τ > 0 an upper bound for ab , || u 0 || L 1 (Ω) , || v 0 || W 2 ,α (Ω) can be pre-scribed in a such a way that ( u,v ) is bounded and H¨older continuous beyond τ ; - for all ( u 0 ,v 0 ), and suﬃciently small ratio ab , there exists a T > 0 such that ( u,v ) is bounded and H¨older continuous beyond T . Finally, we illustrate the range of dynamics present within the chemotaxis system in one, two and three dimensions by means of numerical simulations.


Introduction and motivations.
According to the logistic model in population dynamics (Pierre-François Verhulst, 1838), the self-limiting growth of a biological u t = ∇ · ( ∇u − χu∇v), v t = ∆v − κv + u. (1) In system (1), the parameters , χ and κ are positive constants. Although we will consider such cross diffusive terms in terms of chemotactic movement, they are able to appear for a variety of reasons, including growth of the underlying solution space [14,43,44,5] and inhomogeneities in the underlying environment [4].
In the present case we prescribe the following description to the equations. Chemoattractant, v, spreads diffusively, decays with rate κ and is also produced by the bacteria with rate 1. The bacteria diffuse with mobility and drift in the direction of the gradient of concentration of the chemoattractant with velocity χ|∇v|; χ is called chemosensitivity. Hence, once the initial cells distribution and chemical concentration (that is u 0 (x) = u(x, 0) and v 0 (x) = v(x, 0)) are given, under zero-flux boundary conditions on both u and v, the previous problem describes the chemotactic dynamics of a cells population in a totally insulated domain.
Real observations show that this movement may eventually lead to aggregation processes, in which the density of the cells spatially coalesces and grows without bound (chemotactic collapse). Mathematically, this collapse implies that (possibly) u becomes unbounded at one or more points of its domain at a certain instant (blowup time). It is known that, in a one-dimensional domain, all the solutions of (1) are global and uniformly bounded in time (see [22]), while that in the n-dimensional setting, with n ≥ 2, unbounded solutions to the same problem have been detected (see, for instance, [8] and [39]).
In line with the chemotactic scenario, in [10] for radial solutions and in [20] for non-radial, the authors prove that under suitable assumptions the bacteria concentration blows up in finite time, for certain domains of R 2 and in the cases in which the second differential equation of (1) is replaced by 0 = ∆v − v + u (parabolicelliptic case). Moreover, for the classical parabolic-parabolic (or fully parabolic) case, estimates from below and numerical computations for the blow-up time of unbounded solutions to (1) are derived in [24] and [7], respectively (see also [16] for a more general analysis).
Furthermore, a number of interesting results concerning properties of solutions to chemotaxis-systems have been also attained for a broader class of problems, in which the first equation of (1) reads u t = ∇ · (S(u)∇u) − ∇ · (T (u)∇v). Precisely, bounded or unbounded solutions of the corresponding problem is determined by the asymptotic behaviour of the ratio T (u)/S(u), especially in terms of the space dimension; we refer, for instance, to [6] and [40] for the parabolic-elliptic case and to [9,17,27,28,38] for the parabolic-parabolic case.
As an approach towards the model of self-organizing behaviour of cells populations, it seems coherent to adapt the original Keller-Segel formulation to the case in which the temporal evolution of a cells distribution may be perturbed by the proliferation and the death of the cells themselves. Mathematically it is possible by adding a linear combination of power functions depending on u and, possibly, on |∇u| to the first equation of system (1) (see details in [1,18,32] for pure chemotaxissystems, but also in [31] for weakly coupled systems).
Conforming to the previous paragraph, this investigation focuses on fully parabolic chemotaxis-systems which are complemented by logistic-type effects. To the best of our knowledge, the following are the most recent and partial results in this regard; under Neumann boundary conditions and in a convex smooth and bounded domain Ω of R n , n ≥ 1: i) For the problem the existence of global weak solutions is proven for any nonnegative and sufficient regular initial data (u 0 , v 0 ) and arbitrarily small values of b > 0. Moreover, if n = 3 and a is not too large, these solutions become classical after some time (see [13]). ii) For the problem where g generalizes the logistic source in (2), and satisfies g(0) ≥ 0 and g(s) ≤ a − bs 2 , for s ≥ 0, and with a ≥ 0, b, χ, τ positive constants, it is proved in [37] that if b is big enough, for all sufficiently smooth and nonnegative initial data, (u 0 , v 0 ), the problem possesses a unique bounded and global-in-time classical solution. Furthermore, even though [21] yields global classical solutions to (3) for any (not necessarily large) b > 0, which remain bounded in a (convex smooth and bounded) domain of R 2 , the same conclusion is not clear to occur for n ≥ 3. iii) For the same problem (3), but with source term g controlled, respectively from below and above, by −c 0 (s + s α ) and a − bs α , for s ≥ 0, and with some α > 1, a ≥ 0 and b, c 0 > 0, global existence of very weak solutions is attained in [33]. Moreover, for n = 3, sufficient conditions on the initial data and the coefficients of the source g which ensure the boundedness of such solutions are discussed in [34]. Additionally, analogous conclusions dealing with parabolic-elliptic versions of models related to (2) or (3) are also available. For instance, in [29] it is proven that weak solutions exist for arbitrary b > 0; moreover they are smooth and globally classical if b > (n − 2)/n. Finally, with source term g as in the above item iii), global existence of very weak solutions and their boundedness and eventual smoothness properties are established in [35].
Our starting point are the contributions [33] and [34], where these partials results are, respectively, presented: 1 ) existence of global very weak solutions: for any n ≥ 1 and α ∈ (1, 2) satisfying α > 2 − 1 n , the global existence of very weak solutions (u, v) to the system is shown for any nonnegative initial data, (u 0 , v 0 ) ∈ C 0 (Ω) × C 2 (Ω), and under zero-flux boundary condition on v 0 ; 2 ) boundedness of very weak solutions: in the most realistic three dimensional setting, these very weak solutions derived in 1 ) are bounded. More precisely, if the ratio a/b does not exceed a certain value and the initial data are such that ||u 0 || L p (Ω) and ∇v 0 L 4 (Ω) are small enough then, for appropriate 9/5 < p < α < 2, (u, v) is uniformly-in-time bounded in (0, ∞). A natural and complementary question connected to the above results 1 ) and 2 ) is to show that singularities of solutions to (4), possibly arising after some finite time, disappear. In other words, our main objective is to analyse whether solutions to (4), with "unclear" behaviour over a certain period, become eventually bounded and smooth beyond some time; additionally, it is also interesting to investigate and characterize their long time behaviour. Specifically, we address the following issues, which are directly related to each other.
-For any fixed time τ > 0, is it possible to claim that the very weak solutions of (4) improve their regularity and become bounded beyond such τ ? According to Theorem 2.1 below, by imposing suitable smallness conditions on u 0 and v 0 , measured in proper norms, the very weak solutions are eventually bounded and Hölder continuous provided the ratio a/b does not exceed a certain value.
-Is it possible to claim that any of the very weak solutions of (4) improve their regularity and become bounded beyond some time, regardless the initial sizes of u 0 and v 0 ? Again under smallness assumption on the ratio a/b, this question is positively shown in the forthcoming Theorem 2.2. It is proved that it is always possible to find a T > 0 such that the very weak solution are bounded and Hölder continuous beyond T , independently by some norm of the initial data u 0 and v 0 ; clearly, the time beyond which it occurs depends on such initial norm. -Is it possible to characterize the asymptotics of bounded solutions, i.e. their behaviour for t → ∞? If on the one hand all the solutions of the Pierre-François Verhulst model converge to the constant steady state K (as shown in §1), then the chemotaxis-diffusion-growth models may lead to a spatially uniform steady state, or to a spatially heterogeneous steady state, as well as irregular spatiotemporal solutions, possibly defined by time-periodic or timeirregular pattern formations (see [23,2]). The theoretical analysis dealing with the behaviour of the solutions to (4) when the time increases is currently unclear and goes beyond the scope of this paper; here, we present an important number of numerical simulations (in one, two and three dimensions) concerning the long time behaviour of bounded solutions and also the blow-up scenario of unbounded ones (see §6).
3. Preliminaries and definition of suitable solution. The following results are herein formulated according to our purposes, and they are used through the paper to prove the claimed theorems. For the sake of clarity, we also close the section with the definition of very weak solutions to (4).
Then, for any concave function ϕ : R → R this inequality holds 1 in Ω ⊂ R n , n ≥ 1.
Then the operator (−∆ + 1) is sectorial in L p (Ω) and for any ρ ≥ 0 possesses fractional powers (−∆ + 1) ρ , with dense domain D((−∆ + 1) ρ ). Moreover, there exist positive constants C S and µ 1 such that -for all t > 0 and p ≤ q < ∞ the following L p -L q estimates hold -for all t > 0 and 1 < p ≤ q < ∞ the operator e t∆ ∇· possesses a uniquely determined extension to an operator from L p (Ω) to L q (Ω) obeying this L p -L q estimate Proof. See Section 2. of [9] and Lemma 1.3 of [36].
We also make use of these elementary results.
Lemma 3.3. Let y be a positive real number satisfying y ≤ c(y l + 1) for some c > 0 and 0 < l < 1. Then y ≤ max{1, (2c) Proof. Let us suppose y ≥ 1 (if y ≤ 1 there is nothing left to show): y l ≥ 1 so that y ≤ c(y l + 1) ≤ 2cy l and hence y ≤ (2c) Lemma 3.4. For any A, B ≥ 0 and p > 1 we have Proof. Let us prove relation (12) for both A and B strictly positive; otherwise the inequality is obvious. If Finally, we give these definitions.
is said to be a very weak subsolution of (4) in Ω × (0, T ) if g(u) and u∇v belong to is said to be a weak γ-entropy supersolution of (4) in for all nonnegative ψ ∈ C ∞ 0 (Ω × [0, T )). Definition 3.7. Let T > 0. A pair (u, v) of functions is called very weak solution for problem (4) in Ω × (0, T ) if it is both a very weak subsolution and a weak γ-entropy supersolution of (4) in Ω × (0, T ), in the sense of Definitions 3.5 and 3.6.
A global very weak solution of (4) is a pair (u, v) of functions defined in Ω×(0, ∞) which is a very weak solution of (4) in Ω × (0, T ) for all T > 0. 4. Approximate problem and existence of very weak solutions. In preparation for the main estimates, by means of a parameter ε ∈ (0, 1), we define a perturbed chemotaxis-system, properly constructed to approximate and hence solve and analyze the original problem (4). Precisely, we rely on these specific results.
We also present this result concerning both the existence and boundedness properties of global very weak solutions to system (4).

Proposition 2.
Let Ω be a convex smooth and bounded domain of R n , with n ≥ 1. For some α ∈ (1, 2), with α > 2 − 1 n , let us assume g ∈ C 1 ([0, ∞)), such that g(0) ≥ 0, and verifying both assumptions (H1 α ) and (H2 α ). Then for all nonnegative functions (u 0 (x), v 0 (x)) ∈ C 0 (Ω)×C 2 (Ω), with ∂v0 ∂ν = 0 on ∂Ω, problem (4) admits at least one global very weak solution (u, v), according to Definition 3.7. More precisely, this solution is the limit of a sequence of globally bounded couples of functions (u ε , v ε ) which classically solve the approximate problem (13), in the sense that, for any t > 0 and 1 < p < α < 2, as ε → 0 we have Additionally, for n = 3, there exists a positive real δ > 0 such that if a and b are such that a/b < δ then for all 9/5 < p < α < 2 it is possible to find a λ > 0 with the following property: Proof. For the existence question, see the proof of Theorem 2.1 in [33]. For the boundedness, Theorem 2.1 in [34]. Remark 1. By taking into account the regularity properties that the limit (u, v) provided by Proposition 2 inherits from the sequence (u ε , v ε ), we have that for n = 1 the Sobolev inequalities imply that both W 1,2 (Ω) and W 1, α α−1 (Ω) are compactly embedded in C(Ω). As a consequence, in the one dimensional setting, for any t > 0 the solution (u(·, t), v(·, t)) to problem (4) belongs to L ∞ (Ω) for any choice of the initial distribution (u 0 , v 0 ) and the parameters a and b which determine the behaviour of function g.

5.
A priori estimates and proof of the theorems. In this section our principal objective is to infer some ε-independent and uniform-in-time estimates for both u ε and v ε components of the solution (u ε , v ε ) to (13). In this sense, the following lemma includes some inequalities which are strongly employed with this aim.

GIUSEPPE VIGLIALORO AND THOMAS E. WOOLLEY
In addition, since p/α < 1, the function t → Ω u α ε p α is concave; an application of the Jensen inequality (8) where we have also used relations (14c) and (12) with p = α/p > 1.
On the other hand, since 9/5 < α < 2, the Sobolev embedding W 2,α → W 1,4 infers a positive constantĈ such that The linear (Cauchy) problem extracted from the second equation of (13), that is , has a unique classical and global solution, so that we can apply relation (12) of Proposition 2 of [33] in the interval [t 1 , t 2 ]. Precisely, in line with the nomenclature used in such a proposition, setting p = q = α > 1, using that for any A, B ≥ 0 the relation (A + B) α ≤ 2 α (A α + B α ) holds and estimating the interpolation norm v ε (·, t 1 ) α,1− 1 α with v ε (·, t 1 ) W 2,α (Ω) , we can find a positive constantC, independent on ε, such that t2 t1 Ω In this way, by integrating (20) between t 1 and t 2 , and in view of (14c), relation where C Ω =ĈC. The sum of the expressions (18) and (22), once the bound (19) is also taken into account, concludes the proof.

EVENTUAL SMOOTHNESS AND ASYMPTOTICS OF SOLUTIONS 3033
Proof. For any t 1 ≥ 0, let us define the function By retracing the proof of Lemma 5.2 in [34], and using the same constants therein introduced, one can show that Φ ε satisfies this absorptive differential inequality for any t > t 1 Through a continuous dependence argument, regardless the size of the positive root of Θ 0 (ξ), it is always possible to find a value M lim > 0 with the property that the equation Θ M lim (ξ) = 0 possesses two positive roots (let us say ξ 1 and ξ 2 , with ξ 1 < ξ 2 ) such that ξ lim = ξ 1 < 1. In this sense S M lim ≡ {ξ lim , ξ 2 }. Hence for any ξ ≥ 0 we have Additionally, Φ ε (t) ≡ ξ lim satisfies this initial problem Successively, for any τ > 0 let us define From now on our aim is to justify the existence of at ∈ (t 1 , t 1 + τ /2) and show that Φ ε (t) satisfies the initial problem Indeed, in view of (26), assumption (23) would imply relations (24), which in turn would allow us to apply an ODE comparison principle to problems (25) and (27) and conclude that Φ ε (t) ≤ ξ min for all t >t, ξ min being some zero of Θ Mε (ξ), with 0 < ξ min < ξ lim < 1.
The average theorem establishes the existence of a timet belonging to (t 1 , t 1 + τ 2 ) such that In this way, an application of (15) from Lemma 5.1 with Through (23) and (26) we derive that so that, in particular, 1 2 Ω |∇v ε (·,t)| 4 α 4 < 1. Subsequently, in view of α/4 < 1, we finally get As claimed, this implies that Φ ε (t) ≤ ξ min for all t >t, for some ξ min = ξ min (τ ) < 1; thereafter, since ||u ε || L p (Ω) ≤ p Now, the next conclusion proves that the boundedness of ||u ε || L p (Ω) is sufficient to show that for some 3 < q < 3p/(p − 3) the W 1,q (Ω)-norm of the component v ε is controlled by a certain ε-independent and uniform-in-time positive constant; more exactly, we show that if u ε ∈ L p (Ω) for some t ≥ t 1 , with t 1 ≥ 0, than v ε belongs to some W 1,∞ (Ω) for t beyond t 1 . We will use also some ideas presented in [3].
In addition, for all t ≥ t 1 + τ * , we have
After these preparations, the proof of our main results consists in demonstrating higher regularity for (u ε , v ε ); precisely, we first apply to the classical solution (u ε , v ε ) of problem (13) regularity results which enable us to derive ε-independent bounds of this solution in some Hölder space. Then, by properly interpreting the two equations of the same (13) as Neumann boundary value problems with Hölder and bounded sources and coefficients, we gain higher regularity for (u ε , v ε ). Subsequently, passing to the limit and organizing the statements, our claims follow.
6. Numerical simulations. In this section we numerically test the presented results by simulating the chemotaxis systems in one, two and three dimensions. Further, we investigate whether the global solutions are bounded and stationary, or whether they have complex temporal dynamics, such as moving peaks, oscillations, or chemotactic blow-up. Specifically, we use finite element methods to simulate on domains of different dimension, employing Neumann conditions on the domain boundaries, where applicable. Further, (unless otherwise stated) the initial conditions are uniformly randomly generated with mean equal to the homogeneous steady state, (a/b) 1/α and with range (a/b) 1/α /100. Domain sizes, domain discretisations and parameter values are given in the captions of each figure. Figure 1 illustrates the one-dimensional evolution of system (45) over time and over multiple values of α. As highlighted in Remark 1 of Proposition 2, we expect that when 1 < α < 2 the solution converges to a bounded solution, which we note is heterogeneous and stationary (see Figure 1(a), and, in particular, compare the figures as times 50 and 100). However, although the theory ensures existence of global (and bounded) solutions only for α > 1, we can see that this is not a bifurcation point since in Figure 1(b) the system also converges to a bounded, heterogeneous, stationary solution when α = 0.9.  In Figure 2 we further illustrate the coherence of Remark 1 in this paper as it suggests that for all 1 < α < 2 the one-dimensional simulations should be bounded no matter the initial conditions on u and no matter the size of b. Figure 2(a) illustrates that even after increasing the average initial value of u to 100 the solution still tends to a stable, heterogeneous, bounded solution. Equally, in Figure 2(b), where the parameter b has been reduced to 0.2 we see that the solution is bounded. However, here we see a new dynamic of continuously evolving peaks. Namely, peaks appear approximately in the centre of the domain and then travel towards the boundaries x = 0, or x = 10. The direction of travel depends on which side of x = 5 the peak first appears. Further, the peaks appear to alternate in directions, with the first peak travelling to the left.  Figure 1(a). However, each simulation involves a single parameter change. Specifically, in (a) a larger initial condition for u was used (100 was added to the mean); in (b) the parameter b was reduced to 0.2; Finally, in (c) the spatial solution domain has been reduced from 10 to 1.
Up until now all of the solutions have presented bounded, spatially heterogeneous solutions. However, stationary, bounded, uniform solutions are also possible and these are illustrated in Figure 2(c). Specifically, using spectral analysis it is well known [41,19] that concentration heterogeneity arises through a symmetry breaking of the balance between diffusion and chemotaxis. Hence, patterning of the solution domain requires a minimal spatial scale before instability of the uniform steady state can occur [42,15]. Since the solution domain in Figure 2(c) has been reduced below this critical spatial scale the system simply tends to a uniform steady state. Figure 3 illustrates the two-dimensional evolution of system (45) for α = 1.6 and 1.1. The solutions are simulated on a circular domain of radius 5. The two simulations contain values of α either side of the critical value of 3/2 (recall again Proposition 2), which states that when α is above the critical value the solutions exist, are global and (possibly) bounded. In the specific case of Figure 3(a) the solutions are, once again, bounded, heterogeneous and stationary. However, for shows the full extent of the peak, which is growing without bound. α = 1.1, which is below this critical value, we are able to produce solutions that are prone to blow up. Note that this difference in convergence is purely a property of domain dimension because, apart for this difference, figures 1(a) and 3(b) have the same parameters values. Further, we note that the exploding solution appears to occur foremost on the boundary and that the blow-up occurs very quickly since the maximum value of the solution is over 10 3 (see inset of Figure 3(b)) in just over 17.4 time points. By t ≈ 17.43 the peak is over 10 12 and the solution fails to converge. Finally, Figure 4 illustrates the three-dimensional evolution of system (45) for α = 1.8 and 1.3, which are values either side of the critical value of 5/3. The solutions are simulated within a solid sphere of radius 7. Once again, we see that Proposition 2 holds in Figure 4(a) as, for suitable initial data, the system evolves to a stationary solution, which is both bounded and heterogeneous. On the other hand, we see that the simulation within Figure 4(b) blows up within 1.44 time units. In both cases, we can see that the balls become darker towards their centre. This means that the population u becomes denser there because the isosurfaces corresponding to darker colours correspond to larger values of u (see caption of Figure 4 for more details). Whilst Figure 4(a) remains finite over the entire region, Figure 4(b) contains a region that begins to undergo chemotactic collapse (see the dark spot in Figure 4(b)). At the illustrated time point the density is over 10 6 and grows to over 10 12 , before the simulation fails to converge. In this case, in contrast to Figure 3(b), the unstable growth occurs within the domain, rather than on the boundary. In summary, these simulations illustrate the veracity of the proof contained within this paper. Specifically, global boundedness of a chemotactic system depends on the spatial dimension we are considering, as well as the kinetic parameters of the system.