# Entropy, Duality and Cross Diffusion

Laurent Desvillettes, Thomas Lepoutre, Ayman Moussa

► To cite this version:

Laurent Desvillettes, Thomas Lepoutre, Ayman Moussa. Entropy, Duality and Cross Diffusion. SIAM Journal on Mathematical Analysis, 2014, 46 (1), pp.820-853. 10.1137/130908701. hal-00785379

HAL Id: hal-00785379

https://inria.hal.science/hal-00785379v1

Submitted on 22 Feb 2026

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_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Entropy, duality and Cross Diffusion

L. DESVILLETTES

CMLA, ENS Cachan & CNRS

61 Av. du Pdt. Wilson, 94235 Cachan Cedex, France

E-mail: desville@cmla.ens-cachan.fr

TH. LEPOUTRE

INRIA Rhône Alpes

Project-team DRACULA

Université de Lyon

CNRS UMR 5208

Université Lyon 1

Institut Camille Jordan

43 blvd. du 11 novembre 1918

F-69622 Villeurbanne cedex France

A. MOUSSA

UPMC Université Paris 06 & CNRS

UMR 7598, LJLL, F-75005, Paris, France

& INRIA Paris-Rocquencourt, équipe-Projet REO

BP 105, F-78153 Le Chesnay Cedex, France

E-mail : moussa@ann.jussieu.fr

November 15, 2013

Abstract

This paper is devoted to the use of the entropy and duality methods for the ex-
istence theory of reaction-cross diffusion systems consisting of two equations, in any
dimension of space. Those systems appear in population dynamics when the diffu-
sion rates of individuals of two species depend on the concentration of individuals
of the same species (self diffusion), or of the other species (cross diffusion).

# 1 Introduction

We are interested in populations consisting of individuals belonging to two distinct species (or two classes
in the same species) of (typically) animals, and interacting through diffusion and competition.

The unknowns are the concentration (number density) of individuals of the first species $u_1 := u_1(t, x) \ge 0$ and of the second species $u_2 := u_2(t, x) \ge 0$.

In absence of competition, the respective populations $u_1$ and $u_2$ would grow at a (positive) rate $r_1$ and $r_2$.

The competition is taken into account through logistic-type terms, in such a way that the growth
rate becomes $r_1 - s_{11}(u_1) - s_{12}(u_2)$ for the first species, and $r_2 - s_{22}(u_2) - s_{21}(u_1)$ for the second one.
The (nonnegative) terms $s_{11}(u_1)$ and $s_{22}(u_2)$ are called intraspecific competition, while the (nonnegative)
terms $s_{12}(u_2)$ and $s_{21}(u_1)$ are by definition the interspecific competition. Since we are interested here
in situations when there is no cooperation (or predator-prey type interaction), we shall assume in the
sequel that all functions $s_{ij}$ are nonnegative and increasing (see [26], p.94).

We also assume that the individuals of the two species diffuse with a rate that depends on the
number density of both species: we denote $d_{11}(u_1) + a_{12}(u_2)$ the diffusion rate of the first species, and
$d_{22}(u_2) + a_{21}(u_1)$ the diffusion rate of the second species. The diffusion terms related to $d_{11}(u_1)$ and
$d_{22}(u_2)$ are called self diffusion, those related to $a_{12}(u_2)$ and $a_{21}(u_1)$ are called cross diffusion. We
are interested in this paper in the case when all the functions $d_{ii}$, $a_{ij}$ are nondecreasing (that is, all individuals try to leave the points where the competition is highest). This type of model was introduced in [28].

We are led to write down the following system (on $\mathbb{R}_+ \times \Omega$, where $\Omega$ is a smooth bounded open subset of $\mathbb{R}^M$ with outward unit normal vector $n(x)$ at point $x \in \partial\Omega$):

$$ \partial_t u_1 - \Delta [(d_{11}(u_1) + a_{12}(u_2)) u_1] = u_1(r_1 - s_{11}(u_1) - s_{12}(u_2)) \quad (1) $$

$$ \partial_t u_2 - \Delta [(d_{22}(u_2) + a_{21}(u_1)) u_2] = u_2(r_2 - s_{22}(u_2) - s_{21}(u_1)) \quad (2) $$

with the homogeneous Neumann boundary conditions

$$ \nabla u_1 \cdot n(x) = \nabla u_2 \cdot n(x) = 0 \quad \text{on} \quad \mathbb{R}_+ \times \partial\Omega, \quad (3) $$

and the initial condition

$$ u_1(0, \cdot) = u_1^{in} \geq 0, \quad u_2(0, \cdot) = u_2^{in} \geq 0. \quad (4) $$

The case treated in this paper corresponds to a situation where the reaction terms are strictly sub-quadratic or dominated by the self diffusion, and where the cross diffusions are subquadratic.

We detail below the set of mathematical assumptions that will be imposed on the coefficients $r_i$, $d_{ii}$, $a_{ij}$, $s_{ij}$, and which correspond to the case described above:

**H1** For $i, j \in \{1, 2\}$, $r_i \in \mathbb{R}_+$ and $s_{ij}$ is a nonnegative continuous function on $\mathbb{R}_+$ which is either strictly sublinear: $\lim_{x \to +\infty} \frac{s_{ij}(z)}{z} = 0$; or dominated by the self diffusion in the following sense: $\lim_{z \to +\infty} \frac{s_{ii}(z)}{z d_{ii}(z)} = 0$ and (for $i \neq j$) $\lim_{z \to +\infty} \frac{s_{ij}(z)}{z \sqrt{d_{jj}(z) + a_{ij}(z)}} = 0$;

**H2** For $i \neq j \in \{1, 2\}$, $a_{ij}$ is continuous on $\mathbb{R}_+$ and belongs to $\mathcal{C}^2([0, +\infty])$. It is also nonnegative, nondecreasing, concave and vanishes at point 0. Furthermore, there exists $\alpha \in [0, 1[$, $C > 0$, such that

$$ \forall z \in [0, +\infty[, \quad z^{\alpha} a'_{ij}(z) \geq C. $$

**H3** The self diffusion rate $d_{ii}$ is continuous on $\mathbb{R}_+$ and belongs to $\mathcal{C}^1([0, +\infty[)$. It is also nonnegative, nondecreasing and such that $d_{ii}(0) > 0$ (In the proof of our main Theorem, we use $d_{ii} \geq 1$ for the sake of readability).

The union of the assumptions **H1**, **H2** and **H3** on the parameters $[\text{for } i, j \in \{1, 2\}]$ $r_i$, $d_{ii}$, $s_{ij}$, $a_{ij}$ will be called in the sequel the **H** assumptions.

Since we consider the homogeneous Neumann boundary condition, it is useful to introduce the following notation:

**Definition 1.1.** for any space of functions defined on $\Omega$ whose gradient has a well-defined trace on the boundary $\partial\Omega$ (such as $\mathcal{C}^\infty(\bar{\Omega})$ or $\mathbf{H}^2(\Omega)$ for instance), we add the subscript $\nu$ (the former spaces hence become $\mathcal{C}_\nu^\infty(\bar{\Omega})$ or $\mathbf{H}_\nu^2(\Omega)$) to describe the subspace of functions satisfying the homogeneous Neumann boundary condition.

Our main Theorem reads:

**Theorem 1.2.** Let $\Omega$ be a smooth $(\mathcal{C}^2)$ bounded open subset of $\mathbb{R}^M$ ($M \geq 1$). Let $u^0 := (u_1^{in}, u_2^{in}) \in \mathbf{L}^2(\Omega)^2$ be a couple of nonnegative functions and assume that assumptions **H** are satisfied on the coefficients of the system.

Then, for any $T > 0$, there exists a couple $u := (u_1, u_2)$ of nonnegative weak solutions to (1) – (3) on $[0, T]$ in the following sense: for $i = 1, 2$ ($j \neq i$),

$$ \int_0^T \int_{\Omega} [d_{ii}(u_i(t, x)) + a_{ij}(u_j(t, x))] |u_i(t, x)|^2 dx dt \leq D_T(u^0), \quad (5) $$
and for any $\theta \in \mathcal{D}([0, T]; \mathcal{C}_{\nu}^{\infty}(\bar{\Omega}))$, we have the weak formulation

$$ 
\begin{aligned}
& - \int_{\Omega} u_1^{in}(x) \theta(0, x) dx - \int_{0}^{\infty} \int_{\Omega} u_1(t, x) \partial_t \theta(t, x) dx dt \\
& - \int_{0}^{\infty} \int_{\Omega} \Delta\theta(t, x) [d_{11}(u_1(t, x)) + a_{12}(u_2(t, x))] u_1(t, x) dx dt \\
& = \int_{0}^{\infty} \int_{\Omega} \theta(t, x) u_1(t, x) (r_1 - s_{11}(u_1(t, x)) - s_{12}(u_2(t, x))) dx dt,
\end{aligned}
$$

and

$$ 
\begin{aligned}
& - \int_{\Omega} u_2^{in}(x) \theta(0, x) dx - \int_{0}^{\infty} \int_{\Omega} u_2(t, x) \partial_t \theta(t, x) dx dt \\
& - \int_{0}^{\infty} \int_{\Omega} \Delta\theta(t, x) [d_{22}(u_2(t, x)) + a_{21}(u_1(t, x))] u_2(t, x) dx dt \\
& = \int_{0}^{\infty} \int_{\Omega} \theta(t, x) u_2(t, x) (r_2 - s_{22}(u_2(t, x)) - s_{21}(u_1(t, x))) dx dt.
\end{aligned}
$$

Moreover, denoting $Q_T := [0, T] \times \Omega$, the following bounds hold:

$$ \sup_{t \in [0, T]} \|u_i(t, \cdot)\|_{L^1(\Omega)} \le e^{r_i T} \|u_i^{in}\|_{L^1(\Omega)}, $$

$$ \left\| \frac{1}{u_1} a'_{21}(u_1) \nabla u_1 \right\|_{L^2(Q_T)} + \left\| \frac{1}{u_2} a'_{12}(u_2) \nabla u_2 \right\|_{L^2(Q_T)} \le K_T (1 + \|u^0\|_{L^2(\Omega)}), $$

$$ \left\| \nabla \sqrt{a_{21}(u_1)} \right\|_{L^2(Q_T)} + \left\| \nabla \sqrt{a_{12}(u_2)} \right\|_{L^2(Q_T)} + \left\| \nabla \sqrt{a_{21}(u_1) a_{12}(u_2)} \right\|_{L^2(Q_T)} \le K_T (1 + \|u^0\|_{L^2(\Omega)}). $$

The positive constants $K_T$ and $D_T(u^0)$ used above only depend on $T$, $\Omega$ and the data of the equations $(r_i, d_{ii}, a_{ij}, s_{ij})$, and the latter also depends on $\|u^0\|_{L^2(\Omega)}$.

Note that all the terms in the weak formulation above are well-defined thanks to assumptions **H**. Indeed, remembering that $d_{ii}$ is bounded below, estimate (5) implies $u_i \in L^2([0, T] \times \Omega)$ so that the reaction terms are integrable if sublinear (cf. **H1**), or, if dominated by self diffusion, thanks to (5) and the identity

$$ s_{12}(u_2) u_1 = \left( \frac{s_{12}(u_2)}{u_2 \sqrt{d_{22}(u_1) + a_{12}(u_2)}} \right) \left( \frac{u_2 \sqrt{d_{22}(u_2) + a_{12}(u_2)}}{\sqrt{a_{12}(u_2)}} \right) u_1 \sqrt{a_{12}(u_2)}. $$

The terms coming out of the cross diffusion are also well defined due to concavity assumptions **H2**. Those coming out of self diffusion are integrable thanks to (5).

We now comment the assumptions on the initial data and coefficients of system (1)–(4). The requirement that $u^0 \in (L^2(\Omega))^2$ can certainly be relaxed to $u^0 \in (H^{-1}(\Omega))^2$ together with $u^0 \log u^0 \in (L^1(\Omega))^2$. It is most likely possible to relax assumption **H1**, provided that **H2** and **H3** are reinforced. Finally, there is some hope of treating the special case when reaction terms are exactly quadratic (and there is no self diffusion), thanks to recent improvements in the theory of duality Lemmas [6].

We restricted in this paper our study to the case when both $a_{ij}$, $i \neq j$ are concave, whereas our methods should adapt in the case when one of them is concave and the other one is convex. The corresponding theory is then quite different and will be left to a future work. Note that when both $a_{ij}$, $i \neq j$ are convex, our feeling is that the existence theory is probably quite different (it is not clear if weak global solutions exist in general).

Our opinion is that for systems involving cross diffusion and consisting of more than two equations, the type of Lyapounov functional that we build in the sequel can exist only for a very small class of cross diffusion terms (that is, strong algebraic constraints have to be assumed on the cross diffusion coefficients).

Finally, we do not treat the case when one cross diffusion term is missing (sometimes called the “triangular case”), since the Lyapounov functional that we introduce in the sequel degenerates in this case.
In order to make our result easier to read, we provide in next section a complete description of the case when all functions appearing in the system $(d_{ii}, a_{ij}, s_{ij})$ are power laws.

Let us now describe how our work fits in the existing literature.

The question of local and global existence of classical solutions of systems like (1)–(4) was treated in many particular cases. Most of them deal either with the case when the system is in fact parabolic (that is, when the diffusion matrix is elliptic), which amounts to assume that self diffusion is dominant w.r.t. cross diffusion (cf. [30]), or in the triangular case (when $a_{21}(u_1) = 0$, cf. [24]). The question of local existence is usually treated using Amann’s theorem [1], and extension to global existence requires additional structure to preserve boundedness of solutions see [9, 14] for instance.

Our work is concerned with situations that are neither parabolic, nor triangular. It extends to a very large class of problems the result of [8], dedicated to the cross diffusion model for population dynamics introduced by Shigesada, Kawasaki and Teramoto in [28], where all coefficients ($d_{ii}$, $a_{ij}$, $s_{ij}$) are linear functions. Note that this system can be seen as a limiting case of the equations that we treat.

Our results rely on two main ingredients: entropy structure and duality Lemmas. We show that our systems possess a hidden entropy-like structure, strongly reminiscent of the entropy structure exhibited in [8]. In general, this structure however gives less estimates than in [8] because the cross-diffusion rates that we consider are (possibly strictly) concave. We therefore need another ingredient in order to recover existence of weak global solutions, namely a duality Lemma: we recall that duality Lemmas enable to recover $L^2$ type estimate for solutions to linear singular parabolic equations (with variable coefficients) when the diffusion rate is inside the Laplacian. This is how estimate (5) is derived. Our approximation procedure is a generalization of the one used in [8]. Let us mention that in [2], the authors used the same entropy structure as in [8], together with another (somewhat less involved) approximation procedure which unfortunately doesn’t seem adapted to our generalized setting.

For the use of an underlying entropy structure, its link with symmetrization, and its applications to existence of weak solutions, we refer to [22, 10]. The possibility to use such a structure in the case of cross diffusion was first noticed in [16], and exploited in [8, 7, 20, 18, 2]. The duality estimate that we use comes from [27], and was applied together with entropy methods in the framework of reaction-diffusion systems in [12]. Is was also applied in the framework of cross diffusion or similar models in [3, 4, 5].

We finally quote some works dealing with other aspects of cross diffusion models. For modeling issues, we refer to [19, 11, 4, 3]. For the analysis of equilibria, we cite [19, 23] for instance.

Unfortunately, as often in papers dealing with cross diffusion, the process of approximation enabling to make the estimates and structures rigorous is quite involved, and gives rise to various difficulties which explain the length and technicality of the proofs. As a consequence, we propose in next Section a complete and rigorous study of the *a priori* estimates and weak stability (thus avoiding the problems related to approximation) of a system which is representative of (1)–(4), but in which all the functions ($d_{ii}$, $a_{ij}$, $s_{ij}$) are power laws, in order to make explicit the limitations corresponding to assumptions **H**, and the use of both the entropy and duality Lemmas.

Our paper is structured as follows: in Section 2, we propose the study of *a priori* estimates and weak stability of the model with coefficients given by specific power laws. We also give at the end of this Section a road map for the (unfortunately quite technical) Proof of the main Theorem. In Section 3, we introduce notations which are used in the proof of our main Theorem, especially those related to the Lyapounov functional that we systematically use ; we also present some preliminary Lemmas used in the sequel. Then, Section 4 is devoted to the proof of existence of a solution to a finite-dimensional (Galerkin) approximation of a discrete time version of our (smoothed) system. *A priori* estimates (and their dependence with respect to the various approximations) are provided for such solutions. We let the dimension of the Galerkin approximation go to infinity in Section 5. The duality estimate is presented and proven in Section 6. The last section (that is, Section 7) is devoted to the relaxation of all remaining approximations (that is, the discrete time variable becomes continuous, and all regularizations are removed), which leads to the Proof of our main Theorem.
# 2 A priori estimates on an example and road map for the proof

## 2.1 A specific case of power law rates

In order to keep this Section completely self contained, we recall that in the sequel, $\Omega$ will always be a bounded subset of $\mathbb{R}^M$ with a $\mathcal{C}^2$ boundary. We fix a constant $T > 0$ and denote by $Q_T$ the set $[0, T] \times \Omega$.

In this Section, we also denote by $C(T, \Omega)$ a positive constant depending on $T$ and $\Omega$, but also on the $L^2$ norms of the initial data, and the coefficients of the system. It may change from line to line.

We consider a specific case of power law rates, though our approach applies to a considerably wider set of parameters (as soon as the $H$ assumptions are satisfied). For instance, we did not take into account the self-diffusion terms (which clearly tend to facilitate the study of our system), neither the mixed reaction terms (which may be handled thanks to the duality Lemma). Hopefully this particular case will appear less involved than a (too much) general presentation.

We start with a Proposition stating the a priori estimates for system (1) – (4).

**Proposition 2.1.** Consider $d, \beta > 0$.

Let $(u_1, u_2)$ be a strong (that is, belonging to $\mathcal{C}^2([0, T] \times \bar{\Omega})$) positive (that is, both $u_1$ and $u_2$ are positive on $Q_T$) solution to the reaction-cross diffusion equation:

$$ \partial_t u_1 - \Delta \left[ u_1 \left( d + \sqrt{u_2} \right) \right] = u_1 \left( 1 - u_1^\beta \right), \quad (6) $$

$$ \partial_t u_2 - \Delta \left[ u_2 \left( d + \sqrt{u_1} \right) \right] = u_2 \left( 1 - u_2^\beta \right), \quad (7) $$

with Neumann boundary conditions

$$ \forall x \in \partial\Omega, \qquad \nabla u_1 \cdot n(x) = 0, \qquad \nabla u_2 \cdot n(x) = 0. \quad (8) $$

Then the following a priori estimates hold for $i \neq j$:

$$ \sup_{t \in [0, T]} \|u_i(t, \cdot)\|_{L^1(\Omega)} \le e^T \|u_i(0, \cdot)\|_{L^1(\Omega)}, \quad (9) $$

$$ \int_0^T \int_{\Omega} u_i^{1+\beta} \le (1 + T e^T) \|u_i(0, \cdot)\|_{L^1(\Omega)}; \quad (10) $$

$$ \int_0^T \int_{\Omega} u_i^2 (d + \sqrt{u_j}) \le C(T, \Omega); \quad (11) $$

$$ E(u_1, u_2)(T) + 4 \sum_{i=1}^{2} d \| \nabla \sqrt[4]{u_i} \|_{L^2(Q_T)}^2 + 4 \| \nabla \sqrt[4]{u_1 u_2} \|_{L^2(Q_T)}^2 \le E(u_1, u_2)(0) + C(T, \Omega), \quad (12) $$

where

$$ E(u_1, u_2) := \int_{\Omega} \left( u_1 - 2\sqrt{u_1} + 1 \right) + \int_{\Omega} \left( u_2 - 2\sqrt{u_2} + 1 \right). $$

**Remark 2.2.** Even in low dimension, the use of Sobolev inequalities does not easily provide estimates better than what is obtained thanks to the duality Lemma.

Indeed, if we observe that thanks to estimate (12), $\nabla(\sqrt[4]{u_i}) \in L^2_{t,x}$, we see by the Sobolev injection $H_x^1 \hookrightarrow L_x^{M^*}$ (with $M^* = \frac{2M}{M-2}$ if $M > 2$, $2^* = \infty$ if $M = 1$ and any $p < \infty$ if $M = 2$), we get $\sqrt[4]{u_i} \in L^2_t(L^2_x^*)$. Interpolating with the estimate $\sqrt[4]{u_i} \in L^\infty_t(L^4_x)$, we end up with $u_i \in L_{t,x}^{\frac{M+4}{2M}}$ if $M > 2$, $u_i \in L_{t,x}^{3/2}$ if $M = 1$ and $u_i \in L_{t,x}^p$ for any $p < 3/2$ for $M = 2$.

We see that in all cases, we end up with an estimate which is not as good as what is provided by the duality Lemma (that is, estimate (11)). This estimate turns out to be crucial in our proof, since it allows to avoid any concentration phenomenon in the non-linearity $u_i\sqrt{u_j}$, a term which may not be controlled with the only use of estimates (9)–(10)–(12).
We now present a weak stability estimate which Proof contains the same ingredients as the Proof of our main Theorem, but in which do not appear the technical difficulties related to the approximation processes.

**Proposition 2.3.** Consider the assumptions of Proposition 2.1. We also assume that $\beta < 1$.

Let $(u_{1,n}, u_{2,n})_{n \in \mathbb{N}}$ be a sequence of strong (that is, belonging to $\mathcal{C}^2([0, T] \times \bar{\Omega})$) positive solutions to the reaction-cross diffusion equation (6) – (8). We also assume that $(u_{1,n}(0, \cdot), u_{2,n}(0, \cdot))_{n \in \mathbb{N}}$ converges weakly in $(L^2(\Omega))^2$ towards some functions $(u_1^{in}, u_2^{in})$ (in $(L^2(\Omega))^2$).

Then the sequence $(u_{i,n})_{n \in \mathbb{N}}$ ($i = 1, 2$) converges up to extraction a.e. and weakly in $L^2(Q_T)$ towards some function $u_i \in L^2(Q_T)$.

Moreover, $(u_1, u_2)$ satisfies estimates (9) – (11), and is a weak solution of eq. (6) – (8) in the following sense: For all $\phi \in \mathcal{C}_c^2([0, T[ \times \bar{\Omega})$ satisfying the Neumann boundary conditions, and $i \neq j = 1, 2$,

$$ - \int_{\Omega} u_i^{in} \phi(0, \cdot) - \int_{0}^{T} \int_{\Omega} \partial_t \phi u_i - \int_{0}^{T} \int_{\Omega} \Delta \phi [u_i(d + \sqrt{u_j})] \\ = \int_{0}^{T} \int_{\Omega} \phi u_i (1 - u_i^{\beta}). \quad (13) $$

Note that all terms in the weak formulation above are well-defined thanks to the assumptions on the coefficients of the equation and thanks to the a priori estimates.

*Proof of Proposition 2.1.* Integrating on $[0, t[ \times \Omega$ eq. (6) – (7), we obtain estimate (9) thanks to Gronwall's Lemma. Integrating the same system on $[0, T[ \times \Omega$ and using this time (9), we get (10).

In order to prove (11), we will use the following duality estimate, similar to the classical duality Lemma (see [27] for instance):

**Lemma 2.4.** We consider $\Omega$ an open bounded subset of $\mathbb{R}^M$ with a $\mathcal{C}^2$ boundary, and $T > 0$. We also consider $K > 0$ a constant and $\mu := \mu(t, x)$ a $\mathcal{C}^2([0, T] \times \bar{\Omega})$ positive function.

Let $u \ge 0$ be a strong (that is, belonging to $\mathcal{C}^2([0, T] \times \bar{\Omega})$) solution of the inequality

$$ \partial_t u - \Delta(\mu u) \le K, $$

on $Q_T$, with Neumann boundary condition.

Then

$$ \|\sqrt{\mu}u\|_{L^2(Q_T)} \le 2(KT + |\Omega|^{-1/2} \|u(0, \cdot)\|_{L^2(\Omega)}) \|\sqrt{\mu}\|_{L^2(Q_T)} + C_{\Omega}\|u(0, \cdot)\|_{L^2(\Omega)}, \quad (14) $$

where $C_{\Omega}$ is the Poincaré-Wirtinger constant.

*Proof of Lemma 2.4.* We introduce, for any smooth $(C^2([0, T] \times \bar{\Omega}))$ nonpositive function $\phi := \phi(t, x)$, the nonnegative solution $v := v(t, x)$ of the backward-in-time dual parabolic problem

$$ \partial_t v + \mu \Delta v = \sqrt{\mu} \phi, $$

with Neumann boundary condition and final condition $v(T, \cdot) = 0$.

Multiplying by $v$ the inequality satisfied by $u$, we end up with the identity

$$ - \int_{0}^{T} \int_{\Omega} u \sqrt{\mu} \phi \le K \int_{0}^{T} \int_{\Omega} v + \int_{\Omega} u(0, \cdot) v(0, \cdot). \quad (15) $$

Then, multiplying by $\Delta v$ the equation satisfied by $v$ and integrating on $[0, T]$ we see that

$$ \frac{1}{2} \int_{\Omega} |\nabla v(t, \cdot)|^2 + \int_{t}^{T} \int_{\Omega} \mu (\Delta v)^2 = \int_{t}^{T} \int_{\Omega} \sqrt{\mu} \phi \Delta v, $$

so that using Young's inequality, we end up with

$$ \frac{1}{2} \int_{\Omega} |\nabla v(0, \cdot)|^2 + \frac{1}{2} \int_{0}^{T} \int_{\Omega} \mu (\Delta v)^2 \le \frac{1}{2} \int_{0}^{T} \int_{\Omega} \phi^2. \quad (16) $$
Then integrating the equation on $v$ on $[t, T] \times \Omega$ and using (16) we get for all $t \in [0, T]$

$$ \int_{\Omega} v(t, \cdot) \le 2 \left( \int_{0}^{T} \int_{\Omega} \phi^2 \right)^{1/2} \left( \int_{0}^{T} \int_{\Omega} \mu \right)^{1/2}, $$

so that

$$ \begin{aligned} \left| \int_0^T \int_\Omega v \right| &\le 2T \left( \int_0^T \int_\Omega \phi^2 \right)^{1/2} \left( \int_0^T \int_\Omega \mu \right)^{1/2}, \\ \left( \int_\Omega v(0, \cdot)^2 \right)^{1/2} &\le C_\Omega \left( \int_0^T \int_\Omega \phi^2 \right)^{1/2} + 2 |\Omega|^{-1/2} \left( \int_0^T \int_\Omega \phi^2 \right)^{1/2} \left( \int_0^T \int_\Omega \mu \right)^{1/2}, \end{aligned} $$

where $C_\Omega$ is the Poincaré-Wirtinger constant. Going back to (15), we hence have

$$ \begin{aligned} \left| \int_0^T \int_\Omega u \sqrt{\mu} \phi \right| &\le 2 K T \|\sqrt{\mu}\|_{L^2(Q_T)} \|\phi\|_{L^2(Q_T)} + \|u(0, \cdot)\|_{L^2(\Omega)} \|v(0, \cdot)\|_{L^2(\Omega)} \\ &\le 2(KT + |\Omega|^{-1/2} \|u(0, \cdot)\|_{L^2(\Omega)} \|\sqrt{\mu}\|_{L^2(Q_T)} \|\phi\|_{L^2(Q_T)} + C_\Omega \|u(0, \cdot)\|_{L^2(\Omega)} \|\phi\|_{L^2(Q_T)} \end{aligned} $$

which is exactly the expression of (14) by duality so that Lemma 2.4 is obtained. □

Now notice that

$$ \partial_t u_i - \Delta(\mu_j u_i) \le 1, $$

where $\mu_j := d + \sqrt{u_j}$ ($i \ne j$) is a positive $\mathcal{C}^2([0, T] \times \bar{\Omega})$ function. The previous Lemma hence applies and we have

$$ \int_0^T \int_\Omega u_i^2 (d + \sqrt{u_j}) \le C(T, \Omega) + C(T, \Omega) \int_0^T \int_\Omega (d + \sqrt{u_j}), $$

which, using (9) and $\sqrt{u_j} \le 1 + u_j$ gives (11).

It remains to prove (12). Recall

$$ E(u_1, u_2) := \int_{\Omega} (u_1 - 2\sqrt{u_1} + 1) + \int_{\Omega} (u_2 - 2\sqrt{u_2} + 1). $$

Differentiating this energy along the flow of eq. (6) – (8), we get

$$ \begin{aligned} \frac{dE}{dt} + \sum_{i=1}^{2} \int_{\Omega} (\sqrt{u_i} + u_i^{1+\beta}) &= \sum_{i=1}^{2} \int_{\Omega} (u_i + u_i^{\beta+1/2}) - \sum_{i=1}^{2} \frac{1}{2} \int_{\Omega} d u_i^{-3/2} |\nabla u_i|^2 \\ &\quad - \sum_{i \ne j} \frac{1}{2} \int_{\Omega} u_j^{-3/2} \nabla u_j \cdot \left\{ \sqrt{u_i} \nabla u_j + \frac{1}{2} u_i^{-1/2} u_j \nabla u_i \right\}. \end{aligned} $$

Integrating w.r.t. $t \in [0, T]$, we end up with

$$ \begin{aligned} E(T) + \sum_{i=1}^{2} \frac{1}{2} \int_{0}^{T} \int_{\Omega} d u_{i}^{-3/2} |\nabla u_{i}|^{2} &+ \int_{0}^{T} \int_{\Omega} \left[ \frac{1}{|\nabla u_{1}|} \right]^{2} + \left[ \frac{1}{|\nabla u_{2}|} \right]^{2} + \frac{\nabla u_{1} \cdot \nabla u_{2}}{u_{1}} \\ &+ \sum_{i=1}^{2} \int_{0}^{T} \int_{\Omega} (\sqrt{u_{i}} + u_{i}^{1+\beta}) \\ &= E(0) + \sum_{i=1}^{2} \int_{0}^{T} \int_{\Omega} (u_{i} + u_{i}^{\beta+1/2}). \end{aligned} \quad (17) $$

It can be checked easily that the integrand in the last integral of the first line of identity (17) is larger than $|\nabla \sqrt{u_1 u_2}|^2$. Using $u_i^{\beta+1/2} \le 1 + u_i^{\beta+1}$ and estimates (9)–(10), we get estimate (12). This ends the Proof of Proposition 2.1. □
Proof of Proposition 2.3. We first observe that since we assumed the sequences $(u_{i,n}(0, \cdot))_{n \in \mathbb{N}}$ ($i = 1, 2$) to be weakly converging in $L^2(\Omega)$, they also are bounded in this space, so that according to Proposition 2.1, estimate (11) holds uniformly w.r.t. $n \in \mathbb{N}$, and the sequences $(u_{i,n})_{n \in \mathbb{N}}$ ($i = 1, 2$) are, up to some extraction, converging weakly in $L^2(Q_T)$ towards some function $u_i \in L^2(Q_T)$.

In order to pass to the limit (up to an extraction) in the nonlinear terms in eq. (13) it is sufficient to show on one hand that (for $i \neq j$) the sequences $(u_{i,n}^{1+\beta})_{n \in \mathbb{N}}$, $(u_{i,n} \sqrt{u_{j,n}})_{n \in \mathbb{N}}$ are bounded in $L^p(Q_T)$ for some $p > 1$, and on the other hand, that the sequence $(u_{i,n}(t, x))_{n \in \mathbb{N}}$ ($i = 1, 2$) converges for a.e. $(t, x) \in Q_T$ towards $u_i(t, x)$.

Since $\beta_i < 1$ we first deduce from the bound of $(u_{i,n})_{n \in \mathbb{N}}$ in $L^2(Q_T)$, a bound for $(u_{i,n}^{1+\beta})_{n \in \mathbb{N}}$ in some $L^p(Q_T)$ with $p > 1$.

Thanks to estimate (11), we see that $(u_{i,n} \sqrt[4]{u_{j,n}})_{n \in \mathbb{N}}$ is bounded in $L^2(Q_T)$ and $(\sqrt[4]{u_{j,n}})_{n \in \mathbb{N}}$ is clearly bounded in $L^4(Q_T)$ so that the product of these two sequences (that is the last nonlinearity to handle) is bounded in $L^{4/3}(Q_T)$, thanks to Hölder inequality.

We now turn to the a.e. convergence of $(u_{j,n})_{n \in \mathbb{N}}$ ($j = 1, 2$). Since $(u_{j,n})_{n \in \mathbb{N}}$ is bounded in $L^2(Q_T)$, we also have its boundedness in $L^{3/2}(Q_T)$. Estimate (12) holds uniformly w.r.t. $n$, so that $(\nabla \sqrt[4]{u_{j,n}})_{n \in \mathbb{N}}$ is bounded in $L^2(Q_T)$. Noticing that for any positive regular function $f$

$$ \nabla f = 4f^{3/4} \nabla \sqrt[4]{f}, $$

we get by the Cauchy-Schwarz inequality the boundedness of $(\nabla u_{j,n})_{n \in \mathbb{N}}$ in $L^1(Q_T)$. We see that $(u_{j,n})_{n \in \mathbb{N}}$ is bounded in $L^1([0, T]; W^{1,1}(\Omega))$. Moreover, thanks to eq. (6), (7) and the already mentioned estimates on the sequences $(u_{i,n}^{1+\beta})_{n \in \mathbb{N}}$, $(u_{i,n} \sqrt{u_{j,n}})_{n \in \mathbb{N}}$, we see that $(\partial_t u_i)_{n \in \mathbb{N}}$ is bounded in $L^1([0, T]; W^{-2,1}(\Omega))$. Using Aubin-Simon's Lemma (see [29] for instance) we obtain the convergence for a.e. $(t, x) \in Q_T$ of $(u_i(t, x))_{n \in \mathbb{N}}$ towards $u_i(t, x)$. □

## 2.2 Road map for the Proof of the main Theorem

The approximation processes used in order to prove the existence of solutions to cross diffusion equations are unfortunately often quite complex. In order to help the readers who wish to understand all the details of these processes, we present here a short road map:

* In preliminary Section 3 we give the general version of the entropy estimate (12), which becomes (20). Several technical results linked to this entropy structure are given, together with two standard results that will be useful in the sequel.

* In Section 4 the existence theory is focused on a new set of variables. We consider a time discretization of the system, together with a Galerkin type projection in the spirit of [8] and some smooth truncation of the reaction terms and the self-diffusion. This procedure preserves the entropy structure in its approximate form, this is the goal of Subsection 4.2.

* Section 5 handles the passage from finite to infinite dimension (for the space variable only). Since the time discretization is fixed, compactness in space is sufficient to pass to the limit, and it is mainly obtained thanks to (22) – (23). Of course, since the truncations of the reaction and the self-diffusion terms are fixed in this step, the (possible) concentration phenomena are easily handled, and the duality estimate (goal of the next Section) is not yet required.

* Section 6 is devoted to the Proof of the duality estimate which is the counterpart of the previous Lemma 2.4. The difficulties are twofold. First, we are looking at a discretized (in time) version of the parabolic inequality appearing in Lemma 2.4, we hence have to check carefully all the performed estimates. Secondly (and this is the trickiest part), the function governing the diffusion ($\mu$ in Lemma 2.4) is here a priori not a regular function (it is given by the approximation procedure itself). In order to recover the estimate given in Lemma 2.4, we hence have to consider a smooth approximation of the diffusion coefficients (that is, the cross and the self diffusion): this the goal of Lemma 6.7.

* The last section (that is, Section 7) is devoted to the relaxation of all remaining approximations (that is, the discrete time variable becomes continuous, and all truncations are removed), which leads to the Proof of our main Theorem. The duality estimate proved in Section 6 is crucial in this step (to avoid any concentration phenomena).
# 3 Preliminaries

In this Section, we introduce the entropy structure of our problem. The computations and estimates which are presented in the four first Subsections correspond to (part of) the establishment of eq. (12) in the simplified case studied in Section 2.

## 3.1 Entropy structure

We first introduce some notations which enable to rewrite our system under a form in which a Lyapounov functional naturally appears.

**Definition 3.1.** We define

$$a_{ii}(u_i) = u_i d_{ii}(u_i). \qquad (18)$$

For given cross diffusion parameters $a_{12}$ and $a_{21}$ (satisfying assumption **H2**), we introduce the functions, $\psi_i$ (for $i = 1, 2$), as the only elements of $\mathcal{C}^2([0, +\infty])$ verifying

$$\psi_1''(z) = \frac{a_{21}'(z)}{z}, \quad \psi_1(1) = \psi_1'(1) = 0,$$

$$\psi_2''(z) = \frac{a_{12}'(z)}{z}, \quad \psi_2(1) = \psi_2'(1) = 0.$$

We also define, for $(t, x) \in \Omega_T$,

$$w_i(t, x) := \psi_i'(u_i(t, x)).$$

One can then rewrite the system in a symmetric form:

$$\partial_t \begin{pmatrix} u_1 \\ u_2 \end{pmatrix} - \operatorname{div} \left( \underbrace{\begin{pmatrix} a_{11}'(u_1) + a_{12}(u_2) \\ a_{21}'(u_1) \end{pmatrix}}_{u_1 u_2} u_1 \quad \underbrace{\begin{pmatrix} u_1 u_2 \\ \frac{a_{22}'(u_2) + a_{21}(u_1)}{a_{12}'(u_2)} \end{pmatrix}}_{u_2} \right) \begin{pmatrix} \nabla w_1 \\ \nabla w_2 \end{pmatrix} = \begin{pmatrix} u_1(r_1 - s_{11}(u_1) - s_{12}(u_2)) \\ u_2(r_2 - s_{22}(u_2) - s_{21}(u_1)) \end{pmatrix}. \quad (19)$$

The terms $\nabla w_i$ have to be considered as scalars for the matrix product, the divergence being understood line by line (after the matrix product).

Multiplying the equations of the system by $w_1, w_2$ and integrating in space we obtain formally the entropy identity

$$\frac{d}{dt} \int_{\Omega} (\psi_1(u_1) + \psi_2(u_2)) + \int_{\Omega} (\nabla w_1, \nabla w_2) A(u_1, u_2) \begin{pmatrix} \nabla w_1 \\ \nabla w_2 \end{pmatrix} = (w_1, w_2) \begin{pmatrix} u_1(r_1 - s_{11}(u_1) - s_{12}(u_2)) \\ u_2(r_2 - s_{22}(u_2) - s_{21}(u_1)) \end{pmatrix}. \quad (20)$$

This identity motivates the following definition

**Definition 3.2.** We introduce the entropy (for some vector function $u : \Omega^2 \to \mathbb{R}_+^2$):

$$\mathcal{E}(u) := \sum_{i=1}^{2} \int_{\Omega} \psi_i(u_i)(x) dx.$$

Due to technical difficulties explained in Subsection 3.4, an approximated version of this entropy will be introduced later, see Definition 3.8.

## 3.2 Properties of the symmetric matrix

For $u_1, u_2 > 0$, $u := (u_1, u_2)$, the matrix $A(u)$ is defined by

$$A(u_1, u_2) := \begin{pmatrix} \frac{a_{11}'(u_1) + a_{12}(u_2)}{a_{21}'(u_1)} & u_1 u_2 \\ u_1 u_2 & \frac{a_{22}'(u_2) + a_{21}(u_1)}{a_{12}'(u_2)} \end{pmatrix}, \quad (21)$$
and we denote the associate quadratic form $Q(u)$. For a function $u = (u_1, u_2) : \Omega \to ]0, +\infty]^2$, we defined in the previous subsection $w = (w_1, w_2) = (\psi'_1(u_1), \psi'_2(u_2))$. We will have to deal with the expression

$$Q(u)(\nabla w) = [^t\nabla w][A(u)\nabla w],$$

since this term naturally appears in the entropy estimate. We establish the following Proposition, that will be useful later.

**Proposition 3.3.** Under the H assumptions, the application defined by (21) $A : ]0, +\infty]^2 \to \mathcal{M}_2(\mathbb{R})$ ($2 \times 2$ real matrices) belongs to $\mathcal{C}^0(]0, +\infty]^2)$ and takes its values in $S_2^{++}(\mathbb{R})$ (space of positive symmetric matrices). We have furthermore the two following estimates on the quadratic form $Q$

$$Q(u)(\nabla w) \geq \frac{1}{u_1} a'_{21}(u_1) |\nabla u_1|^2 + \frac{1}{u_2} a'_{12}(u_2) |\nabla u_2|^2, \qquad (22)$$

$$Q(u)(\nabla w) \geq 4 \left\{ \left| \nabla \sqrt{a_{21}(u_1)} \right|^2 + \left| \nabla \sqrt{a_{12}(u_2)} \right|^2 + \left| \nabla \sqrt{a_{21}(u_1) a_{12}(u_2)} \right|^2 \right\}. \qquad (23)$$

*Proof.* The fact that $A \in \mathcal{C}^0(]0, +\infty]^2)$ is an easy consequence of the regularity of the functions $a_{ij}$. As for the positiveness of the matrix, we just decompose $A$ between self (clearly positive) and cross diffusion:

$$A(u_1, u_2) = \overbrace{\left( \begin{array}{cc} \overbrace{\begin{pmatrix} a'_{11}(u_1) & 0 \\ 0 & a'_{22}(u_2) \end{pmatrix}}^{B(u_1, u_2)} & 0 \end{array} \right)}^{C(u_1, u_2)} + \overbrace{\left( \begin{array}{cc} u_1 u_2 & 0 \\ u_1 u_2 & a'_{12}(u_2) \end{array} \right)}^{C(u_1, u_2)},$$

We see that due to the assumptions on $a_{12}$ and $a_{21}$, the cross diffusion matrix $C(u_1, u_2)$ is nonnegative. Indeed, for $u_1, u_2 > 0$

$$\frac{a_{12}(u_2)}{a'_{21}(u_1)} u_1 \times \frac{a_{21}(u_1)}{a'_{12}(u_2)} u_2 \geq u_1 u_2 \times u_1 u_2,$$

is equivalent to

$$a_{12}(u_2)a_{21}(u_1) \geq u_1u_2a'_{12}(u_2)a'_{21}(u_1),$$

which is true by concavity. We get

$$Q(u)(\nabla w) = [^t\nabla w][A(u)\nabla w] \geq [^t\nabla w][B(u)\nabla w] = \frac{a'_{11}(u_1)}{a'_{21}(u_1)}u_1|\nabla w_1|^2 + \frac{a'_{22}(u_2)}{a'_{12}(u_2)}u_2|\nabla w_2|^2.$$

Since $w_i = \psi'_i(u_i)$, with $\psi''_i(x) = \frac{a'_{ji}(x)}{x}$, we end up with

$$Q(u)(\nabla w) \geq \frac{a'_{11}(u_1)}{u_1} a'_{21}(u_1) |\nabla u_1|^2 + \frac{a'_{22}(u_2)}{u_2} a'_{12}(u_2) |\nabla u_2|^2,$$

which, using assumption H3, leads to the first lower bound (22). On the other hand, expanding directly $Q(u)(\nabla w)$ with the definitions of $A(u)$ and $w$, we get

$$Q(u)(\nabla w) = [a'_{11}(u_1) + a_{12}(u_2)] \frac{a'_{21}(u_1)}{u_1} |\nabla u_1|^2 + [a'_{22}(u_2) + a_{21}(u_1)] \frac{a'_{12}(u_2)}{u_2} |\nabla u_2|^2 \\ + 2 a'_{21}(u_1) a'_{12}(u_2) \langle \nabla u_1, \nabla u_2 \rangle.$$

The concavity property $a_{ji}(u_i) \geq a'_{ji}(u_i)u_i$ used before can be written (here for $i = 1, j = 2$)

$$\frac{a'_{21}(u_1)}{u_1} \geq \frac{a'_{21}(u_1)^2}{a_{21}(u_1)},$$
and together we the assumption H3, we get

$$
\begin{align*}
Q(u)(\nabla w) &\ge a'_{11}(u_1) \frac{a'_{21}(u_1)^2}{a_{21}(u_1)} |\nabla u_1|^2 + a'_{22}(u_2) \frac{a'_{12}(u_2)^2}{a_{12}(u_2)} |\nabla u_2|^2 \\
&\quad + a_{12}(u_2) \frac{a'_{21}(u_1)^2}{a_{21}(u_1)} |\nabla u_1|^2 + a_{21}(u_1) \frac{a'_{12}(u_2)^2}{a_{12}(u_2)} |\nabla u_2|^2 + 2a'_{21}(u_1)a'_{12}(u_2) \langle \nabla u_1, \nabla u_2 \rangle \\
&= 4|\nabla \sqrt{a_{21}(u_1)}|^2 + 4|\nabla \sqrt{a_{12}(u_2)}|^2 \\
&\quad + 4|\sqrt{a_{12}(u_2)}\nabla \sqrt{a_{21}(u_1)}|^2 + 4|\sqrt{a_{21}(u_1)}\nabla \sqrt{a_{12}(u_2)}|^2 \\
&\quad + 4 \times 2\sqrt{a_{21}(u_1)}\sqrt{a_{12}(u_2)} \langle \nabla \sqrt{a_{21}(u_1)}, \nabla \sqrt{a_{12}(u_2)} \rangle,
\end{align*}
$$

and we thus end up with eq. (23).

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

3.3 Properties of the functions ψ$_{i}$

We shall need in the sequel the following elementary result:

**Lemma 3.4.** *Take h, ℓ ∈ ℂ⁰(R⁺) ∩ ℂ¹([0, +∞[) with h concave and ℓ nonnegative and convex, with
ℓ'(z) > 0 for all z large enough.*

Then there exists a constant $A_{h,\ell} > 0$ such that $h(x) \le A_{h,\ell}(1 + \ell(z))$ for all $z \in \mathbb{R}_+$.

Proof. If h is bounded from above, then $A_{h,l} = \sup h$ works. Otherwise, $h' > 0$ on $\mathbb{R}_+$. Say that $0 < \ell'(z)$ for $z > Z$. Then for all $z \in \mathbb{R}_+$, $z - Z \le \frac{\ell(z) - \ell(Z)}{\ell'(Z)}$. Then since $h' > 0$ and $h$ is concave, we can write

$$
h(z) \le h(Z) + h'(Z)(z - Z) \le h(Z) + \frac{h'(Z)}{\ell'(Z)}(\ell(z) - \ell(Z)) \le h(Z) + \frac{h'(Z)}{\ell'(Z)}\ell(z),
$$

so that $A_{h,l} = \max \left( h(Z), \frac{h'(Z)}{\ell'(Z)} \right)$ works. This concludes the Proof of Lemma 3.4. $\square$

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

**Lemma 3.5.** We assume **H2** on the coefficients $a_{ij}$ ($i \neq j$). Then the $\psi_i$ (cf. **Definition 3.1**) are convex functions with increasing derivatives. Moreover (for $i \neq j$):

(i) For all $z, y \in \mathbb{R}_+$, $\psi'_i(z)(z - y) \ge \psi_i(z) - \psi_i(y)$.

(ii) $\psi'_i(z) = o_{z \to 0^+}(1/z)$ and hence $z\psi'_i(z) \ge B$, for some constant $B < 0$ for all $z \in ]0, +\infty[$.

(iii) $\psi_i$ has a limit at point $0^+$ ($\psi_1(0) = a_{21}(1)$), furthermore $\psi_i$ is positive on $\mathbb{R}_+$.

(iv) There exists a constant $D > 0$ such that, for all $z \in \mathbb{R}_+$

$$
\begin{align*}
\forall \alpha \in [0, 1], \quad z^{\alpha} + a_{ji}(z) &\le D(1 + \psi_i(z)), \\
z\psi'_i(z) &\le D(1 + \psi_i(z)), \\
\psi_i(z) &\le D(1 + z \ln z - z).
\end{align*}
$$

Proof. Let us treat only the case $i=1$, the other one being similar.

(i) $\psi_1$ is convex.

(ii) We have

$$
\psi'_{1}(z) = \int_{1}^{z} \frac{a'_{21}(t)}{t} dt = \left[ \frac{a_{21}(t)}{t} \right]_{1}^{z} + \int_{1}^{z} \frac{a_{21}(t)}{t^2} dt = \int_{z \to 0^{+}} (1/z),
$$

since $a_{21} \in \mathscr{C}^0(\mathbb{R}_+)$ and $a_{21}(0) = 0$. The function $z \mapsto z\psi_1'(z)$ is increasing after $z=1$ and bounded near 0, hence lower bounded.

(iii) Note that $\psi_1''(t) = a_{21}'(t)$. We have

$$
a_{21}(1) - a_{21}(z) = \int_{z}^{1} t \psi_{1}''(t) dt = \left[ t \psi_{1}'(t) \right]_{z}^{1} - \int_{z}^{1} \psi_{1}'(t) dt,
$$
then recalling that $\psi_i(1) = \psi'_i(1) = a_{21}(0) = 0$, we have

$$a_{21}(1) - a_{21}(z) = \psi_{1}(z) - z\psi'_{1}(z).$$

hence the previous point (ii) gives the limit near 0. For the positiveness, just notice that $\psi'_1$ is nonpositive on $[0, 1]$.

(iv) We use Lemma 3.4 with $h(z) = z + a_{ji}(z)$, $\ell(z) = \psi_1(z)$ and $z^\alpha \le 1 + z$ in order to obtain the first inequality. For the second inequality, we use the same Lemma with $h(z) = z\psi'_1(z) - 2\psi_1(z)$, and $\ell(z) = \psi_1(z)$. Since $a'_{21}$ is bounded, integrating twice the formula in Definition 3.1, we get the third inequality. □

## 3.4 A small perturbation

Since $\psi'_i$ may not be one to one, we have to use in the approximation problem a small perturbation of this function. We introduce the following definition:

**Definition 3.6.** Let us assume **H2** (and recall Definition 3.1) on the coefficients $a_{ij}$ ($i \neq j$), and **H3** on the coefficients $d_{ii}$.

For all $\varepsilon > 0$ (small enough), we introduce

$$\psi_i^\varepsilon(z) := \psi_i(z) + \varepsilon z \ln(z) - \varepsilon z,$$

and, equivalently, (for $i \neq j$) $a_{ij}^\varepsilon := a_{ij} + \varepsilon z$. We have hence $\psi_i^{\varepsilon'}(z) := \psi_i'(z) + \varepsilon \ln(z)$.

We also introduce $a_{ii}^\varepsilon(z) = z d_{ii}^\varepsilon(z)$, with $d_{ii}^\varepsilon = \gamma_\varepsilon(d_{ii})$, where $(\gamma_\varepsilon)_{\varepsilon>0}$ is a smooth approximation of $x \mapsto \min(x, \varepsilon^{-1})$ on $\mathbb{R}_+$, uniformly converging to the identity on compact sets.

Finally, we denote by $A^\varepsilon$ the matrix $A$ (defined by (21)), where the coefficients $a_{ij}$ are replaced by $a_{ij}^\varepsilon$ (for all $i, j \in \{1, 2\}$), and $Q^\varepsilon$ the corresponding quadratic form.

The crucial point is the following: Proposition 3.3 is still true when one replaces the functions $a_{ij}$ and $A$ by their $\varepsilon$-approximations, $a_{ij}^\varepsilon$ and $A^\varepsilon$, (assumption **H2** and **H3** hold for the $a_{ij}^\varepsilon$), and if one tries to reproduce Lemma 3.5 with these new functions, all the inequalities remain the same, the constants being a little bit changed but uniformly bounded in $\varepsilon$. We write the following Lemma, which summarizes the situation. We skip its Proof since it is very close to the one of Lemma 3.5.

**Lemma 3.7.** We assume that **H2** and **H3** holds on the coefficients $a_{ij}$ (for all $i, j \in \{1, 2\}$), and use the notations of Definition 3.6.

Then, the coefficients $a_{ij}^\varepsilon$ also satisfy **H2** and **H3** (with constants independent of $\varepsilon$ for $\varepsilon > 0$ small enough). Also $B$ and $D$ (the constants of Lemma 3.5) may be changed in order to have, for $i = 1, 2$ and $0 < \varepsilon < 1$:

(i) For all $z, y \in \mathbb{R}_+$, $\psi_i^{\varepsilon'}(z)(z - y) \ge \psi_i^{\varepsilon}(z) - \psi_i^{\varepsilon}(y)$.

(ii) $\psi_i^{\varepsilon'}(z) = o_{z \to 0^+}(1/z)$ and furthermore, for all $z \in ]0, +\infty[$, $z\psi_i^{\varepsilon'}(z) \ge B - \varepsilon e^{-1}$.

(iii) $\psi_i^\varepsilon$ has a limit at point $0^+$. For $\varepsilon$ small enough, $\psi_i^\varepsilon$ is positive on $\mathbb{R}_+$.

(iv) For all $z \in \mathbb{R}_+$

$$\begin{aligned} \forall \alpha \in [0, 1], \quad & z^\alpha + a_{ji}^\varepsilon(z) \le D(1 + \varepsilon)(1 + \psi_i^\varepsilon(z)), \\ & z\psi_i^{\varepsilon'}(z) \le D(1 + \varepsilon)(1 + \psi_i^\varepsilon(z)), \\ & \psi_i^\varepsilon(z) \le D(1 + \varepsilon)(1 + z \ln z - z), \end{aligned}$$

where $D$ is the constant defined above in Lemma 3.5.

We end up this Subsection by introducing (bearing in mind the notations of Definition 3.6) the:

**Definition 3.8.** We define the (approximate) entropy (for some vector function $u : \Omega^2 \to \mathbb{R}_+^2$):

$$\mathcal{E}_\varepsilon(u) := \sum_{i=1}^{2} \int_{\Omega} \psi_i^\varepsilon(u_i)(x) \, dx.$$

**Remark 3.9.** Let us mention the following useful fact : because of points (iv) of Lemmas 3.5 and 3.7 we see that for any vector function $u : \Omega^2 \to \mathbb{R}_+^2$, $\mathcal{E}^\varepsilon(u)$ and $\mathcal{E}(u)$ are both finite as soon as $u$ is square integrable: the assumptions on the initial data in our main Theorem are therefore sufficient to ensure the finiteness of the entropy.
## 3.5 Two standard results

We complete this Section devoted to preliminaries by the statement of two results taken from the existing literature (the Proof is not presented here), and that will be used in the sequel.

**Lemma 3.10 (Discrete Gronwall).** Consider two nonnegative sequences $(v_n, w_n)_{n \in \mathbb{N}}$, satisfying for some constant $\theta \in ]0, 1[$,

$$ \forall n \in \mathbb{N}^*, \quad v_n \le v_{n-1} + \theta v_n + w_n. $$

Then, for all $n \in \mathbb{N}^*$:

$$ v_n \le e^{n\lambda_\theta} v_0 + \sum_{k=0}^{n-1} e^{k\lambda_\theta} w_{n-k} \le e^{\lambda_\theta} \left[ v_0 + \sum_{k=1}^{n} e^{(-k+1)\lambda_\theta} w_k \right], $$

with $\lambda_\theta := \theta/(1-\theta)$.

Moreover, if $w_n = C$ is constant, we have

$$ v_n \le e^{n\lambda_\theta} \left[ v_0 + \frac{C}{\theta} \right]. $$

The following Theorem can be found in [17] (in the more general case of an infinite dimensional Banach space, and provided that an extra compactness assumption is assumed) where it is presented as the “Leray-Schauder Theorem” (p.286, Theorem 11.6):

**Theorem 3.11.** Let $(E, \|\cdot\|)$ be a Euclidian vector space and $T: [0, 1] \times E \to E$ a continuous function satisfying $T(0, \cdot) \equiv 0$. Suppose furthermore the existence of $R > 0$ such that for any $s \in [0, 1]$, the following a priori estimate holds for the fixed points of $T(s, \cdot)$:

$$ T(s, x) = x \implies \|x\| < R. $$

Then $T(1, \cdot): E \to E$ has at least one fixed point in $B(0, R)$.

# 4 Approximate system of finite dimension

This Section is devoted to the establishment of the existence of a solution to an approximate problem (corresponding to (1) – (4)). This approximate problem is discretized in time and space (we use an implicit scheme for the time discretization, and a Galerkin approach for the space discretization). Moreover, some parameters (self and cross diffusion rates, reactive terms) are smoothed and truncated.

## 4.1 Notations

We start with a definition related to the discretization w.r.t. time, as in [8].

**Definition 4.1.** We decompose the time interval, $]0, T] = \bigcup_{k=1}^{N} ](k-1)\tau, k\tau]$, where $N \in \mathbb{N}^*$ and $\tau := T/N$, and introduce the finite difference operator : $\partial_\tau u^k := \frac{u^k - u^{k-1}}{\tau}$.

We also introduce new definitions related to the reaction terms:

**Definition 4.2.** We define

$$ R^{+}(u) = \begin{pmatrix} R_{12}^{+}(u) \\ R_{21}^{+}(u) \end{pmatrix} := \begin{pmatrix} r_1 & 0 \\ 0 & r_2 \end{pmatrix} \begin{pmatrix} u_1 \\ u_2 \end{pmatrix}, \quad R^{-}(u) = \begin{pmatrix} R_{12}^{-}(u) \\ R_{21}^{-}(u) \end{pmatrix} := \begin{pmatrix} u_1 & 0 \\ 0 & u_2 \end{pmatrix} \begin{pmatrix} s_{11}(u_1) + s_{12}(u_2) \\ s_{22}(u_2) + s_{21}(u_1) \end{pmatrix}, $$

and

$$ R(u) = \begin{pmatrix} R_{12}(u) \\ R_{21}(u) \end{pmatrix} = R^{+}(u) - R^{-}(u). $$
We introduce the following approximation for this reaction term

$$R^{\epsilon}(u) = R^{+}(u) - R^{-, \epsilon}(u),$$

with

$$R^{-, \epsilon}(u) := \begin{pmatrix} u_1 & 0 \\ 0 & u_2 \end{pmatrix} \begin{pmatrix} \gamma_{\epsilon}(s_{11}(u_1) + s_{12}(u_2)) \\ \gamma_{\epsilon}(s_{22}(u_2) + s_{21}(u_1)) \end{pmatrix},$$

where $\gamma_{\epsilon}$ is the truncation function that we used in Definition 3.6.

We then define our Galerkin approximation:

**Definition 4.3.** Let us consider a sequence $(V_n)_{n \in \mathbb{N}}$ of subspaces of $\mathcal{C}_{\nu}^{\infty}(\bar{\Omega})$, such that for all $n$, $V_n$ is $n$-dimensional, $V_n \subset V_{n+1}$, and $\cup_{n \in \mathbb{N}} V_n$ is dense in $L^2(\Omega)$. We assume that $V_1 = \langle \mathbf{1}_{\Omega} \rangle_{\mathbb{R}}$ (so that the constant function $\mathbf{1}_{\Omega}$ lies in all subspaces $V_n$).

For the sake of clarity, recall how we denote the two component vectors:

$$\chi := \begin{pmatrix} \chi_1 \\ \chi_2 \end{pmatrix}, u := \begin{pmatrix} u_1 \\ u_2 \end{pmatrix}, w := \begin{pmatrix} w_1 \\ w_2 \end{pmatrix} \quad \text{and} \quad R(u) := \begin{pmatrix} R_{12}(u_1, u_2) \\ R_{21}(u_2, u_1) \end{pmatrix}.$$

As previously, all linear operators (such as $\nabla, \partial_{\tau} \dots$) have to be understood line by line in the previous expressions. Any dot product $\langle \cdot, \cdot \rangle_{\mathbb{E}}$ on some space $\mathbb{E}$ of functions defined on $\Omega$ has to be understood as $\langle u, \chi \rangle_{\mathbb{E}} = \langle u_1, \chi_1 \rangle_{\mathbb{E}} + \langle u_2, \chi_2 \rangle_{\mathbb{E}}$.

We are going to define recursively solutions to a discrete-time problem,

all lying in $V_n$ ($n$ is fixed for the moment). Recall that the initial condition $u^0$ belongs to $L^2(\Omega)^2$. Given $u^{k-1} \in L^2(\Omega)^2$, we introduce the following problem:

**Definition 4.4.** For fixed $\epsilon > 0$, $\sigma \in [0, 1]$, $k \ge 1$, problem $P_{\sigma}^{\epsilon}(k, u^{k-1})$ consists in finding $w^k \in V_n^2$ such that $[\text{denoting } u^k = (\psi_{\epsilon}')^{-1}(w^k)]$ we have for all $\chi \in V_n^2$,

$$\sigma \left[ \langle \chi, \partial_{\tau} u^k \rangle_{L^2(\Omega)} + \langle \nabla \chi, A^{\epsilon}(u^k) \nabla w^k \rangle_{L^2(\Omega)} - \langle \chi, R^{\epsilon}(u^k) \rangle_{L^2(\Omega)} \right] = -\epsilon \langle \chi, w^k \rangle_{H^1(\Omega)}. \quad (24)$$

Since $w^k \in V_n \subset L^{\infty}(\Omega)$, $u^k$ takes its values in some compact set of $]0, +\infty]^2$, hence the coefficients of $A^{\epsilon}(u^k)$ all belong to $L^{\infty}(\Omega)$, and this ensures that all the previous brackets are well-defined. Furthermore, $u^k$ clearly belongs to $L^{\infty}(\Omega)^2 \subset L^2(\Omega)^2$, so that the previous inductive definition is consistent.

The collection of problems $(P_{\sigma}^{\epsilon}(k, u^{k-1}))_{k=1...N}$ is hence a discrete-in-time Galerkin-in-space version of system (1) – (4) in which an $\epsilon$-perturbation and a parameter $\sigma$ have been added.

## 4.2 A priori estimate

The first step in the Proof of existence for the problem defined in Definition 4.4 consists in writing the entropy estimate for this problem. This estimate is analogous to estimate (12) of the simplified problem of Section 2.

**Proposition 4.5.** Under the assumptions of Theorem 1.2 and keeping in mind the Definitions of Sub-section 4.1 [that is, Definitions 4.1 – 4.4], and Definitions 3.6, 3.8, there exists a constant $K_T$ depending only on $T$ and the data of the equation $(r_i, d_{ii}, a_{ij}, s_{ij})$ such that for any $\tau > 0$, $\epsilon > 0$ (small enough), and $\sigma \in [0, 1]$ and any sequence (of length less than $N$) $(w^j)_{1 \le j \le k}$, with $w^j$ solving $P_{\sigma}^{\epsilon}(j, u^{j-1})$ (see Definition 4.4), the following entropy inequality holds:

$$\sigma \mathcal{E}_{\epsilon}(u^k) + \tau \sigma \sum_{j=1}^{k} \int_{\Omega} Q^{\epsilon}(u^j)(\nabla w^j) dx + \epsilon \tau \sum_{j=1}^{k} \|w^j\|_{H^1(\Omega)}^2 \le K_T (\mathcal{E}_{\epsilon}(u^0) + 1). \quad (25)$$
Proof. Plug $w^j$ in $P_{\sigma}^{\varepsilon}(j, u^{j-1})$ (see Definition 4.4) to obtain

$$ \sigma \int_{\Omega} w^j \cdot \partial_{\tau} u^j \, dx + \sigma \int_{\Omega} Q^{\varepsilon}(u^j)(\nabla w^j) \, dx + \varepsilon \|w^j\|_{H^1(\Omega)}^2 = \sigma \int_{\Omega} w^j \cdot R^{\varepsilon}(u^j) \, dx. $$

Using (i) of Lemma 3.7, one gets, since $w_i = \psi_i^{\varepsilon'}(u_i)$,

$$ \int_{\Omega} w^j \cdot \partial_{\tau} u^j \, dx = \frac{1}{\tau} \int_{\Omega} w^j \cdot [u^j - u^{j-1}] \, dx \geq \frac{1}{\tau} (\mathcal{E}(u^j) - \mathcal{E}(u^{j-1})). $$

For the reaction term, we have

$$ \int_{\Omega} w^j \cdot R^{\varepsilon}(u^j) \, dx = \int_{\Omega} w_1^j R_{12}^{\varepsilon}(u_1^j, u_2^j) \, dx + \int_{\Omega} w_2^j R_{21}^{\varepsilon}(u_2^j, u_1^j) \, dx. $$

For the sake of clarity, we avoid writing the superscript j for a few lines. Let us focus on the first term of the right-hand side (the second one will be similar):

$$ \begin{aligned} \int_{\Omega} w_1 R_{12}^{\varepsilon}(u_1, u_2) \, dx &= \int_{\Omega} \psi_1^{\varepsilon'}(u_1) u_1 (r_1 - \gamma_{\varepsilon}[s_{11}(u_1) + s_{12}(u_2)]) \, dx \\ &= r_1 \int_{\Omega} \psi_1^{\varepsilon'}(u_1) u_1 \, dx - \int_{\Omega} \psi_1^{\varepsilon'}(u_1) u_1 \gamma_{\varepsilon}[s_{11}(u_1) + s_{12}(u_2)] \, dx. \end{aligned} $$

From (iv) of Lemma 3.7, we get the existence of a constant D such that

$$ r_1 \int_{\Omega} \psi_1^{\varepsilon'}(u_1) u_1 \, dx \leq r_1 D(1 + \varepsilon) (\mu(\Omega) + \mathcal{E}_{\varepsilon}(u)). $$

For the two other terms, we use the fact that $\psi_i^{\varepsilon'} \geq 0$ on $[1, +\infty[$, so that

$$ - \int_{\Omega} \psi_1^{\varepsilon'}(u_1) u_1 \gamma_{\varepsilon}[s_{11}(u_1) + s_{12}(u_2)] \, dx \leq (-B + \frac{1}{e}) \int_{\Omega} s_{11}(u_1) \, dx + (-B + \frac{1}{e}) \int_{\Omega} s_{12}(u_2) \, dx. $$

Finally, we use $s_{ij}(u) \leq K(1+u)$ to get

$$ - \int_{\Omega} \psi_1^{\varepsilon'}(u_1) u_1 \gamma_{\varepsilon}[s_{11}(u_1) + s_{12}(u_2)] \, dx \leq C \left( \mu(\Omega) + \int_{\Omega} [u_1 + u_2] \right). $$

By Lemma 3.7 (iv), we have

$$ \int_{\Omega} [u_1 + u_2] \leq D(1 + \varepsilon)(\mu(\Omega) + \mathcal{E}_{\varepsilon}(u)). $$

Finally we have, recovering the superscripts that we omitted before:

$$ \frac{\sigma}{\tau}(\mathcal{E}_{\varepsilon}(u^j) - \mathcal{E}_{\varepsilon}(u^{j-1})) + \sigma \int_{\Omega} Q^{\varepsilon}(u^j)(\nabla w^j) \, dx + \varepsilon \|w^j\|_{H^1(\Omega)}^2 \leq \sigma K(1 + \mathcal{E}_{\varepsilon}(u^j)), \quad (26) $$

for some constant K independent of j, n, ε, τ, σ. Hence if one takes τ small enough (such that $1 - \tau K \geq 1/2$), one has in particular:

$$ \mathcal{E}_{\varepsilon}(u^j) \leq \mathcal{E}_{\varepsilon}(u^{j-1}) + \tau K \mathcal{E}_{\varepsilon}(u^j) + \tau K, $$

which implies by Lemma 3.10, with $\theta = C := \tau K$,

$$ \mathcal{E}_{\varepsilon}(u^j) \leq e^{2j\theta} [\mathcal{E}_{\varepsilon}(u^0) + 1] \leq e^{2j\tau K} (\mathcal{E}_{\varepsilon}(u^0) + 1) \leq e^{2TK} (\mathcal{E}_{\varepsilon}(u^0) + 1). $$

Plugging this last inequality in (26), we get for some constant $C_T > 0$,

$$ \sigma(\mathcal{E}_{\varepsilon}(u^j) - \mathcal{E}_{\varepsilon}(u^{j-1})) + \tau \sigma \int_{\Omega} Q^{\varepsilon}(u^j)(\nabla w^j) \, dx + \tau \varepsilon \|w^j\|_{H^1(\Omega)}^2 \leq C_T \tau (1 + \mathcal{E}_{\varepsilon}(u^0)), $$

which, after summation over $j \in \{1, \dots, k\}$ and using $k\tau \leq T$, ends the Proof of Proposition 4.5. $\square$
## 4.3 Existence and estimates

The a priori estimate proven in the previous Subsection leads to the following Proposition about existence for the problem of Definition 4.4. We also indicate what becomes of those a priori estimates: in next Proposition, estimate (27) is analogous to (12) in the simplified problem of Section 2, and estimates (28) and (29) are analogous to (9) and (10) in the simplified problem of Section 2.

**Proposition 4.6.** We consider the assumptions of Theorem 1.2 and keep in mind the Definitions of Subsection 4.1 (that is, Definitions 4.1 – 4.4), and Definitions 3.6, 3.8.

For fixed $\tau(=T/N)$ (small enough), $n \in \mathbb{N}^*$, $w^0 \in V_n^2$ and $\varepsilon > 0$, there exists a sequence $(w^k)_{1\le k\le N}$ such that $w^k$ solves $P_1^\varepsilon(k, u^{k-1})$ (see Definition 4.4) for all $k \in \{1, \dots, N\}$, and which furthermore satisfies the estimate (with $K_T$ depending only on $T$ and the data of the equation $(r_i, d_{ii}, a_{ij}, s_{ij})$):

$$ \mathcal{E}_{\varepsilon}(u^k) + \tau \sum_{j=1}^{k} \int_{\Omega} Q^{\varepsilon}(u^j)(\nabla w^j) dx + \varepsilon \tau \sum_{j=1}^{k} \|w^j\|_{H^1(\Omega)}^2 \le K_T(\mathcal{E}_{\varepsilon}(u^0) + 1), \quad (27) $$

$$ \|u^k\|_{L^1(\Omega)} \le \left[ \|u^0\|_{L^1(\Omega)} + O(\sqrt{\varepsilon}) \right] e^{[r^\tau]T}, \quad (28) $$

$$ \tau \sum_{j=1}^{k} \|R^{-,\varepsilon}(u^j)\|_{L^1(\Omega)} \le \left[ \|u^0\|_{L^1(\Omega)} + O(\sqrt{\varepsilon}) \right] (1 + T r e^{[r^\tau]T}), \quad (29) $$

where the sequence of positive scalars $(r^\tau)_\tau$ goes to $r = \max(r_1, r_2)$ when $\tau \to 0$.

**Proof. Step 1: Estimates**

We first notice that since $1_\Omega \in V_n$, it is an admissible test function for the problem $P_1^\varepsilon(j, u^{j-1})$ (see Definition 4.4) for all $j \in \{1, \dots, k\}$, $k \in \{1, \dots, N\}$, and we hence get

$$ \|u^j\|_{L^1(\Omega)} - \|u^{j-1}\|_{L^1(\Omega)} + \varepsilon \tau \int_{\Omega} w^j + \tau \|R^{-,\varepsilon}(u^j)\|_{L^1(\Omega)} = \tau \|R^+(u^j)\|_{L^1(\Omega)}. \quad (30) $$

It follows easily from the definition of $R^+$ that the right-hand side of the previous inequality is not larger than $\tau r \|u^j\|_{L^1(\Omega)}$. Thus,

$$ \|u^j\|_{L^1(\Omega)} \le \|u^{j-1}\|_{L^1(\Omega)} + \tau r \|u^j\|_{L^1(\Omega)} - \varepsilon \tau \int_{\Omega} [w_1^j + w_2^j]. $$

We notice then that,

$$ -\varepsilon\tau \int_{\Omega} (w_1^j + w_2^j) = \sqrt{\varepsilon\tau} \left( |\Omega| + \varepsilon \|w^j\|_{L^2(\Omega)}^2 \right). $$

This leads to

$$ \|u^j\|_{L^1(\Omega)} \le \|u^{j-1}\|_{L^1(\Omega)} + \tau r \|u^j\|_{L^1(\Omega)} + \sqrt{\varepsilon\tau} \left( |\Omega| + \varepsilon \|w^j\|_{L^2(\Omega)}^2 \right), $$

which implies by Lemma 3.10

$$ \|u^j\|_{L^1(\Omega)} \le \left[ \|u^0\|_{L^1(\Omega)} + \sqrt{\varepsilon} \left( j\tau|\Omega| + \varepsilon\tau \sum_{k=1}^{j} \|w^j\|_{L^2(\Omega)}^2 \right) \right] \exp \left\{ j\tau \frac{r}{1 - \tau r} \right\}. $$

Using the entropy estimate, one has

$$ \varepsilon\tau \sum_{k=1}^{j} \|w^j\|_{L^2(\Omega)}^2 \le K_T(\mathcal{E}_{\varepsilon}(u^0) + 1) \quad \text{and} \quad j\tau|\Omega| \le |\Omega|T. $$

Therefore, we obtain

$$ \|u^j\|_{L^1(\Omega)} \le \left[ \|u^0\|_{L^1(\Omega)} + \sqrt{\varepsilon} (T|\Omega| + K_T(\mathcal{E}_{\varepsilon}(u^0) + 1)) \right] e^{r^\tau T}, $$

with $r(\tau)$ having the mentioned property, that is precisely (28).
Quite similarly, we sum up in (30) and obtain after using the bound on $\|w\|_1$:

$$ \|u^k\|_{L^1(\Omega)} + \tau \sum_{j=1}^{k} \|R^{-,\varepsilon}(u^j)\|_{L^1(\Omega)} = \tau \sum_{j=1}^{k} \|R^{+}(u^j)\|_{L^1(\Omega)} + \sqrt{\varepsilon} (T|\Omega| + K_T(\mathcal{E}_{\varepsilon}(u^0) + 1)), $$

and we conclude the estimate using the previous bound.

### Step 2: Existence

In the rest of the Proof, we will work in the finite dimensional Hilbert space $E = V_n^2$ with the dot product $\langle \cdot, \cdot \rangle_{H^1(\Omega)}$ and the associated norm. We will proceed by induction and only prove the first iteration $u^0 \in L^2(\Omega)^2 \Rightarrow \exists w^1$ solving $P_\sigma^\varepsilon(1, u^0)$ (see Definition 4.4), all the other induction steps will be similar.

First notice that the problem $P_\sigma^\varepsilon(1, u^0)$ can be seen in the following way

$$ \text{Find } w^1 \in E \text{ such that } \forall \chi \in E, \sigma L_{u^0, w^1}(\chi) = -\varepsilon \langle \chi, w^1 \rangle_{H^1(\Omega)}, $$

where, for $w, \chi \in E$, $L_{u^0, w}$ is defined by

$$ L_{u^0, w}(\chi) := \frac{1}{\tau} \langle \chi, (\psi'_{\varepsilon})^{-1}(w) - u^0 \rangle_{L^2(\Omega)} + \langle \nabla \chi, A^{\varepsilon}((\psi'_{\varepsilon})^{-1}(w)) \nabla w \rangle_{L^2(\Omega)} - \langle \chi, R^{\varepsilon}((\psi'_{\varepsilon})^{-1}(w)) \rangle_{L^2(\Omega)}. $$

As noticed before, for $w \in E \subset L^\infty(\Omega)^2$, $(\psi'_{\varepsilon})^{-1}(w)$ takes its values in some compact set of $]0, +\infty]^2$, so that the coefficients of $A^{\varepsilon}((\psi'_{\varepsilon})^{-1}(w))$ lie in $L^\infty(\Omega)$, whence $L_{u^0, w^1} \in V_n^*.$

Let us now define a map

$$ T: [0, 1] \times E \to E \\ (\sigma, v) \mapsto w, $$

such that:

$$ \forall \chi \in E \quad \sigma L_{u^0, v}(\chi) = -\varepsilon \langle \chi, w \rangle_{H^1(\Omega)}. $$

Such a map is well-defined because of the usual representation theorem for finite dimensional spaces. The proof of the following lemma is rather standard, so we do not detail it here. It may be found in [13].

### Lemma 4.7. The map T is continuous.

It is now time to use Proposition 4.5 to get an *a priori* estimate on any fixed point of $T(\sigma, \cdot)$, for any $\sigma \in [0, 1]$. In fact, the case $k=1$ of this Proposition exactly tells us that all this fixed points are in the ball of center 0 and radius $K_T(\mathcal{E}(u^0) + 1)/\varepsilon\tau$. Since clearly $T(0, \cdot) \equiv 0$, we now can apply Theorem 3.11 to see that $T(1, \cdot)$ has a fixed point, which is exactly the existence of $w^1$ and hence the first step of our induction machinery. Inequality (27) is then a direct consequence of (25) (with $\sigma = 1$). $\square$

The previous Proposition shows that (for fixed $\tau(= T/N)$ small enough, $n \in \mathbb{N}^*$, $u^0 \in L^2(\Omega)^2$, and $\varepsilon > 0$) there exists a sequence $(w^k)_{1 \le k \le N}$ such that, for all $k \in \{1, \dots, N\}$, denoting $u^k := (\psi'_{\varepsilon})^{-1}(w^k)$, we have for all $\chi \in V_n^2$:

$$ \langle \chi, \partial_{\tau} u^k \rangle_{L^2(\Omega)} + \langle \nabla \chi, A^{\varepsilon}(u^k) \nabla w^k \rangle_{L^2(\Omega)} + \varepsilon \langle \chi, w^k \rangle_{H^1(\Omega)} = \langle \chi, R^{\varepsilon}(u^k) \rangle_{L^2(\Omega)}, \quad (31) $$

that we also can write $(i \ne j = 1, 2)$

$$ \partial_{\tau} \mathbb{P}_n u_i^k - \mathbb{P}_n \Delta [a_{ii}^{\varepsilon}(u_i^k) + u_i^k a_{ij}(u_j^k) + \varepsilon u_i^k u_j^k] + \varepsilon w_i^k - \varepsilon \Delta w_i^k = \mathbb{P}_n [R_{ij}^{\varepsilon}(u_i^k, u_j^k)], $$

where $\mathbb{P}_n$ is the $L^2$-orthonormal projection on $V_n$.
# 5 Asymptotic with respect to $n$

This Section is devoted to the passage to the limit in the space discretization (that is, with the notations of the previous Section, $n \to +\infty$). It can be summarized by the following Proposition, in which existence and a priori estimates are established for an approximate problem (corresponding to (1) – (4)) in which the space is now continuous (but the time is still discretized, and some parameters are smoothed or truncated). In next Proposition, estimates (33) and (34) correspond to (12) in the simplified problem of Section 2, and estimates (35) and (36) are analogous to (9) and (10) in the simplified problem of Section 2.

**Proposition 5.1.** We consider the assumptions of Theorem 1.2 and keep in mind Definitions 4.1, 4.2, and Definitions 3.6, 3.8. For fixed $\tau(= T/N)$ (small enough), and $\varepsilon > 0$, there exists a sequence $(u^k)_{1 \le k \le N}$ of $L^p(\Omega)^2$ for some $p > 1$ (and a corresponding sequence $(w^k)_{1 \le k \le N}$) such that for all $\chi \in \mathcal{C}_\nu^\infty(\Omega)^2$ (and for all $k \in \{1, \dots, N\}$),

$$ 
\begin{aligned}
& \frac{1}{\tau} \langle \chi_i, u_i^k \rangle_{L^2(\Omega)} - \frac{1}{\tau} \langle \chi_i, u_i^{k-1} \rangle_{L^2(\Omega)} - \langle \Delta \chi_i, [a_{ii}^\varepsilon(u_i^k) + u_i^k a_{ij}(u_j^k) + \varepsilon u_i^k u_j^k] \rangle_{L^2(\Omega)} \\
& \qquad + \varepsilon \langle \chi_i - \Delta \chi_i, w_i^k \rangle_{L^2(\Omega)} = \langle \chi_i, R_{ij}^\varepsilon(u_i^k, u_j^k) \rangle_{L^2(\Omega)},
\end{aligned} 
\quad (32) $$

and which furthermore satisfies the estimates (with $K_T$ depending only on $T$ and the data of the equation $(r_i, d_{ii}, a_{ij}, s_{ij})$):

$$ \mathcal{E}_{\varepsilon}(u^k) + \varepsilon\tau \sum_{j=1}^{k} \|w^j\|_{H^1(\Omega)}^2 \leq K_T(\mathcal{E}_{\varepsilon}(u^0) + 1); \quad (33) $$

$$ \tau \sum_{j=1}^{k} \int_{\Omega} \left| \nabla \left( [u^j]^{\frac{1-\alpha}{2}} \right) \right|^2 dx \leq K_T (\mathcal{E}_{\varepsilon}(u^0) + 1), \quad (34) $$

$$ \|u^k\|_{L^1(\Omega)} \le \|u^0\|_{L^1(\Omega)} e^{T r_\tau}; \quad (35) $$

$$ \varepsilon\tau \sum_{j=1}^{k} \|w^j\|_{L^1(\Omega)} + \tau \sum_{j=1}^{k} \|R^{-,\varepsilon}(u^j)\|_{L^1(\Omega)} \le \|u^0\|_{L^1(\Omega)} [1 + \text{Tr} e^{r(\tau)T}]. \quad (36) $$

*Proof.* We first recall the bounds that hold on $u_n^k$ (and $w_n^k$), whose existence is given by Proposition 4.6. Using (22), for all $k \in \{1, \dots, N\}$ and $n \in \mathbb{N}^*$, we have

$$ Q^{\varepsilon}(u_n^k)(\nabla w_n^k) \ge \frac{1}{u_1^{k,n}} a_{21}^{\varepsilon}'(u_1^{k,n}) |\nabla u_1^{k,n}|^2 + \frac{1}{u_2^{k,n}} a_{12}^{\varepsilon}'(u_2^{k,n}) |\nabla u_2^{k,n}|^2. $$

Since $a_{ij}^\varepsilon(x) = a_{ij}(x) + \varepsilon x$, assumption **H2** leads to

$$ Q^{\varepsilon}(u_n^k)(\nabla w_n^k) \ge \frac{4}{(1-\alpha)^2} \left| \nabla[u_1^{k,n}]^{\frac{1-\alpha}{2}} \right|^2 + \frac{4}{(1-\alpha)^2} \left| \nabla[u_2^{k,n}]^{\frac{1-\alpha}{2}} \right|^2 + 4\varepsilon \left| \nabla\sqrt{u_1^{k,n}} \right|^2 + 4\varepsilon \left| \nabla\sqrt{u_2^{k,n}} \right|^2, \quad (37) $$

We also have, because of (28), the boundedness of $(u_n^k)_{n \in \mathbb{N}^*}$ in $L^1(\Omega)$. Since the asymptotics that we are studying is only w.r.t. $n$ (that is, $\varepsilon$ and $\tau$ are fixed), we see finally that $(\sqrt{u_1^{k,n}})_{n \in \mathbb{N}^*}$ and $(\sqrt{u_2^{k,n}})_{n \in \mathbb{N}^*}$ are bounded in $H^1(\Omega) \hookrightarrow L^{2^*}(\Omega)$, with $2^* = 2M/(M-2) > 2$, and hence $(u_n^k)_{n \in \mathbb{N}^*}$ is eventually bounded in some $L^p(\Omega)$ space with $1 < p < \infty$.

On the other hand, estimate (23) gives us

$$ Q^{\varepsilon}(u_n^k)(w_n^k) \ge 4 \left| \nabla \sqrt{a_{21}^{\varepsilon}(u_1^{k,n}) a_{12}^{\varepsilon}(u_2^{k,n})} \right|^2, $$

which, together with (27), leads to

$$ \left\| \nabla \sqrt{a_{21}^{\varepsilon}(u_1^{k,n}) a_{12}^{\varepsilon}(u_2^{k,n})} \right\|_{L^2(\Omega)}^2 \le \frac{1}{4\tau} K_T (\mathcal{E}_{\varepsilon}(u^0) + 1). $$
Using Poincaré-Wirtinger, Cauchy-Schwarz and Young inequalities, we get, for some constant $C_\Omega$ depending only on $\Omega$,

$$ \left\| \sqrt{a_{21}^{\epsilon}(u_1^{k,n}) a_{12}^{\epsilon}(u_2^{k,n})} \right\|_{L^2(\Omega)} \le C_{\Omega} \left[ \|a_{21}^{\epsilon}(u_1^{k,n})\|_{L^1(\Omega)} + \|a_{12}^{\epsilon}(u_2^{k,n})\|_{L^1(\Omega)} \right] + \frac{C_{\Omega}}{2\sqrt{\tau}} \sqrt{K_T(\mathcal{E}_{\epsilon}(u^0) + 1)}, $$

and, again because point (iv) of Lemma 3.7, for some constant $D$ ($\epsilon < 1$),

$$ a_{ji}^{\epsilon}(u_i^{k,n}) \le D(2 + \psi_i^{\epsilon}(u_i^{k,n})), $$

so we finally have

$$ \left\| \sqrt{a_{21}^{\epsilon}(u_1^{k,n}) a_{12}^{\epsilon}(u_2^{k,n})} \right\|_{L^2(\Omega)} \le C_{\Omega} \left[ 8D\mu(\Omega) + 2D\mathcal{E}_{\epsilon}(u_n^k) \right] + \frac{1}{2\sqrt{\tau}} \sqrt{K_T(\mathcal{E}_{\epsilon}(u^0) + 1)}, $$

and hence thanks to (27), the sequence $\left(\sqrt{a_{21}^{\epsilon}(u_1^{k,n}) a_{12}^{\epsilon}(u_2^{k,n})}\right)_{n \in \mathbb{N}^*}$ is bounded in $H^1(\Omega) \hookrightarrow L^{2^*}(\Omega)$, with $2^* = 2M/(M-2) > 2$. The previous continuous injection implies then that $(a_{21}^{\epsilon}(u_1^{k,n}) a_{12}^{\epsilon}(u_2^{k,n}))_{n \in \mathbb{N}^*}$ is bounded in some $L^p(\Omega)$ space with $1 < p < \infty$ and since $a_{ji}^{\epsilon}(u_i^{k,n}) = a_{ji}(u_i^{k,n}) + \epsilon u_i^{k,n}$, we eventually get that $(u_1^{k,n} u_2^{k,n})_{n \in \mathbb{N}^*}$, $(u_1^{k,n} a_{12}(u_2^{k,n}))_{n \in \mathbb{N}^*}$ and $(u_2^{k,n} a_{21}(u_1^{k,n}))_{n \in \mathbb{N}^*}$ are bounded in the same $L^p(\Omega)$.

Summing the bounds already obtained [we recall that they hold for a given $\epsilon$ and $\tau$], we see that (for all $k \in \{1, \dots, N\}$):

* $(\sqrt{u_1^{k,n}})_{n \in \mathbb{N}^*}$ and $(\sqrt{u_2^{k,n}})_{n \in \mathbb{N}^*}$ are bounded in $H^1(\Omega)$,

* $(\sqrt{a_{21}^{\epsilon}(u_1^{k,n}) a_{12}^{\epsilon}(u_2^{k,n})})_{n \in \mathbb{N}^*}$ is bounded in $H^1(\Omega)$,

* $(u_1^{k,n})_{n \in \mathbb{N}^*}$, $(u_2^{k,n})_{n \in \mathbb{N}^*}$, $(u_1^{k,n} u_2^{k,n})_{n \in \mathbb{N}^*}$, $(u_1^{k,n} a_{12}(u_2^{k,n}))_{n \in \mathbb{N}^*}$ and $(u_2^{k,n} a_{21}(u_1^{k,n}))_{n \in \mathbb{N}^*}$ are bounded in some $L^p(\Omega)$ space, with $1 < p < \infty$,

* and obviously $(w_n^k)_{n \in \mathbb{N}^*}$ is bounded in $H^1(\Omega)$, because of estimate (27).

Since we are only dealing with a finite number of values for $k \in \{1, \dots, N\}$ (at this point, the functions are not time-dependent), we shall only detail the study of $(u_n^1)_{n \in \mathbb{N}^*}$, the other values of $k$ being similar (one only has to extract a finite number of subsequences). For every test function $\chi = (\chi_1, \chi_2) \in V_n^2$, the weak formulation (31) may be written $(i \neq j \in \{1, 2\})$

$$ \frac{1}{\tau} \langle \chi_i, u_i^{1,n} \rangle_{L^2(\Omega)} - \frac{1}{\tau} \langle \chi_i, u_i^{0,n} \rangle_{L^2(\Omega)} - \langle \Delta \chi_i, [a_{ii}^{\epsilon}(u_i^{1,n}) + u_i^{1,n} a_{ij}(u_j^{1,n}) + \epsilon u_i^{1,n} u_j^{1,n}] \rangle_{L^2(\Omega)} \\ + \epsilon \langle \chi_i - \Delta \chi_i, w_i^{1,n} \rangle_{L^2(\Omega)} = \langle \chi_i, R_{ij}^{\epsilon}(u_i^{1,n}, u_j^{1,n}) \rangle. $$

First we extract (but do not change the index) a subsequence of $(w_n^1)_{n \in \mathbb{N}^*}$ converging in $L^2(\Omega)$ and almost everywhere to some element $w^1 \in H^1(\Omega)^2$. The almost everywhere convergence is transmitted to $(u_n^1)_{n \in \mathbb{N}^*}$, since $u_n^1 = (\psi'_\epsilon)^{-1}(w_n^1)$, the $\psi'_\epsilon$ functions being homeomorphisms. Because of the $L^p(\Omega)$ $(1 < p < \infty)$ bound for $(u_n^1)_{n \in \mathbb{N}^*}$, we can extract a subsequence (of the previous subsequence) converging weakly in $L^p(\Omega)$ to some function $v^1 \in L^1(\Omega)^2$. It follows by a classical argument that $v^1$ is almost everywhere equal to $u^1 := (\psi'_\epsilon)^{-1}(w^1)$, that is the almost everywhere limit of $(u_n^1)_{n \in \mathbb{N}^*}$.

The same argument holds for $(u_1^{1,n} u_2^{1,n})_{n \in \mathbb{N}^*}$ and $(u_i^{1,n} a_{ij}(u_j^{1,n}))_{n \in \mathbb{N}^*}$, $i \neq j = 1, 2$, since we have the same $L^p(\Omega)$ bound. We can also assume<fnref n="1" /> that $(u^{0,n})_{n \in \mathbb{N}^*}$ converges weakly in $L^p(\Omega)$ to some function $u^0 \in L^p(\Omega)$. Eventually, the cutoff perturbation introduced on $a_{ii}^\epsilon$ and the (superlinear) reaction terms ensures the weak convergence, in $L^p(\Omega)$, of $(a_{ii}^\epsilon(u_i^{1,n}))_{n \in \mathbb{N}^*}$ and $(R^\epsilon(u^{1,n}))_{n \in \mathbb{N}^*}$, respectively to $a_{ii}^\epsilon(u_i^1)$ and $R^\epsilon(u^1)$.

<footnote n="1">¹In fact we have $u^{0,n} = u^0$! This line is just to justify the handling of the same term in the other time steps.</footnote>
We now fix a test function $\chi$ in $\cup_{n \in \mathbb{N}^*} V_n^2$, say $\chi \in V_m^2$, $m \in \mathbb{N}^*$. Then, for $n$ greater than $m$, since the sequence of spaces $V_n$ is increasing, the previous weak formulation passes to the limit thanks to all the extractions that we made, and we get eq. (32).

It remains to prove the bounds announced in the Proposition.

Since we have almost everywhere convergence of $(u_n^k)_{n \in \mathbb{N}}$ and weak $H^1(\Omega)$ convergence of $(w_n^k)_{n \in \mathbb{N}}$, Fatou's Lemma and the classical estimate for the weak limits give directly, for all $k \in \{1, \dots, N\}$, estimate (33).

We then use (37) that we rewrite here:

$$Q^{\varepsilon}(u_n^k)(\nabla w_n^k) \geq \frac{4}{(1-\alpha)^2} \left|\nabla[u_1^{k,n}]^{\frac{1-\alpha}{2}}\right|^2 + \frac{4}{(1-\alpha)^2} \left|\nabla[u_2^{k,n}]^{\frac{1-\alpha}{2}}\right|^2 + 4\varepsilon \left|\nabla\sqrt{u_1^{k,n}}\right|^2 + 4\varepsilon \left|\nabla\sqrt{u_2^{k,n}}\right|^2,$$

We notice that for all $j \in \{1, \dots, k\}$, the sequence $([u_n^j]^{\frac{1-\alpha}{2}})_{n \in \mathbb{N}}$ was, precisely in view of the previous inequality, bounded in $H^1(\Omega)$, so that we can (adding another extraction) assume that we also had the convergence $([u_n^k]^{\frac{1-\alpha}{2}})_{n \in \mathbb{N}} \to [u^k]^{\frac{1-\alpha}{2}}$ in $H^1(\Omega)$ (using the uniqueness of the weak limit) and we hence have estimate (34) using Fatou's Lemma.

As for (28) and (29), we notice that we had a $L^p(\Omega)$ ($p > 1$) bound for both $(u^{j,n})_{n \in \mathbb{N}^*}$ and $(R^{-,\varepsilon}(u^{j,n}))_{n \in \mathbb{N}^*}$, associated with almost everywhere convergence. This is sufficient (using Egoroff's theorem) to get the (strong) convergence of these two sequences in $L^1(\Omega)$. Estimates (28) and (29) give hence (35), and (36) after taking the limit in $n$. $\square$

# 6 Duality Estimate

This (quite technical) Section is devoted to the Proof of an extra a priori estimate (named the duality estimate) for the solutions of the approximate problem obtained in Proposition 5.1. This estimate is analogous to estimate (11) in the simplified problem of Section 2.

The rather simple Proof of the duality Lemma presented in Section 2 cannot unfortunately be directly used, particularly because of the time discretization and the lack of regularity of the diffusion coefficients, and a specific statement and Proof has to be introduced.

## 6.1 Notations

**Definition 6.1.** For a given family $h := (h^k)_{1 \le k \le N}$ of functions defined on $\Omega$, we denote by $\underline{h}^\tau$ the step (in time) function defined on $\mathbb{R} \times \Omega$ by

$$\underline{h}^{\tau}(t, x) := \sum_{k=1}^{N} h^{k}(x) \mathbf{1}_{[k-1)\tau, k\tau]}(t).$$

We then have by definition, for all $p, q \in [1, \infty[$,

$$\|\underline{h}^{\tau}\|_{L^q([0,T];L^p(\Omega))} = \left( \sum_{k=1}^{N} \tau \|h^k\|_{L^p(\Omega)}^q \right)^{1/q},$$

and in particular

$$\|\underline{h}^{\tau}\|_{L^p(Q_T)} = \left( \sum_{k=1}^{N} \tau \int_{\Omega} |h^k(x)|^p dx \right)^{1/p}.$$

## 6.2 Duality estimate: abstract result

This Subsection is devoted to the establishment of a discretized version of the duality estimates devised for singular parabolic equations in [27], [25]. We start with the following standard Lemma, which is a consequence of Lax-Milgram Theorem together with the maximum principle.
**Lemma 6.2.** Consider a smooth function $b \in \mathcal{C}^\infty(\Omega)$ such that $b \ge \gamma > 0$ for some constant $\gamma$. For a given $\Psi \in L^2(\Omega)$, the variational problem associated to the equation

$$-\Phi + b\Delta\Phi = \Psi \qquad (38)$$

is well-posed in $H_\nu^2(\Omega)$ (with $L^2(\Omega)$ test functions). If $\Psi$ is assumed to be nonpositive, then the solution of the previous equation is nonnegative.

Using the previous Lemma, we obtain the following result, which can be seen as a particular instance of the so-called Rothe method (see [21] for a general framework). Since this method is quite standard, we skip the proof of this Lemma; It may nevertheless be found in detail in [13].

**Lemma 6.3.** Consider a real number $r > 0$ such that $1 - 2r\tau > 0$ and two families $b := (b^k)_{1\le k\le N}$, $F := (F^k)_{1\le k\le N}$ of $\mathcal{C}^\infty(\Omega)$ functions. Assume that for all $k \in \{1, \dots, N\}$, $b^k \ge 1$ and $F^k \le 0$ (pointwise).

Then there exists a family $\Phi := (\Phi^k)_{1\le k\le N} \in H_\nu^2(\Omega)^N$ of nonnegative functions such that, (defining $\Phi^{N+1} := 0$),

$$\forall k \in \{1, \dots, N\}, \quad \frac{\Phi^{k+1} - \Phi^k}{\tau} + b^k \Delta \Phi^k = \sqrt{b^k} F^k - r \Phi^k, \qquad (39)$$

where eq. (39) has to be understood weakly, against $L^2(\Omega)$ test functions. This family satisfies furthermore

$$\forall j \in \{1, \dots, N\}, \quad \|\nabla \Phi^j\|_{L^2(\Omega)} \le e^{r(\tau)T} \|\underline{F}^\tau\|_{L^2(Q_T)}, \qquad (40)$$

$$\| \sqrt{b^\tau} \Delta \Phi^\tau \|_{L^2(Q_T)} \le e^{r(\tau)T} \| \underline{F}^\tau \|_{L^2(Q_T)}, \qquad (41)$$

$$\left| \int_{\Omega} \Phi^j dx \right| \le \left[ \left\| \sqrt{b^\tau} \right\|_{L^2(Q_T)} + e^{r(\tau)T} \right] \left\| \underline{F}^\tau \right\|_{L^2(Q_T)} e^{r(\tau)T}, \qquad (42)$$

$$\forall j \in \{1, \dots, N\}, \quad \|\Phi^j\|_{H^1(\Omega)} \le e^{r(\tau)T} C_{\Omega} \left[ e^{r(\tau)T} + \|\sqrt{b^\tau}\|_{L^2(Q_T)} \right] \|\underline{F}^\tau\|_{L^2(Q_T)}, \qquad (43)$$

where the sequence of positive scalars $(r(\tau))_\tau$ goes to $r$ when $\tau \to 0$.

If we want to use functions of the previous family $\Phi$ as test functions in the weak formulation (32), we have to show that they actually belong to $W^{2,p'}(\Omega)$. This can be done using the following lemma

**Lemma 6.4.** Let $r$ be a positive real number. If $w \in H_\nu^2(\Omega)$ satisfies $w \ge 0$ and $-\Delta w \le f + rw$ almost everywhere, for some $f$ in $L^q(\Omega)$, $q > \max(M/2, 2)$, then

$$\|w\|_{L^\infty(\Omega)} \le C_{\Omega,q,r} (\|f\|_{L^q(\Omega)} + \|w\|_{L^1(\Omega)}).$$

*Proof.* $C_{\Omega,q}$ and $C_q$ denotes constants that may vary from line to line. Lemma 6.2 gives us the existence of $g \in H_\nu^2(\Omega)$, unique solution of

$$g - \Delta g = (r+1)w + f.$$

For $p \ge 1$ define $p^\bullet$ as the exponent of the Sobolev embedding $W^{2,p}(\Omega) \hookrightarrow L^{p^\bullet}(\Omega)$. Since $w \in H_\nu^2(\Omega)$, we have $w \in L^{2^\bullet}(\Omega)$. Let us show that $w \in L^q(\Omega)$ by an iteration argument. Indeed, if $2^\bullet < q$ then the r.h.s. of the previous equation lies in $L^{2^\bullet}(\Omega)$. By elliptic regularity of $(\mathrm{Id} - \Delta)$ and a Sobolev embedding (cf. [17]), we hence get

$$\|g\|_{L^{2^\bullet\bullet}(\Omega)} \le C_{\Omega} \|g\|_{W^{2,2^\bullet}(\Omega)} \le C_{\Omega} \|(r+1)w + f\|_{L^{2^\bullet}(\Omega)}.$$

But we also have

$$w - \Delta w \le (r+1)w + f,$$

that we can write

$$(w - g) - \Delta (w - g) \le 0,$$
so that by the weak maximum principle in $H_\nu^1(\Omega)$, we have

$$0 \le w \le g \in L^{2^{\bullet\bullet}}(\Omega).$$

Then $w \in L^{2^{\bullet\bullet}}(\Omega)$, and since the sequence $2^{\bullet\bullet\cdots}$ is strictly increasing, we eventually see that $w \in L^q(\Omega)$ by iteration. Using the elliptic regularity mentioned before, we hence have

$$\|g\|_{W^{2,q}(\Omega)} \le C_{\Omega}\|(r+1)w + f\|_{L^q(\Omega)}.$$

Since $q > M/2$, we get the $L^\infty$ estimate

$$\|g\|_{L^\infty(\Omega)} \le C_{\Omega,q}\|(r+1)w + f\|_{L^q(\Omega)},$$

that is again transmitted to $w$, thanks to the weak maximum principle :

$$\|w\|_{L^\infty(\Omega)} \le C_{\Omega,q}\|f\|_{L^q(\Omega)} + C_{\Omega,q}(r+1)\|w\|_{L^q(\Omega)}.$$

Since

$$\|w\|_{L^q(\Omega)} \le \|w\|_{L^1(\Omega)}^{1/q} \|w\|_{L^\infty(\Omega)}^{1/q'}$$

we get by Young's inequality

$$C_{\Omega,q}(r+1)\|w\|_{L^q(\Omega)} \le C_{q,r}\|w\|_{L^1(\Omega)} + \frac{1}{2}\|w\|_{L^\infty(\Omega)}$$

and we are eventually able to conclude that

$$\|w\|_{L^\infty(\Omega)} \le C_{\Omega,q,r} (\|f\|_{L^q(\Omega)} + \|w\|_{L^1(\Omega)}) .$$

□

Using the elliptic regularity of $(\mathrm{Id} - \Delta)$ and the two previous Lemmas, we are able to prove the following result (see [13] for more details):

**Lemma 6.5.** Under the assumptions of Lemma 6.3, the sequence of functions $\Phi$ satisfies in fact, for all $k \in \{1, \dots, N\}$

$$\|\Phi^k\|_{L^\infty(\Omega)} + \|\Delta\Phi^k\|_{L^\infty(\Omega)} \le C_{\tau,\Omega} \left[ 1 + \|\sqrt{b^\tau}\|_{L^2(Q_T)} \right] \|F^\tau\|_{L^\infty(Q_T)},$$

the constant $C_{\tau,\Omega}$ depending only on $\tau$ and $\Omega$, but (severely) blowing up as $\tau \to 0$. In particular we get that $\Phi^k \in W^{2,q}(\Omega)$ for all $q \in [1, \infty[$.

## 6.3 Duality estimate: application to the system

We now can present the main result of this Section:

**Proposition 6.6.** We consider the assumptions of Theorem 1.2 and keep in mind the Definitions 4.1 – 4.2, and Definition 3.6. Then, any weak solution of eq. (32) satisfies the following estimate:

$$\left\| \underline{u}_1^\tau \sqrt{d_{11}^\varepsilon(\underline{u}_1^\tau) + a_{12}^\varepsilon(\underline{u}_2^\tau)} \right\|_{L^2(Q_T)} \le D_T(u^0), \quad (44)$$

where we denote by $\underline{u}_1^\tau$ and $\underline{u}_2^\tau$ the step (in time) functions (see Definition 6.1) associated to the two families $(u_1^k)_{1 \le k \le N}$ and $(u_2^k)_{1 \le k \le N}$ (excluding hence $u_1^0$ and $u_2^0$), and where $D_T(u^0)$ is a constant depending only on $T, \Omega, \|u^0\|_{L^2(\Omega)}$, and the data of the equation $(r_i, d_{ii}, a_{ij}, s_{ij})$.

The same estimate holds when the subscripts 1 and 2 are exchanged.
Proof. Since the original system is symmetric in $u_1$ and $u_2$, let us focus on $u_1$ for now. As in Lemma 6.3, we consider two families $b := (b^k)_{1\le k\le N}$ and $F := (F^k)_{1\le k\le N}$ of $\mathcal{C}^\infty(\Omega)$ functions such that $b^k \ge 1$ and $F^k \le 0$, and the associated sequence $\Phi := (\Phi^k)_{1\le k\le N}$, with $r := r_1$. As shown in Lemma 6.5, the functions $\Phi^k$ belong to $W^{2,p'}(\Omega)$ (for all $p > 1$) and are hence admissible in the weak formulation (32). We may therefore write, for all $k \in \{1, \dots, N\}$, taking $\Phi^k$ for test function in the weak formulation on $u_1^k$:

$$ \frac{1}{\tau} \langle \Phi^k, u_1^k \rangle_{L^2(\Omega)} - \frac{1}{\tau} \langle \Phi^k, u_1^{k-1} \rangle_{L^2(\Omega)} - \left\langle \Delta \Phi^k, \left[ a_{11}^{\varepsilon}(u_1^k) + u_1^k a_{12}^{\varepsilon}(u_2^k) \right] \right\rangle_{L^2(\Omega)} \\ + \varepsilon \langle \Phi^k - \Delta \Phi^k, w_1^k \rangle_{L^2(\Omega)} = \langle \Phi^k, R_{12}^{\varepsilon}(u_1^k, u_2^k) \rangle. $$

Thanks to the regularity of $\Phi^k$, one can take $u_1^k$ as a test function in the $k$-th equation (39), so that

$$ \frac{1}{\tau} \langle \Phi^{k+1}, u_1^k \rangle_{L^2(\Omega)} - \frac{1}{\tau} \langle \Phi^k, u_1^k \rangle_{L^2(\Omega)} + \langle \Delta \Phi^k, b^k u_1^k \rangle_{L^2(\Omega)} = \langle \sqrt{b^k} F^k, u_1^k \rangle_{L^2(\Omega)} - r_1 \langle \Phi^k, u_1^k \rangle_{L^2(\Omega)}. $$

We hence have, for all $k \in \{1, \dots, N\}$, adding the two previous equations,

$$ \frac{1}{\tau} \langle \Phi^{k+1}, u_1^k \rangle_{L^2(\Omega)} - \frac{1}{\tau} \langle \Phi^k, u_1^{k-1} \rangle_{L^2(\Omega)} + \langle \Delta \Phi^k, [b^k - d_{11}^{\varepsilon}(u_1^k) - a_{12}^{\varepsilon}(u_2^k)] u_1^k \rangle_{L^2(\Omega)} \\ + \varepsilon \langle \Phi^k - \Delta \Phi^k, w_1^k \rangle_{L^2(\Omega)} = \langle \sqrt{b^k} F^k, u_1^k \rangle_{L^2(\Omega)} - \langle \Phi^k, R_{12}^{-\varepsilon}(u_1^k, u_2^k) \rangle, $$

recalling that $a_{12}^{\varepsilon}(x) := a_{12}(x) + \varepsilon x$, $a_{11}^{\varepsilon}(x) := x d_{11}^{\varepsilon}(x)$, $d_{11}^{\varepsilon} := \gamma_{\varepsilon}(d_{11})$, and the decomposition $R_{12}^{\varepsilon} = R_{12}^{+} - R_{12}^{-\varepsilon}$. Since $\Phi^k$ is nonnegative, we get, if we denote $c^k := b^k - d_{11}^{\varepsilon}(u_1^k) - a_{12}^{\varepsilon}(u_2^k)$, after summing over $k \in \{1, \dots, N\}$ (and recalling that $\Phi^{N+1} = 0$),

$$ -\frac{1}{\tau} \langle \Phi^1, u_1^0 \rangle_{L^2(\Omega)} + \sum_{k=1}^{N} \langle \Delta \Phi^k, c^k u_1^k \rangle_{L^2(\Omega)} + \sum_{k=1}^{N} \varepsilon \langle \Phi^k, w_1^k \rangle_{H^1(\Omega)} \le \sum_{k=1}^{N} \langle \sqrt{b^k} F^k, u_1^k \rangle_{L^2(\Omega)}. \quad (45) $$

The following approximation Lemma can easily be deduced from standard measure theory results (for the measures $(1 + u_1^k) dx$), the prove of which may be found in [13], and will help us to handle the sequence $c := (c^k)_{1\le k\le N}$:

**Lemma 6.7.** There exists a sequence of families $(b_m)_{m \in \mathbb{N}} := (b_m^1, b_m^2, \dots, b_m^N)_{m \in \mathbb{N}}$ of $\mathcal{C}^\infty(\Omega)$ functions such that

$$ \forall (k, m) \in \{1, \dots, N\} \times \mathbb{N}, \quad b_m^k \ge 1, \quad (46) $$

$$ \forall k \in \{1, \dots, N\}, \quad \int_{\Omega} |b_m^k - d_{11}^{\varepsilon}(u_1^k) - a_{12}^{\varepsilon}(u_2^k)| (1 + u_1^k) dx \xrightarrow{m \to +\infty} 0, \quad (47) $$

$$ \forall k \in \{1, \dots, N\}, \quad \int_{\Omega} \left| \sqrt{b_m^k} - \sqrt{d_{11}^{\varepsilon}(u_1^k) + a_{12}^{\varepsilon}(u_2^k)} \right| (1 + u_1^k) dx \xrightarrow{m \to +\infty} 0, \quad (48) $$

$$ \left\| \sqrt{\underline{b}_m^\tau} \right\|_{L^2(Q_T)}^2 \xrightarrow{m \to +\infty} \sum_{k=1}^N \tau \int_{\Omega} |d_{11}^{\varepsilon}(u_1^k) + a_{12}^{\varepsilon}(u_2^k)| dx = \left\| \sqrt{d_{11}^{\varepsilon}(\underline{u}_1^\tau) + a_{12}^{\varepsilon}(\underline{u}_2^\tau)} \right\|_{L^2(Q_T)}^2. \quad (49) $$

We now fix a family $F := (F^k)_{1\le k\le N}$ of $\mathcal{C}^\infty(\Omega)$ nonpositive functions, and for each family $b_m$ of the sequence $(b_m)_{m \in \mathbb{N}}$ defined in Lemma 6.7, we define the corresponding family $\Phi_m := (\Phi_m^1, \Phi_m^2, \dots, \Phi_m^N)$, using Lemma 6.3. The previous estimate (45) can now be written

$$ -\langle \Phi_m^1, u_1^0 \rangle_{L^2(\Omega)} + \sum_{k=1}^{N} \tau \langle \Delta \Phi_m^k, c_m^k u_1^k \rangle_{L^2(\Omega)} + \sum_{k=1}^{N} \tau \varepsilon \langle \Phi_m^k, w_1^k \rangle_{H^1(\Omega)} \le \sum_{k=1}^{N} \tau \langle \sqrt{b_m^k} F^k, u_1^k \rangle_{L^2(\Omega)}, $$
where $c_m^k := b_m^k - d_{11}^\epsilon(u_1^k) - a_{12}^\epsilon(u_2^k)$. Since the right-hand side is nonpositive, we may write

$$ \mathbf{1}_m := \left| \sum_{k=1}^{N} \tau \langle \sqrt{b_m^k} F^k, u_1^k \rangle_{L^2(\Omega)} \right| \le \underbrace{\langle \Phi_m^1, u_1^0 \rangle_{L^2(\Omega)}}_{\mathbf{2}_m} + \underbrace{\sum_{k=1}^{N} \tau \left| \langle \Delta \Phi_m^k, c_m^k u_1^k \rangle_{L^2(\Omega)} \right|}_{\mathbf{3}_m} + \underbrace{\sum_{k=1}^{N} \tau \epsilon \left| \langle \Phi_m^k, u_1^k \rangle_{H^1(\Omega)} \right|}_{\mathbf{4}_m}. $$

Using (43) of Lemma 6.3 and (49), we get (for a given $\tau, \epsilon$),

$$ \overline{\lim}_{m \to +\infty} |\mathbf{2}_m| \le C_\Omega \|u^0\|_{L^2(\Omega)} \left[ e^{r_1(\tau)T} + \left\| \sqrt{d_{11}^\epsilon(\underline{u}_1^\tau)} + a_{12}^\epsilon(\underline{u}_2^\tau) \right\|_{L^2(Q_T)} \right] \|\underline{F}^\tau\|_{L^2(Q_T)} e^{r_1(\tau)T}. $$

Using Lemma 6.5, we get on the other hand

$$ \begin{aligned} |\mathbf{3}_m| &\le \tau N \sup_{1 \le k \le N} \left\{ \|\Delta \Phi_m^k\|_{L^\infty(\Omega)} \|c_m^k u_1^k\|_{L^1(\Omega)} \right\} \\ &\le \tau N C_{\tau, \Omega} \left[ 1 + \left\| \sqrt{\underline{b}_m^\tau} \right\|_{L^2(Q_T)} \right] \|\underline{F}^\tau\|_{L^\infty(Q_T)} \sup_{1 \le k \le N} \|c_m^k u_1^k\|_{L^1(\Omega)} \xrightarrow{m \to +\infty} 0, \end{aligned} $$

using (47) and (49) for the convergence. In order to treat $\mathbf{4}_m$, notice that we know, from inequality (33):

$$ \varepsilon \tau \sum_{k=1}^{N} \|w_1^k\|_{H^1(\Omega)}^2 \le K_T(\mathcal{E}_\varepsilon(u^0) + 1), $$

hence, using again (43) of Lemma 6.3, and (49) we get $(\tau N = T)$

$$ \overline{\lim}_{m \to +\infty} |\mathbf{4}_m| \le e^{r_1(\tau)T} C_\Omega \left[ e^{r_1(\tau)T} + \left\| \sqrt{d_{11}^\epsilon(\underline{u}_1^\tau)} + a_{12}^\epsilon(\underline{u}_2^\tau) \right\|_{L^2(Q_T)} \right] \|\underline{F}^\tau\|_{L^2(Q_T)} \sqrt{\varepsilon T} \sqrt{K_T(\mathcal{E}_\varepsilon(u^0) + 1)}. $$

Finally, because of (48), we can see that $\mathbf{1}_m$ converges, as $m$ goes to $+\infty$, to

$$ \left| \sum_{k=1}^{N} \tau \int_{\Omega} u_1^k \sqrt{d_{11}^\epsilon(u_1^k) + a_{12}^\epsilon(u_2^k)} F^k dx \right| = \left| \int_{Q_T} \underline{u}_1^\tau \sqrt{d_{11}^\epsilon(\underline{u}_1^\tau) + a_{12}^\epsilon(\underline{u}_2^\tau)} \underline{F}^\tau dx dt \right|. $$

All the previous estimates give hence, denoting $\underline{h}^\tau := \underline{u}_1^\tau \sqrt{d_{11}^\epsilon(\underline{u}_1^\tau) + a_{12}^\epsilon(\underline{u}_2^\tau)}$,

$$ \begin{aligned} \left| \int_{Q_T} \underline{h}^\tau \underline{F}^\tau dx dt \right| &\le e^{r_1(\tau)T} C_\Omega \left[ e^{r_1(\tau)T} + \left\| \sqrt{d_{11}^\epsilon(\underline{u}_1^\tau)} + a_{12}^\epsilon(\underline{u}_2^\tau) \right\|_{L^2(Q_T)} \right] \\ &\quad \times \left[ \sqrt{\varepsilon T} \sqrt{K_T(\mathcal{E}_\varepsilon(u^0) + 1)} + \|u^0\|_{L^2(\Omega)} \right] \|\underline{F}^\tau\|_{L^2(Q_T)}. \end{aligned} $$

Since $\underline{h}^\tau$ is a step (in time) nonnegative function and the previous holds true for all nonpositive smooth (in $x$) step (in time) functions, we have then by duality (switching to $L^1$ norms)

$$ \| \underline{u}_1^2 \left[ d_{11}^\epsilon(\underline{u}_1^\tau) + a_{12}^\epsilon(\underline{u}_2^\tau) \right] \|_{L^1(Q_T)} \le D_T(u^0) \left[ 1 + \| d_{11}^\epsilon(\underline{u}_1^\tau) \|_{L^1(\Omega)} + \| a_{12}^\epsilon(\underline{u}_2^\tau) \|_{L^1(Q_T)} \right], $$
where the constant $D_T(u^0)$ only depends on $T, \Omega$, $\|u^0\|_{L^2(\Omega)}$ and the data of the equation (we used here $\mathcal{E}_\varepsilon(u^0) \le C + \|u^0\|_{L^2(\Omega)}^2$, see Remark 3.9). Because of the concavity of $a_{12}^\varepsilon$ (and since it vanishes in 0), we have, with some constant $C_1, C_2$ **not depending** on $\varepsilon$,

$$a_{12}^{\varepsilon}(x) \le C_1 + C_2 \psi_2^{\varepsilon}(x),$$

so that, from (33), we have

$$\| \underline{u}_1^{\tau^2} [d_{11}^{\varepsilon}(\underline{u}_1^{\tau}) + a_{12}^{\varepsilon}(\underline{u}_2^{\tau})] \|_{L^1(Q_T)} \le D_T(u^0) \left[ 1 + \|d_{11}^{\varepsilon}(\underline{u}_1^{\tau})\|_{L^1(\Omega)} + K_T(\mathcal{E}_{\varepsilon}(u^0) + 1) \right],$$

that we can write (changing the definition of $D_T(u^0)$)

$$\| \underline{u}_1^{\tau^2} [d_{11}^{\varepsilon}(\underline{u}_1^{\tau}) + a_{12}^{\varepsilon}(\underline{u}_2^{\tau})] \|_{L^1(Q_T)} \le D_T(u^0) \left[ 1 + \|d_{11}^{\varepsilon}(\underline{u}_1^{\tau})\|_{L^1(\Omega)} \right].$$

But we have then, since $d_{11}^\varepsilon \le d_{11}$ which is a nondecreasing function, (simply denoting $D_T(u^0)$ by $D$)

$$\begin{aligned} \| \underline{u}_1^{\tau^2} [d_{11}^{\varepsilon}(\underline{u}_1^{\tau}) + a_{12}^{\varepsilon}(\underline{u}_2^{\tau})] \|_{L^1(Q_T)} &\le D \left[ 1 + \| \mathbf{1}_{\underline{u}_1^{\tau} \ge \sqrt{2D}} d_{11}^{\varepsilon}(\underline{u}_1^{\tau}) \|_{L^1(\Omega)} + \| \mathbf{1}_{\underline{u}_1^{\tau} < \sqrt{2D}} d_{11}^{\varepsilon}(\underline{u}_1^{\tau}) \|_{L^1(\Omega)} \right] \\ &\le D \left[ 1 + \frac{1}{2D} \| \underline{u}_1^{\tau^2} d_{11}^{\varepsilon}(\underline{u}_1^{\tau}) \|_{L^1(\Omega)} + d_{11}(\sqrt{2D})\mu(\Omega) \right], \end{aligned}$$

so that

$$\| \underline{u}_1^{\tau^2} \left[ \frac{1}{2} d_{11}^{\varepsilon}(\underline{u}_1^{\tau}) + a_{12}^{\varepsilon}(\underline{u}_2^{\tau}) \right] \|_{L^1(Q_T)} \le D \left[ 1 + d_{11}(\sqrt{2D})\mu(\Omega) \right],$$

that is to say, eq. (44). □

# 7 Proof of the main Theorem

This section is dedicated to the Proof of our main Theorem. We recall that at this point, we have thanks to Proposition 5.1 a family of approximate solutions (corresponding to (1) – (4)) satisfying the entropy estimates (33), (34) and “easy” estimates (35), (36). This family also satisfies the duality estimate (44) of Proposition 6.6.

The Proof of our main Theorem consists in relaxing the time discretization and in letting the regularization parameter $\varepsilon$ go to 0, in the approximate problem obtained in Proposition 5.1. This is done thanks to the use of the *a priori* estimates (33) – (36). This Proof is analogous to the Proof of Proposition 2.3 in the simplified case of Section 2.

*Proof.* Thanks to Proposition 5.1, we can consider a solution $\underline{u}_1^{\tau^\varepsilon}, \underline{u}_2^{\tau^\varepsilon}$ to system (32). For the sake of clarity, we shall denote this solution by $\underline{u}_1^\tau, \underline{u}_2^\tau$. Let us keep in mind that since we are dealing with the limit $(\tau, \varepsilon) \to (0, 0)$, we shall use bounds that are uniform w.r.t. these two parameters. Hence in the sequel, if not mentioned, the term “bounded” will always have to be understood as “uniformly w.r.t. $\tau$ and $\varepsilon$”. We shall only work with $\underline{u}_1^\tau$, the study of $\underline{u}_2^\tau$ being exactly identical.

We first establish compactness in the $x$ variable. Since for any smooth function $v$ we have

$$\frac{1-\alpha}{2} v^{-\frac{\alpha+1}{2}} \nabla v = \nabla v^{\frac{1-\alpha}{2}},$$

using (34), we have for all $k \in \{1, \dots, N\}$:

$$\tau \sum_{j=1}^{k} \int_{\Omega} |u_{1}^{j}|^{-\alpha-1} |\nabla u_{1}^{j}|^{2} dx \leq K_{T}(\mathcal{E}_{\varepsilon}(u^{0}) + 1).$$
We know that

$$|\nabla \underline{u}_1^\tau|(t, x) = \sum_{k=1}^{N} \mathbf{1}_{[(k-1)\tau, k\tau]}(t) |u_1^k|^{\frac{\alpha+1}{2}} |u_1^k|^{-\frac{\alpha+1}{2}} |\nabla u_1^k|.$$

Hence, using the dual estimate (44), we see that $\nabla \underline{u}_1^\tau$ is bounded in $L^{q_\alpha}(Q_T)$, where $q_\alpha$ is defined by the equality $(\alpha \in ]0, 1[)$

$$\frac{\alpha+1}{4} + \frac{1}{2} = \frac{1}{q_\alpha},$$

that is $q_\alpha = 8/(2\alpha + 6) > 1$.

As for the behavior in the time variable, let us denote by $\sigma_\tau$ the translation operator $g(\cdot, x) \mapsto g(\cdot - \tau, x)$, and introduce the space $E_\tau := L^1([\tau, T]; W^{2,\infty}(\Omega))$, having hence

$$\|\underline{u}_1^\tau - \sigma_\tau \underline{u}_1^\tau\|_{E_\eta} = \tau \sum_{k=2}^{N} \|u_1^k - u_1^{k-1}\|_{W^{2,\infty}(\Omega)}.$$

Now recall that, for all $1 \le k \le N$,

$$\partial_\tau u_1^k - \Delta([d_{11}^\varepsilon(u_1^k) + a_{12}^\varepsilon(u_2^k)]u_1^k) + \varepsilon(w_1^k - \Delta w_1^k) = R_{12}^\varepsilon(u^k),$$

weakly. We then have, since $H^1(\Omega) \hookrightarrow H^{-1}(\Omega) \hookrightarrow W^{2,\infty}(\Omega)$, and $L^1(\Omega) \hookrightarrow W^{2,\infty}(\Omega)$

$$\frac{1}{\tau} \|u_1^k - u_1^{k-1}\|_{W^{2,\infty}(\Omega)} \le \|[d_{11}^\varepsilon(u_1^k) + a_{12}^\varepsilon(u_2^k)]u_1^k\|_{L^1(\Omega)} + \varepsilon C_\Omega \|w_1^k\|_{H^1(\Omega)} + \|R_{12}^\varepsilon(u^k)\|_{L^1(\Omega)}$$

from which we deduce easily, multiplying by $\tau$ and summing up over $2 \le k \le N$

$$\frac{1}{\tau} \| \underline{u}_1^\tau - \sigma_\tau \underline{u}_1^\tau \|_{E_\tau} \le \| [d_{11}^\varepsilon(\underline{u}_1^\tau) + a_{12}^\varepsilon(\underline{u}_2^\tau)] \underline{u}_1^\tau \|_{L^1([\tau, T] \times \Omega)} + C_\Omega \varepsilon \tau \sum_{k=2}^N \|w_1^k\|_{H^1(\Omega)} + \tau \sum_{k=2}^N \|R_{12}^\varepsilon(u^k)\|_{L^1(\Omega)}. \quad (50)$$

The two last terms in the r.h.s. of (50) are bounded (uniformly in $\tau, \varepsilon$) thanks to estimates (33), (35) and (36). Let us handle the remaining term of the r.h.s of (50). with the help of the duality estimate. We have (with $\| \cdot \|$ meaning the $L^1([\tau, T] \times \Omega)$ norm)

$$\begin{aligned} \| [d_{11}^\varepsilon(\underline{u}_1^\tau) + a_{12}^\varepsilon(\underline{u}_2^\tau)] \underline{u}_1^\tau \| &\le \| [d_{11}^\varepsilon(\underline{u}_1^\tau) + a_{12}^\varepsilon(\underline{u}_2^\tau)] \mathbf{1}_{\underline{u}_1^\tau \le 1} \| + \| [d_{11}^\varepsilon(\underline{u}_1^\tau) + a_{12}^\varepsilon(\underline{u}_2^\tau)] \underline{u}_1^{\tau 2} \| \\ &\le \mu(\Omega) \sup_{x \in [0,1]} d_{11}^\varepsilon(x) + \|a_{12}^\varepsilon(u_2^{k+1})\|_{L^1(\Omega)} + \| [d_{11}^\varepsilon(\underline{u}_1^\tau) + a_{12}^\varepsilon(\underline{u}_2^\tau)] \underline{u}_1^{\tau 2} \|, \end{aligned}$$

so that using $d_{11}^\varepsilon \le d_{11}$, $a_{12}^\varepsilon \le D(2 + \psi_2^\varepsilon)$ (see Lemma 3.7) and estimates (33) and (44), we get the boundedness of the r.h.s. of (50).

At the end of the day, we obtained

* $(\underline{u}_1^\tau)_\tau$ is bounded in $L^{q_\alpha}([0, T]; W^{1, q_\alpha}(\Omega))$,

* $\|\sigma_\tau \underline{u}_1^\tau - \underline{u}_1^\tau\|_{E_\tau} \le C\tau$, with independent of $\tau$ and $\varepsilon$.

We know that $W^{1,q_\alpha}(\Omega) \hookrightarrow L^r(\Omega) \hookrightarrow W^{2,\infty}(\Omega)$, where the first injection is compact, the second one is continuous, and $r$ denotes any real number of the interval $[1, q_\alpha^*[$. We hence can use Theorem 1 of [15] to get the compactness of $(\underline{u}_1^\tau)_\tau$ in $L^1([0, T]; L^r(\Omega))$.

Let us now check that the weak formulation passes to the limit. It may be easily checked that

$$\frac{\underline{u}_1^\tau - \sigma_\tau \underline{u}_1^\tau}{\tau}(t, x) = \sum_{k=1}^{N} \mathbf{1}_{[(k-1)\tau, k\tau]}(t) \partial_\tau u_1^k(x) + \frac{1}{\tau} \mathbf{1}_{[0,\tau]}(t) u_1^0(x) - \frac{1}{\tau} \mathbf{1}_{[T,T+\tau]}(t) u_1^N(x),$$
where $\partial_{\tau}$ is the finite difference operator (see Definition 4.1). Now let us fix a test function $\theta \in \mathcal{D}([0, T]; \mathcal{C}_{\nu}^{\infty}(\Omega))$. We have hence

$$ \sum_{k=1}^{N} \int_{(k-1)\tau}^{k\tau} \langle \partial_{\tau} u_{1}^{k}, \theta(t) \rangle \, dt = \frac{1}{\tau} \int_{0}^{\tau} \langle u_{1}^{0}, \theta(t) \rangle \, dt - \int_{\mathbb{R}} \int_{\Omega} \frac{u_{1}^{\tau} - \sigma_{\tau} u_{1}^{\tau}}{\tau} \cdot \theta(t, x) \, dx \, dt, $$

where we denote by $\langle \cdot, \cdot \rangle$ the duality bracket $\mathcal{D}'(\Omega)/\mathcal{D}(\Omega)$ (which is simply the integration on $\Omega$ here ...) and a change of variable leads to

$$ \sum_{k=1}^{N} \int_{(k-1)\tau}^{k\tau} \langle \partial_{\tau} u_{1}^{k}, \theta(t) \rangle \, dt + \frac{1}{\tau} \int_{0}^{\tau} \langle u_{1}^{0}, \theta(t) \rangle \, dt = \int_{0}^{T} \int_{\Omega} \underline{u}_{1}^{\tau} \cdot \frac{\theta - \sigma_{-\tau}\theta}{\tau}(t, x) \, dx \, dt. \quad (51) $$

Since for all $k \in \{1, \dots, N\}$ we have (in the weak sense)

$$ \partial_{\tau} u_{1}^{k} - \Delta \left( \left[ d_{11}^{\varepsilon}(u_{1}^{k}) + a_{12}^{\varepsilon}(u_{2}^{k}) \right] u_{1}^{k} \right) + \varepsilon(w_{1}^{k} - \Delta w_{1}^{k}) = R_{12}^{\varepsilon}(u^{k}), $$

we also can write, for all $t \in [0, T[$,

$$ \langle \partial_{\tau} u_{1}^{k}, \theta(t) \rangle - \langle a_{11}^{\varepsilon}(u_{1}^{k}) + u_{1}^{k} a_{12}(u_{2}^{k}) + \varepsilon u_{1}^{k} u_{2}^{k}, \Delta \theta(t) \rangle + \varepsilon \langle w_{1}^{k}, \theta(t) - \Delta \theta(t) \rangle = \langle R_{12}^{\varepsilon}(u^{k}), \theta(t) \rangle, $$

so that integrating on $](k-1)\tau, k\tau]$ and summing over $k \in \{1, \dots, N\}$, we get

$$ \sum_{k=1}^{N} \int_{(k-1)\tau}^{k\tau} \langle \partial_{\tau} u_{1}^{k}, \theta(t) \rangle \, dt - \int_{0}^{T} \int_{\Omega} \left[ a_{11}^{\varepsilon}(\underline{u}_{1}^{\tau}) + \underline{u}_{1}^{\tau} a_{12}(\underline{u}_{2}^{\tau}) + \varepsilon \underline{u}_{1}^{\tau} \underline{u}_{2}^{\tau} \right] \cdot \Delta \theta(t, x) \, dx \, dt \\ + \varepsilon \sum_{k=1}^{N} \int_{(k-1)\tau}^{k\tau} \langle w_{1}^{k}, \theta(t) - \Delta \theta(t) \rangle \, dt = \int_{0}^{T} \int_{\Omega} R_{12}^{\varepsilon}(\underline{u}^{\tau})(t, x) \theta(t, x) \, dx \, dt. $$

Using (51), we eventually get

$$ \int_{0}^{T} \int_{\Omega} \underline{u}_{1}^{\tau} \cdot \frac{\theta - \sigma_{-\tau}\theta}{\tau}(t, x) \, dx \, dt - \int_{0}^{T} \int_{\Omega} \left[ a_{11}^{\varepsilon}(\underline{u}_{1}^{\tau}) + \underline{u}_{1}^{\tau} a_{12}(\underline{u}_{2}^{\tau}) + \varepsilon \underline{u}_{1}^{\tau} \underline{u}_{2}^{\tau} \right] \cdot \Delta \theta(t, x) \, dx \, dt \\ + \varepsilon \sum_{k=1}^{N} \int_{(k-1)\tau}^{k\tau} \langle w_{1}^{k}, \theta(t) - \Delta \theta(t) \rangle \, dt \\ = \frac{1}{\tau} \int_{0}^{\tau} \langle u_{1}^{0}, \theta(t) \rangle \, dt + \int_{0}^{T} \int_{\Omega} R_{12}^{\varepsilon}(\underline{u}^{\tau})(t, x) \theta(t, x) \, dx \, dt. \quad (52) $$

We now can study the limit $(\tau, \varepsilon) \to (0, 0)$. Since $\theta$ is smooth, $(\theta - \sigma_{-\tau}\theta)\tau^{-1}$ uniformly converges to $-\partial_{t}\theta$ as $\tau \to 0$, and the first term in the third line of (52) goes to $\langle u_{1}^{0}, \theta(0) \rangle$. The term in the second line of (52) is going to 0 with $(\tau, \varepsilon)$ because thanks to (33)

$$ \left| \varepsilon \sum_{k=1}^{N} \int_{(k-1)\tau}^{k\tau} \langle w_{1}^{k}, \theta(t) - \Delta\theta(t) \rangle \, dt \right| \le \|\theta - \Delta\theta\|_{L^{\infty}([0,T]; L^{2}(\Omega))} \sqrt{\varepsilon\sqrt{T}K_{T}(\mathcal{E}_{\varepsilon}(u^{0}) + 1)} \underset{(\tau, \varepsilon) \to (0,0)}{\longrightarrow} 0. $$

We also have

$$ \left| \int_{0}^{T} \int_{\Omega} \varepsilon \underline{u}_{1}^{\tau} \underline{u}_{2}^{\tau} \Delta \theta(t, x) \, dx \, dt \right| \le \varepsilon \| \underline{u}_{1}^{\tau} \|_{L^{2}(Q_{T})} \| \underline{u}_{2}^{\tau} \|_{L^{2}(Q_{T})} \| \Delta \theta \|_{L^{\infty}(Q_{T})}, $$

which goes to zero with $(\tau, \varepsilon)$ thanks to the dual estimate (44), so that it remains just to handle the nonlinearities in order to pass to the limit in (52).

Thanks to the compactness result that we proved above, we get the existence of $u := (u_1, u_2) \in L^1([0, T]; L^r(\Omega))$ (with $r < q_{\alpha}^*$) such that, up to a subsequence, we have

$$ (\underline{u}_1^\tau)_{\tau, \varepsilon} \underset{(\tau, \varepsilon) \to (0,0)}{\longrightarrow} u_1, $$

$$ (\underline{u}_2^\tau)_{\tau, \varepsilon} \underset{(\tau, \varepsilon) \to (0,0)}{\longrightarrow} u_2, $$
in L¹([0, T]; L²(Ω)), and also almost everywhere on $Q_T$. Because of the dual estimate (44) $(\underline{u}_1^\tau \sqrt{a_{12}(\underline{u}_2^\tau)})_{\tau,\varepsilon}$ is bounded in L²($Q_T$). Using assumption **H2**, we see that $a_{12}$ is at most linearly growing, so that using again the dual estimate (44) (but inverting the subscripts 1 and 2), $(a_{12}(\underline{u}_2^\tau))_{\tau,\varepsilon}$ is bounded in L²($Q_T$). Writing

$$ \underline{u}_1^\tau a_{12}(\underline{u}_2^\tau) = \underline{u}_1^\tau \sqrt{a_{12}(\underline{u}_2^\tau)} \sqrt{a_{12}(\underline{u}_2^\tau)}, $$

we see eventually that $(\underline{u}_1^\tau a_{12}(\underline{u}_2^\tau))_{\tau,\varepsilon}$ is bounded in L⁴/³($Q_T$) and we may thus extract a subsequence converging weakly in this space, and whose limit has to be equal to $u_1 a_{12}(u_2)$ (because of the previous almost everywhere convergence).

As for the self-diffusion, the dual estimate (44) ensures that $(\underline{u}_1^\tau \sqrt{d_{11}^\varepsilon(\underline{u}_1^\tau)})_{\tau,\varepsilon}$ is bounded in L²($Q_T$). Since $d_{11}^\varepsilon$ is a nondecreasing function, we infer the boundedness of $(a_{11}^\varepsilon(\underline{u}_1^\tau))_{\tau,\varepsilon} = (\underline{u}_1^\tau d_{11}^\varepsilon(\underline{u}_1^\tau))_{\tau,\varepsilon}$ in L¹($Q_T$), so that the weak convergence of this sequence can be deduced from its uniform integrability (Dunford-Pettis theorem). Since the sequence is bounded in L¹($Q_T$), it remains just to prove that $(E \subset Q_T)$

$$ \sup_{\tau, \varepsilon} \int_E \underline{u}_1^\tau d_{11}^\varepsilon(\underline{u}_1^\tau) \xrightarrow{|E| \to 0} 0. $$

For this purpose, just write

$$ \sup_{\tau, \varepsilon} \int_E \underline{u}_1^\tau d_{11}^\varepsilon(\underline{u}_1^\tau) \le \frac{1}{Z} \sup_{\tau, \varepsilon} \int_{Q_T} \underline{u}_1^{\tau^2} d_{11}^\varepsilon(\underline{u}_1^\tau) + d_{11}(Z) \sup_{\tau, \varepsilon} \int_E \underline{u}_1^\tau, $$

where in the r.h.s. the first term vanishes with 1/Z (because of estimate (44)) and the second one vanishes with $|E|$ (with Z fixed) because of the uniform integrability of $(\underline{u}_1^\tau)_{\tau,\varepsilon}$ (bounded in L²($Q_T$) thanks again to estimate (44)).

For the reactions terms, since the nonlinearities are always strictly sublinear or dominated by the self diffusion, one easily manages (using the dual estimate) to use Dunford Pettis criterion. Indeed, first, since $(\gamma_\varepsilon)_\varepsilon$ is increasing to the identity, the uniform integrability of $R_{12}^{\varepsilon}(\underline{u}^\tau)$ reduces to check this property for both $s_{11}(\underline{u}_1^\tau)\underline{u}_1^\tau$ and $s_{12}(\underline{u}_2^\tau)\underline{u}_1^\tau$. Since $d_{11}$ is nondecreasing, we always have $\lim_{z \to +\infty} \frac{s_{11}(z)}{zd_{11}(z)} = 0$, so that we can proceed as for the self-diffusion term to get the uniform integrability, using (44) again (or a modified version replacing $d_{11}^\varepsilon$ by $d_{11}$ using Fatou's lemma). For $s_{12}(\underline{u}_2^\tau)\underline{u}_1^\tau$, the same method applies but with a slight difference : using assumption **H1** we know that for $z \ge Z$ large enough $s_{12}(z) \le \delta z \sqrt{d_{22}(z)} + a_{12}(z)$, then write for $E \subset Q_T$

$$ \int_E s_{12}(\underline{u}_2^\tau) \underline{u}_1^\tau \le \delta \sup_{\tau, \varepsilon} \int_{Q_T} \underline{u}_1^\tau \underline{u}_2^\tau \sqrt{d_{22}(\underline{u}_2^\tau) + a_{12}(\underline{u}_2^\tau)} + s_{12}(Z) \sup_{\tau, \varepsilon} \int_E \underline{u}_1^\tau \\ \le \delta \sup_{\tau, \varepsilon} \|\underline{u}_1^\tau\|_{L^2(Q_T)} \left\|\underline{u}_2^\tau \sqrt{d_{22}(\underline{u}_2^\tau)}\right\|_{L^2(Q_T)} + \delta \sup_{\tau, \varepsilon} \|\underline{u}_2^\tau\|_{L^2(Q_T)} \left\|\underline{u}_1^\tau \sqrt{a_{12}(\underline{u}_2^\tau)}\right\|_{L^2(Q_T)} \\ + s_{12}(Z) \sup_{\tau, \varepsilon} \int_E \underline{u}_1^\tau, $$

and use again estimate (44) to conclude. We get then

$$ \int_{Q_T} u_1 \left\{ -\partial_t \theta - [d_{11}(u_1) + a_{22}(u_2)] \Delta \theta \right\} dx dt = \int_{Q_T} R_{12}(u_1, u_2) \theta dx dt + \int_{\Omega} u_1^0(x) \theta(0, x) dx, $$

that is the weak formulation on $Q_T$ of the equation

$$ \partial_t u_1 - \Delta [a_{11}(u_1) + u_1 a_{12}(u_2)] = R_{12}(u_1, u_2), $$

with Neumann boundary conditions and initial data $u_1(0, x) = u_1^0(x)$. A similar (symmetric) formulation holds for $u_2$.
We end up the proof of our Theorem with the passage to the limit in the duality estimates. For this purpose just notice that since $(u_1^\tau \sqrt{d_{11}^\xi(u_1^\tau) + a_{12}^\xi(u_2^\tau))}_{(\epsilon, \tau)}$ converges almost everywhere to $u_1 \sqrt{d_{11}(u_1) + a_{12}(u_2)}$, we get by the classical weak estimate

$$ \| u_1 \sqrt{d_{11}(u_1) + a_{12}(u_2)} \|_{L^2(Q_T)} \le D_T(u^0). $$

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

The funding by the ANR project Kibord is acknowledged in the preparation of the final version of this work.

## References

[1] H. Amann. Erratum: "Dynamic theory of quasilinear parabolic systems. III. Global existence" [Math. Z. **202** (1989), no. 2, 219-250; MR1013086 (90i:35125)]. Math. Z., 205(2):331, 1990.

[2] J. W. Barrett and J. F. Blowey. Finite element approximation of a nonlinear cross-diffusion population model. Numer. Math., 98(2):195-221, 2004.

[3] M. Bendahmane, T. Lepoutre, A. Marrocco, and B. Perthame. Conservative cross diffusions and pattern formation through relaxation. Journal de Mathématiques Pures et Appliquées, 92(6):651-667, 2009.

[4] D. Bothe, M. Pierre, and G. Rolland. Cross-diffusion limit for a reaction-diffusion system with fast reversible reaction. Communications in Partial Differential Equations, 37(11):1940-1966, 2012.

[5] N. Boudiba and M. Pierre. Global existence for coupled reaction-diffusion systems. J. Math. Anal. Appl., 250:1-12, 2000.

[6] J. Canizo, L. Desvillettes, and K. Fellner. Improved duality estimates and applications to reaction-diffusion equations. Accepted for publication in Comm. Partial Differential Equations, preprint of CMLA, ENS Cachan (# 2013-07).

[7] L. Chen and A. Jüngel. Analysis of a multidimensional parabolic population model with strong cross-diffusion. SIAM J. Math. Anal., 36(1):301-322 (electronic), 2004.

[8] L. Chen and A. Jüngel. Analysis of a parabolic cross-diffusion population model without self-diffusion. J. Differential Equations, 224(1):39-59, 2006.

[9] Y. S. Choi, R. Lui, and Y. Yamada. Existence of global solutions for the Shigesada-Kawasaki-Teramoto model with strongly coupled cross-diffusion. Discrete Contin. Dyn. Syst., 10(3):719-730, 2004.

[10] P. Degond, S. Génieys, and A. Jüngel. Symmetrization and entropy inequality for general diffusion equations. C. R. Acad. Sci. Paris Sér. I Math., 325(9):963-968, 1997.

[11] L. Desvillettes and F. Conforto. Rigorous passage to the limit in a system of reaction-diffusion equations towards a system including cross diffusions. Accepted for publication in Commun. Math. Sciences, preprint of CMLA, ENS Cachan (# 2009-34).

[12] L. Desvillettes, K. Fellner, M. Pierre, and J. Vovelle. Global existence for quadratic systems of reaction-diffusion. Adv. Nonlinear Stud., 7(3):491-511, 2007.

[13] L. Desvillettes, T. Lepoutre, and A. Moussa. Entropy, Duality and Cross Diffusion. ArXiv e-prints, Feb. 2013.

[14] P. Deuring. An initial-boundary-value problem for a certain density-dependent diffusion system. Mathematische Zeitschrift, 194:375-396, 1987.

[15] M. Dreher and A. Jüngel. Compact families of piecewise constant functions in $L^p(0, T; B)$. Nonlinear Anal., 75(6):3072-3077, 2012.

[16] G. Galiano, M. L. Garzón, and A. Jüngel. Semi-discretization in time and numerical convergence of solutions of a nonlinear cross-diffusion population model. Numer. Math., 93(4):655-673, 2003.
[17] D. Gilbarg and N. S. Trudinger. *Elliptic partial differential equations of second order*, volume 224 of *Grundlehren der Mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences]*. Springer-Verlag, Berlin, second edition, 1983.

[18] S. Hittmeir and A. Jüngel. Cross diffusion preventing blow-up in the two-dimensional Keller-Segel model. *SIAM J. Math. Anal.*, 43(2):997–1022, 2011.

[19] M. Iida, M. Mimura, and H. Ninomiya. Diffusion, cross-diffusion and competitive interaction. *J. Math. Biol.*, 53(4):617–641, 2006.

[20] A. Jüngel and I. V. Stelzer. Entropy structure of a cross-diffusion tumor-growth model. *Math. Models Methods Appl. Sci.*, 22(7):1250009, 26, 2012.

[21] J. Kačur. *Method of Rothe in evolution equations*, volume 80 of *Teubner-Texte zur Mathematik [Teubner Texts in Mathematics]*. BSB B. G. Teubner Verlagsgesellschaft, Leipzig, 1985. With German, French and Russian summaries.

[22] S. Kawashima and Y. Shizuta. On the normal form of the symmetric hyperbolic-parabolic systems associated with the conservation laws. *Tohoku Math. J. (2)*, 40(3):449–464, 1988.

[23] Y. Lou and S. Martinez. Evolution of cross-diffusion and self-diffusion. *Journal of Biological Dynamics*, 3(4):410–429, 2009. PMID: 22876941.

[24] Y. Lou, W.-M. Ni, and Y. Wu. On the global existence of a cross diffusion system. *Discrete Contin. Dynam. Systems*, 4(2):193–203, 1998.

[25] R. H. Martin, Jr. and M. Pierre. Nonlinear reaction-diffusion systems. In *Nonlinear equations in the applied sciences*, volume 185 of *Math. Sci. Engrg.*, pages 363–398. Academic Press, Boston, MA, 1992.

[26] J. D. Murray. *Mathematical biology. II*, volume 18 of *Interdisciplinary Applied Mathematics*. Springer-Verlag, New York, third edition, 2003. Spatial models and biomedical applications.

[27] M. Pierre and D. Schmitt. Blowup in reaction-diffusion systems with dissipation of mass. *SIAM J. Math. Anal.*, 28(2):259–269, 1997.

[28] N. Shigesada, K. Kawasaki, and E. Teramoto. Spatial segregation of interacting species. *J. Theoret. Biol.*, 79(1):83–99, 1979.

[29] J. Simon. Compact sets in the space $L^p(0, T; B)$. *Ann. Mat. Pura Appl. (4)*, 146:65–96, 1987.

[30] A. Yagi. Global solution to some quasilinear parabolic system in population dynamics. *Nonlinear Anal.*, 21(8):603–630, 1993.