Global Renormalised Solutions and Equilibration of Reaction–Diffusion Systems with Nonlinear Diffusion

In this paper, we consider the evolution of concentrations \(u = (u_1, \ldots , u_I)\), \(I \in \mathbb N\), in a bounded domain \(\Omega \subset {\mathbb {R}}^d\) with Lipschitz boundary \(\partial \Omega \) subject to nonlinear diffusion and reactions as modelled by the parabolic system

Here and below, we set \(Q_T \mathrel {\mathop :}=\Omega \times (0, T)\) and \(\Gamma _T \mathrel {\mathop :}=\partial \Omega \times (0, T)\) for any \(T > 0\), and we employ \(\nu \) to denote the unit outward normal vector to \(\partial \Omega \). The diffusion constants \(d_i\) are assumed to satisfy \(d_i > 0\) throughout the paper. The admissible range of the nonlinear diffusion exponents \(m_i\) depends on the specific situation and will be clearly stated at every occurrence below. For instance, for proving the existence of global renormalised solutions, we demand \(m_i \in (0, 2)\), while the exponential equilibration of renormalised solutions is shown for \(m_i > \frac{(d-2)_+}{d}\) with \((d-2)_+ \mathrel {\mathop :}=\max \{0, d-2\}\). We will comment on the bounds on \(m_i\) later on. Moreover, we generally assume \(u_{i,0} \in L^1(\Omega )\), \(u_{i,0} \ge 0\), and impose the following assumptions on the reaction terms \(f_i(u)\), \(i=1,\ldots , I\):

1. (F1)

   (Local Lipschitz continuity) \(f_i: {\mathbb {R}}_+^I \rightarrow {\mathbb {R}}\) is locally Lipschitz continuous.

2. (F2)

   (Quasi-positivity) \(f_i(u)\ge 0\) for all \(u = (u_1, \ldots , u_I)\in {\mathbb {R}}_+^I\) with \(u_i = 0\).

3. (F3)

   (Entropy inequality) There exist constants \((\mu _i)_{i=1,\ldots , I} \in {\mathbb {R}}^I\) satisfying

   $$\begin{aligned} \sum _{i=1}^{I}f_i(u)(\log u_i+\mu _i) \le C \, \sum _{i=1}^{I}\left( 1 + u_i\log u_i\right) \end{aligned}$$

   (1.2)

   for all \(u\in {\mathbb {R}}_+^I\), where \(C>0\) is independent of u.

The second assumption (F2) guarantees non-negativity of solutions provided that the initial data are non-negative, which is a natural assumption since we consider here \(u_i\) as concentrations or densities. Assumption (F3), on the other hand, implies existence and control of the entropy

of the system, and (F3) is only slightly stronger than assuming control of the mass, i.e.

As it is one of the goals of this article to provide a generalisation of the results in Fischer (2015) to nonlinear diffusion models, we use the same entropy structure as in Fischer (2015), where Boltzmann entropy–dissipating systems are considered. Observing that (1.2) itself is just a condition on the involved reactions \(f_i(u)\), it is reasonable to employ the same entropy for both the linear and the nonlinear diffusive case. Other options for an entropy were shown to be feasible if one focuses on different kinds of models such as cross-diffusion population Chen et al. (2018); Chen and Jüngel (2019) or energy–reaction–diffusion systems Fischer et al. (2022); Hopf (2022).

In general, assumptions (F1), (F2), and either (F3) or (1.3) of the nonlinearities \(f_i\) are not enough to obtain suitable a priori estimates which allow to extend local strong solutions globally, see, e.g. Pierre and Schmitt (2000). On the one hand, for the case of linear diffusion (i.e. \(m_i=1\) for all \(i=1,\ldots , I\)), the global existence for (1.1) under extra assumptions on nonlinearities has been extensively investigated. For instance, global bounded solutions were shown under (F3) for the case \(d=1\) with cubic nonlinearities and \(d=2\) with quadratic nonlinearities, see Goudon and Vasseur (2010); Tang (2018). Recent works Souplet (2018); Caputo et al. (2019); Fellner et al. (2020) showed under either (F3) or (1.3) that bounded solutions can be obtained for (slightly super-)quadratic nonlinearities in all dimensions. We refer the reader to the extensive review Pierre (2010) for related results concerning global existence of bounded or weak solutions. It is remarked that in all such results, additional assumptions on nonlinearities must be imposed. A breakthrough has been made in Fischer (2015) where global renormalised solutions were obtained without any extra assumptions on the nonlinearities. The notion of a renormalised solution was initially introduced for Boltzmann’s equation in DiPerna and Lions (1989), and it has been applied subsequently to many other problems, see, e.g. Villani (1996); Dal Maso et al. (1999); Desvillettes et al. (2007). On the other hand, the case of nonlinear, possibly degenerate diffusion is much less understood. Up to our knowledge, there are only a few works (Laamri and Pierre 2017; Laamri and Perthame 2020 and Fellner et al. 2020) which showed global existence of (very) weak or bounded solutions for porous medium type of diffusion, i.e. \(m_i>1\), under (1.3) and some restricted growth conditions on the nonlinearities. However, renormalised solutions to a system as in (1.1) but imposing more regularity on the initial data and the underlying domain are provided in Lankeit and Winkler (2022). Instead of an entropy condition as in (1.2), the authors of Lankeit and Winkler (2022) assume an additional bound on cross-absorptive reaction terms. Furthermore, Fellner et al. (2020) proved the convergence to equilibrium for (1.1) modelling a single reversible chemical reaction. Recent results on the global existence of renormalised solutions also include reaction–cross-diffusion models in population dynamics Chen and Jüngel (2019) and energy–reaction–diffusion systems Fischer et al. (2022).

The present study shows the existence of global renormalised solutions to (1.1) featuring nonlinear diffusion of porous medium type or fast diffusion without any extra conditions on the nonlinearities (up to assuming (F1)–(F3)). In addition, we show that these solutions, in case (1.1) models general complex balanced chemical reaction networks, converge exponentially to equilibrium with explicitly computable rates. Our paper seems to be the first extensive contribution to the nonlinear diffusion-type system (1.1) establishing—for specific parameter regimes—the existence of global renormalised (or even weak and bounded) solutions as well as exponential convergence to equilibrium in various \(L^p\) spaces. A brief summary of our results is given in Table 1 at the end of this section.

The first part of this paper is concerned with the global existence of renormalised solutions to (1.1). We consider the following definition of renormalised solutions.

Definition 1.1

(Renormalised solutions) Let \(m_i \in [0,2]\) for all \(1 \le i \le I\) and \(u_0 = (u_{i,0})_i \in L^1(\Omega )^I\) be (componentwise) non-negative. A non-negative function \(u = (u_1, \ldots , u_I): \Omega \times (0,\infty ) \rightarrow {\mathbb {R}}_+^I\) is called a renormalised solution to (1.1) with initial data \(u_0\) if the following holds:

• \(u_i\in L^\infty _\textrm{loc}([0, \infty ), L^1(\Omega ))\) and \(\nabla {{u_i}^{\frac{m_i}{2}}}\in L^2_\textrm{loc}([0, \infty ), L^2(\Omega )^I)\),

• for any smooth function \(\xi : {\mathbb {R}}_+^I \rightarrow {\mathbb {R}}\) with compactly supported derivative \(D\xi \), any function \(\psi \in C^{\infty }(\overline{\Omega } \times [0,\infty ))\), and a.e. \(T>0\), we have

  $$\begin{aligned}{} & {} \int _{\Omega }\xi (u(\cdot , T))\psi (\cdot , T) \, \textrm{d}x - \int _{\Omega }\xi (u_0)\psi (\cdot , 0) \, \textrm{d}x - \int _{Q_T}\xi (u)\partial _t\psi \, \textrm{d}x \, \textrm{d}t\nonumber \\{} & {} \quad = -\sum _{i,j=1}^{I} d_i m_i \int _{Q_T}\psi \partial _i\partial _j\xi (u)u_i^{m_i-1}\nabla u_i\cdot \nabla u_j \, \textrm{d}x \, \textrm{d}t \nonumber \\{} & {} \qquad - \sum _{i=1}^{I} d_i m_i \int _{Q_T} \partial _i\xi (u) u_i^{m_i-1}\nabla u_i \cdot \nabla \psi \, \textrm{d}x \, \textrm{d}t \nonumber \\{} & {} \qquad + \sum _{i=1}^{I}\int _{Q_T}\partial _i\xi (u)f_i(u)\psi \, \textrm{d}x \, \textrm{d}t, \end{aligned}$$

  (1.4)

• for a.e. \(0\le s < T\), we have

  $$\begin{aligned} \begin{aligned} \sum _{i=1}^I\int _{\Omega }u_i(\log u_i +\mu _i - 1) \, \textrm{d}x\bigg |_s^T \le&-\sum _{i=1}^I\frac{4d_i}{m_i^2}\int _s^T\int _\Omega |\nabla (u_i)^\frac{m_i}{2}|^2 \, \textrm{d}x \, \textrm{d}t \\&+ \sum _{i=1}^I\int _s^T\int _{\Omega }f_i(u)(\log u_i + \mu _i) \, \textrm{d}x \, \textrm{d}t, \end{aligned} \end{aligned}$$

  (1.5)

• and for a.e. \(t>0\), we have

  $$\begin{aligned} \sum _{i=1}^I \int _{\Omega } q_iu_i(x,t) \, \textrm{d}x = \sum _{i=1}^I \int _{\Omega } q_iu_{i,0}(x) \, \textrm{d}x, \end{aligned}$$

  (1.6)

  provided there exist some numbers \(q_i \in {\mathbb {R}}\) such that \(\sum _{i=1}^I q_if_i(z) = 0\) for all \(z\in {\mathbb {R}}^I\).

Remark 1.2

• We notice that all integrals in (1.4) are well defined due to the compact support of \(D\xi \) and as the integrals relating to gradients of solutions can be rewritten as

  $$\begin{aligned}{} & {} \int _{Q_T}\psi \partial _i\partial _j\xi (u)u_i^{m_i-1}\nabla u_i\cdot \nabla u_j \, \textrm{d}x \, \textrm{d}t\\{} & {} \quad = \frac{4}{m_im_j}\int _{Q_T}\psi \partial _i\partial _j\xi (u)u_i^{\frac{m_i}{2}}u_j^{1-\frac{m_j}{2}}\nabla u_i^{\frac{m_i}{2}} \cdot \nabla u_j^{\frac{m_j}{2}} \, \textrm{d}x \, \textrm{d}t \end{aligned}$$

  and

  $$\begin{aligned} \int _{Q_T}\partial _i\xi (u)u_i^{m_i-1}\nabla u_i\nabla \psi \, \textrm{d}x \, \textrm{d}t = \frac{2}{m_i} \int _{Q_T}\partial _i\xi (u)u_i^{\frac{m_i}{2}}\nabla u_i^{\frac{m_i}{2}} \nabla \psi \, \textrm{d}x \, \textrm{d}t \end{aligned}$$

  which are both finite thanks to the boundedness of \(\nabla u_i^{\frac{m_i}{2}}\) in \(L^2(Q_T)\) and the restriction \(m_i \in [0,2]\) for all \(1 \le i \le I\). Note that \(m_i < 0\) or \(m_i > 2\) is not admissible as the integrals above may degenerate in this case even for small values of \(u_i\).

• In the case of linear diffusion, the weak entropy law (1.5) and conservation laws (1.6) in Definition 1.1 are redundant since they can be proved using (1.4), see Fischer (2017). Unfortunately, this seems to be impossible in our present case due to the degeneracy of the diffusion. Nevertheless, if the renormalised solution turns out to be bounded, one can prove these laws directly (see the Proof of Theorem 1.7). Proving (1.5) and (1.6) using only (1.4) in the case of nonlinear diffusion remains an interesting open problem.

The first main result of this article is the following theorem. We stress that even though Definition 1.1 allows for \(m_i \in [0,2]\), the construction of a renormalised solution demands \(m_i \in (0,2)\) for all \(i = 1, \ldots , I\).

Theorem 1.3

(Global existence of renormalised solutions) Let \(\Omega \subset {\mathbb {R}}^d\) be a bounded domain with Lipschitz boundary \(\partial \Omega \). Assume \(m_i \in (0,2)\) for all \(i=1,\ldots , I,\) and that conditions (F1)–(F3) hold. Then, for any non-negative initial data \(u_0 = (u_{i,0})_i \in L^1(\Omega )^I\) subject to

there exists a global renormalised solution to (1.1).

Remark 1.4

• Relation (1.6) usually corresponds to mass conservation laws in chemical reactions, see, e.g. Fellner and Tang (2018).

• The uniqueness of renormalised solutions is widely open. For a weaker notion called weak–strong uniqueness, i.e. the renormalised solution is unique as long as a strong solution exists, we refer the interested reader to Fischer (2017).

• The fact that (1.5) and (1.6) are satisfied by any renormalised solution plays an important role in the next part of this paper where we show that all renormalised solutions, instead of only some of them, converge to equilibrium.

• Assumption (1.7) on the initial data \(u_0 \in L^1(\Omega )^I\) is used on a technical level at the end of the proof of Lemma 2.4 when showing that the candidate \(u = (u_i)_i\) for a renormalised solution satisfies \({u_i}^\frac{m_i}{2} \in L^2(0, T; H^1(\Omega ))\) for all \(i \in \{1, \dotsc , I\}\). Moreover, the additional regularity (1.7) carries over to the renormalised solution (cf. (1.5)) and, in particular, guarantees that the initial entropy (cf. (1.13)) is finite.

• A similar situation is considered in Laamri and Pierre (2017) with homogeneous Dirichlet conditions. In fact, (Laamri and Pierre 2017, Theorem 2.6) is comparable to our Theorem 1.3 but instead of (1.7), the authors assume that an a priori \(L^1\) estimate is available.

The strategy for constructing a global renormalised solution follows the ideas in Fischer (2015). Given a sequence of approximate solutions \(u^\varepsilon \), one proves compactness of the family \(u^\varepsilon \) by deriving bounds on truncations \(\varphi _i^E(u^\varepsilon )\), \(E \in \mathbb N\), \(i = 1, \dotsc , I\), and subsequently employing an Aubin–Lions lemma. The smooth functions \(\varphi _i^E\) defined in (2.6) essentially truncate the mappings \(u \mapsto u_i\) if u becomes too large (measured in terms of E) while leaving \(u \mapsto u_i\) unchanged for sufficiently small u. The main difficulty in our current situation is caused by the different diffusion exponents \(m_i\) for each species \(u_i\), which makes the truncations \(\varphi _i^E\) provided in Fischer (2015) not applicable. We overcome this issue by modifying the functions \(\varphi _i^E\) in such a way that the truncation is not decided by \(\sum _j u_j\) but by a weighted sum \(\sum _j E^{-\alpha ^i_j} u_j\) with suitable constants \(\alpha ^i_j \ge 1\) chosen depending on the diffusion exponents \(m_i\). We then pass to the limit \(\varepsilon \rightarrow 0\) in the equation for \(\varphi _i^E(u^\varepsilon )\) at the cost of an additional defect measure which, nevertheless, vanishes in the subsequent limit \(E \rightarrow \infty \). Finally, an equation for \(\xi (\varphi _i^E(u))\) is derived, where \(\xi \) is subject to the same assumptions as in Definition 1.1, and the limit \(E \rightarrow \infty \) is performed resulting in the desired equation (1.4) for \(\xi (u)\).

The restriction \(m_i \in (0,2)\) for \(i=1,\ldots , I\) is, on the one hand, necessary within the construction of the truncations \(\varphi _i^E\), see, e.g. (2.7), (2.9), and (2.10). On the other hand, it is crucial in showing the vanishing of the defect measure when establishing the renormalised solution. We do not know whether this restriction is purely technical or due to a deeper reason. It is remarked that the same condition is heavily used in Laamri and Pierre (2017), where global weak solutions to (1.1) are investigated. In particular, models with logarithmic diffusion \(m_i = 0\) and ultra-fast diffusion processes are not included in the subsequent results. Systems of porous medium type are, nevertheless, covered to a large extent within Theorem 1.7.

In the second part of this paper, we study the large-time behaviour of reaction–diffusion systems of type (1.1) modelling chemical reaction networks subject to the complex balance condition. The main vocabulary is introduced below but we refer to Desvillettes et al. (2017) for a more detailed discussion of the involved concepts. We assume that there are I chemicals \(S_1, \ldots , S_I\) reacting via the following R reactions

where \(y_{r,i}, y_{r,i}' \in \{0\}\cup [1,\infty )\) are stoichiometric coefficients, and \(k_r>0\) are reaction rate constants. Utilising the notation \(y_r = (y_{r,i})_{i=1,\ldots , I}\) and \(y_r' = (y_{r,i}')_{i=1,\ldots , I}\), we can rewrite the reactions in (1.8) as

Denote by \(u_i(x,t)\) the concentration of \(S_i\) at position \(x\in \Omega \) and time \(t\ge 0\). The reaction network (1.8) results in the following reaction–diffusion system with nonlinear diffusion for \(u = (u_1, \ldots , u_I)\):

in which the reaction term \(f_i(u)\) is determined using the mass action law,

It is obvious that \(f_i(u)\) is locally Lipschitz continuous; hence, (F1) is satisfied. For (F2), it is easy to check that \(f_i(u) \ge 0\) holds true for all \(u\in {\mathbb {R}}_{+}^I\) satisfying \(u_i = 0\) since we assume \(k_r>0\) and \(y\in \mathbb (\{0\}\cup [1,\infty ))^I\) for all \(y\in \{y_r, y_r' \}_{r=1,\ldots , R}\). Before verifying (F3), we recall the definition of a complex balanced equilibrium. A constant state \(u_\infty = (u_{i,\infty })_i \in [0,\infty )^I\) is called a complex balanced equilibrium for (1.10)–(1.11) if

Intuitively, this condition means that for any complex the total outflow from the complex and the total inflow into the complex are balanced at such an equilibrium. Using complex balanced equilibria, one can show that the nonlinearities in (1.11) satisfy (F3) by choosing \(\mu _i = -\log u_{i,\infty }\) with a strictly positive complex balanced equilibrium \(u_\infty \) (see Desvillettes et al. 2017, Proposition 2.1). Therefore, one can apply Theorem 1.3 to obtain global renormalised solutions to (1.10) and (1.11).

In general, there might exist infinitely many solutions to (1.12). To uniquely identify a positive equilibrium, we need a set of conservation laws. More precisely, we consider the Wegscheider matrix (or stoichiometric coefficient matrix)

and define \(m = \dim (\ker (W))\). If \(m>0\), then we can choose a matrix \(\mathbb Q \in {\mathbb {R}}^{m\times I}\), where the rows form a basis of \(\ker (W)\). Note that this implies \(\mathbb Q\, (y_r' - y_r) = 0\) for all \(r=1,\ldots , R\). Therefore, by recalling \(f(u) = (f_i(u))_{i=1,\ldots , I} = \sum _{r=1}^{R}k_ru^{y_r}(y_r' - y_r)\) we obtain

By using the homogeneous Neumann boundary condition, we can then (formally) compute that

and, consequently, for a.e. \(t > 0\),

In other words, systems (1.10) and (1.11) possess m linearly independent mass conservation laws. By the vocabulary used in the literature on chemical reaction network theory, for any fixed non-negative mass vector \(M\in {\mathbb {R}}^m\), the set \(\{u\in {\mathbb {R}}^I_+: \mathbb Q\,u=M\}\) is called the compatibility class corresponding to M. In case that for each strictly positive vector \(M \in {\mathbb {R}}^m\), there exists a strictly positive complex balanced equilibrium, one says that the chemical reaction network is complex balanced. This terminology is well defined since it was proved in Horn (1972/73) that if one equilibrium is complex balanced, then all equilibria are complex balanced. It was further shown that for each non-negative initial mass \(M\in {\mathbb {R}}^m\) with \(M \ne 0\), there exists a unique strictly positive complex balanced equilibrium \(u_{\infty } = (u_{i,\infty })_i \in {\mathbb {R}}^I_+\) in the compatibility class corresponding to M (cf. Horn and Jackson 1972; Feinberg 1979).

We stress that additionally there might exist many so-called boundary equilibria, which are complex balanced equilibria lying on \(\partial {\mathbb {R}}_+^I\). It is remarked that there exists a large class of complex balanced chemical reaction networks, called concordant networks, which do not have boundary equilibria, see Shinar and Feinberg (2013). For the sake of brevity, from now on we call the unique strictly positive complex balanced equilibrium simply the complex balanced equilibrium. Note that all the previous considerations—although established for ODE models for chemical reaction networks—also apply to our PDE setting thanks to the homogeneous Neumann boundary conditions, which allow for spatially homogeneous equilibria.

Theorem 1.5

(Exponential equilibration of renormalised solutions) Let \(\Omega \subset {\mathbb {R}}^d\) be a bounded domain with Lipschitz boundary \(\partial \Omega \). Assume that the reaction network (1.9) is complex balanced, suppose that (F1)–(F3) hold, and let any non-negative initial data \(u_0 = (u_{i,0})_i \in L^1(\Omega )^I\) with \(\sum _{i=1}^I\int _{\Omega }u_{i,0}\log u_{i,0} \, \textrm{d}x < \infty \) be given.

• If \(m_i \in (0,2)\) for all \(1 \le i \le I\), then there exists a global renormalised solution to (1.10)–(1.11).

• If \(m_i > \frac{(d-2)_+}{d}\) for all \(1 \le i \le I\) and assuming that there exist no boundary equilibria, all renormalised solutions to (1.10)–(1.11) converge exponentially to the positive equilibrium \(u_\infty \in {\mathbb {R}}^I_+\) in the same compatibility class as \(u_0\) with a rate which can be explicitly computed up to a finite dimensional inequality, i.e.

  $$\begin{aligned} \sum _{i=1}^I\Vert u_i(t) - u_{i,\infty }\Vert _{L^1(\Omega )} \le Ce^{-\lambda t} \quad \text { for all } \quad t\ge 0, \end{aligned}$$

  for positive constants \(C>0, \lambda > 0\) depending on the initial entropy, the domain \(\Omega \), the diffusion coefficients \(d_i\) and exponents \(m_i\), and the vectors of stoichiometric coefficients \(\{y_r, y_r'\}_{r=1,\ldots , R}\).

We emphasise that the convergence to equilibrium in Theorem 1.5 is proved for all renormalised solutions, not only the ones obtained via an approximation procedure. This is possible because the proof of the convergence uses only the weak entropy–entropy dissipation law (1.5) and the conservation laws (1.6) which are satisfied by all renormalised solutions. It is noted that, when a system possesses boundary equilibria, the convergence to the positive equilibrium (more precisely, the instability of boundary equilibria) is very subtle and strongly connected to the famous Global Attractor Conjecture, see, for instance, Desvillettes et al. (2017); Fellner and Tang (2018); Craciun (2015).

Theorem 1.5, up to our knowledge, is the first result of trend to equilibrium for general complex balanced reaction networks with nonlinear diffusion. A special case of a single reversible reaction was considered recently in Fellner et al. (2020) in which the authors utilised a so-called indirect diffusion effect (see, e.g. Einav et al. 2020), which seems difficult to be generalised to general systems. The proof of Theorem 1.5 uses the relative entropy

and its corresponding entropy dissipation

where \(\Psi (x;y) = x\log (x/y) - x + y\). It holds, at least formally, that

along any trajectory of (1.10) and (1.11). The degeneracy of the nonlinear diffusion makes the classical logarithmic Sobolev inequality not applicable. Our main idea here is to utilise some generalised logarithmic Sobolev inequalities (see Lemmas 3.3 and 3.4), which are suited for nonlinear diffusion, and the established results in Fellner and Tang (2018) for the case of linear diffusion, to firstly derive an entropy–entropy dissipation inequality of the form

where \(\alpha = \max _{i=1,\ldots , I}\{1,m_i\}\). Note that this functional inequality is proved for all non-negative functions \(u:\Omega \rightarrow {\mathbb {R}}_+^I\) satisfying the conservation laws \(\mathbb Q\int _{\Omega }u \, dx = |\Omega |\mathbb Q u_\infty \), and therefore is suitable for renormalised solutions, which have very low regularity. If \(\alpha = 1\) in (1.15), one immediately gets exponential convergence of the relative entropy to zero and, consequently, exponential convergence of the solutions to equilibrium in \(L^1\) thanks to a Csiszár–Kullback–Pinsker-type inequality. If \(\alpha > 1\), we first obtain an algebraic decay of the relative entropy to zero. Thanks to this, we can explicitly compute a finite time \(T_0>0\) from which onwards (in time) the averages of concentrations are strictly bounded below by a positive constant. This helps to compensate the degeneracy of the diffusion and, therefore, to show that solutions with such lower bounds satisfy the linear entropy–entropy dissipation inequality, i.e.

which recovers exponential convergence to equilibrium. \(\square \)

The use of renormalised solutions allows to deal with a large class of nonlinearities, but on the other hand it restricts Theorem 1.5 to the case \(\frac{(d-2)_+}{d}< m_i < 2\). When the nonlinearities are of polynomial type, it was shown in Laamri and Pierre (2017); Fellner et al. (2020), under the assumption of the mass dissipation (1.3) instead of the entropy inequality (1.2), that one can get global weak or even bounded solutions if the porous medium exponents \(m_i\) are large enough. In the following theorem, we show an analogue result under the entropy inequality condition (1.2). Moreover, the solution is shown to converge exponentially to the positive complex balanced equilibrium.

Definition 1.6

(Weak solutions) A vector of non-negative functions \(u = (u_1,\ldots , u_I): \Omega \times (0,\infty ) \rightarrow {\mathbb {R}}_+^I\) is called a (global) weak solution to (1.1) if

and

for all \(T > 0\), \(i=1, \ldots , I\), and test functions \(\psi \in L^{\infty }(Q_T)\cap L^{2(m_i+1)} (0,T;W^{1,2(m_i+1)}(\Omega ) )\).

Theorem 1.7

(Existence and convergence of weak and bounded solutions) Let \(\Omega \subset {\mathbb {R}}^d\) be a bounded domain with Lipschitz boundary \(\partial \Omega \). Assume that the assumptions (F1)–(F3) hold and that the nonlinearities are bounded by polynomials, i.e. there exist \(C>0\) and \(\rho _j \ge 1\) for all \(j = 1, \ldots , I\) such that

• If

  $$\begin{aligned} m_i \ge \rho _i - 1 \quad \text { for all } \quad i=1,\ldots , I, \end{aligned}$$

  (1.17)

  then system (1.1) has a global non-negative weak solution for each non-negative initial data \(u_0 = (u_{i,0})_{i=1,\ldots , I}\) such that \(u_{i,0}\log u_{i,0} \in L^1(\Omega )\cap H^{-1}(\Omega )\) for all \(i=1,\ldots , I\).

• Under an even stronger condition with \(\rho \mathrel {\mathop :}=\max _i \rho _i\), namely

  $$\begin{aligned} \begin{aligned} \min _{i=1,\ldots , I} m_i&> \max \{\rho -1, 1\} \quad \text { for } \quad d \le 2, \\ \min _{i=1,\ldots , I} m_i&> \max \Big \{\rho - \frac{4}{d+2}, 1\Big \} \quad \text { for } \quad d\ge 3, \end{aligned} \end{aligned}$$

  (1.18)

  the weak solution to (1.1) subject to bounded non-negative initial data \(u_0\in L^\infty (\Omega )^I\) is locally bounded in time, i.e.

  $$\begin{aligned} \Vert u_i\Vert _{L^\infty (\Omega \times (0,T))} \le C_T \quad \text { for all } \quad i=1,\ldots , I, \end{aligned}$$

  where \(C_T\) is a constant growing at most polynomially in \(T>0\).

Assume now that the reaction network (1.9) is complex balanced and define

where \(|y_{r}| = \sum _{j=1}^Iy_{r,j}\). Assume moreover that (1.9) has no boundary equilibria.

• If (1.17) holds and \(0\le u_0\) such that \({u_{i,0}\log {u_{i,0}} \in L^1(\Omega )\cap H^{-1}(\Omega )}\) for all \(i=1,\ldots , I\), all bounded weak solutions (locally in time) to (1.10) and (1.11) converge exponentially to the complex balanced equilibrium in \(L^1(\Omega )\), i.e.

  $$\begin{aligned} \sum _{i=1}^{I}\Vert u_{i}^{\varepsilon }(t) - u_{i,\infty }\Vert _{L^1(\Omega )} \le Ce^{-\lambda t} \quad \text { for all } \quad t \ge 0, \end{aligned}$$

  for some \(C>0, \lambda > 0\) depending on the initial entropy, the domain \(\Omega \), the diffusion coefficients \(d_i\) and exponents \(m_i\), and the vectors of stoichiometric coefficients \(\{y_r, y_r'\}_{r=1,\ldots , R}\).

• If (1.18) holds and \(0\le u_0\in L^\infty (\Omega )^I\), weak solutions to (1.10) and (1.11) are bounded locally in time and converge exponentially to equilibrium in any \(L^p\) norm with \(p \in [1,\infty )\), i.e.

  $$\begin{aligned} \sum _{i=1}^I\Vert u_i(t) - u_{i,\infty }\Vert _{L^p(\Omega )} \le C_pe^{-\lambda _p t} \quad \text { for all } \quad t \ge 0, \end{aligned}$$

  (1.20)

  with positive constants \(C_p, \lambda _p>0\) depending on p, the initial entropy, the domain \(\Omega \), the diffusion coefficients \(d_i\) and exponents \(m_i\), and the vectors of stoichiometric coefficients \(\{y_r, y_r'\}_{r=1,\ldots , R}\).

Remark 1.8

In addition to Remark 1.4, we point out that another reason for assuming \(u_{i,0} \log u_{i,0} \in L^1(\Omega ) \cap H^{-1}(\Omega )\) for all \(i \in \{1, \dotsc , I\}\) (rather than demanding just \(u_0 \in (L^1(\Omega ) \cap H^{-1}(\Omega ))^I\) as in Laamri and Perthame 2020) is the supposed entropy inequality (1.2). More precisely, we need this hypothesis on the initial data in the proof of Lemma 3.6, which is itself a crucial ingredient of the proof of Theorem 1.7.

As one can see from Theorems 1.5 and 1.7, the global existence and trend to equilibrium for (1.10) and (1.11) are well established for either \(\frac{(d-2)_+}{d}< m_i < 2\) for all \(i=1,\ldots , I\), or \(m_i\) large enough in the sense of (1.17) or (1.18) for all \(i=1,\ldots , I\). There exists, therefore, a gap in which the global existence of any kind of solution remains open. Remarkably, the proof of the convergence to equilibrium in this paper does not rely on the restriction \(m_i < 2\), and it therefore is applicable to any global solution to (1.10) and (1.11) as long as it satisfies the entropy law (1.5) and the mass conservation laws (1.6). We, thus, arrive at the following result on the convergence to equilibrium for approximating systems of (1.10) and (1.11) uniformly in the approximation parameter.

Theorem 1.9

(Uniform convergence of approximate solutions) Let \(\Omega \subset {\mathbb {R}}^d\) be a bounded domain with Lipschitz boundary \(\partial \Omega \). Assume that (1.10) and (1.11) subject to (F1)–(F3), non-negative initial data \(u_0 = (u_{i,0})_i \in L^1(\Omega )^I\) with \(\sum _{i=1}^I\int _{\Omega }u_{i,0}\log u_{i,0} \, dx < \infty \), and \(m_i > \frac{(d-2)_+}{d}\) for all \(i = 1, \dotsc , I\) admits a complex balanced equilibrium but no boundary equilibria. For any \(\varepsilon >0\), let \(u^\varepsilon = (u^{\varepsilon }_i)_{i=1,\ldots , I}\) be the solution to the approximate system

where \(u_{i,0}^\varepsilon \in L^\infty (\Omega )\) satisfies

Moreover, the approximated initial data \(u_0^\varepsilon \) are chosen such that, for all \(\varepsilon >0\),

(This implies that system (1.21) admits a unique positive complex balanced equilibrium \(u_\infty \in {\mathbb {R}}_+^I\), which is independent of \(\varepsilon >0\).)

Then, \(u^\varepsilon \) converges exponentially to \(u_\infty \) with a rate which is uniform in \(\varepsilon \), i.e.

where \(C>0, \lambda >0\) are constants depending on the initial entropy, the domain \(\Omega \), the diffusion coefficients \(d_i\) and exponents \(m_i\), and the vectors of stoichiometric coefficients \(\{y_r, y_r'\}_{r=1,\ldots , R}\) but independent of \(\varepsilon >0\).

To prove Theorem 1.9, we make use of the same relative entropy as in (1.13),

where \(u^\varepsilon \) is the solution to (1.21). The corresponding non-negative entropy dissipation can be calculated as

Since there is no uniform-in-\(\varepsilon \) \(L^\infty \) bound for \(u^\varepsilon \) available, the factor \(\frac{1}{1+\varepsilon |f(u^{\varepsilon })|}\) is not bounded below uniformly in \(\varepsilon >0\). Therefore, it seems to be impossible to use the entropy dissipation (for limit solutions) in (1.14) as a lower bound for \(\mathcal D[u^{\varepsilon }]\). We overcome this issue by using the ideas in Fellner and Tang (2018). Roughly speaking, we deal with the second sum in \(\mathcal D[u^{\varepsilon }]\) by estimating it below by a term involving only spatial averages of \(u^{\varepsilon }\) (rather than \(u^{\varepsilon }\) pointwise as above), and then exploiting the fact that these averages are bounded uniformly in \(\varepsilon >0\).

We summarise the global existence and convergence to equilibrium for the mass action system (1.10) and (1.11) in Table 1 (assuming the existence of a unique strictly positive complex balanced equilibrium and the absence of boundary equilibria). Note that in our framework the global existence of any kind of solution remains open as soon as there exist \(i, j \in \{1,\ldots ,I\}\) such that the corresponding diffusion exponents satisfy \(m_i \ge 2\) and \(m_j < \rho - 1\), where \(\rho \) is defined in (1.19). In order to keep the presentation simple, we assume within Table 1 that all diffusion exponents \(m_i = m\) coincide.

The rest of this paper is organised as follows. In Sect. 2, we show the global existence of renormalised solutions for the general system (1.1). The proofs of Theorems 1.5, 1.7, and 1.9 on the convergence to equilibrium of chemical reaction networks are presented in Sect. 3. Some technical proofs of auxiliary results are postponed to Appendix.

2.1 Existence of Approximate Solutions

Let \(m_i \in (0,2)\) for all \(i \in \{1,\dotsc ,I\}\) and \(\varepsilon >0\) and consider the following approximating system for \(u^{\varepsilon } = (u_1^{\varepsilon }, \ldots , u_I^{\varepsilon })\),

where \(u_{i,0}^{\varepsilon }\in L^{\infty }(\Omega )\) is non-negative and \(u_{i,0}^{\varepsilon } \xrightarrow {\varepsilon \rightarrow 0} u_{i,0}\) in \(L^1(\Omega )\). Moreover, we demand

With this approximation, it is easy to check that the approximated nonlinearities \(f_i^\varepsilon \) still satisfy assumptions (F1)–(F3).

Lemma 2.1

Provided \(m_i \in (0,2)\) for all \(i \in \{1,\dotsc ,I\}\) and \(\varepsilon >0\), there exists a global weak solution \(u^{\varepsilon } = (u_i^{\varepsilon })_{i} \in L^\infty _\textrm{loc}([0, \infty ), L^\infty (\Omega ))^I\) to the approximate system (2.1) with \(\nabla (u^{\varepsilon }_i)^{\frac{m_i}{2}} \in L^2_\textrm{loc}([0,\infty ), L^2(\Omega ))\) and \(\partial _tu_i^\varepsilon \in L^2_\textrm{loc}([0, \infty ), (H^1(\Omega ))')\). In detail,

for all \(\psi \in L^2_\textrm{loc}([0,\infty ), H^1(\Omega ))\), \(i \in \{1,\dotsc ,I\}\), and a.e. \(0 \le t_0 < t_1\). Moreover, each \(u_i^\varepsilon \) is a.e. non-negative and

for a.e. \(0 \le t_0 < t_1\).

Proof

For each \(\varepsilon >0\), we see by recalling \(|f(u^{\varepsilon })| = \sum _{i=1}^{I}|f_i(u^{\varepsilon })|\) that

The existence of bounded weak solutions to (2.1) is therefore standard. However, since we are not able to find a precise reference, a proof is given in Appendix 1. By multiplying (2.1) by \(\log (u_i^{\varepsilon }) + \mu _i\) (or more rigorously by \(\log (u_i^\varepsilon +\delta ) + \mu _i\) for some \(\delta >0\), then let \(\delta \rightarrow 0\)), summing the resultants over \(i=1,\ldots , I\), and then integrating over \(\Omega \), we obtain

where we used (F3) and \(x\le \delta x\log x + C_\delta \) for all \(x\ge 0\) and any \(\delta > 0\) at the last step. Hence, by integrating (2.5) over (t, T) and using Gronwall’s inequality, we obtain the desired estimate (2.4). The uniform bound on \(u_i^{\varepsilon }\log u_i^{\varepsilon }\) in \(L^{\infty }(0,T;L^1(\Omega ))\) and \(|\nabla (u_i^{\varepsilon })^{\frac{m_i}{2}}|\) in \(L^2(0,T;L^2(\Omega ))\) follow immediately from (2.4). \(\square \)

2.2 Existence of Renormalised Solutions

As it can be seen from Lemma 2.1, the a priori estimates of \(u_i^{\varepsilon }\) are not enough to extract a convergent subsequence. Following the idea from Fischer (2015), we consider another approximation of \(u_i^{\varepsilon }\) by defining

for \(E \in \mathbb N\), and a smooth function \(\xi : {\mathbb {R}} \rightarrow [0,1]\) satisfying \(\xi \equiv 1\) on \((-\infty , 0)\) and \(\xi \equiv 0\) on \((1,\infty )\). The constants \(\alpha ^i_j\) are given by (recall that \(m_j \in (0,2)\) for all \(j \in \{1,\dotsc ,I\}\))

Remark 2.2

The \(\alpha ^i_j\) are chosen in (2.7) for the sake of simplicity. In fact, we can choose any \(\alpha ^i_j\) such that

See the proof of (E2) in Lemma 2.3.

Lemma 2.3

The smooth truncations \(\varphi _i^E\) defined in (2.6) have the following properties:

1. (E1)

   \(\varphi _i^E\in C^2({\mathbb {R}}_+^I)\).

2. (E2)

   For any \(i\in \{1,\ldots , I\}\), there exists a constant \(K_i>0\) such that

   $$\begin{aligned} \max _{1\le j,k\le I}\sup _{E\ge 1}\sup _{v\in {\mathbb {R}}_+^I}v_j^\frac{m_j}{2} v_{k}^{1 - \frac{m_k}{2}} |\partial _j \partial _k \varphi _i^E(v)| \le K_i. \end{aligned}$$

   (2.8)
3. (E3)

   For every E, the set \(\textrm{supp}\ D\varphi _i^E\) is bounded.

4. (E4)

   For all \(j\in \{1,\ldots , I\}\) and all \(v\in {\mathbb {R}}_+^{I}\), there holds \(\lim _{E\rightarrow \infty }\partial _j\varphi _i^E(v)= \delta _{ij}\).

5. (E5)

   There exists a constant \(K>0\) such that

   $$\begin{aligned} \max _{1\le j\le I}\sup _{E\ge 1}\sup _{v\in {\mathbb {R}}_+^I}|\partial _j\varphi _i^E(v)| \le K. \end{aligned}$$
6. (E6)

   \(\varphi _i^E(v) = v_i\) for any \(v\in {\mathbb {R}}_+^I\) with \(\sum _{j=1}^{I} v_j E^{-\alpha ^i_j} \le 1\).

7. (E7)

   For every \(K>0\) and every j, k, there holds

   $$\begin{aligned} \lim _{E\rightarrow \infty }\sup _{|v| \le K}|\partial _j\partial _k\varphi _i^E(v)| = 0. \end{aligned}$$
8. (E8)

   \(\varphi _i^E(v) = v_i\) for any \(v\in {\mathbb {R}}_+^I\) with \(\sum _{j=1}^{I} \varphi _j^E(v) E^{-\alpha ^i_j} \le 1\).

Proof

Properties (E1), (E3), (E4), and (E6) are immediate. To prove (E2), we first compute

for \(j \ne i\), where for brevity we write \((\cdots )\) instead of \(\big (\sum _{j=1}^I E^{-\alpha ^i_j} v_j - 1 \big )\). From that, one further gets the second derivatives

for \(j \ne i\) and \(k \ne i\).

To show (2.8), we will consider the following cases:

• When \(i = j = k\), we have

  $$\begin{aligned} v_i^{\frac{m_i}{2}}v_i^{1-\frac{m_i}{2}}|\partial _i^2 \varphi _i^E(v)|&= v_i|\partial _i^2\varphi _i^E(v)|\\&\le |v_i||v_i-3E|E^{-2\alpha _i^i}|\xi ''(\cdots )| + 2v_iE^{-\alpha _i^i}|\xi '(\cdots )|. \end{aligned}$$

  Employing the identities \(\xi '(\cdots ) = \xi ''(\cdots ) = 0\) for \(v_i \ge 2 E^{\alpha _i^i}\) as well as the bound \(|\xi '(\cdots )| + |\xi ''(\cdots )| \le C\), we can further estimate

  $$\begin{aligned} v_i^{\frac{m_i}{2}}v_i^{1-\frac{m_i}{2}}|\partial _i^2 \varphi _i^E(v)|&\le CE^{-\alpha _i^i}|E^{\alpha _i^i} + E| + C. \end{aligned}$$

  Therefore, by choosing \(\alpha _i^i = 1\), we have the desired estimate (2.8) for \(i=j=k\).

• When \(k=i\) and \(j\ne i\), we estimate using \(\alpha _i^i = 1\), \(\xi '(\cdots ) = \xi ''(\cdots ) = 0\) when \(v_j \ge 2E^{\alpha ^i_j}\) or \(v_i \ge 2E\), and \(|\xi '(\cdots )| + |\xi ''(\cdots )| \le C\),

  $$\begin{aligned} \qquad&v_j^{\frac{m_j}{2}}v_i^{1-\frac{m_i}{2}}|\partial _j\partial _i\varphi _i^E(v)|\\&\le v_j^{\frac{m_j}{2}}v_i^{1-\frac{m_i}{2}}|\xi '(\cdots )|E^{-\alpha ^i_j} + v_j^{\frac{m_j}{2}}v_i^{1-\frac{m_i}{2}}|\xi ''(\cdots )||v_i - 3E|E^{-\alpha ^i_j - 1}\\&\le CE^{\frac{m_j}{2}\alpha ^i_j}E^{1-\frac{m_i}{2}}|\xi '(\cdots )|E^{-\alpha ^i_j} + CE^{\frac{m_j}{2}\alpha ^i_j}E^{1-\frac{m_i}{2}}|\xi ''(\cdots )||E + 3E|E^{-\alpha ^i_j - 1}\\&\le CE^{\alpha ^i_j(\frac{m_j}{2}-1) + (1-\frac{m_i}{2})}. \end{aligned}$$

  Thus, (2.8) is proved in the case \(k=i\ne j\) if we choose \(\alpha ^i_j\) such that

  $$\begin{aligned} \alpha ^i_j(\frac{m_j}{2}-1) + (1-\frac{m_i}{2}) \le 0 \quad \text { or equivalently } \quad \alpha ^i_j \ge \frac{2-m_i}{2-m_j}. \end{aligned}$$

  (2.9)

• When \(j=i\) and \(k\ne i\), we estimate similarly to the previous cases

  $$\begin{aligned} v_i^{\frac{m_i}{2}}v_k^{1-\frac{m_k}{2}}|\partial _i\partial _k\varphi _i^E(v)| \le CE^{\frac{m_i}{2} - \alpha ^i_k\frac{m_k}{2}}. \end{aligned}$$

  If we choose

  $$\begin{aligned} \frac{m_i}{2} - \alpha ^i_k\frac{m_k}{2} \le 0 \quad \text { or equivalently } \quad \alpha ^i_k \ge \frac{m_i}{m_k}, \end{aligned}$$

  (2.10)

  then obviously (2.8) holds true for \(j=i\ne k\).

• Finally, when \(j\ne i\) and \(k\ne i\), we estimate

  $$\begin{aligned} v_j^{\frac{m_j}{2}}v_k^{1-\frac{m_k}{2}}|\partial _j\partial _k\varphi _i^E(v)|&\le CE^{\alpha ^i_j \frac{m_j}{2}}E^{\alpha ^i_k(1-\frac{m_k}{2})}|\xi ''(\cdots )|E^{-\alpha ^i_j - \alpha ^i_k}E\\&\le CE^{\alpha ^i_j(\frac{m_j}{2}-1) - \alpha ^i_k \frac{m_k}{2}+1}. \end{aligned}$$

  It is easy to see that from (2.9) and (2.10) it follows

  $$\begin{aligned} \alpha ^i_j\Big (\frac{m_j}{2}-1\Big ) - \alpha ^i_k \frac{m_k}{2}+1 \le 0, \end{aligned}$$

  and hence (2.8) is proved in this case.

From (2.9) and (2.10), we see that \(\alpha ^i_j > 1\) for all \(j\ne i\); hence, (E5) and (E7) follow.

For (E8), it is enough to show that \(\sum _{j=1}^{I} \varphi _j^E(v) E^{-\alpha ^i_j} \le 1\) yields \(\sum _{j=1}^{I} v_j E^{-\alpha ^i_j} \le 1\). Set

From \(\varphi _j^E(v) = (v_j - 3E)\xi (y) + 3E\) and \(\sum _{j=1}^{I} \varphi _j^E(v) E^{-\alpha ^i_j} \le 1\), we deduce

Due to \(\xi \le 1\) and \(E z \ge 1\), we are now led to

This entails \(y \le 0\) and, thus, completes the proof of the lemma. \(\square \)

Lemma 2.4

Consider non-negative functions \(u_0 = (u_{i,0})_i \in L^1(\Omega )^I\) which satisfy

Let \(u^\varepsilon = (u_1^\varepsilon , \dotsc , u_I^\varepsilon )\) for \(\varepsilon \rightarrow 0\) be the sequence of solutions to the regularised problems as stated in Lemma 2.1.

Then, there exists a subsequence \(u^\varepsilon \) converging a.e. on \(\Omega \times [0,\infty )\) to a limit \(u \in L^\infty _\textrm{loc}([0,\infty ), L^1(\Omega ))^I\) with \(u_i \log u_i \in L^\infty _\textrm{loc}([0,\infty ), L^1(\Omega ))\) for all \(i \in \{1, \dotsc , I\}\). Furthermore, \((u_i^\varepsilon )^\frac{m_i}{2} \rightharpoonup {u_i}^\frac{m_i}{2}\) weakly in \(L^2(0,T; H^1(\Omega ))\) for all \(i \in \{1, \dotsc , I\}\) and \(T>0\).

Proof

Due to the lack of uniform-in-\(\varepsilon \) estimates of the nonlinearities, it is difficult to show directly, for instance, by means of an Aubin–Lions lemma, that \(u_i^{\varepsilon }\) has a convergent subsequence. Following the ideas from Fischer (2015), we first prove that \(\varphi ^E_i(u^{\varepsilon })\) converges (up to a subsequence) to \(z_i^E\) as \(\varepsilon \rightarrow 0\), and then that \(z_i^E\) converges to \(u_i\) (up to a subsequence) as \(E\rightarrow \infty \). In combination with the convergence of \(\varphi _i^E(u^{\varepsilon })\) to \(u^{\varepsilon }_i\) for \(E\rightarrow \infty \), this leads to the desired result.

Due to the bound \(\varphi _i^E \le 3E\), it follows at once that \(\{\varphi _i^E(u^\varepsilon )\}\) is bounded in \(L^2(0,T; L^2(\Omega ))\) uniformly in \(\varepsilon >0\) for each \(E \in \mathbb N\). Next, by the chain rule we have

Thus,

thanks to the compact support of \(D\varphi _i^E\), \(m_j< 2\), and the \(L^2(0,T;L^2(\Omega ))\) bound on \(\nabla (u_j^{\varepsilon })^{\frac{m_j}{2}}\) from Lemma 2.1. As a consequence, we know that \(\{\varphi _i^E(u^\varepsilon )\}\) is bounded (uniformly in \(\varepsilon \)) in \(L^{2}(0,T;H^{1}(\Omega ))\).

To apply the Aubin–Lions Lemma, we need an estimate concerning the time derivative of \(\varphi _i^{E}(u^{\varepsilon })\). Using \(\partial _j \varphi _i^{E}(u^{\varepsilon }) \psi \) as a test function in (2.3) and summing over \(j\in \{1,\ldots , I\}\), we obtain for almost all \(t_2 > t_1 \ge 0\),

The third term on the right-hand side is clearly bounded uniformly in \(\varepsilon \) for each fixed \(E \in \mathbb N\) since \(D\varphi _i^E\) has a compact support. The first and second terms are bounded uniformly in \(\varepsilon \) thanks to the boundedness of \(\nabla (u_k^\varepsilon )^{\frac{m_k}{2}}\) in \(L^2(0,T;L^2(\Omega ))\) for all \(k=1,\ldots , I\), and properties (E2), (E3) in Lemma 2.3. It follows then that \(\partial _t\varphi _i^{E}(u^{\varepsilon })\) is bounded uniformly (w.r.t. \(\varepsilon >0\)) in \(L^1(0,T;(W^{1,\infty }(\Omega ))')\) for each fixed \(E \in \mathbb N\).

Therefore, by applying an Aubin–Lions lemma to the sequence \(\{\varphi _i^E(u^\varepsilon )\}_{\varepsilon >0}\), for fixed \(E \in \mathbb N\), there exists a subsequence (not relabeled) of \(\varphi _i^{E}(u^{\varepsilon })\) converging strongly in \(L^2(0,T; L^2(\Omega ))\) and, thus, almost everywhere as \(\varepsilon \rightarrow 0\). Using a diagonal sequence argument, one can extract a further subsequence such that \(\varphi _i^{E}(u^{\varepsilon })\) converges a.e. in \(Q_T\) to a measurable function \(z_i^{E}\) for all \(E \in \mathbb N\) and \(i = 1, \dotsc , I\).

We next prove that \(z_i^E\) converges a.e. to some measurable function \(u_i\) as \(E\rightarrow \infty \). First, since \(\sum _{j=1}^{I}u_j^{\varepsilon }\log u_j^{\varepsilon }\) is uniformly bounded w.r.t. to \(\varepsilon >0\) in \(L^{\infty }(0,T;L^1(\Omega ))\), it follows that \(\varphi _i^E(u^{\varepsilon })\log \varphi _i^{E}(u^{\varepsilon })\) is uniformly bounded w.r.t. to \(\varepsilon >0\) and \(E \in \mathbb N\) in \(L^{\infty }(0,T;L^1(\Omega ))\). This is trivial in case \(E \le 1\), \(\varphi _i^E(u^\varepsilon ) \le 1\), or \(\sum _{j=1}^I u_j^\varepsilon E^{-\alpha ^i_j} \le 1\). Otherwise, there exists some j such that \(u_j^\varepsilon E^{-\alpha ^i_j} > 1/I\) holds true. And as \(\alpha ^i_j > 1\), we derive

The bound \(\varphi _i^E(u^\varepsilon ) > 1\) ensures \(\varphi _i^E(u^\varepsilon ) \log \varphi _i^E(u^\varepsilon ) \le 3 I u_j^\varepsilon \log ( 3 I u_j^\varepsilon )\) and the claim follows. Thus, \(z_i^E\log z_i^E\) is uniformly bounded in \(L^{\infty }(0,T;L^1(\Omega ))\) w.r.t. \(E \in \mathbb N\) by Fatou’s Lemma. Secondly, \(z_j^E(x,t) = z_j^{\widetilde{E}}(x,t)\) for all \(\widetilde{E}>E\), if \(\sum _{j=1}^{I}z_j^E(x,t)E^{-\alpha _j^i} < 1\). Indeed,

guarantees that \(\sum _{j=1}^{I}\varphi _j^{E}(u^{\varepsilon }(x,t))E^{-\alpha _j^i} < 1\) holds true for sufficiently small \(\varepsilon >0\). Thanks to property (E8) in Lemma 2.3, it follows that \(\varphi _i^{E}(u^{\varepsilon }(x,t)) = u^{\varepsilon }_i(x,t) = \varphi _i^{\widetilde{E}}(u^{\varepsilon }(x,t))\) for small enough \(\varepsilon \) and all \(\widetilde{E} > E\). Therefore, due to \(z_i^E(x,t) = \lim _{\varepsilon \rightarrow 0}\varphi _i^E(u^{\varepsilon }(x,t))\), we obtain \(z_i^{E}(x,t) = z_i^{\widetilde{E}}(x,t)\) for all \(\widetilde{E}>E\) as desired. As a result, if \(\sum _{j=1}^{I}z_j^E(x,t)E^{-\alpha _j^i} < 1\) for some (x, t), then \(\lim _{E\rightarrow \infty }z_i^E(x,t)\) exists and is finite for all \(i=1,\ldots , I\). Using the fact that \(\sum _{j=i}^{I}z_j^E\log z_j^E\) is bounded uniformly w.r.t. \(E \in \mathbb N\) in \(L^{\infty }(0,T;L^1(\Omega ))\), we have

where \(\mathcal L^{n+1}\) is the Lebesgue measure in \({\mathbb {R}}^{n+1}\). Hence, the limit \(u_i(x,t) = \lim _{E\rightarrow \infty }z_i^E(x,t)\) exists for a.e. \((x,t)\in Q_T\). Moreover, \(u_i\log u_i \in L^{\infty }(0,T;L^1(\Omega ))\) due to Fatou’s lemma and \(z_i^E\log z_i^E\in L^{\infty }(0,T;L^1(\Omega ))\). Since \(\sum _{j=1}^{I}u_j^{\varepsilon }\) is uniformly bounded in \(L^{1}(Q_T)\), we find

and this limit is uniform in \(\varepsilon >0\). Now, we estimate for \(\delta >0\)

As we have proved that \(\varphi _i^E(u^{\varepsilon }) \rightarrow z_i^E\) a.e. in \(Q_T\) for \(\varepsilon \rightarrow 0\) and fixed \(E \in \mathbb N\) and that \(z_i^{E} \rightarrow u_i\) a.e. in \(Q_T\) as \(E\rightarrow \infty \), we infer the convergence in measure of \(u_i^\varepsilon \) to \(u_i\) for \(\varepsilon \rightarrow 0\), and convergence a.e. of another subsequence.

The uniform bound on \(u_i^\varepsilon \log u_i^\varepsilon \) in \(L^\infty (0, T; L^1(\Omega ))\) and the convergence \(u_i^{\varepsilon } \rightarrow u_i\) a.e. ensure that \(u_i^\varepsilon \rightarrow u_i\) strongly in \(L^p(0,T; L^1(\Omega ))\) for all \(p \ge 1\). One can easily prove this by truncating \(u_i^\varepsilon \) at a sufficiently large threshold. By the strong convergence of \(u_i^\varepsilon \) to \(u_i\) in \(L^1(0,T; L^1(\Omega ))\), we are able to prove the distributional convergence of \((u_i^\varepsilon )^\frac{m_i}{2}\) to \((u_i)^\frac{m_i}{2}\) by noting that \(m_i<2\). The uniform \(L^2(0,T; H^1(\Omega ))\) bound on \((u_i^\varepsilon )^\frac{m_i}{2}\) then leads to a subsequence \((u_i^\varepsilon )^\frac{m_i}{2}\) which converges to \((u_i)^\frac{m_i}{2}\) weakly in \(L^2(0,T; H^1(\Omega ))\). Thanks to the weak lower semicontinuity of the \(L^2(0,T; L^2(\Omega ))\) norm, we also deduce

for all \(i \in \{1,\dotsc ,I\}\), where the right-hand side is bounded via the uniform \(L^2(0,T; H^1(\Omega ))\) bound on \((u_i^\varepsilon )^\frac{m_i}{2}\). This completes the proof of the lemma. \(\square \)

At this point, we expect the function \(u = (u_1, \ldots , u_I)\) in Lemma 2.4 to be a renormalised solution to (1.1). We will first derive an equation admitting \(\varphi _i^E(u)\) as a solution. This equation is already “almost” matching the formulation for a renormalised solution (1.4) except for a “defect measure”.

Lemma 2.5

Let \(u = (u_1, \ldots , u_I)\) be the functions constructed in Lemma 2.4. Then, for any \(\psi \in C^{\infty }([0,T], C_0^{\infty }(\overline{\Omega }))\), the truncated functions \(\varphi _i^E(u)\), \(i \in \{1, \dotsc , I\}\), satisfy

where \(\mu _i^E\) is a signed Radon measure satisfying

for all \(T>0\) and \(i \in \{1, \ldots , I\}\).

Proof

Choosing \(t_1 = 0\), \(t_2 = T\), and using the convergence \(u^{\varepsilon } \rightarrow u\) a.e. in \(Q_T\), the weak convergence \(\nabla (u_i^\varepsilon )^\frac{m_i}{2} \rightharpoonup \nabla {u_i}^\frac{m_i}{2}\) in \(L^2(0, T; L^2(\Omega ))\), and the fact that \(D \varphi _i^E(u^{\varepsilon })\) vanishes when \(u^{\varepsilon }\) is too large, we straightforwardly obtain the convergence of the left-hand side and the last two terms on the right-hand side of (2.11). It remains to establish

with a signed Radon measure \(\mu _i^E\) satisfying (2.13). By denoting

we can use property (E2) of the truncation function \(\varphi _i^E\) and the \(\varepsilon \)-uniform bound of \(\nabla (u_i^{\varepsilon })^{\frac{m_i}{2}}\) in \(L^2(0,T;L^2(\Omega ))\), to obtain

for all \(\varepsilon >0\). Therefore, by passing to a subsequence, we know that \(\mu _{i,\varepsilon }^E\) weak-\(*\) converges on \(\overline{Q}_T\) to a signed Radon measure \(\mu _i^{E}\) as \(\varepsilon \rightarrow 0\). It remains to prove (2.13). Due to Young’s inequality, we have

where

We stress that \(\nu _{l,K}^{\varepsilon }\) is uniformly bounded w.r.t. l, K, and \(\varepsilon \). Consequently, we may pass to a subsequence \(\nu _{l,K}^{\varepsilon }\) which converges weak-\(*\) on \(\overline{Q}_T\) to a Radon measure \(\nu _{l,K}\) as \(\varepsilon \rightarrow 0\). Together with the weak-\(*\) convergence of \(\mu _{i,\varepsilon }^E\) to \(\mu _i^E\) on \(\overline{Q}_T\), we derive

Moreover,

which is bounded uniformly in \(\varepsilon \). Fatou’s Lemma applied to the counting measure on \(\mathbb N\) now entails

Therefore, employing the dominated convergence theorem, we can finally estimate

thanks to properties (E2) and (E7) of the truncation function \(\varphi _i^E\). \(\square \)

To prove Theorem 1.3, we use the following technical lemma whose proof can be found in Fischer (2015).

Lemma 2.6

(Fischer 2015, Lemma 4)(A weak chain rule for the time derivative) Let \(\Omega \) be a bounded domain with Lipschitz boundary. Assume that \(T>0\), \(v\in L^2(0,T;H^1(\Omega )^I)\), and \(v_0\in L^1(\Omega )^I\). Let \(\nu _i\) be a Radon measure on \(\overline{Q}_T\), \(w_i \in L^1(Q_T)\), and \(z_i\in L^2(0,T;L^2(\Omega )^I)\) for \(1\le i\le I\). Assume moreover that for any \(\psi \in C^{\infty }(\overline{Q}_T)\) with compact support we have

Let \(\xi :{\mathbb {R}}^I \rightarrow {\mathbb {R}}\) be a smooth function with compactly supported first derivatives. Then, we have for all \(\psi \in C^{\infty }(\overline{Q}_T)\) with compact support

We are now ready to prove the existence of global renormalised solutions.

Proof of Theorem 1.3

By applying Lemma 2.6 to (2.12), we obtain an approximate relation for the weak time derivative of \(\xi (\varphi ^E(u))\), which leads to the desired equation for \(\xi (u)\) when passing to the limit \(E \rightarrow \infty \). In detail, we set \(v_i \mathrel {\mathop :}=\varphi _i^E(u)\), \((v_0)_i \mathrel {\mathop :}=\varphi _i^E(u_0)\), \(\nu _i \mathrel {\mathop :}=-\mu _i^E\),

and obtain for any smooth function \(\xi \) with compactly supported first derivatives

We now want to pass to the limit \(E\rightarrow \infty \) in (2.14) to obtain (1.4). Note first that due to (2.13), we obtain in the limit an equality instead of just an estimate. Since \(\xi \) has compactly supported first derivatives, \(\partial _j\varphi _i^E(u)\) is bounded (see (E5)), and \(\varphi _i^E(u) \rightarrow u_i\) a.e. as \(E\rightarrow \infty \), we directly obtain the following limits for the first line in (2.14):

By employing \(m_j u_j^{m_j-1} \nabla u_j = 2 {u_j}^\frac{m_j}{2} \nabla {u_j}^\frac{m_j}{2}\) and \({u_j}^\frac{m_j}{2} \in L^2(0, T; H^1(\Omega ))\), we also derive

It remains to ensure the convergence of the third and fourth line of (2.14). To this end, we recall \({u_j}^\frac{m_j}{2} \in L^2(0, T; H^1(\Omega ))\), and we utilise the following observation: there exists a constant \(E_0>0\) such that for all \(E > E_0\) the inequality \(\sum _{i=1}^I u_i \ge E_0\) implies \(\partial _i \xi (\varphi ^E(u)) = \partial _i \xi (u) = 0\) and \(\partial _i \partial _k \xi (\varphi ^E(u)) = \partial _i \partial _k \xi (u) = 0\) for all \(i,k \in \{1,\dotsc ,I\}\). The convergence of the third and fourth line above is now a consequence of this auxiliary result, as the derivatives of \(\xi \) are zero provided \(\max _i u_i\) is larger than \(E_0\):

The previous observation readily follows by choosing \(E_0>0\) such that \(\textrm{supp}(D\xi ) \subset B_{E_0/\sqrt{I}}(0)\). Given \(E > E_0\) and \(\sum _{i=1}^I u_i \ge E_0\), one gets \(\sum _{i=1}^{I} \varphi _i^E(u) \ge E_0\) from the definition of \(\varphi _i^E\) in (2.6). In particular, \(u, \varphi ^E(u) \notin B_{E_0/\sqrt{I}}(0)\), which proves the claim. \(\square \)

Finally, we show that the constructed renormalised solution satisfies the weak entropy law (1.5) and, when it is applicable, the conservation laws (1.6).

Proof of (1.5) and (1.6)

Since \(u^{\varepsilon }_i \rightarrow u_i\) a.e. in \(Q_T\) and \(f_i\) are continuous, we have

We rewrite (2.5) as

We now integrate this relation on (s, T) and use (2.15) and Fatou’s lemma for the first integral on the right-hand side, and the uniform-in-\(\varepsilon \) in Lemma 3.6 for the second integral to obtain

which is (1.5). To show (1.6), we take \(t_0 = 0\) and \(\psi = q_i\) in (2.3) and sum over i to get

Letting \(\varepsilon \rightarrow 0\) proves (1.6). \(\square \)

3.1 Entropy–Entropy Dissipation Inequality

Without loss of generality, we assume that the domain \(\Omega \) has unit volume, i.e. \(|\Omega | = 1\). We employ the following notation:

To show the convergence to equilibrium for (1.10), we exploit the so-called entropy method. Assuming the complex balanced condition, system (1.10) possesses the relative entropy functional

and the corresponding entropy dissipation function

with \(\Psi (x;y) = x\log (x/y) - x + y\) (see (Desvillettes et al. 2017, Proposition 2.1) for a derivation of \(\mathcal D[u]\)). Formally, we have

along the trajectory of (1.10). The following lemma shows that the \(L^1\) norm can be bounded by the relative entropy.

Lemma 3.1

Assume that \(\mathcal {E}[u|u_{\infty }]<+\infty \). Then,

for \(i=1,\ldots , I\).

Proof

The elementary inequalities \(x\log (x/y) - x + y \ge (\sqrt{x} - \sqrt{y})^2 \ge \frac{1}{2} x - y\) ensure that

Taking the non-negativity of the solution u into account allows to conclude. \(\square \)

Lemma 3.2

(Fellner and Tang 2017, Lemma 2.3) (Csiszár–Kullback–Pinsker-type inequality) Fix a non-negative mass \(M \in {\mathbb {R}}_{+}^m\). For all measurable functions \(u: \Omega \rightarrow {\mathbb {R}}_+^I\) satisfying the mass conservation \(\mathbb Q\, \overline{u} = M\), the following inequality holds,

in which \(C_\textrm{CKP} > 0\) is a constant depending only on \(\Omega \) and \(u_{\infty }\).

The following two lemmas provide generalisations of the classical logarithmic Sobolev inequality, which are suited for nonlinear diffusion.

Lemma 3.3

Dolbeault et al. (2008)(\(L^m\)-logarithmic Sobolev inequality) For any \(m \ge 1\), there exists a constant \(C(\Omega ,m)>0\) such that we have

Lemma 3.4

Mielke and Mittnenzweig (2018) (Generalised logarithmic Sobolev inequality) For any \(m > \frac{(d-2)_+}{d}\), there exists a constant \(C(\Omega , m)>0\) such that

The cornerstone of the entropy method is the entropy–entropy dissipation inequality which is proved in the following lemma.

Lemma 3.5

Let \(K>0\) and \(m_i > \frac{(d-2)_+}{d}\) for all \(1 \le i \le I\). There exist constants \(C>0\) and \(\alpha \ge 1\) depending on K such that for any measurable function \(u: \Omega \rightarrow {\mathbb {R}}_+^I\) subject to \(\mathbb Q\overline{u} = \mathbb Q u_\infty \) and \(\mathcal {E}[u|u_\infty ]\le K\), we have

Proof

First, similar to Lemma 3.1, we have the following \(L^1\) bounds:

For all i with \(\frac{(d-2)_+}{d}< m_i < 1\), we apply Lemma 3.4 to have

since \(m_i < 1\) and \(\overline{u}_i \le M_1\). Therefore, by applying Lemma 3.3, we can estimate the entropy dissipation from below as follows:

For simplicity, we denote \(\alpha = \max _{i=1,\ldots , I}\{1, m_i\}\),

If \(D_i \ge 1\) for some \(i \in \{1,\ldots , I\}\) or \(F\ge 1\), then we have

using \(\mathcal E[u|u_{\infty }]\le K\). It remains to consider the case \(D_i \le 1\) for all \(i=1,\ldots , I\) and \(F \le 1\). Here, we find

By applying now the entropy–entropy production estimate for the case of linear diffusion (cf. Fellner and Tang 2018, Theorem 1.1), we have

and, thus,

From (3.4) and (3.6), we obtain the desired estimate. \(\square \)

3.2 Convergence of Renormalised Solutions

Proof of Theorem 1.5

The global existence of renormalised solutions for (1.10) follows from Theorem 1.3. Indeed, assumptions (F1) and (F2) can be easily checked using (1.11), and (F3) follows from the identity (cf. Desvillettes et al. 2017, Proposition 2.1)

recalling \(\Psi (x;y) = x\log (x/y) - x + y\).

Thanks to (1.6) and (1.5), we have the conservation laws

and

We now apply Lemma 3.5 to get

If \(\alpha = 1\), we immediately get the exponential convergence

for some constant \(\lambda >0\). For \(\alpha > 1\), we first obtain the algebraic decay

Indeed, we define

It follows from (3.7) that

Differentiating \(\varphi (s)\) with respect to s gives

and consequently,

Applying a nonlinear Gronwall inequality yields

thanks to \(\psi (t) + C\varphi (0) \le \psi (0)\) and \(\alpha >1\). Setting \(s = t\) and using \(\varphi (t) = 0\) leads us to the desired estimate (3.8). Owing to the Csiszár–Kullback–Pinsker inequality, it follows

Therefore, there exists an explicit \(T_0>0\) such that

and, thus,

Using this property, we estimate \(\mathcal D[u(t)]\) for \(t\ge T_0\) as follows. If \(m_i<1\), then it is similar to (3.3) that

and for \(m_i\ge 1\),

since \(\overline{u}_i(t) \ge \frac{1}{2} u_{i,\infty }\) for all \(t\ge T_0\). Therefore, for all \(t\ge T_0\), we can estimate \(\mathcal D[u(t)]\) as

Using again the entropy–entropy production estimate in the case of linear diffusion (cf. Fellner and Tang 2018, Theorem 1.1), we have

which leads to the exponential convergence

Since \(T_0\) can be explicitly computed, we in fact get the exponential convergence for all \(t\ge 0\), i.e.

for some constants \(C>0\) and \(\lambda >0\). Thanks to the Csiszár–Kullback–Pinsker inequality, we finally obtain the desired exponential convergence to equilibrium. \(\square \)

3.3 Convergence of Weak and Bounded Solutions

We consider the approximate \(u^{\varepsilon }\) to the regularised system (2.1). From Lemma 2.4, we know that, up to a subsequence,

The next lemma establishes uniform-in-\(\varepsilon \) a priori bounds.

Lemma 3.6

Under assumptions (F1)–(F3), \(u \ge 0\) with \(u_{i,0}\log u_{i,0}\in L^1(\Omega )\cap H^{-1}(\Omega )\), and \(m_i > 0\) for all \(i=1,\ldots , I\), there exists a constant \(C_T>0\) which depends on \(T>0\) but not on \(\varepsilon >0\), such that

Proof

The proof uses the idea from Laamri and Perthame (2020). For completeness, we present it in Appendix 2. \(\square \)

Proof of Theorem 1.7

By testing (2.1) with \(\psi \in L^{\infty }(Q_T)\cap L^{2(m_i+1)} (0,T;W^{1,2(m_i+1)}(\Omega ) )\), we have

From (3.11) and (3.12), we can use Vitali’s lemma to deduce

This combined with Lemma 2.1 gives

Therefore,

From the definition of \(f_i^\varepsilon \) in (2.1) as well as (3.11), we have

We show that the set \(\{f_i^\varepsilon (u^{\varepsilon })\}_{\varepsilon >0}\) is uniformly integrable. Indeed, let \(K \subset Q_T\) be a measurable and compact set with its measure being denoted by |K|. From (1.19) and (1.17), we have

It is easy to see that for any \(\delta >0\), there exists some \(C_{\delta }>0\) such that for all \(x\ge 0\),

It follows from (3.16) and Lemma 3.6 that

Therefore, for any \(\tau > 0\), we first choose \(\delta >0\) fixed such that \(C\delta < \tau /2\) and then \(K \subset Q_T\) arbitrarily according to \(C_\delta |K| < \tau /2\) to ultimately obtain

which is exactly the uniform integrability of \(\{f_i^\varepsilon (u^{\varepsilon })\}_{\varepsilon >0}\). From this and (3.15), we can apply Vitali’s lemma again to get \(f_i^\varepsilon (u^{\varepsilon }) \rightarrow f_i(u)\) strongly in \(L^1(Q_T)\) and, therefore,

Finally, using the boundedness of \(\{u^{\varepsilon }_i\}_{\varepsilon >0}\) in \(L^{m_i+1}(Q_T)\), of \(\big \{\nabla (u^{\varepsilon }_i)^\frac{m_i}{2} \big \}_{\varepsilon >0}\) in \(L^2(Q_T)\), and of \(\{f_i^{\varepsilon }(u^{\varepsilon })\}_{\varepsilon >0}\) in \(L^1(Q_T)\), it follows that

and, consequently,

From (3.14), (3.17), and (3.18), we obtain that \(u = (u_i)_{i=1,\ldots , I}\) is a global weak solution to (1.1). It is remarked that a weak solution is not automatically a renormalised solution; therefore, the weak entropy–entropy dissipation relation (1.5) and the conservation laws (1.6) are not immediate consequences. On the one hand, we can use Definition 1.6 to show directly that (1.6) also holds for weak solutions. Indeed, by choosing \(\psi = q_i\) as the test function for the equation of \(u_i\) and summing the resultants using \(\sum _{i=1}^I q_if_i(u) = 0\), we obtain the conservation laws (1.6).

On the other hand, for verifying the entropy law (1.5), we need to restrict ourselves to bounded weak solutions in order to make the subsequent arguments rigorous. We start with a sequence of functions \(\{u^{\gamma }_i\}_{\gamma >0} \subset L^{2(m_i + 1)} (0, T; W^{1, 2(m_i + 1)})\) such that \(u^{\gamma }_i \rightarrow u_i\) a.e. in \(Q_T\) and \(\nabla (u^{\gamma }_i)^{m_i/2} \rightharpoonup \nabla u_i^{m_i/2}\) in \(L^2(Q_T)\) as \(\gamma \rightarrow 0\). We may further assume that the functions \(\{u^{\gamma }_i\}_{\gamma >0}\) are uniformly bounded, and satisfying \(u^{\gamma }_i \ge u_i\) a.e. in \(Q_T\). For \(\delta >0\), we test (1.16) with \(\log (u^{\gamma }_i + \delta ) + \mu _i\) and sum over \(i = 1,\ldots , I\) to get

Thanks to \(u_i^{m_i/2}\nabla u_i^{m_i/2} \cdot \frac{\nabla (u^{\gamma }_i)}{u^{\gamma }_i+\delta } = \frac{2}{m_i} \frac{u_i^{m_i/2} (u_i^\gamma )^{1 - m_i/2}}{u_i^\gamma + \delta } \nabla u_i^{m_i/2} \cdot \nabla (u^{\gamma }_i)^{m_i/2}\) and \(\frac{u_i^{m_i/2} (u_i^\gamma )^{1 - m_i/2}}{u_i^\gamma + \delta } \le 1\), we conclude the convergence of the corresponding integral above taking the convergence properties of \(u^{\gamma }_i\) and \(\nabla u^{\gamma }_i\) into account. As the dominated convergence theorem applies to the integral on the right-hand side, we can continuously extend the linear functional \(\partial _t u_i\) to \(\log (u_i+\delta )\), which is a priori not contained in the domain of definition of \(\partial _t u_i\). Therefore,

The convergence of the first two sums is now straightforward thanks to the regularity of weak solutions. To pass to the limit for the last sum, we rewrite

where \(u + \delta = (u_i+\delta )_{i=1,\ldots , I}\) and \(C > 0\) is the constant from (1.2). The convergence of (II) is straightforward:

For (I), we use (F3) and Fatou’s lemma to obtain

Finally, from (1.11) it follows that \(|Df_i(z)|\le C \big ( 1+\sum _{j=1}^I|z_j|^{\rho -1} \big )\) and, hence,

thanks to (3.12) and the dominated convergence theorem. Inserting the estimates of (I), (II), and (III) into (3.20), then letting \(\delta \rightarrow 0\) in (3.19), we conclude that any bounded weak solution to (1.1) satisfies the weak entropy–entropy dissipation law (1.5).

Concerning bounded solutions, we first obtain from Lemma 3.6 that

Under assumptions (1.18) on the porous medium exponent, we can apply (Fellner et al. 2020, Lemma 2.2) to get

Note that the constant \(C_T\) grows at most polynomially in T. This implies the exponential convergence to equilibrium in \(L^p(\Omega )\) for any \(1\le p < \infty \) in (1.20). Indeed, from the exponential convergence in \(L^1(\Omega )\), we can use the interpolation inequality to have

for some \(0<\lambda _p < \lambda \), since \(C_T\) grows at most polynomially in T. \(\square \)

3.4 Uniform Convergence of Approximate Solutions

We remark that all the involved constants in the following results do not depend on \(\varepsilon \).

Lemma 3.7

There exists a constant \(L_0>0\) such that

Proof

It follows from

that

Thanks to (1.22), \(\sup _{\varepsilon >0}\max _{i=1,\ldots , I}\mathcal E[u^{\varepsilon }_0|u_{\infty }] < +\infty \). The \(L^1\) bound (3.21) follows with the same arguments as in Lemma 3.1. \(\square \)

Lemma 3.8

There exists a constant \(L_2>0\) such that

Proof

We introduce the shorthand notation

and \(g(v^{\varepsilon }) = f((v^{\varepsilon })^2) = f(u^{\varepsilon })\). Moreover, for each \(i\in \{1,\ldots , I\}\), we set

The desired inequality in Lemma 3.8 becomes

Using Lemma 3.7, we observe by Hölder’s inequality, noting that \(|\Omega | = 1\), that

and

Therefore, there exists a constant \(K_7>0\) such that

We decompose \(\Omega \) into \(\Omega = \Omega _1 \cup \Omega _2\) where

Note that on \(\Omega _1\) we have

Therefore, there exists a constant \(K_8>0\) such that

By using Taylor’s expansion, we have for all \(x\in \Omega _1\),

for some \(K_9>0\). From (3.25) and (3.26), we can estimate the second sum on the left-hand side of (3.23) as

On the other hand, the first sum on the left-hand side of (3.23) is estimated as

By combining (3.27) and (3.28), we have for any \(\theta \in (0,1]\),

Thus, by choosing \(\theta = \min \left\{ 1; \frac{1}{2K_8K_9R} \right\} \), we finally obtain the desired inequality (3.23) with \(L_2 = \min \left\{ \frac{\theta K_8}{2}; \frac{1}{2K_7}\right\} \). \(\square \)

Lemma 3.9

There exists a constant \(L_3>0\) such that

Proof

The proof of this lemma follows by the same arguments as in (Fellner and Tang 2018, Lemma 2.7), so we omit it here. \(\square \)

We are now ready to prove Theorem 1.9.

Proof of Theorem 1.9

We proceed similarly to the proof of Lemma 3.5. If \(\frac{(d-2)_+}{d}< m_i < 1\), we have, thanks to Lemma 3.4 and bound (3.21),

with a constant \(C = C(\Omega ,m_i) > 0\). If \(m_i\ge 1\), we apply Lemma 3.3 to get

Therefore, by using \(\Psi (x,y) = x\log (x/y) - x + y \ge \left( \sqrt{x} - \sqrt{y}\right) ^2\) and noting that

there exists a constant \(K_{10}>0\) such that

We define

and we consider the following two cases.

Case 1. If \(\widehat{F} \ge 1\) or if there exists some \(i\in \{1,\ldots , I\}\) such that \(\widehat{D}_i\ge 1\), we estimate from (3.30) that

Case 2. If \(\widehat{F}\le 1\) and \(\widehat{D}_i\le 1\) for all \(i=1,\ldots , I\), then with \(\alpha = \max _{i=1,\ldots ,I}\{1,m_i\}\), we have from (3.30) that

for some \(K_{11} = K_{11}(K_{10},\alpha )\). Thanks to Lemmas 3.8 and 3.9, we have

where

From here, we can use (Fellner and Tang 2018, Lemma 2.5) to estimate

and (Fellner and Tang 2018, Eq. (11)) to get

It follows from (3.32), (3.33), and (3.34) that

Owing to (3.31) and (3.35), there exists some \(L_4>0\) such that

The rest of this proof follows exactly from that of Theorem 1.5, which helps us to eventually obtain the exponential convergence to equilibrium

where \(C, \lambda >0\) are independent of \(\varepsilon \). \(\square \)