Global regularity and convergence to equilibrium of reaction-diffusion systems with nonlinear diffusion

We study the boundedness and convergence to equilibrium of weak solutions to reaction-diffusion systems with nonlinear diffusion. The nonlinear diffusion is of porous medium type and the nonlinear reaction terms are assumed to grow polynomially and to dissipate (or conserve) the total mass. By utilising duality estimates, the dissipation of the total mass and the smoothing effect of the porous medium equation, we prove that if the exponents of the nonlinear diffusion terms are high enough, then weak solutions are bounded, locally H\"older continuous and their $L^{\infty}(\Omega)$-norm grows in time at most polynomially. In order to show convergence to equilibrium, we consider a specific class of nonlinear reaction-diffusion models, which describe a single reversible reaction with arbitrarily many chemical substances. By exploiting a generalised Logarithmic Sobolev Inequality, an indirect diffusion effect and the polynomial in time growth of the $L^{\infty}(\Omega)$-norm, we show an entropy entropy-production inequality which implies exponential convergence to equilibrium in $L^p(\Omega)$-norm, for any $1\leq p<\infty$, with explicit rates and constants.


Introduction and Main results
In this article, we study the boundedness and convergence to equilibrium of weak solutions to reactiondiffusion systems with nonlinear diffusion x ∈ Ω, t > 0, i = 1, . . . , S, d i ∇(u mi i ) · − → n = 0, x ∈ ∂Ω, t > 0, i = 1, . . . , S, u i (x, 0) = u i,0 (x), x ∈ Ω, i = 1, . . . , S, with the unknown functions u = (u 1 , . . . , u S ) and u i : Ω × R + → R, the positive diffusion coefficients d i > 0, the porous medium exponents m i > 1 and where Ω ⊂ R d denotes a bounded domain with sufficiently smooth boundary ∂Ω (e.g. ∂Ω is of class C 2+ǫ for some ǫ > 0) with outward unit normal − → n on ∂Ω. Moreover, the conditions imposed on the nonlinear reaction terms f i (u) and the nonnegative initial data u i,0 will be specified later.
The first part of this paper considers weak solutions to system (S). Our aim is to provide sufficient conditions on the porous medium exponents m i and on the nonlinearities f i (u), under which weak solutions are indeed bounded in L ∞ (and thus locally Hölder-continuous) for all times and grow at most polynomially in time. More precisely, we assume the following conditions on the nonlinearities: (i) The nonlinearities f i : R S → R are locally Lipschitz functions and satisfy |f i (u)| ≤ C(1 + |u| ν ), ∀u = (u 1 , . . . , u S ) ∈ R S , ∀i = 1, . . . , S, where R ∋ ν ≥ 1 is the maximal growth exponent of the reaction terms. (ii) There exist positive constants λ 1 , . . . , λ S > 0 such that: which formally implies the following mass dissipation law (iii) The nonlinearities are assumed quasi-positive, that is for all i = 1, . . . , S, holds f (u 1 , . . . , u i−1 , 0, u i+1 , . . . , u S ) ≥ 0, ∀u 1 , . . . , u S ≥ 0. (P) The quasi-positivity condition (P) ensures global nonnegativity of solutions subject to nonnegative initial data, see e.g. [Pie10,LP17].
The existence of global weak solutions to (S) subject to homogeneous Dirichlet boundary conditions and under the assumptions (G)-(M)-(P) was recently obtained in [LP17]. The proof of the following Theorem 1.1 on the existence of weak solutions to (S) subject to Neumann boundary conditions uses similar arguments to [LP17] and is postponed to Section 5. where the constant C depends on the L 2 -norm of the initial data, the constants λ i in (M), the diffusion coefficients d i > 0 and the domain Ω.
Remark 1.1. With a more careful analysis, it seems possible to generalise Theorem 1.1 and consider initial data u i,0 ∈ L 1 (Ω). We refer the interested reader to [PR16] for the case of systems with quadratic nonlinearities and L 1 initial data.
Given the weak solutions of Theorem 1.1, our aim is to establish their boundedness and a polynomially in time growing L ∞ -estimate under stronger assumptions on the porous medium exponents m i : First, we recall the a-priori estimate u i ∈ L mi+1 (Q T ) of Theorem 1.1 and the growth condition (G) imply f i (u) ∈ L 1+ǫ (Q T ) for some ǫ > 0, which also justifies the definition of weak solutions in Theorem 1.1. In fact, the L 1+ǫ integrability guarantees uniform integrability of nonlinearities in a suitable approximating scheme (see the proof of Theorem 1.1 in the Section 5).
Intuitively, Theorem 1.1 states that larger exponents m i yield higher integrability of the nonlinearities f i (u). Moreover, the functions u i solve a porous medium equation with the right hand side having higher integrability. Thus, by quantifying the smoothing effect from the porous medium equation, this allows to start a bootstrap argument, which eventually leads to boundedness of u i in L ∞ . In particular, it is of importance that our argument allows to show that the growth in time of the L ∞ -norms is at most polynomial. The first main result of this article is the following theorem.

Theorem 1.2 (Global bounded weak solutions).
Let Ω ⊂ R d be bounded with sufficiently smooth boundary. Let the initial data 0 ≤ u i,0 ∈ L ∞ (Ω), assume the conditions (G),(M) and (P) and m i > max{ν − 1; 1} for all i = 1 . . . S as required by Theorem 1.1. Finally, in dimensions d ≥ 3, we additionally assume Then, any weak solution of (S) obtained in Theorem 1.1 is bounded in L ∞ (Ω) and grows in time at most polynomially in the sense that, for any T > 0, where C T is a constant which depends at most polynomially on time. Consequently, these solutions are locally (in Q T ) Hölder continuous, see e.g. [Vaz07].
Remark 1.2 (Weakened assumptions on mass dissipation and initial data). If one is only interested in the boundedness of solutions but not in the polynomial growth of the L ∞ -norm, then the mass dissipation condition (M) can in fact be weakened to Also the assumed initial regularity u i,0 ∈ L ∞ (Ω) is not optimal and could be relaxed to L p integrability for sufficiently large p according to the details of the proof yet at the price of the readability of the Theorem.
Theorem 1.2 contributes to the large literature on global existence and boundedness of solutions to reaction-diffusion systems, which nevertheless poses still many open questions due to the lack of a unified approach (maximum principles do not hold for general systems). The largest part of the available literature, however, considers the case of linear diffusion, i.e. m i = 1 in system (S). We refer the reader to the extensive review of Michel Pierre [Pie10] and the references therein, in particular [Ba94, BR10, BP00, CV09, DFPV07, HLV98, HMP87, KK00, Laa11, Mas83, Mor89, Pie03,PS97] The case of nonlinear diffusion, on the other hand, is much less investigated. Most of the existing results considered special systems with special structures, see e.g. [Smo94,Leu09] . Up to the best of our knowledge, system (S) under the general structural assumptions (G)-(M)-(P) was only studied very recently in [LP17], where the authors showed the global existence of weak solutions. Therefore, the present paper serves as the first result to show the boundedness of weak solutions by assuming stronger conditions on porous medium exponents. Moreover, our proof allows to estimate explicitly the growth in time of the L ∞ -norm, which turns out to be essential in studying the large time behaviour of solutions in the following second part of the paper.
The second main result of this paper proves exponential convergence to equilibrium for a class of reactiondiffusion systems with porous media diffusion of the form (S), where the nonlinearities model the following reversible reaction with arbitrarily many chemical substances Here α i , β i ∈ [1, +∞) are the stoichiometric coefficients of the M + N involved substances A 1 , . . . , A M , B 1 , . . . , B N and k f , k b > 0 are the forward and backward reaction rate constants. For simplicity, yet without loss of generality, we assume k f = k b = 1. By applying mass action kinetics to (2) and by using the short notation a = (a 1 , . . . , a M ), b = (b 1 , . . . , b N ), α = (α 1 , . . . , α M ), β = (β 1 , . . . , β N ), we study the following reaction-diffusion system: x ∈ Ω. (R) Here d i , h j > 0 are diffusion coefficients and m i , p j > 1 are nonlinear diffusion exponents. It is clear that (R) is a special case of (S). It is also straightforward to verify condition (P), while condition (G) is satisfied by choosing, Finally condition (M) is a consequence from noting that After having the conditions (P), (G) and (M) verified, Theorem 1.1 implies the existence of global weak nonnegative solutions of system (R) provided Moreover by Theorem 1.2, these solutions are bounded in dimensions d = 1, 2, or in dimensions d ≥ 3 when additionally assuming for all i = 1 . . . M, j = 1 . . . N.
By multiplying the equations for a i and b j with β j and α i , respectively, and by adding the resulting terms, integration by parts with the homogeneous Neumann boundary conditions implies that these solutions satisfies the following mass conservation laws: amongst which exactly M + N − 1 linearly independent conservation laws ought to be selected and only the corresponding M + N − 1 components of the initial mass vector M ij need to be calculated from the initial data. System (R) possesses for each fixed positive initial mass vector (M ij ) a unique positive detailed balanced equilibrium (a ∞ , b ∞ ) = (a 1,∞ , . . . , a M,∞ , b 1,∞ , . . . , b N,∞ ) ∈ (0, ∞) M+N , which is the solutions of the following equilibrium equations: where we recall that the second line constitutes of only M + N − 1 linearly independent conditions.
To study the convergence to equilibrium for (R), we will use the so-called entropy method, which recently proved a highly suitable tool in the analysis of the large-time-behaviour of dissipative PDE systems. With respect to reaction-diffusion systems with linear diffusion, we refer in particular to [DF06,DF07,DF08,MHM15,DFT16,FT17a,FT17].
The key entropy functional (or in this case the free energy functional) of system (R) is defined by which dissipates according to the nonnegative entropy production functional, that is formally In the case of linear diffusion, i.e. m i = p j = 1 for all i = 1 . . . M, j = 1 . . . N , the convergence to equilibrium of solutions of (R) (or some special cases) was recently studied in e.g. [DF06,DF08,MHM15,FT17a,PSZ16]. Let us briefly review the entropy method used in the case of linear diffusion and then highlight the difficulties to be overcome in the current paper when dealing with nonlinear diffusion. In the case of linear diffusion, the entropy production writes as and the entropy method consists in establishing a functional inequality of the form for all functions a = (a i ), b = (b j ) satisfying the conservation laws (3). In order to do that, one first uses an additivity property of the relative entropy to calculate The term I 1 is controlled in terms of the entropy production D lin [a, b] thanks to the Logarithmic Sobolev Inequality (LSI) The remain term I 2 only involves the averages of the concentrations a i , b j and can be controlled by D lin [a, b] through lengthly, technical, but constructive estimates (see e.g. [FT17a,PSZ16] for more details). Note that this entropy approach applies successfully to more complex chemical reaction networks than (R), see [MHM15,DFT16,FT17,Mie]. We emphasised that the Logarithmic Sobolev Inequality (5) is not only used to control the term I 1 but also plays an important role in the estimates controlling the term I 2 .
In the case of nonlinear diffusion as here considered, we need a generalisation of the LSI (5) to exponents m i , p j ≥ 1. In this paper, we utilise the following generalisation (see e.g. [MM17]): for any m > (d − 2) + /d with (d − 2) + = max{d − 2; 0}, there exists a constant C(Ω, m) > 0 such that When m = 1, this coincides with the classical Logarithmic Sobolev inequality (5). For system (R), we have in particular Note that if we assume the averages a i and b j to be bounded below by a positive constant, then one can apply the same strategy as for the linear diffusion case in order to obtain the convergence to equilibrium. However, there is no chemical/physical reason for such a lower bound to hold in the transient behaviour of system (R) subject to general initial data. There are even perfectly admissible initial conditions, where some averages are zero since the corresponding species have not yet been formed.
To overcome this difficulty, we first observe that the mass conservation laws (3) subject to a positive mass vector M i,j > 0 implies that the averages a i and b j cannot be simultaneously small. Thus, at any fixed time at least one of the inequalities in (6) is useful, since either a i ≥ ε or b j ≥ ε for some suitably chosen ε > 0 depending on M i,j > 0. Secondly, we are able to compensate the still lacking lower bounds in (6) by a phenomena which can be called "indirect diffusion effect" and which means in our context that the reversible reaction (2) transfers diffusion from a species a i (with strictly positive diffusion bound in (6) due to a i ≥ ε) to other species b j (with lacking positive lower diffusion bound) in terms of a functional inequality, see Lemma 3.2 below.
Examples of indirect diffusion effect inequalities were already derived in e.g. [DF07,FLT17,FPT17], yet typically with a proof which requires uniform in time L ∞ -bounds on the solutions, which is a severe technical restriction as L ∞ -bounds for general reaction-diffusion systems are often unknown due to the lack of comparison principles. Note that also the L ∞ -bounds of Theorem 1.1 would be insufficient since polynomially growing and not uniform in time.
In this work, we are able to prove an indirect diffusion functional inequality without using any L ∞bounds on solutions but instead by exploiting the special structure of (R), see Lemma 3.2. Nevertheless, in the remaining part of applying the entropy method, the polynomial growth in time of the L ∞ -norm of Theorem 1.2 is still needed in one estimate concerning the relative entropy, yet the L ∞ -norm appears only within a logarithm. While it is unclear to us whether this is essential or just technical necessary in our approach, it allows to derive a time-dependent entropy-entropy production inequality (as a generalisation of the functional inequality (4)) of the form where the function Θ : R + → R + is of order 1/ ln(1 + T ) and satisfies +∞ 0 Θ(τ )dτ = +∞. Thus, a classical Gronwall argument implies explicit algebraic decay of E[a(T ), b(T )] − E[a ∞ , b ∞ ] to zero and thus, algebraic convergence to equilibrium in relative entropy.
To obtain exponential from algebraic decay, we show that after some sufficiently large time T 0 > 0, the averages a i (T ) and b j (T ) are bounded below by a positive constant for all T ≥ T 0 (since the equilibrium (a ∞ , b ∞ ) consists of positive constants). Hence, for T ≥ T 0 , we can use the inequalities (6) like in the case for systems with linear diffusion and obtain accordingly exponential convergence to equilibrium. Finally, since T 0 can be explicitly estimated, one recovers global exponential convergence to equilibrium (i.e. for all T ≥ 0) at the price of a smaller, yet explicit constant. Hence, the second main result of this paper is the following theorem.
Moreover, in dimensions d ≥ 3, we additionally assume Finally, consider a positive initial mass vector M ij > 0, which uniquely determines a positive equilibrium (a i∞ , b j∞ ) of system (R). Then, the bounded global weak solutions of Theorem 1.2 converge exponentially to where the constant C > 0 and the convergence rate λ p > 0 can be computed explicitly.

Notation:
• We denote by · the usual norm of L 2 (Ω). For other 1 ≤ p < +∞, we write · p as the norm of L p (Ω).
• For any T > 0, Q T = Ω × (0, T ) and L p (Q T ) =: L p (0, T ; L p (Ω)). The space-time norm is defined as usual • Throughout this work, we will denote by C T a generic positive constant which depends on certain parameters, and more importantly C T grows at most polynomially, i.e. there exists a polynomial P (x) such that C T ≤ P (T ) for all T > 0.
Organisation of the paper: Section 2 states the proof of Theorem 1.2. The proof of Theorem 1.3 is detailed in Section 3. This proof uses also a previously proven entropy-entropy production estimate for reaction-diffusion systems with linear diffusion, which is recalled in Section 4 for the sake of completeness. Finally, the existence of global weak solution is stated in Section 5.

Boundedness and local continuity of weak solutions
In this section, we prove for sufficiently large diffusion exponents m i that the weak solutions obtained in Theorem 1.1 are actually bounded in L ∞ and thus locally Hölder continuous. In Lemma 2.1, we device a bootstrap argument for the inhomogeneous porous media equation which proves that if the porous media exponents m i and the initial integrability are high enough, then the weak solutions of Theorem 1.1 satisfy an improve integrability in a space L s (Q T ) and the L s -norm grows at most polynomially in time T .

Lemma 2.1 (Smoothing effect of porous medium equation).
and subject to initial data u 0 ∈ L ∞ (Ω). Then, u satisfies d+2−2p0 , if p 0 < d+2 2 , and with a constant C T , which only depends on q, d, m, Ω and at most polynomially on T .
Remark 2.1. In the linear case m = 1 Lemma 2.1 recovers the corresponding regularity estimates of the heat equation, see [CDF14]. While the smoothing effect stated in Lemma 2.1 is certainly well-known, our main contribution here lies in the polynomial growth in time of the norms, which will be crucial in Section 3.
Proof. The idea of the proof of this lemma follows [CDF14, Lemma 3.3] and is divided into several steps.
It follows that Next, we construct a sequence p n ≥ 1 based on the estimate (14) and (19) such that and T 0 u rn sn dt ≤ D T,n with r n = m + p n − 1 and in which C T,n and D T,n are constants growing at most polynomially in T .

Like in
Step 2 in case m > 1 and γ < 1, we have by Young's inequality, analog to (17) while the case m = 1 and r n+1 /p n+1 = 1 follows without interpolation and the first term on the right-hand-side above. Combining (30), (31) and (32) yields T,n+1 T,n+1 Step 5. Passing to the limit as n → ∞. Considering the iteration (27), the only possible fixed point p ∞ of the sequence p n is Hence p ∞ < 0 if and only if p 0 > d+2 2 . In particular, it is straightforward to check that the sequence p n define by (27) is strictly monotone increasing if and only if either p n < p ∞ in the case p 0 < d+2 2 or p n > p ∞ in the case p 0 > d+2 2 when p ∞ < 0 holds or p 0 = d+2 2 where p ∞ = +∞. Therefore, we have as n → ∞ Step 6 (Interpolation). From (20) and (21) and by using the interpolation N pn+m−1 (Q T ) we get u ∈ L r (Q T ) for all r < ∞ in the case p 0 ≥ d+2 2 . In the case p 0 < d+2 2 , we obtain u ∈ L s (Q T ) for all This completes the proof of Lemma 2.1.
Lemma 2.2. Let u be a weak solution to (S) and where m = min{m i : i = 1 . . . S} and ν is defined in (G).
Then, it follows that u i L ∞ (QT ) ≤ C T for all i = 1 . . . S.
Proof. From u i ∈ L q0 (Q T ) for all i = 1, . . . , S, we have f i (u) ∈ L q0/ν (Q T ). Moreover note that the quasipositivity assumption (P) ensure non-negative solutions u for non-negative initial data u i,0 . Hence, the concentrations u i satisfy the (non-sign-changing) porous media equation since m ≤ m i . We then construct a sequence q n (equally for all i = 1, . . . , S) such that q n+1 = (md + 2)q n ν(d + 2) − 2q n for n ≥ 0.
Hence with q 0 > d(ν−m)+2(ν−1) 2 , after finitely many steps we arrive at q n > (d+2)ν 2 . From u i ∈ L s (Q T ) for all s < q n , we have in particular 2 (Q T ) for i = 1, . . . , S. By applying Lemma 2.1 once more we obtain u i ∈ L r (Q T ) ∩ L ∞ (0, T ; L q (Ω)) for all r, q < ∞. Thus, where the constant C T depends polynomially on T .
Remark 2.2. This regularity is well known for porous medium equation. However, since we were unable to find a reference, which includes the polynomial grow of the constant C T , we provide in the following a proof for the sake of completeness.
Hence, we estimate the solution of the inhomogeneous equation by using Duhamel's formula for all 0 ≤ t ≤ T , since α p < 1 by the assumption on p guarantees the convergence of the last integral on the right hand side. This finishes the proof.
Now we are ready to prove the boundedness of solutions to (S): Proof of Theorem 1.2. Assuming m i > ν − 1, the existence of weak solutions follows similar to [LP17,PR16] and is proven in Section 5 in detail. By the duality estimates in Lemma 5.1, we have u i ∈ L mi+1 (Q T ) for all i = 1, . . . , S.
Because m i > ν − 4 2+d it follows that Therefore, Lemma 2.2 yields u i ∈ L ∞ (Q T ) and u i L ∞ (QT ) ≤ C T for arbitrary T > 0, which shows that the weak solutions are bounded and the L ∞ (Ω) norms grows at most polynomially in time.
The local Hölder continuity of the bounded weak solutions is a classical result, see e.g. [DF85] or [Vaz07, Theorem 7.17].

Convergence to equilibrium
In this section, we prove exponential convergence to equilibrium of solutions to (R) by using the entropy method. We start by recalling the entropy (free energy) functional and its non-negative entropy production (free energy dissipation) functional D where we have used the short hand notation Moreover, the following additivity property of the relative entropy holds The first Lemma 3.1 of this section states the generalisation of the Logarithmic Sobolev Inequality, which shall use in our approach. where u = Ω udx.
Proof. The first inequality follows from [MM17]. The second estimate follows from an elementary inequality: The estimates in Lemma 3.1 constitute a generalisation of the Logarithmic Sobolev Inequality (5), which is recovered by setting m = 1 and for which the pre-factor u m−1 vanishes. In the case of porous media diffusion m > 1, the pre-factor u m−1 causes the lower bounds in Lemma 3.1 to degenerate for small spatial averages u. In particular, we have by Lemma 3.1 the following lower bound for the entropy production The problem of degeneracy appears when some averages a i or b j do not satisfy a positive lower bound. To overcome this problem, we first observe that due to the mass conservation laws (3) not all spatial averages can be small at the same time. If, for instance, a particular a i is sufficiently small (w.r.t. M ij ) then another b j can't be arbitrarily small because of a mass conservation law (3) connecting these two species, i.e.
The following crucial Lemma 3.2 shows functional inequalities, which quantity the so-called "indirect diffusion effect" and allows to compensate the lacking lower bounds for the species, whose spatial averages do not satisfy a lower bound.
We first introduce some convenient notations:

Moreover,
The conservation laws are now rewritten as Lemma 3.2 ("Indirect diffusion transfer" functional inequality). Let A i , B j : Ω → R + with i = 1 . . . M and j = 1 . . . N be nonnegative functions satisfying the conservation laws (37) and ε > 0 be a constant to be determined later. Assume that for some J ∈ {1, . . . , N }, B 2 j ≤ ε for all j = 1 . . . J. Then, there exists a constant K 1 which depends on ε such that: Remark 3.1. Note that when the last term on the left hand side A α − B β 2 diverges, the inequality holds trivially. Therefore, in the proof we only consider the case when it is finite.
Proof. Due to the mass conservation laws (37), we have the following natural bounds, . . , M, ∀j = 1, . . . , N for some constant M 0 > 0. Therefore, by Jensen's inequality, recalling that |Ω| = 1, From these bounds we get an upper bound for the right hand side of (38) We consider the following two cases. Case 1: If there exists i ∈ {1, . . . , M } such that δ i 2 ≥ ε or there exists a j ∈ {J + 1, . . . , N } such that η j 2 ≥ ε, we have: hence, the desired inequality (38) holds with K 1 = ε M 2 0 J . Case 2: Assume δ i 2 ≤ ε for all i ∈ {1, . . . , M } and η j 2 ≤ ε for all j ∈ {J + 1, . . . , N }, which together with the above assumption B 2 j ≤ ε and η 2 j ≤ B 2 j for all j = 1 . . . J implies η j 2 ≤ ε for all j ∈ {1, . . . , N }., Let λ > 0 and denote by Similarly we get, By Taylor's expansion, we have where the remainder terms R depends polynomially on A i and δ i . Note that |R( where we used δ i 2 ≤ ε in the last inequality. In order to estimate further, we use again Taylor's expansion where again, Q depends polynomially on B j , η j , which implies |Q(B j , η j )| ≤ C 1 (M 0 ) on G. Therefore, where we used that η j 2 ≤ ε for all j = 1, . . . , N . Combining these two estimates, we arrive at By Jensen's inequality and the assumption of the Lemma, we have On the other hand B j ≤ B 2 j ≤ M 0 , ∀j = J + 1, . . . , N . Thus, the conservation law (37) and δ i 2 ≤ ε yield Hence, by using |G| ≥ 1 2 we get from (39) that Because the right hand side of the above inequality converges to 1 Mi1 β1 αi as ε → 0, we can choose ε > 0 small enough, but still explicit, such that which implies the desired inequality (38) with the constant

Lemma 3.3 (An time-dependent entropy-entropy production estimate).
Let (a, b) = (a 1 , . . . , a M , b 1 , . . . , b N ) with a i , b j : Q T → R + be nonnegative functions, which satisfy the conservation laws (3). Moreover, Then, there exists a constant K 2 > 0 such that for all T > 0, Proof. Let ε > 0 be a small constant chosen in Lemma 3.2. We will consider two cases and for convenience we will drop T in a i (T ) and b j (T ) when there is no confusion. Case 1. Assume a i ≥ ε for all i = 1, . . . , M and b j ≥ ε for all j = 1, . . . , N . By applying Lemma 3.1, we have Using an entropy-entropy production inequality in case of system (R) with linear diffusion, see Lemma 4.1 below, we know that for an explicit constant K 4 > 0. Therefore, Case 2. Suppose either a i ≤ ε for some i ∈ {1, . . . , M } or b j ≤ ε for some j = 1, . . . , N . Due to the mass conservation laws β j a i + α i b j = M ij , it cannot happen that a i ≤ ε and b j ≤ ε simultaneously for a sufficiently small ε, e.g. ε < Mij 2 min 1 βj ; 1 αi . Therefore, without loss of generality, we can assume that b j ≤ ε ∀j = 1, . . . , J and b j ≥ ε ∀j = J + 1, . . . , N for some J ∈ {1, . . . , N }. Moreover, by mass conservation laws Thus, we can apply Lemma 3.1 to D[a, b] and estimate Applying Lemma 3.2 yields By using another functional inequality, which was already proven in the case of linear diffusion, see (48) in Section 4, we have Now, we estimate E[a, b] − E[a ∞ , b ∞ ] from above. Consider the two variables function which is continuous in (0, ∞) 2 and Φ(·, y) is increasing for each fixed y > 0. It holds that where in the last inequality we have used the estimates a i L ∞ (QT ) ≤ C T and b j L ∞ (QT ) ≤ C T and that C T is a constant growing at most polynomially w.r.t. T .
It's obvious that Q(A i ) ≥ 0 and moreover Therefore, and similarly Hence it follows from (41) that (42) A combination of (40) and (42) yields Finally, from Case 1 and Case 2, we can conclude the proof of Lemma 3.3 with Remark 3.2. The assumptions a i L ∞ (QT ) ≤ C T and b j L ∞ (QT ) ≤ C T in Lemma 3.3 are only needed to estimate E[a, b] − E[a ∞ , b ∞ ] above as in (41). In the case of linear diffusion, it is possible to avoid these L ∞ -bounds by using the additivity of the relative entropy (see also the proof of Lemma 4.1 in Section 4), i.e.
However, while for linear diffusion, the Logarithmic Sobolev Inequality controls to first part , such an estimate is unclear in the case of porous media diffusion, where the generalised Logarithmic Sobolev Inequality in Lemma 3.1 degenerates for states without lower bounds on the spatial averages.
We need also the following Csiszár-Kullback-Pinsker type inequality. The proof is standard and can be found in e.g. [DFT16,FT17a].
Lemma 3.4. There exists a constant C CKP > 0 such that for any measurable nonnegative functions a i , b j : Ω → R + satisfying the mass conservation (36), there holds We are ready to prove Theorem 1.3. By applying Lemma 3.3 this yields Moreover, due to the boundedness of solutions, we have the entropy-entropy production relation A classical Gronwall's inequality leads to By direct calculations Hence, and therefore thanks to the Csiszár-Kullback-Pinsker inequality in Lemma 3.4 which implies algebraic convergence to equilibrium of solutions to (R).
We will now show that from this it is possible to recover exponential convergence. Since the right hand side of (44) tends to zero as T → ∞, we can choose which implies for all t ≥ T 0 and thus Therefore, for all t ≥ T 0 , we can apply these lower bounds on the spatial averages bounds and Lemma 3.1 to estimate the entropy-entropy production as follows By applying again Lemma 4.1, we obtain which in a combination with the classical Gronwall's inequality yields for all t ≥ T 0 , where we used (43) for the second inequality. On the other hand, it follows from (43) that for all 0 ≤ t < T 0 , Due to the explicitness of T 0 in (45), we eventually get the exponential convergence with the constant C 2 = e λC1T0 and the rate λ = λC 1 . Note that C 2 is explicit since T 0 is explicit (see (45)).
With another application of the Csiszár-Kullback-Pinsker inequality in Lemma 3.4, this yields . Finally, by combining the above exponential L 1 -convergence with the at most polynomial grow L ∞ a-priori estimates a i L ∞ (QT ) , b j L ∞ (QT ) ≤ C T , interpolation yields for any 1 < p < ∞, for some 0 < λ p < λ(1 − θ) since C T grows at most polynomially in T , and similarly This concludes the proof of Theorem 1.3.

Entropy-entropy production Inequality
Lemma 4.1 (Entropy-entropy production estimate). Let a ∞ ∈ (0, ∞) M and b ∞ ∈ (0, ∞) N satisfy Then, there exists an explicit constant λ > 0 depending on a ∞ , b ∞ , α, β and the domain Ω, such that for any nonnegative functions a = (a i ) : for all i = 1, . . . , M, j = 1, . . . , N, the following entropy-entropy production inequality holds Remark 4.1. The above entropy-entropy production inequality was first proved in [FT17a] in a constructive way with explicit bounds on the constant λ. The proof stated here follows the line of a significantly simplified version presented in [FT17].
Proof. First, by the additivity of the relative entropy, we have It is straightforward that (I) can be controlled by D[a, b], i.e.
It remains to control (II). To do that, we first introduce the following useful notations and definitions In order to bound to estimate the right-hand-side of (46) with an upper bound of (II), we first observe from the conservation laws that there exists a constant M 0 > 0 such that Next, we note that the two variables function is continuous on (0, ∞) 2 and Φ(·, y) is increasing for each fixed y. Then, the term (II) is estimated as From (46) and (47), it remains to show that for some constant C 0 > 0. By using Lemma 4.2, we have with A = (A 1 , . . . , A M ) and B = (B 1 , . . . , B N ) for some constant C 1 > 0. Using the ansatz the right hand side of (48) writes as RHS of (48) = C 0 Moreover, the bounds for some constant M 1 > 0. From the ansatz (50) (and similar to the proof of Lemma 3.3), we have Next, we use Taylor expansion to estimate in which the Lagrange remainder term Q i = Q(µ i , δ i ) is uniformly bounded above by a constant for all admissible values of µ i and δ i thanks to the boundedness of µ i and δ i ≤ A 2 i ≤ M 0 . Similarly, with uniformly bounded remainder R j (ζ j , η j ). Thus is also uniformly bounded. Thus, by using (x+y) 2 ≥ 1 2 x 2 −y 2 and A α and the Cauchy-Schwarz inequality, Hence, for any δ ∈ (0, 1) holds (1 + ζ j ) βj 2 by choosing δ small enough such that 1 ≥ δ|Θ| 2 (M + N ) 2 since Θ is uniformly bounded above. This leads in combination with (49) to a lower bound of the left hand side of (48) From (51) and (54), it is sufficient to prove In order to do so, we note that the conservation laws rewritten in terms of the ansatz (50), i.e.
Lemma 5.1 (Duality estimates and uniform a-priori estimates for the approximating solutions, cf. [LP17]). Let u ε = (u 1,ε , . . . , u S,ε ) be the nonnegative solutions to the approximating system (58). Then, u i,ε L m i +1 (QT ) ≤ C for all T > 0 and i = 1, . . . , S, where the ε-independent constant C only depends on the L 2 -norm of the initial data, the positive constants λ i of assumption (M), the positive diffusion coefficients d i and the domain Ω. Moreover, we have f i,ε (u ε ) L 1+δ (QT ) ≤ C for some δ > 0, where the constant C depends only on the L 2 -norm of u i,0,ε , the positive constants λ i of assumption (M), the diffusion coefficients d i , the exponents m i and the domain Ω.
Proof. The proof follows [LP17] with straightforward changes due to the considered Neumann (instead of Dirichlet) boundary conditions. By setting Then, integration over (0, t) and multiplication with W (t) in L 2 (Ω) (due to the regularity of the approximative solutions) leads after integration over Ω to Therefore, by integrating (60) with respect to t on (0, T ), we obtain Moreover, we note that due to the nonnegativity of functions u i and the constant λ i . To estimate the right hand side of (61) in terms of the L 2 -norm of Z(0), we first notice from ∂ t Z − ∆W ≤ 0 that Multiplying this inequality with θ 0 in L 2 (Ω), where θ 0 ≥ 0 solves −∆θ 0 = Z(0), θ 0 | ∂Ω = 0, and using integration by parts − Ω θ 0 ∆ By inserting (62) and (63) into (61), we obtain which completes the proof of the first a-priori estimate of Lemma 5.1. Concerning the second uniform a-priori estimate for the nonlinearities, we have where C does not depend on ε. By the assumption m i > ν − 1 and the estimate of u i,ε L m i +1 (QT ) , we obtain f i,ε (u ε ) L 1+δ (QT ) ≤ C.
The following compactness lemma allows to extract a converging subsequence from the approximating system. Proof of Theorem 1.1. Thanks to the uniform bounds of the nonlinearities in Lemma 5.1 and the compactness Lemma 5.2, there exists a subsequence (not relabeled) {u i,ε } ε which converges in L 1 (Q T ) to limit functions u i ∈ L 1 (Q T ). From the L mi+1 -bound in Lemma 5.1, it holds in fact that u i,ε (up to another subsequence) converges strongly to u i in L mi (Q T ). For the nonlinearities, we first notice from Lemma 5.1 that the sequence {f i,ε (u ε )} is uniformly integrable. Moreover, for another subsequence u i,ε → u i a.e. in Q T it follows that Therefore, we can apply Vitali's Lemma, see e.g. [Sch05,Chapter 16], to obtain f i,ε (u ε ) → f i (u i ) strongly in L 1 (Q T ). All this allows to pass to the limit in the weak formulation (59) for any test function ψ ∈ |∇u mi i | β dxdt ≤ C(T, u i,0 1 , f i (u) L 1 (QT ) ) for all 1 ≤ β < 1 + 1 1 + m i d .