%~Mouliné par MaN_auto v.0.27.3 2023-08-18 09:12:21
\documentclass[AHL,Unicode,longabstracts,published]{cedram}


%
\makeatletter
\def\editors#1{%
\def\editor@name{#1}
\if@francais
\def\editor@string{Recommandé par les éditeurs \editor@name.}
\else
\def\editor@string{Recommended by Editors \editor@name.}
\fi}   
\makeatother

%
%\AtBeginDocument{
%%\newtheorem{myname}[cdrthm]{My beautiful theorem} }
%
\newcommand{\bbR}{\mathbb{R}}
\newcommand{\bbN}{\mathbb{N}}
\newcommand{\clE}{\mathcal{E}}
\newcommand{\divv}{\mathrm{div}}
\newcommand{\bfn}{\mathbf{n}}
\newcommand{\rmd}{\mathrm{d}}
\newcommand{\clP}{\mathcal{P}}
\newcommand{\bfM}{\mathbf{M}}
\newcommand{\bSym}{\mathbf{Sym}}
\newcommand{\bSPD}{\mathbf{SPD}}
\newcommand{\clV}{\mathcal{V}}
\newcommand{\trr}{\mathrm{tr}}



\newcounter{stepcount}
\newenvironment{step}{\refstepcounter{stepcount}\begin{proof}[\textbf{Step~\thestepcount}]}{\let\qed\relax\end{proof}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\graphicspath{{./figures/}}

\newcommand*{\mk}{\mkern -1mu}
\newcommand*{\Mk}{\mkern -2mu}
\newcommand*{\mK}{\mkern 1mu}
\newcommand*{\MK}{\mkern 2mu}

%\hypersetup{urlcolor=purple, linkcolor=blue, citecolor=red}

\newcommand*{\romanenumi}{\renewcommand*{\theenumi}{\roman{enumi}}}
\newcommand*{\Romanenumi}{\renewcommand*{\theenumi}{\Roman{enumi}}}
\newcommand*{\alphenumi}{\renewcommand*{\theenumi}{\alph{enumi}}}
\newcommand*{\Alphenumi}{\renewcommand*{\theenumi}{\Alph{enumi}}}


\newcommand*{\Penumi}{\renewcommand*{\theenumi}{P\arabic{enumi}}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\title[Bounded weak solutions to a class of degenerate cross-diffusion systems]{Bounded weak solutions to a class of degenerate cross-diffusion systems}
\alttitle{Solutions faibles bornées pour une classe de systèmes à diffusion croisée dégénérés}

\subjclass{35K65, 35K51, 37L45, 35B65}

\keywords{Degenerate parabolic system, cross-diffusion, boundedness, Liapunov functionals, global existence}


\author[\initial{Ph.} \lastname{Laurençot}]{\firstname{Philippe} \lastname{Laurençot}}
\address{Institut de Mathématiques de Toulouse,\\ UMR~5219, Université de Toulouse, CNRS \\ F--31062 Toulouse Cedex 9, France}
\email{laurenco@math.univ-toulouse.fr}
\thanks{Ph. L. gratefully acknowledges the hospitality and support of the Fakult\"at f\"ur Mathematik, Universit\"at Regensburg, where part of this work was done.}

\author[\initial{B.-V.} \lastname{Matioc}]{\firstname{Bogdan-Vasile} \lastname{Matioc}}
\address{Fakult\"at f\"ur Mathematik, \\
Universit\"at Regensburg \\
D--93040 Regensburg, Deutschland}
\email{bogdan.matioc@ur.de}

\begin{abstract}
Bounded weak solutions are constructed for a degenerate parabolic system with a full diffusion matrix, which is a generalized version of the thin film Muskat system. Boundedness is achieved with the help of a sequence $(\clE_n)_{n\,\ge\,2}$ of Liapunov functionals such that~$\clE_n$ is equivalent to the $L_n$-norm for each $n\ge 2$ and $\clE_n^{1/n}$ controls the $L_\infty$-norm in the limit~${n\to\infty}$. Weak solutions are built by a compactness approach, special care being needed in the construction of the approximation in order to preserve the availability of the above-mentioned Liapunov functionals.
\end{abstract}

\begin{altabstract}
Des solutions faibles bornées sont construites pour un système parabolique dégénéré avec une matrice de diffusion pleine, qui est une version généralisée d'une approximation de type \og~film mince\fg{} du système de Muskat. Le caractère borné des solutions est obtenu à l'aide d'une suite $(\clE_n)_{n\,\ge\,2}$ de fonctionnelles de Liapunov avec les propriétés suivantes~: $\clE_n$ est équivalente à la norme $L_n$ pour chaque $n\ge 2$ et~$\clE_n^{1/n}$ contrôle la norme $L_\infty$ dans la limite~${n\to\infty}$. Les solutions faibles sont construites par une méthode de compacité, la construction des approximations requérant une attention particulière afin d'être compatibles avec les fonctionnelles de Liapunov mentionnées ci-dessus.
\end{altabstract}

\datereceived{2022-01-17}
\daterevised{2023-06-06}
\dateaccepted{2023-06-09}

\editors{F. Hubert and N. Séguin}
\begin{DefTralics}
\newcommand{\clE}{\mathcal{E}}

\end{DefTralics}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%	
\dateposted{2023-10-02}
\begin{document}
\maketitle



\section{Introduction}\label{sec1}



Let $\Omega$ be a bounded domain of $\bbR^N$, $N\ge 1$, with smooth boundary $\partial\Omega$ and let $R$ and $\mu$ be two positive real numbers. In a recent paper~\cite{LM22}, we noticed that there is an infinite family $(\clE_n)_{n\,\ge\,1}$ of Liapunov functionals associated with the thin film Muskat system
\begin{align*}
\partial_t f & = \divv\left(f \nabla\left[(1+R)f + R g \right] \right) \;\text{ in }\; (0,\infty)\times \Omega\,, \\
\partial_t g & = \mu R\, \divv\left(g \nabla\left[f + g \right] \right) \;\text{ in }\; (0,\infty)\times \Omega\,,
\end{align*}
supplemented with homogeneous Neumann boundary conditions and initial conditions, with the following properties: for all $n\ge 2$, there are $0<c_n < C_n$ such that
\begin{equation*}
c_n \|f+g\|_n^n \le \clE_n(f,g) \le C_n \|f+g\|_n ^n\,, \qquad (f,g)\in L_{n,+}\left(\Omega,\bbR^2\right)\,,
\end{equation*}
and there are $0<c_\infty < C_\infty$ such that
\begin{equation*}
c_\infty \|f+g\|_\infty \le \liminf_{n\,\to\,\infty} \clE_n(f,g)^{1/n} \le \limsup_{n\,\to\,\infty} \clE_n(f,g)^{1/n} \le C_\infty \|f+g\|_\infty
\end{equation*}
for $(f,g)\in L_{\infty,+}(\Omega,\bbR^2)$, where $L_{p,+}(\Omega,\bbR^m)$ denotes the positive cone of $L_p(\Omega,\bbR^m)$ for $m\ge 1$ and~${p\in [1,\infty]}$.

On the one hand, the thin film Muskat system being of cross-diffusion type (i.e., featuring a diffusion matrix with no zero entry), the availability of such a family of Liapunov functionals is rather seldom within this class of systems and paves the way towards the construction of bounded weak solutions, a result that we were only able to show in one space dimension $N=1$ in~\cite{LM22}. On the other hand, it is tempting to figure out whether this property is peculiar to the thin film Muskat system or extends to the generalization thereof
\begin{subequations}\label{PB}
\begin{align}
\partial_t f & = \divv\left(f \nabla\left[af + b g \right] \right) \;\text{ in }\; (0,\infty)\times \Omega\,, \label{PB1a}
\\
\partial_t g & = \divv\left(g \nabla\left[cf + d g \right] \right) \;\text{ in }\; (0,\infty)\times \Omega\,, \label{PB1b}
\end{align}
with $(a,\, b,\, c,\, d) \in (0,\infty)^4$, supplemented with homogeneous Neumann boundary conditions
\begin{equation}\label{PB1c}
\nabla f\cdot \bfn = \nabla g\cdot \bfn = 0 \;\text{ on }\; (0,\infty)\times \partial\Omega\,, 
\end{equation}
and non-negative initial conditions
\begin{equation}\label{PB1d}
(f,g)(0) = \left(f^{in},g^{in}\right) \;\text{ in }\; \Omega\,, 
\end{equation}
\end{subequations}
which is proposed in~\cite[Section~4]{BGHP85} to describe the dispersal of two interacting population species and is also a particular case of a model of interacting particles derived in~\cite{GS14}. Obviously, the thin film Muskat system is a particular case of~\eqref{PB1a}-\eqref{PB1b}, corresponding to the choice
\begin{equation*}
(a,\, b,\, c,\, d)=(1+R, \, R,\, \mu R,\, \mu R)\,.
\end{equation*}

It is worth mentioning at this point that the existence of global weak solutions to several cross-diffusion systems relies on the availability of a Liapunov functional or an entropy and we refer to~\cite[Chapter~4]{J2016} and the references therein for results in that direction. In the most favourable cases, an a priori $L_\infty$-bound can even be retrieved from the structure of the entropy functional, see~\cite{BFPS10} for instance. In contrast, the cornerstone of our approach is the construction of countably infinitely many Liapunov functionals, leading to $L_\infty$-bounds after performing a suitable limiting process. Our contribution is somewhat closer in spirit to~\cite{JM06}, where an algorithmic method for the construction of Liapunov functionals is developed. Let us also mention~\cite{AM23}, where a system of two coupled degenerate parabolic equations (without cross-diffusion) is studied which also features an infinite family of Liapunov functionals. This family of functionals provides $L_n$-estimates for all~${n\geq 1}$, but no~$L_\infty$-bound as in~\cite{LM22} and herein.

Coming back to~\eqref{PB}, the main result of this paper is to show that, for any quadruple $(a,\, b,\, c,\, d)$ satisfying
\begin{equation}\label{condabcd}
(a,\, b,\, c,\, d) \in(0,\infty)^4 \qquad\text{and}\qquad ad>bc\,,
\end{equation}
we can associate a countably infinite family of Liapunov functionals with~\eqref{PB} and prove the global existence of bounded non-negative weak solutions to~\eqref{PB}, whatever the dimension $N\ge 1$. More precisely, given a quadruple $(a,\, b,\, c,\, d)$ satisfying~\eqref{condabcd}, we define a sequence $(\Phi_n)_{n\,\ge\,1}$ of functions as follows. Setting $L(r):=r\ln{r}-r+1\ge 0$ for~$r\ge 0$, we first define the function $\Phi_1$ by the relation
\begin{equation}\label{p4d}
\Phi_1(X) := L(X_1) + \frac{b^2}{ad}L(X_2)\,, \qquad X=(X_1,X_2)\in [0,\infty)^2\,. 
\end{equation}
Next, for each integer $n\geq 2$, let $\Phi_n$ be the homogeneous polynomial of degree $n$ defined by
\begin{equation}\label{p4b}
\Phi_n(X) := \sum_{j=0}^n a_{j,n} X_1^j X_2^{n-j}\,, \qquad X=(X_1,X_2)\in \bbR^2\,, 
\end{equation}
with $a_{0,n}:=1$ and
\begin{equation}\label{p4c}
a_{j,n}:=\binom{n}{j} \prod_{k=0}^{j-1}\frac{ak+c(n-k-1)}{bk+d(n-k-1)}> 0\,,\qquad 1\leq j\leq n\,.
\end{equation}
We then define, for $n\ge 1$, the functional
\begin{equation}\label{p4}
\clE_n(u):=\int_{\Omega}\Phi_n(u(x))\, dx,\qquad u=(f,g)\in L_{\max\{2,n\},+}\left(\Omega,\bbR^2\right)\,.
\end{equation}
We finally observe that~\eqref{condabcd} guarantees that
\begin{equation}\label{constants}
\Theta_1:=\frac{b(ad+bc)}{2ad}>0\qquad\text{and}\qquad \Theta_2:=\frac{b^2(ad-bc)(3ad+bc)}{4a^2d^2}>0\,.
\end{equation}

%\pagebreak

With this notation, the main result of this paper is the following:


\begin{theo}\label{ThBWS}
Assume~\eqref{condabcd} and let $ u^{in} := (f^{in},g^{in})\in L_{\infty,+}(\Omega,\bbR^2)$ be given. Then, there is a bounded weak solution $u=(f,g)$ to~\eqref{PB} such that:
\begin{enumerate}\romanenumi
\item for each $T>0$,
\begin{equation}\label{p1}
\begin{split}
(f,g)\in L_{\infty,+}\left((0,T)\times \Omega,\bbR^2\right) &\cap L_2\left((0,T),H^1\left(\Omega,\bbR^2\right)\right)\\
&\cap W_2^1\left((0,T),H^1 \left(\Omega,\bbR^2\right)'\right)\,;
\end{split} 
\end{equation}
\item for all $\varphi\in H^1(\Omega)$ and $t\ge 0$,
\begin{subequations}\label{p2}
\begin{multline}\label{p2a}
\int_\Omega \left(f(t,x)-f^{in}(x)\right) \varphi(x)\ \rmd x \\
+\int_0^t \int_\Omega f(s,x) \nabla[a f + b g](s,x) \cdot \nabla\varphi(x)\, \rmd x\rmd s =0 
\end{multline}
and
\begin{multline}\label{p2b}
\int_\Omega \left(g(t,x)-g^{in}(x)\right) \varphi(x)\ \rmd x \\
+ \int_0^t \int_\Omega g(s,x) \nabla[cf + dg](s,x) \cdot \nabla\varphi(x)\, \rmd x\rmd s =0\,; 
\end{multline}
\end{subequations}
\item for all $t\ge 0$,
\begin{equation}\label{p3}
\begin{aligned}
&\clE_1(u(t))+ \frac{1}{a}\int_0^t \int_\Omega\big[|\nabla (af+\Theta_1 g)|^2+\Theta_2 |\nabla g|^2\big](s,x)\ \rmd x\rmd s\le \clE_1(u^{in})\,,
\end{aligned}
\end{equation}
where the positive constants $\Theta_1$ and $\Theta_2$ are defined in~\eqref{constants};
\item for all $n\ge 2$ and all $t\ge 0$,
\begin{equation}
\clE_n(u(t)) \le \clE_n(u^{in})\,; \label{p4a}
\end{equation}
\item for $t\ge 0$,
\begin{equation}
\|f(t)+g(t)\|_\infty\le \frac{d}{b}\frac{\max\{a,\,b\}}{\min\{c,\, d\}} \left\|f^{in}+g^{in}\right\|_\infty\,. \label{p5}
\end{equation}
\end{enumerate}
\end{theo}


Let us first mention that Theorem~\ref{ThBWS} improves~\cite{LM22} in two directions: on the one hand, it shows that the structural properties~\eqref{p3}, \eqref{p4a}, and~\eqref{p5}, uncovered there for the thin film Muskat system, are also available for the whole class~\eqref{PB}. On the other hand, it provides the existence of non-negative bounded weak solutions to~\eqref{PB} in all space dimensions, a result which was only established in one space dimension in~\cite{LM22}. Theorem~\ref{ThBWS} may also be viewed as a partial extension of~\cite{GS14}, where global weak solutions to~\eqref{PB} are constructed in space dimensions $N\in\{1,\, 2,\, 3\}$ when the coefficients~$a,\,b,\, c,\, d$ are non-negative bounded functions which satisfy a more restrictive condition than~\eqref{condabcd}, namely the inequality~${4ad-(b+c)^2>\lambda}$ for some positive constant $\lambda$. Whether the analysis performed below could be adapted to non-constant coefficients is yet unclear. Let us also mention that global weak solutions to the thin film Muskat system are also constructed in~\cite{AIJM2018, BGB2019, ELM2011, LM2013, LM2017,ACCL2019}, but their boundedness is an open question, to which an affirmative answer is only provided in~\cite{BGB2019}. The latter however requires some smallness condition on the initial data, in contrast to Theorem~\ref{ThBWS}. Finally, the local well-posedness of the thin film Muskat system in the classical sense is investigated in~\cite{EMM2012}.

%\bigskip

We next outline the main steps of the proof of Theorem~\ref{ThBWS}. As in~\cite{LM22}, the starting point is to notice that, introducing the mobility matrix
\begin{equation}\label{x1}
M(X) = \left(m_{jk}(X)\right)_{1\,\le\,j,\,k\,\le \,2} :=
\begin{pmatrix}
a X_1& b X_1 \\
c X_2 & d X_2
\end{pmatrix}\,, \qquad X=(X_1,X_2)\in \bbR^2\,, 
\end{equation}
and $u:=(f,g)$, an alternative formulation of the system~\eqref{PB1a}-\eqref{PB1b} is
\begin{equation}\label{x2}
\partial_t u = \sum_{i=1}^N\partial_i (M(u)\partial_i u) \;\text{ in }\; (0,\infty)\times \Omega\,. 
\end{equation}
Then, given $\Phi\in C^2(\bbR^2,\bbR)$, it readily follows from~\eqref{x2}, the homogeneous Neumann boundary conditions~\eqref{PB1c}, and the symmetry of the Hessian matrix $D^2\Phi$ that
\begin{equation}\label{x3}
\frac{\rmd }{\rmd t} \int_\Omega \Phi(u)\ \rmd x + \sum_{i=1}^N\int_\Omega \left\langle D^2\Phi(u) M(u) \partial_i u, \partial_i u \right\rangle\ \rmd x = 0\,, 
\end{equation}
where $\langle \cdot,\cdot \rangle$ stands for the scalar product on $\bbR^2$. As a straightforward consequence of~\eqref{x3} we note that $\int_\Omega \Phi(u)\ \rmd x$ is a Liapunov functional for~\eqref{x2} when the matrix~${D^2\Phi(u) M(u)}$ is positive semidefinite. We shall then show in Appendix~\ref{secA} that, for all $n\ge 2$, it is possible to construct an homogeneous polynomial $\Phi_n\in \bbR[X_1,X_2]$ of degree $n$ which is convex on $[0,\infty)^2$ and such that the matrix $D^2\Phi_n(X) M(X)$ is positive semidefinite for all $X\in [0,\infty)^2$. A closed form formula is actually available for the polynomial $\Phi_n$, see~\eqref{p4b} and~\eqref{p4c}.

We next construct weak solutions to~\eqref{x2} by a compactness method. It is here of utmost importance to construct approximations which do not alter the inequalities~\eqref{x3} for $\Phi=\Phi_n$ and~${n\ge 1}$. As a first step, it is well-known that implicit time discrete schemes are well-suited in that direction. Thus, given $\tau>0$, we shall first prove the existence of a sequence $(u^\tau_l)_{l\,\ge\,0}$ which satisfies~${u^\tau_0=u^{in}:=(f^{in},g^{in})}$ and, for $l\ge 0$,
\begin{equation}\label{x4}
u^\tau_{l+1} - \tau \sum_{i=1}^N\partial_i \Big(M(u^\tau_{l+1})\partial_i u^\tau_{l+1}\Big) = u^\tau_l \;\;\text{ in }\;\; \Omega\,, 
\end{equation}
supplemented with homogeneous Neumann boundary conditions. Furthermore, the sequence $(u^\tau_l)_{l\,\ge\,0}$ has the property that, for $n\ge 1$ and $l\ge 0$,
\begin{equation}\label{x5}
\clE_n\left(u_{l+1}^\tau\right) + \tau \sum_{i=1}^N\int_\Omega \left\langle D^2\Phi_n\left(u_{l+1}^\tau\right) M\left(u_{l+1}^\tau\right) \partial_i u_{l+1}^\tau, \partial_i u_{l+1}^\tau \right\rangle\ \rmd x \le \clE_n(u_{l}^\tau)\,, 
\end{equation}
so that the structural property~\eqref{x3} is indeed preserved by the time discrete scheme. The existence of a solution to~\eqref{x4} is achieved by a compactness method relying on an approximation of the matrix $M(\cdot)$ by bounded ones. This step is actually the more delicate one, as we have to construct matrices approximating $M(\cdot)$ which do not alter~\eqref{x5}. To this end, a two-parameter approximation procedure is required and it is detailed in Section~\ref{sec2.2}. The existence of a weak solution to~\eqref{x4} satisfying~\eqref{x5} is shown in Section~\ref{sec2.4}, building upon preliminary and intermediate results established in Section~\ref{sec2.1} and Section~\ref{sec2.3}.


\begin{rema}\label{rem.W}
A common feature of system~\eqref{PB} is that it has, at least formally, a gradient flow structure for the functional $\clE_2$ with respect to the $2$-Wasserstein distance in the space $\clP_2(\Omega,\bbR^2)$ of probability measures with finite second moments, as pointed out in~\cite{LM2013,ACCL2019} for the thin film Muskat system. In particular, there is a natural variational structure associated with~\eqref{PB} which is suitable to construct weak solutions. However, the connection between this variational structure and the whole family $(\clE_n)_{n\,\ge\,2}$ of Liapunov functionals is yet unclear.
\end{rema}


%\bigskip

\begin{nota}
For $p\in [1,\infty]$, we denote the $L_p$-norm in $L_p(\Omega)$ by $\|\cdot\|_p$ and set
\begin{equation*}
L_p\left(\Omega,\bbR^2\right) := L_p(\Omega)\times L_p(\Omega)\,, \quad H^1\left(\Omega,\bbR^2\right):= H^1(\Omega)\times H^1(\Omega)\,.
\end{equation*}
The positive cone of a Banach lattice $E$ is denoted by $E_+$. The space of $2\times 2$ real-valued matrices is denoted by ${\bfM_2(\bbR)}$, while ${\bSym_2(\bbR)}$ is the subset of ${\bfM_2(\bbR)}$ consisting of symmetric matrices and~${\bSPD_2(\bbR)}$ is the set of symmetric and positive definite matrices in $\bfM_2(\bbR)$. The positive part of a real number $r\in\bbR$ is given by~${r_+:=\max\{r,0\}}$ and, for $X=(X_1,X_2)\in\bbR^2$, we define the positive part of~$X$ componentwise; that is, ${X_+:=(X_{1,+}, X_{2,+})}$. Finally, $\langle \cdot,\cdot \rangle$ is the scalar product on~$\bbR^2$.
\end{nota}


\section{A time discrete scheme}\label{sec2}



In order to construct bounded non-negative global weak solutions to the evolution problem~\eqref{PB}, we employ a compactness approach, paying special attention to preserve as much as possible the structural properties~\eqref{p3}, \eqref{p4a}, and~\eqref{p5} in the design of the approximation. It turns out that implicit time discrete schemes are well-suited for that purpose and we thus establish in this section the existence of solutions to the implicit time discrete scheme associated with~\eqref{PB}, see~\eqref{ex1a}-\eqref{ex1b}.


\begin{prop}\label{P:1}
Given $\tau>0$ and $U=(F,G)\in L_{\infty,+}(\Omega,\bbR^2)$, there is a solution
\begin{equation*}
u=(f,g) \in H^1\left(\Omega,\bbR^2\right)\cap L_{\infty,+}\left(\Omega,\bbR^2\right)
\end{equation*}
to
\begin{subequations}\label{ex1}
\begin{align}
\int_\Omega \big(f \varphi + \tau f \nabla\left[a f + b g \right] \cdot\nabla\varphi \big)\ \rmd x & = \int_\Omega F \varphi\ \rmd x\,, \quad \varphi\in H^1(\Omega)\,, \label{ex1a}
\\
\int_\Omega \big(g \psi + \tau g \nabla\left[c f +d g \right]\cdot \nabla\psi \big)\ \rmd x & = \int_\Omega G \psi\ \rmd x\,, \quad \psi\in H^1(\Omega)\,, \label{ex1b}
\end{align}
\end{subequations}
which also satisfies
\begin{equation}\label{ex2}
\clE_n(u) \le \clE_n(U)\qquad\text{for $n\ge 2$}
\end{equation}
and
\begin{equation}\label{ex2b}
\clE_1(u) + \frac{\tau}{a} \int_\Omega \big[|\nabla (af+\Theta_1 g)|^2 + \Theta_2|\nabla g|^2 \big]\ \rmd x \le \clE_1(U) \,, 
\end{equation}
recalling that, see~\eqref{constants},
\begin{equation*}
\Theta_1 = \frac{b(ad+bc)}{2ad}>0 \;\;\text{ and }\;\; \Theta_2 = \frac{b^2(ad-bc)(3ad+bc)}{4a^2 d^2}>0\,.
\end{equation*}
\end{prop}


As already mentioned, several steps are involved in the proof of Proposition~\ref{P:1}. We begin with the existence of bounded weak solutions to an auxiliary elliptic system which shares the same structure with~\eqref{ex1}, but has bounded coefficients instead of linearly growing ones, see Section~\ref{sec2.1}. As a next step, we introduce in Section~\ref{sec2.2} the approximation to~\eqref{ex1} which is derived from~\eqref{ex1} by replacing the matrix $M(\cdot)$ defined in~\eqref{x1} by a suitable invertible and bounded matrix $M_\varepsilon^\rho(\cdot)$ with~$(\varepsilon,\rho)\in (0,1)\times (1,\infty)$. We emphasize here once more that the matrix $M_\varepsilon^\rho(\cdot)$ is designed in such a way that the inequalities~\eqref{ex2} and~\eqref{ex2b} are not significantly altered. Passing to the limit, first as $\rho\to\infty$, and then as $\varepsilon\to 0$, is performed in Section~\ref{sec2.3} and Section~\ref{sec2.4}, respectively, this last step completing the proof of Proposition~\ref{P:1}.

Throughout this section, $C$ and $(C_l)_{l\,\ge\,0}$ denote various positive constants depending only on~$N$,~$\Omega$, and $(a,\, b,\, c,\, d)$. Dependence upon additional parameters will be indicated explicitly.



\subsection{An auxiliary elliptic system}\label{sec2.1}



Let $A=(a_{jk})_{1\,\le\,j,\,k\,\le\,2}$ and~$B=(b_{jk})_{1\,\le\,j,\,k\,\le\,2}$ be chosen such that~${A\in \bSPD_2(\bbR)}$ and~${B\in BC(\bbR^2,\bfM_2(\bbR))}$, with $AB(X)\in \bSPD_2(\bbR)$ for all $X\in \bbR^2$. Moreover, we assume that there is $\delta_1>0$ such that
\begin{equation}\label{ap1}
\langle AB(X)\xi,\xi \rangle \ge \delta_1 |\xi|^2\,, \qquad (X,\xi)\in \bbR^2\times\bbR^2\,. 
\end{equation}
Since $A\in \bSPD_2(\bbR)$, there is also $\delta_2>0$ such that
\begin{equation}\label{ap2}
\langle A\xi,\xi \rangle \ge\delta_2 |\xi|^2\,, \qquad \xi\in\bbR^2\,. 
\end{equation}


\begin{lemm}\label{lem.ap1}
Given $\tau>0$ and $U=(U_1,U_2)\in L_2(\Omega,\bbR^2)$, there exists a solution~$u=(u_1,u_2)\in H^1(\Omega,\bbR^2)$ to the nonlinear equation
\begin{equation}\label{ap3}
\int_\Omega \left[\langle u, v \rangle + \tau \sum_{i=1}^N\left\langle B(u) \partial_i u, \partial_i v \right\rangle \right]\ \rmd x = \int_\Omega \langle U, v \rangle\ \rmd x\,, \qquad v\in H^1\left(\Omega,\bbR^2\right)\,. 
\end{equation}
Additionally:
\begin{enumerate}\romanenumi
\item\label{lemm2.2.i} If
\begin{equation}\label{ap10}
\begin{split}
b_{11}(X) \geq b_{12}(X) & = 0\,, \qquad X\in (-\infty,0)\times \bbR\,, \\
b_{22}(X) \geq b_{21}(X) & = 0\,, \qquad X\in \bbR\times (-\infty,0)\,,
\end{split} 
\end{equation}
and if $U(x)\in [0,\infty)^2$ for a.a. $x\in\Omega$, then $u(x)\in [0,\infty)^2$ for a.a. $x\in\Omega$.
\item\label{lemm2.2.ii} If there exists $\rho>0$ such that
\begin{equation}\label{ap10'}
\begin{split}
b_{11}(X) \geq b_{12}(X) & = 0\,, \qquad X\in (\rho,\infty)\times \bbR\,, \\
b_{22}(X) \geq b_{21}(X) & = 0\,, \qquad X\in \bbR\times (\rho,\infty)\,,
\end{split} 
\end{equation}
and if $\max\{U_1, U_2\} \leq \rho$ a.e. in $\Omega$, then $\max\{u_1, u_2\}\leq \rho$ a.e. in $\Omega$.
\end{enumerate}
\end{lemm}


\begin{proof}
The proof of Lemma~\ref{lem.ap1} is rather classical and it is actually similar to that of~\cite[Lemma~B.1]{LM22}. We nevertheless sketch it below for the sake of completeness.

\begin{step}
To set up a fixed point scheme, we consider $u\in L_2(\Omega,\bbR^2)$ and define a bilinear form $b_u$ on~${H^1(\Omega,\bbR^2)}$ by
\begin{equation*}
b_u(v,w) := \int_\Omega \left[\langle Av, w \rangle + \tau \sum_{i=1}^N\left\langle AB(u) \partial_i v, \partial_i w \right\rangle \right]\ \rmd x
\end{equation*}
for $ (v,w)\in H^1(\Omega,\bbR^2)\times H^1(\Omega,\bbR^2)$. Owing to~\eqref{ap1} and~\eqref{ap2},
\begin{equation}\label{ap4}
b_u(v,v) \ge \delta_0 \|v\|_{H^1}^2\,, \qquad v\in H^1\left(\Omega,\bbR^2\right)\,, 
\end{equation}
where $\delta_0:=\min\{\tau \delta_1, \delta_2\},$ while the boundedness of $B$ guarantees that
\begin{equation*}
|b_u(v,w)| \le b_* \|v\|_{H^1} \|w\|_{H^1}\,, \qquad (v,w)\in H^1\left(\Omega,\bbR^2\right)\times H^1\left(\Omega,\bbR^2\right) \,,
\end{equation*}
with
\begin{equation*}
b_* := 2\max_{1\,\le\,j,\,k\,\le\,2}\left\{\left|a_{jk}\right|\right\} \left(1 + 2\tau \max_{1\,\le\,j,\,k\,\le\,2}\left\{\left\|b_{jk}\right\|_\infty\right\} \right)\,.
\end{equation*}
We then infer from Lax--Milgram's theorem that there is a unique $\clV[u]\in H^1(\Omega,\bbR^2)$ such that
\begin{equation}\label{ap5}
b_u(\clV[u],w) = \int_\Omega \langle AU, w \rangle\ \rmd x\,, \qquad w\in H^1\left(\Omega,\bbR^2\right)\,.
\end{equation}
An immediate consequence of~\eqref{ap4}, \eqref{ap5} (with $w=\clV[u]$), and H\"older's inequality is the following estimate:
\begin{equation*}
\delta_0 \|\clV[u]\|_{H^1}^ 2 \le b_u\left(\clV[u],\clV[u]\right) \le \|AU\|_2 \|\clV[u]\|_2 \le \|AU\|_2 \|\clV[u]\|_{H^1}\,.
\end{equation*}
Hence
\begin{equation}\label{ap6}
\|\clV[u]\|_{H^1} \le \frac{\|AU\|_2}{\delta_0}\,. 
\end{equation}

We next argue as in the proof of~\cite[Lemma~B.1]{LM22} to show that the map $\clV$ is continuous and compact from $L_2(\Omega,\bbR^2)$ to itself, the proof relying on~\eqref{ap6}, the compactness of the embedding of~$H^1(\Omega,\bbR^2)$ in $L_2(\Omega,\bbR^2)$, and the continuity and boundedness of $B$.

Consider now $\theta\in [0,1]$ and a function $u\in L_2(\Omega,\bbR^2)$ satisfying $u = \theta\clV[u]$. Then we have $u\in H^1(\Omega,\bbR^2)$ and, in view of~\eqref{ap6},
\begin{equation*}
\|u\|_2 = \theta \|\clV[u]\|_2 \le \|\clV[u]\|_2 \le \|\clV[u]\|_{H^1} \le \frac{\|AU\|_2}{\delta_0}\,.
\end{equation*}
Thanks to the above bound and the continuity and compactness properties of the map $\clV$ in~$L_2(\Omega,\bbR^2)$, we are in a position to apply Leray--Schauder's fixed point theorem, see~\cite[Theorem~11.3]{GT2001} for instance, and conclude that the map $\clV$ has a fixed point $u\in L_2(\Omega,\bbR^2)$. Since $\clV$ ranges in $H^1(\Omega,\bbR^2)$, the function~$u$ actually belongs to $H^1(\Omega,\bbR^2)$ and satisfies
\begin{equation*}
b_u(u,w) = \int_\Omega \langle AU, w \rangle\ \rmd x\,, \qquad w\in H^1\left(\Omega,\bbR^2\right)\,.
\end{equation*}
Finally, given $v\in H^1(\Omega,\bbR^2)$, the function $w=A^{-1}v$ also belongs to $H^1(\Omega,\bbR^2)$ and we infer from the above identity and the symmetry of $A$ that
\begin{align*}
\int_\Omega \langle U, v \rangle\ \rmd x = \int_\Omega \langle AU, w \rangle\ \rmd x & = b_u(u,w) =b_u\left(u,A^{-1}v\right) \\
& = \int_\Omega \left[\langle u, v \rangle + \tau \sum_{i=1}^N\left\langle B(u) \partial_i u, \partial_i v \right\rangle \right]\ \rmd x \,.
\end{align*}
We have thus constructed a weak solution $u\in H^1(\Omega,\bbR^2)$ to~\eqref{ap3}.
\end{step}

%\noindent\textbf{Step~2.}
\begin{step}
We now turn to the sign-preserving property~\eqref{lemm2.2.i} and assume $U(x)\in [0,\infty)^2$ for a.a.~${x\in\Omega}$. Let $u\in H^1(\Omega,\bbR^2)$ be a weak solution to~\eqref{ap3} and set~$\varphi:=-u$. The function $(\varphi_{1,+},\varphi_{2,+})$ then belongs to $H^1(\Omega,\bbR^2)$ and it follows from~\eqref{ap3} that
\begin{multline}\label{ap11}
 \int_\Omega \left[\varphi_1 \varphi_{1,+} + \varphi_2 \varphi_{2,+} + \tau \sum_{i=1}^N\sum_{j,k=1}^2 b_{jk}(u) \partial_i \varphi_k \partial_i (\varphi_{j,+}) \right]\ \rmd x \\
 = - \int_\Omega \left(U_1 \varphi_{1,+} + U_2 \varphi_{2,+} \right)\ \rmd x \le 0\,.
\end{multline}
We now infer from~\eqref{ap10} that, for $1\le i \le N$,
\begin{align*}
b_{11}(u) \partial_i\varphi_1 \partial_i\varphi_{1,+} & = b_{11}(u) \mathbf{1}_{(-\infty,0)}(u_1) |\partial_i u_1|^2\geq 0\,, \\
b_{12}(u) \partial_i\varphi_2 \partial_i\varphi_{1,+} & = b_{12}(u) \mathbf{1}_{(-\infty,0)}(u_1) \partial_i u_1 \partial_i u_2 = 0\,, \\
b_{21}(u) \partial_i\varphi_1 \partial_i\varphi_{2,+} & = b_{21}(u) \mathbf{1}_{(-\infty,0)}(u_2) \partial_i u_1 \partial_i u_2 = 0\,, \\
b_{22}(u) \partial_i\varphi_2 \partial_i\varphi_{2,+} & = b_{22}(u) \mathbf{1}_{(-\infty,0)}(u_2) |\partial_i u_2|^2 \geq0\,,
\end{align*}
so that the second term on the left-hand side of~\eqref{ap11} is non-negative. Consequently,~\eqref{ap11} gives
\begin{equation*}
\int_\Omega \left[|\varphi_{1,+}|^2 + |\varphi_{2,+}|^2 \right]\ \rmd x \le 0\,,
\end{equation*}
which implies that $\varphi_{1,+}=\varphi_{2,+}=0$ a.e. in $\Omega$. Hence, $u(x)\in [0,\infty)^2$ for a.a. $x\in\Omega$ as claimed.
\end{step}
%\medskip

%\noindent\textbf{Step~3.}
\begin{step}
It remains to prove~\eqref{lemm2.2.ii}. We thus assume that $\max\{U_1,U_2\} \le \rho$ a.e. in $\Omega$ and consider a weak solution $u\in H^1(\Omega,\bbR^2)$ to~\eqref{ap3}. As $v=((u_1-\rho)_+, (u_2-\rho)_+)$ belongs to $H^1(\Omega,\bbR^2)$, we deduce from~\eqref{ap3} that
\begin{align*}
\int_\Omega \left[\sum_{j=1}^2 (u_j-U_j) (u_j-\rho)_+ + \tau \sum_{i=1}^N\sum_{j,k=1}^2 b_{jk}(u) \partial_i u_k \partial_i (u_j-\rho)_+ \right]\ \rmd x=0\,.
\end{align*}
On the one hand,
\begin{equation*}
u_j-U_j \ge u_j - \rho \;\text{ a.e. in}\; \Omega\,, \qquad j=1,2\,,
\end{equation*}
so that
\begin{equation*}
\left(u_j-U_j\right) \left(u_j-\rho\right)_+ \ge \left(u_j-\rho\right) \left(u_j-\rho\right)_+ = \left(u_j-\rho\right)_+^2 \;\text{ a.e. in }\; \Omega\,, \qquad j=1,2\,.
\end{equation*}
On the other hand, we infer from~\eqref{ap10'} that, for $1\le i \le N$,
\begin{align*}
b_{11}(u) \partial_i u_1 \partial_i(u_1-\rho)_+ & = b_{11}(u) \mathbf{1}_{(\rho,\infty)}(u_1) |\partial_i u_1|^2\geq 0\,, \\
b_{12}(u) \partial_i u_2 \partial_i(u_1-\rho)_+ & = b_{12}(u) \mathbf{1}_{(\rho,\infty)}(u_1) \partial_i u_1 \partial_i u_2 = 0\,, \\
b_{21}(u) \partial_i u_1 \partial_i(u_2-\rho)_+ & = b_{21}(u) \mathbf{1}_{(\rho,\infty)}(u_2) \partial_i u_1 \partial_i u_2 = 0\,, \\
b_{22}(u) \partial_i u_2 \partial_i(u_2-\rho)_+ & = b_{22}(u) \mathbf{1}_{(\rho,\infty)}(u_2) |\partial_i u_2|^2 \geq0\,.
\end{align*}
Therefore,
\begin{equation*}
\sum_{j=1}^2 \int_\Omega (u_j-\rho)_+^2\, \rmd x\leq 0\,,
\end{equation*}
from which we deduce that $\max\{u_1,u_2\}\le\rho$ a.e. in $\Omega$.\qedhere
\end{step}\let\qed\relax
\end{proof}



\subsection{A regularised system}\label{sec2.2}



We now introduce the two-parameter approximation of~\eqref{ex1} on which the subsequent analysis relies. Specifically, given $\rho>1$, we define
\begin{equation*}
\alpha_\rho(z):=
\left\{
\begin{array}{clll}
0, &z\leq 0,\\
z, &0\leq z\leq \rho-1,\\
(\rho-1)(\rho-z), &\rho-1\leq z\leq \rho,\\
0, &z\geq\rho,
\end{array}
\right.
\end{equation*}
and observe that $\alpha_\rho\in BC(\bbR)$ with
\begin{equation*}
0\le \alpha_\rho(z)\le \min\{\rho, z_+\}\,, \qquad z\in\bbR\,.
\end{equation*}
Next, for $\varepsilon\in (0,1)$ and $X\in\bbR^2$, we set
\begin{equation*}
M_{\varepsilon}^\rho(X) = \left(m_{\varepsilon, jk}^\rho(X)\right)_{1\,\le\,j,\,k\,\le\,2}:= \varepsilon I_2 + \lambda_\varepsilon(X_{+}) M^\rho(X),
\end{equation*}
where
\begin{equation}\label{momar}
M^\rho(X) = \left(m_{jk}^\rho(X)\right)_{1\,\le\, j,\,k\,\le\,2} :=
\begin{pmatrix}
a\alpha_\rho(X_1) & b\alpha_\rho(X_1) \\
c\alpha_\rho(X_2) & d\alpha_\rho(X_2)
\end{pmatrix}\,, \qquad X\in\bbR^2\,,
\end{equation}
and
\begin{equation*}
\lambda_\varepsilon(X) := \frac{2}{1+\exp{[\varepsilon (X_1+X_2)}]}\,, \qquad X\in\bbR^2\,.
\end{equation*}
Note that $(M^\rho)_{\rho\,>\,1}$ converges to $M$, defined in~\eqref{x1}, locally uniformly in $[0,\infty)^2$ as~$\rho\to\infty$, while~$(\lambda_\varepsilon)_{\varepsilon\,\in\,(0,1)}$ converges to $1$ locally uniformly in $\bbR^2$ as $\varepsilon\to 0$. In fact, for $R>0$,
\begin{equation}\label{UL}
|\lambda_\varepsilon(X) - 1| \le 2 R \varepsilon\,, \qquad X\in [-R,R]^2\,. 
\end{equation}

The outcome of this section is that, given $\tau>0$, $\varepsilon\in (0,1)$, $\varrho>1$, and a function~$U\in L_{\infty,+}(\Omega,\bbR^2)$, there is a weak solution $u_\varepsilon^\rho\in H^1(\Omega,\bbR^2)\cap L_{\infty,+}(\Omega,\bbR^2)$ to
\begin{equation*}
u_\varepsilon^\rho - \tau \sum_{i=1}^N\partial_i \big(M_\varepsilon^\rho(u_\varepsilon^\rho) \partial_i u_\varepsilon^\rho \big) = U \;\text{ in }\; \Omega\,,
\end{equation*}
which satisfies an appropriate weak version of~\eqref{ex2}, as stated below. The next lemma is actually the building block of the proof of Proposition~\ref{P:1}.


\begin{lemm}\label{lem.ex2}
Given $\tau>0$, $U=(F,G)\in L_{\infty,+}(\Omega,\bbR^2)$, $\varepsilon\in (0,1)$, and a real number~${\rho\geq \max\{1,\|F\|_\infty,\|G\|_\infty\}}$, there exists a weak solution $u_{\varepsilon}^\rho=(u_{\varepsilon,1}^\rho, u_{\varepsilon,2}^\rho)$ in~${H^1(\Omega,\bbR^2)\cap L_{\infty,+}(\Omega,\bbR^2)}$ to
\begin{multline}\label{ex12a}
\int_\Omega \left[\left\langle u_{\varepsilon}^\rho, v \right\rangle + \tau \sum_{i=1}^N\big\langle M_{\varepsilon }^\rho(u_{\varepsilon}^\rho) \partial_i u_{\varepsilon}^\rho, \partial_i v \big\rangle \right]\ \rmd x\\ = \int_\Omega \langle U, v \rangle\ \rmd x\,, \quad v\in H^1\left(\Omega,\bbR^2\right)\,, 
\end{multline}
which additionally satisfies
\begin{align}
\max\left\{\|u_{\varepsilon,1}^\rho\|_\infty,\,\|u_{\varepsilon,2}^\rho\|_\infty\right\}&\leq \rho\label{ex12a'}
,\\[1ex]
\left\|u_{\varepsilon}^\rho\right\|_2&\leq C_0\|U\|_2,\label{ex12d}
\\[1ex]
\left\|\nabla u_{\varepsilon}^\rho\right\|_2&\leq \frac{C_1}{\sqrt{\tau\varepsilon}} \|U\|_2\,.\label{ex12e}
\end{align}

Moreover, given $n\ge 2$, there exists a constant $C(n)$ such that
\begin{equation}\label{ex12c}
\clE_n(u_\varepsilon^\rho)\le \tau C(n) \frac{\rho^{n-1}}{e^{\varepsilon\rho}}\|\nabla u_\varepsilon^\rho\|_2^2\ + \clE_n(U)\,. 
\end{equation}
\end{lemm}


\begin{proof}
Let $\varepsilon\in (0,1)$ and $\rho\geq \max\{1,\|F\|_\infty,\|G\|_\infty\}$. To deduce the existence result stated in Lemma~\ref{lem.ex2} from the already established Lemma~\ref{lem.ap1}, we first recast~\eqref{ex12a} in the form~\eqref{ap3}. First, owing to the definition of the function $\alpha_\rho$, the matrix $M_\varepsilon^\rho$ lies in $BC(\bbR^2,\bfM_2(\bbR))$ and satisfies
\begin{subequations}\label{ex11}
\begin{equation}\label{ex11a}
0 \le m_{\varepsilon,jk}^\rho(X) \le \varepsilon + 2 \rho\max\{a,\, b,\, c,\, d\}\,, \qquad 1\le j,k\le 2\,, \ X\in\bbR^2\,, 
\end{equation}
as well as
\begin{equation}\label{ex11b}
\begin{split}
m_{\varepsilon,11}^\rho(X) \geq m_{\varepsilon,12}^\rho(X) & = 0\,, \qquad X\in (-\infty,0)\times \bbR\,, \\
m_{\varepsilon,22}^\rho(X) \geq m_{\varepsilon,21}^\rho(X) & = 0\,, \qquad X\in \bbR\times (-\infty,0)\,.
\end{split} 
\end{equation}
and
\begin{equation}\label{ex11b'}
\begin{split}
m_{\varepsilon,11}^\rho(X) \geq m_{\varepsilon,12}^\rho(X) & = 0\,, \qquad X\in (\rho,\infty)\times \bbR\,, \\
m_{\varepsilon,22}^\rho(X) \geq m_{\varepsilon,21}^\rho(X) & = 0\,, \qquad X\in \bbR\times (\rho,\infty)\,.
\end{split} 
\end{equation}

Next, according to~\cite{DGJ1997}, it is natural to use the Hessian matrix of the convex function $\Phi_2$ to symmetrize~\eqref{ex12a}. We thus set
\begin{equation*}
S:=\frac{bd}{2}D^2\Phi_2 =
\begin{pmatrix}
ac& bc \\[1ex] bc & bd
\end{pmatrix}
\end{equation*}
and observe that $S$ is symmetric and positive definite by~\eqref{condabcd}. In addition, for all~${X\in\bbR^2}$,
\begin{equation*}
SM_\varepsilon^\rho(X) = \varepsilon S + \lambda_\varepsilon(X_{+}) SM^\rho(X)
\end{equation*}
with
\begin{equation*}
SM^\rho(X) =
\begin{pmatrix}
a^2 c \alpha_\rho(X_1) + b c^2 \alpha_\rho(X_2) & abc \alpha_\rho(X_1) + bcd \alpha_\rho(X_2) \\
& \\
abc \alpha_\rho(X_1) + bcd \alpha_\rho(X_2) & b^2 c \alpha_\rho(X_1) + bd^2 \alpha_\rho(X_2)
\end{pmatrix} \in {\bSym_2(\bbR)}\,.
\end{equation*}
Since $\trr(SM^\rho(X))\ge 0$ and
\begin{equation*}
\det\left(S M^\rho(X)\right) = \det(S) \det\left(M^\rho(X)\right) = bc (ad-bc)^2 \alpha_\rho(X_1) \alpha_\rho(X_2)\ge 0
\end{equation*}
by~\eqref{condabcd}, the matrix $SM^\rho(X)$ is positive semidefinite, so that the matrix $SM_{\varepsilon}^\rho(X)$ belongs to ${\bSPD_2(\bbR)}$ for all~${X\in\bbR^2}$ with
\begin{equation}\label{ex11c}
\left\langle SM_{\varepsilon}^\rho(X)\xi,\xi\right\rangle \ge \varepsilon \langle S\xi,\xi \rangle \ge \varepsilon \frac{\det(S)}{\trr(S)} |\xi|^2 = \varepsilon \frac{bc(ad-bc)}{ac+bd} |\xi|^2\,, \qquad \xi\in\bbR^2\,.
\end{equation}
\end{subequations}
According to the properties~\eqref{ex11}, we are now in a position to apply Lemma~\ref{lem.ap1} (with~$A=S$ and~${B=M_{\varepsilon}^\rho}$) and deduce that there is~${u_{\varepsilon}^\rho\in H^1(\Omega,\bbR^2)\cap L_{\infty,+}(\Omega,\bbR^2)}$ which solves~\eqref{ex12a} and satisfies~\eqref{ex12a'}. Moreover, it follows from the integral identity~\eqref{ex12a} (with $v=Su_\varepsilon^\rho\in H^1(\Omega,\bbR^2)$), \eqref{ex11c}, and the positive definiteness of~$S$,
\begin{equation*}
\langle S\xi,\xi\rangle \ge \frac{bc(ad-bc)}{ac+bd} |\xi|^2\,, \qquad \xi\in\bbR^2\,,
\end{equation*}
that
\begin{align*}
\|SU\|_2 \|u_\varepsilon^\rho\|_2 \ge \int_\Omega \langle SU, u_\varepsilon^\rho \rangle\ \rmd x & = \int_\Omega \left[\left\langle u_\varepsilon^\rho, S u_\varepsilon^\rho \right\rangle + \tau \sum_{i=1}^N\big\langle M_\varepsilon^\rho\left(u_\varepsilon^\rho\right) \partial_i u_\varepsilon^\rho, \partial_i S u_\varepsilon^\rho \big\rangle \right]\ \rmd x \\
& = \int_\Omega \left[\left\langle Su_\varepsilon^\rho, u_\varepsilon^\rho \right\rangle + \tau \sum_{i=1}^N\big\langle S M_\varepsilon^\rho\left(u_\varepsilon^\rho\right) \partial_i u_\varepsilon^\rho, \partial_i u_\varepsilon^\rho \big\rangle \right]\ \rmd x \\
& \ge \frac{bc(ad-bc)}{ac+bd} \left(\left\|u_\varepsilon^\rho\right\|_2^2 + \tau \varepsilon \left\|\nabla u_\varepsilon^\rho\right\|_2^2 \right)\,.
\end{align*}
Owing to~\eqref{condabcd}, we conclude that the estimates~\eqref{ex12d} and~\eqref{ex12e} are satisfied.

It remains to establish the estimate~\eqref{ex12c}. Let therefore $n\ge 2$. Since $u_{\varepsilon}^\rho$ belongs to~$H^1(\Omega,\bbR^2)\cap L_\infty(\Omega,\bbR^2)$, the vector field $D\Phi_n(u_\varepsilon^\rho)$ lies in $H^1(\Omega,\bbR^2)$ and we infer from~\eqref{ex12a} (with $v=D\Phi_n(u_\varepsilon^\rho)$) that
\begin{equation}\label{pro1}
\int_\Omega \left[\big\langle u_{\varepsilon}^\rho -U, D\Phi_n(u_\varepsilon^\rho) \big\rangle + \tau \sum_{i=1}^N\big\langle M_{\varepsilon }^\rho(u_{\varepsilon}^\rho) \partial_i u_{\varepsilon}^\rho, \partial_i D\Phi_n(u_\varepsilon^\rho)\big\rangle \right]\ \rmd x = 0\,.
\end{equation}
On the one hand, the convexity of $\Phi_n$ implies that
\begin{equation}\label{pro2}
\int_\Omega \big\langle u_{\varepsilon}^\rho -U, D\Phi_n(u_\varepsilon^\rho) \big\rangle\ \rmd x \ge \int_\Omega \left[\Phi_n(u_{\varepsilon}^\rho) - \Phi_n(U)\right]\ \rmd x =\clE_n(u_\varepsilon^\rho) - \clE_n(U)\,.
\end{equation}
On the other hand, using the symmetry and the positive semidefiniteness of the matrix $ D^2\Phi_n(u_\varepsilon^\rho)$, see Lemma~\ref{lem.p2}, we have
\begingroup
\allowdisplaybreaks
\begin{equation}\label{pro3}
\begin{split}
\tau &\sum_{i=1}^N\int_\Omega \big\langle M_{\varepsilon}^\rho(u_{\varepsilon}^\rho) \partial_i u_{\varepsilon}^\rho, \partial_i D\Phi_n(u_\varepsilon^\rho)\big\rangle\ \rmd x \\[1ex]
&=\tau \sum_{i=1}^N\int_\Omega \big\langle M_{\varepsilon }^\rho(u_{\varepsilon}^\rho) \partial_i u_{\varepsilon}^\rho, D^2\Phi_n(u_\varepsilon^\rho)\partial _i u_\varepsilon^\rho\big\rangle\ \rmd x \\[1ex]
&=\tau \sum_{i=1}^N\int_\Omega \left\langle D^2\Phi_n(u_\varepsilon^\rho) M_{\varepsilon }^\rho(u_{\varepsilon}^\rho) \partial_i u_{\varepsilon}^\rho, \partial _i u_\varepsilon^\rho\right\rangle\ \rmd x
\\
&=\tau\varepsilon \sum_{i=1}^N\int_\Omega \left\langle D^2\Phi_n(u_\varepsilon^\rho) \partial_i u_{\varepsilon}^\rho, \partial _i u_\varepsilon^\rho\right\rangle\ \rmd x\\[1ex]
&\quad + \tau \sum_{i=1}^N\int_\Omega \lambda_\varepsilon(u_\varepsilon^\rho) \left\langle D^2\Phi_n(u_\varepsilon^\rho) M^\rho(u_\varepsilon^\rho) \partial_i u_{\varepsilon}^\rho, \partial _i u_\varepsilon^\rho\right\rangle\ \rmd x\\[1ex]
&\geq \tau \sum_{i=1}^N\int_\Omega\lambda_\varepsilon(u_\varepsilon^\rho) \left\langle D^2\Phi_n(u_\varepsilon^\rho) M^\rho(u_\varepsilon^\rho) \partial_i u_{\varepsilon}^\rho, \partial _i u_\varepsilon^\rho\right\rangle\ \rmd x\,. 
\end{split}
\end{equation}
\endgroup
Since $S_n(u_\varepsilon^\rho) := D^2\Phi_n(u_\varepsilon^\rho) M(u_\varepsilon^\rho)$ is positive semidefinite by Lemma~\ref{lem.p3}, we further have
\begin{equation}\label{pro4}
\begin{split}
\tau &\sum_{i=1}^N\int_\Omega \lambda_\varepsilon(u_\varepsilon^\rho) \left\langle D^2\Phi_n(u_\varepsilon^\rho) M^\rho(u_\varepsilon^\rho) \partial_i u_{\varepsilon}^\rho, \partial _iu_\varepsilon^\rho\right\rangle\ \rmd x\\[1ex]
& = \tau \sum_{i=1}^N\int_\Omega \lambda_\varepsilon(u_\varepsilon^\rho) \left\langle D^2\Phi_n(u_\varepsilon^\rho) M(u_\varepsilon^\rho) \partial_i u_{\varepsilon}^\rho, \partial _iu_\varepsilon^\rho\right\rangle\ \rmd x \\[1ex]
&\qquad + \tau \sum_{i=1}^N\int_\Omega \lambda_\varepsilon(u_\varepsilon^\rho) \left\langle D^2\Phi_n(u_\varepsilon^\rho)\big[M^\rho(u_\varepsilon^\rho) - M(u_\varepsilon^\rho)\big] \partial_i u_{\varepsilon}^\rho, \partial _i u_\varepsilon^\rho\right\rangle\ \rmd x\\[1ex]
&\geq \tau \sum_{i=1}^N\int_\Omega \lambda_\varepsilon(u_\varepsilon^\rho) \left\langle D^2\Phi_n(u_\varepsilon^\rho) \big[M^\rho(u_\varepsilon^\rho) - M(u_\varepsilon^\rho)\big] \partial_i u_{\varepsilon}^\rho, \partial _i u_\varepsilon^\rho\right\rangle\ \rmd x\,. 
\end{split}
\end{equation}
Taking now advantage of the fact that $0 \le u_{\varepsilon,j}^\rho\leq \rho$ a.e. in $\Omega$ for $j=1,\, 2$ by~\eqref{ex12a'}, we further have
\begin{multline*}
\left| \tau \sum_{i=1}^N\int_\Omega \lambda_\varepsilon(u_\varepsilon^\rho) \left\langle D^2\Phi_n(u_\varepsilon^\rho) \big[M^\rho(u_\varepsilon^\rho) - M(u_\varepsilon^\rho)\big] \partial_i u_{\varepsilon}^\rho, \partial _i u_\varepsilon^\rho\right\rangle\ \rmd x\right|\\[1ex]
\begin{aligned}
& \le 2\tau \max\{a,\,b,\,c,\,d\} \left\|D^2\Phi_n\right\|_{L_\infty\left((0,\rho)^2\right)} \sum_{j=1}^2 \int_\Omega \lambda_\varepsilon(u_\varepsilon^\rho) \left|\alpha_\rho\left(u_{\varepsilon,j}^\rho\right) - u_{\varepsilon,j}^\rho\right|\, \left|\nabla u_\varepsilon^\rho\right|^2 \ \rmd x\\[1ex]
& \leq 8\tau \max\{a,\,b,\,c,\,d\}\kappa_n \rho^{n-2} \sum_{j=1}^2 \int_{\left\{\rho-1\,\leq\, u_{\varepsilon,j}^\rho\,\leq\,\rho\right\}} \frac{\left|\alpha_\rho\left(u_{\varepsilon,j}^\rho\right) - u_{\varepsilon,j}^\rho\right|}{1+\exp\left(\varepsilon u_{\varepsilon,j}^\rho\right)} \left|\nabla u_\varepsilon^\rho\right|^2 \ \rmd x\,,
\end{aligned}
\end{multline*}
where $\kappa_n\in\bbR$ is a positive constant such that
\begin{equation*}
\left|D^2\Phi_n(X)\right|\leq \kappa_n\left(X_1^{n-2}+X_2^{n-2}\right)\qquad\text{for all }X\in[0,\infty)^2.
\end{equation*}
Owing to the definition of $\alpha_\rho$, we further obtain
\begin{multline}\label{pro5}
\left| \tau \sum_{i=1}^N\int_\Omega \lambda_\varepsilon(u_\varepsilon^\rho) \left\langle D^2\Phi_n(u_\varepsilon^\rho) \big[M^\rho(u_\varepsilon^\rho) - M(u_\varepsilon^\rho)\big] \partial_i u_{\varepsilon}^\rho, \partial _i u_\varepsilon^\rho\right\rangle\ \rmd x\right| \\[1ex] 
\begin{aligned}
&\leq 8 \tau \max\{a,\,b,\,c,\,d\}\kappa_n \rho^{n-2} \sum_{j=1}^2 \int_{\left\{\rho-1\,\leq\,u_{\varepsilon,j}^\rho\,\leq\,\rho\right\}}\frac{\rho}{1+e^{\varepsilon(\rho-1)}} |\nabla u_\varepsilon^\rho|^2\ \rmd x\\[1ex]
&\leq 16e \tau \max\{a,\,b,\,c,\,d\}\kappa_n \rho^{n-1} e^{-\varepsilon\rho} \left\|\nabla u_\varepsilon^\rho\right\|_2^2\,.
\end{aligned}
\end{multline}
The desired estimate~\eqref{ex12c} is now a straightforward consequence of~\eqref{pro1}-\eqref{pro5}.
\end{proof}



\subsection{A regularised system: \texorpdfstring{$\rho\to\infty$}{rho to infty}}\label{sec2.3}



We next study the cluster points of the family~${\{u_{\varepsilon}^\rho\,:\, \rho\geq \max\{1,\|F\|_\infty,\|G\|_\infty\}\}}$ provided in Lemma~\ref{lem.ex2}, as $\rho\to\infty$, the parameter $\varepsilon\in (0,1)$ being held fixed.


\begin{lemm}\label{lem.ex3}
Given $\tau>0$, $U=(F,G)\in L_{\infty,+}(\Omega,\bbR^2)$, and $\varepsilon\in (0,1)$, there exist a sequence~$(\rho_l)_{l\,\ge\,1}$ and a function~${u_{\varepsilon} = (u_{\varepsilon,1}, u_{\varepsilon,2}) \in H^1(\Omega,\bbR^2)\cap L_{\infty,+}(\Omega,\bbR^2)}$ such that $\rho_l\to\infty$ and
\begin{align}
u_\varepsilon^{\rho_l}&\to u_\varepsilon\qquad\text{in }L_p\left(\Omega,\bbR^2\right) \text{for all }p\in [1,\infty)\text{ and pointwise a.e. in }\Omega\,,\label{conv1}
\\[1ex]
\nabla u_\varepsilon^{\rho_l}&\rightharpoonup \nabla u_\varepsilon\quad\, \text{in }L_2\left(\Omega,\bbR^{2N}\right)\,.\label{conv3}
\end{align}
Moreover, $u_\varepsilon$ solves the equation
\begin{align}
\int_\Omega \bigg[\langle u_{\varepsilon}, v \rangle + \tau \sum_{i=1}^N\big\langle M_{\varepsilon}(u_{\varepsilon}) \partial_i u_{\varepsilon}, &\partial_i v \big\rangle \bigg]\ \rmd x \label{ex12aa}\\
&=\int_\Omega \langle U, v \rangle\ \rmd x\,, \quad v\in H^1\left(\Omega,\bbR^2\right)\,, \nonumber
\end{align}
where
\begin{equation*}
M_{\varepsilon}(X) = \left(m_{\varepsilon, jk}(X)\right)_{1\,\le\,j,\,k\,\le\,2}:= \varepsilon I_2 + \lambda_\varepsilon(X_{+}) M(X),
\end{equation*}
with $M(X)$ defined in~\eqref{x1}, and, for each $n\ge 2$, we have
\begin{equation}\label{ex12ba}
\clE_n(u_\varepsilon)\le \clE_n(U)\,. 
\end{equation}
Furthermore,
\begin{equation}\label{ex12bac}
\min\left\{ 1, \frac{c}{d}\right\} \|u_{\varepsilon,1}+u_{\varepsilon,2}\|_\infty \le \max\left\{1,\frac{a}{b} \right\}\|F+G\|_\infty\,. 
\end{equation}
\end{lemm}


\begin{proof}
Recalling~\eqref{ex12d}-\eqref{ex12e}, we deduce that $(u_\varepsilon^\rho)_\rho$ is bounded in $H^1(\Omega,\bbR^2)$. Moreover, since
\begin{equation}\label{x6}
\frac{\varepsilon^n z^n}{n!} \le e^{\varepsilon z}\,, \qquad z\in [0,\infty)\,, \quad n\ge 1\,, 
\end{equation}
the estimates~\eqref{ex12e} and~\eqref{ex12c}, along with Lemma~\ref{lelfb}, ensure that $(u_\varepsilon^\rho)_\rho$ is bounded in $L_n(\Omega,\bbR^2)$ for any integer~$n\geq 2$ (with an $\varepsilon$-dependent bound). We may then use a Cantor diagonal process, together with Rellich--Kondrachov' theorem and an interpolation argument, to deduce the convergences~\eqref{conv1} and~\eqref{conv3} along a sequence~$\rho_l\to\infty$, as well as the componentwise non-negativity of~$u_\varepsilon$.

Since $\Phi_n$ is convex on $[0,\infty)^2$ for all $n\geq 2$, see Lemma~\ref{lem.p2}, it follows from the relations~\eqref{ex12e}, \eqref{ex12c}, \eqref{conv1}, and~\eqref{x6} that~\eqref{ex12ba} holds true. Using once more Lemma~\ref{lelfb}, we infer from~\eqref{ex12ba} that
\begin{equation*}
\left\|c u_{\varepsilon,1} + d u_{\varepsilon, 2}\right\|_n \le \frac{d}{b} \|a F + b G\|_n
\end{equation*}
for all $n\ge 2$. Passing to the limit $n\to\infty$ in the above inequality, we deduce that the function $u_\varepsilon\in L_\infty(\Omega,\bbR^2)$ satisfies~\eqref{ex12bac}.

Let us now consider $v\in H^1(\Omega,\bbR^2)$. Since~\eqref{conv1} and~\eqref{conv3} imply that
\begin{equation*}
\lim_{l\,\to\,\infty} \int_\Omega \left\langle u_{\varepsilon}^{\rho_l}, v \right\rangle\ \rmd x = \int_\Omega \left\langle u_{\varepsilon}, v \right\rangle\ \rmd x
\quad\text{and}\quad \lim_{l\,\to\,\infty} \int_\Omega \left\langle \partial_i u_{\varepsilon}^{\rho_l}, \partial_i v \right\rangle \ \rmd x = \int_\Omega \left\langle \partial_i u_{\varepsilon}, \partial_i v \right\rangle \ \rmd x
\end{equation*}
for $1\le i \le N$, the identity~\eqref{ex12aa} is satisfied provided that
\begin{equation}\label{ex12ac}
\lim_{l\,\to\,\infty} \int_\Omega \lambda_\varepsilon(u_{\varepsilon}^{\rho_l}) \big\langle M^{\rho_l}(u_{\varepsilon}^{\rho_l}) \partial_i u_{\varepsilon}^{\rho_l}, \partial_i v \big\rangle \ \rmd x = \int_\Omega \lambda_\varepsilon(u_{\varepsilon})\left\langle M(u_{\varepsilon}) \partial_i u_{\varepsilon}, \partial_i v \right\rangle\ \rmd x 
\end{equation}
for each $1\leq i\leq N$. To prove~\eqref{ex12ac}, we observe that, for $1\leq i\leq N$ and $j\in\{1,2\}$,
\begin{equation}\label{x7}
\int_\Omega \lambda_\varepsilon(u_\varepsilon^{\rho_l}) \big\langle M^{\rho_l}(u_{\varepsilon}^{\rho_l}) \partial_i u_{\varepsilon}^{\rho_l}, \partial_i v \big\rangle \ \rmd x = \int_\Omega \lambda_\varepsilon(u_\varepsilon^{\rho_l}) \left\langle M^{\rho_l}(u_{\varepsilon}^{\rho_l})^t \partial_i v, \partial_i u_{\varepsilon}^{\rho_l} \right\rangle \ \rmd x 
\end{equation}
with
\begin{align*}
\left| \lambda_\varepsilon(u_\varepsilon^{\rho_l}) \sum_{k=1}^2 m_{kj}^{\rho_l}(u_{\varepsilon}^{\rho_l}) \partial_i v_k \right| & \le 2\max\{a,b,c,d\} \frac{u_{\varepsilon,1}^{\rho_l} + u_{\varepsilon,2}^{\rho_l}}{1+ \exp\left[\varepsilon\left(u_{\varepsilon,1}^{\rho_l} + u_{\varepsilon,2}^{\rho_l}\right)\right]} |\partial_i v| \\
&\le \frac{2\max\{a,b,c,d\} }{\varepsilon} |\partial_i v| \qquad\text{a.e. in $\Omega$}\,,
\end{align*}
by the definition of $\lambda_\varepsilon$ and~\eqref{x6} (with $n=1$), and
\begin{equation*}
\lim_{l\,\to\,\infty} \lambda_\varepsilon(u_\varepsilon^{\rho_l}) \sum_{k=1}^2 m_{kj}^{\rho_l}(u_{\varepsilon}^{\rho_l}) \partial_i v_k = \lambda_\varepsilon(u_\varepsilon) \sum_{k=1}^2 m_{ kj}(u_{\varepsilon}) \partial_i v_k \quad\text{a.e. in $\Omega$}\,,
\end{equation*}
by~\eqref{momar}, the pointwise almost everywhere convergence in $\Omega$ established in~\eqref{conv1}, and the properties of $\alpha_{\rho_l}$. Lebesgue's dominated convergence theorem then guarantees that
\begin{equation*}
\lim_{l\,\to\,\infty} \left\| \lambda_\varepsilon(u_\varepsilon^{\rho_l}) \sum_{k=1}^2 m_{kj}^{\rho_l}(u_{\varepsilon}^{\rho_l}) \partial_i v_k - \lambda_\varepsilon(u_\varepsilon) \sum_{k=1}^2 m_{kj}(u_{\varepsilon}) \partial_i v_k \right\|_2 = 0\,.
\end{equation*}
Combining the above convergence with~\eqref{conv3} allows us to pass to the limit $l\to\infty$ in~\eqref{x7} and find
\begin{align*}
\lim_{l\,\to\,\infty} \int_\Omega \lambda_\varepsilon(u_\varepsilon^{\rho_l}) \left\langle M^{\rho_l}(u_{\varepsilon}^{\rho_l}) \partial_i u_{\varepsilon}^{\rho_l}, \partial_i v \right\rangle \ \rmd x & = \int_\Omega \lambda_\varepsilon(u_\varepsilon) \left\langle M(u_{\varepsilon})^t \partial_i v, \partial_i u_{\varepsilon} \right\rangle \ \rmd x \\
& = \int_\Omega \lambda_\varepsilon(u_\varepsilon) \left\langle M(u_{\varepsilon}) \partial_i u_{\varepsilon}, \partial_i v \right\rangle \ \rmd x
\end{align*}
for $1\le i \le N$, which proves~\eqref{ex12ac}. We have thus shown that $u_\varepsilon$ solves~\eqref{ex12aa} and thereby completed the proof of Lemma~\ref{lem.ex3}.
\end{proof}

We next show that the entropy functional $\clE_1$ evaluated at the function~$u_\varepsilon$ identified in Lemma~\ref{lem.ex3} is dominated by $\clE_1(U)$ and that the associated dissipation term~${\clE_1(U)-\clE_1(u_\varepsilon)}$ provides a control on the gradient of $u_\varepsilon$ which is essential when considering the limit $\varepsilon\to 0$.


\begin{lemm}\label{lem.ex4}
Let $\tau>0$, $U=(F,G)\in L_{\infty,+}(\Omega,\bbR^2)$, and $\varepsilon\in (0,1)$. The function
\begin{equation*}
u_{\varepsilon}=(u_{\varepsilon,1},u_{\varepsilon,2})\in H^1\left(\Omega,\bbR^2\right)\cap L_{\infty,+}\left(\Omega,\bbR^2\right)
\end{equation*}
identified in Lemma~\ref{lem.ex3} satisfies
\begin{equation*}
\clE_1(u_{\varepsilon}) + \frac{ \tau}{a} \int_\Omega \lambda_\varepsilon(u_{\varepsilon}) \big[\left|\nabla (a u_{\varepsilon,1} + \Theta_1 u_{\varepsilon,2})\right|^2 + \Theta_2 \left|\nabla u_{\varepsilon,2}\right|^2 \big]\ \rmd x \le \clE_1(U)\,.
\end{equation*}
\end{lemm}


\begin{proof}
Let $\eta\in (0,1)$. Then $\left(\ln{(u_{\varepsilon,1}+\eta)}, (b^2/ad)\ln{(u_{\varepsilon,2}+\eta) } \right)\in H^1(\Omega,\bbR^2)$ and we infer from~\eqref{ex12aa} that
\begin{equation}\label{ex19}
0 = \int_\Omega \Big[(u_{\varepsilon,1} - U_1) \ln{(u_{\varepsilon,1}+\eta)} + \frac{b^2}{ad} (u_{\varepsilon,2} - U_2)\ln{(u_{\varepsilon,2}+\eta)} \Big]\ \rmd x +D(\eta)\,,
\end{equation}
where
\begin{align*}
D(\eta)&:= \tau \int_\Omega \sum_{i=1}^N\big(m_{\varepsilon,11}(u_{\varepsilon}) \partial_i u_{\varepsilon,1} + m_{\varepsilon,12}(u_{\varepsilon}) \partial_i u_{\varepsilon,2} \big)
\frac{\partial_i u_{\varepsilon,1}}{u_{\varepsilon,1}+\eta}\ \rmd x \\[1ex]
& \qquad + \frac{\tau b^2}{ad} \int_\Omega \sum_{i=1}^N\big(m_{\varepsilon,21}(u_{\varepsilon}) \partial_i u_{\varepsilon,1} + m_{\varepsilon,22}(u_{\varepsilon}) \partial_i u_{\varepsilon,2} \big)
\frac{\partial_i u_{\varepsilon,2}}{u_{\varepsilon,2}+\eta}\ \rmd x \,.
\end{align*}
Since $L(r)=r\ln r-r+1$ is convex on $[0,\infty)$ with $L'(r)=\ln{r}$, the first term on the right-hand side of~\eqref{ex19} can be estimated as follows:
\begin{multline*}
\int_\Omega \left[(u_{\varepsilon,1} - U_1) \ln{(u_{\varepsilon,1}+\eta)} + \frac{b^2}{ad} (u_{\varepsilon,2} - U_2)\ln{(u_{\varepsilon,2}+\eta)} \right]\ \rmd x \\[1ex]
\begin{aligned}
& \ge \int_\Omega \left[\big(L(u_{\varepsilon,1}+\eta) - L(U_1+\eta)\big) + \frac{b^2}{ad}\big(L(u_{\varepsilon,2}+\eta) - L(U_2+\eta)\big) \right]\ \rmd x \\[1ex]
& = \clE_1((u_{\varepsilon,1}+\eta,u_{\varepsilon,2}+\eta))- \clE_1((U_{1}+\eta,U_{2}+\eta))\,.
\end{aligned}
\end{multline*}
Using the continuity of $\Phi_1$ and the boundedness of $u_\varepsilon$, see~\eqref{ex12bac}, we deduce that
\begin{multline}\label{ex20}
\liminf_{\eta\,\to\,0} \int_\Omega \left[(u_{\varepsilon,1} - U_1) \ln{(u_{\varepsilon,1}+\eta)} + \frac{b^2}{ad} (u_{\varepsilon,2} - U_2)\ln{(u_{\varepsilon,2}+\eta)} \right]\ \rmd x \\
\ge \clE_1(u_{ \varepsilon})- \clE_1(U)\,. 
\end{multline}
Next, recalling the definition of the matrix $M_\varepsilon$, see Lemma~\ref{lem.ex3}, we have
\begin{align*}
D(\eta) & = \tau\varepsilon\int_\Omega
\left(\frac{|\nabla u_{\varepsilon,1}|^2}{u_{\varepsilon,1}+\eta} +\frac{b^2}{ad} \frac{|\nabla u_{\varepsilon,2}|^2}{u_{\varepsilon,2}+\eta} \right)\ \rmd x \\
&\qquad + \tau \int_\Omega \lambda_\varepsilon(u_{\varepsilon}) \left(\frac{u_{\varepsilon,1}}{u_{\varepsilon,1}+\eta} - 1 + 1 \right) \nabla u_{\varepsilon,1}\cdot \nabla(a u_{\varepsilon,1} +b u_{\varepsilon,2}) \ \rmd x \\
&\qquad + \frac{\tau b^2}{ad} \int_\Omega \lambda_\varepsilon(u_{\varepsilon}) \left(\frac{ u_{\varepsilon,2}}{u_{\varepsilon,2}+\eta} - 1 + 1 \right)
\nabla u_{\varepsilon,2}\cdot \nabla(c u_{\varepsilon,1} +d u_{\varepsilon,2}) \ \rmd x \\
&=\tau\varepsilon\int_\Omega \left(\frac{|\nabla u_{\varepsilon,1}|^2}{u_{\varepsilon,1}+\eta} +\frac{b^2}{ad} \frac{|\nabla u_{\varepsilon,2}|^2}{u_{\varepsilon,2}+\eta} \right)\ \rmd x \\
& \qquad + \frac{\tau}{a} \int_\Omega\lambda_\varepsilon(u_{\varepsilon}) \big[|\nabla(a u_{\varepsilon,1}+\Theta_1 u_{\varepsilon,2})|^2 + \Theta_2|\nabla u_{\varepsilon,2} |^2 \big]\ \rmd x\\
& \qquad - J_1(\eta) - J_2(\eta)\,,
\end{align*}
where
\begin{align*}
J_1(\eta) & := \tau \int_\Omega \frac{\eta\lambda_\varepsilon(u_{\varepsilon}) }{u_{\varepsilon,1}+\eta} \nabla u_{\varepsilon,1} \cdot \nabla(a u_{\varepsilon,1} + b u_{\varepsilon,2}) \ \rmd x\,, \\
J_2(\eta) & := \frac{\tau b^2}{ad} \int_\Omega \frac{\eta\lambda_\varepsilon(u_{\varepsilon}) }{u_{\varepsilon,2}+\eta} \nabla u_{\varepsilon,2} \cdot \nabla(c u_{\varepsilon,1} + d u_{\varepsilon,2}) \ \rmd x\,.
\end{align*}
Since $u_{\varepsilon}\in H^1(\Omega,\bbR^2)$ satisfies $\nabla u_{\varepsilon,j}=0$ a.e. on the level set $\{x\in\Omega\ : \ u_{\varepsilon,j}=0\}$ for~$j\in\{1,2\}$, we have
\begin{align*}
\lim_{\eta\,\to\,0} \frac{\eta\lambda_\varepsilon(u_{\varepsilon}) }{u_{\varepsilon,j}+\eta} \nabla u_{\varepsilon,j}=0\qquad &\text{ a.e. in $\Omega$}\,, \\
 \left| \frac{\eta\lambda_\varepsilon(u_{\varepsilon}) }{u_{\varepsilon,j}+\eta} \nabla u_{\varepsilon,j} \right|
\leq |\nabla u_{\varepsilon,j}|\qquad &\text{ a.e. in $\Omega$}\,.
\end{align*}
Lebesgue's dominated convergence theorem ensures now that
\begin{equation*}
\lim_{\eta\,\to\,0} \left(J_1(\eta) + J_2(\eta) \right) = 0\,.
\end{equation*}
This shows that
\begin{equation}
\liminf_{\eta\,\to\,0} D(\eta) \ge \frac{\tau}{a} \int_\Omega\lambda_\varepsilon(u_{\varepsilon}) \big[\left|\nabla \left(a u_{\varepsilon,1}+\Theta_1 u_{\varepsilon,2}\right)\right|^2 + \Theta_2\left|\nabla u_{\varepsilon,2}\right|^2 \big]\ \rmd x\,. \label{ex23}
\end{equation}
Passing to the limit $\eta\to 0$ in~\eqref{ex19}, we get the desired estimate in view of~\eqref{ex20} and~\eqref{ex23}.
\end{proof}



\subsection{A regularised system: \texorpdfstring{$\varepsilon\to 0$}{varepsilon to 0}}\label{sec2.4}



We complete this section with the proof of Proposition~\ref{P:1}.

\begin{proof}[Proof of Proposition~\ref{P:1}]
Consider $\tau>0$ and $U=(F,G)\in L_{\infty,+}(\Omega,\bbR^2)$. Given~$\varepsilon\in (0,1)$, let
\begin{equation*}
u_\varepsilon =(u_{\varepsilon,1},u_{\varepsilon,2}) \in H^1\left(\Omega,\bbR^2\right)\cap L_{\infty,+}\left(\Omega,\bbR^2\right)
\end{equation*}
denote the weak solution to~\eqref{ex12aa} provided by Lemma~\ref{lem.ex3}. According to~\eqref{ex12bac},
\begin{equation}\label{x8}
\max\{\|u_{\varepsilon,1}\|_\infty, \| u_{\varepsilon,2}\|_\infty \} \le\|u_{\varepsilon,1}+ u_{\varepsilon,2}\|_\infty \le R_0 := \frac{d}{b} \frac{\max\{a,b\}}{\min\{c,d\}} \|F+G\|_\infty\,. 
\end{equation}
Hence,
\begin{equation*}
\lambda_\varepsilon(u_\varepsilon) \ge \frac{2}{1+e^{R_0}}\,,
\end{equation*}
a lower bound which, together with Lemma~\ref{lem.ex4} and the non-negativity of $\clE_1$, ensures that
\begin{equation}\label{x9}
(\nabla u_{\varepsilon})_{\varepsilon\,\in\,(0,1)} \text{ is bounded in }L_2\left(\Omega,\bbR^{2N}\right). 
\end{equation}
We now infer from~\eqref{x8}, \eqref{x9}, Rellich--Kondrachov' theorem, an interpolation argument, and a Cantor diagonal process that there exist a function
\begin{equation*}
u=(f,g)\in H^1\left(\Omega,\bbR^2\right)\cap L_{\infty,+}\left(\Omega,\bbR^2\right)
\end{equation*}
and a sequence $(\varepsilon_l)_{l\,\ge\,1}$, with $\varepsilon_l\to 0$, such that
\begin{align}
u_{\varepsilon_l}&\to u\qquad\text{in }L_p\left(\Omega,\bbR^2\right)\text{ for all }p\in [1,\infty)\,,\label{conv1'}
\\[1ex] u_{\varepsilon_l}&\overset{*}\rightharpoonup u\qquad \text{in }L_\infty\left(\Omega,\bbR^2\right)\,,\label{conv2'}
\\[1ex]
\nabla u_{\varepsilon_l}&\rightharpoonup \nabla u\quad\: \text{in }L_2\left(\Omega,\bbR^{2N}\right)\,.\label{conv3'}
\end{align}
An immediate consequence of~\eqref{ex12ba} and~\eqref{conv1'} is the estimate~\eqref{ex2}. As~${\sqrt{\lambda_{\varepsilon_l}(u_{\varepsilon_l})}\to \!1}$ in~${L_\infty(\Omega)}$ by~\eqref{UL} and~\eqref{x8}, we conclude together with~\eqref{conv3'} that
\begin{align*}
\sqrt{\lambda_{\varepsilon_l}(u_{\varepsilon_l})} \nabla \big(a u_{\varepsilon_l,1} + \Theta_1 u_{\varepsilon_l,2} \big) &
\rightharpoonup \nabla\big (a u_1 + \Theta_1 u_2\big) \qquad \text{in }L_2\left(\Omega,\bbR^{N}\right)\,, \\[1ex]
\sqrt{\Theta_2 \lambda_{\varepsilon_l}(u_{\varepsilon_l})} \nabla u_{\varepsilon_l,2} & \rightharpoonup \sqrt{\Theta_2}\nabla u_2 \mkern87mu \text{in }L_2\left(\Omega,\bbR^{N}\right)\,.
\end{align*}
Moreover, the $L_\infty$-bound~\eqref{x8} and the convergence~\eqref{conv1'} imply that
\begin{equation*}
\liminf_{l\,\to\,\infty} \clE_1(u_{\varepsilon_l}) \geq \clE_1(u)\,,
\end{equation*}
and the estimate~\eqref{ex2b} is now obtained by passing to $\liminf$ in the inequality reported in Lemma~\ref{lem.ex4} (with $\varepsilon$ replaced by $\varepsilon_l$).

Finally,~\eqref{conv1'}, along with~\eqref{x8} and the convergence property
\begin{equation*}
\lim_{\varepsilon\,\to\,0} \left| m_{\varepsilon,jk}(X) - m_{jk}(X) \right| = 0\,,
\end{equation*}
which is uniform with respect to $X\in [0,R_0]^2 $ and $1 \le j,k\le 2$, enables us to use Lebesgue's dominated convergence theorem to show that, for $v=(\varphi,\psi)\in H^1(\Omega,\bbR^2)$,
\begin{equation*}
\lim_{l\,\to\,\infty} \big\| M_{\varepsilon_l}(u_{\varepsilon_l})^t \partial_i v - M(u)^t \partial_i v \big\|_2 = 0\,, \qquad 1\le i \le N\,.
\end{equation*}
Together with~\eqref{conv1'} and~\eqref{conv3'}, the above convergence allows us to let $\varepsilon_l\to 0$ in~\eqref{ex12aa} and conclude that $u=(f,g)$ satisfies~\eqref{ex1}. This completes the proof of Proposition~\ref{P:1}.
\end{proof}



\section{Existence of bounded weak solutions}\label{sec3}



This section is devoted to the proof of Theorem~\ref{ThBWS}, which relies on rather classical arguments, besides the estimates derived in Proposition~\ref{P:1}, and proceeds along the lines of the proof of~\cite[Theorem~1.2]{LM22}. As a first step, we use Proposition~\ref{P:1} to construct a family of piecewise constant functions $(u^\tau)_{\tau\,\in\,(0,1)}$ starting from the initial condition $(f^{in},g^{in})\in L_{\infty,+}(\Omega,\bbR^2)$. More precisely, for~${\tau\in (0,1)}$, we set $u^\tau(0):=u_0^\tau $ and
\begin{equation}\label{x01}
u^\tau(t)= u^\tau_{l}\,, \qquad t\in ((l-1)\tau, l\tau]\,, \qquad l\in\bbN\setminus\{0\}\,,
\end{equation}
where the sequence $(u_{l}^\tau)_{l\,\geq\,0}$ is defined as follows:
\begin{equation}
\begin{split}
&u_0^\tau = u^{in} := (f^{in},g^{in})\in L_{\infty,+}\left(\Omega,\bbR^2\right)\,, \\
&u_{l+1}^\tau =(f_{l+1}^\tau,g_{l+1}^\tau)\in H^1\left(\Omega,\bbR^2\right)\cap L_{\infty,+}\left(\Omega,\bbR^2\right) \;\text{is the solution to~\eqref{ex1}} \\
&\text{with } U=u_l^\tau=(f_l^\tau,g_l^\tau) \text{ constructed in Proposition~\ref{P:1} for }l\ge 0\,.
\end{split}\label{x02}
\end{equation}
In order to establish Theorem~\ref{ThBWS}, we show that the family $(u^\tau)_{\tau\,\in\,(0,1)}$ defined in~\eqref{x02} converges along a subsequence~${\tau_j\to0}$ towards a pair $u=(f,g)$ which fulfills all the requirements of Theorem~\ref{ThBWS}.

%\medskip

Below, $C$ and $(C_l)_{l\,\ge\,0}$ denote various positive constants depending only on~${(a, b, c, d)}$ and $u^{in}$. Dependence upon additional parameters will be indicated explicitly.

\begin{proof}[Proof of Theorem~\ref{ThBWS}]
Let $\tau\in (0,1)$ and let $u^\tau$ be defined in~\eqref{x01}-\eqref{x02}. Given~${l\ge 0}$, we infer from Proposition~\ref{P:1} that
\begin{subequations}\label{x03}
\begin{multline}\label{x03a}
\int_\Omega \Big(f_{l+1}^\tau \varphi + \tau f_{l+1}^\tau \nabla\left[af_{l+1}^\tau + b g_{l+1}^\tau\right] \cdot\nabla \varphi \Big)\ \rmd x \\
= \int_\Omega f_{l}^\tau \varphi\ \rmd x\,, \quad \varphi\in H^1(\Omega)\,, 
\end{multline}
\begin{multline}\label{x03b}
\int_\Omega \Big(g_{l+1}^\tau \psi + \tau g_{l+1}^\tau \nabla\left[cf_{l+1}^\tau + dg_{l+1}^\tau\right] \cdot\nabla\psi \Big)\ \rmd x \\
= \int_\Omega g_{l}^\tau \psi\ \rmd x\,, \quad \psi\in H^1(\Omega)\,. 
\end{multline}
\end{subequations}
Moreover,
\begin{equation}\label{x04}
\clE_n(u_{l+1}^\tau) \le \clE_n(u_{l}^\tau)\qquad\text{for }n\ge 2, 
\end{equation}
and we also have
\begin{equation}\label{x05}
\clE_1(u_{l+1}^\tau) +\frac{ \tau}{a} \int_\Omega \left[\left|\nabla \left(a f_{l+1}^\tau + \Theta_1g_{l+1}^\tau\right)\right|^2+\Theta_2\left|\nabla g_{l+1}^\tau\right|^2\right]\ \rmd x \le\clE_1(u_{l}^\tau)\,. 
\end{equation}
It readily follows from~\eqref{x01}, \eqref{x02}, \eqref{x04}, and~\eqref{x05} that, for $t>0$,
\begin{equation}\label{x06}
\clE_n(u^\tau(t)) \le \clE_n (u^{in}) \,, \quad \ n\ge 2\,, 
\end{equation}
and
\begin{equation}\label{x07}
\clE_1(u^\tau(t)) +\frac{ 1}{a} \ \int_0^t \int_\Omega \big[\left|\nabla \left(a f^\tau + \Theta_1g^\tau\right)\right|^2+\Theta_2\left|\nabla g^\tau\right|^2\big]\ \rmd x\rmd s \le \clE_1(u^{in})\,. 
\end{equation}

An immediate consequence of~\eqref{x06} and Lemma~\ref{lelfb} is the estimate
\begin{equation*}
\left\|f^\tau(t)+g^\tau(t)\right\|_n \le \frac{d}{b}\frac{\max\{a,\,b\}}{\min\{c,\, d\}}\left\|f^{in}+g^{in}\right\|_n\,, \quad n\ge 2\,, \ t>0\,.
\end{equation*}
Letting $n\to\infty$ in the above inequality gives
\begin{equation}\label{x08}
\left\|f^\tau(t)+g^\tau(t)\right\|_\infty \le C_2 := \frac{d}{b}\frac{\max\{a,\,b\}}{\min\{c,\, d\}} \left\|f^{in}+g^{in}\right\|_\infty\,, \quad t>0\,.
\end{equation}
Also, taking advantage of the non-negativity of $\clE_1$, we deduce from~\eqref{x07} that
\begin{equation}\label{x09}
\int_0^t \big[\left\|\nabla f^\tau(s)\right\|_2^2 + \left\|\nabla g^\tau(s)\right\|_2^2 \big]\ \rmd s\le C_3
:= \frac{a^2+2\left(\Theta_2+\Theta_1^2\right)}{a\Theta_2} \clE_1(u^{in})
\end{equation} 
for $t>0$.

Next, for $l\ge 1$ and $t\in ((l-1)\tau,l\tau]$, we deduce from~\eqref{x03a}, \eqref{x08}, and H\"older's inequality that, for $\varphi\in H^1(\Omega)$,
\begin{align*}
\left| \int_\Omega \big(f^{\tau}(t+\tau) - f^\tau(t) \big) \varphi\ \rmd x \right| & =
\left| \int_{l\tau}^{(l+1) \tau} \int_\Omega f_{l+1}^\tau \nabla\left[af_{l+1}^\tau + b g_{l+1}^\tau\right] \cdot\nabla \varphi\ \rmd x\rmd s \right| \\
& \le \int_{l\tau}^{(l+1) \tau} \left\| f^\tau(s)\right\|_\infty \left\|\nabla\left[af^\tau(s) + b g^\tau(s)\right]\right\|_2 \|\nabla\varphi\|_2\ \rmd s \\
& \le C_2 \|\nabla\varphi\|_2 \int_{l\tau}^{(l+1) \tau} \left\|\nabla\left[af^\tau(s) + b g^\tau(s)\right]\right\|_2\ \rmd s \,.
\end{align*}
A duality argument then gives
\begin{equation*}
\left\| f^{\tau}(t+\tau) - f^\tau(t) \right\|_{(H^1)'} \le C_2 \int_{l\tau}^{(l+1) \tau} \left\|\nabla\left[af^\tau(s) + b g^\tau(s)\right]\right\|_2\ \rmd s
\end{equation*}
for $t\in ((l-1)\tau,l\tau]$ and $l\ge 1$. Now, for $l_0\ge 2$ and $T\in ((l_0-1)\tau,l_0\tau]$, the above inequality, along with H\"older's inequality, entails that
\begin{align*}
\int_0^{T-\tau} \left\| f^{\tau}(t+\tau) - f^\tau(t) \right\|_{(H^1)'}^2\ \rmd t & \le \int_0^{(l_0-1)\tau} \left\| f^{\tau}(t+\tau) - f^\tau(t)\right\|_{(H^1)'}^2\ \rmd t \\
& =\sum_{l=1}^{l_0-1} \int_{(l-1)\tau}^{l\tau} \left\| f^{\tau}(t+\tau) - f^\tau(t) \right\|_{(H^1)'}^2\ \rmd t \\
& \le C_2^2 \tau \sum_{l=1}^{l_0-1} \left(\int_{l\tau}^{(l+1) \tau} \left\|\nabla\left[af^\tau(s) + b g^\tau(s)\right]\right\|_2\ \rmd s \right)^2\\
& \le C_2^2 \tau^2 \sum_{l=1}^{l_0-1} \int_{l\tau}^{(l+1) \tau} \left\|\nabla\left[af^\tau(s) + b g^\tau(s)\right]\right\|_2^2\ \rmd s \\
& \le C_2^2 \tau^2 \int_{0}^{l_0 \tau}\left\|\nabla\left[af^\tau(s) + b g^\tau(s)\right]\right\|_2^2\ \rmd s\,.
\end{align*}
We then use~\eqref{x09} (with $t=l_0\tau$) and Young's inequality to obtain
\begin{multline}\label{x10}
\int_0^{T-\tau} \left\| f^{\tau}(t+\tau) - f^\tau(t) \right\|_{(H^1)'}^2\ \rmd t\\
\begin{aligned}
& \le C_2^2 \tau^2
\int_{0}^{l_0\tau} \left(2 a^2 \left\|\nabla f^\tau(s)\right\|_2^2 + 2b^2 \left\|\nabla g^\tau(s)\right\|_2^2 \right)\ \rmd s \\
& \le C_4 \tau^2 \,, 
\end{aligned}
\end{multline}
with $C_4 := 2 (a^2+b^2)^2 C_2^2 C_3$. Similarly,
\begin{equation}\label{x11}
\int_0^{T-\tau} \left\| g^{\tau}(t+\tau) - g^\tau(t) \right\|_{(H^1)'}^2\ \rmd t \le C_5 \tau^2 \,, 
\end{equation}
with $C_5 := 2 (c^2+d^2) C_2^2 C_3$.

According to Rellich--Kondrachov' theorem, $H^1(\Omega,\bbR^2)$ is compactly embedded in~${L_2(\Omega,\bbR^2)}$, while~${L_2(\Omega,\bbR^2)}$ is continuously and compactly embedded in~$H^1(\Omega,\bbR^2)'$. Gathering~\eqref{x08}-\eqref{x11}, we infer from~\cite[Theorem~1]{DJ2012} that, for any~${T>0}$,
\begin{equation}\label{x12}
(u^\tau)_{\tau\,\in\,(0,1)} \;\text{ is relatively compact in }\; L_2\left((0,T)\times\Omega,\bbR^2\right)\,. 
\end{equation}
Owing to~\eqref{x08}, \eqref{x09}, and~\eqref{x12}, we may use a Cantor diagonal argument to find a function
\begin{equation*}
u=(f,g)\in L_{\infty,+}\left((0,\infty)\times\Omega,\bbR^2\right)
\end{equation*}
and a sequence $(\tau_m)_{m\,\ge\,1}$, $\tau_m\to 0$, such that, for any $T>0$ and $p\in [1,\infty)$,
\begin{equation}\label{x13}
\begin{split}
u^{\tau_m} & \rightarrow u \quad\text{in}\quad L_p\left((0,T)\times\Omega,\bbR^2\right)\,, \\
u^{\tau_m} & \stackrel{*}{\rightharpoonup} u \quad\text{in}\quad L_\infty\left((0,T)\times\Omega,\bbR^2\right)\,, \\
u^{\tau_m} & \rightharpoonup u \quad\text{in}\quad L_2\left((0,T),H^1\left(\Omega,\bbR^2\right)\right)\,.
\end{split}
\end{equation}
In addition, the compact embedding of $L_2(\Omega,\bbR^2)$ in $H^1(\Omega,\bbR^2)'$, along with~\eqref{x06} for $n=2$, \eqref{x10}, and~\eqref{x11}, allows us to apply once more~\cite[Theorem~1]{DJ2012} to conclude that
\begin{equation}\label{x14}
u\in C\left([0,\infty),H^1\left(\Omega,\bbR^2\right)'\right)\,. 
\end{equation}

Let us now identify the equations solved by the components $f$ and $g$ of $u$. To this end, let~${\chi\in W^1_\infty([0,\infty))}$ be a compactly supported function and $\varphi\in C^1(\overline{\Omega})$. In view of~\eqref{x03a}, classical computations give
\begin{multline*}
\int_0^\infty \int_\Omega \frac{\chi(t+\tau)-\chi(t)}{\tau} f^\tau(t) \varphi\ \rmd x\rmd t + \left(\frac{1}{\tau} \int_0^\tau \chi(t)\ \rmd t \right) \int_\Omega f^{in}\varphi\ \rmd x \\
= \int_0^\infty \int_\Omega \chi(t) f^\tau(t) \nabla\left[af^\tau(t) + b g^\tau(t)\right] \cdot\nabla\varphi\ \rmd x\rmd t\,.
\end{multline*}
Taking $\tau=\tau_m$ in the above identity, it readily follows from~\eqref{x13} and the regularity of $\chi$ and $\varphi$ that we may pass to the limit as $m\to\infty$ and conclude that
\begin{multline}\label{x15}
\int_0^\infty \int_\Omega \frac{d\chi}{dt}(t) f(t,x) \varphi(x)\ \rmd x\rmd t + \chi(0) \int_\Omega f^{in}(x)\varphi(x)\ \rmd x \\
= \int_0^\infty \int_\Omega \chi(t) f(t,x) \nabla\left[af+ b g\right](t,x) \cdot\nabla \varphi(x)\ \rmd x\rmd t\,.
\end{multline}
Since $f\nabla f$ and $f\nabla g$ belong to $L_2((0,T)\times\Omega)$ for all $T>0$ by~\eqref{x13}, a density argument ensures that the identity~\eqref{x15} is valid for any $\varphi\in H^1(\Omega)$. We next use the time continuity~\eqref{x14} of~$f$ and a classical approximation argument to show that~$f$ solves~\eqref{p2a}. A similar argument allows us to derive~\eqref{p2b} from~\eqref{x03b}.

Finally, combining~\eqref{x13}, \eqref{x14}, and a weak lower semicontinuity argument, we may let $m\to\infty$ in~\eqref{x06}, \eqref{x07}, and~\eqref{x08} with $\tau=\tau_m$ to show that $u=(f,g)$ satisfies~\eqref{p3}, \eqref{p4a}, and~\eqref{p5}, thereby completing the proof of Theorem~\ref{ThBWS}.
\end{proof}



\appendix
\section{The polynomials \texorpdfstring{$\Phi_n$}{Phin}, \texorpdfstring{$n\ge 2$}{n>= 2}}\label{secA}



Let $n\ge 2$. According to the discussion in the introduction, we look for an homogeneous polynomial $\Phi_n$ of degree $n$ such that:
\begin{enumerate}\Penumi
\item\label{P1} $\Phi_n$ is convex on $[0,\infty)^2$;
\item\label{P2} the matrix $S_n(X) := D^2\Phi_n(X) M(X)$ is symmetric and positive semidefinite for~${X\in [0,\infty)^2}$.
\end{enumerate}
We recall that the mobility matrix $M(X)$ is given by
\begin{equation*}
M(X) = \left(m_{jk}(X)\right)_{1\,\le\,j\,,\,k\,\le\,2} :=
\begin{pmatrix}
a X_1 & b X_1 \\
c X_2 & d X_2
\end{pmatrix}\,, \qquad X\in\bbR^2\,,
\end{equation*}
see~\eqref{x1}. Specifically, we set
\begin{equation}\label{Polen}
\Phi_n(X) := \sum_{j=0}^n a_{j,n} X_1^j X_2^{n-j}\,, \qquad X=(X_1,X_2)\in \bbR^2\,, 
\end{equation}
with $a_{j,n}$, $0\leq j\leq n$, to be determined in order for properties~\eqref{P1}-\eqref{P2} to be satisfied. We recall that the parameters $(a,\, b,\, c,\, d)$ are assumed to satisfy~\eqref{condabcd}.


\begin{lemm}\label{lem.p1}
Set $a_{0,n}:=1$ and
\begin{equation}\label{ajn}
a_{j,n}:= \prod_{k=0}^{j-1}\frac{(n-k)[ak+c(n-k-1)]}{(k+1)[bk+d(n-k-1)]} = \binom{n}{j} \prod_{k=0}^{j-1}\frac{ak+c(n-k-1)}{bk+d(n-k-1)}
\end{equation}
for $1\leq j\leq n$. Then $a_{j,n}>0$ for $0\le j \le n$ and $S_n(X) = D^2\Phi_n(X) M(X)\in {\bSym_2(\bbR)}$ for~$X\in\bbR^2$.
\end{lemm}


\begin{proof}
Given $X\in\bbR^2$, we compute
\begin{align*}
\partial_1^2 \Phi_n(X) & = \sum_{j=1}^{n-1} j(j+1) a_{j+1,n} X_1^{j-1} X_2^{n-j-1} = \sum_{j=0}^{n-2} (j+1)(j+2) a_{j+2,n} X_1^j X_2^{n-j-2}, \\
\partial_1 \partial_2 \Phi_n(X) & = \sum_{j=1}^{n-1} j(n-j) a_{j,n} X_1^{j-1} X_2^{n-j-1} = \sum_{j=0}^{n-2} (j+1)(n-j-1) a_{j+1,n} X_1^j X_2^{n-j-2}, \\
\partial_2^2 \Phi_n(X) & = \sum_{j=0}^{n-2} (n-j)(n-j-1) a_{j,n} X_1^j X_2^{n-j-2}.
\end{align*}
It then follows that
\begin{align*}
[S_n(X)]_{11} & = a X_1 \partial_1^2 \Phi_n(X) + c X_2 \partial_1 \partial_2 \Phi_n(X)\\
& = a \sum_{j=1}^{n-1} j(j+1) a_{j+1,n} X_1^j X_2^{n-j-1} + c \sum_{j=0}^{n-2} (j+1)(n-j-1) a_{j+1,n} X_1^j X_2^{n-j-1}, \\
[S_n(X)]_{12} & = bX_1 \partial_1^2\Phi_n(X) + d X_2 \partial_1 \partial_2 \Phi_n(X) \\
& = b \sum_{j=1}^{n-1} j(j+1) a_{j+1,n} X_1^j X_2^{n-j-1} + d \sum_{j=0}^{n-2} (j+1)(n-j-1) a_{j+1,n} X_1^j X_2^{n-j-1}\\
& = bn(n-1)a_{n,n}X_1^{n-1} +\sum_{j=1}^{n-2} (j+1)[bj+d(n-j-1)] a_{j+1,n} X_1^j X_2^{n-j-1} \\
& \qquad + d(n-1) a_{1,n} X_2^{n-1},\\
[S_n(X)]_{21} & = aX_1 \partial_1 \partial_2 \Phi_n(X) + c X_2 \partial_2^2 \Phi_n(X)\\
& = a \sum_{j=1}^{n-1} j(n-j) a_{j,n} X_1^{j} X_2^{n-j-1} +c \sum_{j=0}^{n-2} (n-j)(n-j-1) a_{j,n} X_1^j X_2^{n-j-1}\\
& = a (n-1)a_{n-1,n}X_1^{n-1}+\sum_{j=1}^{n-2} (n-j)[aj+c(n-j-1)] a_{j,n} X_1^{j} X_2^{n-j-1} \\
& \qquad + c n (n-1)a_{0,n} X_2^{n-1},\\
[S_n(X)]_{22} & = b X_1 \partial_1 \partial_2 \Phi_n(X) + d X_2 \partial_2^2 \Phi_n(X)\\
& = b \sum_{j=1}^{n-1} j(n-j) a_{j,n} X_1^{j} X_2^{n-j-1} + d \sum_{j=0}^{n-2} (n-j)(n-j-1) a_{j,n} X_1^j X_2^{n-j-1}.
\end{align*}
Hence, $S_n(X)$ is symmetric provided that
\begin{equation*}
(j+1)[bj+d(n-j-1)]a_{j+1,n}= (n-j)[aj+c(n-j-1)] a_{j,n},\quad 0\leq j\leq n-1\,,
\end{equation*}
or, equivalently,
\begin{equation}\label{recursive}
a_{j+1,n} = \frac{(n-j)[aj+c(n-j-1)]}{(j+1)[bj+d(n-j-1)]} a_{j,n}\,,\quad 0\leq j\leq n-1\,.
\end{equation}
Since $a_{0,n}=1$, the closed form formula~\eqref{ajn} readily follows from~\eqref{recursive} and we deduce from~\eqref{ajn} and the positivity of $(a,\,b,\,c,\, d)$ that $a_{j,n}>0$ for all $0\leq j\leq n$.
\end{proof}

We next show that ${D^2\Phi_n(X)}$ is positive definite for $X\in[0,\infty)^2\setminus\{(0,0)\}$. This property implies in particular that ${D^2\Phi_n(X)}$ is positive semidefinite for $X\in[0,\infty)^2$.


\begin{lemm}\label{lem.p2}
Let $\Phi_n$ be the polynomial defined by~\eqref{Polen} and~\eqref{ajn}. Then we have~$D^2\Phi_n(X)\in \bSPD_2(\bbR)$ for~${X\in [0,\infty)^2\setminus\{(0,0)\}}$.
\end{lemm}


\begin{proof}
Given $X\in [0,\infty)^2$, the positivity of the coefficients $a_{j,n}$, $0\leq j\leq n$, of~$\Phi_n$ ensures that
\begin{equation*}
\trr(D^2\Phi_n(X)) := \partial_1^2 \Phi_n(X) + \partial_2^2 \Phi_n(X) \ge 0\,, \qquad X\in [0,\infty)^2\,.
\end{equation*}
It remains to show that the determinant $\det(D^2\Phi_n(X))$ is also non-negative. To this end we compute
\begin{equation}\label{lf07}
\begin{split}
\det(D^2\Phi_n(X)) & = \partial_1^2 \Phi_n(X) \partial_2^2 \Phi_n(X) - [\partial_1 \partial_2 \Phi_n(X)]^2 \\
& = \sum_{j=0}^{n-2} \sum_{k=0}^{n-2} (j+1)(n-k-1) A_{j,k} X_1^{j+k} X_2^{2n-j-k-4}, 
\end{split}
\end{equation}
where
\begin{equation*}
A_{j,k} := (j+2)(n-k) a_{j+2,n} a_{k,n} - (n-j-1)(k+1) a_{j+1,n} a_{k+1,n}\,, \qquad 0 \le j,k \le n-2\,.
\end{equation*}
Using~\eqref{recursive}, we express $a_{j+2,n}$ and $a_{k+1,n}$ in terms of $a_{j+1,n}$ and $a_{k,n}$, respectively, to arrive at the following formula
\begin{multline}\label{Ajk}
A_{j,k}\\
\begin{aligned}
&= (n-k)(n-j-1) \left[\frac{a(j+1)+c(n-j-2)}{b(j+1)+d(n-j-2)} - \frac{ak+c(n-k-1)}{bk+d(n-k-1)}\right] a_{j+1,n} a_{k,n}\\[1ex]
&= (ad-bc)(n-k)(n-j-1) \\[1ex] &\qquad\times\frac{(j+1)(n-k-1)-k(n-j-2)}{\left[b(j+1)+d(n-j-2)\right] \left[bk+d(n-k-1)\right]} a_{j+1,n} a_{k,n} \\[1ex]
&= (ad-bc)\frac{(n-1)(n-k)(n-j-1)(j+1-k)}{\alpha_{j+1,n}\alpha_{k,n}} a_{j+1,n} a_{k,n}\,,
\end{aligned}
\end{multline}
where $\alpha_{k,n}$ denotes the positive number
\begin{equation*}
\alpha_{k,n}:=bk+d(n-k-1)\,,\quad 0\leq k\leq n-1\,.
\end{equation*}
In particular,
\begin{equation}\label{lf09}
A_{k-1,j+1} = - A_{j,k}\,, \qquad 0\le j\le n-3\,, \ 1\le k\le n-2\,. 
\end{equation}
It then follows from~\eqref{lf07} that
\begingroup
\allowdisplaybreaks
\begin{align*}
2 \det(D^2\Phi_n(X))& = \sum_{j=0}^{n-2} \sum_{k=0}^{n-2} (j+1)(n-k-1) A_{j,k} X_1^{j+k} X_2^{2n-j-k-4} \\
& \qquad + \sum_{l=1}^{n-1} \sum_{i=-1}^{n-3} l(n-i-2) A_{l-1,i+1} X_1^{i+l} X_2^{2n-i-l-4} \\
& = \sum_{j=0}^{n-2} \sum_{k=0}^{n-2} (j+1)(n-k-1) A_{j,k} X_1^{j+k} X_2^{2n-j-k-4} \\
& \qquad + \sum_{j=-1}^{n-3} \sum_{k=1}^{n-1} k(n-j-2) A_{k-1,j+1} X_1^{j+k} X_2^{2n-j-k-4}\\
& = \sum_{j=0}^{n-3} \sum_{k=1}^{n-2} (j+1)(n-k-1) A_{j,k} X_1^{j+k} X_2^{2n-j-k-4} \\
& \qquad + \sum_{k=0}^{n-2} (n-1)(n-k-1) A_{n-2,k} X_1^{n-2+k} X_2^{n-k-2} \\
& \qquad + \sum_{j=0}^{n-3} (j+1)(n-1) A_{j,0} X_1^{j} X_2^{2n-j-4} \\
& \qquad + \sum_{j=0}^{n-3} \sum_{k=1}^{n-2} k(n-j-2) A_{k-1,j+1} X_1^{j+k} X_2^{2n-j-k-4} \\
& \qquad + \sum_{k=1}^{n-1} k(n-1) A_{k-1,0} X_1^{k-1} X_2^{2n-k-3} \\
& \qquad + \sum_{j=0}^{n-3} (n-1)(n-j-2) A_{n-2,j+1} X_1^{j+n-1} X_2^{n-j-3}.
\end{align*}
\endgroup
According to~\eqref{condabcd} and~\eqref{Ajk},
\begin{align*}
A_{l,0} & = (ad-bc) \frac{n(n-1)(n-1-l)(l+1)}{\alpha_{0,n}\alpha_{l+1,n}}>0\,, \qquad 0\le l\le n-2\,,\\
A_{n-2,l} & = (ad-bc) \frac{(n-1)(n-l)(n-1-l)}{\alpha_{n-1,n}\alpha_{l,n}}>0\,, \qquad 0\le l\le n-2\,.
\end{align*}
In particular, all the terms in the above identity involving a single sum are non-negative. Therefore, using the symmetry property~\eqref{lf09} and retaining in the last two sums only the terms corresponding to $k=1$ and $j=n-3$, respectively, we get
\begin{align*}
2 \det\left(D^2\Phi_n(X)\right)& \ge \sum_{j=0}^{n-3} \sum_{k=1}^{n-2} \left[(j+1)(n-k-1) - k(n-j-2) \right] A_{j,k} X_1^{j+k} X_2^{2n-j-k-4} \\
& \qquad + (n-1) A_{n-2,n-2} X_1^{2n-4} + (n-1) A_{0,0} X_2^{2n-4}
\end{align*}
\begin{align*}
& = \sum_{j=0}^{n-3} \sum_{k=1}^{n-2} (n-1)(j+1-k) A_{j,k} X_1^{j+k} X_2^{2n-j-k-4} \\
& \qquad + (n-1) A_{n-2,n-2} X_1^{2n-4} + (n-1) A_{0,0} X_2^{2n-4}.
\end{align*}
Observing that
\begin{equation*}
(n-1)(j+1-k) A_{j,k} = (ad-bc) \frac{(n-1)^2(n-k)(n-j-1)(j+1-k)^2}{\alpha_{j+1,n}\alpha_{k,n}} a_{j+1,n} a_{k,n} \ge 0
\end{equation*}
for $0\le j, k\le n-2$, we conclude that
\begin{equation}\label{lf10}
2\det\left(D^2\Phi_n(X)\right) \ge (n-1) A_{n-2,n-2} X_1^{2n-4} + (n-1) A_{0,0} X_2^{2n-4} 
\end{equation}
for $X\in [0,\infty)^2$. Since $A_{0,0}>0$ and $A_{n-2,n-2}>0$, we have thus established that the symmetric matrix $D^2\Phi_n(X)$ has non-negative trace and positive determinant, so that it is positive definite for each $X\in [0,\infty)^2\setminus\{(0,0)\}$.
\end{proof}

We next turn to the positive definiteness of $S_n = D^2\Phi_n M$.


\begin{lemm}\label{lem.p3}
Let $\Phi_n$ be defined by~\eqref{Polen} and~\eqref{ajn}. Then
\[
S_n(X)=D^2\Phi_n(X)M(X)\in \bSPD_2(\bbR)\text{ for }X\in (0,\infty)^2\,.
\]
\end{lemm}


\begin{proof}
Let $X\in (0,\infty)^2$. On the one hand, by~\eqref{condabcd}, \eqref{lf10}, and the positivity of~$A_{0,0}$ and $A_{n-2,n-2}$,
\begin{align*}
2\det(S_n(X)) & = 2(ad-bc)X_1 X_2 \det\left(D^2\Phi_n(X)\right) \\
& \geq(ad-bc)X_1X_2(n-1) \Big[A_{n-2,n-2} X_1^{2n-4} + A_{0,0} X_2^{2n-4} \Big]>0\,.
\end{align*}
On the other hand, the positivity of $a_{j,n}$ for $0\le j \le n$ and~\eqref{condabcd} imply that
\begin{equation*}
\trr(S_n(X)) = [S_n(X)]_{11} + [S_n(X)]_{22}> 0\,.
\end{equation*}
Consequently, $S_n(X)$ has positive trace and positive determinant, and is thus positive definite as claimed.
\end{proof}

We end up this section with useful upper and lower bounds for $\Phi_n$.


\begin{lemm}\label{lelfb}
Let $\Phi_n$ be defined by~\eqref{Polen} and~\eqref{ajn}. Then
\begin{equation}\label{LUB}
\frac{(c X_1 + d X_2)^n}{d^n} \le \Phi_n(X) \le \frac{(a X_1 + b X_2)^n}{b^n}\,, \qquad X\in [0,\infty)^2\,.
\end{equation}
\end{lemm}


\begin{proof}
Since the function
\begin{equation*}
\chi (z) := \frac{ (a-c) z+c}{(b-d) z+d}\,, \qquad z\in [0,1]\,,
\end{equation*}
is increasing and positive, we deduce from~\eqref{ajn} that, for $1\le j \le n$,
\begin{align*}
a_{j,n} &= \binom{n}{j} \prod_{k=0}^{j-1}\chi\left(\frac{k}{n-1} \right)\leq \binom{n}{j} [\chi(1)]^j=\binom{n}{j} \left(\frac{a}{b} \right)^j\\
\intertext{and}
a_{j,n} &= \binom{n}{j} \prod_{k=0}^{j-1}\chi\left(\frac{k}{n-1} \right)\geq \binom{n}{j} [\chi(0)]^j=\binom{n}{j} \left(\frac{c}{d} \right)^j\,.
\end{align*}
The upper and lower bounds in~\eqref{LUB} are direct consequences of the above inequalities.
\end{proof}


\subsection*{Acknowledgments}
We would like to express our thanks to the anonymous referee for the careful reading of the manuscript and valuable comments.

\bibliography{laurencot}
\end{document}
