Blow-up of solutions to semi-discrete parabolic-elliptic Keller-Segel models

The existence of weak solutions and upper bounds for the blow-up time for time-discrete parabolic-elliptic Keller-Segel models for chemotaxis in the two-dimensional whole space are proved. For various time discretizations, including the implicit Euler, BDF, and Runge-Kutta methods, the same bounds for the blow-up time as in the continuous case are derived by discrete versions of the virial argument. The theoretical results are illustrated by numerical simulations using an upwind finite-element method combined with second-order time discretizations.


Introduction
This paper is concerned with the derivation of upper bounds for the blow-up time in semi-discrete Keller-Segel systems in R 2 . The (Patlak-) Keller-Segel model describes the evolution of the cell density of bacteria or amoebae that are attracted by a chemical substance [23,34]. Under simplifying assumptions, the cell density n(x, t) and the density of the chemo-attractant c(x, t) satisfy the parabolic-elliptic equations (1) ∂ t n = div(∇n − n∇c), −∆c + αc = n in R 2 , where α ≥ 0 measures the degradation rate of the chemical substance. Denoting by B α the Bessel potential if α > 0 and the Newton potential if α = 0 (see the appendix for the definitions), this system can be formulated more compactly as a single equation, (2) ∂ t n = div ∇n − n(∇B α * n) in R 2 with the initial condition n(·, 0) = n 0 in R 2 .
1.1. Blow-up properties. The coupling in equations (1) is a positive feedback: The cells produce a signal that attracts other cells. The cell aggregation is counterbalanced by diffusion, but if the cell density is sufficiently large, the chemical interaction dominates diffusion and may lead to finite-time blow-up of the cell density. System (1) exhibits a dichotomy. Consider first the case without degradation, α = 0. If the initial mass satisfies M := R 2 n 0 dx < 8π, no aggregation takes place and the cells just diffuse for all times [6]. On the other hand, if the mass is supercritical, M > 8π, and the second moment I 0 := R 2 n 0 |x| 2 dx is finite, the solutions blow up in finite time. In the limit case M = 8π, there is blow-up in infinite time with constant second moment [5]. We remark that the critical space is L d/2 (R d ) in space dimensions d ≥ 2 [15]. When the signal degradates, α > 0, a similar criterium holds: The solutions exist for all time for subcritical initial masses M < 8π, and they blow up in finite time for supercritical masses M > 8π and sufficiently small second moment I 0 [9]. The blow-up time T max can be estimated from above by (3) T max ≤ T * : The upper bound for the system with degradation becomes larger than that one without degradation, indicating that blow-up happens later than without degradation, which is biologically reasonable. The idea of the proof is a virial argument: The second moment I(t) = R 2 n(·, t)dx solves the differential inequality If α = 0, this expression becomes an identity. Now, if t > T * , we infer that I(t) < 0, contradicting n(·, t) ≥ 0. Very detailed numerical simulations of the collapse phenomenon in the parabolic-parabolic system have been performed in, e.g., [8]. The asymptotic profile of blow-up solutions in the parabolic-elliptic model was studied in [16]. Numerical blow-up times were computed, for instance, in [18] using a kind of H 2 norm indicator; in [8] using the moving-mesh method; and in [17] using discontinuous Galerkin approximations.
The aim of this paper is to prove the existence of weak solutions to various versions of time-discrete equations and to derive discrete versions of inequality (4). We analyze the implicit Euler and higher-order schemes, including BDF (Backward Differentiation Formula) and Runge-Kutta schemes. Our goal is not to design efficient numerical schemes but to continue our program to "translate" mathematical techniques from continuous to discrete situations [27,28,29].
The proof of inequality (4) is based on three properties: nonnegativity of n(x, t), mass conservation, and a symmetry argument in the nonlocal term in (2). The implicit Euler scheme satisfies all these properties, while the nonnegativity cannot be proved for higherorder schemes, although it seems to hold numerically (see Section 4). We discuss this issue below.
1.2. State of the art. The literature on the analysis and numerical approximation of Keller-Segel models is enormous. Therefore, we review only papers concerned with the analysis of numerical schemes and the blow-up behavior of the discrete solutions, in particular possessing structure-preserving features. We do not claim completeness and refer to the introduction of [1] for more references.
Most of the numerical schemes for the Keller-Segel model utilize the implicit Euler method for the time discretization and aim to preserve some properties of the continuous equations, like (local) mass conservation, positivity (more precisely, nonnegativity) preservation, or energy dissipation. These schemes use (semi-implicit) finite-difference methods [12,31,38]; an upwind finite-element discretization [36]; an Eulerian-Lagrangian scheme based on the characteristics method [39]; a Galerkin method with a diminishing flux limiter [42]; and finite-volume methods [19,45]. A finite-volume scheme was also studied in [1] but with a first-order semi-exponentially fitted time discretization. Finally, a mass-transport steepest descent scheme was analyzed in [4]. All these schemes are based on first-order discretizations.
Only few results are concerned with higher-order time integrators. Strongly A(θ)-stable Runge-Kutta finite-element discretizations were analyzed in [33], and the convergence of the discrete solution was shown. However, mass conservation or positivity preservation was not verified. As A-stable time integrators are computationally very costly, often splitting methods are used. A third-order strong stability preserving (SSP) Runge-Kutta time discretization for the advection term and the second-order Krylov IIF (implicit integration factor) method for the reaction-diffusion term, together with a positivity-preserving discontinuous Galerkin approximation in space, was suggested in [44], but without any analysis.
A number of papers are concerned with the preservation of the nonnegativity of the cell density. An example is [13], where a semi-discrete central-upwind scheme was proposed. Moreover, in [14], a hybrid finite-volume finite-difference method was combined with SSP explicit Runge-Kutta schemes (for the parabolic-parabolic case) or the explicit Euler scheme (for the parabolic-elliptic case). Clearly, a CFL condition is needed to ensure the stability of the explicit schemes. Another approach was used in [43] for a related tumor-angiogenesis model, where a Taylor series method in time allows for higher-order but still explicit schemes.
We remark that there do not exist SSP implicit Runge-Kutta or multistep methods of higher order [21,Section 6]. Moreover, SSP for such methods is guaranteed only under some finite time step condition [7,40] (also see the recent work [24]). In view of these results, positivity preservation of our higher-order schemes cannot be expected.
We do not aim to preserve the free energy of the Keller-Segel system since such schemes usually destroy the symmetry property needed for the blow-up argument; see Remark 3 for details. Some estimates on the discrete free energy in an implicit Euler finite-volume approximation were shown in [45]. A numerical scheme that dissipates the free energy numerically was suggested in [11] using a gradient-flow formulation of the energy functional with respect to a quadratic transportation distance. It is shown in [38] that the dissipation of the discrete free energy may fail in an upwind finite-difference scheme. The dissipation of the discrete entropy in a finite-volume modified Keller-Segel system was proved in [3].
The novelty of this paper is the derivation of the critical initial mass and an upper bound for the blow-up time for various time discretizations, where the values are the same as in the continuous case. This shows that the "continuous" methods carry over to the semi-discrete case. Because of the virial argument, the extension to spatial discretizations is more delicate. A possible direction is to extend the virial argument to a discrete Bessel potential, but we leave details to a future work.

Main results.
Our main results can be sketched as follows (details will be given in sections 2 and 3). Let n k be an approximation of n(·, kτ ), where τ > 0 is the time step and k ∈ N. We recall the definitions M = R 2 n 0 dx of the initial mass and I 0 = R 2 n 0 |x| 2 dx of the initial second moment.
• Existence of solutions to the implicit Euler scheme (Theorem 1): For given n k−1 ≥ 0 and sufficiently small time step τ , there exists a unique weak solutions n k to the semi-discrete equations. Moreover, the scheme preserves the positivity, conserves the mass, and has finite second moment. for the implicit Euler scheme, the solution blows up and kτ ≤ T * . The same result holds for the implicit midpoint and trapezoidal rule if α > 0. The interesting fact is that the upper bound T * is the same for both the continuous and semi-discrete equations. Observe that the existence results do not require the condition M < 8π, since they are local. In particular, the smallness condition on τ is natural, and the time step needs to be chosen smaller and smaller when the blow-up time is approached.
The paper is organized as follows. Section 2 is concerned with the analysis of the implicit Euler scheme, while some higher-order schemes (BDF, Runge-Kutta) are investigated in section 3. In section 4 we provide some numerical examples to illustrate our theoretical statements. Finally, the appendix collects some auxiliary results.

Implicit Euler scheme
We show the existence of weak solutions to the implicit Euler approximation of the Keller-Segel system (2) and their finite-time blow-up, where Here, τ > 0 is the time step and B α is the Bessel potential if α > 0 and the Newton potential if α = 0 (see the appendix). First, we study the existence of solutions. For this, let X : Theorem 1 (Existence for the implicit Euler scheme). Let n k−1 ∈ X and τ < 1 (π + 1/2) 2 1 ( n k−1 X + 1) 4 .
Then there exists a unique weak solution n k ∈ X ∩ H 1 (R 2 ) to (5) such that • nonnegativity: if n k−1 ≥ 0 then n k ≥ 0 in R 2 , • conservation of mass: The strategy is to prove first the existence of a unique very weak solution n k ∈ X to the truncated problem The second step is to show that n k ∈ H 1 (R 2 ). Then we can use n − k = min{0, n k } as a test function in the weak formulation and prove that n k ≥ 0.
Step 1: Solution to (6). To simplify the notation, we write n := n k and n 0 := n k−1 . We introduce the operator ∇c : X → L ∞ (R 2 ) 2 , ∇c[n] = ∇B α * n. We claim that this operator is well-defined and continuous. Indeed, set for α ≥ 0, Then, after the substitution s = |x| 2 /(4t), we can write and it follows that where b := 1 + 1/(2π). This shows the continuity of ∇c.
Next, for given n ∈ X, the linear problem has a unique solution in H 1 (R 2 ). Indeed, since we have n + ∈ L 2 (R 2 ), ∇c[ n] ∈ L ∞ (R 2 ) 2 , and therefore f := n + ∇c[ n] ∈ L 2 (R 2 ) 2 , we can apply Lemma 15 in the appendix yielding the unique solvability of the linear problem. The solution is given by Hence, we can define the fixed-point operator T : X → X by T [ n] := n, and n is given by (10). Clearly, any fixed point of T is a solution to (6). We apply the Banach fixed-point theorem to T on the set S := {n ∈ X : n X ≤ n 0 X + 1}. For this, we need to show that T : S → S is a contraction.
For the proof of T (S) ⊂ S, we use Lemma 13, the Young inequality with p = q = r = 1 (see the appendix), and (9): Combining the previous two estimates, we conclude that To show the contraction property, let n, m ∈ S. Then, estimating as above, Since πτ 1/2 b( n 0 X + 1) < 1, T is a contraction. By the Banach fixed-point theorem, T has a fixed point n k ∈ X, which is a solution to (6).
Step 2: Regularity of the solution to (6). We prove that for any α ≥ 0, n k ∈ H 1 (R 2 ) solves for some positive constant C 0 which depends on τ and n k−1 X . Again, we set n := n k and n 0 := n k−1 . Since n ∈ X ⊂ L 2 (R 2 ), we have n + ∇c[n] ∈ L 2 (R 2 ) 2 . Therefore, by Lemma 15, n ∈ H 1 (R 2 ) is the unique solution to (11). We take n as a test function in that equation and use the Young inequality: By (9) and n X ≤ n 0 X + 1, we find that Then, since β := 1 − b 2 ( n 0 X + 1) 2 τ > 0, we infer that and the claim follows as a test function in (11): since n k−1 n − k ≤ 0 and the last integral on the right-hand side vanishes. This shows that n k ≥ 0 in R 2 , and n + k = n k in (11).
Step 4: Mass conservation. The statement follows immediately if we could use φ(x) = 1 as a test function in (6). Since this function is not integrable, we need to approximate. As in [25], we introduce the radially symmetric cut-off function The following properties hold: for some constants R as a test function in (6) and insert (8): By symmetry and (14), the second integral can be estimated as These estimates allow us to apply the dominated convergence theorem, which leads, in the limit ε → 0, to The same estimates as before show that both integrals on the right-hand side can be estimated by a multiple of 1/R 2 such that the limit R → ∞ leads to which gives mass conservation.
Step 5: Control of the second moment. Similarly as in step 4, we approximate |x| 2 by setting ψ . Using ψ ε R as a test function in (6) and passing to the limit ε → 0 and then R → ∞, it follows that Young's inequality and estimate (9) for ∇c[n k ] = ∇B α * n k show that Therefore, with n k X ≤ n k−1 X + 1, (15) gives Since 1 − τ b( n k−1 X + 1) > 0, we infer that the second moment of n k is bounded if the second moment of n k−1 does so.
Next, we turn to the finite-time blow up of semi-discrete solutions. Set, for α > 0, Theorem 2 (Blow-up for the implicit Euler scheme). Assume that be a sequence of nonnegative weak solutions to (5). Then this sequence is finite with maximal index k max , where, if α = 0, .
In case α > 0, if additionally I 0 ≤ I * and τ < τ * then Proof. Let first α = 0 and let n k be a weak solution to (5) i.e., we assume that (17) does not hold. We set I k = R 2 n k |x| 2 dx. Then, by (15), for any j ≤ k, where we used the conservation of mass and the definition of ∇B 0 * n k . A symmetry argument leads to Summing this identity over j = 1, . . . , k and taking into account the choice of k gives which is a contradiction to n k ≥ 0. Next, let α > 0. For the proof, we follow the lines of [9, Section 6] but the end of the proof is different. Let n k ≥ 0 be a weak solution to (5) such that (18) is not true. Similarly as above, we find that where we recall definition (7) of g α .
We need to estimate 1 − g α (r). For this, let z ∈ R 2 , r = |z| ∈ (0, 1/ √ α). We compute where K = 2π sup |x|<1 |x|B 1 (|x|). It is known that B 1 (x) behaves asymptotically as − log |x| as |x| → 0, so K is finite. A numerical computation shows that sup |x|<1 |x|B 1 (|x|) ≈ 0.0742 and K ≈ 0.4662. We conclude that This bound, together with 1 − g α (|z|) ≤ 1, shows that the last integral in (19) can be estimated as follows: where we have applied the Cauchy-Schwarz inequality in the last step. We infer from (19) that Now, the argument differs from that one used in [9]. Set β = √ αM 3/2 /π and γ = M (M − 8π)/(2π). Then we need to solve the recursive inequality By definition of I * , we have f (I * ) = 0. Since f is increasing and I 0 ≤ I * , it holds that f (I 0 ) ≤ 0. We proceed by induction. Let f (I k−1 ) ≤ 0. We suppose that f (I k ) > 0 and show that this leads to a contradiction. Inequality (21) is equivalent to Since f (I k−1 ) ≤ 0, the first factor is larger than or equal to one, and taking into account which contradicts the smallness condition τ < τ * = γ/β 2 . We infer that f (I k ) ≤ 0.
Remark 3 (Semi-discrete energy dissipation). It is possible to design semi-discrete schemes that dissipate the discrete free energy (and conserve the mass and preserve the positivity). An example, taken from [38], is the semi-implicit scheme Indeed, by the convexity of s → s log s, we obtain Furthermore, translating the computation in [35, Section 5.2.1] to the semi-discrete case, Subtracting the latter from the former expression, we conclude that Unfortunately, scheme (22) does not allow us to apply the symmetrization argument used in the proof of Theorem 2, since the drift part depends on two different time steps. The question whether (22) admits solutions that blow up in finite time remains open.

Higher-order schemes
We investigate BDF and general Runge-Kutta schemes.
3.1. BDF-2 scheme. The scheme reads as for k ≥ 2, where n 0 is given and n 1 is computed from n 0 using the implicit Euler scheme.
Lemma 4 (Existence for the BDF-2 scheme). Let α ≥ 0, n k−2 , n k−1 ∈ L 1 (R 2 ) ∩ L ∞ (R 2 ), and τ ≤ 3 2 Then there exists a weak solution n k ∈ L 1 (R 2 )∩L ∞ (R 2 )∩H 1 (R 2 ) to (23) with the following properties: • conservation of mass: The proof is exactly as for Theorem 1 since we can formulate scheme (23) as and the first term on the right-hand side plays the role of n k−1 in the implicit Euler scheme. The only difference to the proof of Theorem 1 is that we replace n + k in (6) by n k , since the truncation was only needed to show the nonnegativity of n k , which we are not able to show for the BDF-2 scheme. Let (n k ) ⊂ L 1 (R 2 ) ∩ H 1 (R 2 ) be a sequence of nonnegative weak solutions to (23). Then this sequence is finite with maximal index k max , where k max is bounded from above according to (17) (if α = 0) or (18) (if α > 0 and additionally I 0 ≤ I * and τ ≤ τ * , where I * and τ * are defined in (16)).
Proof. The proof is similar to that one of Theorem 2 but the iteration argument is different. First, let α = 0. We know from the proof of Theorem 2 that (24) where we recall that γ = M (M − 8π)/(2π). Using the same approximation of |x| 2 as in step 5 of the proof of Theorem 2, we can justify the weak formulation (see (15)) The last integral can be calculated as in the proof of Theorem 2 and we end up with We iterate this identity and insert (24): As in the proof of Theorem 2, this leads to a contradiction for large values of k.
Next, let α > 0. As the first step is computed with the implicit Euler scheme, the proof of Theorem 2 gives the estimate Moreover, since f (I 0 ) ≤ 0, we know that I 1 ≤ I 0 , and this gives f (I 1 ) ≤ 0.
For the following time steps, we obtain Let us assume, by induction, that I k−1 ≤ I k−2 and f (I k−1 ) ≤ 0 for k ≥ 2. We will prove that I k ≤ I k−1 and f (I k ) ≤ 0. Assume by contradiction that f (I k ) > 0, which is equivalent to I 1/2 k > γ/β. Then, using I k−1 − I k−2 ≤ 0, we reformulate (25) as 3 .
Since f (I k−1 ) ≤ 0, the first factor is larger than or equal to one, so and τ ≤ γ/β 2 leads to We conclude that f (I k ) ≤ 0 and therefore, by (25), I k ≤ I k−1 ≤ 0, showing the claim.
3.2. BDF-3 scheme. The finite-time blow-up can be also shown for solutions to higherorder BDF schemes, at least in the case α = 0. As an example, let us consider the BDF-3 scheme where n k−1 , n k−2 , and n k−3 are given. The existence of solutions can be shown as in Lemma 4. First, we prove that the scheme preserves the mass.
Lemma 6 (Conservation of mass). Let n 0 , n 1 , and n 2 be given and having the same mass M . Then the solution n k has the same mass, R 2 n k dx = M , for k ≥ 3. Proof. We proceed by induction. Employing the mollified version of the cut-off function (13) as a test function in (26) and passing to the limits ε → 0 and R → ∞, we arrive at since n 2 , n 1 , and n 0 have the same mass M . For the induction step, if n k−1 , n k−2 , and n k−3 for k > 4 have the same mass, the same argument as above shows that R 2 n k = M . We recall that I k = R 2 n k |x| 2 dx for k ∈ N 0 and γ = M (M − 8π)/(2π). Theorem 7 (Blow-up for the BDF-3 scheme). Assume that α = 0, I 2 −I 1 = I 1 −I 0 = −τ γ, and Let (n k ) ⊂ L 1 (R 2 ) ∩ H 1 (R 2 ) be a sequence of nonnegative weak solutions to (26). Then this sequence is finite with maximal index k max , where k max is bounded from above according to (17) (if α = 0) or (18) (if α > 0 and additionally I 0 ≤ I * and τ ≤ τ * , where I * and τ * are defined in (16)).
Proof. We claim that I k − I k−1 = −τ γ. To prove this, we proceed by induction. Let k = 3. We take an approximation of |x| 2 as a test function in (26). Then, arguing as in the previous sections, we find that Since I 2 − I 1 = I 1 − I 0 = −τ γ, it follows that For the induction step, we assume that I k−1 − I k−2 = I k−2 − I k−3 = −τ γ for k > 3. Then, as above, which shows that I k − I k−1 = −τ γ. As in the proof of Theorem 2, this leads to a contradiction for large values of k.
The previous proof can be generalized to all BDF-m methods where a i ∈ R satisfy m i=0 a i = 1. Note, however, that only the BDF-m schemes with m ≤ 6 are A(α)-stable, while they are instable for m > 6.

3.3.
Runge-Kutta schemes. The Runge-Kutta scheme reads as follows: where s ∈ N is the number of stages, b i ≥ 0 are the weights, and a ij are the Runge-Kutta coefficients. We assume that s i=1 b i = 1. The existence of solutions is only shown for two particular Runge-Kutta schemes; see below.
First, we claim that the mass is conserved in the following sense.
Lemma 8 (Conservation of mass). Let n k ∈ L 1 (R 2 ) be a solution to (27) such that m i ∈ L 1 (R 2 ) and Note, however, that we do not know whether m i ≥ 0 in R 2 .
Proof. Using the mollifier ρ ε and the cut-of function φ R (x) = φ(|x|/R), where φ is defined in (13), as a test function in the weak formulation of the equation for K i and performing the limit ε → 0, we find that According to (14), the first term on the right-hand side can be estimated as For the second term, we use formulation (8) of ∇B α , the symmetry argument, and the Lipschitz estimate (14) for ∇φ R , which leads to We deduce that for R → ∞, which concludes the proof.
We are able to show finite-time blow-up for all Runge-Kutta schemes if α = 0.
Theorem 9 (Blow-up for Runge-Kutta schemes). Let α = 0. Assume that be a sequence of nonnegative weak solutions to (27). Then this sequence is finite with maximal index k max defined in (17).
Proof. Using an approximation of |x| 2 as a test function in (27) and passing to the deregularization limit (see step 5 of the proof of Theorem 2), we find that By Lemma 8, the symmetry argument, and s i=1 b i = 1, Now, we argue as in the proof of Theorem 2 to conclude.
The case α > 0 is more delicate since m i ≥ 0 is generally not guaranteed. Indeed, it follows that (see (19) and (20)) where we recall that γ = M (M − 8π)/(2π). By the Cauchy-Schwarz inequality, , and this cannot be estimated further as m i ≥ 0 may not hold. However, for the midpoint and trapezoidal rule, we are able to give a result.
Proof. We set n := n k and n 0 := n k−1 . For given n ∈ X, we solve the linear problem By Lemma 15, this problem has a unique solution n ∈ H 1 (R 2 ), and it can be represented by writing ∇c[n] = ∇B α * n as in section 2. This defines the fixed-point operator T : S → S, where S = {n ∈ X : n X ≤ C B } and It holds T (S) ⊂ S since, using similar arguments as in the proof of Theorem 1 and the smallness condition on τ , We claim that T : S → S is a contraction. Indeed, let n, m ∈ S. Then and we have πb √ τ (C B + 1)/2 < 1. The Banach fixed-point theorem now implies that there exists a unique fixed point n ∈ X.
By the same arguments used in step 2 of the proof of Theorem 1, we infer that n k ∈ H 1 (R 2 ). Steps 4 and 5 show the conservation of mass and the finiteness of the second moment.
It remains to show that if n k−1 ∈ Y then n k has the same regularity. By Lemma 13, n k ∈ H 1 (R 2 ) implies that ∇B α * n k ∈ H 2 (R 2 ). Therefore, Elliptic regularity then gives n k ∈ H 2 (R 2 ). We bootstrap this argument to find that n k ∈ H 3 (R 2 ) → W 1,∞ (R 2 ). Taking the L 1 norm of the gradient of n = n k in (29) shows that ∇n k L 1 (R 2 ) can be estimated in terms of the H 3 norms of n k , n k−1 , and c[n k ]. We conclude that n k ∈ W 1,1 (R 2 ), finishing the proof.
Lemma 11 (Blow-up for the midpoint scheme). Let α > 0. Assume that be a sequence of nonnegative weak solutions to (28). Suppose that I 0 ≤ I * and τ ≤ 2τ * (see (16)). Then this sequence is finite with maximal index k max defined in (18).
Proof. Approximating |x| 2 as in step 5 of the proof of Theorem 2 and using the nonnegativity of n k and n k−1 , we can estimate as Setting again f (s) = β √ s − γ, it follows that Again, since f (I * ) = 0 and f is increasing, we have f (I 0 ) ≤ 0. Let f (I k−1 ) ≤ 0. Then and we can proceed as in the proof of Theorem 2.
Trapezoidal rule. The (implicit, two-stage) trapezoidal rule is defined by s = 2, a 11 = a 12 = 1 2 , b 1 = b 2 = 1 2 , which gives the scheme The existence of weak solutions can be shown exactly as in the proof of Lemma 10, therefore we leave the details to the reader.
Proposition 12 (Finite-time blow-up for the trapezoidal rule). Let α > 0. Assume that Let (n k ) be a sequence of nonnegative weak solutions to (30). Suppose that I 0 ≤ I * and τ ≤ τ * (see (16)). Then this sequence is finite with maximal index k max defined in (18).
Proof. Arguing as in the previous blow-up proofs, we obtain k−1 . Now, we proceed as in the proof of Proposition 11.

Numerical examples
The numerical experiments are performed by using the finite-element method introduced by Saito in [36] and analyzed in [37]. In contrast to [36], we choose higher-order temporal approximations. The scheme uses a first-order upwind technique for the drift term, the lumped mass method, and a decoupling procedure. We take α = 1 in (1) and consider bounded domains only. Equations (1) are supplemented with no-flux boundary conditions. In the first example, the domain is large enough to avoid effects arising from boundary conditions. The second example, on the other hand, illustrates blow-up at the boundary.
where v h ∈ X h , and the mass-lumped inner product is given by where (·, ·) 2 is the L 2 inner product.
For the discretization of the drift term, we define the discrete Green operator G h : (Recall that α = 1.) The drift term is approximated by the trilinear form b h : where Λ i is the set of vertices P j that share an edge with P i . For the definition of β ± ij , we first introduce the set S ij h of all elements K ∈ T h such that P i , P j ∈ K and the exterior normal vector ν ij to ∂D i ∩ ∂D j with respect to D i . Then where [s] ± = max{0, ±s} for s ∈ R. It is explained in [36] that the trilinear form approximates the integral Ω v∇(Gu) · ∇wdx, where Gu is the Green operator associated with −∆ + 1 on L 2 . Equation (2) is solved in a semi-implicit way. This means, for the BDF-2 scheme and for given u h , that we solve the linear problem (31) 1 τ for all w h ∈ X h . Here, n k h is an approximation of n(·, τ k). This defines the solution operator N (u h ) = n h and the scheme is completed by chosing u h . Saito has taken u h = n k−1 h , giving the usual semi-implicit scheme of first order. For higher-order schemes, we need to iterate. For this, we introduce the iteration u where (x 0 , y 0 ) ∈ (0, 1) 2 , M > 0, and θ > 0. Clearly, the mass of W x 0 ,y 0 equals M . For the first example, we choose θ = 1/500, M = 6π, and n 0 = W 0.33,0.33 + W 0.33,0.66 + W 0.66,0.33 + W 0.66,0.66 .
The initial mass is 24π > 8π and thus, we expect the solutions to blow up in finite time. Figure 1 shows the cell density n(x, t) at various time instances computed from the BDF-2 scheme. As expected, the solution blows up in finite time in the center of the domain. Note that the numerical solution is always nonnegative and conserve the total mass. A nonsymmetric situation is given by the initial data taking θ = 1/500, M = 6π, and the same numerical parameters as above. The total mass of n 0 is 11π, so we expect again blow up of the solutions. This is illustrated in Figure 2.
The solution aggregates, moves to the boundary and blows up. Again, the discrete solution stays nonnegative and conserves the mass.

Convergence rate.
To calculate the temporal convergence rates and to show that the schemes are indeed of second order, we compute a reference solution n ref with the very small time step τ = 10 −6 and compare it in various L p norms with the solutions n τ using larger time step sizes τ . We choose the same initial datum as in the first example with M = 24π. Figure 3 shows the L p error where the end time T = 0.01 is chosen such that the density already started to aggregate but blow up still did not happen. As expected, the L p errors are approximately of second order.

4.4.
Numerical blow-up. We demonstrate that the bound for the blow-up time T * = τ k max derived in the time-discrete situation can serve as a bound for the numerical blow up. It is well known that the computation of the numerical blow-up time is rather delicate. For instance, Chertock et al. [14] use the L ∞ norm of the density as a measure of the numerical blow-up time, since n L ∞ (Ω) is proportional to h 2 (recall that h is the spatial grid size). Numerical blow-up may be reached when the numerical solution becomes negative [13,41] or when the seond moment becomes negative [22]. However, since our scheme conserves the mass and the grid is finite, the numerical solution cannot blow up in the L ∞ norm. Instead, the solution converges to a state where the mass concentrates at certain points and no further growth is possible. Moreover, as the scheme conserves the nonnegativity numerically, the second moment cannot become negative, but it will be small near blowup. A lower bound for the blow-up time was derived in, for instance, [ For the numerical test, we choose the initial datum n 0 = W 1,1 on the domain Ω = (0, 2) 2 with parameters θ = 1/500, M = 30π, and τ = 10 −5 . The grid sizes are h = 0.02, 0.04, 0.08. Figure 4 illustrates the evolution of n k L ∞ (Ω) and I k . The vertical line marks the bound k max from (18). We observe that the L ∞ norm and the second moment reach a limit close to k max . The initial density and the density at k max are displayed in Figure 5 for comparison.

Appendix A. Some auxiliary results
We recall that the function   defined for x ∈ R n , is called the Newton potential if α = 0 and the Bessel potential if α = 0. We need the following properties of the Bessel potential.
Lemma 15 (Elliptic problem). Let τ > 0, f ∈ L 2 (R 2 ) 2 , and g ∈ L 2 (R 2 ). Then there exists a unique weak solution n ∈ H 1 (R 2 ) to (34) − ∆n + τ −1 (n − g) = − div f in R 2 , and this solution can be represented as Equation (34) correspond to the implicit Euler discretization of a parabolic problem with n being the solution at the actual time step and g being the solution at the previous time step. Although the result is standard, we give proof for the sake of completeness.
Taking the test function n k − n in the difference of the weak formulations for n k , n corresponding to f k , f , respectively, it follows that .
Since (f k ) is a Cauchy sequence, (n k ) is a Cauchy sequence in H 1 (R 2 ) and hence there exists a n ∈ H 1 (R 2 ), such that n k → n strongly in H 1 (R 2 ) as k → ∞. Therefore, we can perform the limit k → ∞ in the weak formulation of (36) leading to (34). It remains to show (35). Let n = (1/τ )B 1/τ * g − ∇B 1/τ * f . Then, by Lemma 14 and (37), The right-hand side can be made arbitrarily small by choosing k sufficiently large. This shows that n = n in R 2 .