Global classical solutions for mass-conserving, (super)-quadratic reaction-diffusion systems in three and higher space dimensions

This paper considers quadratic and super-quadratic reaction-diffusion systems for reversible chemistry, for which all species satisfy uniform-in-time $L^1$ a-priori estimates, for instance, as a consequence of suitable mass conservation laws. A new result on the global existence of classical solutions is proved in three and higher space dimensions by combining regularity and interpolation arguments in Bochner spaces, a bootstrap scheme and a weak comparison argument. Moreover, provided that the considered system allows for entropy entropy-dissipation estimates proving exponential convergence to equilibrium, we are also able to prove that solutions are bounded uniformly-in-time.


Introduction
In this article, we first consider the following quadratic reaction-diffusion system: on Ω T , on Ω T , n(x)·∇ x u i (t, x) = 0, on ∂Ω T , where u i := u i (t, x) ≥ 0 for i = 1, 2, 3, 4 denotes non-negative concentrations at time t and position x of four species A i subject to non-negative initial concentrations u i0 (x) ≥ 0. Moreover, d i > 0 are the corresponding positive and constant diffusion coefficients. We suppose x ∈ Ω, where Ω is a bounded domain of IR N (N ≥ 1) with sufficiently smooth boundary ∂Ω ∈ C 2+α , α > 0. Finally, n(x) is the outer normal unit vector at point x of ∂Ω and we denote Ω T = [0, T ] × Ω and ∂Ω T = [0, T ] × ∂Ω for any T > 0.
The above system (1) constitutes a mass-action law model of the evolution of a mixture of four diffusive species A i , i = 1, 2, 3, 4 undergoing the single reversible reaction until a unique positive detail balance equilibrium is reached in the large time behaviour, see Section 2. For a derivation of system (1) or related mass-action law reaction-diffusion systems from kinetic models or fast reaction limits, we refer to [2,3,4,8].
Note that for the sake of readability, we have set the forward and backward reaction rates in (1) equal to one (the general case can be treated without any additional difficulty). Also, without loss of generality, we shall also assume that Ω is normalised (i.e. |Ω| = 1), which can always be achieved by rescaling the spatial variable.
The quadratic reaction-diffusion system (1) and in particular the question of global-in-time solutions has lately received a lot of attention, see e.g. [5,12,13].
In [12], for instance, a duality argument in terms of entropy density variables was used to prove the existence of global, non-negative L 2 -weak solutions in any space dimension, see also Section 2.
While the existence of global classical solutions in 1D follows already from Amann, see e.g. [1], the existence of global classical solutions in 2D was recently shown by Goudon and Vasseur in [13] by using De Giorgi's method. For higher space dimensions the existence of classical solutions constitutes in general an open problem, for which the Hausdorff dimension of possible singularities was characterised in [13].
The (technical) criticality of quadratic nonlinearities was underlined by Caputo and Vasseur in [6], where smooth solutions were shown to exist in any space dimension for systems with nonlinearities of sub-quadratic power law type, see also e.g. [1].
A further related result by Hollis and Morgan [14] showed that if blow-up (here that is a concentration phenomenon since the total mass is conserved) occurs in system (1) in one concentration u i (t, x) at some time t and position x, then at least one more concentration u j =i has also to blow-up (i.e. concentrate) at the same time and position. The proof of this result is based on a duality argument.
Recently in [5], an improved duality method allowed to show global classical solutions of system (1) in 2D in a significantly shorter and less technical way than De Giorgi's method.
In this work, we will prove in higher space dimensions N ≥ 3 and under a dimension-dependent closeness condition of the diffusion coefficients d i that the global, L 2 -weak solutions to (1) (see [12]) are in fact global classical solutions. More precisely, we shall denote by the maximal distance between the diffusion rates appearing in (1) and prove the following: Theorem 1.1 (Global classical solutions and uniform-in-time bounds). Consider N ≥ 3 and assume that δ, a and b as given in (2) satisfy where C a+b 2 ,2Γ is defined in Proposition 2.1. Assume initial data u i0 ∈ L 3N 4 (Ω) (which also implies with 3N 4 > 2Γ, assumption (3), estimate (17) of Proposition 2.1 and the theory of global, weak solutions in [12], the existence of global, weak solutions u i ∈ L 2Γ (Ω T ) to (1) for any T > 0).
Then, any global, weak solutions u i ∈ L 2Γ (Ω T ) for any T > 0 is a global, classical solutions to (1). Moreover, these solutions are bounded uniformly-in-time in the sense that for all τ > 0, there exists a positive constant C = C(τ ) such that Remark 1.2. The actual novelty of Theorem 1.1 concerns the range Γ ∈ ( N +2 2 − 2N N +2 , N +2 2 ) because when Γ ≥ N +2 2 , then the quadratic right hand side terms of (1) is bounded in L N +2 2 (Ω T ) and standard parabolic regularity estimates (which reflect the convolution with the heat kernel being in L N +2 N −ε , for all ε > 0) implies u i ∈ L ∞ (Ω T ) and thus classical solutions, see e.g. Proposition 2.3 from [5], where the more stringent δ-smallness condition (21) was assumed, which yields directly u i ∈ L N +2 (Ω T ).
In this paper, we are able to show classical solutions under the weaker δ-smallness condition (3) by combining regularity and bootstrap arguments in Bochner spaces with interpolation of Bochner spaces (see Lemma 4.1) and a further bootstrap based on a weak comparison argument (Lemma 3.6).
For example, for N = 3, our approach lowers the initially required regularity u i ∈ L 2Γ (Ω T ) from the seemingly natural assumption Γ = 5 2 (as assumed in Proposition 2.3 from [5]) down to Γ > 13 10 . Since global L 2 -weak solutions to (1) exist in all space dimensions, Theorem 1.1 leaves the gap 1 < Γ ≤ 13 10 in the case N = 3, for which it remains an open problem if global weak solutions u i ∈ L 2Γ (Ω T ) are also global classical solutions. This corresponds to systems (1), for which Remark 1.3. We remark that the assumption on the initial data u i0 ∈ L 3N 4 (Ω) is not optimal, but chosen for sake of the readability of the proof of Theorem 1.1. In fact, the proof of Theorem 1.1 would work equally for u i0 ∈ L µ0 (Ω) with µ 0 > max ν, (ν−1)N 2 . Remark 1.4. We further remark that the uniform-in-time L ∞ -bound (4) follows from an interpolation argument between the exponential convergence to equilibrium of system (1) (see Proposition 2.2 below) and the fact that the proof of Theorem 1.1 only involves regularity constants, which grow not faster than polynomially in time.
Remark 1.5. There is an alternative argument to derive the initially required regularity u i ∈ L 2Γ (Ω T ) provided a suitable δ-smallness condition of the diffusion coefficients d i of system (1). System (1) implies, for instance that Then, by a duality estimate (see e.g. the proof of [18,Lemma 3.4]), it follows directly for any p ∈ (1, +∞) and some C = C(p, T ) that Thus, with u 3 L p (ΩT ) ≤ u 1 + u 3 L p (ΩT ) , the required estimate for u 3 L p (ΩT ) holds provided that We emphasise that our proof of Theorem 1.1 only relies on the closeness assumption (3), uniform-intime L 1 a-priori estimates (which follows typically from mass conservation laws) and the quadratic order of the non-linearity of the reaction terms on the r.h.s. of (1). Thus our results can be similarly applied to related mass-conserving reaction-diffusion systems with quadratic non-linearities. Moreover, one can also easily generalise the proof to systems with super-quadratic reaction terms provided the more stringent closeness condition (7), see the following Corollary. Corollary 1.6. Let N ≥ 3 and ν ≥ 2 and consider the following non-linear reaction-diffusion systems: on Ω T , where u = (u 1 , . . . , u k ) ≥ 0 denotes k ∈ IN + non-negative concentrations and for all i ∈ {1, . . . , k} the nonlinearities f i (u) are positivity-preserving, polynomials of order ν, i.e. f i (u) ≤ C k i=1 u ν i for a constant C > 0 and satisfy conservation laws of the following form: For all i ∈ {1, . . . , k}, there exists a constant, positive vector 0 ≤ (a i j ) j∈{1,...,k} ∈ IR k with a i i > 0 such that k j=1 a i j f j (u) ≤ 0 and thus where C a+b 2 ,2Γ is defined in Proposition 2.1. Assume initial data u i0 ∈ L (2ν−1)N 4 (Ω) (which also implies with (2ν−1)N 4 > 2Γ, assumption (7), estimate (17) of Proposition 2.1 and a theory of global, weak solutions analog to [12], the existence of global, weak solutions u i ∈ L 2Γ (Ω T ) to (1) for any T > 0).
Then, any global, weak solutions u i ∈ L 2Γ (Ω T ) for any T > 0 is a global, classical solutions to (1).
Remark 1.7. We remark that in contrast to Theorem 1.1, Corollary 1.6 does not state a uniform-in-time bound analog to (4). In fact, general systems (5) will not feature an entropy functional like system (1). Without the entropy method proving exponential convergence to equilibrium, our method only shows that global classical solutions to (5) grow (at most) polynomially-in-time. More precisely, it is due to the uniform-in-time L 1 -bounds (6) that we are able to show global solutions with polynomially growing bounds for quadratic and super-quadratic systems.
For the existence of classical solutions to systems with sub-quadratic nonlinearities (i.e. ν < 2), we refer to Caputo and Vasseur [6] and the references therein.
Remark 1.8. Again, like in Remark 1.2, the novel range is Γ ∈ ( (ν−1)N 2 +4 2(N +2) , N +2 2 ) because when Γ ≥ N +2 2 , then L ∞ -bounds and the existence of classical solutions follows from standard arguments. We remark that depending on the polynomial order ν, Corollary 1.6 only constitutes an improvement under to following condition on the dimension N : which is false for ν large. Thus, our approach does not lead to improvements in the case of strongly super-quadratic nonlinearities. However, for N = 3, we are able to lower the threshold, which is sufficient to show the existence of classical solutions, for systems with third order and slightly higher nonlinearities.
Notation: Besides standard notations, we shall denote by L p,q (Ω T ) the Bochner space with space-time norm: Idea of the proofs: For the convenience of the reader we present a short outline of the proofs of Theorem 1.1 and Corollary 1.6: First, we use an improved duality estimate as proven in [5] in order to derive the starting a-priori estimate u i ∈ L 2Γ (Ω T ) (or u i ∈ L νΓ (Ω T ) respectively). Then, we apply Sobolev's embedding theorem, interpolation with the uniform-in-time L 1 -bound and other classical inequalities in order to derive L p -estimates for some p to be appropriately chosen. Next, we use these estimates as a starting point for a bootstrap scheme, which allows to improve the regularity of the solutions to u i ∈ L pn,∞ , i.e. u i is L ∞ in time with values L pn in space. Then, after bootstrapping p n sufficiently large, the regularity u i ∈ L pn,∞ allows to use the weak comparison Lemma 3.6, which leads to global-in-time L ∞ and, thus, classical solution to system (1). Finally, uniform-in-time bounds for solutions of systems (1) follow from an interpolation argument with the exponential convergence to equilibrium.
Outline: In Section 2, we shall first state basic properties of the systems (1) and recall previous existence results and methods. Then, in Section 3, we shall first prove several Lemmas required for the proof of Theorem 1.1 and finally prove Corollary 1.6. In the Appendix 4, we provide a proof of Lemma 4.1 for the convenience of the reader.

Preliminaries
Non-negativity, mass conservation laws and equilibrium The chemical reaction terms on the r.h.s. of (1) satisfy the quasi-positivity property and thus ensure that solutions to (1) propagate the non-negativity of the initial data u i0 ≥ 0, i.e. for all T > 0, we have on Ω T , for all i = 1, 2, 3, 4.
Moreover, system (1) subject to homogeneous Neumann boundary conditions satisfies the following mass conservation laws: where j = {1, 2} and k = {3, 4}. Notice that only three of these four conservation laws are linearly independent.
The non-negativity of the solutions together with the mass conservation laws (8) provide natural uniform-in-time a-priori L 1 -estimates for the concentrations, i.e.
where M denotes the total initial mass M = M 13 + M 24 = M 14 + M 23 .
System (1) features a unique positive detail balance equilibrium. Due to the homogeneous Neumann boundary conditions this equilibrium (u i,∞ ) i=1,..,4 consists of the unique positive constants balancing the reversible reaction u 1,∞ u 2,∞ = u 3,∞ a u,∞ and satisfying the conservation laws u j,∞ + u k,∞ = M jk for (j, k) ∈ ({1, 2}, {3, 4}), that is: Duality estimates, entropy variables and global weak solutions The system (1) can also be rewritten in terms of the entropy density variables z i := u i log(u i ) − u i (compare with the entropy functional (18) below). By introducing the sum z := 4 i=1 z i , it holds (with a and b defined as in (2)) that Then, by a duality argument (see e.g. [12,14,19] and the references therein), the parabolic problem (11) satisfies for all T > 0 and for all space dimensions N ≥ 1 the following a-priori estimate where C is a constant independent of T , see [5,12]. Thus, given (u i0 ) i=1,..,4 ∈ L 2 (log L) 2 (Ω), we have (u i ) i=1,..,4 ∈ L 2 (log L) 2 (Ω T ) and the quadratic nonlinearities on the right hand side of (1) are uniformly integrable, which allows to prove the existence of global L 2 -weak solutions in all space dimensions N ≥ 1, see [12]. Moreover, in 2D and in higher space dimension under the assumption of sufficiently strong δ-smallness condition, the following duality Lemma allows to show global classical solutions: Let Ω be a bounded domain of IR N with smooth (e.g. C 2+α , α > 0) boundary ∂Ω and T > 0. We consider a coefficient function Then, any function u satisfying and (with p ′ < 2 denoting the Hölder conjugate of p) satisfies for all T > 0 Here, the constant C T depends polynomially on T and the constant C m,q > 0 in (14) is defined for m > 0, Specifically for system (1), we can rewrite the system in terms of z = u 1 + u 2 + u 3 + u 4 (it is not even necessary to consider entropy density variables) where a and b are defined as in (2) and obtain for all T > 0 ∀i = 1, 2, 3, 4 : Entropy functional, entropy method and exponential convergence to equilibrium The detail balance structure of system (1) ensures also the following non-negative entropy functional E((u i ) i=1,..,4 ) and the associated entropy dissipation functional D( It is easy to verify that the following entropy dissipation law holds (for sufficiently regular solutions (u i ) i=1,..,4 of (1)) for all t ≥ 0 In [11], exponential convergence in L 1 towards the unique constant equilibrium (10) (with explicitly computable rates) was shown for global L 2 -weak solutions in all space dimensions N . The proof of this statement was based on the so called entropy method, where a quantitative entropy entropy-dissipation estimate of the form with an explicitly computable constant C was established, which uses only natural a-priori bounds of the system and thus significantly improved the results of [10].
Let Ω be a bounded domain with sufficiently smooth boundary such that Poincaré's and Sobolev's inequalities and thus the Logarithmic Sobolev inequality hold. (8)).
Then, any global solution (u i ) i=1,..,4 of (1), which satisfies the entropy dissipation law (2) (this is true for any weak or classical solutions as shown to exist in [5,12]) decays exponentially towards the positive equilibrium state (u i,∞ ) i=1,..,4 > 0 defined by (10): for all t ≥ 0 and for constants C 1 and C 2 , which can be explicitly computed.
The following Proposition 2.3 from [5] shows the consequences of the Propositions 2.1 and 2.2 specific to the system (1).

Existence of Global Classical Solutions
In this Section, we first prove the existence of global classical solutions of the system (1) under the δ-smallness assumption (3), where δ (as defined in (2)) measures the maximal distance of the diffusion coefficients of systems (1) on a domain Ω ⊂ IR N for a given space dimension N ≥ 3 (recall that global classical solutions of system (1) are known for N = 1, 2). As a consequence of the δ-smallness assumption (3), estimate (17) of Proposition 2.1 provides an u i ∈ L 2Γ (Ω T )-estimate for the concentration u i of (1).
In particular, we are interested in the range since for Γ ≥ N + 2 and the considered quadratic nonlinearities f ≤ |u 1 u 2 − u 3 u 4 | ∈ L Γ (Ω T ), a standard parabolic regularity bootstrap argument for the heat equation with Neumann boundary conditions, i.e.
implies for a right-hand-side f ∈ L N +2 2 (Ω T ) that the solution satisfies the following L ∞ estimates where C T grows at most polynomially w.r.t. T , (see [5] for the polynomial dependence on T ). In particular, the range Γ ≥ N +2 2 was considered in Proposition 2.3, where δ was assumed to satisfy the more stringent δ-smallness assumption (21), which yields f ∈ L N +2 2 (Ω T ).
Thus, we consider here the δ-smallness assumption (3) and estimate (17) of Proposition 2.1 yields a-priori estimate of the form where C T is a constant depending polynomially on the time T and N +2 2 − 2N N +2 ≥ 13 10 for N ≥ 3.
As a first step in proving Theorem 1.1, we are faced with the fact that duality estimates like (25) address spaces L 2Γ,2Γ (Ω T ), which feature equal time-and space-integrability. However, as we shall see in the following, bootstrapping the problem (24) naturally leads to Bochner spaces with higher time-than space-integrability.
In this Section, we present a bootstrap, which allows to conclude from u i ∈ L 2Γ,2Γ to u i ∈ L 3N 4 ,∞ . Furthermore, we shall take care to ensure that all involved constants grow at most polynomially in T . (for instance as a consequence of a suitable mass conservation law) and consider initial data u 0 ∈ L p (Ω) for a p > 1 to be chosen. Then, (by using Sobolev's embedding theorem, the uniform-in-time bound (26) as well as Hölder's and Young's inequality) u satisfies the following L p estimate for all T > 0: where q = q(p, N ) := p N N −2 > p for N ≥ 3 and the constants C ′ andC only depend on p, M , d, N and C s , which is the constant from the Sobolev's embedding theorem. Moreover, γ, γ ′ , µ and µ ′ are Hölder conjugates to be chosen.
Furthermore, by using once more the uniform-in-time bound (26) and provided that 1 ≤ µ ′ (p − 1) ≤ q holds, u also satisfies Proof. By testing (24) with p |u| p−1 sgn(u) (more precisely by testing with a smoothed version of the modulus |u| and its derivative sgn(u) and letting then the smoothing tend to zero), we obtain by integrationby-parts and with the constant C 0 (p, d) :

Sobolev's embedding theorem
After adding on both sides C 0 Ω |u| p dx, we apply Sobolev's embedding theorem for N ≥ 3, i.e.
with Sobolev constant C s and and q(p) > p for all N ≥ 3.
Therefore, we get:

Interpolation and the mass conservation property
Here we remark that integration of (29) by means of a Gronwall type argument would lead to global, yet exponentially growing estimates of u p p . However, such exponential growing estimates can be avoided (and should be avoided in order to retain the possibility of interpolating a-priori estimates with exponential convergence to equilibrium) by using the mass conservation laws (8), which provides a uniform-in-time L 1 (Ω) bound of the form (26). More precisely, we interpolate the L p term u p p on the right hand side of (29) between L 1 and L q , i.e. 1 p = θ 1 + 1−θ q such that θ = q−p q−1 1 p ∈ (0, 1) since 1 < p < q. Thus, due to the uniform-in-time bound u(t, ·) 1 ≤ M for all t ≥ 0, we obtain with an exponent Now, since p > s, the last term on the right hand side of (30) can be controlled by the term u p q on the left hand side of (30) as in the following step.

Young's and Hölders's inequality, time-integration and further Hölder inequality
Next, we apply Young's inequality to the last term on the right hand side of (30) to derive for a δ > 0 to be chosen sufficiently small In the following, we estimate the first term on the right hand side with Hölders's inequality. In fact, this shall first be done for a general Hölder exponent 1 where µ can be chosen in order to match the necessary assumptions.
Remark 3.2. We remark at this point that we will mostly chose µ such that µ ′ (p − 1) = q, i.e. we shall choose µ according to the maximal space integrability, which can be controlled by the left hand side term u p q . Next, we integrate over (0, T ) and apply Hölder again: which proves the estimate (27).

Further interpolation with the uniform-in-time L 1 -bound
We first remind that In order to further exploit the mass conservation property in cases such that 1 ≤ µ ′ (p − 1) ≤ q, we will interpolate L µ ′ (p−1) ֒→ L q ∩ L 1 with a parameter θ ∈ [0, 1] that satisfies .
The estimate (27) still allows to choose µ and γ in our goal of deriving a-priori estimates for u. At this point, there are two basic options: First, we could aim to control the term u p−1 µ ′ (p−1),γ ′ (p−1) entirely by the duality estimate (25), which we have anyway already used to ensure the integrability f ∈ L Γ,Γ and u i ∈ L 2Γ,2Γ . Thus, in this first case, we choose µ and γ such that If we consider p = 2, then (32) requires to assume the δ-closeness condition in Proposition 2.1 small enough such that the duality estimate ensures u i ∈ L 3,3 and f ∈ L 3 2 , 3 2 . In the case N = 3, under such a δ-closeness assumption, (32) then implies directly u i ∈ L 2,∞ and we can apply Lemma 3.6 in order to establish global classical solutions in 3D, see also [16]. In order to show classical solutions in higher space dimensions N > 3, we can proceed similar by considering (32) for suitable exponents p > 2. However, the idea to directly control the term u p−1 µ ′ (p−1),γ ′ (p−1) via a duality estimate (25) does not lead to such general results as stated in Theorem 1.1. This which will become clear from the following considerations.
From now onwards, we are interested in the second case, where the term u p−1 µ ′ (p−1),γ ′ (p−1) cannot be controlled by the duality estimate (25), but has to be controlled by the term u p q,p on the left hand side of (27). In order to constitute an improvement of the argument described in the previous paragraph, we are interested in cases when holds.
Remark 3.3. In fact, we will see in the following that our method naturally considers µ ≤ γ, i.e. that the time-integrability is larger (or at least equal) to the space-integrability. This leads to two possible subcases where the first subcase would in principle still allow to partially use the duality estimate u 2Γ,2Γ ≤ C T in a bootstrap argument. However, by anticipating some details from the below calculations, in particular by considering the starting value of a bootstrap exponent p 0 = Γ(N +2)−2 N −2 and γ ′ (p 0 − 1) = p 0 , we see that the below arguments are able to successfully treat some cases (in dimensions N = 3, 4), where 2Γ < γ ′ (p 0 − 1) < µ ′ (p 0 − 1).
This means that the most general results as stated in Theorem 1.1 are derived from entirely controlling the term u p−1 µ ′ (p−1),γ ′ (p−1) by the term u p q,p on the left hand side of (27).
The following two Lemmata 3.4 and 3.5 develop further the estimates of Lemma 3.1 for two different choices of µ and γ: Lemma 3.4 takes into account that the duality estimates (25) only allow to control f in spaces with the same regularity in space and in time, i.e. f ∈ L Γ,Γ (Ω T ). Therefore, Lemma 3.4 sets µ = γ: N +2 and u 0 ∈ L p (Ω) for some p > 1.
Then, we have with q = pN with a constant C, which grows at most polynomially in time T .
Proof. By recalling (28) from Lemma 3.1 and setting γ = µ, we have especially Then, the first term on the right hand side can be controlled by the term u p q,p on the left hand side of (28) by setting if we verify that the corresponding Hölder exponent satisfies µ ′ > 1 (such that 1 µ ′ < 1) and that the interpolation L µ ′ (p−1) ֒→ L q ∩ L 1 used in Lemma 3.1 is admissible. This is done by recalling that q = pN N −2 and straightforward calculations show the corresponding values of µ ′ , µ and θ to be for all p > 1. Therefore, the estimate (34) becomes: Next, since µ ′ , µ > 1, we estimate u p q,p 1 µ ′ with Young's inequality and obtain for a sufficiently small enough ǫ > 0 and u 0 ∈ L p ⇒ u ∈ L p,∞ ∩ L q,p , where q = pN N −2 and the involved constants grow at most polynomially in time.
The following Lemma 3.5 will allow to perform a bootstrap argument build on estimate (33) of Lemma 3.4 in order to improve the integrability of f .
We point out that estimate (33) leads naturally to a-priori estimates for u which exhibit larger timethan space-integrability (see e.g. Lemma 4.1 for the interpolation of the space L p,∞ and L q,p ). Therefore, the following Lemma 3.5 provides a suitable variant of Lemma 3.4 based on the assumption that f features higher time-than space-integrability. More precisely, we shall assume f ∈ L r,p where r = r(p, N ) is a specific space-integrability exponent with r(p, N ) < p.
Lemma 3.5 (µ < γ estimate). Consider N ≥ 3. Assume for some p > 1 that u 0 ∈ L p and f ∈ L r,p with where as above q = q(p, N ) = pN N −2 > p. Moreover, by using interpolation in Bochner spaces (see Lemma 4.1), we have for any interpolation exponent θ ∈ (0, 1): Proof. We return to the proof of Lemma 3.1 and consider µ < γ in order to derive a-priori estimates on u based on f featuring higher time-rather than space-integrability. In particular, we can set µ ′ (p − 1) = q in (31), having in mind that the L q norm is the highest possible space, which can be controlled by the the left hand side of (31), see also Remark 3.2. Thus, we get: where the index denotes the minimal space-integrability requirement on f such that the corresponding term u p−1 q can still be controlled by the left hand side term u p q . Next, we integrate (37) over (0, T ) for any T > 0 and apply on the right hand side Hölder's inequality with 1 = 1 γ ′ + 1 γ to obtain Thus, we set γ ′ (p − 1) = p and γ = p in (38) and apply once more Young's inequality to the right hand side term u p−1 q,p such that for a sufficiently small ε > 0 u(T ) p p +Ĉ u p q,p ≤ u 0 whereĈ =C − ε > 0.
Altogether, if f r,p is bounded, i.e. f ∈ L r,p and u 0 ∈ L p , we have f ∈ L r,p and u 0 ∈ L p ⇒ u ∈ L p,∞ ∩ L q,p .
We remark that in f ∈ L r,p , since we have set γ(p − 1) = p and γ ′ = p in (38), the time-integrability index p represents the minimal required time-regularity such that estimate (39) holds.
In the last step of the proof of Lemma 3.5, we use the Bochner space interpolation Lemma (4.1) (see Appendix) and get L p,∞ ∩ L q,p ֒→ L σ,τ , where σ and τ are defined depending on an interpolation parameter θ ∈ (0, 1) by Finally, the following Lemma 3.6 provides a bootstrap argument based on ideas of [15,16,21], which combines semigroup estimates with weak comparison arguments and allows to bootstrap polynomially in time growing L ∞ -bounds. Lemma 3.6 (Weak comparison bootstrap argument). Consider ν ≥ 2 and let u i (t, ·) be weak solutions to the systems (1) or (5) for i = 1, . . . , k. In particular we have ν = 2 and k = 4 for solutions u i (t, ·) of system (1). Assume that for an exponent and for all i = 1, . . . , k holds where C(t) > 0 is a continuous function in time, which grows polynomially. Then, the solution u i (·, t) satisfies where C(t) > 0 is again a (different) continuous function in time on (0, ∞), which grows polynomially.
Proof. We remind that each ( In particular, the concentrations (u i ) i=1...k are non-negative and the considered nonlinearities |f i | ≤ C k j=1 u ν j are super-quadratic in the u i of order ν ≥ 2. This motivates us, by neglecting the negative terms of f , to consider: where u ν = k j=1 u ν j for simplicity and d i is the diffusion rate of u i . Now, we define g i = λu i + Cu ν for λ > 0 and all i = 1, . . . , k and consider the linear comparison problem: on Ω T , where g i (t, ·) is considered a given function satisfying g(t, ·) L µ 0 and C 1 (t) has the same properties as C(t).
Thus, by the comparison principle of weak solutions of parabolic equations with values in the reflexive Banach space L µ 0 ν (Ω) for µ0 ν > 1 (see e.g. [7, Theorem 11.9]), we have that 0 ≤ u ≤ u. Next, we refer e.g. to Rothe [20] for the following semi-group estimates of the Laplace operator subject to Neumann boundary conditions with C(q, r) > 0, which implies for L i = d i ∆ − λ subject to Neumann boundary condition (e tLi u)(t, ·) r ≤ C(q, r)e −λt max{1, ( By taking the r-norm of (43) and applying the above semi-group estimate with q = µ0 ν , we obtain and at this point, in order to have the above right hand side integrable near the singularity s = t, we need while for µ0 ν = N 2 , we can chose any r < ∞ and for µ0 ν > N 2 , we can set r = ∞. In all these cases, we can estimate eq. (44) as and the r-norm is bounded by a polynomial of order n + 1, if C 1 is a polynomial of order n. Recalling that 0 ≤ u ≤ u, we conclude that where C 2 (t) grows polynomially in time.
In order to bootstrap the estimates (43)-(45), we require µ 1 := r > µ 0 > ν, which is possible provided that Thus, we find that assumption (40) allows to start a bootstrap argument, which leads after finitely many steps to where µn ν ≥ N 2 and C n+1 (t) grows polynomially in time. Thus, after one or two (if µn ν = N 2 ) further bootstrap steps, we reach where for all τ > 0 the bound C(t) grows polynomially in time.
We can now present the proof of our main Theorem.
Proof of Theorem 1.1 (Global existence of classical solutions to system (1)) The first part of the proof establishes a bootstrap argument based on the Lemmata 3.4 and 3.5 in order to prove the following gain of regularity: Provided that Γ > N +2 2 − 2N N +2 and initial data u i0 ∈ L 3N 4 (Ω), then Once we know that u ∈ L 3N 4 ,∞ , then Lemma 3.6 implies u ∈ L ∞,∞ with constants, which grow at most polynomially in time. In a final step, we shall then interpolate polynomially growing a-priori estimates with the exponentially convergence of solutions towards equilibrium in order to prove uniform-in-time boundedness of the solutions.

Initial
Step: A-priori estimate via duality Proposition 2.1 In order to start our argument, we recall assumption (3) and estimate (17) of Proposition 2.1, which ensures that u ∈ L 2Γ,2Γ for Γ > N +2 2 − 2N N +2 > 1 for N ≥ 3. We remark again that for Γ ≥ N +2 2 , the statement of the Theorem 1.1 follows from well known parabolic regularity results, which directly ensure u ∈ L ∞ (Ω T ) for all T > 0, see e.g. [5], which also proves the polynomial growth of the constants in T .
Thus, Gagliadro-Nirenberg interpolation between L 1 and H k for k > N 2 yields a uniform-in-time L ∞ bound, i.e. there exists a θ ∈ (0, 1) such that for all T > 0 where we have used the exponential convergence to equilibrium stated in Proposition 2.2 and the polynomial-in-time dependence of the constant C T to obtain the uniform-in-time L ∞ -bound (4). This concludes the proof of Theorem 1.1. In fact, we expect in general that σ 0 < 2Γ < τ 0 .

Proof Corollary 1.6
Proof. First, we remark that thanks due the conservation laws assumed in Corollary 1.6, we observe that for every i = 1, . . . , k the sum z i := k j=1 a i j u j where a i i > 0 satisfies a problem of the form Thus, the δ-smallness condition (7) and Proposition 2.1 ensure for all T > 0 the formal a-priori estimate u i ∈ L νΓ (Ω T ) and thus for all i = 1, . . . , k, f i (u) ∈ L Γ and it is easy verified that Γ > 1. Thus, the existence of weak global solutions to system (5) can be shown as the limit of global solutions to suitable approximating systems: By using, for instance, approximating systems, which truncate the nonlinearities f i (u) following the lines of e.g. Michel Pierre [17], we can pass to the limit in system (5) and in particular in the nonlinear terms on the right hand side f i (u) ∈ L Γ for Γ > 1. This will yield global weak solutions of (5) satisfying f i (u) ∈ L Γ .
In the following, we shall prove that these global weak solutions are indeed classical solutions. This is done analog to the proof of Theorem 1.1: We start by proving that u ∈ L νΓ,νΓ with initial data u i0 ∈ L (2ν−1)N 4 (Ω) can be bootstrapped to u ∈ L (2ν−1)N 4 ,∞ provided that Γ > (ν − 1)N 2 + 4 2(N + 2) .