# On a class of triangular cross-diffusion systems and its fast reaction approximation

Elisabetta Brocchieri, Lucilla Corrias

<figure><img src="image_1.png" /></figure>

## To cite this version:

Elisabetta Brocchieri, Lucilla Corrias. On a class of triangular cross-diffusion systems and its fast reaction approximation. 2025. hal-04983668v3

HAL Id: hal-04983668

https://hal.science/hal-04983668v3

Preprint submitted on 16 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# On a class of triangular cross-diffusion systems and its fast reaction approximation

E. Brocchieri¹, L. Corrias²

4th June 2025

## Abstract

The purpose of this article is to investigate the emergence of cross-diffusion in the time evolution of two slow-fast species in competition. A class of triangular cross-diffusion system is obtained as the singular limit of a fast reaction-diffusion system. We first prove the convergence of the unique strict solution of the fast reaction-diffusion system towards a (weak, strong) solution of the cross-diffusion system, as the reaction rate $\varepsilon^{-1}$ goes to $+\infty$. Furthermore, under the assumption of small cross-diffusion, we obtain a convergence rate as well as the influence of the initial layer, due to initial data, on the convergence rate itself. Both results are obtained through energy functionals that handle the fast reaction terms uniformly in $\varepsilon$.

*Keywords.* Cross-diffusion, singular limits, dynamical systems, slow-fast manifold. 2010 Mathematics Subject Classification. Primary: 35B25, 35B40, 35K57, 35Q92, 92D25. Secondary 35B45, 35K45.

## 1 Introduction

This article deals with the emergence of cross-diffusion in the singular limit of a reaction-diffusion system with multiple time scales. The system models two species, say **u** and **v**, in competition in a bounded region of $\mathbb{R}^N$, $N \ge 1$, with reflecting boundary. Due to the inter et intra competition, the individuals of the species **u** may switch between two different states, **a** and **b**, with switching rate of order $\varepsilon^{-1}$, $\varepsilon > 0$, and change the diffusivity together with their state. Hence, the density $u^\varepsilon$ of the population **u** writes as $u^\varepsilon = u_a^\varepsilon + u_b^\varepsilon$, where $u_a^\varepsilon \ge 0$, $u_b^\varepsilon \ge 0$ are the densities of the two subpopulations (one for each state). The slow population **v** has density $v^\varepsilon$ and the global dynamic is modelled by the fast reaction system below

$$ \left\{ \begin{aligned} \partial_t u_a^\varepsilon - d_a \Delta u_a^\varepsilon &= f_a(u_a^\varepsilon, u_b^\varepsilon, v^\varepsilon) + \varepsilon^{-1} Q(u_a^\varepsilon, u_b^\varepsilon, v^\varepsilon), & \text{in } (0, \infty) \times \Omega, \\ \partial_t u_b^\varepsilon - d_b \Delta u_b^\varepsilon &= f_b(u_a^\varepsilon, u_b^\varepsilon, v^\varepsilon) - \varepsilon^{-1} Q(u_a^\varepsilon, u_b^\varepsilon, v^\varepsilon), & \text{in } (0, \infty) \times \Omega, \\ \partial_t v^\varepsilon - d_v \Delta v^\varepsilon &= f_v(u_a^\varepsilon, u_b^\varepsilon, v^\varepsilon), & \text{in } (0, \infty) \times \Omega, \end{aligned} \right. \quad (1.1) $$

where $d_a, d_b, d_v > 0$, $d_a \neq d_b$ and $\Omega$ is a smooth bounded domain of $\mathbb{R}^N$, $N \ge 1$.
System (1.1) is supplemented with homogeneous Neumann boundary conditions

$$ \nabla u_a^\varepsilon \cdot \vec{n} = \nabla u_b^\varepsilon \cdot \vec{n} = \nabla v^\varepsilon \cdot \vec{n} = 0, \quad \text{in } (0, \infty) \times \partial\Omega, \qquad (1.2) $$

where $\vec{n}$ is the unit outward normal vector, and nonnegative initial data

$$ u_a^\varepsilon(0) = u_a^{\text{in}} \ge 0, \quad u_b^\varepsilon(0) = u_b^{\text{in}} \ge 0, \quad v^\varepsilon(0) = v^{\text{in}} \ge 0, \quad \text{in } \Omega. \qquad (1.3) $$

The competitive dynamics are given by the reaction functions

$$ \begin{aligned} f_a(u_a, u_b, v) &:= \eta_a u_a (1 - au_a - cv) - \gamma_a u_a u_b, \\ f_b(u_a, u_b, v) &:= \eta_b u_b (1 - bu_b - dv) - \gamma_b u_b u_a, \\ f_v(u_a, u_b, v) &:= \eta'_v v (1 - au_a - cv) + \eta''_v v (1 - bu_b - dv), \end{aligned} \qquad (1.4) $$

with $a, b > 0$, $c, d \in \mathbb{R}_+$, $\eta_a, \eta_b > 0$, $\eta'_v, \eta''_v, \in \mathbb{R}_+$, $(c\eta'_v, d\eta''_v) \neq (0, 0)$, $\gamma_a, \gamma_b \in \mathbb{R}_+$, while the switching dynamic between the states **a** and **b** is modelled by

$$ \begin{aligned} Q(u_a, u_b, v) &:= q(u_a, u_b, v) / \Lambda(u_a, u_b, v), \\ q(u_a, u_b, v) &:= \phi(b_{\text{fast}} u_b + d_{\text{fast}} v) u_b - \psi(a_{\text{fast}} u_a + c_{\text{fast}} v) u_a, \\ \Lambda(u_a, u_b, v) &:= \phi(b_{\text{fast}} u_b + d_{\text{fast}} v) + \psi(a_{\text{fast}} u_a + c_{\text{fast}} v), \end{aligned} \qquad (1.5) $$

where $\psi, \phi$ are chosen so that $\Lambda > 0$ and $a_{\text{fast}}, b_{\text{fast}}, c_{\text{fast}}, d_{\text{fast}} \in \mathbb{R}_+$, $(a_{\text{fast}}, c_{\text{fast}}) \neq (0, 0)$, $(b_{\text{fast}}, d_{\text{fast}}) \neq (0, 0)$.

We are interested in the case where the increase in density $u_b^\varepsilon$ (respectively $u_a^\varepsilon$) pushes individuals of species **u** to migrate towards the state **a** (respectively **b**). Therefore, it is natural to consider increasing transition functions $\psi, \phi : \mathbb{R}_+ \mapsto \mathbb{R}_+$. In addition, this choice gives us relative satisfaction measures

$$ \begin{aligned} \Lambda_a(u_a, u_b, v) &:= \psi(a_{\text{fast}} u_a + c_{\text{fast}} v) / \Lambda(u_a, u_b, v), \\ \Lambda_b(u_a, u_b, v) &:= \phi(b_{\text{fast}} u_b + d_{\text{fast}} v) / \Lambda(u_a, u_b, v), \end{aligned} \qquad (1.6) $$

increasing with respect to $u_a$ and $u_b$, respectively, as well as a family of energy functionals well fitted to handle the fast reaction term $\varepsilon^{-1}Q$. Furthermore, as power laws are meaningful from the biological point of view, and general enough to allow to extend our results to any pair $(\psi, \phi)$ behaving like power functions as $x \searrow 0$ and $x \nearrow \infty$, we choose

$$ \psi(x) = (A + x)^\alpha, \quad \phi(x) = (B + x)^\beta, \quad x \ge 0, \quad \alpha, \beta > 0. \qquad (1.7) $$

Due to the symmetry of (1.1)-(1.7), i.e. to the interchangeable role of the subpopulations with densities $u_a^\varepsilon$ and $u_b^\varepsilon$, we can assume without loss of generality that $\alpha \le \beta$. As the transition rates $\alpha, \beta$ can be different from each other, we have an additional slow-fast intra-dynamic that deeply affects the mathematical analysis of the system. Hence, in order to have well defined relative satisfaction measures (1.6) and well defined and manageable energy functionals, we need to assume the following

$$ 0 < \alpha \le \beta, \quad A > 0, \quad B \ge 0 \quad \text{and} \quad B > 0 \quad \text{if} \quad \beta < 1, \qquad (H1) $$

$$ 0 \le \beta - \alpha < 2(\alpha + 3), \qquad (H2) $$

$$ a_{\text{fast}} \le a, \quad b_{\text{fast}} \le b, \quad c_{\text{fast}} \le c, \quad d_{\text{fast}} \le d. \qquad (H3) $$
Note that assumption (H2) is written in a way that highlights the difference between the transition rates, since this gap will be crucial in the further analysis. Obviously, if $\beta \le \alpha$, one has to switch the role of A and B and the role of $\alpha$ and $\beta$ in (H1), (H2). Furthermore, assumption (H3) beyond being relevant from the mathematical point of view, it is also biologically meaningful since it implies that the switching between the two states of the individuals of the species u cannot occur too fast.

The fast reaction-diffusion system system (1.1)-(1.7) is the natural generalisation of the system introduced in [3] to investigate the impact of dietary diversity of populations in competition. However, in [3] only the subpopulation with density $u_b^\varepsilon$ is in direct competition with the population v and no intra-specific competition in the population u is taken into account, i.e. $c = \gamma_a = \gamma_b = \eta_v' = 0$, $b = d$, $\eta_v'' = \eta_b$, so that $u_b^{-1}f_b = v^{-1}f_v$. Moreover, u is in direct competition with the species v, i.e. $a = b$, $c = d$, $\gamma_a = \gamma_b$, $\eta_a = \eta_b$ and $\eta_v' = \eta_v''$, so that $u_a^{-1}f_a = u_b^{-1}f_b$. Finally, in [3], $a_{\text{fast}}$ and $b_{\text{fast}}$ are zero, meaning that the fast switching between individuals of population u is uniquely determined by v.

The aim of the paper is to investigate the singular limit of (1.1)-(1.7). More specifically, we are concerned with the convergence analysis of the problem

$$ egin{cases} 
\partial_t u_b^\varepsilon - d_b \Delta u_b^\varepsilon = f_b(u_a^\varepsilon, u_b^\varepsilon, v^\varepsilon) - \varepsilon^{-1} Q(u_a^\varepsilon, u_b^\varepsilon, v^\varepsilon), & \text{in } (0, \infty) \times \Omega, \\ 
\partial_t u^\varepsilon - \Delta(d_a u_a^\varepsilon + d_b u_b^\varepsilon) = f_u(u_a^\varepsilon, u_b^\varepsilon, v^\varepsilon), & \text{in } (0, \infty) \times \Omega, \\ 
\partial_t v^\varepsilon - \Delta(d_v v^\varepsilon) = f_v(u_a^\varepsilon, u_b^\varepsilon, v^\varepsilon), & \text{in } (0, \infty) \times \Omega, \\ 
\nabla u_b^\varepsilon \cdot \vec{n} = \nabla u^\varepsilon \cdot \vec{n} = \nabla v^\varepsilon \cdot \vec{n} = 0, & \text{in } (0, \infty) \times \partial \Omega \\ 
u_b^\varepsilon = u_b^{\text{in}}, \quad u^\varepsilon(0) = u^{\text{in}} := u_a^{\text{in}} + u_b^{\text{in}}, \quad v^\varepsilon(0) = v^{\text{in}}, & \text{in } \Omega 
\end{cases} 
\quad (1.8) $$

as $\varepsilon \to 0$ and the relative rate of convergence, where we denote

$$ f_u(u_a, u_b, v) := f_a(u_a, u_b, v) + f_b(u_a, u_b, v). \quad (1.9) $$

The class of system obtained in the limit and for which we investigate the existence and uniqueness issue of global solutions, reads as the cross-diffusion system

$$ egin{cases} 
\partial_t u - \Delta(A(u, v)) = f_u(u_a^*(u, v), u_b^*(u, v), v), & \text{in } (0, \infty) \times \Omega, \\ 
\partial_t v - \Delta(d_v v) = f_v(u_a^*(u, v), u_b^*(u, v), v), & \text{in } (0, \infty) \times \Omega, \\ 
\nabla A(u, v) \cdot \vec{n} = \nabla v \cdot \vec{n} = 0, & \text{in } (0, \infty) \times \partial \Omega, \\ 
u(0) = u^{\text{in}}, \quad v(0) = v^{\text{in}}, & \text{in } \Omega 
\end{cases} 
\quad (1.10) $$

where

$$ A(u, v) := d_a u_a^*(u, v) + d_b u_b^*(u, v) \quad (1.11) $$

and $(u_a^*, u_b^*)$ is the $C^1$ maps from $\mathbb{R}_+^2$ to $\mathbb{R}_+^2$ such that, for all $(\tilde{u}, \tilde{v}) \in \mathbb{R}_+^2$, the pair $(u_a^*(\tilde{u}, \tilde{v}), u_b^*(\tilde{u}, \tilde{v}))$ is the unique nonnegative solution of the nonlinear system

$$ egin{cases} 
u_a^* + u_b^* = \tilde{u} \\ 
Q(u_a^*, u_b^*, \tilde{v}) = 0, 
\end{cases} 
\quad (1.12) $$

with $Q$ defined in (1.5).

Note that it is not necessary to consider initial data $u_a^{\text{in}}, u_b^{\text{in}}, v^{\text{in}}$ belonging to the critical manifold given by (1.12), i.e. initial data satisfying $Q(u_a^{\text{in}}, u_b^{\text{in}}, v^{\text{in}}) = 0$. Moreover,

$$ u_a^*(u, v) = \Lambda_b^*(u, v)u, \quad u_b^*(u, v) = \Lambda_a^*(u, v)u, $$
where $\Lambda_a^*, \Lambda_b^*$ are the satisfaction measures (1.6) evaluated at $(u_a^*(u, v), u_b^*(u, v), v)$. Therefore, the reaction terms in (1.10) still write as Lotka-Volterra competitive reactions, with coefficients depending on $(u, v)$.

The emergence of cross-diffusion as a fast reaction singular limit has been observed in several mathematical models for ecology, biology or chemistry, and more specifically in the context of competitive interactions [14], predator-prey interactions [8, 17, 26], dietary diversity and starvation [3] and enzyme reaction [27]. On the other hand, cross-diffusion systems have attracted significant interest, at least since the seminal paper of Shigesada-Kawasaki-Teramoto [25], because cross-diffusion can induce instability and thus explain pattern formation whereas linear diffusion cannot (see also [2, 15, 18, 23] and the references therein). The mathematical analysis of cross-diffusion systems is delicate and has received a lot of attention in the last two decades. A fundamental contribution to the existence of solutions issue was given by Chen and Jüngel [6, 7], for cross-diffusion systems with entropy structure, including the SKT model in [25]. These entropy methods were shown to be robust enough to treat generalisations of the SKT system [12, 13]. Afterwards, the relation between the structure of systems involving cross-diffusion and the existence of an entropy functional has been deeply investigated (see e.g. [9, 20]).

In this article, the authors aim to give a contribution to the aforementioned tusk, showing that there are biologically relevant cross-diffusion systems, that, even though they do not have an entropy structure, they have a family of underlying energies helpful to obtain a positive answer to the existence issue. Moreover, we also address the issue of the convergence rate, as $\epsilon$ goes to 0, of $(u^\epsilon, v^\epsilon)$ towards $(u, v)$ and of $(u_a^\epsilon, u_b^\epsilon)$ towards $(u_a^*(u, v), u_b^*(u, v))$. This problem is strictly linked with the crucial topic of invariant slow manifolds of slow-fast dynamical systems in finite dimension. The underlying theory has been extended to the case of slow-fast PDE systems, e.g. fast reaction-diffusion system, (see [11] and the reference therein). Here, employing purely analytical tools, different from those used in [11], we obtain the convergence rate and quantify its dependence on the initial layer due to initial data not lying on the critical manifold given by (1.12).

Triangular cross-diffusion systems, similar to (1.10), has been considered also in [4], where the authors investigate the existence and weak-strong stability issues, using regularisation techniques and fixed-point arguments. However, system (1.10)–(1.12) don’t fit in the class analysed in [4].

The article is organised as follows. In Section 2 we state the main results and a set of notations. In Section 3 we introduce the family of energy functionals for (1.1)–(1.7). The time evolution of the energies along the trajectories of the unique strict solution $(u_a^\epsilon, u_b^\epsilon, v^\epsilon)$ of (1.1)–(1.7) is analysed in Sections 4–6. This analysis will provide us a priori estimates on $(u_a^\epsilon, u_b^\epsilon, v^\epsilon)$, uniform in $\epsilon$, stated in Lemma 2.2. The proof of the existence result for the cross-diffusion system (1.10)–(1.12), stated in Theorem 2.5, is given in Section 7, while Section 8 is devoted to the proof of the uniqueness of bounded solutions. We conclude the article with Section 9 where we obtain a convergence rate for the singular limit $\epsilon \to 0$, stated in Theorem 2.7. For the sake of completeness, in Appendix A we give the existence result of the unique strict solution $(u_a^\epsilon, u_b^\epsilon, v^\epsilon)$ of the fast reaction-diffusion system, stated in Theorem 2.1. Finally, Appendix B is devoted to the proof of the solvability of the nonlinear system (1.12).

**Acknowledgment** The authors warmly thank Laurent Desvillettes and Helge Dietert for the fruitful discussions they had all along the preparation of the article.
## 2 Statements and main results

The starting results are the existence of a unique nonnegative strict solution of system (1.1)–(1.7) together with basic estimates (independent of $\varepsilon$) on the solution. All of this is stated in Theorem 2.1 and obtained applying the classical theory of analytic semigroups, the maximal (parabolic) regularity and taking advantage of the competition dynamics (1.4). We refer mainly to [22] and we sketch the proof in Appendix A, for the reader convenience.

Let us define

$$D_p := \{w \in W^{2,p}(\Omega) : \nabla w \cdot \vec{n} = 0 \text{ on } \partial\Omega\}, \quad p \in (1, +\infty), \quad (2.1)$$

$$\bar{\eta} := \eta_a \lor \eta_b, \quad \eta := a\eta_a \land b\eta_b, \quad \eta_v := \eta'_v + \eta''_v, \quad r_v := c\eta'_v + d\eta''_v. \quad (2.2)$$

**Theorem 2.1** (Well-posedness of the fast reaction-diffusion system). Let $\Omega \subset \mathbb{R}^N, N \ge 1$, be a bounded open set with $C^2$ boundary $\partial\Omega$ and let $\varepsilon > 0$. Assume (1.4)–(1.7), (H1) and let

$$u_a^{in}, u_b^{in}, v^{in} \in \bigcap_{1<p<\infty} D_p \quad (2.3)$$

be non-negative initial data, with $u_a^{in}, u_b^{in}$ not identically zero. It follows that there exists a triplet

$$(u_a^{\varepsilon}, u_b^{\varepsilon}, v^{\varepsilon}) \in C^1([0, \infty); (L^p(\Omega))^3) \cap C^0([0, \infty); D_p^3), \quad \forall p \in (1, +\infty) \quad (2.4)$$

with $u_a^{\varepsilon}, u_b^{\varepsilon} > 0$ and $v^{\varepsilon} \ge 0$ on $(0, \infty) \times \Omega$, which is the unique strict solution of (1.1) with boundary conditions (1.2) and initial conditions (1.3). In addition, the solution satisfies the following estimates (independent of $\varepsilon > 0$)

$$\|u_a^{\varepsilon} + u_b^{\varepsilon}\|_{L^{\infty}((0,\infty);L^1(\Omega))} \le \max\{\|u_a^{in} + u_b^{in}\|_{L^1(\Omega)}, 2|\Omega|\bar{\eta}\eta^{-1}\} =: K_1, \quad (2.5)$$

$$\|v^{\varepsilon}\|_{L^{\infty}((0,\infty)\times\Omega)} \le \max \left\{ \|v^{in}\|_{L^{\infty}(\Omega)}, \frac{\eta_v}{r_v} \right\} =: K_{\infty}, \quad (2.6)$$

and, for all $T > 0$,

$$\|u_a^{\varepsilon}\|_{L^2((0,T)\times\Omega)}^2 + \|u_b^{\varepsilon}\|_{L^2((0,T)\times\Omega)}^2 \le \eta^{-1}\|u_a^{in} + u_b^{in}\|_{L^1(\Omega)} + \bar{\eta}\eta^{-1}K_1 T =: K_2, \quad (2.7)$$

$$\|\partial_t v^{\varepsilon}\|_{L^2((0,T)\times\Omega)} + \sum_{i,j} \|\partial_{x_i} x_j v^{\varepsilon}\|_{L^2((0,T)\times\Omega)} \le C_1(K_2, K_{\infty}, T, |\Omega|), \quad (2.8)$$

$$\|\nabla v^{\varepsilon}\|_{L^4((0,T)\times\Omega)} \le C_2(K_2, K_{\infty}, N, T, |\Omega|). \quad (2.9)$$

As mentioned in the introduction, the main results of the article are the existence and uniqueness of a global (weak, strong) solution $(u, v)$ of the cross-diffusion system (1.10)–(1.12), obtained when $\varepsilon \to 0$ in (1.8), and the rate of convergence.

The existence result for the cross-diffusion system, stated in Theorems 2.5 and proved in Section 7, needs the proof of further (uniform in $\varepsilon$) estimates. For this purpose, we construct in Sections 3 a well fitted family of energy functionals (3.1)–(3.3). The analysis of the time evolution of the energies is developed in Sections 4–6 and leads to Lemma 2.2 below. For the proof of this Lemma, we need to define the exponents increasing in $\alpha, \beta$

$$\begin{aligned} q(p) &:= p + \alpha(p - 1) = (\alpha + 1)(p - 1) + 1, \quad (2.10) \\ r(p) &:= p + \beta(p - 1) = (\beta + 1)(p - 1) + 1, \end{aligned}$$
and the critical values

$$p_{\alpha} := 1 + \frac{1}{1 + \alpha} \in (1, 2), \quad p_{\beta} := 1 + \frac{1}{1 + \beta} \in (1, 2). \qquad (2.11)$$

Note that $p_{\beta} \le p_{\alpha}$ and $q(p) \le r(p)$, since $0 < \alpha \le \beta$, and that

$$q(p_{\alpha}) = r(p_{\beta}) = 2. \qquad (2.12)$$

Furthermore, writing $r(p) = q(p) + (\beta - \alpha)(p - 1)$, one see that the gap between $r(p)$ and $q(p)$ is controlled by the gap $\beta - \alpha$ between the transition rates. Therefore, the size of $\beta - \alpha$ will be crucial in the bootstrap procedure performed in the proof of Lemma 2.2. To carry out the bootstrap, we define the decreasing family of intervals

$$I_n := (2(\alpha + 1), 2(\alpha + 1) + \frac{4}{(\alpha + 1)^n}), \quad n \in \mathbb{N} \cup \{0\}, \qquad (2.13)$$

so that the admissible set $[0, 2(\alpha+3))$ for $\beta - \alpha$ in (H2) reads as $[0, 2(\alpha+1)] \cup (\cup_n I_n)$. Moreover, if $\beta - \alpha > 2(\alpha + 1)$, we also denote

$$n_{\alpha, \beta} \text{ the largest integer such that } \beta - \alpha \in I_{n_{\alpha, \beta}}. \qquad (2.14)$$

**Lemma 2.2 (Energy estimates).** Under hypothesis of Theorem 2.1 and assuming in addition (H2), (H3), for all $T > 0$, it holds

(i) there exists $C(T) > 0$ such that, for all $p \in [p_{\beta}, p_{\alpha}]$ and $\varepsilon > 0$,

$$\begin{aligned} & \|u_a^{\varepsilon}\|_{L^{\infty}(0,T;L^{q(p)}(\Omega))} + \|u_b^{\varepsilon}\|_{L^{\infty}(0,T;L^{r(p)}(\Omega))} \\ & + \|\nabla(u_a^{\varepsilon})^{q(p)/2}\|_{L^2((0,T)\times\Omega)} + \|\nabla(u_b^{\varepsilon})^{r(p)/2}\|_{L^2((0,T)\times\Omega)} \\ & + \|u_a^{\varepsilon}\|_{L^{q(p)+1}((0,T)\times\Omega)} + \|u_b^{\varepsilon}\|_{L^{r(p)+1}((0,T)\times\Omega)} \le C(T); \end{aligned} \qquad (2.15)$$

(ii) if $\beta - \alpha \in [0, 2(\alpha + 1)]$, for all $p \in [2, +\infty)$, there exists $C(p, T) > 0$ such that, for all $\varepsilon > 0$,

$$\begin{aligned} & \|u_a^{\varepsilon}\|_{L^{\infty}(0,T;L^{q(p)}(\Omega))} + \|u_b^{\varepsilon}\|_{L^{\infty}(0,T;L^{r(p)}(\Omega))} \\ & + \|\nabla(u_a^{\varepsilon})^{q(p)/2}\|_{L^2((0,T)\times\Omega)} + \|\nabla(u_b^{\varepsilon})^{r(p)/2}\|_{L^2((0,T)\times\Omega)} \le C(p,T); \end{aligned} \qquad (2.16)$$

(iii) if $\beta - \alpha \in (2(\alpha + 1), 2(\alpha + 3)) = \cup_n I_n$, there exists $C(n_{\alpha, \beta}, T) > 0$ such that, for all $p \in [2, 1 + (\alpha + 1)^{n_{\alpha, \beta}}]$ and $\varepsilon > 0$,

$$\begin{aligned} & \|u_a^{\varepsilon}\|_{L^{\infty}(0,T;L^{q(p)}(\Omega))} + \|u_b^{\varepsilon}\|_{L^{\infty}(0,T;L^{r(p)}(\Omega))} \\ & + \|\nabla(u_a^{\varepsilon})^{q(p)/2}\|_{L^2((0,T)\times\Omega)}^2 + \|\nabla(u_b^{\varepsilon})^{r(p)/2}\|_{L^2((0,T)\times\Omega)}^2 \\ & + \|u_a^{\varepsilon}\|_{L^{q(p)+1}((0,T)\times\Omega)}^{q(p)+1} + \|u_b^{\varepsilon}\|_{L^{r(p)+1}((0,T)\times\Omega)}^{r(p)+1} \le C(n_{\alpha, \beta}, T). \end{aligned} \qquad (2.17)$$

(iv) there exists $C(T) > 0$ such that, for all $\varepsilon > 0$,

$$\varepsilon^{-\frac{1}{2}} \|\Lambda^{1/2} Q\|_{L^2((0,T)\times\Omega)} \le C(T)(1 + \|u_a^{\varepsilon}\|_{L^2((0,T)\times\Omega)} + \|u_b^{\varepsilon}\|_{L^2((0,T)\times\Omega)}). \qquad (2.18)$$
Finally, for all $T > 0$ and $p \in (1, +\infty)$ if $\beta - \alpha \in [0, 2(\alpha+1)]$, or $p \in (1, 2+(\alpha+1)^{n_{\alpha,\beta}+1})$ if $\beta - \alpha \in (2(\alpha+1), 2(\alpha+3))$, there exist $C_1(p, T) > 0$ and $C_2(p, T) > 0$ such that, for all $\varepsilon > 0$, it holds

$$ \|\partial_t v^\varepsilon\|_{L^p((0,T)\times\Omega)} + \sum_{i,j} \|\partial_{x_i x_j} v^\varepsilon\|_{L^p((0,T)\times\Omega)} \le C_1 \quad (2.19) $$

$$ \|\nabla v^\varepsilon\|_{L^{2p}((0,T)\times\Omega)} \le C_2. \quad (2.20) $$

**Remark 1.** Thanks to identity (2.12), we obtain the $L^2((0,T) \times \Omega)$ uniform estimate of $\nabla u_a^\varepsilon$ and $\nabla u_b^\varepsilon$ from (2.15), taking $p = p_\alpha$ and $p = p_\beta$, respectively.

Next, the cross-diffusion term in (1.10) is due to the convergence of the pair $(u_a^\varepsilon, u_b^\varepsilon)$ towards the unique solution of the nonlinear system (1.12). Estimate (2.23) in Lemma 2.4 (proved in Section 7) is the key tool to obtain this convergence. The solvability of (1.12) and the regularity of the map $(u_a^*, u_b^*) : \mathbb{R}_+^2 \to \mathbb{R}_+^2$, in turn, are straightforward consequences of the regularity of $Q$ and of the implicit function theorem. These properties are resumed in Lemma 2.3 below and the proof is given for completeness in Appendix B.

**Lemma 2.3** (Existence and regularity of the map $(u_a^*, u_b^*)$). Assume (1.5), (1.7), (H1). For all $(\tilde{u}, \tilde{v}) \in \mathbb{R}_+^2$, there exists a unique nonnegative solution $(u_a^*(\tilde{u}, \tilde{v}), u_b^*(\tilde{u}, \tilde{v}))$ of (1.12). Moreover, $(u_a^*, u_b^*) \in (C^1(\mathbb{R}_+^2))^2$, with

$$ \partial_{\tilde{u}} u_a^*(\tilde{u}, \tilde{v}), \partial_{\tilde{u}} u_b^*(\tilde{u}, \tilde{v}) \in (0, 1) \quad (2.21) $$

and, assuming $a_{fast}$ and $b_{fast}$ both non-zero,

$$ -\frac{c_{fast}}{a_{fast}} \le \partial_{\tilde{v}} u_a^*(\tilde{u}, \tilde{v}) = -\partial_{\tilde{v}} u_b^*(\tilde{u}, \tilde{v}) \le \frac{d_{fast}}{b_{fast}}. \quad (2.22) $$

**Lemma 2.4** (Convergence of $(u_a^\varepsilon, u_b^\varepsilon)$). Under assumptions of Lemma 2.2, let $(u_a^\varepsilon, u_b^\varepsilon, v^\varepsilon)$ be the unique strict solution of (1.1)-(1.7), given by Theorem 2.1, and let $u^\varepsilon = u_a^\varepsilon + u_b^\varepsilon$. For all $T > 0$, there exists $C(T) > 0$ such that, for all $\varepsilon > 0$, it holds

$$ \|u_b^\varepsilon - u_b^*(u^\varepsilon, v^\varepsilon)\|_{L^2((0,T)\times\Omega)} \le \|Q(u_a^\varepsilon, u_b^\varepsilon, v^\varepsilon)\|_{L^2((0,T)\times\Omega)} \le A^{-\frac{\alpha}{2}} C(T) \sqrt{\varepsilon}. \quad (2.23) $$

As $|u_a^\varepsilon - u_a^*(u^\varepsilon, v^\varepsilon)| = |u_b^\varepsilon - u_b^*(u^\varepsilon, v^\varepsilon)|$, the same inequality holds for $u_a^\varepsilon - u_a^*(u^\varepsilon, v^\varepsilon)$.

Having these results at hand, we can state the existence of a global (weak, strong) solution $(u, v)$ of the cross-diffusion system. When $\beta - \alpha \in (2(\alpha+1), 2(\alpha+3))$, we will: (i) use the best integrability for $u_a^\varepsilon$ given in (2.17), namely the $L^{q(p)+1}$ integrability with $p = 1 + (\alpha + 1)^{n_{\alpha,\beta}}$, that gives $q(p) + 1 = 2 + (\alpha + 1)^{n_{\alpha,\beta}+1}$; (ii) add a dimension depending condition on $\beta - \alpha$.

**Theorem 2.5** (Existence for the cross-diffusion system). Let $\Omega \subset \mathbb{R}^N$, $N \ge 1$, be a bounded open set with $C^2$ boundary $\partial\Omega$ and assume (H1), (H2), (H3). Let $u_a^{in}$, $u_b^{in}$, $v^{in}$ be non-negative initial data satisfying (2.3), with $u_a^{in}$, $u_b^{in}$ not identically zero. Furthermore, if $N > 6$ and $\beta - \alpha \in (2(\alpha+1), 2(\alpha+3))$, assume that $(\alpha+1)^{n_{\alpha,\beta}+1} \ge 2$. There exists a pair of nonnegative measurable functions $(u, v) : (0, \infty) \times \Omega \to \mathbb{R}_+^2$ such that, for all $p$ satisfying

$$ \left\{ \begin{array}{ll} p \in (2, +\infty), & \text{if } \beta - \alpha \in [0, 2(\alpha + 1)], \\ p = 2 + (\alpha + 1)^{n_{\alpha, \beta} + 1}, & \text{if } \beta - \alpha \in (2(\alpha + 1), 2(\alpha + 3)), \end{array} \right. \quad (2.24) $$

(see (2.13), (2.14)), $s = \frac{p}{2} \wedge 2$ and all $T > 0$, it holds
(i) $u \in L^2(0, T; H^1(\Omega)) \cap L^p((0, T) \times \Omega) \cap L^\infty(0, T; L^{p-1}(\Omega))$ and
$\partial_t u \in L^s(0, T; H^{-1}(\Omega))$,

(ii) $v \in W^{1,p}((0, T) \times \Omega) \cap L^p(0, T; W^{2,p}(\Omega)) \cap L^\infty((0, \infty) \times \Omega)$,

(iii) up to the extraction of a subsequence from the strict solutions $(u_a^\varepsilon, u_b^\varepsilon, v^\varepsilon)_\varepsilon$ \\
$$ (u_a^\varepsilon, u_b^\varepsilon) \to (u_a^*(u, v), u_b^*(u, v)), \quad \text{a.e. in } (0, \infty) \times \Omega, \text{ as } \varepsilon \to 0, $$ \\
where $(u_a^*(u, v), u_b^*(u, v))$ is the unique solution of (1.12),

(iv) $u_a^*(u, v), u_b^*(u, v) \in L^2(0, T; H^1(\Omega)) \cap L^p((0, T) \times \Omega) \cap L^\infty(0, T; L^{p-1}(\Omega)),$

(v) $(u, v)$ is a global (weak, strong) solution of (1.10), (1.11), (1.12), i.e. for all $T > 0$ and all $w \in C^1([0, T]; H^1(\Omega))$ such that $w(T) = 0$, it holds

$$
\left\{
\begin{array}{l}
- \int_0^T \int_\Omega u \partial_t w \, dxdt + \int_0^T \int_\Omega \nabla(d_a u_a^*(u, v) + d_b u_b^*(u, v)) \cdot \nabla w \, dxdt \\
\qquad = \int_\Omega u^{in} w(0) \, dx + \int_0^T \int_\Omega f_u(u_a^*(u, v), u_b^*(u, v), v) w \, dxdt \\
\partial_t v = d_v \Delta v + f_v(u_a^*(u, v), u_b^*(u, v), v) \quad \text{in } L^p((0, T) \times \Omega) \quad (2.25) \\
\sum_i \gamma_0(\partial_{x_i} v) n_i = 0, \quad \text{a.e. in } (0, T) \times \partial \Omega \\
u(0) = u^{in}, \quad v(0) = v^{in} \quad \text{a.e. in } \Omega
\end{array}
\right.
$$

One can show that the global (weak, strong) solution $(u, v)$ obtained in Theorem 2.5 enjoys additional regularity properties. However, exploring the regularity of $(u, v)$ lies beyond the scope of the paper. Instead, we conclude the analysis by proving the uniqueness of solutions of (1.10)–(1.12) that are bounded in both components. Note that the boundedness of the component $u$ is required only to handle the reaction. Theorem 2.6 is proved in Section 8.

**Theorem 2.6** (Uniqueness and stability). Under the assumptions of Theorem 2.5, let $(u_i, v_i)$, $i = 1, 2$, be two solutions with initial data $(u_i^{in}, v_i^{in})$. Assume in addition that, for all $T > 0$, $u_i \in L^\infty((0, T) \times \Omega)$, $i = 1, 2$, and $a_{fast}$, $b_{fast}$ are both non-zero. It follows that, for all $T > 0$, there exists $C(T, \|u_i\|_{L^\infty((0,T)\times\Omega)}, \|v_i\|_{L^\infty((0,\infty)\times\Omega)}) > 0$ such that

$$
\begin{align}
& \|u_1 - u_2\|_{L^2((0,T) \times \Omega)}^2 + \|v_1 - v_2\|_{L^2((0,T) \times \Omega)}^2 \nonumber \\
& \le C \left( \|u_1^{in} - u_2^{in}\|_{L^2(\Omega)}^2 + \|v_1^{in} - v_2^{in}\|_{L^2(\Omega)}^2 \right). \tag{2.26}
\end{align}
$$

Concerning the rate of convergence issue, we analyse the time evolution of the $L^2(\Omega)$ norms of $u^\varepsilon - u$, $v^\varepsilon - v$, $u_b^\varepsilon - u_b^*(u, v)$ via the functional (9.6) which include the ad hoc sub-functional (9.7)–(9.8) designed to handle the fast reaction. The result is established under additional regularity assumptions of the solutions and small cross-diffusion, i.e. small $|d_b - d_a|$. A careful reading of the proof makes it clear that the diffusivity coefficient $d_v$ can only help to optimise the smallness condition (2.30), but cannot remove it. The obtained estimates also illustrate how the initial layer

$$ \varepsilon_{\text{in}} := \|u_b^{\text{in}} - u_b^*(u^{\text{in}}, v^{\text{in}})\|_{H^1(\Omega)} \quad (2.27) $$

slows down the convergence rate.
**Theorem 2.7 (Rate of convergence).** Under the hypothesis of Theorem 2.5, let $(u_a^\epsilon, u_b^\epsilon, v^\epsilon)$ be the unique nonnegative global strict solution of (1.1)-(1.7). Assume in addition $d_b > d_a$, $A, B > 0$ for all $\alpha, \beta$ and that $u^\epsilon = u_a^\epsilon + u_b^\epsilon$ is uniformly bounded locally in time, i.e. for all $T > 0$ there exists $M_T > 0$ such that, for all $\epsilon \in (0, 1)$, it holds

$$ \|u^\epsilon\|_{L^\infty((0,T)\times\Omega)} \le M_T. \qquad (2.28) $$

Let $(u, v)$ be a nonnegative global classical solution of (1.10)-(1.12) such that, for all $T > 0$,

$$ u, v \in C^0([0, T]; C^3(\bar{\Omega})) \cap C^1([0, T]; C^1(\bar{\Omega})). \qquad (2.29) $$

Then, for all $T > 0$, there exists a constant $C_1(\alpha, \beta, A, B, T) > 0$ such that, if

$$ 1 < \frac{d_b}{d_a} < 1 + C_1(\alpha, \beta, A, B, T), \qquad (2.30) $$

there exists $C_2(T) > 0$ such that, for all $\epsilon \in (0, 1)$ and for $\epsilon_{in}$ defined in (2.27), it holds

$$ \|u^\epsilon - u\|_X + \|v^\epsilon - v\|_X \le C_2(T)(\epsilon + \epsilon^{\frac{1}{2}}\epsilon_{in}), \qquad (2.31) $$

$$ \|u_a^\epsilon - u_a^*(u, v)\|_Y + \|u_b^\epsilon - u_b^*(u, v)\|_Y \le C_2(T) (\epsilon + \epsilon^{\frac{1}{2}} \epsilon_{in}), \qquad (2.32) $$

where $\|\cdot\|_X := \|\cdot\|_{L^2(0,T;H^1(\Omega))} + \|\cdot\|_{L^\infty(0,T;L^2(\Omega))}$ and $\|\cdot\|_Y := \|\cdot\|_{L^2(0,T;H^1(\Omega))}$.

A by-product of the above convergence result is again the uniqueness of smooth solutions, without the requirements $a_{fast} \ne 0$ and $b_{fast} \ne 0$.

**Notations.** Hereafter, $\nabla$ and $\Delta$ will denote the gradient with respect to the spatial variable $x$ and the Laplacian, while $D$ and $D^2$ will denote the gradient and the Hessian with respect to non spatial variables. Moreover, $\partial_i$ and $\partial_{ij}$ will denote respectively the partial derivative with respect to the $i$-th variable and the partial derivative with respect to the $j$-th and $j$-th variables, whatever they are. For the sake of simplicity, in the computations of Sections 3–6 we will keep explicit only the constants depending on $p \ge 2$, and we will write $x \lesssim y$ meaning that there exists a universal constant $C > 0$, not depending on $p$, such that $x \le Cy$. Finally, when there is no risk of confusion, we will denote $\Omega_T$, $T > 0$, the usual open cylinder $(0, T) \times \Omega$ and we will omit the $\epsilon$ superscript.

## 3 A family of energy functionals $\mathcal{E}_p(u_a, u_b, v)$

System (1.1)-(1.7) is naturally endowed with the following family of energy functionals

$$ \mathcal{E}_p(u_a, u_b, v) := \int_{\Omega} h_p(u_a, u_b, v) dx, \quad u_a, u_b, v \in \mathbb{R}_+, \quad p \ge 1, \qquad (3.1) $$

with the total energy density $h_p$ given by

$$ h_p(u_a, u_b, v) := h_{a,p}(u_a, v) + h_{b,p}(u_b, v) \qquad (3.2) $$

and the partial energy densities $h_{a,p}$, $h_{b,p}$ defined as

$$ h_{a,p}(u_a, v) := \int_0^{u_a} \psi^{p-1}(a_{\text{fast}}z + c_{\text{fast}}v)z^{p-1}dz, \qquad (3.3) $$

$$ h_{b,p}(u_b, v) := \int_0^{u_b} \phi^{p-1}(b_{\text{fast}}z + d_{\text{fast}}v)z^{p-1}dz. $$
Moreover, setting

$$ \theta(z, v) := A + a_{\text{fast}}z + c_{\text{fast}}v, \quad \omega(z, v) := B + b_{\text{fast}}z + d_{\text{fast}}v, \quad (3.4) $$

by (1.7), the partial energy densities in (3.3) rewrite as

$$ \begin{aligned} h_{a,p}(u_a, v) &= \int_0^{u_a} \theta^{\alpha(p-1)}(z, v) z^{p-1} dz, \\ h_{b,p}(u_b, v) &= \int_0^{u_b} \omega^{\beta(p-1)}(z, v) z^{p-1} dz. \end{aligned} \quad (3.5) $$

**Remark 2.** It is worth noticing that in definition (3.3) the transition functions $\psi, \phi$ are not renormalised by $\Lambda(u_a, u_b, v)$ unlike in $Q(u_a, u_b, v)$ (see (1.5)). This choice will be clear later.

The interest in the family of energies (3.1)–(3.5) is threefold. It allows us to obtain further a priori estimates on the densities $u_a^\varepsilon, u_b^\varepsilon$ and their gradients in Lebesgue spaces, to handle easily the contribution due to the fast reaction $\varepsilon^{-1}Q$ in the aforementioned estimates, and to obtain the convergence of $\Lambda^{\frac{1}{2}}Q$ towards 0, as $\varepsilon \to 0$. Indeed, on the one hand, using for all $z, v \ge 0$,

$$ \theta(z, v) \ge a_{\text{fast}} z \quad \text{and} \quad \omega(z, v) \ge b_{\text{fast}} z, \quad (3.6) $$

it is easily seen from (3.2), (3.5), that, for $p \ge 1$,

$$ a_{\text{fast}}^{\alpha(p-1)} \frac{u_a^{q(p)}}{q(p)} + b_{\text{fast}}^{\beta(p-1)} \frac{u_b^{r(p)}}{r(p)} \le h_p(u_a, u_b, v) \le a_{\text{fast}}^{-p} \frac{\theta(u_a, v)^{q(p)}}{q(p)} + b_{\text{fast}}^{-p} \frac{\omega(u_b, v)^{r(p)}}{r(p)} $$

implying that $\mathcal{E}_p$ is well defined, for all $p \ge 1$, along the trajectories of the solution $(u_a^\varepsilon, u_b^\varepsilon, v^\varepsilon)$ obtained in Theorem 2.1, and also that, for all $T > 0$,

$$ \|u_a^\varepsilon\|_{L^\infty(0,T;L^{q(p)}(\Omega))}^q + \|u_b^\varepsilon\|_{L^\infty(0,T;L^{r(p)}(\Omega))}^r \le \frac{r(p)}{a_{\text{fast}}^{\alpha(p-1)} \wedge b_{\text{fast}}^{\beta(p-1)}} \mathcal{E}_p(T). \quad (3.7) $$

On the other hand, denoting $\mathcal{F} := (f_a, f_b, f_v)^T$, the evolution of $\mathcal{E}_p$ along the solution is described by the differential equation

$$ \begin{aligned} \mathcal{E}'_p(t) &= \frac{d}{dt} \int_{\Omega} h_p(u_a^\varepsilon, u_b^\varepsilon, v^\varepsilon) dx = \int_{\Omega} (\partial_1 h_p \partial_t u_a^\varepsilon + \partial_2 h_p \partial_t u_b^\varepsilon + \partial_3 h_p \partial_t v^\varepsilon) dx \\ &= \int_{\Omega} (d_a \partial_1 h_p \Delta u_a^\varepsilon + d_b \partial_2 h_p \Delta u_b^\varepsilon + d_v \partial_3 h_p \Delta v^\varepsilon) dx \\ &\quad + \int_{\Omega} Dh_p \cdot \mathcal{F} dx + \frac{1}{\varepsilon} \int_{\Omega} (\partial_1 h_p - \partial_2 h_p) Q dx \\ &=: I_{\text{diff}}^p + I_{\text{rea}}^p + I_{\text{fast}}^p. \end{aligned} \quad (3.8) $$

Then, we see that, for all $p \ge 1$, it holds

$$ \begin{aligned} I_{\text{fast}}^p &:= \frac{1}{\varepsilon} \int_{\Omega} (\partial_1 h_p - \partial_2 h_p) Q dx \\ &= -\frac{1}{\varepsilon} \int_{\Omega} \left[ (\phi(b_{\text{fast}}u_b^\varepsilon + d_{\text{fast}}v^\varepsilon)u_b^\varepsilon)^{p-1} - (\psi(a_{\text{fast}}u_a^\varepsilon + c_{\text{fast}}v^\varepsilon)u_a^\varepsilon)^{p-1} \right] Q dx. \end{aligned} $$
As $x \mapsto x^{p-1}$ is an increasing function on $\mathbb{R}_+$, the latter and (1.5) gives

$$I_{\text{fast}}^{p} \leq 0, \quad \forall p \geq 1, \qquad (3.9)$$

so that $I_{\text{fast}}^p$ can be neglected in the evolution equation (3.8), whenever it is useless, i.e. $p \neq 2$. When $p=2$, $I_{\text{fast}}^2$ reads as

$$I_{\text{fast}}^2 = -\frac{1}{\varepsilon} \int_{\Omega} q(u_a^{\varepsilon}, u_b^{\varepsilon}, v^{\varepsilon}) Q(u_a^{\varepsilon}, u_b^{\varepsilon}, v^{\varepsilon}) dx = -\frac{1}{\varepsilon} \int_{\Omega} \Lambda Q^2(u_a^{\varepsilon}, u_b^{\varepsilon}, v^{\varepsilon}) dx, \quad (3.10)$$

thus allowing to obtain the convergence of $\|\Lambda^{1/2} Q\|_{L^2(\Omega_T)}$ towards zero, as $\varepsilon \to 0$, with rate $\varepsilon^{1/2}$, under assumption (H2), (see (2.18)). This convergence result will be crucial in proving the convergence of the solution of the fast reaction system towards the solution of cross-diffusion system.

Finally, when $p=1$ the total energy density (3.2) reduce to $h_1(u_a, u_b, v) = u_a + u_b$, so that $I_{\text{diff}}^1 = I_{\text{fast}}^1 = 0$. It follows the uniform control of the densities $u_a^{\varepsilon}$, $u_b^{\varepsilon}$ in the Lebesgue spaces $L^{\infty}(0, \infty; L^1(\Omega))$ and $L^2(\Omega_T)$ obtained in Theorem 2.1. In order to estimate $I_{\text{rea}}^p$ and $I_{\text{diff}}^p$ in (3.8) with $p > 1$ and bootstrap the above $L^1$ and $L^2$ estimates to $L^p$ estimates, $p > 2$, (see Lemma 2.2), we need instead a suitable analysis of the hessian matrix Hess($h_p$). Particular attention must be paid to the critical case $p = p_\beta$ (see(2.11)), which requires to assume $A > 0$ to control the Hess($h_{a,p}$), (see (3.16)). This analysis is done in the rest of this section.

**Remark 3.** The energy functionals (3.1) are reminiscent of the functionals introduced in [14, 3]. It is worth noticing that they are not the sum of functionals of the single densities $u_a, u_b$ and $v$.

## 3.1 The gradient of the total energy density $h_p$

Let $p > 1$. From (3.2),(3.4),(3.5), the gradient $Dh_p$ reads as

$$Dh_p(u_a, u_b, v) = (\partial_1 h_{a,p}(u_a, v), \partial_1 h_{b,p}(u_b, v), \partial_2 h_{a,p}(u_a, v) + \partial_2 h_{b,p}(u_b, v)),$$

where

$$\partial_1 h_{a,p}(u_a, v) = \theta(u_a, v)^{\alpha(p-1)} u_a^{p-1}, \quad \partial_1 h_{b,p}(u_b, v) = \omega(u_b, v)^{\beta(p-1)} u_b^{p-1}, \quad (3.11)$$

and

$$\begin{aligned} \partial_2 h_{a,p}(u_a, v) &= c_{\text{fast}} \alpha(p-1) \int_0^{u_a} \theta(z, v)^{\alpha(p-1)-1} z^{p-1} dz, \\ \partial_2 h_{b,p}(u_b, v) &= d_{\text{fast}} \beta(p-1) \int_0^{u_b} \omega(z, v)^{\beta(p-1)-1} z^{p-1} dz. \end{aligned} \qquad (3.12)$$

The derivatives $\partial_1 h_{a,p}$ and $\partial_1 h_{b,p}$ in (3.11) are well defined for all $u_a, u_b, v \ge 0$ and all $A, B \ge 0$. The same holds true for $\partial_2 h_{a,p}$ and $\partial_2 h_{b,p}$ since the integrals in (3.12) are finite. Indeed, by (3.6) and (2.10), we have

$$\begin{aligned} 0 \le \theta(z, v)^{\alpha(p-1)-1} z^{p-1} &\le a_{\text{fast}}^{1-p} \theta(z, v)^{(1+\alpha)(p-1)-1} = a_{\text{fast}}^{1-p} \theta(z, v)^{q(p)-2}, \\ 0 \le \omega(z, v)^{\beta(p-1)-1} z^{p-1} &\le b_{\text{fast}}^{1-p} \omega(z, v)^{(1+\beta)(p-1)-1} = b_{\text{fast}}^{1-p} \omega(z, v)^{r(p)-2}, \end{aligned}$$

where $q(p), r(p) > 1$ as $p > 1$. Hence, by integration, we end up with

$$\begin{aligned} 0 \le \partial_2 h_{a,p}(u_a, v) &\lesssim a_{\text{fast}}^{-p} \theta(u_a, v)^{q(p)-1}, \\ 0 \le \partial_2 h_{b,p}(u_b, v) &\lesssim b_{\text{fast}}^{-p} \omega(u_b, v)^{r(p)-1}. \end{aligned} \qquad (3.13)$$
## 3.2 The Hessian of the total energy density $h_p$

Let $p > 1$. The Hessian matrix of $h_p$ is

$$ \text{Hess}(h_p) = \begin{pmatrix} \partial_{11}h_{a,p} & 0 & \partial_{12}h_{a,p} \\ 0 & 0 & 0 \\ \partial_{21}h_{a,p} & 0 & \partial_{22}h_{a,p} \end{pmatrix} + \begin{pmatrix} 0 & 0 & 0 \\ 0 & \partial_{11}h_{b,p} & \partial_{12}h_{b,p} \\ 0 & \partial_{21}h_{b,p} & \partial_{22}h_{b,p} \end{pmatrix}, $$

where

$$ \begin{aligned} \partial_{11}h_{a,p}(u_a, v) &= a_{\text{fast}}\alpha(p-1)\theta(u_a, v)^{\alpha(p-1)-1}u_a^{p-1} + (p-1)\theta(u_a, v)^{\alpha(p-1)-1}u_a^{p-2}, \\ \partial_{12}h_{a,p}(u_a, v) &= \partial_{21}h_{a,p}(u_a, v) = c_{\text{fast}}\alpha(p-1)\theta(u_a, v)^{\alpha(p-1)-1}u_a^{p-1}, \quad (3.14) \\ \partial_{22}h_{a,p}(u_a, v) &= c_{\text{fast}}^2\alpha(p-1)(\alpha(p-1)-1)\int_0^{u_a} \theta(z, v)^{\alpha(p-1)-2}z^{p-1}dz, \end{aligned} $$

and

$$ \begin{aligned} \partial_{11}h_{b,p}(u_b, v) &= b_{\text{fast}}\beta(p-1)\omega(u_b, v)^{\beta(p-1)-1}u_b^{p-1} + (p-1)\omega(u_b, v)^{\beta(p-1)-1}u_b^{p-2}, \\ \partial_{12}h_{b,p}(u_b, v) &= \partial_{21}h_{b,p}(u_b, v) = d_{\text{fast}}\beta(p-1)\omega(u_b, v)^{\beta(p-1)-1}u_b^{p-1}, \quad (3.15) \\ \partial_{22}h_{b,p}(u_b, v) &= d_{\text{fast}}^2\beta(p-1)(\beta(p-1)-1)\int_0^{u_b} \omega(z, v)^{\beta(p-1)-2}z^{p-1}dz. \end{aligned} $$

### 3.2.1 Hess($h_{a,p}$)

From (3.14), we see that, for $p \in (1, 2)$, we need a strictly positive density $u_a$ in order to define $u_a^{p-2}$ in $\partial_{11}h_{a,p}$. More precisely, this is necessary when $p = p_\alpha < 2$ and $p = p_\beta < 2$. The strict positivity of $u_a^\epsilon$ is given by Theorem 2.1.

Furthermore, the following cases will be considered to control the term appearing in both $\partial_{11}h_{a,p}$ and $\partial_{12}h_{a,p}$, i.e. $\theta(u_a, v)^{\alpha(p-1)-1}u_a^{p-1}$.

(a1) If $1 < p < p_\alpha$, then $\frac{1}{p-1} - \alpha > 1$ and, for all $u_a, v \ge 0$, as $A > 0$ (see (H1)), it holds

$$ 0 \le \theta(u_a, v)^{\alpha(p-1)-1} u_a^{p-1} = \left( \frac{u_a}{\theta(u_a, v)^{\frac{1}{p-1}-\alpha}} \right)^{p-1} \le C_A(p). \quad (3.16) $$

Estimate (3.16) will be crucial when $p = p_\beta$, since $p_\beta < p_\alpha$, if $\alpha < \beta$.

(a2) If $p \ge p_\alpha$, then $\alpha(p-1)-1 \ge -(p-1)$ and, for all $u_a, v \ge 0$, by (3.6) it holds

$$ 0 \le \theta(u_a, v)^{\alpha(p-1)-1} u_a^{p-1} \le a_{\text{fast}}^{1-p} \theta(u_a, v)^{(1+\alpha)(p-1)-1} = a_{\text{fast}}^{1-p} \theta(u_a, v)^{q(p)-2}, \quad (3.17) $$

with $q(p) \ge q(p_\alpha) = 2$, (see (2.12)).

Finally, as $A > 0$, $\partial_{22}h_{a,p}$ in (3.14) is well defined. In order to simplify the computations, we will handle $\partial_{22}h_{a,p}$ only for $p > p_\alpha$ and the cases to be analysed are the following.

(a3) If $p_\alpha < p < 1 + \frac{1}{\alpha}$, then $-(p-1) < \alpha(p-1) - 1 < 0$ and by (3.6), we have

$$ 0 \le \theta(z, v)^{\alpha(p-1)-2} z^{p-1} \le a_{\text{fast}}^{1-p} \theta(z, v)^{(1+\alpha)(p-1)-2} = a_{\text{fast}}^{1-p} \theta(z, v)^{q(p)-3}. $$
Hence

$$ \begin{gathered} |\partial_{22}h_{a,p}(u_a, v)| \lesssim \alpha(p-1) C(\alpha, p) a_{\text{fast}}^{-p} \theta(u_a, v)^{q(p)-2}, \\ C(\alpha, p) := \frac{1 - \alpha(p-1)}{(\alpha+1)(p-1)-1} > 0. \end{gathered} \tag{3.18} $$

(a4) If $p \ge 1 + \frac{1}{\alpha}$, then $\alpha(p-1) - 1 \ge 0$, so that $\partial_{22}h_{a,p}$ is positive and gives a negative term (see (5.3)) in the evolution equation (3.8) that will be neglected.

### 3.2.2 Hess($h_{b,p}$)

A similar analysis can be done for $\partial_{11}h_{b,p}$ and $\partial_{12}h_{b,p}$ in (3.15). As before, we see that, for $p \in (1, 2)$, we need a strictly positive density $u_b$, in order to define $u_b^{p-2}$. On the other hand, the critical case corresponding to the (a1) case above will not appear, since we will not consider any $p \in (1, p_\beta)$. This is one of the reason why, we do not need to assume $\phi$ strictly positive for all $\beta > 0$ (i.e. $B > 0$ for all $\beta > 0$) as we do for $\psi$, (see (H1)). However, we need to make sure that $\partial_{22}h_{b,p}$ is well defined. To this end, we proceed as follows.

First, concerning $\partial_{11}h_{b,p}$ and $\partial_{12}h_{b,p}$, we observe, as in the (a2) case, that

(b1) if $p \ge p_\beta$, then $\beta(p-1) - 1 \ge -(p-1)$ and, for all $u_b, v \ge 0$, it holds

$$ 0 \le \omega(u_b, v)^{\beta(p-1)-1} u_b^{p-1} \le b_{\text{fast}}^{1-p} \omega(u_b, v)^{(1+\beta)(p-1)-1} = b_{\text{fast}}^{1-p} \omega(u_b, v)^{r(p)-2}, \tag{3.19} $$

with $r(p) \ge r(p_\beta) = 2$, (see (2.12)).

Next, as $B \ge 0$ when $\beta \ge 1$, $\partial_{22}h_{b,p}$ is not well defined when $p = p_\beta$, since the integrability of the function $\omega(z, v)^{\beta(p-1)-2}z^{p-1}$ in the neighbourhood of $z = 0^+$ is not guaranteed, as it holds

$$ \omega(z, v)^{\beta(p_\beta-1)-2} z^{p_\beta-1} = \omega(z, v)^{-1-\frac{1}{\beta+1}} z^{\frac{1}{\beta+1}}. $$

Therefore, we will avoid this criticality and the following cases are considered.

(b2) If $p_\beta < p < 1 + \frac{1}{\beta}$ then $-(p-1) < \beta(p-1) - 1 < 0$ and we have by (3.6)

$$ 0 \le \omega(z, v)^{\beta(p-1)-2} z^{p-1} \le b_{\text{fast}}^{1-p} \omega(z, v)^{(1+\beta)(p-1)-2} = b_{\text{fast}}^{1-p} \omega(z, v)^{r(p)-3}. $$

Hence, as $r(p) > r(p_\beta) = 2$, $\partial_{22}h_{b,p}$ is well defined and

$$ \begin{gathered} |\partial_{22}h_{b,p}(u_b, v)| \lesssim \beta(p-1) C(\beta, p) b_{\text{fast}}^{-p} \omega(u_b, v)^{r(p)-2}, \\ C(\beta, p) := \frac{1 - \beta(p-1)}{(\beta+1)(p-1)-1} > 0. \end{gathered} \tag{3.20} $$

(b3) If $p \ge 1 + \frac{1}{\beta}$, then $\beta(p-1) - 1 \ge 0$ and $\partial_{22}h_{b,p}$ is positive and gives a negative term (see (5.3)) in the evolution equation (3.8) that will be neglected.

## 4 The reaction contribution $I_{\text{rea}}^p$ to $\mathcal{E}_p'(t)$ in (3.8)

In the sequel, it will be useful to employ the following elementary interpolation inequality: for all $x \ge 0$, $C > 0$ and $\gamma, \gamma_1 \in \mathbb{R}$ such that $0 < \gamma < \gamma_1$, it holds

$$ x^{\gamma} \le C^{\gamma} + x^{\gamma_1} C^{\gamma - \gamma_1}. \tag{4.1} $$
Let $p > 1$. From (3.8) we have that $I_{\text{rea}}^p = \int_{\Omega} Dh_p \cdot \mathcal{F} dx$, i.e.

$$ 
\begin{aligned}
I_{\text{rea}}^p &= \int_{\Omega} \partial_1 h_{a,p}(u_a, v) f_a(u_a, u_b, v) dx + \int_{\Omega} \partial_1 h_{b,p}(u_b, v) f_b(u_a, u_b, v) dx \\
&\quad + \int_{\Omega} \partial_2 h_{a,p}(u_a, v) f_v(u_a, u_b, v) dx + \int_{\Omega} \partial_2 h_{b,p}(u_b, v) f_v(u_a, u_b, v) dx \\
&:= J_1^p + J_2^p + J_3^p + J_4^p.
\end{aligned} 
\quad (4.2) $$

The competitive expression of the reaction functions $f_a$ and $f_b$ enable us to obtain estimates of $\|u_a\|_{L^{q(p)+1}(\Omega_T)}$ and $\|u_b\|_{L^{r(p)+1}(\Omega_T)}$ from $J_1^p$ and $J_2^p$, respectively (see (4.5) and (4.6)) and to absorb some terms arising from the diffusion. Indeed, using (3.11), the definitions of $f_a$ and $\theta$ in (1.4), (3.4), assumption (H3) and neglecting the non-negative intra competition term $u_a u_b$ in $f_a$, it holds

$$ 
\begin{aligned}
J_1^p &= \int_{\Omega} \partial_1 h_{a,p}(u_a, v) f_a(u_a, u_b, v) dx = \int_{\Omega} \theta(u_a, v)^{\alpha(p-1)} u_a^{p-1} f_a(u_a, u_b, v) dx \\
&\le \eta_a \int_{\Omega} \theta(u_a, v)^{\alpha(p-1)} u_a^{p-1} u_a(1+A - (A + a_{\text{fast}}u_a + c_{\text{fast}}v)) dx \\
&= \eta_a \int_{\Omega} u_a^p \left[ \theta(u_a, v)^{\alpha(p-1)} (1+A - \theta(u_a, v)) \right] dx.
\end{aligned} 
\quad (4.3) $$

Next, taking in (4.1), $\gamma_1 = \gamma + 1$ and $C > 1$, we have the inequality

$$ x^{\gamma}(1-x) \le C^{\gamma} - (1-C^{-1})x^{\gamma+1}. \quad (4.4) $$

Applying (4.4) with $x = \frac{\theta(u_a, v)}{1+A}$, $\gamma = \alpha(p-1)$ and $C = 1+A > 1$, to the function in the square brackets in (4.3), we have

$$ \theta(u_a, v)^{\alpha(p-1)}(1+A - \theta(u_a, v)) = (1+A)^{\gamma+1} x^{\gamma}(1-x) \le C^{2\gamma+1} - \frac{A}{1+A} \theta(u_a, v)^{\gamma+1}. $$

Then using (3.6) and (2.10), we end up with

$$ J_1^p \lesssim (1+A)^{2\alpha(p-1)+1} \|u_a\|_{L^p(\Omega)}^p - a_{\text{fast}}^{\alpha(p-1)+1} \frac{A}{1+A} \|u_a\|_{L^{q(p)+1}(\Omega)}^{q(p)+1}. \quad (4.5) $$

Similarly, for $J_2^p$ we obtain that

$$ 
\begin{aligned}
J_2^p &= \int_{\Omega} \omega(u_b, v)^{\beta(p-1)} u_b^{p-1} f_b(u_a, u_b, v) dx \\
&\le \eta_b \int_{\Omega} u_b^p \left[ \omega(u_b, v)^{\beta(p-1)} (1 + B - \omega(u_b, v)) \right] dx.
\end{aligned} 
$$

However, as $B \ge 0$ when $\beta \ge 1$, we choose an arbitrary $\sigma > 0$, we replace $B$ with $B \lor \sigma$ and, proceeding as before, we have

$$ J_2^p \lesssim (1 + B \lor \sigma)^{2\beta(p-1)+1} \|u_b\|_{L^p(\Omega)}^p - b_{\text{fast}}^{\beta(p-1)+1} \frac{B \lor \sigma}{1 + B \lor \sigma} \|u_b\|_{L^{r(p)+1}(\Omega)}^{r(p)+1}. \quad (4.6) $$

The terms $J_3^p$ and $J_4^p$ in (4.2) cannot give similar estimates since they contain the interaction between $\theta$ and $\omega$, through the reaction function $f_v$, so that (4.4) can
not be applied. Therefore, we simply use (3.13) and neglect all the negative terms in $f_v$, to end up with

$$J_3^p = \int_{\Omega} \partial_2 h_{a,p}(u_a, v) f_v(u_a, u_b, v) dx \lesssim a_{\text{fast}}^{-p} \int_{\Omega} v \theta(u_a, v)^{q(p)-1} dx, \quad (4.7)$$

and

$$J_4^p = \int_{\Omega} \partial_2 h_{b,p}(u_b, v) f_v(u_a, u_b, v) dx \lesssim b_{\text{fast}}^{-p} \int_{\Omega} v \omega(u_b, v)^{r(p)-1} dx. \quad (4.8)$$

Finally, plugging (4.5)-(4.8) into (4.2), we have for all $p > 1$

$$\begin{aligned} I_{\text{rea}}^p \lesssim & (1+A)^{2\alpha(p-1)+1} \|u_a\|_{L^p(\Omega)}^p - a_{\text{fast}}^{\alpha(p-1)+1} \frac{A}{1+A} \|u_a\|_{L^{q(p)+1}(\Omega)}^{q(p)+1} \\ & + (1+B \vee \sigma)^{2\beta(p-1)+1} \|u_b\|_{L^p(\Omega)}^p - b_{\text{fast}}^{\beta(p-1)+1} \frac{B \vee \sigma}{1+B \vee \sigma} \|u_b\|_{L^{r(p)+1}(\Omega)}^{r(p)+1} \\ & + a_{\text{fast}}^{-p} \|v \theta(u_a, v)^{q(p)-1}\|_{L^1(\Omega)} + b_{\text{fast}}^{-p} \|v \omega(u_b, v)^{r(p)-1}\|_{L^1(\Omega)}. \end{aligned} \quad (4.9)$$

# 5 The diffusion contribution $I_{\text{diff}}^p$ to $\mathcal{E}_p'(t)$ in (3.8)

From (3.8) we have

$$I_{\text{diff}}^p = d_a \int_{\Omega} \partial_1 h_p \Delta u_a dx + d_b \int_{\Omega} \partial_2 h_p \Delta u_b dx + d_v \int_{\Omega} \partial_3 h_p \Delta v dx,$$

and, by definition (3.2),

$$\begin{aligned} I_{\text{diff}}^p = & d_a \int_{\Omega} \partial_1 h_{a,p} \Delta u_a dx + d_b \int_{\Omega} \partial_1 h_{b,p} \Delta u_b dx \\ & + d_v \int_{\Omega} \partial_2 h_{a,p} \Delta v dx + d_v \int_{\Omega} \partial_2 h_{b,p} \Delta v dx. \end{aligned} \quad (5.1)$$

As it is not possible to have a priori estimates on $\Delta u_a$ and $\Delta u_b$ uniform in $\varepsilon$, we have to apply Green's formula to the first and second integral in the right hand side of (5.1). Assumption (H1) appears to be fundamental here, since we need to control $\partial_{11}h_{a,p}$ and $\partial_{12}h_{a,p}$, for $p < p_\alpha$ (see (3.14) and (3.16)).

On the other hand, as we do not have assumed the strict positivity of the transition function $\phi$ when $\beta \ge 1$, we cannot bound $\partial_{22}h_{b,p}$ for $p = p_\beta$, when $\beta \ge 1$, (see Subsection 3.2.2). Hence, the third and fourth integral in the right hand side of (5.1) are left as they are for $p \le p_\alpha$ and $p \le p_\beta$, respectively, and the bound on $\Delta v$ used (see (2.8), (2.19)). Green's formula will be used for these two terms when $p > p_\alpha$ and $p > p_\beta$, respectively.

To resume, using the boundary conditions (1.2), and the Heaviside functions

$$\chi_{\beta}(p) := \begin{cases} 1, & p > p_{\beta}, \\ 0, & p \le p_{\beta}, \end{cases} \quad \chi_{\alpha}(p) := \begin{cases} 1, & p > p_{\alpha}, \\ 0, & p \le p_{\alpha}, \end{cases} \quad (5.2)$$
$I_{\text{diff}}^p$ rewrites as

$$
\begin{align}
I_{\text{diff}}^{p} = & -d_a \int_{\Omega} \partial_{11} h_{a,p} |\nabla u_a|^2 dx - (d_a + \chi_{\alpha}(p) d_v) \int_{\Omega} \partial_{12} h_{a,p} \nabla u_a \cdot \nabla v \, dx \nonumber \\
& - d_b \int_{\Omega} \partial_{11} h_{b,p} |\nabla u_b|^2 dx - (d_b + \chi_{\beta}(p) d_v) \int_{\Omega} \partial_{12} h_{b,p} \nabla u_b \cdot \nabla v \, dx \nonumber \\
& - d_v \chi_{\alpha}(p) \int_{\Omega} \partial_{22} h_{a,p} |\nabla v|^2 dx - d_v \chi_{\beta}(p) \int_{\Omega} \partial_{22} h_{b,p} |\nabla v|^2 dx \tag{5.3} \\
& + d_v (1 - \chi_{\alpha}(p)) \int_{\Omega} \partial_2 h_{a,p} \Delta v \, dx + d_v (1 - \chi_{\beta}(p)) \int_{\Omega} \partial_2 h_{b,p} \Delta v \, dx \nonumber \\
& := K_1^p + K_2^p + K_3^p + K_4^p + K_5^p + K_6^p. \nonumber
\end{align}
$$

In the rest of the section, we will estimate each of the $K_i^p$ terms above.

**Estimate of $K_1^p$.** From (3.14), $K_1^p$ reads as

$$
\begin{align}
K_1^p = & -d_a \int_{\Omega} \partial_{11} h_{a,p} |\nabla u_a|^2 dx - (d_a + \chi_{\alpha}(p) d_v) \int_{\Omega} \partial_{12} h_{a,p} \nabla u_a \cdot \nabla v \, dx \nonumber \\
= & -d_a a_{\text{fast}} \alpha(p-1) \int_{\Omega} \theta(u_a, v)^{\alpha(p-1)-1} u_a^{p-1} |\nabla u_a|^2 \, dx \nonumber \\
& - d_a(p-1) \int_{\Omega} \theta(u_a, v)^{\alpha(p-1)} u_a^{p-2} |\nabla u_a|^2 \, dx \nonumber \\
& - (d_a + \chi_{\alpha}(p) d_v) c_{\text{fast}} \alpha(p-1) \int_{\Omega} \theta(u_a, v)^{\alpha(p-1)-1} u_a^{p-1} \nabla u_a \cdot \nabla v \, dx. \tag{5.4}
\end{align}
$$

Then, by Young's inequality applied to the third integral in the right hand side of (5.4), we have

$$
\begin{align}
K_1^p \lesssim & -\alpha(p-1) \int_{\Omega} \theta(u_a, v)^{\alpha(p-1)-1} u_a^{p-1} |\nabla u_a|^2 \, dx \nonumber \\
& - d_a(p-1) \int_{\Omega} \theta(u_a, v)^{\alpha(p-1)} u_a^{p-2} |\nabla u_a|^2 \, dx \tag{5.5} \\
& + \alpha(p-1) \int_{\Omega} \theta(u_a, v)^{\alpha(p-1)-1} u_a^{p-1} |\nabla v|^2 \, dx. \nonumber
\end{align}
$$

Next, we use (3.6) in the second integral in the right hand side of (5.5) to have

$$
\theta(u_a, v)^{\alpha(p-1)} u_a^{p-2} \geq a_{\text{fast}}^{\alpha(p-1)} u_a^{(\alpha+1)(p-1)-1} = a_{\text{fast}}^{\alpha(p-1)} u_a^{q(p)-2}.
$$

Hence, neglecting the first integral, we obtain

$$
\begin{align}
K_1^p \lesssim & -d_a(p-1)a_{\text{fast}}^{\alpha(p-1)} \int_{\Omega} u_a^{q(p)-2} |\nabla u_a|^2 \, dx \nonumber \\
& + \alpha(p-1) \int_{\Omega} \theta(u_a, v)^{\alpha(p-1)-1} u_a^{p-1} |\nabla v|^2 \, dx. \tag{5.6}
\end{align}
$$

Finally, for the second integral in the right hand side of (5.6), we proceed according the value of $p$

(i) if $1 < p < p_\alpha$, we employ (3.16) to obtain

$$
K_1^p \lesssim -d_a \frac{4(p-1)}{q(p)^2} a_{\text{fast}}^{\alpha(p-1)} \|\nabla u_a^{q(p)/2}\|_{L^2(\Omega)}^2 + \alpha(p-1) C_A(p) \|\nabla v\|_{L^2(\Omega)}^2; \quad (5.7)
$$
(ii) if $p \ge p_\alpha$, we employ (3.17) and we have

$$ K_1^p \lesssim -d_a \frac{4(p-1)}{q(p)^2} a_{\text{fast}}^{\alpha(p-1)} \|\nabla u_a^{q(p)/2}\|_{L^2(\Omega)}^2 + \alpha(p-1) a_{\text{fast}}^{1-p} \|\theta(u_a, v)^{q(p)/2-1} \nabla v\|_{L^2(\Omega)}^2. \quad (5.8) $$

**Estimate of $K_2^p$.** From (3.15), $K_2^p$ reads as

$$ \begin{aligned} K_2^p &= -d_b \int_{\Omega} \partial_{11} h_{b,p} |\nabla u_b|^2 dx - (d_b + \chi_\beta(p) d_v) \int_{\Omega} \partial_{12} h_{b,p} \nabla u_b \cdot \nabla v dx \\ &= -d_b b_{\text{fast}} \beta(p-1) \int_{\Omega} \omega(u_b, v)^{\beta(p-1)-1} u_b^{p-1} |\nabla u_b|^2 dx \\ &\quad - d_b(p-1) \int_{\Omega} \omega(u_b, v)^{\beta(p-1)} u_b^{p-2} |\nabla u_b|^2 dx \\ &\quad - (d_b + \chi_\beta(p) d_v) d_{\text{fast}} \beta(p-1) \int_{\Omega} \omega(u_b, v)^{\beta(p-1)-1} u_b^{p-1} \nabla u_b \cdot \nabla v dx. \end{aligned} $$

Again by Young's inequality we have

$$ \begin{aligned} K_2^p \lesssim & -\beta(p-1) \int_{\Omega} \omega(u_b, v)^{\beta(p-1)-1} u_b^{p-1} |\nabla u_b|^2 dx \\ & - d_b(p-1) \int_{\Omega} \omega(u_b, v)^{\beta(p-1)} u_b^{p-2} |\nabla u_b|^2 dx \\ & + \beta(p-1) \int_{\Omega} \omega(u_b, v)^{\beta(p-1)-1} u_b^{p-1} |\nabla v|^2 dx. \end{aligned} $$

Then, as before, we neglect the first integral, we use (3.6) in the second integral and (3.19) in the third one, to end up with, for all $p \ge p_\beta$,

$$ K_2^p \lesssim -d_b \frac{4(p-1)}{r(p)^2} b_{\text{fast}}^{\beta(p-1)} \|\nabla u_b^{r(p)/2}\|_{L^2(\Omega)}^2 + \beta(p-1) b_{\text{fast}}^{1-p} \|\omega(u_b, v)^{r(p)/2-1} \nabla v\|_{L^2(\Omega)}^2. \quad (5.9) $$

**Estimate of $K_3^p$ and $K_4^p$.** It is sufficient to estimate $K_3^p$ for $p_\alpha < p < 1 + \alpha^{-1}$, since out of this range of $p$, $K_3^p$ is either zero or nonpositive (see (3.14)) and can be neglected. From (3.18) it holds

$$ \begin{aligned} 0 < K_3^p &= -d_v \chi_\alpha(p) \int_{\Omega} \partial_{22} h_{a,p} |\nabla v|^2 dx \\ &\lesssim \alpha(p-1) C(\alpha, p) a_{\text{fast}}^{-p} \|\theta(u_a, v)^{q(p)/2-1} \nabla v\|_{L^2(\Omega)}^2, \quad p_\alpha < p < 1 + \alpha^{-1}. \end{aligned} \quad (5.10) $$

Similarly, $K_4^p$ is strictly positive for $p_\beta < p < 1 + \beta^{-1}$ (see (3.15)), and in this range, by (3.20), it holds

$$ \begin{aligned} 0 < K_4^p &= -d_v \chi_\beta(p) \int_{\Omega} \partial_{22} h_{b,p} |\nabla v|^2 dx \\ &\lesssim \beta(p-1) C(\beta, p) b_{\text{fast}}^{-p} \|\omega(u_b, v)^{r(p)/2-1} \nabla v\|_{L^2(\Omega)}^2, \quad p_\beta < p < 1 + \beta^{-1}. \end{aligned} \quad (5.11) $$

**Estimate of $K_5^p$ and $K_6^p$.** The term $K_5^p$ is not trivial for $p \le p_\alpha$ and it will be employed only for $p = p_\beta$ and $p = p_\alpha$. Thus, using (3.13) and (2.10), we obtain

$$ \begin{aligned} 0 < K_5^{p_\beta} &= d_v \int_{\Omega} \partial_2 h_{a,p_\beta} \Delta v \, dx \lesssim a_{\text{fast}}^{-p_\beta} \int_{\Omega} \theta(u_a, v)^{\frac{\alpha+1}{\beta+1}} |\Delta v| \, dx \\ &\lesssim \|\theta(u_a, v)^{q(p_\beta)-1}\|_{L^2(\Omega)} \|\Delta v\|_{L^2(\Omega)}, \end{aligned} \quad (5.12) $$
and

$$0 < K_5^{p_\alpha} \lesssim \|\theta(u_a, v)\|_{L^2(\Omega)} \|\Delta v\|_{L^2(\Omega)}. \quad (5.13)$$

Finally, the term $K_6^p$ is not trivial for $p \le p_\beta$ and it will be employed only for $p = p_\beta$. Using (3.13), we have

$$0 < K_6^{p_\beta} = d_v \int_{\Omega} \partial_2 h_{b,p_\beta} \Delta v \, dx \lesssim \|\omega(u_b, v)\|_{L^2(\Omega)} \|\Delta v\|_{L^2(\Omega)}. \quad (5.14)$$

# 6 Energy estimates: proof of Lemma 2.2

This section is devoted to the proof of Lemma 2.2, based on the computations obtained in Sections 3 to 5, on the maximal regularity (A.17) giving (A.18) and on a bootstrap argument. We recall that $\alpha \le \beta$, so that $p_\beta \le p_\alpha$, and that $q(p_\alpha) = r(p_\beta) = 2$. Therefore, we estimate $\mathcal{E}_p$ along the trajectories of the solution $(u_a^\varepsilon, u_b^\varepsilon, v^\varepsilon)$, starting with $p = p_\beta$, then $p = p_\alpha$ and finally $p \in [2, p_{\alpha,\beta}^0)$ (see (6.21)), using the differential equation below (see (3.8) and (5.3))

$$\mathcal{E}'_p(t) = I_{\text{diff}}^p + I_{\text{rea}}^p + I_{\text{fast}}^p = \sum_{i=1}^6 K_i^p + I_{\text{rea}}^p + I_{\text{fast}}^p. \quad (6.1)$$

## 6.1 Estimates from $\mathcal{E}_{p_\beta}$, $\alpha < \beta$

Assume $\alpha < \beta$. Taking into account that $K_3^{p_\beta} = K_4^{p_\beta} = 0$ (see (5.10), (5.11)), from estimates (5.7), (5.9) where we use $r(p_\beta) = 2$, (5.12) and (5.14), we have

$$I_{\text{diff}}^{p_\beta} \lesssim - \|\nabla u_a^{q(p_\beta)/2}\|_{L^2(\Omega)}^2 - \|\nabla u_b\|_{L^2(\Omega)}^2 + \|\nabla v\|_{L^2(\Omega)}^2 \\ + \|\theta(u_a, v)^{q(p_\beta)-1}\|_{L^2(\Omega)} \|\Delta v\|_{L^2(\Omega)} + \|\omega(u_b, v)\|_{L^2(\Omega)} \|\Delta v\|_{L^2(\Omega)}. \quad (6.2)$$

Using (3.9), plugging (6.2) and estimate (4.9) of $I_{\text{rea}}^{p_\beta}$ into (6.1), rearranging the terms and integrating the obtained inequality over $(0, T)$, we get

$$\begin{aligned} \mathcal{E}_{p_\beta}(T) - \mathcal{E}_{p_\beta}(0) &+ \|\nabla u_a^{q(p_\beta)/2}\|_{L^2(\Omega_T)}^2 + \|\nabla u_b\|_{L^2(\Omega_T)}^2 + \|u_a\|_{L^{q(p_\beta)+1}(\Omega_T)}^{q(p_\beta)+1} + \|u_b\|_{L^3(\Omega_T)}^3 \\ &\lesssim \|\nabla v\|_{L^2(\Omega_T)}^2 + \|\Delta v\|_{L^2(\Omega_T)}^2 + \|\theta(u_a, v)^{q(p_\beta)-1}\|_{L^2(\Omega_T)}^2 + \|\omega(u_b, v)\|_{L^2(\Omega_T)}^2 \\ &+ \|u_a\|_{L^{p_\beta}(\Omega_T)}^{p_\beta} + \|u_b\|_{L^{p_\beta}(\Omega_T)}^{p_\beta} + \|v\theta(u_a, v)^{q(p_\beta)-1}\|_{L^1(\Omega_T)} + \|v\omega(u_b, v)\|_{L^1(\Omega_T)}^2. \end{aligned} \quad (6.3)$$

As $p_\beta < 2$ and $q(p_\beta) - 1 = \frac{\alpha+1}{\beta+1} < 1$, recalling that $\theta$ and $\omega$ are affine functions, the estimates obtained in Theorem 2.1 allow us to control all the terms in the right hand side of (6.3). Hence, using (3.7), we end up with

$$\begin{aligned} \|u_a\|_{L^\infty(0,T;L^{q(p_\beta)}(\Omega))}^{q(p_\beta)} + \|u_b\|_{L^\infty(0,T;L^2(\Omega))}^2 + \|\nabla u_a^{q(p_\beta)/2}\|_{L^2(\Omega_T)}^2 + \|\nabla u_b\|_{L^2(\Omega_T)}^2 \\ + \|u_a\|_{L^{q(p_\beta)+1}(\Omega_T)}^{q(p_\beta)+1} + \|u_b\|_{L^3(\Omega_T)}^3 \lesssim \mathcal{E}_{p_\beta}(0) + C(|\Omega|, T). \end{aligned} \quad (6.4)$$
## 6.2 Estimates from $\mathcal{E}_{p_\alpha}$, $\alpha < \beta$

Assume $\alpha < \beta$. Taking into account that $K_3^{p_\alpha} = K_6^{p_\alpha} = 0$, using (5.8) with $q(p_\alpha) = 2$, (5.9), (5.11), (5.13), we have

$$ I_{\text{diff}}^{p_\alpha} \lesssim -\|\nabla u_a\|_{L^2(\Omega)}^2 - \|\nabla u_b^{r(p_\alpha)/2}\|_{L^2(\Omega)}^2 + \|\nabla v\|_{L^2(\Omega)}^2 \\ + \|\omega(u_b, v)^{r(p_\alpha)/2-1} \nabla v\|_{L^2(\Omega)}^2 + \|\theta(u_a, v)\|_{L^2(\Omega)}\|\Delta v\|_{L^2(\Omega)}. \quad (6.5) $$

Using (3.9), plugging (6.5) and estimate (4.9) of $I_{\text{rea}}^{p_\alpha}$ into (6.1), rearranging the terms and integrating the obtained inequality over $(0, T)$, we get

$$ \begin{aligned} \mathcal{E}_{p_\alpha}(T) - \mathcal{E}_{p_\alpha}(0) &+ \|\nabla u_a\|_{L^2(\Omega_T)}^2 + \|\nabla u_b^{r(p_\alpha)/2}\|_{L^2(\Omega_T)}^2 + \|u_a\|_{L^3(\Omega_T)}^3 + \|u_b\|_{L^{r(p_\alpha)+1}(\Omega_T)}^{r(p_\alpha)+1} \\ &\lesssim \|\nabla v\|_{L^2(\Omega_T)}^2 + \|\omega(u_b, v)^{r(p_\alpha)/2-1} \nabla v\|_{L^2(\Omega_T)}^2 + \|\theta(u_a, v)\|_{L^2(\Omega_T)}^2 + \|\Delta v\|_{L^2(\Omega_T)}^2 \\ &+ \|u_a\|_{L^{p_\alpha}(\Omega_T)}^{p_\alpha} + \|u_b\|_{L^{p_\alpha}(\Omega_T)}^{p_\alpha} + \|v\theta(u_a, v)\|_{L^1(\Omega_T)} + \|v\omega(u_b, v)^{r(p_\alpha)-1}\|_{L^1(\Omega_T)} \end{aligned} \quad (6.6) $$

It is worth noticing that despite $p_\alpha < 2$, $r(p_\alpha)$ can be large without any restriction on $\beta - \alpha$, since from (2.10),(2.11) and $\alpha < \beta$, it holds

$$ r(p_\alpha) = \frac{\beta + 1}{\alpha + 1} + 1 = \frac{\beta - \alpha}{\alpha + 1} + 2 > 2. \quad (6.7) $$

Hence, in order to obtain new a priori estimates on $u_a$ and $u_b$ from (6.6), we need to get rid of the terms

$$ I := \|\omega(u_b, v)^{r(p_\alpha)/2-1} \nabla v\|_{L^2(\Omega_T)}^2 \quad \text{and} \quad J := \|v \omega(u_b, v)^{r(p_\alpha)-1}\|_{L^1(\Omega_T)} \quad (6.8) $$

Let $\delta \in (0, 1)$. Applying Young's inequality into $I$ in (6.8) we have

$$ I \lesssim \delta \int_{\Omega_T} |\nabla v|^6 dx dt + \delta^{-\frac{1}{2}} \int_{\Omega_T} \omega(u_b, v)^{\frac{3}{2}(r(p_\alpha)-2)} dx dt =: I_1 + I_2. \quad (6.9) $$

Then, by (A.18) with $p=3$,

$$ I_1 = \delta \|\nabla v\|_{L^6(\Omega_T)}^6 \lesssim \delta(1+T) + \delta \|u_a\|_{L^3(\Omega_T)}^3 + \delta \|u_b\|_{L^3(\Omega_T)}^3. \quad (6.10) $$

On the other hand, by (6.7) and assumption (H2),

$$ 0 < \frac{3}{2}(r(p_\alpha) - 2) = \frac{3}{2}\frac{\beta - \alpha}{\alpha + 1} < \frac{\beta - \alpha}{\alpha + 1} + 3 = r(p_\alpha) + 1. $$

Using (4.1) with $\gamma = \frac{3}{2}(r(p_\alpha) - 2)$, $\gamma_1 = r(p_\alpha) + 1$ and $C = C(\delta) > 0$ such that $\delta^{-\frac{1}{2}}C(\delta)^{\gamma-\gamma_1} = \delta$, we obtain

$$ I_2 = \delta^{-\frac{1}{2}} \int_{\Omega_T} \omega(u_b, v)^{\frac{3}{2}(r(p_\alpha)-2)} dx dt \le \delta C(\delta)^{r(p_\alpha)+1} |\Omega_T| + \delta \int_{\Omega_T} \omega(u_b, v)^{r(p_\alpha)+1} dx dt. $$

Finally, recalling definition (3.4) of $\omega$, by Jensen's inequality and the boundedness of $v$ (see (2.6)), we have

$$ I_2 \lesssim \delta C(\delta)^{r(p_\alpha)+1} |\Omega_T| + \delta |\Omega_T| + \delta \|u_b\|_{L^{r(p_\alpha)+1}(\Omega_T)}^{r(p_\alpha)+1}. \quad (6.11) $$
Plugging (6.10) and (6.11) into (6.9), we obtain

$$I \lesssim \delta(1 + C(\delta, |\Omega|, T)) + \delta \|u_a\|_{L^3(\Omega_T)}^3 + \delta \|u_b\|_{L^3(\Omega_T)}^3 + \delta \|u_b\|_{L^{r(p_\alpha)+1}(\Omega_T)}^{r(p_\alpha)+1}. \quad (6.12)$$

Next, by (6.7), $r(p_\alpha) - 1 > 1$. Hence, using again Jensen's inequality and the boundedness of $v$ in $J$ defined in (6.8), we have

$$J = \int_{\Omega_T} v \omega(u_b, v)^{r(p_\alpha)-1} dx dt \lesssim \int_{\Omega_T} u_b^{r(p_\alpha)-1} dx dt + |\Omega_T|. \quad (6.13)$$

Using (4.1) with $\gamma = r(p_\alpha) - 1$, $\gamma_1 = r(p_\alpha) + 1$ and $C = \delta^{-\frac{1}{2}}$ so that $C^{\gamma - \gamma_1} = \delta$, we end up with

$$J \lesssim \delta^{\frac{1-r(p_\alpha)}{2}} |\Omega_T| + \delta \|u_b\|_{L^{r(p_\alpha)+1}(\Omega_T)}^{r(p_\alpha)+1} + |\Omega_T|. \quad (6.14)$$

To conclude, (6.12) and (6.14) imply that there exists $C(\delta, |\Omega|, T) > 0$ such that

$$I + J \lesssim C(\delta, |\Omega|, T) + \delta \|u_a\|_{L^3(\Omega_T)}^3 + \delta \|u_b\|_{L^3(\Omega_T)}^3 + \delta \|u_b\|_{L^{r(p_\alpha)+1}(\Omega_T)}^{r(p_\alpha)+1}. \quad (6.15)$$

Choosing $\delta$ small enough, plugging (6.15) into (6.6) and rearranging the terms, we obtain

$$\begin{aligned} \mathcal{E}_{p_\alpha}(T) - \mathcal{E}_{p_\alpha}(0) &+ \|\nabla u_a\|_{L^2(\Omega_T)}^2 + \|\nabla u_b^{r(p_\alpha)/2}\|_{L^2(\Omega_T)}^2 \\ &+ (1-\delta)\|u_a\|_{L^3(\Omega_T)}^3 + (1-\delta)\|u_b\|_{L^{r(p_\alpha)+1}(\Omega_T)}^{r(p_\alpha)+1} \\ &\lesssim C(\delta, |\Omega|, T) + \|\nabla v\|_{L^2(\Omega_T)}^2 + \|\theta(u_a, v)\|_{L^2(\Omega_T)}^2 + \|\Delta v\|_{L^2(\Omega_T)}^2 \\ &+ \|u_a\|_{L^{p_\alpha}(\Omega_T)}^{p_\alpha} + \|u_b\|_{L^{p_\alpha}(\Omega_T)}^{p_\alpha} + \|v\theta(u_a, v)\|_{L^1(\Omega_T)} + \delta \|u_b\|_{L^3(\Omega_T)}^3. \end{aligned} \quad (6.16)$$

Recalling that $p_\alpha < 2$, the estimates obtained in Theorem 2.1 plus the estimate of $\|u_b\|_{L^3(\Omega_T)}$ obtained in (6.4), will allow us to control all the terms in the right hand side of (6.16), so that, using (3.7), we get

$$\begin{aligned} \|u_a\|_{L^\infty(0,T;L^2(\Omega))}^2 + \|u_b\|_{L^\infty(0,T;L^{r(p_\alpha)}(\Omega))}^{r(p_\alpha)} + \|\nabla u_a\|_{L^2(\Omega_T)}^2 + \|\nabla u_b^{r(p_\alpha)/2}\|_{L^2(\Omega_T)}^2 \\ + \|u_a\|_{L^3(\Omega_T)}^3 + \|u_b\|_{L^{r(p_\alpha)+1}(\Omega_T)}^{r(p_\alpha)+1} \lesssim \mathcal{E}_{p_\alpha}(0) + C(|\Omega|, T). \end{aligned} \quad (6.17)$$

## 6.3 Estimates from $\mathcal{E}_{p_\alpha} = \mathcal{E}_{p_\beta}$, $\alpha = \beta$

If $\alpha = \beta$, then $p_\alpha = p_\beta$, $K_3^{p_\alpha} = K_4^{p_\alpha} = 0$ and we can use (5.8), (5.9), (5.13), (5.14), to obtain

$$\begin{aligned} I_{\text{diff}}^{p_\beta} \lesssim & -\|\nabla u_a\|_{L^2(\Omega)}^2 - \|\nabla u_b\|_{L^2(\Omega)}^2 + \|\nabla v\|_{L^2(\Omega)}^2 \\ &+ \|\theta(u_a, v)\|_{L^2(\Omega)} \|\Delta v\|_{L^2(\Omega)} + \|\omega(u_b, v)\|_{L^2(\Omega)} \|\Delta v\|_{L^2(\Omega)}. \end{aligned} \quad (6.18)$$

Employing again (3.9) and estimate (4.9) for $I_{\text{rea}}^{p_\beta} = I_{\text{rea}}^{p_\alpha}$, plugging (6.18) into (6.1) and integrating the obtained inequality over $(0, T)$, we get

$$\begin{aligned} \mathcal{E}_{p_\beta=p_\alpha}(T) - \mathcal{E}_{p_\beta=p_\alpha}(0) &+ \|\nabla u_a\|_{L^2(\Omega_T)}^2 + \|\nabla u_b\|_{L^2(\Omega_T)}^2 + \|u_a\|_{L^3(\Omega_T)}^3 + \|u_b\|_{L^3(\Omega_T)}^3 \\ &\lesssim \|\nabla v\|_{L^2(\Omega_T)}^2 + \|\Delta v\|_{L^2(\Omega_T)}^2 + \|\theta(u_a, v)\|_{L^2(\Omega_T)}^2 + \|\omega(u_b, v)\|_{L^2(\Omega_T)}^2 \\ &+ \|u_a\|_{L^{p_\alpha}(\Omega_T)}^{p_\alpha} + \|u_b\|_{L^{p_\alpha}(\Omega_T)}^{p_\alpha} + \|v(\theta(u_a, v) + \omega(u_b, v))\|_{L^1(\Omega_T)} \end{aligned} \quad (6.19)$$

Then, we see that all the terms in the right hand side of (6.19) are controlled by the estimates obtained in Theorem 2.1, and we get

$$\begin{aligned} \|u_a\|_{L^\infty(0,T;L^2(\Omega))}^2 + \|u_b\|_{L^\infty(0,T;L^2(\Omega))}^2 + \|\nabla u_a\|_{L^2(\Omega_T)}^2 + \|\nabla u_b\|_{L^2(\Omega_T)}^2 \\ + \|u_a\|_{L^3(\Omega_T)}^3 + \|u_b\|_{L^3(\Omega_T)}^3 \lesssim \mathcal{E}_{p_\beta=p_\alpha}(0) + C(|\Omega|, T). \end{aligned} \quad (6.20)$$
## 6.4 Estimates from $\mathcal{E}_p$, $p \in [2, p_{\alpha,\beta}^0)$

To begin with, we fix $p$ in $[2, p_{\alpha,\beta}^0)$, where $p_{\alpha,\beta}^0$ is defined below and $p_{\alpha,\beta}^0 > 2$ by (H2)

$$ p_{\alpha, \beta}^0 := \begin{cases} 1 + \frac{4}{\beta - 3\alpha - 2}, & \text{if } 2(\alpha + 1) < \beta - \alpha < 2(\alpha + 3), \\ +\infty & \text{if } 0 \le \beta - \alpha \le 2(\alpha + 1). \end{cases} \qquad (6.21) $$

Observing that (because of definition (5.2)) the terms $K_5^p$ and $K_6^p$ in (5.3) both vanish for $p \ge 2$, (5.3) reads as $I_{\text{diff}}^p = K_1^p + K_2^p + K_3^p + K_4^p$. Moreover, the term $K_3^p$ in (5.10) (respectively $K_4^p$ in (5.11)) gives a positive contribution to $I_{\text{diff}}^p$ if and only if $p \in (p_{\alpha}, 1 + \frac{1}{\alpha})$ (respectively $p \in (p_{\beta}, 1 + \frac{1}{\beta})$). If $\alpha < 1$ (respectively $\beta < 1$), there are $p \in [2, 1 + \frac{1}{\alpha})$, (respectively $p \in [2, 1 + \frac{1}{\beta})$), and in that case we use the decreasing character of the constant $C(\alpha, p)$ in (5.10), defined in (3.18), to obtain

$$ 0 < C(\alpha, p) \le C(\alpha, 2) = (1 - \alpha)/\alpha. $$

Hence, the upper bound (5.10) of $K_3^p$ can be absorbed by the upper bound (5.8) of $K_1^p$ (respectively $0 < C(\beta, p) \le (1 - \beta)/\beta$ and the upper bound (5.11) of $K_4^p$ can be absorbed by the upper bound (5.9) of $K_2^p$). Therefore, by (5.8),(5.9), it holds

$$ I_{\text{diff}}^p \lesssim -\frac{4(p-1)}{q(p)^2} a_{\text{fast}}^{\alpha(p-1)} \|\nabla u_a^{\frac{q(p)}{2}}\|_{L^2(\Omega)}^2 + (p-1) a_{\text{fast}}^{1-p} \|\theta(u_a, v)^{\frac{q(p)}{2}-1} \nabla v\|_{L^2(\Omega)}^2 \\ - \frac{4(p-1)}{r(p)^2} b_{\text{fast}}^{\beta(p-1)} \|\nabla u_b^{\frac{r(p)}{2}}\|_{L^2(\Omega)}^2 + (p-1) b_{\text{fast}}^{1-p} \|\omega(u_b, v)^{\frac{r(p)}{2}-1} \nabla v\|_{L^2(\Omega)}^2. \quad (6.22) $$

Plugging (6.22) and the estimate (4.9) of $I_{\text{rea}}^p$ into (6.1), rearranging the terms and integrating the obtained inequality over $(0, T)$, we end up with

$$ \begin{aligned} \mathcal{E}_p(T) - \mathcal{E}_p(0) &+ \frac{4(p-1)}{q(p)^2} a_{\text{fast}}^{\alpha(p-1)} \|\nabla u_a^{\frac{q(p)}{2}}\|_{L^2(\Omega_T)}^2 + \frac{4(p-1)}{r(p)^2} b_{\text{fast}}^{\beta(p-1)} \|\nabla u_b^{\frac{r(p)}{2}}\|_{L^2(\Omega_T)}^2 \\ &+ a_{\text{fast}}^{\alpha(p-1)+1} \frac{A}{1+A} \|u_a\|_{L^{q(p)+1}(\Omega_T)}^{q(p)+1} + b_{\text{fast}}^{\beta(p-1)+1} \frac{B \vee \sigma}{1+B \vee \sigma} \|u_b\|_{L^{r(p)+1}(\Omega_T)}^{r(p)+1} \\ &- \int_0^T I_{\text{fast}}^p(t) dt \\ &\lesssim (p-1) a_{\text{fast}}^{1-p} \|\theta(u_a, v)^{\frac{q(p)}{2}-1} \nabla v\|_{L^2(\Omega_T)}^2 \\ &+ (p-1) b_{\text{fast}}^{1-p} \|\omega(u_b, v)^{\frac{r(p)}{2}-1} \nabla v\|_{L^2(\Omega_T)}^2 \\ &+ (1+A)^{2\alpha(p-1)+1} \|u_a\|_{L^p(\Omega_T)}^p + (1+B \vee \sigma)^{2\beta(p-1)+1} \|u_b\|_{L^p(\Omega_T)}^p \\ &+ a_{\text{fast}}^{-p} \|v \theta(u_a, v)^{q(p)-1}\|_{L^1(\Omega_T)} + b_{\text{fast}}^{-p} \|v \omega(u_b, v)^{r(p)-1}\|_{L^1(\Omega_T)} \end{aligned} $$
Next, we observe that, by (3.7) and the above inequality, for all $p \ge 2$, there exists $C_p(\alpha, \beta, A, B, a_{\text{fast}}, b_{\text{fast}}) > 0$ such that

$$
\begin{align}
& \|u_a^{\varepsilon}\|_{L^{\infty}(0,T;L^{q(p)}(\Omega))}^{q(p)} + \|u_b^{\varepsilon}\|_{L^{\infty}(0,T;L^{r(p)}(\Omega))}^{r(p)} + \|\nabla u_a^{\frac{q(p)}{2}}\|_{L^2(\Omega_T)}^{2} + \|\nabla u_b^{\frac{r(p)}{2}}\|_{L^2(\Omega_T)}^{2} \nonumber \\
& \qquad + \|u_a\|_{L^{q(p)+1}(\Omega_T)}^{q(p)+1} + \|u_b\|_{L^{r(p)+1}(\Omega_T)}^{r(p)+1} - \int_0^T I_{\text{fast}}^p(t) dt \nonumber \\
& \le C_p \|\theta(u_a, v)^{\frac{q(p)}{2}-1} \nabla v\|_{L^2(\Omega_T)}^2 + C_p \|\omega(u_b, v)^{\frac{r(p)}{2}-1} \nabla v\|_{L^2(\Omega_T)}^2 \tag{6.23} \\
& \qquad + C_p \|u_a\|_{L^p(\Omega_T)}^p + C_p \|u_b\|_{L^p(\Omega_T)}^p \nonumber \\
& \qquad + C_p \|v \theta(u_a, v)^{q(p)-1}\|_{L^1(\Omega_T)} + C_p \|v \omega(u_b, v)^{r(p)-1}\|_{L^1(\Omega_T)} + C_p \mathcal{E}_p(0) \nonumber \\
& := Z_1^p + Z_2^p + Z_3^p + Z_4^p + Z_5^p + Z_6^p + C_p \mathcal{E}_p(0). \nonumber
\end{align}
$$

We proceed estimating $Z_1^p, Z_2^p$ and $Z_5^p + Z_6^p$ in such a way to be all absorbed by $\|u_a\|_{L^{q(p)+1}(\Omega_T)}^{q(p)+1}$ and $\|u_b\|_{L^{r(p)+1}(\Omega_T)}^{r(p)+1}$ in the left hand side of (6.23). As the computations will be in the line of the computations carried out in the case $p = p_\alpha$ with $\alpha < \beta$ (see Subsection 6.2), redundant details will be omitted. Hereafter, the constant $C_p$ will change from line to line and it may depends also on $|\Omega|, T$ and the constants $C_p^{\text{MR}}, C_p^{\text{ADN}}, C_p^{\text{GN}}$ in (A.13),(A.14),(A.15) respectively.

**Estimate of $Z_1^p$.** Let $\delta > 0$ to be chosen later and note that, if $p \ge 2$, $q(p)-2 \ge \alpha > 0$. Applying Young's inequality, we have

$$
\begin{align}
Z_1^p = C_p \int_{\Omega_T} \theta(u_a, v)^{q(p)-2} |\nabla v|^2 dxdt \tag{6.24} \\
\le \delta C_p \int_{\Omega_T} |\nabla v|^{2(q(p)+1)} dxdt + \delta^{-\frac{1}{q(p)}} C_p \int_{\Omega_T} \theta(u_a, v)^{(q(p)-2)\frac{q(p)+1}{q(p)}} dxdt. \nonumber
\end{align}
$$

By (A.18) and $q(p) \le r(p)$, we obtain

$$
\delta \int_{\Omega_T} |\nabla v|^{2(q(p)+1)} dx dt \le \delta C_p \left( 1 + \|u_a\|_{L^{q(p)+1}(\Omega_T)}^{q(p)+1} + \|u_b\|_{L^{r(p)+1}(\Omega_T)}^{r(p)+1} \right). \quad (6.25)
$$

Next, using (4.1) with $\gamma = (q(p) - 2) \frac{q(p)+1}{q(p)}$, $\gamma_1 = q(p) + 1$ and $C = \delta^{-\frac{1}{2}}$ so that
$\delta^{-\frac{1}{q(p)}} C^{\gamma-\gamma_1} = \delta$, we get

$$
\begin{align}
\delta^{-\frac{1}{q(p)}} \int_{\Omega_T} \theta(u_a, v)^{(q(p)-2)\frac{q(p)+1}{q(p)}} dx dt &\le \delta^{-\frac{q(p)-1}{2}} |\Omega_T| + \delta \int_{\Omega_T} \theta(u_a, v)^{q(p)+1} dx dt \nonumber \\
&\le \delta^{-\frac{q(p)-1}{2}} |\Omega_T| + \delta C_p \left(1 + \|u_a\|_{L^{q(p)+1}(\Omega_T)}^{q(p)+1}\right). \tag{6.26}
\end{align}
$$

Plugging (6.25),(6.26) into (6.24), we end up with

$$
Z_1^p \le C_p \delta \left(1 + \delta^{-\frac{q(p)+1}{2}}\right) + \delta C_p \left( \|u_a\|_{L^{q(p)+1}(\Omega_T)}^{q(p)+1} + \|u_b\|_{L^{r(p)+1}(\Omega_T)}^{r(p)+1} \right). \quad (6.27)
$$

**Estimate of $Z_2^p$.** Let us observe that $r(p) - 2 \ge \beta > 0$, if $p \ge 2$. Hence, proceeding
as above,

$$
\begin{align}
Z_2^p &= C_p \int_{\Omega_T} \omega(u_b, v)^{r(p)-2} |\nabla v|^2 dxdt \nonumber \\
&\le \delta C_p \int_{\Omega_T} |\nabla v|^{2(q(p)+1)} dxdt + \delta^{-\frac{1}{q(p)}} C_p \int_{\Omega_T} \omega(u_b, v)^{(r(p)-2)\frac{q(p)+1}{q(p)}} dxdt \nonumber \\
&\le \delta C_p \left( 1 + \|u_a\|_{L^{q(p)+1}(\Omega_T)}^{q(p)+1} + \|u_b\|_{L^{r(p)+1}(\Omega_T)}^{r(p)+1} \right) \nonumber \\
&\quad + \delta^{-\frac{1}{q(p)}} C_p \int_{\Omega_T} \omega(u_b, v)^{(r(p)-2)\frac{q(p)+1}{q(p)}} dxdt. \tag{6.28}
\end{align}
$$

If $\alpha = \beta$, then $r(p) = q(p)$ and we can proceed as in (6.26) to obtain from (6.28)

$$
Z_2^p \le \delta C_p \left(1 + \delta^{-\frac{r(p)+1}{2}}\right) + \delta C_p \left( \|u_a\|_{L^{q(p)+1}(\Omega_T)}^{q(p)+1} + \|u_b\|_{L^{r(p)+1}(\Omega_T)}^{r(p)+1} \right). \quad (6.29)
$$

If $\alpha < \beta$, using definitions (2.10), as $p \in [2, p_{\alpha,\beta}^0)$, it holds

$$
0 < (r(p) - 2) \frac{q(p) + 1}{q(p)} < r(p) + 1. \tag{6.30}
$$

Therefore, we apply once again (4.1) with $\gamma = (r(p) - 2) \frac{q(p)+1}{q(p)}$, $\gamma_1 = r(p) + 1$ and
$C = C(\delta) > 0$ such that $\delta^{-\frac{1}{q(p)}} C(\delta)^{\gamma-\gamma_1} = \delta$, to get

$$
\begin{align}
\delta^{-\frac{1}{q(p)}} \int_{\Omega_T} \omega(u_b, v)^{(r(p)-2)\frac{q(p)+1}{q(p)}} dxdt &\le \delta C(\delta)^{r(p)+1} |\Omega_T| + \delta \int_{\Omega_T} \omega(u_b, v)^{r(p)+1} dxdt \nonumber \\
&\le \delta C(\delta)^{r(p)+1} |\Omega_T| + \delta C_p \left(1 + \|u_b\|_{L^{r(p)+1}(\Omega_T)}^{r(p)+1}\right). \tag{6.31}
\end{align}
$$

Plugging (6.31) into (6.28) and taking into account that $C(\delta)^{r(p)+1} = \delta^{-\frac{(q(p)+1)(r(p)+1)}{2+3q(p)-r(p)}}$,
we obtain

$$
Z_2^p \le \delta C_p \left(1 + \delta^{-\frac{(q(p)+1)(r(p)+1)}{2+3q(p)-r(p)}}\right) + \delta C_p \left( \|u_a\|_{L^{q(p)+1}(\Omega_T)}^{q(p)+1} + \|u_b\|_{L^{r(p)+1}(\Omega_T)}^{r(p)+1} \right). \quad (6.32)
$$

Note that (6.32) becomes (6.29) when $r(p) = q(p)$ (i.e. $\alpha = \beta$).

Estimate of $Z_5^p + Z_6^p$. Proceeding as in (6.13), (6.14), we have

$$
\begin{align}
Z_5^p + Z_6^p &= C_p \int_{\Omega_T} v \theta(u_a, v)^{q(p)-1} dx dt + C_p \int_{\Omega_T} v \omega(u_b, v)^{r(p)-1} dx dt \nonumber \\
&\le C_p \int_{\Omega_T} u_a^{q(p)-1} dx dt + C_p \int_{\Omega_T} u_b^{r(p)-1} dx dt + C_p |\Omega_T| \tag{6.33} \\
&\le C_p \delta \left( \delta^{-\frac{q(p)+1}{2}} + \delta^{-\frac{r(p)+1}{2}} \right) \nonumber \\
&\quad + \delta C_p \left( \|u_a\|_{L^{q(p)+1}(\Omega_T)}^{q(p)+1} + \|u_b\|_{L^{r(p)+1}(\Omega_T)}^{r(p)+1} \right) + C_p |\Omega_T|. \nonumber
\end{align}
$$

Final energy estimate. Plugg (6.27), (6.32), (6.33) into (6.23) and, for $\delta' \in (0, 1)$,
set $\delta$ so that $3\delta C_p = 1 - \delta'$. Then, $\delta^{-1} \lesssim C_p$ and it holds

$$
\begin{align}
& \|u_a\|_{L^\infty(0,T;L^{q(p)}(\Omega))}^{q(p)} + \|u_b\|_{L^\infty(0,T;L^{r(p)}(\Omega))}^{r(p)} + \|\nabla u_a^{q(p)/2}\|_{L^2(\Omega_T)}^2 + \|\nabla u_b^{r(p)/2}\|_{L^2(\Omega_T)}^2 \nonumber \\
& \qquad + \delta' \left( \|u_a\|_{L^{q(p)+1}(\Omega_T)}^{q(p)+1} + \|u_b\|_{L^{r(p)+1}(\Omega_T)}^{r(p)+1} \right) - \int_0^T I_{\text{fast}}^p(t) dt \nonumber \\
& \le C_p (\|u_a\|_{L^p(\Omega_T)}^p + \|u_b\|_{L^p(\Omega_T)}^p + \mathcal{E}_p(0) + |\Omega_T|) \tag{6.34} \\
& \qquad + (1 + C_p^{\frac{q(p)+1}{2}} + C_p^{\frac{r(p)+1}{2}} + C_p^{\frac{(q(p)+1)(r(p)+1)}{2+3q(p)-r(p)}}). \nonumber
\end{align}
$$

6.5 Bootstrapping and end of the proof

We are now ready to prove Lemma 2.2.

First, recall that $q(p_\alpha) = r(p_\beta) = 2$. So, if $\alpha < \beta$, estimates (6.4) and (6.17)
imply (2.15) by interpolation. If $\alpha = \beta$, (6.20) is exactly (2.15).

Next, we have proved that, if $p \in [2, p_{\alpha,\beta}^0)$ (see (6.21)), (6.34) holds true and gives a bound on $\|u_a\|_{L^{q(p)}(\Omega_T)}$ and $\|u_b\|_{L^{r(p)}(\Omega_T)}$ as soon as we have a bound on the $L^p(\Omega_T)$ norm of $u_a$ and $u_b$. Hence, recalling that $q(p) \le r(p)$ and starting from exponent 2 $= q(p_\alpha)$ and the $L^2(\Omega_T)$ estimates of $u_a$ and $u_b$ in (2.7), we can bootstrap an $L^{q(p)}(\Omega_T)$ bound of $u_a$ and $u_b$ to an $L^{q(q(p))}(\Omega_T)$ bound of $u_a$ and $u_b$ until $q(p) < p_{\alpha,\beta}^0$. The two cases below have to be considered.

(i) $\beta - \alpha \in [0, 2(\alpha + 1)]$. Then, $p_{\alpha,\beta}^0 = +\infty$ and estimate (6.34) implies (2.16).

(ii) $\beta - \alpha \in (2(\alpha + 1), 2(\alpha + 3)) = \cup_n I_n$, (see (2.13)). Then, $p_{\alpha,\beta}^0 \in (2, +\infty)$. In order to set up the bootstrap procedure, we denote

$$
p_{\alpha, \beta}^{n} := 1 + \frac{4}{(\alpha + 1)^n (\beta - 3\alpha - 2)}, \quad n \ge 1,
$$

and we observe that, by (2.10), it holds

$$
p_{\alpha, \beta}^{0} = (\overbrace{q \circ \cdots \circ q}^{n \text{-times}})(p_{\alpha, \beta}^{n}) := q_{n}(p_{\alpha, \beta}^{n}), \quad n \geq 0.
$$

Let $n_{\alpha,\beta} \ge 0$ be the largest integer such that $\beta - \alpha \in I_{n_{\alpha,\beta}}$. We have that
$2 < p_{\alpha,\beta}^{n_{\alpha,\beta}}$ and $q_{n_{\alpha,\beta}}(2) < q_{n_{\alpha,\beta}}(p_{\alpha,\beta}^{n_{\alpha,\beta}}) = p_{\alpha,\beta}^{0}$. Therefore, starting from (2.7),
we can bootstrap (6.34) till $q_{n_{\alpha,\beta}}(2) = (\alpha+1)^{n_{\alpha,\beta}} + 1$ to get estimate (2.17).

Finally, taking $p=2$ in (6.34), by (3.10) we obtain (2.18). The maximal regularity estimate (2.19) and the gradient estimate (2.20) follow by (A.17) and (A.18) respectively, using (2.16) or (2.17), according to the value of $\beta - \alpha$. For that, it is worth recall that when $\beta - \alpha \in (2(\alpha + 1), 2(\alpha + 3))$ and $p = (\alpha + 1)^{n_{\alpha,\beta}} + 1$, then $q(p) + 1 = 2 + (\alpha + 1)^{n_{\alpha,\beta}+1}$.

7 Existence for the cross-diffusion system

This section is devoted to the proofs of Lemma 2.4 and Theorem 2.5. The latter
follows by compactness arguments based on the previous estimates on the unique
global strict solution $(u_a^\epsilon, u_b^\epsilon, v^\epsilon)$ of (1.1)-(1.7). The key point is the identification
of the limit (as $\epsilon \to 0$) of the densities pair $(u_a^\epsilon, u_b^\epsilon)$ with the unique solution of
the nonlinear system (1.12) corresponding to the limit of the densities pair $(u^\epsilon =
u_a^\epsilon + u_b^\epsilon, v^\epsilon)$. This is the object of Lemma 2.4.
## 7.1 Proof of Lemma 2.4

By estimates (2.18), (2.7) there exists $C(T) > 0$ such that, for all $\varepsilon > 0$,

$$ \|(\Lambda^{1/2}Q)(u_a^\varepsilon, u_b^\varepsilon, v^\varepsilon)\|_{L^2(\Omega_T)} \le \sqrt{\varepsilon} C(T). $$

As $\Lambda(u_a^\varepsilon, u_b^\varepsilon, v^\varepsilon) \ge A^\alpha > 0$, the latter implies

$$ \|Q(u_a^\varepsilon, u_b^\varepsilon, v^\varepsilon)\|_{L^2(\Omega_T)} \le A^{-\alpha/2} \sqrt{\varepsilon} C(T). \quad (7.1) $$

Now, following [3], let us define $Q(u, u_b, v) := Q(u - u_b, u_b, v)$, for $(u, u_b, v) \in \mathbb{R}_+^3$ such that $u_b \le u$. Note that, by (1.5), (1.7), for all $(u_a, u_b, v) \in \mathbb{R}_+^3$,

$$ \begin{aligned} \partial_1 Q(u_a, u_b, v) = & -\psi(a_{\text{fast}} u_a + c_{\text{fast}} v) / \Lambda(u_a, u_b, v) \\ & - a_{\text{fast}}(u_a + u_b)\phi(b_{\text{fast}} u_b + d_{\text{fast}} v)\psi'(a_{\text{fast}} u_a + c_{\text{fast}} v) / \Lambda^2(u_a, u_b, v) \end{aligned} \quad (7.2) $$

$$ \begin{aligned} \partial_2 Q(u_a, u_b, v) = & \phi(b_{\text{fast}} u_b + d_{\text{fast}} v) / \Lambda(u_a, u_b, v) \\ & + b_{\text{fast}}(u_a + u_b)\phi'(b_{\text{fast}} u_b + d_{\text{fast}} v)\psi'(a_{\text{fast}} u_a + c_{\text{fast}} v) / \Lambda^2(u_a, u_b, v). \end{aligned} \quad (7.3) $$

Hence,

$$ \partial_2 Q(u, u_b, v) = -\partial_1 Q(u - u_b, u_b, v) + \partial_2 Q(u - u_b, u_b, v) \ge 1. \quad (7.4) $$

Recalling that $(u_a^\varepsilon(u^\varepsilon, v^\varepsilon), u_b^\varepsilon(u^\varepsilon, v^\varepsilon))$ is the unique solution of system (1.12) corresponding to $(u^\varepsilon, v^\varepsilon)$, it holds $Q(u^\varepsilon, u_b^\varepsilon(u^\varepsilon, v^\varepsilon), v^\varepsilon) = 0$. Therefore, for some intermediate value $\xi \in [0, u^\varepsilon]$ between $u_b^\varepsilon \in [0, u^\varepsilon]$ and $u_b^\varepsilon(u^\varepsilon, v^\varepsilon) \in [0, u^\varepsilon]$, it follows

$$ Q(u_a^\varepsilon, u_b^\varepsilon, v^\varepsilon) = Q(u^\varepsilon, u_b^\varepsilon, v^\varepsilon) - Q(u^\varepsilon, u_b^\varepsilon(u^\varepsilon, v^\varepsilon), v^\varepsilon) = \partial_2 Q(u^\varepsilon, \xi, v^\varepsilon)(u_b^\varepsilon - u_b^\varepsilon(u^\varepsilon, v^\varepsilon)). $$

Using (7.4) and the previous identity, we have

$$ |u_b^\varepsilon - u_b^\varepsilon(u^\varepsilon, v^\varepsilon)| \le \partial_2 Q(u^\varepsilon, \xi, v^\varepsilon) |u_b^\varepsilon - u_b^\varepsilon(u^\varepsilon, v^\varepsilon)| = |Q(u_a^\varepsilon, u_b^\varepsilon, v^\varepsilon)|, \quad (7.5) $$

and (2.23) follows by (7.1). The Lemma is proved.

## 7.2 Proof of Theorem 2.5

In the sequel $p$ will satisfies (2.24).

**Convergence of** $(v_\varepsilon)_\varepsilon$. Thanks to estimates (2.6), (2.19), (2.20), for all $T > 0$, the sequence $(v^\varepsilon)_\varepsilon$ is bounded in $W^{1,p}(\Omega_T) \cap L^p(0, T; W^{2,p}(\Omega)) \cap L^\infty((0, \infty) \times \Omega)$. By (2.16), (2.17), the sequence $(f_v(u_a^\varepsilon, u_b^\varepsilon, v^\varepsilon))_\varepsilon$ is bounded in $L^p(\Omega_T)$, (see also (1.4)). Owing to the Aubin-Lions's Lemma [1] and standard weak compactness arguments, it follows that, for any $T > 0$, there exists a subsequence of $(v^\varepsilon)_\varepsilon$ (still denoted $v^\varepsilon$) and $v \in W^{1,p}(\Omega_T) \cap L^p(0, T; W^{2,p}(\Omega)) \cap L^\infty((0, T) \times \Omega)$ such that,

$$ \begin{aligned} v^\varepsilon \to v, & \quad \text{in } L^p(0, T; W^{1,p}(\Omega)) \text{ and a.e. in } \Omega_T, \quad \text{as } \varepsilon \to 0, \\ v^\varepsilon \to v & \quad \text{in } W^{1,p}(\Omega_T) \text{ and in } L^p(0, T; W^{2,p}(\Omega)), \quad \text{as } \varepsilon \to 0. \end{aligned} \quad (7.6) $$

Moreover, $v$ is nonnegative as $v^\varepsilon$ is nonnegative.

**Convergence of** $(u^\varepsilon)_\varepsilon = (u_a^\varepsilon + u_b^\varepsilon)_\varepsilon$. By estimates (2.15)-(2.17), (see also Remark 1), for all $T > 0$, the sequences $(u_a^\varepsilon)_\varepsilon$, $(u_b^\varepsilon)_\varepsilon$ are bounded in $L^2(0, T; H^1(\Omega)) \cap$
$L^p(\Omega_T) \cap L^\infty(0, T; L^{p-1}(\Omega))$. Hence, $(f_u(u_a^\varepsilon, u_b^\varepsilon, v^\varepsilon))_\varepsilon$ is bounded in $L^\frac{p}{2}(\Omega_T)$ (because of the quadratic terms, see (1.9), (1.4)). We denote $s = \frac{p}{2} \wedge 2$ and we continue according the value of $\beta - \alpha$.

If $\beta - \alpha \in [0, 2(\alpha + 1)]$, we can choose $p = 4$ (see (2.24)), which yields $s = 2$ and $(f_u(u_a^\varepsilon, u_b^\varepsilon, v^\varepsilon))_\varepsilon$ is bounded in $L^2(\Omega_T)$. Therefore, the equation satisfied by $u^\varepsilon$ (see (1.8)) implies that $(\partial_t u^\varepsilon)_\varepsilon$ is bounded in $L^2(0, T; H^{-1}(\Omega))$.

If $\beta - \alpha \in (2(\alpha + 1), 2(\alpha + 3))$, then $p > 3$ (see (2.24)), $s \in (\frac{3}{2}, 2]$ and we have to argue according the space dimension. If $N \le 6$, by Sobolev's embedding theorem $L^s(\Omega) \subset H^{-1}(\Omega)$ with continuous embedding, so that $(\partial_t u^\varepsilon)_\varepsilon$ is bounded in $L^s(0, T; H^{-1}(\Omega))$. If $N > 6$, by the assumption $(\alpha + 1)^{n_{\alpha,\beta}+1} \ge 2$, we have that $p \ge 4$, which yields $s = 2$ and again $(\partial_t u^\varepsilon)_\varepsilon$ is bounded in $L^2(0, T; H^{-1}(\Omega))$.

As above, for any $T > 0$, it holds the existence of a subsequence of $(u^\varepsilon)_\varepsilon$ (still denoted $u^\varepsilon$) and $u \in L^2(0, T; H^1(\Omega)) \cap L^p(\Omega_T) \cap L^\infty(0, T; L^{p-1}(\Omega))$, such that

$$ 
\begin{align} 
u^\varepsilon &= u_a^\varepsilon + u_b^\varepsilon \to u, \quad \text{in } L^2(\Omega_T) \text{ and a.e. in } \Omega_T, \quad \text{as } \varepsilon \to 0, \nonumber \\ 
u^\varepsilon &\to u \quad \text{in } L^2(0, T; H^1(\Omega)) \text{ and in } L^p(\Omega_T), \quad \text{as } \varepsilon \to 0, \tag{7.7} \\ 
\partial_t u^\varepsilon &\to \partial_t u \quad \text{in } L^s(0, T; H^{-1}(\Omega)), \quad \text{as } \varepsilon \to 0. \nonumber 
\end{align} 
$$

Furthermore, $u$ is nonnegative as $u^\varepsilon$ is positive.

**Convergence of** $(u_a^\varepsilon)_\varepsilon, (u_b^\varepsilon)_\varepsilon$. By the boundedness of the sequences $(u_a^\varepsilon)_\varepsilon, (u_b^\varepsilon)_\varepsilon$ quoted above, for any $T > 0$, there exists $u_a, u_b \in L^2(0, T; H^1(\Omega)) \cap L^p(\Omega_T) \cap L^\infty(0, T; L^{p-1}(\Omega))$ such that, for subsequences (still denoted $u_a^\varepsilon$ and $u_b^\varepsilon$)

$$ u_a^\varepsilon \to u_a, \quad u_b^\varepsilon \to u_b \quad \text{in } L^2(0, T; H^1(\Omega)) \text{ and in } L^p(\Omega_T), \quad \text{as } \varepsilon \to 0. \quad (7.8) $$

On the other hand, using (2.23), there exists subsequences such that

$$ Q(u_a^\varepsilon, u_b^\varepsilon, v^\varepsilon) \to 0, \quad |u_b^\varepsilon - u_b^*(u^\varepsilon, v^\varepsilon)| = |u_a^\varepsilon - u_a^*(u^\varepsilon, v^\varepsilon)| \to 0, \quad \text{a.e. in } \Omega_T. $$

By the continuity of the map $(u_a^*, u_b^*)$ (see Lemma 2.3) and the a.e. convergences of $(u^\varepsilon, v^\varepsilon)$ towards $(u, v)$ obtained above, it follows that

$$ (u_a^\varepsilon, u_b^\varepsilon) \to (u_a^*(u, v), u_b^*(u, v)), \quad \text{a.e. in } \Omega_T, \quad \text{as } \varepsilon \to 0, \quad (7.9) $$

Therefore, $(u_a^*(u, v), u_b^*(u, v)) = (u_a, u_b)$ a.e. in $\Omega_T$, and

$$ u_a^*(u, v), u_b^*(u, v) \in L^2(0, T; H^1(\Omega)) \cap L^p(\Omega_T) \cap L^\infty(0, T; L^{p-1}(\Omega)). $$

**Diagonal extraction.** Since $T > 0$ is arbitrarily large, we can apply the diagonal extraction argument. It follows that there exists a subsequence $(\varepsilon_k)_k$ and a pair of nonnegative measurable functions $(u, v) : (0, \infty) \times \Omega \to \mathbb{R}_+^2$ satisfying (i), (ii) in the statement of Theorem 2.5 and

$$ u^{\varepsilon k} \to u, \quad v^{\varepsilon k} \to v, \quad (u_a^{\varepsilon k}, u_b^{\varepsilon k}) \to (u_a^*(u, v), u_b^*(u, v)), \quad \text{a.e. in } (0, \infty) \times \Omega, \quad \text{as } k \to \infty, \quad (7.10) $$

and also the convergence in (7.6), (7.7), (7.8), thus giving (iii), (iv).

**Conclusion.** It remains to show that $(u, v)$ is a global (weak, strong) solution according to (2.25). To begin with, we consider the weak formulation of the equation for $u^\varepsilon$ in (1.8) with test functions as in (2.25), i.e.

$$ 
\begin{align} 
& - \int_0^T \int_\Omega u^\varepsilon \partial_t w \, dxdt - \int_\Omega u^{\text{in}} w(0) \, dx + \int_0^T \int_\Omega \nabla(d_a u_a^\varepsilon + d_b u_b^\varepsilon) \cdot \nabla w \, dxdt \nonumber \\ 
& = \int_0^T \int_\Omega f_u(u_a^\varepsilon, u_b^\varepsilon, v^\varepsilon) w \, dxdt. \tag{7.11} 
\end{align} 
$$
It is worth noticing that the term in the right hand side of (7.11) is well defined. Indeed, this is clear if $\beta - \alpha \in [0, 2(\alpha + 1)]$ or if $\beta - \alpha \in (2(\alpha + 1), 2(\alpha + 3))$ and $N > 6$, since then $f_u(u_a^\varepsilon, u_b^\varepsilon, v^\varepsilon) \in L^2(\Omega_T)$, as observed above. On the other hand, if $\beta - \alpha \in (2(\alpha + 1), 2(\alpha + 3))$ and $2 < N \le 6$, the Sobolev's embedding $H^1(\Omega) \subset L^{\frac{2N}{N-2}}(\Omega)$ implies that the quadratic terms in $f_u(u_a^\varepsilon, u_b^\varepsilon, v^\varepsilon)$ belong to $L^1(0, T; L^{\frac{2N}{N-2}}(\Omega))$, while $w \in L^\infty(0, T; L^{\frac{2N}{N-2}}(\Omega))$, and we are able to conclude. If $N = 1, 2$, similar arguments give us the claim.

Hence, thanks to the above convergence properties of $u_a^{\varepsilon^k}$, $u_b^{\varepsilon^k}$, $v^{\varepsilon^k}$, the boundedness of $(f_u(u_a^\varepsilon, u_b^\varepsilon, v^\varepsilon))_\varepsilon$ in $L^\frac{p}{2}(\Omega_T)$, for all $T > 0$, the convergence of $(f_u(u_a^{\varepsilon^k}, u_b^{\varepsilon^k}, v^{\varepsilon^k}))_k$ towards $f_u(u_a^*(u, v), u_b^*(u, v), v)$ a.e. in $(0, \infty) \times \Omega$ as $k \to \infty$, we can pass to the limit $k \to \infty$ in (7.11) and the equation for $u$ in (2.25) holds true.

Furthermore, $u \in W^{1,s}(0, T; H^{-1}(\Omega))$ and $W^{1,s}(0, T; H^{-1}(\Omega))$ is continuously embedded in $C^0([0, T]; H^{-1}(\Omega))$. Hence, the operator $w \to w(0)$ is weakly sequentially continuous from $W^{1,s}(0, T; H^{-1}(\Omega))$ weak to $H^{-1}(\Omega)$ weak. As, for all $T > 0$, $u^{\varepsilon^k} \to u$ in $W^{1,s}(0, T; H^{-1}(\Omega))$ and $u^{\varepsilon^k}(0) = u^{\text{in}}$, we have $u(0) = u^{\text{in}}$ a.e. in $\Omega$.

Next, concerning the $v$ component of the solution, recall that $(f_v(u_a^\varepsilon, u_b^\varepsilon, v^\varepsilon))_\varepsilon$ is bounded in $L^p(\Omega_T)$, for all $T > 0$, and that $f_v(u_a^{\varepsilon^k}, u_b^{\varepsilon^k}, v^{\varepsilon^k}) \to f_v(u_a^*(u, v), u_b^*(u, v), v)$ a.e. in $(0, \infty) \times \Omega$ as $k \to \infty$. Hence, the previous convergence properties of $(v^{\varepsilon^k})_k$, applied to the equation

$$ \partial_t v^{\varepsilon^k} - d_v \Delta v^{\varepsilon^k} = f_v(u_a^{\varepsilon^k}, u_b^{\varepsilon^k}, v^{\varepsilon^k}) $$

yield that

$$ \partial_t v = d_v \Delta v + f_v(u_a^*(u, v), u_b^*(u, v), v), \quad \text{in } \mathcal{D}'((0, \infty) \times \Omega), $$

and then that $v$ satisfies the above equation in $L^p(\Omega_T)$, for all $T > 0$, by the $L^p(\Omega_T)$-integrability of each term in the equation. As $W^{1,p}(\Omega_T) \subset W^{1,p}(0, T; L^p(\Omega)) \subset C^0([0, T]; L^p(\Omega))$ with continuous embeddings, the operator $w \to w(0)$ is weakly sequentially continuous from $W^{1,p}(\Omega_T)$ weak to $L^p(\Omega)$ weak. As, for all $T > 0$, $v^{\varepsilon^k} \to v$ in $W^{1,p}(\Omega_T)$ and $v^{\varepsilon^k}(0) = v^{\text{in}}$, we have $v(0) = v^{\text{in}}$ a.e. in $\Omega$. Finally, the regularity of $v$ implies that $v$ satisfies homogeneous Neumann boundary condition in the sense of traces. Theorem 2.5 is proved.

# 8 Uniqueness for the cross-diffusion system

This section is devoted to the proof of Theorem 2.6. More precisely, we show that, for all $T > 0$, the function

$$ \phi_\lambda(\tau) := \|u_1 - u_2\|_{L^2(\Omega_\tau)}^2 + \lambda \|v_1 - v_2\|_{L^2(\Omega_\tau)}^2, \quad \tau \in [0, T], \quad (8.1) $$

with $\lambda > 0$ large enough, satisfies an integral inequality that gives (2.26) by Gronwall's lemma. As a consequence, solutions of (2.25), whose components are both bounded, are unique. The key tool are test functions introduced by Oleinik (see [28] and the references therein, see also [21]).

Let us denote

$$ f_u^*(u, v) := f_u(u_a^*(u, v), u_b^*(u, v), v), \quad f_v^*(u, v) := f_v(u_a^*(u, v), u_b^*(u, v), v). $$
The weak formulation of (1.10), (1.11), (1.12) for both components $u, v$, writes as

$$ -\iint_{\Omega_T} u \partial_t w_1 \, dxdt - \int_{\Omega} u^{\text{in}} w_1(0) \, dx + \iint_{\Omega_T} \nabla A(u, v) \cdot \nabla w_1 \, dxdt = \iint_{\Omega_T} f_u^*(u, v) w_1 \, dxdt, \quad (8.2) $$

$$ -\iint_{\Omega_T} v \partial_t w_2 \, dxdt - \int_{\Omega} v^{\text{in}} w_2(0) \, dx + \iint_{\Omega_T} \nabla v \cdot \nabla w_2 \, dxdt = \iint_{\Omega_T} f_v^*(u, v) w_2 \, dxdt, \quad (8.3) $$

for all $T > 0$ and all $w_1, w_2 \in C^1([0, T]; H^1(\Omega))$ such that $w(T) = 0$.

We may assume without loss of generality that $d_b > d_a$, so that the function $A(u, v)$ in (1.11) rewrites as

$$ A(u, v) = d_a u_a^*(u, v) + d_b u_b^*(u, v) = d_a u + (d_b - d_a) u_b^*(u, v), \quad (8.4) $$

and, using (2.21), (2.22), it holds

$$ \partial_1 A(u, v) = d_a + (d_b - d_a) \partial_1 u_b^*(u, v) \in (d_a, d_b) \quad (8.5) $$

$$ \partial_2 A(u, v) = (d_b - d_a) \partial_2 u_b^*(u, v) \in \left( -(d_b - d_a) \frac{d_{\text{fast}}}{b_{\text{fast}}}, (d_b - d_a) \frac{c_{\text{fast}}}{a_{\text{fast}}} \right). \quad (8.6) $$

Let $(u_i, v_i)$, $i = 1, 2$, be two solutions with initial data $(u_i^{\text{in}}, v_i^{\text{in}})$ and let us denote $A_i := A(u_i, v_i)$, $i = 1, 2$. As the functions $A_i$ belong to $L^2(0, T; H^1(\Omega))$, the function

$$ w_1^\tau(t, x) = \begin{cases} \int_t^\tau (A_1(s, x) - A_2(s, x)) \, ds, & \text{if } 0 \le t \le \tau, \\ 0, & \text{if } \tau \le t \le T, \end{cases} \quad (8.7) $$

belongs to $H^1([0, T]; H^1(\Omega))$ and $w_1^\tau(T) = 0$, for all $\tau \in [0, T]$. Therefore, by the density of $C^1([0, T]; H^1(\Omega))$ in $H^1([0, T]; H^1(\Omega))$, we can use $w_1^\tau$ as test functions. Hence, testing the equation satisfied by $u_1 - u_2$ against $w_1^\tau$, we obtain (see (8.2))

$$ \iint_{\Omega_{\tau}} (u_1 - u_2)(t, x)(A_1 - A_2)(t, x) dxdt - \int_{\Omega} (u_1^{\text{in}} - u_2^{\text{in}})(x) \int_{0}^{\tau} (A_1 - A_2)(s, x) ds \, dx \\ + \iint_{\Omega_{\tau}} \nabla(A_1 - A_2)(t, x) \cdot \int_{t}^{\tau} \nabla(A_1 - A_2)(s, x) \, ds \, dxdt \quad (8.8) \\ = \iint_{\Omega_{\tau}} (f_u^*(u_1, v_1) - f_u^*(u_2, v_2))(t, x) \int_{t}^{\tau} (A_1 - A_2)(s, x) ds \, dxdt. $$

It is convenient to split the first term in the left hand side of (8.8) as

$$ \begin{aligned} & \iint_{\Omega_{\tau}} (u_1 - u_2)(t, x)(A_1 - A_2)(t, x) dxdt \\ &= \iint_{\Omega_{\tau}} (u_1 - u_2)(t, x)(A(u_1, v_1) - A(u_2, v_1))(t, x) dxdt \\ & \quad + \iint_{\Omega_{\tau}} (u_1 - u_2)(t, x)(A(u_2, v_1) - A(u_2, v_2))(t, x) dxdt \\ &=: I_1 + I_2. \end{aligned} \quad (8.9) $$
Indeed, by (8.5), A is increasing in u with $\partial_1 A(u, v)$ lower bounded, so that $I_1$ is positive and lower bounded

$$I_1 \geq d_a \|u_1 - u_2\|_{L^2(\Omega_\tau)}^2. \qquad (8.10)$$

On the other hand, denoting $c_2 := (d_b - d_a) (\frac{c_{\text{fast}}}{a_{\text{fast}}} \lor \frac{d_{\text{fast}}}{b_{\text{fast}}})$, by (8.6) we obtain for $I_2$

$$|I_2| \leq \sigma d_a \|u_1 - u_2\|_{L^2(\Omega_\tau)}^2 + \frac{c_2^2}{4\sigma d_a} \|v_1 - v_2\|_{L^2(\Omega_\tau)}^2, \qquad (8.11)$$

with $\sigma > 0$ to be chosen later.

Using (8.5), (8.6) again, the second term in the left hand side of (8.8) can be estimated as following, for all $\tau \in [0, T]$,

$$\begin{aligned} \int_{\Omega} (u_1^{\text{in}} - u_2^{\text{in}})(x) \int_{0}^{\tau} (A_1 - A_2)(s, x) ds dx &\leq \frac{\sigma d_a}{2d_b^2} \|A_1 - A_2\|_{L^2(\Omega_\tau)}^2 \\ &\quad + \frac{T d_b^2}{2\sigma d_a} \|u_1^{\text{in}} - u_2^{\text{in}}\|_{L^2(\Omega)}^2 \\ &\leq \sigma d_a \|u_1 - u_2\|_{L^2(\Omega_\tau)}^2 + \sigma \frac{c_2^2 d_a}{d_b^2} \|v_1 - v_2\|_{L^2(\Omega_\tau)}^2 + \frac{T d_b^2}{2\sigma d_a} \|u_1^{\text{in}} - u_2^{\text{in}}\|_{L^2(\Omega_\tau)}^2. \end{aligned} \qquad (8.12)$$

Denoting $y_i(t, x) = \int_t^\tau \partial_{x_i}(A_1 - A_2)(s, x) ds$, the third term in the left hand side of (8.8) turn out to be positive and it can be neglected, since it writes as

$$\begin{aligned} \iint_{\Omega_\tau} \nabla(A_1 - A_2)(t, x) \cdot \int_t^\tau \nabla(A_1 - A_2)(s, x) ds \,dx dt \\ = -\frac{1}{2} \int_{\Omega} \int_0^\tau \partial_t \left( \sum_i y_i^2(t, x) \right) dt dx \qquad (8.13) \\ = \frac{1}{2} \int_{\Omega} \sum_i \left( \int_0^\tau \partial_{x_i}(A_1 - A_2)(s, x) ds \right)^2 dx. \end{aligned}$$

Finally, for the term in the right hand side of (8.8) we have

$$\begin{aligned} \iint_{\Omega_\tau} (f_u^*(u_1, v_1) - f_u^*(u_2, v_2))(t, x) \int_t^\tau (A_1 - A_2)(s, x) ds \,dx dt \qquad (8.14) \\ \leq \int_0^\tau ds \int_0^s dt \int_\Omega dx \|(f_u^*(u_1, v_1) - f_u^*(u_2, v_2))(t, x)\|(A_1 - A_2)(s, x)| \end{aligned}$$

and

$$\begin{aligned} \int_{\Omega} |(f_u^*(u_1, v_1) - f_u^*(u_2, v_2))(t, x)| |(A_1 - A_2)(s, x)| dx \\ \leq \int_{\Omega} |(f_u^*(u_1, v_1) - f_u^*(u_2, v_2))(t, x)| (d_b|(u_1 - u_2)(s, x)| + c_2|(v_1 - v_2)(s, x)|) dx \\ \leq \frac{\sigma d_b^2}{T} \|(u_1 - u_2)(s)\|_{L^2(\Omega)}^2 + \frac{\sigma c_2^2}{T} \|(v_1 - v_2)(s)\|_{L^2(\Omega)}^2 \\ + \frac{T}{2\sigma} \|(f_u^*(u_1, v_1) - f_u^*(u_2, v_2))(t)\|_{L^2(\Omega)}^2. \end{aligned} \qquad (8.15)$$
Plugging (8.15) into (8.14), we end up with

$$ \begin{aligned} & \iint_{\Omega_{\tau}} (f_{u}^{*}(u_{1}, v_{1}) - f_{u}^{*}(u_{2}, v_{2}))(t, x) \int_{t}^{\tau} (A_{1} - A_{2})(s, x) ds \, dxdt \\ & \leq \sigma d_{b}^{2} \|u_{1} - u_{2}\|_{L^{2}(\Omega_{\tau})}^{2} + \sigma c_{2}^{2} \|v_{1} - v_{2}\|_{L^{2}(\Omega_{\tau})}^{2} \\ & \qquad + \frac{T}{2\sigma} \int_{0}^{\tau} \int_{0}^{s} \|(f_{u}^{*}(u_{1}, v_{1}) - f_{u}^{*}(u_{2}, v_{2}))(t)\|_{L^{2}(\Omega)}^{2} dt ds. \end{aligned} \quad (8.16) $$

Now, gathering (8.9)-(8.13) and (8.16), we have that there exists $C_1 = C_1(\sigma, d_a, d_b, c_2) > 0$ such that, for all $T > 0$, it holds

$$ \begin{aligned} & (d_a - \sigma(2d_a + d_b^2)) \|u_1 - u_2\|_{L^2(\Omega_\tau)}^2 \le C_1 \|v_1 - v_2\|_{L^2(\Omega_\tau)}^2 + \frac{T d_b^2}{2\sigma d_a} \|u_1^{\text{in}} - u_2^{\text{in}}\|_{L^2(\Omega)}^2 \\ & \qquad + \frac{T}{2\sigma} \int_0^\tau \int_0^s \|(f_u^*(u_1, v_1) - f_u^*(u_2, v_2))(t)\|_{L^2(\Omega)}^2 dt ds. \end{aligned} \quad (8.17) $$

The same type of computations can be performed for the equation satisfied by $v_1 - v_2$ using the test functions

$$ w_2^\tau(t, x) = \begin{cases} \int_t^\tau (v_1(s, x) - v_2(s, x)) \, ds, & \text{if } 0 \le t \le \tau, \\ 0, & \text{if } \tau \le t \le T. \end{cases} \quad (8.18) $$

Indeed, from (8.3) we have

$$ \begin{aligned} & \iint_{\Omega_{\tau}} (v_1 - v_2)^2(t, x) \, dxdt - \int_{\Omega} (v_1^{\text{in}} - v_2^{\text{in}})(x) \int_{0}^{\tau} (v_1 - v_2)(s, x) \, ds \, dx \\ & \qquad + \iint_{\Omega_{\tau}} \nabla(v_1 - v_2)(t, x) \cdot \int_{t}^{\tau} \nabla(v_1 - v_2)(s, x) \, ds \, dxdt \quad (8.19) \\ & = \iint_{\Omega_{\tau}} (f_v^*(u_1, v_1) - f_v^*(u_2, v_2))(t, x) \int_{t}^{\tau} (v_1 - v_2)(s, x) \, ds \, dxdt. \end{aligned} $$

The first term in the left hand side of (8.19) is left as it is, while all the other terms are estimated similarly as before, to obtain

$$ \begin{aligned} & (1 - 2\sigma) \|v_1 - v_2\|_{L^2(\Omega_\tau)}^2 \le \frac{T}{4\sigma} \|v_1^{\text{in}} - v_2^{\text{in}}\|_{L^2(\Omega)}^2 \\ & \qquad + \frac{T}{4\sigma} \int_0^\tau \int_0^s \|(f_v^*(u_1, v_1) - f_v^*(u_2, v_2))(t)\|_{L^2(\Omega)}^2 dt ds. \end{aligned} \quad (8.20) $$

It remains to choose $\sigma > 0$ small enough to have $d_a - \sigma(2d_a + d_b^2) > 0$ in (8.17) and then $\lambda > 0$ large enough (depending on $C_1(\sigma, d_a, d_b, c_2)$), so that, adding (8.17), (8.20), the function $\phi_\lambda(\tau)$ in (8.1) satisfies, for all $\tau \in [0, T]$,

$$ \begin{aligned} \phi_\lambda(\tau) \le T C_2 \left( \|u_1^{\text{in}} - u_2^{\text{in}}\|_{L^2(\Omega)}^2 + \lambda \|v_1^{\text{in}} - v_2^{\text{in}}\|_{L^2(\Omega)}^2 \right) \\ + T C_3 \left( \int_0^\tau \int_0^s \|(f_u^*(u_1, v_1) - f_u^*(u_2, v_2))(t)\|_{L^2(\Omega)}^2 dt ds \right. \\ \left. + \lambda \int_0^\tau \int_0^s \|(f_v^*(u_1, v_1) - f_v^*(u_2, v_2))(t)\|_{L^2(\Omega)}^2 dt ds \right), \end{aligned} \quad (8.21) $$

where $C_2, C_3$ are positive constants depending only on $\sigma, d_a, d_b, c_2$.
Finally, as $f_u(u_a, u_b, v)$, $f_v(u_a, u_b, v)$ in (1.9), (1.4) are locally Lipschitz continuous (due to the quadratic terms) and recalling that $u_a^*, u_b^*$ are $C^1(\mathbb{R}_+^2)$ with bounded gradient, $f_u^*(u, v)$, $f_v^*(u, v)$ are also locally Lipschitz continuous. Therefore, there exists a positive constant $C_4$, depending on $\sigma$, $d_a$, $d_b$, $c_2$, $\lambda$ and $\|u_i\|_{L^\infty((0,T)\times\Omega)}$, $\|v_i\|_{L^\infty((0,\infty)\times\Omega)}$, such that (8.21) gives us, for all $\tau \in [0, T]$,

$$ \phi_{\lambda}(\tau) \le T C_2 \left( \|u_1^{\text{in}} - u_2^{\text{in}}\|_{L^2(\Omega)}^2 + \lambda \|v_1^{\text{in}} - v_2^{\text{in}}\|_{L^2(\Omega)}^2 \right) + T C_4 \int_0^{\tau} \phi_{\lambda}(s) ds. $$

Gronwall's lemma implies, for all $\tau \in [0, T]$,

$$ \|u_1 - u_2\|_{L^2(\Omega_\tau)}^2 + \lambda \|v_1 - v_2\|_{L^2(\Omega_\tau)}^2 \le T C_2 e^{T^2 C_4} \left( \|u_1^{\text{in}} - u_2^{\text{in}}\|_{L^2(\Omega)}^2 + \lambda \|v_1^{\text{in}} - v_2^{\text{in}}\|_{L^2(\Omega)}^2 \right) $$

and (2.26) follows.

## 9 The rate of convergence: proof of Theorem 2.7

Let $(u, v)$ be a nonnegative global classical solution of (1.10)–(1.12) satisfying (2.29), $(u_a^\epsilon, u_b^\epsilon, v^\epsilon)$ be the unique nonnegative global strict solution of (1.1)–(1.7) and $u^\epsilon = u_a^\epsilon + u_b^\epsilon$. We denote

$$ U^{\epsilon} := u^{\epsilon} - u, \quad V^{\epsilon} := v^{\epsilon} - v, \quad w := u_{b}^{*}(u, v), \quad W^{\epsilon} := u_{b}^{\epsilon} - w, \quad (9.1) $$

and $\delta := d_b - d_a > 0$ so that, by (9.1), (1.11),

$$ \begin{aligned} d_a u_a^\epsilon + d_b u_b^\epsilon &= d_a u^\epsilon + \delta u_b^\epsilon = d_a(U^\epsilon + u) + \delta(W^\epsilon + w), \\ A(u, v) &= d_a u_a^*(u, v) + d_b u_b^*(u, v) = d_a u + \delta w. \end{aligned} \quad (9.2) $$

It is worth noticing that by (9.2) and the definition of $w$ in (9.1), it holds

$$ \nabla A(u, v) = (d_a + \delta \partial_1 u_b^*(u, v)) \nabla u + \delta \partial_2 u_b^*(u, v) \nabla v. $$

As observed in (8.5), $d_a + \delta \partial_1 u_b^*(u, v) \in (d_a, d_b)$. Therefore, the homogeneous Neumann boundary conditions satisfied by $A(u, v)$ and $v$ give homogeneous Neumann boundary conditions for $u$ and consequently for $w$, $U^\epsilon$ and $W^\epsilon$.

Furthermore, under the constraint $u = u_a + u_b$, with $u_a, u_b \in \mathbb{R}_+$, the reaction functions $f_u = f_a + f_b, f_v$ and $Q$, defined in (1.4), (1.5), can be written as

$$ \begin{aligned} f_u(u_a, u_b, v) &= f_u(u - u_b, u_b, v) := \mathcal{F}_U(u, u_b, v) \\ f_v(u_a, u_b, v) &= f_v(u - u_b, u_b, v) := \mathcal{F}_V(u, u_b, v) \\ f_b(u_a, u_b, v) &= f_b(u - u_b, u_b, v) := \mathcal{F}_W(u, u_b, v) \\ Q(u_a, u_b, v) &= Q(u - u_b, u_b, v) := Q(u, u_b, v). \end{aligned} \quad (9.3) $$

Then, it is easily seen that the triplet $(U^\epsilon, V^\epsilon, W^\epsilon)$ satisfies over $(0, T) \times \Omega$ and for all $T > 0$, the fast reaction-diffusion system below

$$ \left\{ \begin{aligned} \partial_t U^\epsilon - \Delta(d_a U^\epsilon + \delta W^\epsilon) &= \mathcal{F}_U(U^\epsilon + u, W^\epsilon + w, V^\epsilon + v) - \mathcal{F}_U(u, w, v) \\ \partial_t V^\epsilon - d_v \Delta V^\epsilon &= \mathcal{F}_V(U^\epsilon + u, W^\epsilon + w, V^\epsilon + v) - \mathcal{F}_V(u, w, v) \\ \partial_t W^\epsilon - d_b \Delta W^\epsilon &= \mathcal{F}_W(U^\epsilon + u, W^\epsilon + w, V^\epsilon + v) \\ & \quad - \epsilon^{-1} Q(U^\epsilon + u, W^\epsilon + w, V^\epsilon + v) - (\partial_t w - d_b \Delta w) \end{aligned} \right. \quad (9.4) $$
together with the initial and boundary conditions

$$ \begin{cases} \nabla U^{\epsilon} \cdot \vec{n} = \nabla V^{\epsilon} \cdot \vec{n} = \nabla W^{\epsilon} \cdot \vec{n} = 0, & \text{in } (0, T) \times \partial\Omega, \\ U^{\epsilon}(0) = V^{\epsilon}(0) = 0, \quad W^{\epsilon}(0) = u_{b}^{\text{in}} - u_{b}^{*}(u^{\text{in}}, v^{\text{in}}), & \text{in } \Omega. \end{cases} \quad (9.5) $$

Theorem 2.7 will be proved estimating (term by term) the time evolution of the functional

$$ \mathcal{L}(t) := \frac{\gamma_1}{2} \|U^{\epsilon}(t)\|_{L^2(\Omega)}^2 + \frac{\gamma_2}{2} \|V^{\epsilon}(t)\|_{L^2(\Omega)}^2 + \frac{\epsilon\gamma_3}{2} \|W^{\epsilon}(t)\|_{L^2(\Omega)}^2 \\ + \frac{\epsilon\delta}{2} \|\nabla W^{\epsilon}(t)\|_{L^2(\Omega)}^2 + E(t) \quad (9.6) $$

where

$$ E(t) := - \int_{\Omega} [\mathcal{P}(U^{\epsilon} + u, w, V^{\epsilon} + v) - \mathcal{P}(u, w, v) - \partial_1 \mathcal{P}(u, w, v) U^{\epsilon} - \partial_3 \mathcal{P}(u, w, v) V^{\epsilon}] dx \quad (9.7) $$

$$ \mathcal{P}(x, y, z) := \int_{0}^{x} Q(\xi, y, z) d\xi, \quad (x, y, z) \in \mathbb{R}_{+}^{3}, \quad (9.8) $$

and $\gamma_1, \gamma_2, \gamma_3$ are strictly positive constants (independent of $\epsilon$) to be chosen later. The functional (9.6)-(9.8) is inspired by [18]. However we require minor properties for $Q$, than in [18]. The key tool employed to handle the functional is simply Taylor's formula.

**Step 1. Preliminaries.** To begin with, note that by assumptions (2.28)-(2.29) and estimate (2.6), for all $T > 0$, there exists $M_T > 0$, such that $u, v, u^{\epsilon}, v^{\epsilon} \in [0, M_T]$, for all $\epsilon \in (0, 1)$. As $0 \le u_b^{\epsilon} \le u^{\epsilon}$ and $0 \le w = u_b^{*}(u, v) \le u$, it follows that $u_b^{\epsilon}, w \in [0, M_T]$, for all $\epsilon \in (0, 1)$. Therefore, in the sequel we can innoque the boundedness over $[0, M_T]^3$ of the reaction functions in (9.3) and their derivatives and of $\partial_i Q, \partial_{ij} Q, \partial_{ijk} Q, i, j, k \in \{1, 2, 3\}$. Indeed, as we have assumed $A, B > 0$, for all $\alpha, \beta$, the function $Q$ belongs to $C^{\infty}(\mathbb{R}_+^3)$. In particular, from (7.2) it is easily seen that $\partial_1 Q < 0$ and that there exists $C(\alpha, \beta, A, B, M_T) > 0$ such that

$$ K_0 := \inf_{[0, M_T]^3} |\partial_1 Q(u_a, u_b, v)| \ge \inf_{[0, M_T]^3} \psi(a_{\text{fast}}u_a + c_{\text{fast}}v)/\Lambda(u_a, u_b, v) \ge C > 0. \quad (9.9) $$

On the other hand, using the inequalities

$$ \psi/\Lambda \le 1, \quad \phi/\Lambda \le 1, \quad a_{\text{fast}}u_a\psi'/\Lambda \le \alpha\psi/\Lambda \le \alpha, \quad \psi'/\Lambda \le \psi'/\psi \le \alpha A^{-1}, $$

it holds

$$ K_1 := \sup_{[0, M_T]^3} |\partial_1 Q(u_a, u_b, v)| \le 1 + \alpha + a_{\text{fast}} M_T \alpha A^{-1}. \quad (9.10) $$

Furthermore, by (7.2), (7.3), we have

$$ \inf_{\mathbb{R}_+^3} (-\partial_1 Q(u_a, u_b, v) + \partial_2 Q(u_a, u_b, v)) \ge 1. \quad (9.11) $$

By (2.29) again, for all $T > 0$, $\partial_t u, \partial_t v, \nabla u, \nabla v, \Delta u, \Delta v, \nabla \partial_t u, \nabla \partial_t v, \nabla \Delta u, \nabla \Delta v$ are bounded over $[0, T] \times \bar{\Omega}$. As a consequence, $\partial_t w, \Delta w, \nabla \partial_t w, \nabla \Delta w$ are also bounded over $[0, T] \times \bar{\Omega}$. Indeed, $w = u_b^*(u, v)$ and the gradient of the map $u_b^* : \mathbb{R}_+^2 \to$
$\mathbb{R}_+^2$ is given by (B.1), (B.2), with $q$ defined in (1.5). As $A, B > 0$, $q \in C^\infty(\mathbb{R}_+^3)$. In particular $\partial_2 q - \partial_1 q \ge A^\alpha + B^\beta > 0$. Therefore, $\partial_t u_b^*, \partial_{ij} u_b^*, \partial_{ijk} u_b^*, i, j, k \in \{1, 2\}$, are locally bounded and this is sufficient since $u, v \in [0, M_T]$.

Hereafter, the $\xi_i$, $i = 1, 2, 3$, appearing when applying Taylor's formula, belong all to $[0, M_T]$. The constants in the estimates will change from line to line and only the dependence on $d_a, d_b, d_v, \delta, \mathcal{K}_0, \mathcal{K}_1$ is kept explicit.

**Step 2. The evolution equation of $\|U^\epsilon\|_{L^2(\Omega)}^2$ and $\|V^\epsilon\|_{L^2(\Omega)}^2$.** We have

$$ 
\begin{aligned}
\frac{1}{2} \frac{d}{dt} \|U^\epsilon\|_{L^2(\Omega)}^2 &= -d_a \|\nabla U^\epsilon\|_{L^2(\Omega)}^2 - \delta \int_{\Omega} \nabla W^\epsilon \cdot \nabla U^\epsilon dx \\
&\quad + \int_{\Omega} [\mathcal{F}_U(U^\epsilon + u, W^\epsilon + w, V^\epsilon + v) - \mathcal{F}_U(u, w, v)] U^\epsilon dx \\
&= -d_a \|\nabla U^\epsilon\|_{L^2(\Omega)}^2 - \delta \int_{\Omega} \nabla W^\epsilon \cdot \nabla U^\epsilon dx \\
&\quad + \int_{\Omega} (D\mathcal{F}_U(\xi_1, \xi_2, \xi_3) \cdot (U^\epsilon, W^\epsilon, V^\epsilon)) U^\epsilon dx,
\end{aligned}
$$

and

$$ 
\begin{aligned}
\frac{1}{2} \frac{d}{dt} \|V^\epsilon\|_{L^2(\Omega)}^2 &= -d_v \|\nabla V^\epsilon\|_{L^2(\Omega)}^2 \\
&\quad + \int_{\Omega} [\mathcal{F}_V(U^\epsilon + u, W^\epsilon + w, V^\epsilon + v) - \mathcal{F}_V(u, w, v)] V^\epsilon dx \\
&= -d_v \|\nabla V^\epsilon\|_{L^2(\Omega)}^2 + \int_{\Omega} (D\mathcal{F}_V(\xi_1, \xi_2, \xi_3) \cdot (U^\epsilon, W^\epsilon, V^\epsilon)) V^\epsilon dx.
\end{aligned}
$$

Hence, the local boundedness of $D\mathcal{F}_U, D\mathcal{F}_V$ and ad hoc Young's inequality, give

$$ 
\begin{aligned}
\frac{1}{2} \frac{d}{dt} \|U^\epsilon\|_{L^2(\Omega)}^2 &\le -\frac{d_a}{2} \|\nabla U^\epsilon\|_{L^2(\Omega)}^2 + \frac{\delta^2}{2d_a} \|\nabla W^\epsilon\|_{L^2(\Omega)}^2 \\
&\quad + C(\|U^\epsilon\|_{L^2(\Omega)} + \|W^\epsilon\|_{L^2(\Omega)} + \|V^\epsilon\|_{L^2(\Omega)}) \|U^\epsilon\|_{L^2(\Omega)} \\
&\le -\frac{d_a}{2} \|\nabla U^\epsilon\|_{L^2(\Omega)}^2 + \frac{\delta^2}{2d_a} \|\nabla W^\epsilon\|_{L^2(\Omega)}^2 \\
&\quad + C(\|U^\epsilon\|_{L^2(\Omega)}^2 + \|V^\epsilon\|_{L^2(\Omega)}^2) + \frac{1}{4} \|W^\epsilon\|_{L^2(\Omega)}^2,
\end{aligned}
\quad (9.12) $$

and similarly

$$ 
\begin{aligned}
\frac{1}{2} \frac{d}{dt} \|V^\epsilon\|_{L^2(\Omega)}^2 &\le -d_v \|\nabla V^\epsilon\|_{L^2(\Omega)}^2 + C(\|U^\epsilon\|_{L^2(\Omega)}^2 + \|V^\epsilon\|_{L^2(\Omega)}^2) + \frac{1}{4} \|W^\epsilon\|_{L^2(\Omega)}^2. 
\end{aligned}
\quad (9.13) $$

**Step 3. The evolution equation of $\|W^\epsilon\|_{L^2(\Omega)}^2$.** We have

$$ 
\begin{aligned}
\frac{1}{2} \frac{d}{dt} \|W^\epsilon\|_{L^2(\Omega)}^2 &= -d_b \|\nabla W^\epsilon\|_{L^2(\Omega)}^2 + I_1 + I_2, \\
I_1 &:= \int_{\Omega} \mathcal{F}_W(U^\epsilon + u, W^\epsilon + w, V^\epsilon + v) W^\epsilon dx - \int_{\Omega} (\partial_t w - d_b \Delta w) W^\epsilon dx, \\
I_2 &:= -\frac{1}{\epsilon} \int_{\Omega} Q(U^\epsilon + u, W^\epsilon + w, V^\epsilon + v) W^\epsilon dx.
\end{aligned}
\quad (9.14) $$
By the local boundedness of $\mathcal{F}_W$ and the boundedness of $\partial_t w, \Delta w$

$$I_1 \leq C(1 + d_b) \|W^\varepsilon\|_{L^1(\Omega)} \quad (9.15)$$

Next, observing that $Q(u, w, v) = Q(u_a^*(u, v), u_b^*(u, v), v) = 0$ because $(u, v)$ is a solution of (1.10)–(1.12), we write $I_2$ as

$$I_2 = -\frac{1}{\varepsilon} \int_{\Omega} [Q(U^\varepsilon + u, W^\varepsilon + w, V^\varepsilon + v) - Q(U^\varepsilon + u, w, V^\varepsilon + v)] W^\varepsilon dx \\ -\frac{1}{\varepsilon} \int_{\Omega} [Q(U^\varepsilon + u, w, V^\varepsilon + v) - Q(u, w, v)] W^\varepsilon dx.$$

The latter allows us to use again Taylor's formula to obtain

$$I_2 = -\frac{1}{\varepsilon} \int_{\Omega} \partial_2 Q(U^\varepsilon + u, \xi_1, V^\varepsilon + v)(W^\varepsilon)^2 dx \\ -\frac{1}{\varepsilon} \int_{\Omega} [\partial_1 Q(\xi_2, w, \xi_3)U^\varepsilon + \partial_3 Q(\xi_2, w, \xi_3)V^\varepsilon]W^\varepsilon dx.$$

As by the definition of $Q$ in (9.3) it holds that $\partial_2 Q = -\partial_1 Q + \partial_2 Q$, using (9.11) together with the local boundedness of $DQ$, we get

$$I_2 \leq -\frac{1}{\varepsilon} \|W^\varepsilon\|_{L^2(\Omega)}^2 + \frac{C}{\varepsilon} (\|U^\varepsilon\|_{L^2(\Omega)} + \|V^\varepsilon\|_{L^2(\Omega)}) \|W^\varepsilon\|_{L^2(\Omega)}. \quad (9.16)$$

Plugging (9.15), (9.16) into (9.14), we end up with the estimate

$$\frac{1}{2} \frac{d}{dt} \|W^\varepsilon\|_{L^2(\Omega)}^2 \leq -d_b \|\nabla W^\varepsilon\|_{L^2(\Omega)}^2 + C(1+d_b) \|W^\varepsilon\|_{L^1(\Omega)} \\ - \frac{1}{\varepsilon} \|W^\varepsilon\|_{L^2(\Omega)}^2 + \frac{C}{\varepsilon} (\|U^\varepsilon\|_{L^2(\Omega)} + \|V^\varepsilon\|_{L^2(\Omega)}) \|W^\varepsilon\|_{L^2(\Omega)}.$$

Finally, we multiply the above inequality by $\varepsilon$ and use Young's inequality to get

$$\begin{aligned} \frac{\varepsilon}{2} \frac{d}{dt} \|W^\varepsilon\|_{L^2(\Omega)}^2 &\leq -d_b \varepsilon \|\nabla W^\varepsilon\|_{L^2(\Omega)}^2 + C(1+d_b)^2 \varepsilon^2 - \frac{1}{2} \|W^\varepsilon\|_{L^2(\Omega)}^2 \\ &\quad + C(\|U^\varepsilon\|_{L^2(\Omega)} + \|V^\varepsilon\|_{L^2(\Omega)}) \|W^\varepsilon\|_{L^2(\Omega)} \\ &\leq -d_b \varepsilon \|\nabla W^\varepsilon\|_{L^2(\Omega)}^2 + C(1+d_b)^2 \varepsilon^2 - \frac{1}{4} \|W^\varepsilon\|_{L^2(\Omega)}^2 + C(\|U^\varepsilon\|_{L^2(\Omega)}^2 + \|V^\varepsilon\|_{L^2(\Omega)}^2). \end{aligned} \quad (9.17)$$

**Step 4.** The evolution equation of $\|\nabla W^\varepsilon\|_{L^2(\Omega)}^2$. Multiplying the equation for $W^\varepsilon$ in (9.4) by $-\Delta W^\varepsilon$, we have

$$\frac{1}{2} \frac{d}{dt} \|\nabla W^\varepsilon\|_{L^2(\Omega)}^2 = -d_b \|\Delta W^\varepsilon\|_{L^2(\Omega)}^2 + J_1 + J_2 + J_3, \quad (9.18)$$

$$J_1 := - \int_{\Omega} \mathcal{F}_W (U^\varepsilon + u, W^\varepsilon + w, V^\varepsilon + v) \Delta W^\varepsilon dx$$

$$J_2 := \int_{\Omega} (\partial_t w - d_b \Delta w) \Delta W^\varepsilon dx$$

$$J_3 := \frac{1}{\varepsilon} \int_{\Omega} Q(U^\varepsilon + u, W^\varepsilon + w, V^\varepsilon + v) \Delta W^\varepsilon dx.$$
Using Taylor's formula into $J_1$ and the local boundedness of $D\mathcal{F}_W$, we have

$$J_1 = - \int_{\Omega} \mathcal{F}_W(u, w, v) \Delta W^{\varepsilon} dx - \int_{\Omega} (D\mathcal{F}_W(\xi_1, \xi_2, \xi_3) \cdot (U^{\varepsilon}, W^{\varepsilon}, V^{\varepsilon})) \Delta W^{\varepsilon} dx \\ \leq C \|\nabla W^{\varepsilon}\|_{L^1(\Omega)} + C(\|U^{\varepsilon}\|_{L^2(\Omega)} + \|W^{\varepsilon}\|_{L^2(\Omega)} + \|V^{\varepsilon}\|_{L^2(\Omega)}) \|\Delta W^{\varepsilon}\|_{L^2(\Omega)}.$$

Moreover, by the boundedness of $\nabla\partial_t w, \nabla\Delta w$, we obtain for $J_2$

$$J_2 \leq C(1 + d_b) \|\nabla W^{\varepsilon}\|_{L^1(\Omega)}.$$

Therefore,

$$\begin{aligned} J_1 + J_2 &\leq C(1 + d_b) \|\nabla W^{\varepsilon}\|_{L^1(\Omega)} \\ &\quad + C(\|U^{\varepsilon}\|_{L^2(\Omega)} + \|W^{\varepsilon}\|_{L^2(\Omega)} + \|V^{\varepsilon}\|_{L^2(\Omega)}) \|\Delta W^{\varepsilon}\|_{L^2(\Omega)} \\ &\leq C(1 + d_b)^2 \varepsilon + \frac{1}{4\varepsilon} \|\nabla W^{\varepsilon}\|_{L^2(\Omega)}^2 \quad (9.19) \\ &\quad + C d_b^{-1} (\|U^{\varepsilon}\|_{L^2(\Omega)}^2 + \|W^{\varepsilon}\|_{L^2(\Omega)}^2 + \|V^{\varepsilon}\|_{L^2(\Omega)}^2) + \frac{d_b}{2} \|\Delta W^{\varepsilon}\|_{L^2(\Omega)}^2. \end{aligned}$$

The term $J_3$ is the more challenging and we proceed as for $I_2$ in Step 4, i.e.

$$\begin{aligned} J_3 &= \frac{1}{\varepsilon} \int_{\Omega} [\mathcal{Q}(U^{\varepsilon} + u, W^{\varepsilon} + w, V^{\varepsilon} + v) - \mathcal{Q}(U^{\varepsilon} + u, w, V^{\varepsilon} + v)] \Delta W^{\varepsilon} dx \\ &\quad + \frac{1}{\varepsilon} \int_{\Omega} [\mathcal{Q}(U^{\varepsilon} + u, w, V^{\varepsilon} + v) - \mathcal{Q}(u, w, v)] \Delta W^{\varepsilon} dx \\ &=: J_3^1 + J_3^2. \end{aligned} \quad (9.20)$$

For $J_3^1$ we write

$$J_3^1 = - \frac{1}{\varepsilon} \int_{\Omega} \nabla [\mathcal{Q}(U^{\varepsilon} + u, W^{\varepsilon} + w, V^{\varepsilon} + v) - \mathcal{Q}(U^{\varepsilon} + u, w, V^{\varepsilon} + v)] \cdot \nabla W^{\varepsilon} dx,$$

so that

$$\begin{aligned} J_3^1 = &-\frac{1}{\varepsilon} \int_{\Omega} [\partial_1 \mathcal{Q}(U^{\varepsilon} + u, W^{\varepsilon} + w, V^{\varepsilon} + v) - \partial_1 \mathcal{Q}(U^{\varepsilon} + u, w, V^{\varepsilon} + v)] \nabla(U^{\varepsilon} + u) \cdot \nabla W^{\varepsilon} dx \\ &-\frac{1}{\varepsilon} \int_{\Omega} \partial_2 \mathcal{Q}(U^{\varepsilon} + u, W^{\varepsilon} + w, V^{\varepsilon} + v) |\nabla W^{\varepsilon}|^2 dx \\ &-\frac{1}{\varepsilon} \int_{\Omega} [\partial_2 \mathcal{Q}(U^{\varepsilon} + u, W^{\varepsilon} + w, V^{\varepsilon} + v) - \partial_2 \mathcal{Q}(U^{\varepsilon} + u, w, V^{\varepsilon} + v)] \nabla w \cdot \nabla W^{\varepsilon} dx \\ &-\frac{1}{\varepsilon} \int_{\Omega} [\partial_3 \mathcal{Q}(U^{\varepsilon} + u, W^{\varepsilon} + w, V^{\varepsilon} + v) - \partial_3 \mathcal{Q}(U^{\varepsilon} + u, w, V^{\varepsilon} + v)] \nabla(V^{\varepsilon} + v) \cdot \nabla W^{\varepsilon} dx. \end{aligned}$$
Moreover, using $\partial_2 Q = -\partial_1 Q + \partial_2 Q$ and (9.11), and rearranging the terms

$$ 
\begin{aligned}
J_3^1 \le & -\frac{1}{\varepsilon} \int_{\Omega} [\partial_1 Q(U^\varepsilon + u, W^\varepsilon + w, V^\varepsilon + v) - \partial_1 Q(U^\varepsilon + u, w, V^\varepsilon + v)] \nabla U^\varepsilon \cdot \nabla W^\varepsilon dx \\
& -\frac{1}{\varepsilon} \int_{\Omega} \partial_{12} Q(U^\varepsilon + u, \xi_1, V^\varepsilon + v) W^\varepsilon (\nabla u \cdot \nabla W^\varepsilon) dx \\
& -\frac{1}{\varepsilon} \|\nabla W^\varepsilon\|_{L^2(\Omega)}^2 \\
& -\frac{1}{\varepsilon} \int_{\Omega} \partial_{22} Q(U^\varepsilon + u, \xi_2, V^\varepsilon + v) W^\varepsilon (\nabla w \cdot \nabla W^\varepsilon) dx \\
& -\frac{1}{\varepsilon} \int_{\Omega} [\partial_3 Q(U^\varepsilon + u, W^\varepsilon + w, V^\varepsilon + v) - \partial_3 Q(U^\varepsilon + u, w, V^\varepsilon + v)] \nabla V^\varepsilon \cdot \nabla W^\varepsilon dx \\
& -\frac{1}{\varepsilon} \int_{\Omega} \partial_{32} Q(U^\varepsilon + u, \xi_3, V^\varepsilon + v) W^\varepsilon (\nabla v \cdot \nabla W^\varepsilon) dx.
\end{aligned}
$$

Next, using (9.10), the boundedness of $\nabla u, \nabla v, \nabla w$ and the local boundedness of $DQ, D^2Q$, we obtain

$$ 
\begin{aligned}
J_3^1 \le & \frac{2K_1}{\varepsilon} \|\nabla U^\varepsilon\|_{L^2(\Omega)} \|\nabla W^\varepsilon\|_{L^2(\Omega)} + \frac{C}{\varepsilon} \|W^\varepsilon\|_{L^2(\Omega)} \|\nabla W^\varepsilon\|_{L^2(\Omega)} \\
& - \frac{1}{\varepsilon} \|\nabla W^\varepsilon\|_{L^2(\Omega)}^2 + \frac{C}{\varepsilon} \|\nabla V^\varepsilon\|_{L^2(\Omega)} \|\nabla W^\varepsilon\|_{L^2(\Omega)} \\
\le & \frac{6K_1^2}{\varepsilon} \|\nabla U^\varepsilon\|_{L^2(\Omega)}^2 - \frac{1}{2\varepsilon} \|\nabla W^\varepsilon\|_{L^2(\Omega)}^2 + \frac{C}{\varepsilon} \|W^\varepsilon\|_{L^2(\Omega)}^2 + \frac{C}{\varepsilon} \|\nabla V^\varepsilon\|_{L^2(\Omega)}^2.
\end{aligned}
\quad (9.21) $$

The term $J_3^2 = \varepsilon^{-1} \int_{\Omega} [Q(U^\varepsilon + u, w, V^\varepsilon + v) - Q(u, w, v)] \Delta W^\varepsilon dx$ will be absorbed by the evolution equation of $E(t)$. Therefore, we let it as it is. Plugging into (9.18) multiplied by $\varepsilon$, the estimates (9.19), (9.21) and the definition of $J_3^2$, we have

$$ 
\begin{aligned}
\frac{\varepsilon}{2} \frac{d}{dt} \|\nabla W^\varepsilon\|_{L^2(\Omega)}^2 \le & -\frac{d_b \varepsilon}{2} \|\Delta W^\varepsilon\|_{L^2(\Omega)}^2 + C(1+d_b)^2 \varepsilon^2 - \frac{1}{4} \|\nabla W^\varepsilon\|_{L^2(\Omega)}^2 \\
& + 6K_1^2 \|\nabla U^\varepsilon\|_{L^2(\Omega)}^2 + C \|\nabla V^\varepsilon\|_{L^2(\Omega)}^2 \\
& + C d_b^{-1} \varepsilon (\|U^\varepsilon\|_{L^2(\Omega)}^2 + \|V^\varepsilon\|_{L^2(\Omega)}^2) + C(d_b^{-1} \varepsilon + 1) \|W^\varepsilon\|_{L^2(\Omega)}^2 \\
& + \int_{\Omega} [Q(U^\varepsilon + u, w, V^\varepsilon + v) - Q(u, w, v)] \Delta W^\varepsilon dx.
\end{aligned}
\quad (9.22) $$

**Step 5. The evolution equation of $E(t)$ defined in (9.7).** Let us write

$$ \frac{d}{dt} E(t) = \sum_{i=1}^{5} L_i(t), \quad (9.23) $$

with

$$ L_1 := - \int_{\Omega} [\partial_1 P(U^\varepsilon + u, w, V^\varepsilon + v) - \partial_1 P(u, w, v)] \partial_t U^\varepsilon dx $$

$$ L_2 := - \int_{\Omega} [\partial_1 P(U^\varepsilon + u, w, V^\varepsilon + v) - \partial_1 P(u, w, v) \\ - \partial_{11} P(u, w, v) U^\varepsilon - \partial_{31} P(u, w, v) V^\varepsilon] \partial_t u dx $$

$$ L_3 := - \int_{\Omega} [\partial_2 P(U^\varepsilon + u, w, V^\varepsilon + v) - \partial_2 P(u, w, v) \\ - \partial_{12} P(u, w, v) U^\varepsilon - \partial_{32} P(u, w, v) V^\varepsilon] \partial_t w dx $$
$$L_4 := - \int_{\Omega} [\partial_3 \mathcal{P}(U^{\varepsilon} + u, w, V^{\varepsilon} + v) - \partial_3 \mathcal{P}(u, w, v)] \partial_t V^{\varepsilon} dx$$

$$L_5 := - \int_{\Omega} [\partial_3 \mathcal{P}(U^{\varepsilon} + u, w, V^{\varepsilon} + v) - \partial_3 \mathcal{P}(u, w, v) \\ - \partial_{13} \mathcal{P}(u, w, v) U^{\varepsilon} - \partial_{33} \mathcal{P}(u, w, v) V^{\varepsilon}] \partial_t v \, dx.$$

The terms $L_2, L_3, L_5$ are easily controlled using second order Taylor's formula applied to $\partial_1 \mathcal{P} = \mathcal{Q}, \partial_2 \mathcal{P}, \partial_3 \mathcal{P}$ respectively, to obtain

$$ \begin{aligned} L_2 &\le C \int_{\Omega} ((U^{\varepsilon})^2 + (V^{\varepsilon})^2) |\partial_t u| \, dx \le C (\|U^{\varepsilon}\|_{L^2(\Omega)}^2 + \|V^{\varepsilon}\|_{L^2(\Omega)}^2), \\ L_3 &\le C \int_{\Omega} ((U^{\varepsilon})^2 + (V^{\varepsilon})^2) |\partial_t w| \, dx \le C (\|U^{\varepsilon}\|_{L^2(\Omega)}^2 + \|V^{\varepsilon}\|_{L^2(\Omega)}^2), \quad (9.24) \\ L_5 &\le C \int_{\Omega} ((U^{\varepsilon})^2 + (V^{\varepsilon})^2) |\partial_t v| \, dx \le C (\|U^{\varepsilon}\|_{L^2(\Omega)}^2 + \|V^{\varepsilon}\|_{L^2(\Omega)}^2). \end{aligned} $$

Next, by (9.8) and the equation for $U^{\varepsilon}$ in (9.4), we have for $L_1$

$$ \begin{aligned} L_1 &= - \int_{\Omega} [\mathcal{Q}(U^{\varepsilon} + u, w, V^{\varepsilon} + v) - \mathcal{Q}(u, w, v)] [d_a \Delta U^{\varepsilon} + \delta \Delta W^{\varepsilon} \\ &\quad + \mathcal{F}_U(U^{\varepsilon} + u, W^{\varepsilon} + w, V^{\varepsilon} + v) - \mathcal{F}_U(u, w, v)] \, dx \quad (9.25) \\ &=: L_1^1 + L_1^2 + L_1^3. \end{aligned} $$

The $L_1^2$ term is fundamental since it allows us to get rid of the term $J_3^2$ defined in (9.20). Indeed, $L_1^2$ reads as

$$L_1^2 = -\delta \int_{\Omega} [\mathcal{Q}(U^{\varepsilon} + u, w, V^{\varepsilon} + v) - \mathcal{Q}(u, w, v)] \Delta W^{\varepsilon} \, dx = -\varepsilon \delta J_3^2. \quad (9.26)$$

The control of $L_1^3$ follows simply as

$$ \begin{aligned} L_1^3 &= - \int_{\Omega} [\mathcal{Q}(U^{\varepsilon} + u, w, V^{\varepsilon} + v) - \mathcal{Q}(u, w, v)] [\mathcal{F}_U(U^{\varepsilon} + u, W^{\varepsilon} + w, V^{\varepsilon} + v) - \mathcal{F}_U(u, w, v)] \, dx \\ &\le C \int_{\Omega} (|U^{\varepsilon}| + |V^{\varepsilon}|)(|U^{\varepsilon}| + |V^{\varepsilon}| + |W^{\varepsilon}|) \, dx \quad (9.27) \\ &\le C (\|U^{\varepsilon}\|_{L^2(\Omega)}^2 + \|V^{\varepsilon}\|_{L^2(\Omega)}^2) + \frac{1}{8} \|W^{\varepsilon}\|_{L^2(\Omega)}^2. \end{aligned} $$

The control of $L_1^1$ follows the same computations done for $J_3^1$ defined in (9.20). First we write

$$ \begin{aligned} L_1^1 &= d_a \int_{\Omega} \nabla [\mathcal{Q}(U^{\varepsilon} + u, w, V^{\varepsilon} + v) - \mathcal{Q}(u, w, v)] \cdot \nabla U^{\varepsilon} \, dx \\ &= d_a \int_{\Omega} \partial_1 \mathcal{Q}(U^{\varepsilon} + u, w, V^{\varepsilon} + v) |\nabla U^{\varepsilon}|^2 \, dx \\ &\quad + d_a \int_{\Omega} [\partial_1 \mathcal{Q}(U^{\varepsilon} + u, w, V^{\varepsilon} + v) - \partial_1 \mathcal{Q}(u, w, v)] \nabla u \cdot \nabla U^{\varepsilon} \, dx \\ &\quad + d_a \int_{\Omega} [\partial_2 \mathcal{Q}(U^{\varepsilon} + u, w, V^{\varepsilon} + v) - \partial_2 \mathcal{Q}(u, w, v)] \nabla w \cdot \nabla U^{\varepsilon} \, dx \\ &\quad + d_a \int_{\Omega} \partial_3 \mathcal{Q}(U^{\varepsilon} + u, w, V^{\varepsilon} + v) \nabla V^{\varepsilon} \cdot \nabla U^{\varepsilon} \, dx \\ &\quad + d_a \int_{\Omega} [\partial_3 \mathcal{Q}(U^{\varepsilon} + u, w, V^{\varepsilon} + v) - \partial_3 \mathcal{Q}(u, w, v)] \nabla v \cdot \nabla U^{\varepsilon} \, dx. \end{aligned} $$
Next, since $\partial_1 Q = \partial_1 Q$, using (9.9) and skipping few details, we have

$$L_1^1 \le -d_a K_0 \|\nabla U^\epsilon\|_{L^2(\Omega)}^2 + d_a C \int_{\Omega} (|U^\epsilon| + |V^\epsilon|)(|\nabla u| + |\nabla w| + |\nabla v|) |\nabla U^\epsilon| dx \\ + d_a C \|\nabla U^\epsilon\|_{L^2(\Omega)} \|\nabla V^\epsilon\|_{L^2(\Omega)},$$

and by the boundedness of $\nabla u, \nabla w, \nabla v$ and ad hoc Young's inequalities, we obtain

$$L_1^1 \le -\frac{d_a K_0}{2} \|\nabla U^\epsilon\|_{L^2(\Omega)}^2 + \frac{C}{K_0} (\|U^\epsilon\|_{L^2(\Omega)}^2 + \|V^\epsilon\|_{L^2(\Omega)}^2) + C \frac{d_a}{K_0} \|\nabla V^\epsilon\|_{L^2(\Omega)}^2. \quad (9.28)$$

Plugging (9.26), (9.27), (9.28) into (9.25) we end up with

$$L_1 \le -\frac{d_a K_0}{2} \|\nabla U^\epsilon\|_{L^2(\Omega)}^2 + C \frac{d_a}{K_0} \|\nabla V^\epsilon\|_{L^2(\Omega)}^2 \\ + C(1 + K_0^{-1})(\|U^\epsilon\|_{L^2(\Omega)}^2 + \|V^\epsilon\|_{L^2(\Omega)}^2) + \frac{1}{8} \|W^\epsilon\|_{L^2(\Omega)}^2 \quad (9.29) \\ - \delta \int_{\Omega} [Q(U^\epsilon + u, w, V^\epsilon + v) - Q(u, w, v)] \Delta W^\epsilon dx.$$

We consider now the term $L_4$ and use the equation for $V^\epsilon$ in (9.4)

$$L_4 = - \int_{\Omega} [\partial_3 P(U^\epsilon + u, w, V^\epsilon + v) - \partial_3 P(u, w, v)] [d_v \Delta V^\epsilon \\ + F_V(U^\epsilon + u, W^\epsilon + w, V^\epsilon + v) - F_V(u, w, v)] dx \quad (9.30) \\ =: L_4^1 + L_4^2.$$

For the $L_4^1$ term we have

$$\begin{aligned} L_4^1 &= d_v \int_{\Omega} \nabla [\partial_3 P(U^\epsilon + u, w, V^\epsilon + v) - \partial_3 P(u, w, v)] \cdot \nabla V^\epsilon dx \\ &= d_v \int_{\Omega} \partial_{31} P(U^\epsilon + u, w, V^\epsilon + v) \nabla U^\epsilon \cdot \nabla V^\epsilon dx \\ &\quad + d_v \int_{\Omega} [\partial_{31} P(U^\epsilon + u, w, V^\epsilon + v) - \partial_{31} P(u, w, v)] \nabla u \cdot \nabla V^\epsilon dx \\ &\quad + d_v \int_{\Omega} [\partial_{32} P(U^\epsilon + u, w, V^\epsilon + v) - \partial_{32} P(u, w, v)] \nabla w \cdot \nabla V^\epsilon dx \\ &\quad + d_v \int_{\Omega} \partial_{33} P(U^\epsilon + u, w, V^\epsilon + v) |\nabla V^\epsilon|^2 dx \\ &\quad + d_v \int_{\Omega} [\partial_{33} P(U^\epsilon + u, w, V^\epsilon + v) - \partial_{33} P(u, w, v)] \nabla v \cdot \nabla V^\epsilon dx, \end{aligned}$$

so that

$$\begin{aligned} L_4^1 \le & \frac{d_a K_0}{4} \|\nabla U^\epsilon\|_{L^2(\Omega)}^2 + C \frac{d_v^2}{d_a K_0} \|\nabla V^\epsilon\|_{L^2(\Omega)}^2 \\ & + C d_v \int_{\Omega} (|U^\epsilon| + |V^\epsilon|)(|\nabla u| + |\nabla w| + |\nabla v|) |\nabla V^\epsilon| dx + C d_v \|\nabla V^\epsilon\|_{L^2(\Omega)}^2 \\ \le & \frac{d_a K_0}{4} \|\nabla U^\epsilon\|_{L^2(\Omega)}^2 + C d_v (1 + d_v (d_a K_0)^{-1}) \|\nabla V^\epsilon\|_{L^2(\Omega)}^2 \\ & + C d_v (\|U^\epsilon\|_{L^2(\Omega)}^2 + \|V^\epsilon\|_{L^2(\Omega)}^2). \end{aligned} \quad (9.31)$$
For the $L_4^2$ term in (9.30), skipping again few details, we have

$$ 
\begin{aligned}
L_4^2 &\le C \int_{\Omega} (|U^\epsilon| + |V^\epsilon|)(|U^\epsilon| + |V^\epsilon| + |W^\epsilon|) dx \\
&\le C (\|U^\epsilon\|_{L^2(\Omega)}^2 + \|V^\epsilon\|_{L^2(\Omega)}^2) + \frac{1}{8} \|W^\epsilon\|_{L^2(\Omega)}^2.
\end{aligned} 
\quad (9.32) $$

Plugging (9.31), (9.32) into (9.30) we obtain

$$ 
\begin{aligned}
L_4 \le & \frac{d_a K_0}{4} \|\nabla U^\epsilon\|_{L^2(\Omega)}^2 + C d_v (1 + d_v (d_a K_0)^{-1}) \|\nabla V^\epsilon\|_{L^2(\Omega)}^2 \\
& + C (1 + d_v) (\|U^\epsilon\|_{L^2(\Omega)}^2 + \|V^\epsilon\|_{L^2(\Omega)}^2) + \frac{1}{8} \|W^\epsilon\|_{L^2(\Omega)}^2.
\end{aligned} 
\quad (9.33) $$

Gathering (9.24), (9.29), (9.33) in (9.23), we end up with the estimate of $E'(t)$

$$ 
\begin{aligned}
\frac{d}{dt} E(t) \le & -\frac{d_a K_0}{4} \|\nabla U^\epsilon\|_{L^2(\Omega)}^2 + C(d_v + d_v^2 (d_a K_0)^{-1} + d_a K_0^{-1}) \|\nabla V^\epsilon\|_{L^2(\Omega)}^2 \\
& + C(1 + d_v + K_0^{-1}) (\|U^\epsilon\|_{L^2(\Omega)}^2 + \|V^\epsilon\|_{L^2(\Omega)}^2) + \frac{1}{4} \|W^\epsilon\|_{L^2(\Omega)}^2 \\
& - \delta \int_{\Omega} [Q(U^\epsilon + u, w, V^\epsilon + v) - Q(u, w, v)] \Delta W^\epsilon dx.
\end{aligned} 
\quad (9.34) $$

**Step 6.** A positive lower bound. Let $T > 0$ and $\gamma_1, \gamma_2 \in \mathbb{R}_+^2$. We claim that, if $\gamma_2 > 0$ is large enough, there exists $c_0 = c_0(\gamma_1, K_0, T) > 0$ such that

$$ E(t) + \frac{\gamma_1}{2} \|U^\epsilon\|_{L^2(\Omega)}^2 + \frac{\gamma_2}{2} \|V^\epsilon\|_{L^2(\Omega)}^2 \ge c_0 (\|U^\epsilon\|_{L^2(\Omega)}^2 + \|V^\epsilon\|_{L^2(\Omega)}^2). \quad (9.35) $$

Indeed, from (9.7) and a second order Taylor's formula, we have

$$ E(t) = - \int_{\Omega} [\partial_{11} \mathcal{P}(\xi_1, w, \xi_2) (U^\epsilon)^2 + 2 \partial_{13} \mathcal{P}(\xi_1, w, \xi_2) U^\epsilon V^\epsilon + \partial_{33} \mathcal{P}(\xi_1, w, \xi_2) (V^\epsilon)^2] dx. $$

Moreover, from (9.3), (9.8) and (9.9), we have

$$ -\partial_{11} \mathcal{P}(\xi_1, w, \xi_2) = -\partial_1 \mathcal{Q}(\xi_1, w, \xi_2) = -\partial_1 \mathcal{Q}(\xi_1, w, \xi_2) \ge K_0 > 0. \quad (9.36) $$

Therefore, using (9.36) and the local boundedness of $DQ$ and $D^2Q$, we obtain

$$ E(t) \ge K_0 \|U^\epsilon\|_{L^2(\Omega)}^2 - C_1 \|U^\epsilon\|_{L^2(\Omega)} \|V^\epsilon\|_{L^2(\Omega)} - C_2 \|V^\epsilon\|_{L^2(\Omega)}^2, $$

and, for $\theta > 0$ to be chosen,

$$ 
\begin{aligned}
E(t) + \frac{\gamma_1}{2} \|U^\epsilon\|_{L^2(\Omega)}^2 + \frac{\gamma_2}{2} \|V^\epsilon\|_{L^2(\Omega)}^2 \ge & (K_0 + \frac{\gamma_1}{2} - \theta) \|U^\epsilon\|_{L^2(\Omega)}^2 \\
& + (\frac{\gamma_2}{2} - C_2 - \frac{C_1^2}{4\theta}) \|V^\epsilon\|_{L^2(\Omega)}^2.
\end{aligned} 
$$

In order to obtain (9.35) it is sufficient to choose $\theta > 0$ small enough so that $K_0 + \frac{\gamma_1}{2} - \theta > 0$ and then $\gamma_2 > 0$ large enough so that $\frac{\gamma_2}{2} - C_2 - \frac{C_1^2}{4\theta} > 0$.

**Step 7.** End of the proof. Recall definition (9.6) of $\mathcal{L}(t)$. By (9.12), (9.13), (9.17), (9.22) and (9.34), neglecting useless negative terms, we have

$$ 
\begin{aligned}
\frac{d}{dt} \mathcal{L}(t) \le & c_1(\gamma_1) \|\nabla U^\epsilon\|_{L^2(\Omega)}^2 + c_2(\gamma_2) \|\nabla V^\epsilon\|_{L^2(\Omega)}^2 + c_3(\gamma_1) \|\nabla W^\epsilon\|_{L^2(\Omega)}^2 \\
& + c_4(\gamma_1, \gamma_2, \gamma_3) (\|U^\epsilon\|_{L^2(\Omega)}^2 + \|V^\epsilon\|_{L^2(\Omega)}^2) \\
& + c_5(\gamma_1, \gamma_2, \gamma_3) \|W^\epsilon\|_{L^2(\Omega)}^2 + C(1 + d_b)^2(\gamma_3 + \delta)\epsilon^2,
\end{aligned} 
\quad (9.37) $$
where

$$
\begin{align}
c_1(\gamma_1) &= -\frac{d_a}{4} \mathcal{K}_0 - \frac{d_a}{2} \gamma_1 + 6\delta \mathcal{K}_1^2, \nonumber \\
c_2(\gamma_2) &= -d_v \gamma_2 + C(\delta + d_v + d_v^2(d_a \mathcal{K}_0)^{-1} + d_a \mathcal{K}_0^{-1}), \nonumber \\
c_3(\gamma_1) &= \frac{\delta^2}{2d_a} \gamma_1 - \frac{\delta}{4}, \nonumber \\
c_5(\gamma_1, \gamma_2, \gamma_3) &= \frac{1}{4}(1 + \gamma_1 + \gamma_2 - \gamma_3) + \delta C(1 + \varepsilon d_b^{-1}). \tag{9.38}
\end{align}
$$

We need now to determine $\gamma_1, \gamma_2, \gamma_3 \in \mathbb{R}_+$ (independent of $\varepsilon$) such that the constant $c_1, c_2, c_3, c_5$ are negative. We will see that this is possible provided that the ratio $d_b/d_a$ is small enough (see (2.30)). Since the constant $c_4$ is positive for all $\gamma_1, \gamma_2, \gamma_3 \in \mathbb{R}_+$, we do not need the explicit formula.

Note that $c_1$ and $c_3$ in (9.38) are both negative if and only if

$$
\frac{12 \kappa_1^2 \delta}{d_a} - \frac{\kappa_0}{2} < \gamma_1 < \frac{d_a}{2\delta}. \qquad (9.39)
$$

If $\frac{12\kappa_1^2\delta}{d_a} - \frac{\kappa_0}{2} \le 0$, i.e. $\frac{d_b}{d_a} \le 1 + \frac{\kappa_0}{24\kappa_1^2}$, then it is sufficient to choose $\gamma_1 \in (0, \frac{d_a}{2\delta})$
to have $c_1, c_3$ negative. On the other hand, if $\frac{12\kappa_1^2\delta}{d_a} - \frac{\kappa_0}{2} > 0$, then it is easily
seen that $\frac{12\kappa_1^2\delta}{d_a} - \frac{\kappa_0}{2} < \frac{d_a}{2\delta}$ if $\frac{d_b}{d_a} < 1 + x_+$, where $x_+$ is the positive root of the
polynomial function $(12\kappa_1^2 x^2 - \frac{\kappa_0}{2} x - \frac{1}{2})$ and satisfies $x_+ > \frac{\kappa_0}{24\kappa_1^2}$. Therefore, it is
again possible to find $\gamma_1 > 0$ satisfying (9.39) and giving negative $c_1$ and $c_3$.

For this $\gamma_1 > 0$, we choose $\gamma_2 > 0$ large enough so that at the same time (9.35)
holds true and $c_2$ is negative, and finally $\gamma_3 > 0$ large enough so that $c_5$ is negative
as well.

Now, by (9.37) and (9.35), $\mathcal{L}(t)$ satisfies

$$
\begin{align*}
\frac{d}{dt} \mathcal{L}(t) &\le c_4 (\|U^\epsilon\|_{L^2(\Omega)}^2 + \|V^\epsilon\|_{L^2(\Omega)}^2) + C(1+d_b)^2 \epsilon^2 (\gamma_3 + \delta) \\
&\le c_4 c_0^{-1} (E(t) + \frac{\gamma_1}{2} \|U^\epsilon\|_{L^2(\Omega)}^2 + \frac{\gamma_2}{2} \|V^\epsilon\|_{L^2(\Omega)}^2) + C(1+d_b)^2 \epsilon^2 (\gamma_3 + \delta) \\
&\le c_4 c_0^{-1} \mathcal{L}(t) + C(1+d_b)^2 (\gamma_3 + \delta) \epsilon^2.
\end{align*}
$$

Integrating the differential inequality above over $(0, t)$, we obtain that there exists $C(T) > 0$ such that, for all $t \in (0, T)$, it holds

$$
E(t) + \frac{\gamma_1}{2} \|U^\epsilon(t)\|_{L^2(\Omega)}^2 + \frac{\gamma_2}{2} \|V^\epsilon(t)\|_{L^2(\Omega)}^2 \le \mathcal{L}(t) \le C(T)(\epsilon^2 + \mathcal{L}(0)).
$$

Observing that (see (9.6), (9.1) and (2.27))

$$
\mathcal{L}(0) = \frac{\varepsilon\gamma_3}{2} \|W^\varepsilon(0)\|_{L^2(\Omega)}^2 + \frac{\varepsilon\delta}{2} \|\nabla W^\varepsilon(0)\|_{L^2(\Omega)}^2 = C(T) \varepsilon \varepsilon_{\text{in}}^2, \quad (9.40)
$$

by (9.35) again, the above inequality implies

$$
\|U^{\epsilon}(t)\|_{L^2(\Omega)}^2 + \|V^{\epsilon}(t)\|_{L^2(\Omega)}^2 \le C(T)(\epsilon^2 + \epsilon \epsilon_{\text{in}}^2), \quad t \in (0, T). \quad (9.41)
$$

Next, we plug (9.41) into (9.17), we neglect useless negative terms and we obtain

$$
\frac{d}{dt} \|W^{\epsilon}(t)\|_{L^2(\Omega)}^2 \le C(T)(\epsilon + \epsilon_{\text{in}}^2) - \frac{1}{2\epsilon} \|W^{\epsilon}(t)\|_{L^2(\Omega)}^2,
$$
i.e.

$$ \|W^{\epsilon}(t)\|_{L^2(\Omega)}^2 \leq \epsilon_{\text{in}}^2 e^{-\frac{1}{2\epsilon}t} + C(T)(\epsilon^2 + \epsilon \epsilon_{\text{in}}^2), \quad t \in (0, T). \quad (9.42) $$

Finally, we plug (9.41), (9.42) into (9.37) to obtain

$$ \begin{aligned} & \frac{d}{dt} \mathcal{L}(t) - c_1 \|\nabla U^{\epsilon}(t)\|_{L^2(\Omega)}^2 - c_2 \|\nabla V^{\epsilon}(t)\|_{L^2(\Omega)}^2 - c_3 \|\nabla W^{\epsilon}(t)\|_{L^2(\Omega)}^2 \\ & \leq \epsilon_{\text{in}}^2 e^{-\frac{1}{2\epsilon}t} + C(T)(\epsilon^2 + \epsilon \epsilon_{\text{in}}^2). \end{aligned} \quad (9.43) $$

We are now able to conclude. Indeed, (2.31) follows by (9.41) and (9.43) integrated over $(0, t)$ and taking into account the positivity of $\mathcal{L}(t)$, the negativity of $c_1, c_2, c_3$ and (9.40). (2.32) follows by (9.42), (9.43) and $u_a^{\epsilon} - u_a^{*}(u, v) = (u^{\epsilon} - u) - W^{\epsilon}$.

## A Proof of Theorem 2.1

Throughout the proof, we omit the $\epsilon$ superscript for the sake of clarity.

First, we truncate the reaction functions in (1.1) in order to obtain globally Lipschitz functions in $L^p(\Omega)$. To do this, fix $M > 0$. Let us denote $U = (u_a, u_b, v)$, $|U| = \max(|u_a|, |u_b|, |v|)$ and let

$$ f_{\nu}^{M}(U) := \begin{cases} f_{\nu}(U), & \text{if } |U| \leq M, \\ f_{\nu}\left(M\frac{U}{|U|}\right), & \text{if } |U| > M, \end{cases} $$

for $\nu = a, b, v$, where $f_\nu$ is defined by (1.4). Moreover, let

$$ \psi^M(x) := \begin{cases} A^\alpha, & \text{if } x < 0, \\ \psi(x), & \text{if } x \in [0, M], \\ \psi(M), & \text{if } x > M, \end{cases} \quad \phi^M(x) := \begin{cases} B^\beta, & \text{if } x < 0, \\ \phi(x), & \text{if } x \in [0, M], \\ \phi(M), & \text{if } x > M, \end{cases} $$

where $\psi, \phi$ are the transiction functions defined by (1.7).

Next, we let $\Lambda_a^M, \Lambda_b^M$ be the relative satisfaction measures defined by (1.6) with $\psi, \phi$ replaced by $\psi^M, \phi^M$ and also $Q^M(U) = \Lambda_b^M g^M(u_b) - \Lambda_a^M g^M(u_a)$, where

$$ g^M(x) := \begin{cases} x, & \text{if } |x| \leq M, \\ M \frac{x}{|x|}, & \text{if } |x| > M. \end{cases} $$

Finally, we define the nonlinear mapping $F^M : \mathbb{R}^3 \to \mathbb{R}^3$ as

$$ F^M(U) := (f_a^M(U) + \epsilon^{-1}Q^M(U), f_b^M(U) - \epsilon^{-1}Q^M(U), f_v^M(U)). $$

Note that the functions $f_\nu^M, \psi^M, \phi^M$ are globally Lipschitz and bounded. Therefore, $F^M$ is globally Lipschitz, bounded and

$$ F^M(U) := (f_a(U) + \epsilon^{-1}Q(U), f_b(U) - \epsilon^{-1}Q(U), f_v(U)), \quad \text{if } |U| \leq M. $$

Given $p \in (1, +\infty)$, we consider the operator $A_p$ on $X_p := (L^p(\Omega))^3$ defined by

$$ \begin{cases} D(A_p) = D_p^3 & (\text{see (2.1)}) \\ A_p U = (d_a \Delta u_a, d_b \Delta u_b, d_v \Delta v) & \text{for } U \in D(A_p), \end{cases} \quad (A.1) $$
and the abstract initial value problem

$$U'(t) = A_p U(t) + F^M(U(t)), \quad t > 0, \quad U(0) = U^{\text{in}} := (u_a^{\text{in}}, u_b^{\text{in}}, v^{\text{in}}). \quad (A.2)$$

We will solve (A.1), (A.2) and then we will get rid of the truncation. The main ingredient is that $A_p: D(A_p) \subset X_p \to X_p$ is a sectorial operator ([22], Theorem 3.1.3). Hence it generates in $X_p$ an analytic semigroup denoted $(e^{tA_p})_{t \ge 0}$ ([22], Chapter 2). Moreover, $A_p$ is closed so that $D(A_p)$, endowed with the graph norm, is a Banach space. $D(A_p)$ being also dense in $X_p$, the semigroup is strongly continuous, i.e. $\lim_{t \to 0} e^{tA_p} U = U$, for all $U \in X_p$. Furthermore, there exists $K_p > 0$ and $\omega_p \in \mathbb{R}$ such that (see [22], Proposition 2.1.1)

$$\|e^{tA}\|_{L(X_p)} \le K_p e^{\omega_p t}, \quad \forall t \ge 0. \quad (A.3)$$

*First step: well-posedness of (A.1), (A.2).* Let $\| \cdot \|_p$ denote the usual norm in $X_p$. We start proving that (A.1), (A.2) has a unique mild solution, i.e. a unique function $U \in C^0([0, \infty), X_p)$ such that

$$U(t) = e^{tA}U^{\text{in}} + \int_{0}^{t} e^{(t-s)A} F^{M}(U(s)) ds, \quad \forall t \ge 0. \quad (A.4)$$

It is easily seen that $F^M$ maps $X_p$ into $X_p$ and

$$\|F^M(U)\|_p \le \|F^M(0)\|_p + L_M \|U\|_p = L_M \|U\|_p, \quad \forall U \in X_p. \quad (A.5)$$

Therefore, (A.4) makes sense since, by assumption (2.3), $U^{\text{in}} \in X_p$ and, for all $U \in C^0([0, \infty), X_p)$ and all $t > 0$, $F^M(U(\cdot)) \in L^1((0, t); X_p)$. Moreover, the Lipschitz property of $F^M$ together with (A.3) and Gronwall's Lemma gives us the uniqueness of (A.4). The same ingredients give us the continuous dependence of $U$ with respect to $U^{\text{in}}$. Therefore, it remains to prove the existence of $U$ and that $U$ belongs to $C^1([0, \infty); X_p) \cap C^0([0, \infty); D_p(A))$ for all $p \in (1, +\infty)$.

Let $\theta > 0$ be such that $\omega_p + \theta > 0$. The existence is proved using the contraction mapping principle in the space

$$E := \{U \in C^0([0, \infty), X_p) : \|U\|_E = \sup_{t \ge 0} e^{-(\omega_p + \theta)t} \|U(t)\|_p < \infty\},$$

that is a Banach space when endowed with the norm $\|U\|_E$. Hence, given $U \in E$, we set

$$\Phi(U)(t) = e^{tA}U^{\text{in}} + \int_{0}^{t} e^{(t-s)A} F^{M}(U(s)) ds, \quad \forall t \ge 0.$$

We claim that $\Phi$ maps $E$ into $E$ and it is a contraction provided that $\theta > K_p L_M$. Indeed, it is clear that $\Phi(U) \in C([0, \infty), X_p)$. Moreover, using (A.3) and (A.5) or the Lipschitz property of $F^M$, we obtain that, for all $U, V \in E$,

$$\|\Phi(U)\|_E \le K_p \|U^{\text{in}}\|_p + K_p \theta^{-1} L_M \|U\|_E,$$

and

$$\|\Phi(U) - \Phi(V)\|_E \le K_p \theta^{-1} L_M \|U - V\|_E.$$

Now, for all $T \in (0, \infty)$, the mild solution $U$ is Lipschitz $[0, T] \to X_p$. Indeed, on the one hand, using (A.3), the Lipschitz property of $F^M$ and Gronwall's Lemma again, we have

$$\|U(t)\|_p \le K e^{(K_p L_M + \omega_p)t} \|U^{\text{in}}\|_p, \quad \forall t \ge 0. \quad (A.6)$$
On the other hand, as $U(t+h)$, $t, h > 0$, is a mild solution of (A.2) with initial data $U(h) \in X_p$, proceeding as above, we obtain

$$ \|U(t+h) - U(t)\|_{p} \le K e^{(K_p L_M + \omega_p)t} \|U(h) - U^{\text{in}}\|_{p}, \quad \forall t \ge 0. \quad (A.7) $$

Next, by [22] Propositions 2.1.1, 2.1.4, and since $U^{\text{in}} \in D_p(A)$ (see (2.3)), it holds

$$ e^{hA}U^{\text{in}} - U^{\text{in}} = A \int_{0}^{h} e^{sA}U^{\text{in}} ds = \int_{0}^{h} A e^{sA}U^{\text{in}} ds = \int_{0}^{h} e^{sA}AU^{\text{in}} ds. $$

Therefore, from (A.3)-(A.5) and the above equality, we get

$$ \|U(h) - U^{\text{in}}\|_{p} \le K_{p}e^{\omega_{p}h}h \left( \|AU^{\text{in}}\|_{p} + L_{M} \sup_{0 \le s \le h} \|U(s)\|_{p} \right), \quad (A.8) $$

and the Lipschitz property follows from (A.6), (A.7), (A.8).

As a consequence, $F^M(U(\cdot))$ is Lipschitz $[0, T] \to X_p$. Therefore, taking also into account that $F^M$ is bounded, and applying Theorem 4.3.1 (ii) and Lemma 4.1.6 in [22], we have that $U \in C^1([0, \infty); X_p) \cap C^0([0, \infty); D(A_p))$.

Finally, note that, since $\Omega$ is bounded, if $p \ge q$, $X_p \subset X_q$ (with continuous embedding), $D(A_p) \subset D(A_q)$ and $A_p U = A_q U$ if $U \in D(A_p)$. Therefore, by assumption (2.3), the above time regularity holds true for all $p \in (1, +\infty)$. We can drop the subscript $p$ in the sequel.

**Second step: well-posedness of (1.1)-(1.7).** Let $U = (u_a, u_b, v)$ be the unique solution of (A.1), (A.2). Then, $u_a, u_b > 0$ and $v \ge 0$ on $(0, +\infty) \times \Omega$, since the semi-group $(e^{tA})_{t \ge 0}$ is strongly positive, the initial data are non-negative, with $u_a^{\text{in}}, u_b^{\text{in}}$ not identically zero, and the nonlinear mapping $F^M$ is quasi-positive (see [16]).

Multiplying the equation for $u_a$ in (A.2) by $u_a^{p-1}$, $p > 1$, using the positivity of the solution in $f_a^M$ and $Q^M$, the fact that the satisfaction measure $\Lambda_b^M$ lies in $[0, 1)$ and the Young inequality, we get

$$ \begin{aligned} \frac{1}{p} \frac{d}{dt} \int_{\Omega} u_a^p dx &\le \eta_a \int_{\Omega} u_a^p dx + \frac{1}{\varepsilon} \int_{\Omega} u_b u_a^{p-1} dx \\ &\le \bar{\eta} \int_{\Omega} u_a^p dx + \frac{1}{\varepsilon p} \int_{\Omega} u_b^p dx + \frac{1}{\varepsilon} \left(1 - \frac{1}{p}\right) \int_{\Omega} u_a^p dx. \end{aligned} $$

Similarly, for $u_b$ we have

$$ \frac{1}{p} \frac{d}{dt} \int_{\Omega} u_b^p dx \le \bar{\eta} \int_{\Omega} u_b^p dx + \frac{1}{\varepsilon p} \int_{\Omega} u_a^p dx + \frac{1}{\varepsilon} \left(1 - \frac{1}{p}\right) \int_{\Omega} u_b^p dx, $$

so that, for all $t > 0$,

$$ \frac{d}{dt} (\|u_a(t)\|_{L^p(\Omega)}^p + \|u_b(t)\|_{L^p(\Omega)}^p) \le p(\bar{\eta} + \frac{1}{\varepsilon})(\|u_a(t)\|_{L^p(\Omega)}^p + \|u_b(t)\|_{L^p(\Omega)}^p). $$

Integrating the above differential inequality over $(0, T)$, we end up with

$$ \sup_{0 \le t \le T} (\|u_a\|_{L^p(\Omega)}^p + \|u_b\|_{L^p(\Omega)}^p) \le (\|u_a^{\text{in}}\|_{L^p(\Omega)}^p + \|u_b^{\text{in}}\|_{L^p(\Omega)}^p)^{\frac{1}{p}} e^{(\bar{\eta} + \frac{1}{\varepsilon})T}, $$

implying, as $p \to \infty$ and for all $T > 0$,

$$ \|u_a\|_{L^\infty(0,T;L^\infty(\Omega))}, \|u_b\|_{L^\infty(0,T;L^\infty(\Omega))} \le (\|u_a^{\text{in}}\|_{L^\infty(\Omega)} \vee \|u_b^{\text{in}}\|_{L^\infty(\Omega)}) e^{(\bar{\eta} + \frac{1}{\varepsilon})T} := M_1. $$
Similarly, we prove that

$$ \|v\|_{L^{\infty}(0,T;L^{\infty}(\Omega))} \le \|v^{\text{in}}\|_{L^{\infty}(\Omega)} e^{\eta_v T} := M_2. $$

We now fix $T > 0$ and we let $M_T = M_1 \cup M_2$. It follows from the above estimates that $|U(t,x)| \le M_T$, for all $(t,x) \in [0,T] \times \Omega$. Therefore, $F^{M_T}(U) := (f_a(U) + \varepsilon^{-1}Q(U), f_b(U) - \varepsilon^{-1}Q(U), f_v(U))$, on $[0,T] \times \Omega$. Thus, we see that $U$ satisfies (1.1) on $[0,T] \times \Omega$. We denote this solution by $U^T$, since at this stage it might depend on $T$. We claim that in fact it does not. Indeed, let $0 < T < S$. On $[0,T]$, both $U^T$ and $U^S$ are bounded. Therefore, both $U^T$ and $U^S$ satisfy the same equation (A.2) provided $M$ is chosen sufficiently large. By uniqueness for the equation (A.2), it follows that $U^T = U^S$ on $[0,T] \times \Omega$. Thus we obtain a solution $U$ of (1.1) on $[0, \infty) \times \Omega$.

Third step: estimates uniform in $\varepsilon$. It is easily seen from (1.1), (1.4) (using the positivity of the solution and (2.2)) that $u_a + u_b$ satisfies

$$ \frac{d}{dt} \|u_a + u_b\|_{L^1(\Omega)} \le \bar{\eta} \|u_a + u_b\|_{L^1(\Omega)} - a\eta_a \|u_a\|_{L^2(\Omega)}^2 - b\eta_b \|u_b\|_{L^2(\Omega)}^2, \quad (A.9) $$

so that $y(t) = \|(u_a + u_b)(t)\|_{L^1(\Omega)}$ satisfies

$$ y'(t) \le \bar{\eta} y(t) - \frac{\eta}{2|\Omega|} y^2(t), \quad \forall t \ge 0. $$

Integrating the above differential inequality, we obtain (2.5). Furthermore, integrating (A.9) over $(0,T)$, we obtain

$$ a\eta_a \|u_a\|_{L^2(0,T;L^2(\Omega))}^2 + b\eta_b \|u_b\|_{L^2(0,T;L^2(\Omega))}^2 \le \|u_a^{\text{in}} + u_b^{\text{in}}\|_{L^1(\Omega)} + \bar{\eta} \int_0^T y(t) dt. $$

Using (2.5) in the above inequality, we get (2.7).

Next, multiplying the equation for $v$ in (1.1) by $v^{p-1}$, $p > 1$, and using again the positivity of the solution in $f_v$, we get

$$ \frac{1}{p} \frac{d}{dt} \int_{\Omega} v^p dx \le \eta_v \int_{\Omega} v^p dx - r_v \int_{\Omega} v^{p+1} dx. \quad (A.10) $$

Plugging into (A.10) the Hölder inequality

$$ \int_{\Omega} v^p dx \le \left( \int_{\Omega} v^{p+1} dx \right)^{\frac{p}{p+1}} |\Omega|^{\frac{1}{p+1}}, $$

we have

$$ \frac{1}{p} \frac{d}{dt} \int_{\Omega} v^p dx \le \eta_v \int_{\Omega} v^p dx - \frac{r_v}{|\Omega|^{1/p}} \left( \int_{\Omega} v^p dx \right)^{1+\frac{1}{p}}. \quad (A.11) $$

Integrating (A.11) over $(0,t)$, $t > 0$, we obtain

$$ \|v(t)\|_{L^p(\Omega)} \le \max \left\{ \|v^{\text{in}}\|_{L^p(\Omega)}, \frac{\eta_v}{r_v} |\Omega|^{\frac{1}{p}} \right\}, \quad (A.12) $$

implying (2.6) as $p \to \infty$.
Fourth step: maximal regularity and further estimates. Let $p = 2$. Then, $(e^{tA})_{t \ge 0}$ is a semigroup of contraction, i.e. $\|e^{tA}\Phi\|_{X_2} \le \|\Phi\|_{X_2}$, for all $\Phi \in X_2$ and all $t \ge 0$ ([5], Theorem 3.1.1). As $-A \ge 0$ and $A$ is self-adjoint, Theorem 1.3.9 in [10] gives

$$ \|e^{tA}\Phi\|_{X_p} \le \|\Phi\|_{X_p}, \quad \forall \Phi \in X_p, p \in [1, +\infty], t \ge 0. $$

It follows from [19] Theorem 1 applied to the component $v$ of the strict solution $(u_a, u_b, v)$ that, for all $T > 0$ and all $p \in (1, +\infty)$, there exists $C_p^{\text{MR}} > 0$ (not depending on $T$) such that

$$ \|\partial_t v\|_{L^p(\Omega_T)} + \|\Delta v\|_{L^p(\Omega_T)} \le C_p^{\text{MR}} (\|\Delta v^{\text{in}}\|_{L^p(\Omega)} + \|f_v(u_a, u_b, v)\|_{L^p(\Omega_T)}). \quad (A.13) $$

Next, by the classical Agmon-Douglis-Nirenberg a priori estimates (see [22], Theorem 3.1.1), there exists $C_p^{\text{ADN}} > 0$ such that, for all $t \ge 0$,

$$ \|v(t)\|_{W^{2,p}(\Omega)} \le C_p^{\text{ADN}} (\|v(t)\|_{L^p(\Omega)} + \|\Delta v(t)\|_{L^p(\Omega)}), \quad (A.14) $$

while, by the Gagliardo-Nirenberg inequality [24], there exists $C_p^{\text{GN}} > 0$ such that, for all $t \ge 0$ and all $i = 1, \dots, N$,

$$ \|\partial_i v(t)\|_{L^{2p}(\Omega)} \le C_p^{\text{GN}} \left( \max_{1 \le i,j \le N} \|\partial_{i,j} v(t)\|_{L^p(\Omega)} \right)^{\frac{1}{2}} \|v(t)\|_{L^\infty(\Omega)}^{\frac{1}{2}} + C_p^{\text{GN}} \|v(t)\|_{L^\infty(\Omega)}. \quad (A.15) $$

Finally, by (2.6), there exists $C(K_\infty, |\Omega|) > 0$, such that

$$ \|f_v(u_a, u_b, v)\|_{L^p(\Omega_T)} \le C(K_\infty, |\Omega|) (\|u_a\|_{L^p(\Omega_T)} + \|u_b\|_{L^p(\Omega_T)} + T^{\frac{1}{p}}). \quad (A.16) $$

Hence, plugging (A.16) into (A.13), and combining the resulting inequality with (A.14) integrated over $(0, T)$, and (2.6), we obtain that there exists a constant $C_1(C_p^{\text{ADN}}, C_p^{\text{MR}}, K_\infty, |\Omega|) > 0$ such that

$$ \|\partial_t v\|_{L^p(\Omega_T)} + \sum_{i,j} \|\partial_{ij} v\|_{L^p(\Omega_T)} \le C_1 (\|\Delta v^{\text{in}}\|_{L^p(\Omega)} + \|u_a\|_{L^p(\Omega_T)} + \|u_b\|_{L^p(\Omega_T)} + T^{\frac{1}{p}}). \quad (A.17) $$

Finally, combining the Gagliardo-Nirenberg inequality above integrated over $(0, T)$ and (A.17), give us the existence of $C_2(C_p^{\text{GN}}, C_1, K_\infty, N)$ such that

$$ \|\nabla v\|_{L^{2p}(\Omega_T)}^{2p} \le C_2^{2p} (\|\Delta v^{\text{in}}\|_{L^p(\Omega)}^p + \|u_a\|_{L^p(\Omega_T)}^p + \|u_b\|_{L^p(\Omega_T)}^p + T). \quad (A.18) $$

(2.8) follows taking $p = 2$ in (A.17),(A.18) and using (2.7).

## B Proof of Lemma 2.3

By definitions (1.5),(1.7) and assumption (H1), $\Lambda(u_a, u_b, v) \ge A^\alpha > 0$, for all $(u_a, u_b, v) \in \mathbb{R}_+^3$, so that $Q(u_a, u_b, v) = 0$ if and only if $q(u_a, u_b, v) = 0$. Moreover, if $\tilde{u} = 0$, the unique nonnegative solution of (1.12) is $(u_a^*, u_b^*) = (0, 0)$, $\forall \tilde{v} \ge 0$.

Let us denote $\Sigma = \{(\tilde{u}, u_b, \tilde{v}) \in \mathbb{R}^3 : 0 \le u_b \le \tilde{u}, \tilde{u} \ge 0, \tilde{v} \ge 0\}$ and

$$ \tilde{q}(\tilde{u}, u_b, \tilde{v}) = q(\tilde{u} - u_b, u_b, \tilde{v}), \quad (\tilde{u}, u_b, \tilde{v}) \in \Sigma. $$
By (1.5), (1.7) and (H1) again, it is easily seen that $\tilde{q}$ is a continuous differentiable function such that, for all $u_b \in (0, \tilde{u}), \tilde{u} > 0, \tilde{v} \ge 0$,

$$ \partial_2 \tilde{q}(\tilde{u}, u_b, \tilde{v}) = -\partial_1 q(\tilde{u} - u_b, u_b, \tilde{v}) + \partial_2 q(\tilde{u} - u_b, u_b, \tilde{v}) \ge A^\alpha > 0 $$

and $\tilde{q}(\tilde{u}, 0, \tilde{v}) = q(\tilde{u}, 0, \tilde{v}) < 0$, $\tilde{q}(\tilde{u}, \tilde{u}, \tilde{v}) = q(0, \tilde{u}, \tilde{v}) > 0$. Therefore, $\forall \tilde{u} > 0, \tilde{v} \ge 0$, there exists a unique $U_b = U_b(\tilde{u}, \tilde{v})$ such that $U_b \in (0, \tilde{u}), \tilde{q}(\tilde{u}, U_b, \tilde{v}) = 0$, i.e. for all $\tilde{u} > 0, \tilde{v} > 0$, there exists a unique solution of (1.12) given by $(u_a^*(\tilde{u}, \tilde{v}), u_b^*(\tilde{u}, \tilde{v})) = (\tilde{u} - U_b, U_b)$.

As it holds $\partial_2 \tilde{q}(U_b, \tilde{u}, \tilde{v}) > 0$, by the implicit function theorem, for all $(\tilde{u}, \tilde{v}) \in (0, +\infty)^2$, there exists a neighbourhood $W$ of $(\tilde{u}, \tilde{v})$ and a unique continuously differentiable map $u_b^* : W \to \mathbb{R}_+$ such that, $\forall (\tilde{u}, \tilde{v}) \in W$, $u_b^*(\tilde{u}, \tilde{v}) = U_b$ and $\tilde{q}(\tilde{u}, u_b^*(\tilde{u}, \tilde{v}), \tilde{v}) = 0$. Furthermore, it is easily seen that $u_b^*$ is defined and continuously differentiable over $(0, +\infty)^2$. Hence, defining $u_b^*(0, \tilde{v}) = 0$ and

$$ u_a^*(\tilde{u}, \tilde{v}) = \tilde{u} - u_b^*(\tilde{u}, \tilde{v}), \quad (\tilde{u}, \tilde{v}) \in [0, +\infty)^2, $$

we have that the pair $(u_a^*(\tilde{u}, \tilde{v}), u_b^*(\tilde{u}, \tilde{v}))$ is the unique solution of (1.12).

Finally, differentiating the identities below with respect to $\tilde{u}$ and $\tilde{v}$, we obtain

$$ \tilde{q}(\tilde{u}, u_b^*(\tilde{u}, \tilde{v}), \tilde{v}) = q(\tilde{u} - u_b^*(\tilde{u}, \tilde{v}), u_b^*(\tilde{u}, \tilde{v}), \tilde{v}) = 0 \quad \text{and} \quad u_a^*(\tilde{u}, \tilde{v}) = \tilde{u} - u_b^*(\tilde{u}, \tilde{v}). $$

$$ \partial_{\tilde{u}} u_b^*(\tilde{u}, \tilde{v}) = \frac{\partial_1 q(u_a^*(\tilde{u}, \tilde{v}), u_b^*(\tilde{u}, \tilde{v}), \tilde{v})}{\partial_1 q(u_a^*(\tilde{u}, \tilde{v}), u_b^*(\tilde{u}, \tilde{v}), \tilde{v}) - \partial_2 q(u_a^*(\tilde{u}, \tilde{v}), u_b^*(\tilde{u}, \tilde{v}), \tilde{v})} \quad (B.1) $$

$$ \partial_{\tilde{v}} u_b^*(\tilde{u}, \tilde{v}) = \frac{\partial_3 q(u_a^*(\tilde{u}, \tilde{v}), u_b^*(\tilde{u}, \tilde{v}), \tilde{v})}{\partial_1 q(u_a^*(\tilde{u}, \tilde{v}), u_b^*(\tilde{u}, \tilde{v}), \tilde{v}) - \partial_2 q(u_a^*(\tilde{u}, \tilde{v}), u_b^*(\tilde{u}, \tilde{v}), \tilde{v})} $$

and

$$ \partial_{\tilde{u}} u_a^*(\tilde{u}, \tilde{v}) = 1 - \partial_{\tilde{u}} u_b^*(\tilde{u}, \tilde{v}), \quad \partial_{\tilde{v}} u_a^*(\tilde{u}, \tilde{v}) = -\partial_{\tilde{v}} u_b^*(\tilde{u}, \tilde{v}). $$

Therefore, (2.21), (2.22) follow taking into account that (see (1.5))

$$ \begin{aligned} \partial_1 q(u_a, u_b, v) &= -\psi(b_{\text{fast}} u_b + d_{\text{fast}} v) - a_{\text{fast}} u_a \psi'(a_{\text{fast}} u_a + c_{\text{fast}} v) \\ \partial_2 q(u_a, u_b, v) &= \phi(b_{\text{fast}} u_b + d_{\text{fast}} v) + b_{\text{fast}} u_b \phi'(b_{\text{fast}} u_b + d_{\text{fast}} v) \end{aligned} \quad (B.2) $$

$$ \partial_3 q(u_a, u_b, v) = d_{\text{fast}} u_b \phi'(b_{\text{fast}} u_b + d_{\text{fast}} v) - c_{\text{fast}} u_a \psi'(a_{\text{fast}} u_a + c_{\text{fast}} v) $$

and the positivity of $\psi, \phi, \psi', \phi'$.

## References

[1] J.-P. AUBIN, *Compact sets in the space Lp(0, T; B)*, Ann. Mat. Pura Appl. (4), 146 (1987), pp. 65–96.

[2] M. BREDEN, C. KUEHN, AND C. SORESINA, *On the influence of cross-diffusion in pattern formation*, Journal of Computational Dynamics, 8 (2021), pp. 213–240.

[3] E. BROCCCHIERI, L. CORRIAS, H. DIETERT, AND Y.-J. KIM, *Evolution of dietary diversity and a starvation driven cross-diffusion system as its singular limit*, J. Math. Biol., 83 (2021), pp. Paper No. 58, 40.

[4] E. BROCCCHIERI, L. DESVILLETTE, AND H. DIETERT, *Study of a class of triangular starvation driven cross-diffusion systems*, Ricerche di Matematica, (2024), pp. 1–27.
[5] T. CAZENAVE AND A. HARAUX, *An introduction to semilinear evolution equations*, vol. 13 of Oxford Lecture Series in Mathematics and its Applications, The Clarendon Press, Oxford University Press, New York, 1998.

[6] L. CHEN AND A. JÜNGEL, *Analysis of a multidimensional parabolic population model with strong cross-diffusion*, SIAM Journal on Mathematical Analysis, 36 (2004), pp. 301–322.

[7] L. CHEN AND A. JÜNGEL, *Analysis of a parabolic cross-diffusion population model without self-diffusion*, Journal of Differential Equations, 224 (2006), pp. 39–59.

[8] F. CONFORTO, L. DESVILLETTE'S, AND C. SORESINA, *About reaction-diffusion systems involving the Holling-type II and the Beddington-DeAngelis functional responses for predator-prey models*, NoDEA Nonlinear Differential Equations Appl., 25 (2018), pp. Paper No. 24, 39.

[9] E. S. DAUS, L. DESVILLETTE'S, AND H. DIETERT, *About the entropic structure of detailed balanced multi-species cross-diffusion equations*, Journal of Differential Equations, 266 (2019), pp. 3861–3882.

[10] E. B. DAVIES, *Heat kernels and spectral theory*, vol. 92 of Cambridge Tracts in Mathematics, Cambridge University Press, Cambridge, 1989.

[11] L. DESVILLETTE'S, C. KUEHN, J.-E. SULZBACH, B. Q. TANG, AND B.-N. TRAN, *Slow manifolds for pde with fast reactions and small cross diffusion*, arXiv preprint arXiv:2501.16775, (2025).

[12] L. DESVILLETTE'S, T. LEPOUTRE, AND A. MOUSSA, *Entropy, duality, and cross diffusion*, SIAM Journal on Mathematical Analysis, 46 (2014), pp. 820–853.

[13] L. DESVILLETTE'S, T. LEPOUTRE, A. MOUSSA, AND A. TRESCASES, *On the entropic structure of reaction-cross diffusion systems*, Comm. Partial Differential Equations, 40 (2015), pp. 1705–1747.

[14] L. DESVILLETTE'S AND A. TRESCASES, *New results for triangular reaction cross diffusion system*, J. Math. Anal. Appl., 430 (2015), pp. 32–59.

[15] G. GAMBINO, M. LOMBARDO, AND M. SAMMARTINO, *Turing instability and traveling fronts for a nonlinear reaction-diffusion system with cross-diffusion*, Mathematics and Computers in Simulation, 82 (2012), pp. 1112–1132. Non-linear Waves: Computation and Theory-IX, WAVES 2009.

[16] A. HARAUX, *A simple characterization of positivity preserving semi-linear parabolic systems*, J. Korean Math. Soc., 54 (2017), pp. 1817–1828.

[17] M. IIDA, H. IZUHARA, AND R. KON, *Cross-diffusion predator-prey model derived from the dichotomy between two behavioral predator states*, Discrete Contin. Dyn. Syst. Ser. B, 28 (2023), pp. 6159–6178.

[18] M. IIDA, M. MIMURA, AND H. NINOMIYA, *Diffusion, cross-diffusion and competitive interaction*, J. Math. Biol., 53 (2006), pp. 617–641.

[19] D. LAMBERTON, *équations d'évolution linéaires associées à des semi-groupes de contractions dans les espaces L<sup>p</sup>*, J. Funct. Anal., 72 (1987), pp. 252–262.
[20] T. LEPOUTRE AND A. MOUSSA, *Entropic structure and duality for multiple species cross-diffusion systems*, Nonlinear Analysis, 159 (2017), pp. 298-315.

[21] J.-L. LIONS AND E. MAGENES, *Non-homogeneous boundary value problems and applications. Vol. II*, vol. Band 182 of Die Grundlehren der mathematischen Wissenschaften, Springer-Verlag, New York-Heidelberg, 1972. Translated from the French by P. Kenneth.

[22] A. LUNARDI, *Analytic semigroups and optimal regularity in parabolic problems*, vol. 16 of Progress in Nonlinear Differential Equations and their Applications, Birkhäuser Verlag, Basel, 1995.

[23] A. MOUSSA, B. PERTHAME, AND D. SALORT, *Backward parabolicity, cross-diffusion and Turing instability*, J. Nonlinear Sci., 29 (2019), pp. 139-162.

[24] L. NIRENBERG, *An extended interpolation inequality*, Ann. Scuola Norm. Sup. Pisa Cl. Sci., 20 (1966), pp. 733-737.

[25] N. SHIGESADA, K. KAWASAKI, AND E. TERAMOTO, *Spatial segregation of interacting species*, J. Theoret. Biol., 79 (1979), pp. 83-99.

[26] C. SORESINA, Q. B. TANG, AND B. N. TRAN, *Fast-reaction limits for predator-prey reaction-diffusion systems: improved convergence*, arXiv.

[27] B. Q. TANG AND B.-N. TRAN, *Rigorous derivation of michaelis-menten kinetics in the presence of slow diffusion*, SIAM Journal on Mathematical Analysis, 56 (2024), pp. 5995-6024.

[28] J. L. VÁZQUEZ, *An introduction to the mathematical theory of the porous medium equation*, in Shape optimization and free boundaries (Montreal, PQ, 1990), vol. 380 of NATO Adv. Sci. Inst. Ser. C: Math. Phys. Sci., Kluwer Acad. Publ., Dordrecht, 1992, pp. 347-389.

Email addresses:

Elisabetta Brocchieri: elisabetta.brocchieri@uni-graz.at

Lucilla Corrias: lucilla.corrias@univ-evry.fr

<footnote n="1">¹ Department of Mathematics and Scientific Computing, University of Graz, 8010 Graz, Austria.</footnote>

<footnote n="2">² Université Paris-Saclay, CNRS, Univ. Evry, Laboratoire de Mathématiques et Modélisation d'Evry, 91037, Evry-Courcironnes, France.</footnote>