\documentclass[reqno]{amsart}
\usepackage{hyperref}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2016 (2016), No. 200, pp. 1--25.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu}
\thanks{\copyright 2016 Texas State University.}
\vspace{8mm}}

\begin{document}
\title[\hfilneg EJDE-2016/200\hfil Existence of bounded solutions]
{Existence of bounded solutions for quasilinear parabolic systems
 with quadratic growth}

\author[R. Souilah \hfil EJDE-2016/200\hfilneg]
{Rezak Souilah}

\address{Rezak Souilah \newline
Universit\'e Dr. Yahia Fares M\'ed\'ea
Facul\'e des sciences \'economiques,
commerciales et de gestions,
M\'ed\'ea, Algeria. \newline
- Laboratoire d'EDP Non Lin\'eaires et Histoires des Math\'ematiques,
 ENS, Kouba, Alger, Algeria}
\email{souilah.rezak@yahoo.com}

\thanks{Submitted February 23, 2016. Published July 27, 2016.}
\subjclass[2010]{35K40, 35K59, 35K51}
\keywords{Quasilinear parabolic systems; quadratic growth; 
\hfill\break\indent upper and lower solutions; bounded solutions}

\begin{abstract}
 Assuming the existence of an upper and a lower solution, we prove
 the existence of at least one bounded solution of a quasilinear parabolic
 systems, with nonlinear second member having a quadratic growth with
 respect to the gradient of the solution.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{remark}[theorem]{Remark}
\newtheorem{definition}[theorem]{Definition}
\allowdisplaybreaks

 \section{Introduction}

 Let $\Omega$ be a bounded open subset of $\mathbb{R}^N$, with boundary
$\partial\Omega$ and let $Q$ be the cylinder $\Omega \times (0,T)$ with
 some given $T>0$. Consider the quasilinear parabolic system
\begin{equation}\label{system}
 \begin{gathered}
 \frac{\partial u^1}{\partial t}
-\operatorname{div}(A(u)\nabla u^1)
=G^1(u,\nabla u)+  F(u,\nabla u)\cdot\nabla u^1\quad\text{in }Q, \\
 \frac{\partial u^2}{\partial t}
-\operatorname{div}(A(u)\nabla u^2)=G^2(u,\nabla u)+
 F(u,\nabla u)\cdot\nabla u^2\quad\text{in }Q, \\
 u^1(x,t)= u^2(x,t)=0\quad\text{on}\quad\Sigma=\partial\Omega \times (0,T), \\
 u^1(x,0)=u^1_0(x),\quad u^2(x,0)=u^2_0(x)\quad\text{in }\Omega.
  \end{gathered}
\end{equation}
 where $u_0^1(x), u_0^2(x)\in L^{\infty}(\Omega) $ and
 $$
\operatorname{div}(A(u)\nabla u^{\gamma})
=\sum_{i,j=1}^{N}\frac{\partial}{\partial
 x_i}\Big( A_{i,j}(u)\frac{\partial u^{\gamma}}{\partial x_j}\Big), \quad
\gamma =1,2
$$
with  $A_{i,j}:Q\times\mathbb{R}^2\to \mathbb{R} $ are
Carath\'eodory functions which satisfy the following assumptions
 \begin{equation}\label{ellipticite}
 \quad\exists\alpha>0,\;\forall\xi\in\mathbb{R}^{N},\;\forall s\in \mathbb{R}^2,
 \sum_{i,j=1}^{N} A_{i,j}(x,t,s)\xi_i\xi_j\geq\alpha|\xi|^2\;\;\text{a.e. }(x,t)\in Q;
 \end{equation}
 and there exists $\varrho>0$ such that for all $s\in \mathbb{R}^2$,
 \begin{equation}\label{uniformborne}
 |A_{i,j}(x,t,s)|\leq\varrho\quad\text{a.e. }(x,t)\in Q.
 \end{equation}
 The functions $G^1, G^2:Q\times\mathbb{R}^2\times\mathbb{R}^{2N}
\to \mathbb{R}$ are Carath\'{e}odory functions which satisfy the quadratic
growth assumptions:
for all $s\in\mathbb{R}^2$ and all $\xi=(\xi^1, \xi^2)\in{\mathbb{R}^{2N}}$,
 \begin{gather}\label{G1croissance}
|G^1(x,t,s,\xi)|\leq  C_0+C_2|\xi^1|^2+\eta|\xi^2|^2\quad
\text{a.e. }(x,t)\in Q; \\
\label{G2croissance}
|G^2(x,t,s,\xi)|\leq  C_0+C_2[|\xi^1|^2+|\xi^2|^2]\quad
\text{a.e. }(x,t)\in Q.
 \end{gather}
$F:Q\times\mathbb{R}^2\times\mathbb{R}^{2N}\to \mathbb{R}^{N}$ is a
Carath\'{e}odory function satisfying the sublinearity assumption:
for all $s\in\mathbb{R}^2$ and all $\xi=(\xi^1, \xi^2)\in{\mathbb{R}^{2N}}$,
 \begin{equation}\label{Fcroissance}
 |F(x,t,s,\xi)|\leq  C_3+C_4|\xi|\quad\text{a.e. }(x,t)\in Q.
 \end{equation}
where $C_0,C_2,C_3,C_4$ and $\eta$ are positive constants,
 $\eta$ being small enough.
 Several papers concern mainly the regularity properties of
solution of elliptic and parabolic system;
 see e.g. [\cite{DK}, \cite{Gia}--\cite{KIKo}, \cite{Per}, \cite{Struwe},
\cite{Wid}--\cite{Wie2}].
In the elliptic case with quadratic growth, in \cite{Sa}, the author studies
a unilalteral problem for $L^1$-data, in which the truncate function is used
instead of upper and lower solutions.
In \cite{Ben1}--\cite{Ben3}, the authors study renormalized or entropic
parabolic systems to overcame the lake of regularities of solutions.

Others articles study the existence and regularity of a solution using as a
main tool some regularity arguments and strong maximum principles see e.g.
\cite{JF1}--\cite{JF8}. Others extend the weak maximum principles to a special
class of systems of parabolic equations, the so-called weakly coupled
systems(it is coupled only through the terms which are not differentiated,
each equation containing derivatives of just one component)
see e.g. \cite{Dickstein,Sleeman,Portnyagin}.
Others articles study the existence of a solution using monotony arguments
(the method of upper and lower solutions) see e.g.
\cite{Aris,Cosner,Leung,Liang,Mckenna,Pao2,Sindler} and the book
\cite{Leung2} and the references therein. In \cite{Abudiab,Abudiab2,Liang}
 the authors have extended the method of classical upper-lower methods
for elliptic and parabolic systems without the assumption of quasi-monotonicity.
The Growing conditions \eqref{G1croissance} and \eqref{G2croissance}
imposed on $G^1$ and $G^2$ and the growth condition \eqref{Fcroissance}
imposed on $F$ are sufficient to have an uniform estimate of $u$ in the
space $\big( L^2(0, T, H^1_0(\Omega))\big) ^2$.
Note that if the condition \eqref{G1croissance} is the same as
\eqref{G2croissance} we need to add a condition of type
$C_2\|u\|_{\infty}<\alpha$ to have a uniform estimate of $u$ in the space
$\big( L^2(0, T, H^1_0(\Omega))\big) ^2$ and also an uniform estimation
in space $\mathcal{C}^{\delta}$ see e.g. \cite{Hil5} for more detailed
on this subject.


When $\eta=0$, the system has a triangular structure with respect to the
quadratic terms, and the system can be decoupled.
When $\eta$ is positive and small, in \cite{AM2} the author establishes
the existence of solutions of the associated elliptic systems.

 \section{Statement of the main result}

\begin{theorem}\label{theor}
  Under hypotheses \eqref{ellipticite}--\eqref{Fcroissance} and the smallness
condition \eqref{Eta} for $\eta$, there exists at least one solution $u$ of
system \eqref{system}.
\end{theorem}

A solution to  \eqref{system} must be interpreted in the weak sense:
$$
u\in \big( L^2(0,T;H_0^1(\Omega))\cap L^{\infty}(Q)\big)^2,\quad
 \frac{\partial u}{\partial t}\in
\big( L^2(0,T;H^{-1}(\Omega))+ L^1(Q)\big)^2,
$$
such that for all $v\in( L^2(0,T;H_0^1(\Omega))\cap L^{\infty}(Q)^2$,
$\frac{\partial v}{\partial t}=\beta_1+\beta_2$ where
$\beta_1\in \big( L^2(0,T;H^{-1}(\Omega))\big)^2$ and
$\beta_2\in \big(L^1(Q)\big)^2$,
\begin{align*}
&\int_{\Omega}u^{\gamma}(T)v^{\gamma}(T)dx
 -\int_{\Omega}u^{\gamma}(0)v^{\gamma}(0)dx
 -\int_0^{T} \langle\,\beta_1^{\gamma}, u^\gamma
 \rangle _{H^{-1}(\Omega)\times H^1_0(\Omega)}dt \\
& -\int_{Q}\beta_2^{\gamma}u^{\gamma}dx\,dt
 +\int_{Q}A(u)\nabla u^{\gamma}\cdot\nabla v^\gamma dx\,dt\\
&=\int_{Q}G^{\gamma}(u,\nabla u)v^\gamma dx\,dt
 +\int_{Q}F(u,\nabla u)\cdot\nabla u^\gamma v^\gamma dx\,dt,\quad \text{for }
\gamma=1,2
\end{align*}


The proof of Theorem \ref{theor}] will be performed in three steps:
firstly prove that the approximated system of \eqref{system} admits at least
one bounded solution, denoted $u_{\varepsilon }$;
 secondly we prove an $(L^2(0,T; H^1_0(\Omega)))^2$-estimate for $u_{\varepsilon }$,
then the strong convergence in $(L^2(0,T; H^1_0(\Omega)))^2$ of $u_{\varepsilon }$;
and finally we pass to the limit in the approximated system of \eqref{system}.

 \section{Approximation}

 According to \cite{BM,AM2,Sa}, we regularize the nonlinear terms to be bounded,
 for that we consider now the approximated system of \eqref{system}:
 \begin{equation}\label{asystem}
 \begin{gathered}
  \frac{\partial u_{\varepsilon }^1}{\partial t}
-\operatorname{div}(A(u_{\varepsilon })\nabla u_{\varepsilon }^1)
=G_{\varepsilon }^1(u_{\varepsilon },\nabla u_{\varepsilon })
+F_{\varepsilon }(u_{\varepsilon },\nabla u_{\varepsilon })
\cdot\nabla u_{\varepsilon }^1\quad\text{in }Q,
 \\
 \frac{\partial u_{\varepsilon }^2}{\partial t}-\operatorname{div}
(A(u_{\varepsilon })\nabla u_{\varepsilon }^2)
=G_{\varepsilon }^2(u_{\varepsilon },\nabla u_{\varepsilon })
+F_{\varepsilon }(u_{\varepsilon },\nabla u_{\varepsilon })
\cdot\nabla u_{\varepsilon }^2\quad\text{in }Q,
 \\
 u_{\varepsilon }^1(x,t)=0,\; u_{\varepsilon}^2(x,t)=0\quad\text{on }\Sigma,
 \\
 u_{\varepsilon }^1(x,0)=u^1_0(x),\;u_{\varepsilon }^2(x,0)=u^2_0(x)
\quad\text{in }\Omega,
 \end{gathered}
 \end{equation}
 where $\varepsilon >0$, and
$G_{\varepsilon }^1(x,t,s,\xi),\;G_{\varepsilon }^2(x,t,s,\xi):
 Q\times\mathbb{R}^2\times\mathbb{R}^{2  N}\to \mathbb{R}$
 and $F_{\varepsilon }(x,s,\xi):Q\times\mathbb{R}^2\times\mathbb{R}^{2N}\to
\mathbb{R}^{N}$ are Carath\'{e}odory functions such  that:
 \begin{equation} \label{G1G2Fdef}
\begin{gathered}
 G_{\varepsilon }^1(x,t,s,\xi)=\frac{G^1(x,t,s,\xi)}{1+\varepsilon |G^1(x,t,s,\xi)|},
 \\
 G_{\varepsilon }^2(x,t,s,\xi)=\frac{G^2(x,t,s,\xi)}{1+\varepsilon |G^2(x,t,s,\xi)|},
 \\
 F_{\varepsilon }(x,t,s,\xi)=\frac{F(x,t,s,\xi)}{1+\varepsilon |
 F(x,t,s,\xi)||\xi|}.
 \end{gathered}
\end{equation}
 Noting that the functions $G_\varepsilon ^1$, $G_\varepsilon ^2$  and
$F_\varepsilon $  satisfy the following conditions:
for all $s\in\mathbb{R}^2$, all $\xi=(\xi^1,\xi^2)\in\mathbb{R}^{2N}$,
and for a.e. $(x,t)\in Q$, we have
 \begin{gather}\label{G1G2Fineq}
|G_{\varepsilon }^1(x,t,s,\xi)|\leq\frac{1}{\varepsilon },\quad
|G_{\varepsilon }^2(x,t,s,\xi)|\leq\frac{1}{\varepsilon },\quad
|F_{\varepsilon }(x,t,s,\xi)\xi^{\gamma}|\leq\frac{1}{\varepsilon}, \\
\label{G1G2Fcroissance}
|G_{\varepsilon }^1(x,t,s,\xi)|\leq|G^1(x,t,s,\xi)|,\quad
|G_{\varepsilon }^2(x,t,s,\xi)|\leq|G^2(x,t,s,\xi)|, \\
|F_{\varepsilon }(x,t,s,\xi)|\leq|F(x,t,s,\xi)|.
 \end{gather}

 Since the right hand side of each equation in \eqref{asystem} is
 bounded by $\frac{2}{\varepsilon}$, then by applying the De Giorgi
 iteration technique \cite[theorem 4.2.1]{JZ},   for  each $\gamma=1,2$,
we have
 \begin{equation}\label{H}
 \sup_{Q}u_{\varepsilon }^{\gamma}\leq\sup_{\partial Q}
 u_{\varepsilon }^{\gamma}+C\|H_{\varepsilon }^{\gamma}\|_{L^{\infty}(Q)}
 \end{equation}
 where $H_{\varepsilon }^\gamma=G_{\varepsilon }^{\gamma}(u_{\varepsilon },\nabla
 u_{\varepsilon })+ F_{\varepsilon }(u_{\varepsilon },\nabla
 u_{\varepsilon })\cdot\nabla u_{\varepsilon }^{\gamma}$ and $C$ is a constant
depending   only on $N$ and $\Omega$.
Hence
 \begin{equation}\label{inff}
 \|u_{\varepsilon }^\gamma\|_{L^{\infty}(Q)}\leq\|u_0^\gamma\|_{L^{\infty}(Q)}+
 \frac{2C}{\varepsilon }=C_{\varepsilon}.
 \end{equation}
 Unfortunately the estimate \eqref{inff} is not sufficient to obtain a
uniform estimate of a possible solution of the regularization
system \eqref{asystem} in the space $\big( L^2(0,T;H_0^1(\Omega))\big)^2$
(see the proof of Lemma \ref{lemma3.1}). In fact, we will need the uniform estimate
 \begin{equation*}
 \|u_{\varepsilon }^\gamma\|_{L^{\infty}(Q)}\leq M
 \quad\text{for } \gamma=1,2,
 \end{equation*}
 where $M$ is positive constant independent of $\varepsilon$, this is the
 main goal of the next step. We will define the upper and lower solution
of regularization system \eqref{asystem} and we will consider an auxiliary
 modified system whose solution is between the upper and lower solutions
and satisfy the system \eqref{asystem}. First, we need to state some notations.
For given $v=(v^1,v^2)$ and $w=(w^1,w^2)$, we say
 $v\le w$ if $v^1\le w^1$ and $v^2\le w^2$. We consider also this notation
$[v]_w^1=(w^1,v^2)$ and  $[v]_w^2=(v^1,w^2)$.

 \begin{definition} \label{def3.1} \rm
 Let $\varphi$ and $\phi\in\left( L^{\infty}(0,T;W^{1,\infty}(\Omega))\right)^2$
such that  $\frac{\partial\varphi}{\partial t}$ and
$\frac{\partial\phi}{\partial t}\in  \left( L^2(0,T;H^{-1}(\Omega))\right)^2$.
Then $\varphi$ and $\phi$ are called ordered coupling weak upper and lower
solution of \eqref{asystem}, if $\phi\le\varphi$ and for all
$ v\in\big( L^2(0,T;H_0^1(\Omega))\big)^2$ such that
 $ \phi\le v\le \varphi$, a.e. in $Q$, they satisfy
 \begin{equation}\label{upper}
 \begin{gathered}
 \frac{\partial \varphi^1}{\partial t}
-\operatorname{div}(A([v]_{\varphi}^1)\nabla\varphi^1)
\ge  G_{\varepsilon }^1 ( [v]_{\varphi}^1,\nabla[v]_{\varphi}^1)
+  F_{\varepsilon }([v]_{\varphi}^1,\nabla[v]_{\varphi}^1)\cdot\nabla\varphi^1
\quad\text{in }Q,
 \\
 \frac{\partial \varphi^2}{\partial t}
-\operatorname{div}(A([v]_{\varphi}^1)\nabla\varphi^2)
\ge  G_{\varepsilon }^2\big( [v]_{\varphi}^2,\nabla[v]_{\varphi}^2\big)
+  F_{\varepsilon }([v]_{\varphi}^2,\nabla[v]_{\varphi}^2)\cdot
\nabla\varphi^2\quad \text{in }Q,
 \\
 \varphi\ge 0,\quad\text{on }\Sigma, \\
 \varphi(0)\ge u_0,\quad\text{in }\Omega,
 \end{gathered}
 \end{equation}
 and
 \begin{equation}\label{lower}
 \begin{gathered}
 \frac{\partial \phi^1}{\partial t}
-\operatorname{div}(A([v]_{\phi}^1)\nabla\phi^1)
\le  G_{\varepsilon }^1([v]_{\phi}^1,\nabla[v]_{\phi}^1)+
 F_{\varepsilon }([v]_{\phi}^1,\nabla[v]_{\phi}^1)\cdot\nabla\phi^1\quad\text{in }Q,
 \\
 \frac{\partial \phi^2}{\partial t}
-\operatorname{div}(A([v]_{\phi}^2)\nabla\phi^2)
\le  G_{\varepsilon }^2([v]_{\phi}^2,\nabla[v]_{\phi}^2)+
 F_{\varepsilon }([v]_{\phi}^2,\nabla[v]_{\phi}^2)\cdot\nabla\phi^2\quad\text{in }
Q,  \\
 \phi\le 0,\quad\text{on }\Sigma, \\
 \phi(0)\le u_0,\quad\text{in }\Omega.
 \end{gathered}
 \end{equation}
\end{definition}

\begin{remark} \label{rmk1} \rm
For applications to stochastic differential games (see for example, the
\cite[condition (2.7)]{JF8}) we can assume that $G^1$ and $G^2$ satisfy
the following conditions:
There exists $K>0$ such that for all $s\in\mathbb{R}^2$ and all
$\xi=(\xi^1, \xi^2)\in{\mathbb{R}^{2N}}$, we have
 \begin{equation}\label{K}
 |G^1(x,t,s,\xi)|_{\xi^1=0}\leq K,\quad
|G^2(x,t,s,\xi)|_{\xi^2=0}\leq K \quad\text{a.e. }(x,t)\in Q.
\end{equation}
According to this  conditions, it is easy to see that we can take as
upper and lower solution the constant functions with respect to the
space variable $\varphi=(\varphi^1,\varphi^2)
=(Kt+\|u_0^1\|_{L^{\infty}(\Omega)},Kt+\|u_0^2\|_{L^{\infty}(\Omega)})$
and $\phi=(\phi^1,\phi^2)=(-Kt-\|u_0^1\|_{L^{\infty}(\Omega)},
-Kt-\|u_0^2\|_{L^{\infty}(\Omega)})$.
\end{remark}

Let us extend the solution of problem  \eqref{asystem} by continuity as
 \begin{equation}\label{troncature}
 \bar{u}_{\varepsilon }=(\bar{u_{\varepsilon }}^1,\bar{u_{\varepsilon }}^2),\quad
 \bar{u_{\varepsilon }}^{\gamma} =
 \begin{cases}
 \varphi^{\gamma} & \text{if $u_{\varepsilon }^{\gamma}\ge \varphi^{\gamma}$},\\
 u_{\varepsilon }^{\gamma} & \text{if $\phi^{\gamma}\leq u_{\varepsilon }^{\gamma}\le \varphi^{\gamma}$},\\
 \phi^{\gamma} & \text{if $u_{\varepsilon }^{\gamma}\le \phi^{\gamma}$}.
 \end{cases}
 \quad \gamma=1,2.
 \end{equation}

 Now we consider the  approximated and the extended system
 \begin{equation}\label{modification}
 \begin{gathered}
  \frac{\partial u_{\varepsilon }^1}{\partial t}
-\operatorname{div}( \widehat{A}(u_{\varepsilon })\nabla u_{\varepsilon }^1)
=\widehat{G}_{\varepsilon }^1(u_{\varepsilon },\nabla u_{\varepsilon })
+ \widehat{F}_{\varepsilon }(u_{\varepsilon },\nabla u_{\varepsilon })
\cdot\nabla\bar{u_{\varepsilon }}^1\quad\text{in }Q, \\
 \frac{\partial u_{\varepsilon }^2}{\partial t}
-\operatorname{div}( \widehat{A}(u_{\varepsilon })\nabla u_{\varepsilon }^2)
=\widehat{G}_{\varepsilon }^2(u_{\varepsilon },\nabla u_{\varepsilon })
+ \widehat{F}_{\varepsilon }(u_{\varepsilon },\nabla u_{\varepsilon })
 \cdot\nabla\bar{u_{\varepsilon }}^2\quad\text{in }Q, \\
 u_{\varepsilon }^1(x,t)=0,\; u_{\varepsilon}^2(x,t)=0\quad\text{on }\Sigma, \\
 u_{\varepsilon }^1(x,0)=u^1_0(x),\;u_{\varepsilon}^2(x,0)=u^2_0(x)\quad
\text{in }\Omega,
 \end{gathered}
 \end{equation}
where
 \begin{equation}\label{hate}
 \widehat{G}_{\varepsilon }^{\gamma}(u_{\varepsilon },\nabla u_{\varepsilon })
= G_{\varepsilon }^{\gamma}(\bar{u_{\varepsilon }},\nabla \bar{u_{\varepsilon }}),\quad
 \widehat{F}_{\varepsilon }(u_{\varepsilon },\nabla u_{\varepsilon })
= F_{\varepsilon }(\bar{u_{\varepsilon }},\nabla \bar{u_{\varepsilon }}),\quad
\widehat{A}(u_{\varepsilon })=A(\bar{u}_{\varepsilon}).
 \end{equation}

\subsection{Existence of solutions of system (3.16)}

 \begin{theorem}\label{theor1}
 If there exist an upper and  lower solutions $ \varphi$ and $ \phi$ of
 system \eqref{asystem}, then for all $\varepsilon>0$, there exists at least one
solution $u_{\varepsilon }$ of system \eqref{modification} which satisfies
\begin{gather*}
u_{\varepsilon }\in \left( L^2(0,T;H_0^1(\Omega))\right)^2,\quad
\frac{\partial u_{\varepsilon }}{\partial t}\in
 \left( L^2(0,T;H^{-1}(\Omega))\right)^2, \\
\phi\leq u_{\varepsilon }\leq \varphi,\quad\text{a.e. in }Q.
\end{gather*}
\end{theorem}



\begin{proof}[Proof of Theorem \ref{theor1}]
In view of \eqref{G1G2Fineq} and \eqref{hate}, an application of Schauder's
fixed point theorem implies that  system \ref{modification} has at least one
solution for $\varepsilon >0$ given. Let now $u_{\varepsilon }$ be a solution
of system \eqref{modification} and let us show that
$ u^{\gamma}_{\varepsilon }\leq\varphi^{\gamma}$ a.e. in $Q$ for  all $\gamma=1,2$.
Using \eqref{upper} and \eqref{modification} (with $v=\bar{u_{\varepsilon }}$)
for  all $\gamma=1,2$ we obtain
 \begin{equation}\label{ineqq}
 \begin{gathered}
  \frac{\partial (u_{\varepsilon }^{\gamma}-\varphi^{\gamma})}{\partial t}
-\operatorname{div}\left(\widehat{A}(u_{\varepsilon })
\nabla(u_{\varepsilon }^{\gamma}-\varphi^{\gamma})\right)
-\operatorname{div}\left(\left[\widehat{A}(u_{\varepsilon })
-A([\bar{u_{\varepsilon }}]_{\varphi}^{\gamma})\right]\nabla\varphi^{\gamma}\right)
 \\
 +\big[G_{\varepsilon }^{\gamma}([\bar{u_{\varepsilon }}
]_{\varphi}^{\gamma},\nabla[\bar{u_{\varepsilon }}]_{\varphi}^{\gamma})
-\widehat{G}_{\varepsilon }^{\gamma}(u_{\varepsilon },\nabla
 u_{\varepsilon })\big]
 \\
 +\big[F_{\varepsilon }([\bar{u_{\varepsilon }}]_{\varphi}^{\gamma},\nabla[\bar{u_{\varepsilon }}]_{\varphi}^{\gamma})\cdot\nabla\varphi^{\gamma}-\widehat{F}_{\varepsilon }(u_{\varepsilon },\nabla
 u_{\varepsilon })\cdot\nabla\bar{u_{\varepsilon }}^{\gamma}\big]\leq 0
\quad\text{in }Q,
 \\
 u_{\varepsilon}^{\gamma}-\varphi^{\gamma}\leq 0,\quad\text{on }\Sigma,
 \\
 (u_{\varepsilon}^{\gamma}-\varphi^{\gamma})(x,0)\leq 0,\quad\text{in }\Omega.
\end{gathered}
 \end{equation}
 We multiply  by $(u_{\varepsilon }^{\gamma}-\varphi^{\gamma})^{+}$,
using \cite[lemma 2.4]{BM} we obtain
 \begin{equation}\label{ineq2}
 \begin{aligned}
&\frac{1}{2}\| (u_{\varepsilon }^{\gamma}-\varphi^{\gamma})^{+}(T)
 \|^2_{L^2(\Omega)}+\int_{Q} \widehat{A}(u_{\varepsilon })
 \nabla(u_{\varepsilon }^{\gamma}
-\varphi^{\gamma})\cdot\nabla(u_{\varepsilon }^{\gamma}-\varphi^{\gamma})^{+}dx\,dt
 \\
 &+\int_{Q} \left[\widehat{A}(u_{\varepsilon })-A([\bar{u_{\varepsilon }}]_{\varphi}^{\gamma})\right]\nabla\varphi^{\gamma}\cdot\nabla(u_{\varepsilon }^{\gamma}-\varphi^{\gamma})^{+}dx\,dt
 \\
&+\int_{Q}\left[G_{\varepsilon }^{\gamma}([\bar{u_{\varepsilon }}]_{\varphi}^{\gamma},\nabla[\bar{u_{\varepsilon }}]_{\varphi}^{\gamma})-{\widehat{G}_{\varepsilon }^{\gamma}(u_{\varepsilon },\nabla
 u_{\varepsilon })}\right](u_{\varepsilon }^{\gamma}-\varphi^{\gamma})^{+}dx\,dt\\
&+\int_{Q} \big[F_{\varepsilon }([\bar{u_{\varepsilon }}]_{\varphi}^{\gamma},
 \nabla[\bar{u_{\varepsilon }}]_{\varphi}^{\gamma})\cdot\nabla\varphi^{\gamma}
 -\widehat{F}_{\varepsilon }(u_{\varepsilon },\nabla  u_{\varepsilon })
\cdot\nabla\bar{u_{\varepsilon }}^{\gamma}\big]
(u_{\varepsilon }^{\gamma}-\varphi^{\gamma})^{+} dx\,dt\leq 0.
 \end{aligned}
 \end{equation}
 At the points where $(u_{\varepsilon }^{\gamma}-\varphi^{\gamma})^{+}$ is not zero,
we have in particular $u_{\varepsilon }^{\gamma}\geq\varphi^{\gamma} $,
then  $\widehat{G}_{\varepsilon }^{\gamma}(u_{\varepsilon },\nabla
 u_{\varepsilon })
=G_{\varepsilon }^{\gamma}([\bar{u_{\varepsilon }}]_{\varphi}^{\gamma},
\nabla[\bar{u_{\varepsilon }}]_{\varphi}^{\gamma})$ and
$F_{\varepsilon }([\bar{u_{\varepsilon }}]_{\varphi}^{\gamma},
\nabla[\bar{u_{\varepsilon }}]_{\varphi}^{\gamma})\cdot\nabla\varphi^{\gamma}
=\widehat{F}_{\varepsilon }(u_{\varepsilon },\nabla
 u_{\varepsilon })\cdot\nabla\bar{u_{\varepsilon }}^{\gamma}$ and
$\widehat{A}(u_{\varepsilon })=A([\bar{u_{\varepsilon }}]_{\varphi}^{\gamma})$.
On the other hand, $\nabla(u_{\varepsilon }^{\gamma}-\varphi^{\gamma})^{+}=0$,
a.e. on the set where $(u_{\varepsilon }^{\gamma}-\varphi^{\gamma})^{+}=0$ and
 \begin{equation}\label{ke}
 \nabla(u_{\varepsilon }^{\gamma}-\varphi^{\gamma})\cdot
\nabla(u_{\varepsilon }^{\gamma}-\varphi^{\gamma})^{+}
=\nabla(u_{\varepsilon }^{\gamma}-\varphi^{\gamma})^{+}\cdot
\nabla(u_{\varepsilon }^{\gamma}-\varphi^{\gamma})^{+}\quad\text{a.e. in } Q.
 \end{equation}
Using \eqref{ellipticite}, \eqref{ineq2} and \eqref{ke} we obtain
 \begin{equation}
 \alpha \|(u_{\varepsilon }^{\gamma}-\varphi^{\gamma})^{+}
\|^2_{L^2(0,T;H_0^1(\Omega))}\leq 0,
 \end{equation}
and so $(u_{\varepsilon }^{\gamma}-\varphi^{\gamma})^{+}=0$
 hence
 $ u_{\varepsilon }^{\gamma}\leq\varphi^{\gamma}$ a.e. in $Q$.
In the same way  we can show that $u_{\varepsilon }^{\gamma}\geq\phi^{\gamma}$
a.e. in $Q$ for  $\gamma=1,2$.
Then, we have
 \begin{equation}\label{phipsi}
 \phi\leq u_{\varepsilon }\leq \varphi\quad\text{a.e. in } Q.
 \end{equation}
Finally, we have prove that \eqref{modification} admits at least one solution
$u_{\varepsilon }$ satisfy \eqref{phipsi}. By the definition of
$\widehat{G}_{\varepsilon }^{\gamma}(u_{\varepsilon },\nabla
 u_{\varepsilon })$, $\widehat{F}_{\varepsilon }(u_{\varepsilon },\nabla
 u_{\varepsilon })$ and $\bar{u_{\varepsilon }}^{\gamma}$ for $\gamma=1,2$,
it is clear that $u_{\varepsilon }$ is solution of \eqref{asystem}.
Using \eqref{phipsi} we obtain
 \begin{equation}\label{Maj}
 \|u_{\varepsilon }^{\gamma}\|_{L^{\infty}(Q)}\leq M,\quad\text{for  }\gamma=1,2.
 \end{equation}
where $M$ is positive constant independent of $\varepsilon $.
 We are now able to specify the smallness of the constant $\eta$ which appears
in the growth condition \eqref{G1croissance}: we will assume that
 \begin{equation}\label{Eta}
 0\leq\eta\leq\frac{C_2}{4\exp\big(\frac{64C_2M}{\alpha}\big)}.
 \end{equation}
\end{proof}

 \section{$(L^2(0,T; H^1_0(\Omega)))^2$-estimate}

 We have show that the regularized system \eqref{asystem} admits
at least one solution $u_{\varepsilon }$ for all $\varepsilon >0$.
In the next step we will establish the sufficient conditions which
allow us to pass to the limit in system \eqref{asystem} to obtain a
solution of \eqref{system}. In this step we will need some needful lemmas.
 Firstly we consider the functions:
 $\varphi: \mathbb{R}\to  \mathbb{R}$ defined by
 \begin{equation}\label{varphii}
 \varphi(\tau)=\exp({\lambda \tau})+\exp({-\lambda \tau})-2,\quad
\forall\tau\in\mathbb{R}
 \end{equation}
 and $\psi: \mathbb{R}^2\to  \mathbb{R}$
 defined by
 \begin{equation}\label{psi}
 \psi(s^1,s^2)=\beta_1\varphi(s^1)+\beta_2\varphi(s^2)
 \quad\forall s=(s^1,s^2)\in\mathbb{R}^2,
 \end{equation}
 and $\lambda$, $\mu$ and $\beta$ are positive constants that we choose as
 \begin{equation}\label{data}
  \lambda=\frac{2C_2}{\alpha},\quad \beta_2=\frac{\beta_1}{2\exp(\lambda M)},\quad
\mu=\frac{C_3^2}{2\theta\alpha}  +\frac{C_4^2}{2\theta\alpha}
\end{equation}
with $\theta$ being a fixed number such that
$ 0<\theta\leq\frac{C_2^2\beta_1}{2\alpha\exp(\lambda M)}$.

 \begin{lemma}\label{lemma3.1}
 For any $u_{\varepsilon }$, such that $|u_{\varepsilon }^{\delta}|\leq M$
for any $\delta=1,2$, we have
 \begin{gather}\label{rrho1}
 \rho_1=\alpha\beta_1\varphi''(u_{\varepsilon }^1)
-C_2\beta_1|\varphi'(u_{\varepsilon }^1)|
-C_2\beta_2|\varphi'(u_{\varepsilon }^2)|-\frac{\theta}{2}\geq\alpha_0,\\
\label{rho2}
 \rho_2= \alpha\beta_2\varphi''(u_{\varepsilon }^2)
-C_2\beta_2|\varphi'(u_{\varepsilon }^2)|
 -\eta\beta_1|\varphi'(u_{\varepsilon }^1)|-\frac{\theta}{2}\geq\alpha_0
 \end{gather}
where $\alpha_0$ is defined by
 \begin{equation}\label{alpha 0}
 \alpha_0=\frac{C_2^2\beta_1}{4\alpha\exp(\lambda M)}.
 \end{equation}
 \end{lemma}


\begin{proof}
 Since
 \begin{gather*}
 \forall\tau, \; |\tau|\leq M ,|\varphi'(\tau)|\leq\lambda\exp(\lambda |\tau|)
 \leq\lambda\exp(\lambda M), \\
 \varphi''(\tau)=\lambda^2\left( \exp(\lambda \tau)+\exp(-\lambda \tau)\right)
\geq\lambda^2\exp(\lambda |\tau|),
 \end{gather*}
we obtain that, for any $u_{\varepsilon }^{\delta}$ with
 $|u_{\varepsilon }^{\delta}|\leq M$, $\delta=1,2$,
 \begin{equation}\label{rho1}
 \rho_1\geq \alpha\beta_1 \lambda^2\exp(\lambda|u_\varepsilon ^1|)-C_2\beta_1
\lambda\exp(\lambda|u_\varepsilon ^1|)
 - C_2\beta_2 \lambda\exp(\lambda M)-\frac{\theta}{2}
 \end{equation}
 Since, from \eqref{data}, we have $\alpha\lambda>C_2$, the infinimum of
the right hand side of \eqref{rrho1} is achieved for $|u_{\varepsilon }^1|=0$.
 We estimate from below the right hand side of \eqref{rho1} by
 \begin{equation}\label{rho111}
 C_2\beta_1\lambda-C_2\beta_2 \lambda\exp(\lambda M)-\frac{\theta}{2}.
 \end{equation}
 In view of the values of $\lambda,\beta_2$ and $\theta$ given by \eqref{data},
the right hand side of  \eqref{rho111} is greater than
 \begin{equation}
 \frac{C_2^2\beta_1}{\alpha}-\frac{C_2^2\beta_1}{4\alpha\exp(\lambda M)}
\geq\frac{C_2^2\beta_1}{4\alpha\exp(\lambda M)}=\alpha_0.
 \end{equation}
Inequality \eqref{rho2} can be proved by same way.
\end{proof}


 \begin{lemma}\label{lemma3.2}
 if $v\in L^2(0,T;H_0^1(\Omega))\cap L^{\infty}(Q)$ and
$\frac{\partial v}{\partial t}\in  L^2(0,T;H^{-1}(\Omega))$,
then there exists a sequence $w_j$ such that
 $w_j\in L^2(0,T;H_0^1(\Omega))$, $\frac{\partial  w_j}{\partial t}\in
 L^2(0,T;H^1(\Omega))$, $w_j$ bounded in $L^{\infty}(Q)$
and
 \begin{gather}
 w_j\to  v\quad \text{strongly in } L^2(0,T;H_0^1(\Omega)), \\
  \frac{\partial w_j}{\partial t}\to \frac{\partial v}{\partial t}\quad
\text{strongly in }L^2(0,T;H^{-1}(\Omega)), \\
  w_j(0)\to  v(0)= \quad \text{strongly in }L^2(\Omega).
 \end{gather}
 \end{lemma}

The proof of the above Lemma is given by \cite[lemma 2.2]{BM}.

 \begin{proposition}\label{proposition1}
 Assume that \eqref{ellipticite}--\eqref{Fcroissance} and
 \eqref{G1G2Fcroissance} hold. If the solutions $u_{\varepsilon }$ of the
 approximated problem \eqref{asystem}  satisfy \eqref{Maj},
then the solution $u_{\varepsilon }$ remains bounded in $(L^2(0,T; H^1_0(\Omega)))^2$
 \end{proposition}

\begin{proof}
 We consider the test functions
$$
v_{\varepsilon }^{\gamma}=\beta_{\gamma}\varphi'(u_{\varepsilon }^{\gamma})
 \exp[\mu\psi(u_{\varepsilon})],\quad\text{for  }\gamma=1,2.
$$
 Noting that
 \begin{equation}\label{eq16}
 \nabla{\psi(u_{\varepsilon})}
=\sum_{\gamma=1}^2\beta_{\gamma}\varphi'(u_{\varepsilon }^{\gamma})
\nabla u_{\varepsilon }^{\gamma}.
 \end{equation}
We set
 \begin{equation*}
 I_{\varepsilon }=\sum_{\gamma=1}^2 \int_0^{T} \langle
\frac{\partial u_{\varepsilon }^{\gamma}}{\partial t},
\beta_{\gamma} \varphi'(u_{\varepsilon }^{\gamma})\exp[\mu{\psi(u_{\varepsilon})}]
\rangle dt.
 \end{equation*}
We use $v_{\varepsilon }^{\gamma}$ as test function in the $\gamma$-th equation
of system \eqref{asystem} and sum up from $\gamma = 1$ to $\gamma = 2$, we obtain
\begin{equation}\label{eq15}
 \begin{aligned}
&I_{\varepsilon } +\sum_{\gamma=1}^2\int_{Q} A(u_{\varepsilon })
 \nabla u_{\varepsilon }^{\gamma}\cdot\nabla u_{\varepsilon }^{\gamma}\beta_{\gamma}
 \varphi''(u_{\varepsilon }^{\gamma})\exp[\mu{
 \psi(u_{\varepsilon})}]dx\,dt \\
&+\mu\sum_{\gamma=1}^2\int_{Q}A(u_{\varepsilon })\nabla u_{\varepsilon }^{\gamma}
 \cdot\nabla{\psi(u_{\varepsilon})}\beta_{\gamma}
 \varphi'(u_{\varepsilon }^{\gamma})
 \exp[\mu{\psi(u_{\varepsilon})}]dx\,dt\\
&=\sum_{\gamma=1}^2\int_{Q}G_{\varepsilon }^{\gamma}(u_{\varepsilon },
 \nabla u_{\varepsilon })  \beta_{\gamma}\varphi'(u_{\varepsilon }^\gamma)
 \exp[\mu{\psi(u_{\varepsilon})}]dx\,dt
 \\
&\quad +\sum_{\gamma=1}^2 \int_{Q}F_{\varepsilon }(u_{\varepsilon },
\nabla u_{\varepsilon })
 \cdot\nabla u_{\varepsilon }^{\gamma}\beta_{\gamma}\varphi'(u_{\varepsilon }^\gamma)
\exp[\mu{\psi(u_{\varepsilon})}]dx\,dt.
 \end{aligned}
 \end{equation}
Firstly, we prove that
 \begin{gather}\label{rez}
 I_{\varepsilon }\geq-\frac{1}{\mu}\int_{\Omega}\exp\left[\mu\psi(u_0)\right] dx
 \end{gather}
Since (for $\gamma=1,2$)
 $u_{\varepsilon }^{\gamma}\in L^2(0,T;H_0^1(\Omega))\cap L^{\infty}(Q)$ and
 $\frac{\partial  u_{\varepsilon }^{\gamma}}{\partial t}\in  L^2(0,T;H^{-1}(\Omega))$,
using lemma \ref{lemma3.2}, there exists a sequence $w_j^{\gamma}$ such that
 $w_j^{\gamma}\in L^2(0,T;H_0^1(\Omega))$, $\frac{\partial
 w_j^{\gamma}}{\partial t}\in  L^2(0,T;H^1(\Omega))$, $w_j^{\gamma}$ bounded
in $L^{\infty}(Q)$ and
 \begin{gather}
 w_j^{\gamma}\to  u_{\varepsilon }^{\gamma} \quad \text{strongly in }
 L^2(0,T;H_0^1(\Omega)), \\
 \frac{\partial w_j^{\gamma}}{\partial t}\to \frac{\partial
 u_{\varepsilon }^{\gamma}}{\partial t}\quad \text{strongly in }
L^2(0,T;H^{-1}(\Omega)), \\
 w_j^{\gamma}(0)\to  u_{\varepsilon }^{\gamma}(0)=u_0^{\gamma} \quad
\text{strongly in }L^2(\Omega).
 \end{gather}
 Then
\begin{equation} \label{wj}
\begin{aligned}
&\sum_{\gamma=1}^2\int_0^{T} \langle\frac{\partial
 w_j^{\gamma}}{\partial t},
 \beta_{\gamma}\varphi'(w_j^{\gamma})\exp[\mu\psi(w_j)]\rangle dt\\
&=\frac{1}{\mu}\int_{\Omega}\exp\left[ {\mu\psi(w_j(T))}\right]dx
 -\frac{1}{\mu}\int_{\Omega}\exp\left[ {\mu\psi(w_j(0))}\right]dx \\
&\geq-\frac{1}{\mu}\int_{\Omega}\exp\left[ {\mu\psi(w_j(0))}\right]dx,
 \end{aligned}
\end{equation}
 then, we obtain \eqref{rez} by letting $j\to  \infty$ in \eqref{wj}.

Using \eqref{rez}, the coercivity condition \eqref{ellipticite} and the growth
conditions \eqref{G1G2Fcroissance}, \eqref{G1croissance} and \eqref{G2croissance}
on $G_{\varepsilon }^1,\;G_{\varepsilon }^2$ we obtain:
 \begin{align}
&\alpha \sum_{\gamma=1}^2\int_{Q}|\nabla u_{\varepsilon }^{\gamma}|^2
 \beta_{\gamma}\varphi''(u_{\varepsilon }^{\gamma})
 \exp[\mu{\psi(u_{\varepsilon})}]dx\,dt \nonumber\\
&+\alpha\mu\int_{Q}|\nabla{\psi(u_{\varepsilon})}|^2
 \exp[\mu{\psi(u_{\varepsilon})}]dx\,dt
\nonumber \\
 &\leq C_0\sum_{\gamma=1}^2\int_{Q}
 \beta_{\gamma}|\varphi'(u_{\varepsilon }^{\gamma})|
 \exp[\mu{\psi(u_{\varepsilon})}]dx\,dt
\nonumber \\
&\quad +\int_{Q}F_{\varepsilon }(u_{\varepsilon },\nabla u_{\varepsilon })
 \cdot\nabla{\psi(u_{\varepsilon})}
 \exp[\mu{\psi(u_{\varepsilon})}]dx\,dt
\nonumber \\
&\quad +C_2\int_{Q}|\nabla u_{\varepsilon }^1|^2\beta_1
 |\varphi'(u_{\varepsilon }^1)|\exp[\mu{\psi(u_{\varepsilon})}]dx\,dt
\nonumber \\
&\quad +\eta\int_{Q} |\nabla u_{\varepsilon }^2|^2
 \beta_1|\varphi'(u_{\varepsilon }^1)|
 \exp[\mu{\psi(u_{\varepsilon})}]dx\,dt
\nonumber \\
&\quad +C_2\int_{Q}|\nabla u_{\varepsilon }^1|^2\beta_2
 |\varphi'(u_{\varepsilon }^2)|\exp[\mu{\psi(u_{\varepsilon})}]dx\,dt
\nonumber \\
&\quad +C_2\int_{Q} |\nabla u_{\varepsilon }^2|^2\beta_2|\varphi'(u_{\varepsilon }^2)|
 \exp[\mu{\psi(u_{\varepsilon})}]dx\,dt
 +\frac{1}{\mu}\int_{\Omega}\exp\left[ {\mu\psi(u_0)}\right]dx.\label{eq17}
 \end{align}
 We estimate the second integral of the right hand side of \eqref{eq17} by using the
 growth conditions \eqref{Fcroissance} and \eqref{G1G2Fcroissance} on
$F_{\varepsilon }$ and Young's inequality,
 we obtain
 \begin{equation}\label{eq18}
 \begin{aligned}
& \int_{Q}F_{\varepsilon }(u_{\varepsilon },\nabla u_{\varepsilon })
 \cdot\nabla{\psi(u_{\varepsilon})}
 \exp[\mu{\psi(u_{\varepsilon})}]dx\,dt\\
&\leq \int_{Q}\left[ C_3+C_4|\nabla
 u_{\varepsilon }|\right]|\nabla{\psi(u_{\varepsilon})}|
 \exp[\mu{\psi(u_{\varepsilon})}]dx\,dt
 \\
&\leq\int_{Q}\big[\frac{\theta}{2}+\frac{C_3^2}{2\theta}|\nabla
 \psi(u_{\varepsilon})|^2+
 \frac{\theta}{2}|\nabla u_{\varepsilon }|^2
+\frac{C_4^2}{2\theta}|\nabla \psi(u_{\varepsilon})|^2\big]
\exp[\mu{\psi(u_{\varepsilon})}]dx\,dt
 \\
&=\int_{Q}\big[\frac{\theta}{2}+(\frac{C_3^2}{2\theta}
 +\frac{C_4^2}{2\theta})|\nabla{\psi(u_{\varepsilon})}|^2+
 \frac{\theta}{2}|\nabla u_{\varepsilon }|^2\big]
 \exp[\mu{\psi(u_{\varepsilon})}]dx\,dt,
 \end{aligned}
 \end{equation}
 using the hypothesis \eqref{data} on $\mu$ and \eqref{eq18},
the inequality \eqref{eq17} becomes
 \begin{align}
&\int_{Q}|\nabla u_{\varepsilon }^1|^2\exp[\mu{\psi(u_{\varepsilon})}]
 \underbrace{[\alpha\beta_1
 \varphi''(u_{\varepsilon }^1) -C_2\beta_1|\varphi'(u_{\varepsilon }^1)|
-C_2\beta_2|\varphi'(u_{\varepsilon }^2)|-\frac{\theta}{2}]}_{=\rho_1}dx\,dt
 \nonumber \\
& +\int_{Q}|\nabla u_{\varepsilon }^2|^2\exp[\mu{\psi(u_{\varepsilon})}]
 \underbrace{[\alpha\beta_2
 \varphi''(u_{\varepsilon }^2) -C_2{\beta_2}|\varphi'(u_{\varepsilon }^2)|
-{\beta_1}\eta|\varphi'(u_{\varepsilon }^1)|-\frac{\theta}{2}]}_{=\rho_2}dx\,dt
 \nonumber  \\
& +\int_{Q}|\nabla{\psi(u_{\varepsilon})}|^2
 \underbrace{[\alpha\mu-\frac{C_3^2}{2\theta}-\frac{C_4^2}{2\theta}]}_{\geq0}
\exp[\mu{\psi(u_{\varepsilon})}]dx\,dt
\nonumber  \\
&\leq C_0\sum_{\gamma=1}^2\int_{Q}
 \beta_{\gamma}|\varphi'(u_{\varepsilon }^{\gamma})
|\exp[\mu{\psi(u_{\varepsilon})}]dx\,dt
 +\int_{Q}\frac{\theta}{2}\exp[\mu{\psi(u_{\varepsilon})}]dx\,dt
\nonumber   \\
&\quad +\frac{1}{\mu}\int_{\Omega}\exp\left[ {\mu\psi(u_0)}\right] dx.  \label{21}
 \end{align}
Employing \eqref{21} and the lemma \eqref{lemma3.1} yields
 \begin{equation} \label{alpha0}
\begin{aligned}
&\alpha_0\sum_{\gamma=1}^2\int_{Q}|\nabla u_{\varepsilon }^{\gamma}|^2
 \exp[\mu{\psi(u_{\varepsilon})}]dx\,dt \\
&\leq\int_{Q}\frac{\theta}{2}\exp[\mu{\psi(u_{\varepsilon})}]dx\,dt
 +C_0\sum_{\gamma=1}^2\int_{Q}  \beta_{\gamma}|\varphi'(u_{\varepsilon }^{\gamma})
 |\exp[\mu{\psi(u_{\varepsilon})}]dx\,dt \\
&\quad +\frac{1}{\mu}\int_{\Omega}\exp\left[ {\mu\psi(u_0)}\right] dx,
 \end{aligned}
\end{equation}
using the facts that
 $\exp[\mu{\psi(u_{\varepsilon})}]\geq1$ and that $u_{\varepsilon }$
satisfy \eqref{Maj}
 (which implies that ${\psi(u_{\varepsilon})}$, is bounded in
 $L^{\infty}(Q)$)
and $u_0\in  (L^{\infty}(\Omega))^2$, implies that $u_{\varepsilon }$ is bounded in
 $(L^2(0,T;H_0^1(\Omega)))^2$.
\end{proof}

 Since by proposition \ref{proposition1} $u_{\varepsilon }$ remains bounded
in $(L^2(0,T;H_0^1(\Omega)))^2$, we can extract a subsequence, still denoted
 by $u_{\varepsilon }$, such that
 \begin{equation}\label{weak convergence}
 u_{\varepsilon }\rightharpoonup u\quad \text{in }(L^2(0,T;H_0^1(\Omega)))^2.
 \end{equation}

 \begin{proposition}\label{proposition2}
 The sequence $u_{\varepsilon }$ is relatively compact in  $(L^2(Q))^2$.
 \end{proposition}

\begin{proof}
For $s>0$ large enough, we have $L^1(\Omega)\subset H^{-s}(\Omega)$.
Then $\frac{\partial u_{\varepsilon }}{\partial t}$ is bounded in
$(L^1(0,T; H^{-s}(\Omega))^2$
and $u_{\varepsilon }$ is bounded in $(L^2(0,T;H_0^1(\Omega)))^2$.

As $H_0^1(\Omega)\subset L^2(\Omega)\subset H^{-s}(\Omega)$,
the injection being compact, proposition \ref{proposition2} follows from
a compacticity lemma of Aubin's type. Such a lemma can be found for example
in  \cite[p. 271]{Temam}  or in \cite[section 8, corollary 4)]{Si}.
Therefore we can extract a subsequence, still denoted by $u_{\varepsilon }$
such that if $\varepsilon \to  0$
 \begin{equation}\label{Hweakconv}
 u_{\varepsilon }\rightharpoonup u\quad \text{in}\;(L^2(0,T;H_0^1(\Omega)))^2.
 \end{equation}
 By possibly extracting a subsequence, we can suppose, without loss of generality
 using proposition \ref{proposition2}, that
 \begin{gather}\label{L2conv}
 u_{\varepsilon }\to  u\quad  \text{in }(L^2(Q))^2, \\
 \label{almostconv}
 u_{\varepsilon }\to  u\quad \text{a.e. in } Q
 \end{gather}
\end{proof}

 \section{Strong convergence in $(L^2(0,T;H_0^1(\Omega)))^2$}

 To pass to the limit as $\varepsilon \to  0$ in the nonlinearities $G_{\varepsilon }^1(u_{\varepsilon }, \nabla u_{\varepsilon } ),
 \\G_{\varepsilon }^2(u_{\varepsilon }, \nabla u_{\varepsilon } )$ and $F_{\varepsilon }(u_{\varepsilon }, \nabla u_{\varepsilon } )$ in system \eqref{asystem}, we need the strong convergence of $u_{\varepsilon }\to  u$ in $(L^2(0,T;H_0^1(\Omega)))^2$. This is our goal in this step.
 \begin{proposition}\label{proposition3}
 Assume that \eqref{ellipticite}--\eqref{Fcroissance} and \eqref{G1G2Fcroissance} hold true. If the
 solutions $u_{\varepsilon }$ of the approximated problem \eqref{asystem}  satisfy \eqref{Maj}  and \eqref{Hweakconv}--\eqref{almostconv}
 then $u_{\varepsilon }$ converges strongly to $u$ in $(L^2(0,T;H_0^1(\Omega)))^2$.
 \end{proposition}
 We consider the functions
$\bar{\varphi}:\mathbb{R}\to \mathbb{R}$ and 
$\bar{\psi}:\mathbb{R}^2\to \mathbb{R}$
 defined by:
 \begin{gather*}\label{varphi}
 \bar{\varphi}(\tau)=e^{\bar{\lambda}\tau}+e^{-\bar{\lambda}\tau}-2,\quad\forall\tau\in\mathbb{R}, \\
 \bar{\psi}(s)=\bar{\beta}_1\bar{\varphi}(s^1)+\bar{\beta}_2\bar{\varphi}(s^2)
 ,\quad\forall s\in\mathbb{R}^2,
 \end{gather*}
 where $\bar{\lambda}$, $\bar{\mu}$, $\bar{\beta}_1$ and $\bar{\beta}_2$ are
 positive constants defined by
 \begin{equation}\label{bardata}
 \bar{\lambda}=\frac{16C_2}{\alpha},\quad
 \bar{\beta}_2 =\frac{\bar{\beta}_1}{2\exp(2\bar{\lambda}M)},\quad
\bar{\mu}=\frac{2}{\alpha}\Big(\frac{C_3^2}{\bar{\theta}}
+\frac{C_4^2}{\bar{\theta}}\Big)
\end{equation}
with $\bar{\theta}$ a fixed number such that
$0<\bar{\theta}\leq\frac{4C_2^2
 \bar{\beta}_1}{\alpha\exp(2\bar{\lambda }M)}$.

\begin{lemma}\label{lemma4.2}
For any $u_{\varepsilon }$ and $u_{\nu}$, such that
$|u_{\varepsilon }^{\delta}-u_{\nu}^{\delta}|\leq 2M$ for any $\delta=1,2$, we have
 \begin{equation}\label{rho1bar}
\begin{aligned}
 \bar{\rho}_1&=\alpha\bar{\beta}_1\bar{\varphi}''(u_{\varepsilon }^1-u_{\nu}^1)\\
&\quad -4\underbrace{\left(2C_2\bar{\beta}_1|\bar{\varphi}'(u_{\varepsilon }^1
-u_{\nu}^1)|  +2C_2\bar{\beta}_2 |\bar{\varphi}'(u_{\varepsilon }^2-u_{\nu}^2)|
 +\bar{\theta}\right)}_{=L_1(u_{\varepsilon},u_{\nu})}\geq\bar{\alpha}_0
\end{aligned}
 \end{equation}
 and
 \begin{equation}\label{rho2bar}
 \bar{\rho}_2=\alpha\bar{\beta}_2\bar{\varphi}''(u_{\varepsilon }^2-u_{\nu}^2)
-4\underbrace{\big(2C_2\bar{\beta}_2|\bar{\varphi}'(u_{\varepsilon }^2-u_{\nu}^2)|
 +2\eta\bar{\beta}_1 |\bar{\varphi}'(u_{\varepsilon }^1-u_{\nu}^1)|
 +\bar{\theta}\big)}_{=L_2(u_{\varepsilon },u_{\nu})}\geq\bar{\alpha}_0
 \end{equation}
 where
 \begin{equation}\label{baralpha}
 \bar{\alpha}_0=\frac{16C_2^2\bar{\beta}_1}
 {\alpha\exp(2\lambda M)}.
 \end{equation}
\end{lemma}

The proof of the above lemma  is the same of the proof of \eqref{lemma3.1}
where $\varphi$, $\lambda$, $\beta_1$, and $\beta_2$ are
replaced by $\bar{\varphi}$, $\bar{\lambda}$, $\bar{\beta_1}$ and $\bar{\beta_2}$;
and where $C_2,\eta,\theta$ and $M$ are replaced by
$8C_2,8\eta,8\theta$ and $2M$.


\begin{proof}[Proof of proposition \ref{proposition3}]
 If $\varepsilon $ and $\nu$ are two parameters, we write \eqref{asystem} as
\begin{equation}\label{nueps}
 \begin{gathered}
 \begin{aligned}
 &\frac{\partial( u_{\varepsilon }^1-u_{\nu}^1)}{\partial t}
-\operatorname{div}(A(u_{\varepsilon })\nabla(u_{\varepsilon }^1-u_{\nu}^1))
+\frac{\partial u_{\nu}^1}{\partial t}-\operatorname{div}(A(u_{\varepsilon })
\nabla u_{\nu}^1)
 \\
&=G_{\varepsilon }^1(u_{\varepsilon },\nabla u_{\varepsilon })+
 F_{\varepsilon }(u_{\varepsilon },\nabla u_{\varepsilon })
\cdot\nabla u_{\varepsilon }^1\quad \text{in }Q,
\end{aligned} \\
\begin{aligned}
& \frac{\partial( u_{\varepsilon }^2-u_{\nu}^2)}{\partial t} -\operatorname{div}
(A(u_{\varepsilon })\nabla(u_{\varepsilon }^2-u_{\nu}^2))
+\frac{\partial u_{\nu}^2}{\partial t}-\operatorname{div}(A(u_{\varepsilon })
\nabla u_{\nu}^2)\\
&=G_{\varepsilon }^2(u_{\varepsilon },\nabla u_{\varepsilon })
+  F_{\varepsilon }(u_{\varepsilon },\nabla u_{\varepsilon })
 \cdot\nabla u_{\varepsilon }^2\quad \text{in }Q,
\end{aligned} \\
 u_{\varepsilon }^1-u_{\nu}^1=0,\quad
u_{\varepsilon }^2-u_{\nu}^2=0,\quad \text{on }\Sigma, \\
(u_{\varepsilon }^1-u_{\nu}^1)(0)=0,\quad
(u_{\varepsilon }^2-u_{\nu}^2)(0)=0,\quad \text{in }\Omega.
 \end{gathered}
 \end{equation}
We consider the test functions
 \begin{equation*}
 \bar{v}_{\varepsilon \nu}^{\gamma}=\bar{\beta}_{\gamma}
\bar{\varphi}'(u_{\varepsilon }^{\gamma}-u_{\nu}^{\gamma})
 \exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}],\quad \gamma=1,2
 \end{equation*}
Noting that
 \begin{equation}\label{Dbarpsi}
 \nabla{\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}=\sum_{\gamma=1}^2
\bar{\beta}_{\gamma}\bar{\varphi}'(u_{\varepsilon }^{\gamma}-u_{\nu}^{\gamma})
\nabla(u_{\varepsilon }^{\gamma}-u_{\nu}^{\gamma})
 \end{equation}

 Using $\bar{v}_{\varepsilon \nu}^{\gamma}$ as test function in the
 $\gamma$-th equation of system \eqref{nueps} and summing from
 $\gamma=1$ to $\gamma=2$, we obtain
 \begin{align}
&\underbrace{\sum_{\gamma=1}^2\int_0^{T}\langle
 \frac{\partial( u_{\varepsilon }^{\gamma}-u_{\nu}^{\gamma})}{\partial t},
\bar{\beta}_{\gamma}
 \bar{\varphi}'(u_{\varepsilon }^{\gamma}
-u_{\nu}^{\gamma})\exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }
-{u}_{\nu})}]\rangle dt}_{=J_1(\varepsilon )} \nonumber\\
& +\sum_{\gamma=1}^2\int_{Q}A(u_{\varepsilon })
 \nabla(u_{\varepsilon }^{\gamma}-u_{\nu}^{\gamma})
 \cdot\nabla(u_{\varepsilon }^{\gamma}-u_{\nu}^{\gamma})
 \bar{\beta}_{\gamma}\bar{\varphi}''(u_{\varepsilon }^{\gamma}
 -u_{\nu}^{\gamma})\exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}]dx\,dt
 \nonumber\\
&+\bar{\mu}\sum_{\gamma=1}^2\int_{Q}A(u_{\varepsilon })
 \nabla(u_{\varepsilon }^{\gamma}-u_{\nu}^{\gamma})
 \cdot\nabla{\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}
 \bar{\beta}_{\gamma}\bar{\varphi}'(u_{\varepsilon }^{\gamma}
 -u_{\nu}^{\gamma})\exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}]dx\,dt
 \nonumber\\
&+\sum_{\gamma=1}^2\int_0^{T}\langle\frac{\partial u_{\nu}^{\gamma}}{\partial t},
 \bar{\beta}_{\gamma}\bar{\varphi}'(u_{\varepsilon }^{\gamma}
 -u_{\nu}^{\gamma})\exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }
 -{u}_{\nu})}]\rangle dt
 \nonumber\\
&+\sum_{\gamma=1}^2\int_{Q}A(u_{\varepsilon })\nabla u_{\nu}^{\gamma}
 \cdot\nabla\left[\bar{\beta}_{\gamma}\bar{\varphi}'(u_{\varepsilon }^{\gamma}
 -u_{\nu}^{\gamma})\exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}]
 \right]dx\,dt
 \nonumber\\
&=\sum_{\gamma=1}^2\int_{Q}G_{\varepsilon }^{\gamma}(u_{\varepsilon },
 \nabla u_{\varepsilon })\bar{\beta}_{\gamma}\bar{\varphi}'
 (u_{\varepsilon }^{\gamma}-u_{\nu}^{\gamma})\exp[\bar{\mu}
 {\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}]dx\,dt
 \nonumber
\\
&\quad +  \int_{Q} F_{\varepsilon }(u_{\varepsilon },\nabla
 u_{\varepsilon })\cdot
\nabla\bar\psi(u_\varepsilon -u_\nu)
%\nabla\bar\psi_{\varepsilon \nu}
 \exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}]dx\,dt \nonumber
\\
 &\quad + \sum_{\gamma=1}^2\int_{Q}
 F_{\varepsilon }(u_{\varepsilon },\nabla
 u_{\varepsilon })\cdot\nabla u_{\nu}^{\gamma}\bar{\beta}_{\gamma}
 \bar{\varphi}'(u_{\varepsilon }^{\gamma}-u_{\nu}^{\gamma})
 \exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}]dx\,dt
 \nonumber\\
 &=I_1(\varepsilon )+I_2(\varepsilon )+I_3(\varepsilon ) \label{int}
 \end{align}

We claim that the first term of \eqref{int} is nonnegative. Indeed
 as
 $(u_{\varepsilon }^{\gamma}-u_{\nu}^{\gamma})\in L^2(0,T;H_0^1(\Omega))\cap L^{\infty}(Q)$ and $\frac{\partial
 (u_{\varepsilon }^{\gamma}-u_{\nu}^{\gamma})}{\partial t}\in
 L^2(0,T;H^{-1}(\Omega))$ for  $\gamma=1,2$, then there exists a sequence
$w_j^{\gamma}$ such
 that
 $w_j^{\gamma}\in L^2(0,T;H_0^1(\Omega))$, $\frac{\partial
 w_j^{\gamma}}{\partial t}\in
 L^2(0,T;H^1(\Omega))$, $w_j^{\gamma}$ bounded in $L^{\infty}(Q)$
and
 \begin{gather}\label{w}
 w_j^{\gamma}\to  u_{\varepsilon }^{\gamma}-u_{\nu}^{\gamma} \quad
\text{strongly in } L^2(0,T;H_0^1(\Omega)), \\
\label{partialw}
 \frac{\partial w_j^{\gamma}}{\partial t}\to \frac{\partial
 (u_{\varepsilon }^{\gamma}-u_{\nu}^{\gamma})}{\partial t}\quad
\text{strongly in }L^2(0,T;H^{-1}(\Omega)), \\
 w_j^{\gamma}(0)\to  (u_{\varepsilon }^{\gamma}-u_{\nu}^{\gamma})(0)=0 \quad
\text{strongly in }L^2(\Omega).
 \end{gather}
 By \eqref{w}, \eqref{partialw} and the the continuous injection of the space
 \begin{equation}
 W(0,T)=\{ v\in L^2(0,T;H_0^1(\Omega)),\;\frac{\partial v}{\partial t}
 \in L^2(0,T;H^{-1}(\Omega))\}
 \end{equation}
in $\mathcal{C}([0,T];L^2(\Omega))$, we obtain
 \begin{gather*}
 w_j^{\gamma}(T)\to  (u_{\varepsilon }^{\gamma}-u_{\nu}^{\gamma})(T)\quad
 \text{strongly in }L^2(\Omega), \\
 \begin{aligned}
&\sum_{\gamma=1}^2\int_0^{T} \langle\frac{\partial
 w_j^{\gamma}}{\partial t},
 \bar{\beta}_{\gamma}\varphi'(w_j^{\gamma})
 \exp[\bar{\mu}\bar{\psi}(w_j)]\rangle  dt\\
&= \frac{1}{\bar\mu}\int_{Q}\exp[{\bar{\mu}\bar{\psi}(w_j(T))}] dx
 -\frac{1}{\bar\mu}\int_{Q}\exp[{\bar{\mu}\bar{\psi}(w_j(0))}] dx
 \end{aligned}
 \end{gather*}
and letting $j\to +\infty$ shows that
 \begin{equation}\label{J1eps}
 J_1(\varepsilon )= \frac{1}{\bar\mu}\int_{Q}\{
 \exp[\bar{\mu}\bar\psi((u_{\varepsilon }^{\gamma}
-u_{\nu}^{\gamma})(T))]-1\} dx\geq0.
 \end{equation}
 Secondly we estimate various terms of the right hand side of \eqref{int}.
For the third term we have by  using the growth conditions \eqref{Fcroissance},
 \eqref{G1G2Fcroissance} on $ F_{\varepsilon }$ and Young's inequality:
 \begin{equation}\label{I3eps}
 \begin{split}
I_3(\varepsilon )&\leq\sum_{\gamma=1}^2\int_{Q}[C_3+C_4|\nabla u_{\varepsilon }|]
|\nabla u_{\nu}^{\gamma}|\bar{\beta}_{\gamma}|\bar{\varphi}'(u_{\varepsilon }^{\gamma}
-u_{\nu}^{\gamma})|
 \exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}]dx\,dt\\
&\leq C_3\sum_{\gamma=1}^2\int_{Q}|
 \nabla u_{\nu}^{\gamma}|\bar{\beta}_{\gamma}|\bar{\varphi}'(u_{\varepsilon }^{\gamma}
-u_{\nu}^{\gamma})|
 \exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}]dx\,dt\\
&\quad +C_4\sum_{\gamma=1}^2\int_{Q}|\nabla u_{\varepsilon }|
 |\nabla u_{\nu}^{\gamma}|\bar{\beta}_{\gamma}|\bar{\varphi}'
 (u_{\varepsilon }^{\gamma}-u_{\nu}^{\gamma})|
 \exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}]dx\,dt\\
&=I_3^1(\varepsilon )+I_3^2(\varepsilon )
 \end{split}
 \end{equation}
Concerning the second term, we use \eqref{Dbarpsi} and the growth conditions
\eqref{Fcroissance},  \eqref{G1G2Fcroissance} on $ F_{\varepsilon }$ and Young's
inequality, we obtain
 \begin{equation}\label{I2eps}
 \begin{aligned}
 I_2(\varepsilon )
&\leq\int_{Q}[C_3+C_4|\nabla u_{\varepsilon }|]
 |\nabla{\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}| 
  \exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}]dx\,dt
 \\
&\leq\int_{Q}\big[\frac{\bar{\theta}}{2}+\frac{C_3^2}{2\bar{\theta}}
 |\nabla{\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}|^2+
 \frac{\bar{\theta}}{2}|\nabla u_{\varepsilon }|^2+\frac{C_4^2}{2\bar{\theta}}|
 \nabla{\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}|^2\big] \\
&\quad\times \exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}]dx\,dt
 \\
&\leq\bar{\theta}\int_{Q}\Big(\frac{1}{2}+|\nabla u_{\nu}|^2\Big)
 \exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}]dx\,dt
 \\
&\quad +\Big(\frac{C_3^2}{2\bar{\theta}}+
 \frac{C_4^2}{2\bar{\theta}}\Big)\int_{Q}
 |\nabla{\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}|^2
\exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}]dx\,dt\\
&\quad +\bar{\theta}\int_{Q}|\nabla(u_{\varepsilon }-u_{\nu})|^2
 \exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}]dx\,dt\\
&=I_2^1(\varepsilon )+I_2^2(\varepsilon )+I_2^{3}(\varepsilon).
 \end{aligned}
 \end{equation}
Then, we estimate the first term by using the growth conditions
\eqref{G1croissance}, \eqref{G2croissance} and \eqref{G1G2Fcroissance} on
$G_{\varepsilon }^1$, $G_{\varepsilon }^2$ to obtain
 \begin{equation}\label{I1eps}
 \begin{aligned}
 I_1(\varepsilon )
&\leq C_0\sum_{\gamma=1}^2\int_{Q}
 \bar{\beta}_{\gamma}|\bar{\varphi}'(u_{\varepsilon }^{\gamma}-u_{\nu}^{\gamma})|
\exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}]dx\,dt\\
&\quad +C_2\int_{Q}|\nabla u_{\varepsilon }^1|^2\bar{\beta}_1|\bar{\varphi}'
 (u_{\varepsilon }^1-u_{\nu}^1)| \exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }
 -{u}_{\nu})}]dx\,dt\\
&\quad +\eta\int_{Q}|\nabla u_{\varepsilon }^2|^2\bar{\beta}_1|\bar{\varphi}'
 (u_{\varepsilon }^1-u_{\nu}^1)| \exp[\bar{\mu}{\bar{\psi}
 ({u}_{\varepsilon }-{u}_{\nu})}]dx\,dt\\
&\quad +C_2\int_{Q}|\nabla u_{\varepsilon }^1|^2\bar{\beta}_2
 |\bar{\varphi}'(u_{\varepsilon }^2-u_{\nu}^2)|
 \exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}]dx\,dt\\
&\quad  +C_2\int_{Q}|\nabla u_{\varepsilon }^2|^2\bar{\beta}_2
 |\bar{\varphi}'(u_{\varepsilon }^2-u_{\nu}^2)|
  \exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}]dx\,dt\\
&=I_1^1(\varepsilon )+I_1^2(\varepsilon )+I_1^{3}(\varepsilon )
 +I_1^{4}(\varepsilon )+I_1^{5}(\varepsilon).
 \end{aligned}
 \end{equation}
 Now we estimate the four last terms of the right hand side of the inequality
 \eqref{I1eps}; for what concerns the second term we have by using the
Young's inequality
 \begin{align*}%\label{I12}
 I_1^2(\varepsilon )
&\leq 2C_2\int_{Q}|\nabla(u_{\varepsilon }^1-u_{\nu}^1)|^2\bar{\beta}_1
 |\bar{\varphi}'(u_{\varepsilon }^1-u_{\nu}^1)|
 \exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}]dx\,dt
 \\
&\quad +2C_2\int_{Q}|\nabla u_{\nu}^1|^2\bar{\beta}_1|\bar{\varphi}'
 (u_{\varepsilon }^1-u_{\nu}^1)| \exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }
 -{u}_{\nu})}]dx\,dt=I_{11}^2(\varepsilon )+I_{12}^2(\varepsilon )
 \end{align*}
concerning the third term we have
 \begin{equation}\label{I13}
 \begin{split}
 I_1^{3}(\varepsilon )
&\leq 2\eta\int_{Q}|\nabla(u_{\varepsilon }^2-u_{\nu}^2)|^2
 \bar{\beta}_1|\bar{\varphi}'(u_{\varepsilon }^1-u_{\nu}^1)
 | \exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}]dx\,dt
 \\
&\quad +2\eta\int_{Q}|\nabla u_{\nu}^2|^2\bar{\beta}_1
 |\bar{\varphi}'(u_{\varepsilon }^1-u_{\nu}^1)
 | \exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}]dx\,dt \\
&=I_{11}^{3}(\varepsilon )+I_{12}^{3}(\varepsilon )
 \end{split}
 \end{equation}
we estimate the fourth term as
 \begin{equation}\label{I14}
 \begin{split}
 I_1^{4}(\varepsilon )
&\leq 2C_2\int_{Q}|\nabla(u_{\varepsilon }^1-u_{\nu}^1)|^2
 \bar{\beta}_2|\bar{\varphi}'(u_{\varepsilon }^2-u_{\nu}^2)|
 \exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}]dx\,dt
 \\
&\quad +2C_2\int_{Q}|\nabla u_{\nu}^1|^2\bar{\beta}_2
 |\bar{\varphi}'(u_{\varepsilon }^2-u_{\nu}^2)|
 \exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}]dx\,dt\\
&=I_{11}^{4}(\varepsilon )+I_{12}^{4}(\varepsilon )
 \end{split}
 \end{equation}
for the last term we have
 \begin{equation}\label{I15}
 \begin{split}
 I_1^{5}(\varepsilon )
&\leq 2C_2\int_{Q}|\nabla(u_{\varepsilon }^2-u_{\nu}^2)|^2
 \bar{\beta}_2|\bar{\varphi}'(u_{\varepsilon }^2-u_{\nu}^2)
| \exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}]dx\,dt
 \\
&\quad +2C_2\int_{Q}|\nabla u_{\nu}^2|^2\bar{\beta}_2
 |\bar{\varphi}'(u_{\varepsilon }^2-u_{\nu}^2)|
 \exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}]dx\,dt \\
&=I_{11}^{5}(\varepsilon )+I_{12}^{5}(\varepsilon).
 \end{split}
 \end{equation}
Set afterwards
 \begin{gather}\label{J2eps}
 J_2(\varepsilon )=\sum_{\gamma=1}^2\int_{Q}A(u_{\varepsilon })
\nabla u_{\nu}^{\gamma}\nabla\left[\bar{\beta}_{\gamma}\bar{\varphi}'
(u_{\varepsilon }^{\gamma}-u_{\nu}^{\gamma})
\exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}]\right]dx\,dt,
\\
\label{J3ep}
 J_3(\varepsilon )=\sum_{\gamma=1}^2\int_0^{T}
 \langle\frac{\partial u_{\nu}^{\gamma}}{\partial t},
 \bar{\beta}_{\gamma}\bar{\varphi}'(u_{\varepsilon }^{\gamma}
 -u_{\nu}^{\gamma})\exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }
 -{u}_{\nu})}]\rangle dt
 \end{gather}
 Now we consider the matrices
 \begin{gather}\label{Bepsilon1}
\begin{aligned}
 B_{\varepsilon \nu}^1
&=\bar{\beta}_1\bar{\varphi}''(u_{\varepsilon }^1
 -u_{\nu}^1)\exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }
 -{u}_{\nu})}]A(u_{\varepsilon }) \\
&\quad - L_1(u_{\varepsilon },u_{\nu})\exp[\bar{\mu}{\bar{\psi}
 ({u}_{\varepsilon }-{u}_{\nu})}]Id,
\end{aligned}\\
\label{Bepsilon2}
\begin{aligned}
 B_{\varepsilon \nu}^2
&=\bar{\beta}_2\bar{\varphi}''(u_{\varepsilon }^2
 -u_{\nu}^2)\exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }
 -{u}_{\nu})}]A(u_{\varepsilon }) \\
&\quad - L_2(u_{\varepsilon },u_{\nu})\exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }
 -{u}_{\nu})}]Id,
\end{aligned}\\
\label{Eepsilon}
 E_{\varepsilon \nu}=\bar{\mu}\exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }
-{u}_{\nu})}]A(u_{\varepsilon })-\Big(\frac{C_3^2}{2\bar{\theta}}
+\frac{C_3^2}{2\bar{\theta}}\Big)\exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }
-{u}_{\nu})}] Id.
 \end{gather}
 where $Id$ denotes the identity matrix.
Inequalities \eqref{int} and \eqref{J1eps}--\eqref{Eepsilon} give
 \begin{equation}\label{inqeps}
 \begin{split}
&J_2(\varepsilon )+J_3(\varepsilon )+\int_{Q}B_{\varepsilon \nu}^1
 \nabla(u_{\varepsilon }^1-u_{\nu}^1)\cdot\nabla(u_{\varepsilon }^1-u_{\nu}^1)dx\,dt\\
&+\int_{Q}B_{\varepsilon \nu}^2
 \nabla(u_{\varepsilon }^2-u_{\nu}^2)\cdot\nabla(u_{\varepsilon }^2-u_{\nu}^2)dx\,dt
 + \int_{Q}E_{\varepsilon \nu}  \nabla{\bar{\psi}({u}_{\varepsilon }
 -{u}_{\nu})}\cdot\nabla{\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}dx\,dt\\
&\leq I_1^1(\varepsilon )+I_{12}^2(\varepsilon )+I_{12}^{3}(\varepsilon )
 +I_{12}^{4}(\varepsilon )+  I_{12}^{5}(\varepsilon )+I_3^1(\varepsilon )
 +I_3^2(\varepsilon )+I_2^1(\varepsilon).
 \end{split}
 \end{equation}
Now we pass to the limit in \eqref{inqeps} for $\varepsilon \to  0$ and $\nu$ fixed.
Using the fact that
$\bar{\beta}_{\gamma}\bar{\varphi}'(u_{\varepsilon }^{\gamma}
-u_{\nu}^{\gamma})\exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}]$
remains bounded in $L^2(0,T;H_0^1(\Omega))$ and converges almost everywhere in
$Q$ towards $\bar{\beta}_{\gamma}\bar{\varphi}'(u^{\gamma}-u_{\nu}^{\gamma})
\exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]$.
 We have for $\gamma=1,2$,
 \begin{equation}\label{weak}
 \bar{\beta}_{\gamma}\bar{\varphi}'(u_{\varepsilon }^{\gamma}
-u_{\nu}^{\gamma})\exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}]
\to  \bar{\beta}_{\gamma}\bar{\varphi}'(u^{\gamma}
-u_{\nu}^{\gamma})\exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]
 \end{equation}
  weakly in $L^2(0,T;H_0^1(\Omega))$, and strongly in $L^2(Q)$.

Using \eqref{uniformborne}, \eqref{almostconv} and Lebesgue's dominated convergence
 theorem, we obtain
 \begin{equation}\label{str}
 A(u_{\varepsilon })\nabla u_{\nu}^{\gamma}\to  A(u)\nabla u_{\nu}^{\gamma}
\quad\text{strongly in }(L^2(Q))^{N}\,.
\end{equation}
 By \eqref{weak} and \eqref{str} we have
 \begin{equation}\label{J2epsi}
 J_2(\varepsilon )=\sum_{\gamma=1}^2\int_{Q}A(u_{\varepsilon })
\nabla u_{\nu}^{\gamma}\cdot\nabla\big[\bar{\beta}_{\gamma}
\bar{\varphi}'(u_{\varepsilon }^{\gamma}-u_{\nu}^{\gamma})
\exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}]\big]dx\,dt\to  J_2
 \end{equation}
 where
 \begin{gather}
 J_2=\sum_{\gamma=1}^2\int_{Q}A(u)\nabla u_{\nu}^{\gamma}\cdot\nabla
\big[\bar{\beta}_{\gamma}\bar{\varphi}'(u^{\gamma}-u_{\nu}^{\gamma})
\exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]\big]dx\,dt, \\
\label{J3eps}
 J_3(\varepsilon )=\sum_{\gamma=1}^2\int_0^{T}\langle
\frac{\partial u_{\nu}^{\gamma}}{\partial t},
\bar{\beta}_{\gamma}\bar{\varphi}'(u_{\varepsilon }^{\gamma}
-u_{\nu}^{\gamma})\exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }
-{u}_{\nu})}]\rangle dt\to  J_3,
 \end{gather}
 where
 \begin{equation}
 J_3=\sum_{\gamma=1}^2\int_0^{T}\langle\frac{\partial u_{\nu}^{\gamma}}{\partial t},
\bar{\beta}_{\gamma}\bar{\varphi}'(u^{\gamma}-u_{\nu}^{\gamma})
\exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]\rangle dt.
 \end{equation}
The matrices $B_{\varepsilon \nu}^1$, $B_{\varepsilon \nu}^2$, and
$E_{\varepsilon \nu}$ are positive definite, because of the coercivity
condition \eqref{ellipticite}, \eqref{rho1bar} and \eqref{rho2bar} imply that
 \begin{equation}
 B_{\varepsilon \nu}^1\geq \bar{\alpha}_0Id,\quad
B_{\varepsilon \nu}^2\geq \bar{\alpha}_0Id.
 \end{equation}
 Using \eqref{ellipticite} and the definition of $\bar{\mu}$ in \eqref{bardata},
we obtain
 \begin{equation}
 E_{\varepsilon \nu}\geq \frac{\bar{\mu}}{2}Id.
 \end{equation}
 Using \eqref{almostconv} we obtain that the matrix
$B_{\varepsilon \nu}^1$ converges a.e. in $Q$ to the matrix $B_{\nu}^1$ defined by
 \begin{equation}\label{Bnu1}
 B_{\nu}^1=\bar{\beta}_1\bar{\varphi}''(u^1-u_{\nu}^1)
\exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]A(u)-
 L_1(u,u_{\nu})\exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]Id,
 \end{equation}
and the matrix $B_{\varepsilon \nu}^2$ converges a.e. in $Q$ to the matrix
$B_{\nu}^2$ defined by
 \begin{equation}\label{Bnu2}
 B_{\nu}^2=\bar{\beta}_2\bar{\varphi}''(u^2-u_{\nu}^2)\exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]A(u)-
 L_2(u,u_{\nu})\exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]Id,
 \end{equation}
and the matrix $E_{\varepsilon \nu}$ converges a.e. in $Q$ to
 the matrix $E_{\nu}$ defined by
 \begin{equation}
 E_{\nu}=\bar{\mu}\exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]A(u)
-\Big(\frac{C_3^2}{2\bar{\theta}}+\frac{C_3^2}{2\bar{\theta}} \Big)
\exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}] Id.
 \end{equation}
Thanks to the positive definiteness of $B_{\varepsilon \nu}^1$,
$B_{\varepsilon \nu}^2$
 and $E_{\varepsilon \nu}$ we obtain
 \begin{gather*}
 \liminf _{\varepsilon \to
 0}\int_{Q}B_{\varepsilon \nu}^1
 \nabla(u_{\varepsilon }^1-u_{\nu}^1)\cdot\nabla(u_{\varepsilon }^1-u_{\nu}^1)dx\,dt
 \geq  \int_{Q}B_{\nu}^1
 (\nabla u^1-u_{\nu}^1)\cdot\nabla(u^1-u_{\nu}^1)dx\,dt
\\
 \liminf _{\varepsilon \to  0}\int_{Q}B_{\varepsilon \nu}^2
 \nabla(u_{\varepsilon }^2-u_{\nu}^2)\cdot\nabla(u_{\varepsilon }^2-u_{\nu}^2)dx\,dt
 \geq  \int_{Q}B_{\nu}^2
 (\nabla u^2-u_{\nu}^2)\cdot\nabla(u^2-u_{\nu}^2)dx\,dt
\\
\begin{aligned}
&\liminf _{\varepsilon \to  0}\int_{Q}E_{\varepsilon \nu}
 \nabla{\bar{\psi}({u}_{\varepsilon }
 -{u}_{\nu})}\cdot\nabla{\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}dx\,dt\\
&\geq  \int_{Q}E_{\nu}\nabla{\bar{\psi}(u-u_{\nu})}\cdot
 \nabla{\bar{\psi}(u-u_{\nu})}dx\,dt.
\end{aligned}
 \end{gather*}
To pass to the limit in the right hand side of inequality \eqref{inqeps},
 we use Lebesgue's dominated convergence theorem which implies that
 \begin{align*}
& I_3^1(\varepsilon )+I_2^1(\varepsilon )+I_1^1(\varepsilon )
 +\sum_{i=2}^{5} I_{12}^{i}(\varepsilon ) \\
&=  C_3\sum_{\gamma=1}^2\int_{Q}|\nabla u_{\nu}^{\gamma}
 |\bar{\beta}_{\gamma}|\bar{\varphi}'(u_{\varepsilon }^{\gamma}-u_{\nu}^{\gamma})|
 \exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}]dx\,dt
 \\
&\quad+  \bar{\theta}\int_{Q}\left(\frac{1}{2}+|\nabla u_{\nu}|^2\right)
\exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}]dx\,dt
 \\
&\quad +C_0\sum_{\gamma=1}^2\int_{Q}\bar{\beta}_{\gamma}|\bar{\varphi}'
 (u_{\varepsilon }^{\gamma}-u_{\nu}^{\gamma})
 |\exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}]dx\,dt
 \\
&\quad +  2C_2\int_{Q}|\nabla u_{\nu}^1|^2\bar{\beta}_1
 |\bar{\varphi}'(u_{\varepsilon }^1
-u_{\nu}^1)| \exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}]dx\,dt
 \\
 &\quad +  2\eta\int_{Q}|\nabla u_{\nu}^2|^2\bar{\beta}_1
 |\bar{\varphi}'(u_{\varepsilon }^1-u_{\nu}^1)|
 \exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}]dx\,dt
 \\
 &\quad + 2C_2\int_{Q}|\nabla u_{\nu}^1|^2\bar{\beta}_2|\bar{\varphi}'
 (u_{\varepsilon }^2-u_{\nu}^2)| \exp[\bar{\mu}{\bar{\psi}
 ({u}_{\varepsilon }-{u}_{\nu})}]dx\,dt
 \\
 &\quad +  2C_2\int_{Q}|\nabla u_{\nu}^2|^2\bar{\beta}_2|
 \bar{\varphi}'(u_{\varepsilon }^2-u_{\nu}^2)|
  \exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}]dx\,dt
 \to  I_3^1+I_2^1+I_1^1+\sum_{i=2}^{5} I_{12}^{i},
 \end{align*}
 where
 \begin{align*}
&I_3^1+I_2^1+I_1^1+\sum_{i=2}^{5} I_{12}^{i} \\
&=  C_3\sum_{\gamma=1}^2\int_{Q}|\nabla u_{\nu}^{\gamma}
 |\bar{\beta}_{\gamma}|\bar{\varphi}'(u^{\gamma}-u_{\nu}^{\gamma})|
 \exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt
 \\
&\quad  +  \bar{\theta}\int_{Q}\left(\frac{1}{2}+|\nabla u_{\nu}|^2\right)
 \exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt
 \\
&\quad +C_0\sum_{\gamma=1}^2\int_{Q}\bar{\beta}_{\gamma}|\bar{\varphi}'
 (u^{\gamma}-u_{\nu}^{\gamma})|\exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt
 \\
&\quad +  2C_2\int_{Q}|\nabla u_{\nu}^1|^2\bar{\beta}_1|\bar{\varphi}'
 (u^1-u_{\nu}^1)| \exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt
 \\
&\quad +  2\eta\int_{Q}|\nabla u_{\nu}^2|^2\bar{\beta}_1|\bar{\varphi}'
 (u^1-u_{\nu}^1)| \exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt
 \\
&\quad +  2C_2\int_{Q}|\nabla u_{\nu}^1|^2\bar{\beta}_2|\bar{\varphi}'
 (u^2-u_{\nu}^2)| \exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt
 \\
&\quad +  2C_2\int_{Q}|\nabla u_{\nu}^2|^2\bar{\beta}_2
 |\bar{\varphi}'(u^2-u_{\nu}^2)| \exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt.
\end{align*}
 Now we pass to the limit in the remaining term of the right hand side
of \eqref{inqeps}. Due to the fact that $|\nabla u_{\varepsilon }|$ is
uniformly bounded in $L^2(Q)$, one can extract a subsequence which
is not relabled such that
 \begin{equation}\label{nablaueps}
 |\nabla u_{\varepsilon }|\rightharpoonup w \quad\text{weakly in }L^2(Q).
 \end{equation}
 Applying Lebesgue's dominated convergence theorem we obtain
 \begin{equation}\label{nablau}
 \begin{split}
& |\nabla u_{\nu}^{\gamma}|\bar{\beta}_{\gamma}|\bar{\varphi}'
 (u_{\varepsilon }^{\gamma}-u_{\nu}^{\gamma})|
 \exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}]\\
&\to
 |\nabla u_{\nu}^{\gamma}|\bar{\beta}_{\gamma}|\bar{\varphi}'
 (u^{\gamma}-u_{\nu}^{\gamma})|
 \exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]
 \quad  \text{strongly in }L^2(Q)\,.
 \end{split}
 \end{equation}
 \eqref{nablaueps} and \eqref{nablau} yield
 \begin{equation}
 I_3^2(\varepsilon )=C_4\sum_{\gamma=1}^2\int_{Q}|\nabla u_{\varepsilon }||\nabla u_{\nu}^{\gamma}|\bar{\beta}_{\gamma}|\bar{\varphi}'(u_{\varepsilon }^{\gamma}-u_{\nu}^{\gamma})|
 \exp[\bar{\mu}{\bar{\psi}({u}_{\varepsilon }-{u}_{\nu})}]dx\,dt\to  I_3^2,
 \end{equation}
 where
 \begin{equation}
 I_3^2=C_4\sum_{\gamma=1}^2\int_{Q}w|\nabla u_{\nu}^{\gamma}|\bar{\beta}_{\gamma}|\bar{\varphi}'(u^{\gamma}-u_{\nu}^{\gamma})|
 \exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt.
 \end{equation}

 Passing to the limit infimum in the inequality \eqref{inqeps} we obtain
 \begin{equation}\label{inqnu}
\begin{aligned}
&J_2+J_3+\int_{Q}B_{\nu}^1
 \nabla(u^1-u_{\nu}^1)\cdot\nabla(u^1-u_{\nu}^1)dx\,dt
 \\
&+\int_{Q}B_{\nu}^2
 \nabla(u^2-u_{\nu}^2)\cdot\nabla(u^2-u_{\nu}^2)dx\,dt+
 \int_{Q}E_{\nu}
 \nabla\bar{\psi}(u-u_{\nu})\cdot\nabla\bar{\psi}(u-u_{\nu})dx\,dt
 \\
&\leq I_1^1+I_{12}^2+I_{12}^{3}+I_{12}^{4}+
 I_{12}^{5}+I_3^1+I_3^2+I_2^1.
\end{aligned}
 \end{equation}
Multiplying the $\gamma$-th equation of  system \eqref{asystem}
relative to the parameter $\nu$ by
$\bar{\beta}_{\gamma}\bar{\varphi}'(u^{\gamma}-u_{\nu}^{\gamma})
\exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]$
 and sum  from $\gamma=1$ to $\gamma=2$ we obtain
 \begin{equation}\label{jequality}
 \begin{aligned}
 J_3&=\sum_{\gamma=1}^2\int_{Q}
 A(u_{\nu})\nabla(u^{\gamma}-u_{\nu}^{\gamma})\cdot\nabla
\left[\bar{\beta}_{\gamma}  \bar{\varphi}'(u^{\gamma}-u_{\nu}^{\gamma})
 \exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]\right]dx\,dt
 \\
 &\quad -\sum_{\gamma=1}^2\int_{Q}
 A(u_{\nu})\nabla u^{\gamma}\cdot\nabla\left[\bar{\beta}_{\gamma}
 \bar{\varphi}'(u^{\gamma}-u_{\nu}^{\gamma})
\exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]\right]dx\,dt
 \\
 &\quad +\sum_{\gamma=1}^2\int_{Q}G_{\nu}^{\gamma}(u_{\nu},
 \nabla u_{\nu})\bar{\beta}_{\gamma}\bar{\varphi}'(u^{\gamma}
 -u_{\nu}^{\gamma})\exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt
 \\
 &\quad -\int_{Q} F_{\nu}(u_{\nu},\nabla
 u_{\nu})\cdot\nabla{\bar{\psi}(u-u_{\nu})}
 \exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt
 \\
 &\quad +\sum_{\gamma=1}^2\int_{Q} F_{\nu}(u_{\nu},\nabla
 u_{\nu})\cdot\nabla u^{\gamma}\bar{\beta}_{\gamma}
 \bar{\varphi}'(u^{\gamma}-u_{\nu}^{\gamma})
 \exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt
 \\
 &=K_1+K_2+K_3+K_4+K_5.
 \end{aligned}
 \end{equation}
 Using \eqref{jequality} and replacing the values of the matrix
$B_{\nu}^1$, $B_{\nu}^2$ and $E_{\nu}$ in \eqref{inqnu} we have
 \begin{equation}\label{remplacement}
 J_2+K_1+S_1+S_2+S_3+S_4
\leq I_1^1+ I_3^1+ I_3^2+I_2^1 +\sum_{i=1}^{5}I_{12}^{i}-\sum_{i=2}^{5}K_i
 \end{equation}
 where
 \begin{gather*}
S_1=-\int_{Q}L_1(u,u_{\nu})|\nabla(u^1-u_{\nu}^1)|^2\exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt,
 \\
 S_2=-\int_{Q}L_2(u,u_{\nu})|\nabla(u^2-u_{\nu}^2)|^2\exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt,
 \\
 S_3=-\Big(\frac{C_3^2}{2\bar{\theta}}+\frac{C_4^2}{2\bar{\theta}}\Big)
 \int_{Q}|\nabla{\bar{\psi}(u-u_{\nu})}|^2
 \exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt,
 \\
 \begin{aligned}
 S_4&=\int_{Q}A(u)\nabla(u^1-u_{\nu}^1)\cdot\nabla(u-u_{\nu}^1)
 \bar{\beta}_1\bar{\varphi}''(u^1-u_{\nu}^1)
 \exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt
 \\
 &\quad +\int_{Q}A(u)\nabla(u^2-u_{\nu}^2)\cdot\nabla(u-u_{\nu}^1)
 \bar{\beta}_2\bar{\varphi}''(u^2-u_{\nu}^2)
 \exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt
 \\
 &\quad +\bar{\mu}\int_{Q}A(u)\nabla{\bar{\psi}(u-u_{\nu})}
 \cdot\nabla{\bar{\psi}(u-u_{\nu})}\exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt
 \\
 &=\sum_{\gamma=1}^2\int_{Q}
 A(u)\nabla(u^{\gamma}-u^{\gamma}_{\nu})\cdot\nabla\left[\bar{\beta}_{\gamma}
 \bar{\varphi}'(u^{\gamma}-u_{\nu}^{\gamma})
 \exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]\right]dx\,dt
 \\
&=\sum_{\gamma=1}^2\int_{Q}
 A(u)\nabla u ^{\gamma}\cdot\nabla\left[\bar{\beta}_{\gamma}
 \bar{\varphi}'(u^{\gamma}-u_{\nu}^{\gamma})\exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]\right]dx\,dt
 \\
&\quad -\sum_{\gamma=1}^2\int_{Q}
 A(u)\nabla u^{\gamma}_{\nu}\cdot\nabla\left[\bar{\beta}_{\gamma}
 \bar{\varphi}'(u^{\gamma}-u_{\nu}^{\gamma})
\exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]\right]dx\,dt \\
&=K_{6}-J_2.
 \end{aligned}
 \end{gather*}
Replacing the value of $S_4$ in \eqref{remplacement} we obtain
 \begin{equation}\label{remplacing}
 K_1+S_1+S_2+S_3
\leq I_1^1+I_3^1+ I_3^2+I_2^1+\sum_{i=1}^{5}I_{12}^{i}-\sum_{i=2}^{6}K_i.
 \end{equation}
By using the growth conditions \eqref{Fcroissance} and \eqref{G1G2Fcroissance}
on $F_{\nu}$ and Young's inequality we obtain
 \begin{equation}\label{K4}
 \begin{aligned}
 -K_4
&\leq\int_{Q}(C_3+C_4|\nabla u_{\nu}|)|\nabla{\bar{\psi}(u-u_{\nu})}|
 \exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt
 \\
&\leq\int_{Q}\Big[\frac{\bar{\theta}}{2}+\frac{C_3^2}{2\bar{\theta}}
 |\nabla{\bar{\psi}(u-u_{\nu})}|^2+
 \frac{\bar{\theta}}{2}|\nabla u_{\nu}|^2+\frac{C_4^2}{2\bar{\theta}}|
 \nabla{\bar{\psi}(u-u_{\nu})}|^2\Big] \\
&\quad\times  \exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt
 \\
&\leq-S_3+\bar{\theta}\int_{Q}|\nabla(u-u_{\nu})|^2
 \exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt+K_4^1,
 \end{aligned}
\end{equation}
 where
 \begin{equation}\label{K41}
 K_4^1=\int_{Q}\bar{\theta}big(\frac{1}{2}+|\nabla
 u|^2\big)\exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt.
 \end{equation}
 Employing the growth conditions \eqref{G1croissance}, \eqref{G2croissance}
and \eqref{G1G2Fcroissance} on $G_{\nu}^{\gamma}$, the term $-K_3$
can be estimated as
 \begin{equation}\label{K3}
 \begin{aligned}
 -K_3
&\leq C_0\sum_{\gamma=1}^2\int_{Q}\bar{\beta}_{\gamma}
 |\bar{\varphi}'(u^{\gamma}-u_{\nu}^{\gamma})|
 \exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt
 \\
&\quad +  C_2\int_{Q}|\nabla u_{\nu}^1|^2\bar{\beta}_1
 |\bar{\varphi}'(u^1-u_{\nu}^1)| \exp[\bar{\mu}{\bar{\psi}({u}
 -{u}_{\nu})}]dx\,dt
 \\
&\quad +  \eta\int_{Q}|\nabla u_{\nu}^2|^2\bar{\beta}_
 1|\bar{\varphi}'(u^1-u_{\nu}^1)| \exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt
 \\
 &\quad +  C_2\int_{Q}|\nabla u_{\nu}^1|^2\bar{\beta}_2
 |\bar{\varphi}'(u^2-u_{\nu}^2)| \exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt
 \\
&\quad +  C_2\int_{Q}|\nabla u_{\nu}^2|^2\bar{\beta}_2|\bar{\varphi}'
 (u^2-u_{\nu}^2)| \exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt
 \\
&\leq I_1^1+\frac{1}{2}\sum_{i=2}^{5}I_{12}^2.
 \end{aligned}
 \end{equation}
Using  Young's inequality the term $\sum_{i=2}^{5}I_{12}^2$ can be controlled by
 \begin{equation}
 \begin{aligned}
 \sum_{i=2}^{5}I_{12}^2
&\leq 4C_2\int_{Q}|\nabla(u^1-u_{\nu}^1)|^2\bar{\beta}_1|\bar{\varphi}'
 (u^1-u_{\nu}^1)| \exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt
 \\
&\quad +4\eta\int_{Q}|\nabla(u^2-u_{\nu}^2)|^2\bar{\beta}_1|\bar{\varphi}'
 (u^1-u_{\nu}^1)| \exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt
 \\
&\quad +4C_2\int_{Q}|\nabla(u^1-u_{\nu}^1)|^2\bar{\beta}_2
|\bar{\varphi}'(u^2-u_{\nu}^2)| \exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt
 \\
&\quad + 4C_2\int_{Q}|\nabla(u^2-u_{\nu}^2)|^2\bar{\beta}_2
 |\bar{\varphi}'(u^2-u_{\nu}^2)| \exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt
 +\sum_{i=2}^{5}R_{12}^2,
 \end{aligned}
\end{equation}
where
 \begin{equation}
 \begin{aligned}
 \sum_{i=2}^{5}R_{12}^2
&=4C_2\int_{Q}|\nabla u^1|^2\bar{\beta}_1|\bar{\varphi}'
 (u^1-u_{\nu}^1)| \exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt
 \\
&\quad +4\eta\int_{Q}|\nabla u^2|^2\bar{\beta}_1|\bar{\varphi}'
 (u^1-u_{\nu}^1)| \exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt
 \\
&\quad +4C_2\int_{Q}|\nabla u^1|^2\bar{\beta}_2|\bar{\varphi}'
 (u^2-u_{\nu}^2)| \exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt
 \\
&\quad +4C_2\int_{Q}|\nabla u^2|^2\bar{\beta}_2|\bar{\varphi}'
 (u^2-u_{\nu}^2)| \exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt.
\end{aligned}
 \end{equation}
We use again Young's inequality to obtain
 \begin{equation}
\begin{aligned}
 I_2^1&\leq 2\bar{\theta}\int_{Q}|\nabla(u-u_{\nu})|^2
 \exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt\\
&\quad +\underbrace{\bar{\theta}\int_{Q}\left(\frac{1}{2}
 +2|\nabla u|^2\right)\exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt}_{=R_2^1}.
\end{aligned}
 \end{equation}
 We employ the coercivity condition \eqref{ellipticite} to estimate the
term $K_1$ from below by
 \begin{equation}\label{K1}
 \begin{aligned}
 K_1&\geq \alpha\int_{Q}|\nabla(u^1-u_{\nu}^1)|^2\bar{\beta}_1
 \bar{\varphi}''(u^1-u_{\nu}^1)
 \exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt
 \\
 &\quad +\alpha\int_{Q}|\nabla(u^2-u_{\nu}^2)|^2
 \bar{\beta}_2\bar{\varphi}''(u^2-u_{\nu}^2)
 \exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt
 \\
&\quad +\alpha\bar{\mu}\int_{Q}|\nabla\bar{\psi}|^2
 \exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt.
 \end{aligned}
 \end{equation}
Using \eqref{K4}--\eqref{K1}, inequality \eqref{remplacing} becomes
 \begin{equation}\label{soui}
 \begin{aligned}
&\sum_{\gamma=1}^2\int_{Q}|\nabla(u^{\gamma}-u_{\nu}^{\gamma})|^2\left\lbrace\alpha\bar{\beta}_{\gamma}\bar{\varphi}''(u^{\gamma}-u_{\nu}^{\gamma})
 -4L_{\gamma}(u,u_{\nu})\right\rbrace\exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt
 \\
&+\int_{Q}|\nabla\bar{\psi}|^2\left[\alpha\bar{\mu}-\left(\frac{C_3^2}{\bar{\theta}}+\frac{C_4^2}{\bar{\theta}}\right)\right]\exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt
 \\
&\leq 2I_1^1+I_3^1+I_3^2+R_2^1-K_2+K_4^1-K_5-K_{6}+\frac{3}{2}
 \sum_{i=2}^{5}R_{12}^{i}\\
&=R_{\nu}+R_2^1+K_4^1.
 \end{aligned}
 \end{equation}
Using \eqref{rho1bar}, \eqref{rho2bar} and the definition of $\bar{\mu}$
in \eqref{bardata} we obtain from \eqref{soui} that
 \begin{equation}
 \begin{aligned}
&\bar{\alpha}_0\sum_{\gamma=1}^2\int_{Q}|\nabla(u^{\gamma}-u_{\nu}^{\gamma})|^2
\exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt \\
&+\underbrace{\frac{\alpha\bar{\mu}}{2}\int_{Q}|\nabla\bar{\psi}|^2
\exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt}_{ \geq0}
 \\
&\leq R_{\nu}+\bar{\theta}\int_{Q}\left(1+3|\nabla u|^2\right)
\exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt,
 \end{aligned}
 \end{equation}
hence
 \begin{equation}\label{inqfin}
 \begin{aligned}
&\bar{\alpha}_0\sum_{\gamma=1}^2\int_{Q}|\nabla(u^{\gamma}-u_{\nu}^{\gamma})|^2
\exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt\\
&\leq R_{\nu} + \bar{\theta}\int_{Q}\left(1+3|\nabla u|^2\right)
\exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt,
 \end{aligned}
 \end{equation}
 such that
\begin{equation}\label{Rnu}
 R_{\nu}=2I_1^1+I_3^1+I_3^2-K_5-\left(K_2+K_{6}\right)+\frac{3}{2}
 \sum_{i=2}^{5}R_{12}^{i},
 \end{equation}
where
 \begin{equation}\label{part1ofRnu}
 \begin{aligned}
&2I_1^1+I_3^1+I_3^2-K_5\\
&=2C_0\sum_{\gamma=1}^2\int_{Q}\bar{\beta}_{\gamma}|\bar{\varphi}'(u^{\gamma}-u_{\nu}^{\gamma})|^2 \exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt
 \\
&\quad + C_3\sum_{\gamma=1}^2\int_{Q}|\nabla u_{\nu}^{\gamma}
 |\bar{\beta}_{\gamma}|\bar{\varphi}'(u^{\gamma}-u_{\nu}^{\gamma})|
 \exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt
 \\
&\quad +C_4\sum_{\gamma=1}^2\int_{Q}w|\nabla u_{\nu}^{\gamma}
 |\bar{\beta}_{\gamma}|\bar{\varphi}'(u^{\gamma}-u_{\nu}^{\gamma})|
 \exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt
 \\
&\quad -\sum_{\gamma=1}^2\int_{Q} F_{\nu}(u_{\nu},\nabla
 u_{\nu})\cdot\nabla u^{\gamma}\bar{\beta}_{\gamma}
 \bar{\varphi}'(u^{\gamma}-u_{\nu}^{\gamma})\exp[\bar{\mu}{\bar{\psi}
 (u-u_{\nu})}]dx\,dt,
\end{aligned}
 \end{equation}
such that
 \begin{equation}\label{part2ofRnu}
\begin{aligned}
&-\left(K_2+K_{6}\right)\\
&=\sum_{\gamma=1}^2\int_{Q}
 \left[A(u_{\nu})-A(u)\right]\cdot\nabla u^{\gamma}\nabla
\left[\bar{\beta}_{\gamma}
 \bar{\varphi}'(u^{\gamma}-u_{\nu}^{\gamma})
\exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]\right]dx\,dt,
\end{aligned}
 \end{equation}
 and
 \begin{equation}\label{part3ofRnu}
 \begin{aligned}
 \frac{3}{2}\sum_{i=2}^{5}R_{12}^{i}
&=6C_2\int_{Q}|\nabla u^1|^2\bar{\beta}_1|\bar{\varphi}'(u^1-u_{\nu}^1)|
\exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt
 \\
&\quad +6\eta\int_{Q}|\nabla u^2|^2\bar{\beta}_1|\bar{\varphi}'
(u^1-u_{\nu}^1)| \exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt
 \\
&\quad +6C_2\int_{Q}|\nabla u^1|^2\bar{\beta}_2|\bar{\varphi}'
(u^2-u_{\nu}^2)| \exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt
 \\
&\quad +6C_2\int_{Q}|\nabla u^2|^2\bar{\beta}_2|\bar{\varphi}'
(u^2-u_{\nu}^2)| \exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt.
 \end{aligned}
 \end{equation}
In view of
 \begin{equation}\label{converge}
 \begin{gathered}
 u_{\nu}\rightharpoonup
 u\quad\text{weakly in }(L^2(0,T;H_0^1(\Omega)))^{m},\\
 u_{\nu}\to  u\quad\text{strongly in }(L^2(Q))^{m},\\
 u_{\nu}\to  u\quad\text{a.e.in}\;Q,\\
 \end{gathered}
 \end{equation}
and the fact that $u-u_{\nu}$ is uniformly bounded in $L^{\infty}(Q)$,
for all $\gamma=1,2$, it holds
 \begin{gather}\label{strongly}
 \bar{\beta}_{\gamma}\bar{\varphi}'(u^{\gamma}-u_{\nu}^{\gamma})
\exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]\to  0
\quad\text{strongly in }L^2(Q), \\
\label{weaklyH10}
 \bar{\beta}_{\gamma}\bar{\varphi}'(u^{\gamma}-u_{\nu}^{\gamma})\exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]\rightharpoonup 0
\quad\text{weakly in }L^2(0,T;H_0^1(\Omega)),
 \end{gather}
this follows from the fact that
$\bar{\beta}_{\gamma}\bar{\varphi}'(u^{\gamma}-u_{\nu}^{\gamma})
\exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]$ remains bounded in
$L^2(0,T;H_0^1(\Omega))$ and converges
almost everywhere in $Q$ towards 0.
By \eqref{converge} and employing Lebesgue's dominated convergence theorem,
we obtain
 \begin{equation}
 \frac{3}{2}\sum_{i=2}^{5}R_{12}^{i}\to  0 \quad\text{as }\nu\to  0.
 \end{equation}
By \eqref{strongly}, Young's inquality and employing Lebesgue's dominated
convergence theorem, we obtain
 \begin{equation}
 2I_1^1+I_3^1+I_3^2\to  0 \quad\text{in }\nu\to  0,
 \end{equation}
 \eqref{converge} and Lebesgue's dominated convergence yield
 \begin{equation}
 \left[ A(u_{\nu})-A(u) \right]\nabla u^{\gamma}\to  0\quad\text{strongly in }
\left(L^2(Q)\right)^{N}.
 \end{equation}
Therefore
 \begin{equation}\label{div}
 -\operatorname{div} \left(\left[ A(u_{\nu})-A(u) \right]\nabla u^{\gamma} \right)\to
 0  \quad\text{strongly in }L^2(0,T;H^{-1}(\Omega)),
 \end{equation}
 \eqref{div} and \eqref{weaklyH10} imply that
 \begin{equation}\label{K2plusK6}
 \begin{split}
& -\left(K_2+K_{6} \right)\\
&= \sum_{\gamma=1}^2\langle -div \left(\left[ A(u_{\nu})-A(u) \right]
\nabla u^{\gamma} \right);
 \bar{\beta}_{\gamma}\bar{\varphi}'(u^{\gamma}
-u_{\nu}^{\gamma})\exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]
 \rangle
 \to  0.
 \end{split}
 \end{equation}
 For the rest terms of \eqref{part1ofRnu}, using the growth conditions
\eqref{Fcroissance} and \eqref{G1G2Fcroissance} on $F_{\nu}$, and
 the fact that $u_{\nu}$ is uniformly bounded in $(L^2(Q))^2$,
by applying the Young's inequality and afterwards Lebesgue's theorem one
can prove that
 \begin{equation}
 -K_5\to  0,\quad\text{when } \nu\to  0.
 \end{equation}
 Finally we have
 $R_{\nu}\to  0,\quad\text{when } \nu\to  0$,
 similarly we have
 \begin{equation}
 \bar{\theta}\int_{Q}\left(1+4|\nabla u|^2\right)
\exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]dx\,dt\to
 \bar{\theta}\int_{Q}\left(1+4|\nabla u|^2\right)dx\,dt.
 \end{equation}
Since $\exp[\bar{\mu}{\bar{\psi}(u-u_{\nu})}]\geq1$,
 we deduce from \eqref{inqfin} that
 \begin{equation}
 \limsup_{\varepsilon \to 0}\bar{\alpha}_0\sum_{\gamma=1}^2
\int_{Q}|\nabla(u^{\gamma}-u_{\nu}^{\gamma})|^2dx\,dt
\leq  \bar{\theta}\int_{Q}\left(1+4|\nabla u|^2\right)dx\,dt.
 \end{equation}
Since $\bar{\alpha}_0$
 is independent of\;$\bar{\theta}$ this implies that
 $u_{\nu}\to  u$ strongly in $(L^2(0,T;H_0^1(\Omega)))^2$.


 \subsection{Passing to the limit}

 As $u_{\varepsilon }$ converges strongly to $u$ in
$\left( L^2(0,T;H_0^1(\Omega))\right)^2$, for all $\gamma=1,2$ we have
 \begin{equation}
 G^{\gamma}(u_{\varepsilon },\nabla u_{\varepsilon })+F(u_{\varepsilon },
\nabla u_{\varepsilon })\cdot\nabla u_{\varepsilon }^{\gamma}\to
G^{\gamma}(u,\nabla u)+F(u,\nabla u)\cdot\nabla u^{\gamma}
\end{equation}
strongly in  $L^1(Q)$.
 Therefore, applying Vitali's theorem,
 \begin{equation}
 G^{\gamma}_{\varepsilon }(u_{\varepsilon },\nabla u_{\varepsilon })
+F_{\varepsilon }(u_{\varepsilon },\nabla u_{\varepsilon })
\cdot\nabla u_{\varepsilon }^{\gamma}\to  G^{\gamma}(u,\nabla u)
+F(u,\nabla u)\cdot\nabla u^{\gamma}
 \end{equation}
strongly in $L^1(Q)$.

 On the other hand,
 \begin{equation}
 -\operatorname{div}\big(A(u_{\varepsilon })\nabla u_{\varepsilon }^{\gamma}\big)
\to  -\operatorname{div}\big(A(u)\nabla u^{\gamma}\big)
 \end{equation}
strongly in $L^2(0,T;H^{-1}(\Omega))$.

 Thus $\frac{\partial u_{\varepsilon }^{\gamma}}{\partial t}$ converges strongly to
$\frac{\partial u^{\gamma}}{\partial t}$ in $L^2(0,T;H^{-1}(\Omega))+L^1(Q)$
 which implies
 \begin{equation*}
 \frac{\partial u^{\gamma}}{\partial t} \in L^2(0,T;H^{-1}(\Omega))+L^1(Q),
 \end{equation*}
 and
 \begin{gather*}
 \frac{\partial u^1}{\partial t} -\sum_{i,j=1}^{N}\frac{\partial}{\partial
 x_i}(a_{i,j}(u)\frac{\partial u^1}{\partial x_j}=G^1(u,\nabla u)+
 F(u,\nabla u)\cdot\nabla u^1\quad\text{in }Q, \\
 \frac{\partial u^2}{\partial t} -\sum_{i,j=1}^{N}\frac{\partial}{\partial
 x_i}(a_{i,j}(u)\frac{\partial u^2}{\partial x_j}=G^2(u,\nabla u)+
 F(u,\nabla u)\cdot\nabla u^2\quad\text{in }Q.
 \end{gather*}
For $s$ large enough,  for all ($ \gamma=1,2$), we have
$\frac{\partial u_{\varepsilon }^{\gamma}}{\partial t}$ converges strongly
to $\frac{\partial u^{\gamma}}{\partial t}$ in $L^1(0,T;H^{-s}(\Omega))$.
Then  $u_{\varepsilon }^{\gamma}$ converges strongly to $u^{\gamma}$ in
$\mathcal{C}([0,T];H^{-s}(\Omega))$, hence
 $u_{\varepsilon }^{\gamma}(0)$ converges strongly to $u^{\gamma}(0)$
in $H^{-s}(\Omega)$.
 Therefore, $u$ satisfies \eqref{system}.
 We thus have proved the existence of at least one solution of system
\eqref{system}. this solution $u$ is bounded because $u_{\varepsilon }$
is uniformly bounded in $\big( L^{\infty}(Q)\big)^2$ and $u_{\varepsilon }\to  u$
a.e. in $Q$. This completes the proof.
\end{proof}

\begin{thebibliography}{99}

 \bibitem{Abudiab} M. Abudiab, I. Ahn, L. Li;
\emph{A theorem on upper-lower solutions for nonlinear elliptic systems
 and its applications}, J. Math. Anal. Appl., 340 (2008), 175-182.

 \bibitem{Abudiab2} M. Abudiab, I. Ahn, L. Li;
\emph{Upper-lower solutions for nonlinear parabolic systems and their applications},
J. Math. Anal. Appl., 378 (2011), 620-633.

 \bibitem{Aris} R. Aris, K. Zygourakis;
\emph{Weakly coupled systems of nonlinear elliptic boundary
 value problems}, Nonlinear Analysis, Theory Methods, and Appi., 6 (1982), 555-569.

\bibitem{Ben} M. Bendahmane,  K. H. Karlsen;
\emph{Nonlinear anisotropic elliptic and parabolic equations with variable
exponents and $L^1$-data}. Commun. Pure Appl. Anal., 12 (2013) 1201-1220.
,
\bibitem{Ben1} M. Bendahmane,  K. H. Karlsen;
\emph{Renormalized entropy solutions for quasi-linear anisotropic degenerate
 parabolic equations}, SIAM journal on mathematical analysis, 36(2) (2004), 405-422.

\bibitem{Ben2} M. Bendahmane, M. Langlais, M. Saad;
\emph{On some anisotropic reaction-diffusion systems with L 1-data modeling
the propagation of an epidemic disease}. Nonlinear Analysis, Theory, Methods
and Applications, 54(4) (2003), 617-636.

\bibitem{Ben3} M. Bendahmane,  K. H. Karlsen;
\emph{Analysis of a class of degenerate reaction-diffusion
systems and the bidomain model of cardiac tissue}, Netw. Heterog.
Media., 1(1) (2006) 185-218.

\bibitem{BM} L. Boccardo, F. Murat, J. P. Puel;
\emph{Existence results for some quasilinear
 parabolic equations}, Non Lin. Anal. TMA, 13 (1988) 373-392.

 \bibitem{Cosner} C. Cosner, F. Schindler;
Upper and lower solutions for systems of second order equations with nonnegative
 characteristic form and discontinuous nonlinearities,
Rocky Mt. J. Math. 14 (1984) 549-557.

\bibitem{DK} Z. Dahmani, S. Kerbal;
\emph{Solution to nonlinear gradient dependent systems with a balance law, Electron. J. Differ.
Equ., 158 (2007) 1–7 .}

\bibitem{Dickstein} F. Dickstein, M. Escobedo;
 A maximum principle for semilinear parabolic systems and applications,
 Nonlinear Anal., 45 (2001) 825-837.

 \bibitem{JF1} J. Frehse;
\emph{Existence and perturbation theorems for nonlinear elliptic systems}. Nonlinear
Partial Diferential Equations and Their Applications, Coll\`{e}ge de France Seminar,
Vol. 4, Research Notes Math., 84 (1983) 87-110.

 \bibitem{JF3} J. Frehse, A. Bensoussan, M.Bulucek;
\emph{Existence and compactness for weak solutions to Bellman systems with
 critical growth, Discrete Contin}. Dyn. Syst. Ser. B, 17(6) (2012) 1729-1750.

 \bibitem{JF8} J. Frehse, A. Bensoussan;
\emph{Smooth solutions of systems of quasilinear parabolic equations},
ESAIM: Control, Optimisation and Calculus of Variations, 8 (2002) 169-193.

\bibitem{Gia3} M. Giaquinta;
\emph{Multiple Integrals in the Calculus of Variations and Nonlinear Elliptic
Systems}. Annals of Mathematics Studies 105. Princeton University Press, Princeton,
1983.

 \bibitem{Gia} M. Giaquinta, J. Necas;
\emph{On the regularity of weak solutions to non
 linear elliptic systems of partial differential equations}, Journal f\"ur die
reine und angewandte Mathematik, 316 (1980), 140-159.

 \bibitem{Gia2} M. Giaquinta, M. Struwe;
\emph{On the partial regularity of weak solutions of nonlinear parabolic systems},
Math. Z., 179 (1982) 437-451.

 \bibitem{Hil1} S. Hildebrandt, K. O. Widman;
\emph{Some regularity results for quasilinear elliptic systems of second order}.
 Math. Z., 142 (1975), 67-86.
Non. Lin. Ana. Great Britain. A 127 (1989), 373-392.

 \bibitem{Hil82} S. Hildebrandt;
\emph{Nonlinear elliptic systems and harmonic mappings}, Proc,
Beijing Symposium on diff. geom. and diff. Equ., Beijing (1980).

 \bibitem{Hil4} S. Hildebrandt;
\emph{Quasilinear elliptic systems in diagonal form}, in Systems of nonlinear
partial differential equations, 173-217. Editor: J.M. Ball. D. Reidel
Publishing Company 1983.

 \bibitem{Hil5} S. Hildebrandt, K. O. Widman;
\emph{On the H\"{o}lder continuity of weak solutions of quasilinear
elliptic systems of  second order}, Ann. Sc. Norm. Sup. Pisa, 4 (1977), 145-178.

\bibitem{Ki} M. Kirane;
\emph{Global bounds and asymptotic for a system of reaction-diffusion equations},
J. of Math. Anal. and Appl.,138 (1989), 328-342.

\bibitem{KIKo} M. Kirane, S. Kouachi;
\emph{Global solutions to a system of strongly coupled reaction-diffusion equations},
 Nonlin. Anal. Theo. Methods and Appl., 26 (1996), 1387-1396.

 \bibitem{Leung} A. Leung;
\emph{System of Nonlinear Partial Differential Equations: Applications to
Biology and Engineering}, Kluwer Acad. Publ., Norwell, MA,
 1989.

 \bibitem{Leung2} A. Leung;
\emph{Nonlinear systems of partial differential equations}.
 Applications to life and physical sciences.
 World Scientific Publishing Co. Pte. Ltd., Hackensack, NJ, 2009.

 \bibitem{Liang} J. Liang;
\emph{The reaction-diffusion system without quasi-monotone conditions},
Acta Sci. Natur. Univ. Pekinensis, 3 (1987), 9-20.

 \bibitem{L} J. L. Lions;
\emph{Quelques m\'ethodes de resolution des  probl\`emes aux limites non
 lin\'eaires}, Dunod, Gauthier-Villars, Paris (1969).

 \bibitem{Sleeman} G. Lu, B.D. Sleeman;
\emph{Maximum principles and comparison theorems for semi-
 linear parabolic systems and their applications}, Proc. Roy. Soc. Edinburgh Sect. A
 123 (1993) 857-885.

 \bibitem{Chen} R. Ma, R. Chen, Y. Lu;
\emph{Method of Lower and Upper Solutions for Elliptic Systems with Nonlinear
Boundary Condition and Its Applications}, Journal of Applied Mathematics,
vol. 2014 (2014).

 \bibitem{Mckenna} P.J. Mckenna, W. Walter;
\emph{On the Dirrichlet problem for elliptic systems},
Appl. Anal., 21 (1986) 201-224.

 \bibitem{AM2} A. Mokrane;
\emph{Existence for quasilinear elliptic
 system with quadratic growth having a particular structure},
 Portugaliae Mathematica, 55. Fasc. 1 (1998).

 \bibitem{Pao2} C. V. Pao;
\emph{Nonlinear Parabolic and Elliptic Equations}, Plenum Press, 1992.

 \bibitem{Per} I. Per-Anders;
\emph{On Quasilinear Elliptic Systems of Diagonal Form},
Mathematische Zeitschrift, 170(3) (1980), 283-286.

 \bibitem{Portnyagin} D. Portnyagin;
\emph{A generalization of the maximum principle to nonlinear parabolic systems},
Annales Pol. Math., 81(3) (2003), 217-236.

\bibitem{Sa} M. Saad;
\emph{A unilateral problem for elliptic equations with quadratic growth
and for $L^1-$ data}, Nonlinear Analysis, 44 (2001), 217-238.

 \bibitem{Si} J. Simon;
\emph{Compact sets in the space $L^p(0,T;B)$},
 Ann. Mat. Pura Appl, 146 (1987), 65-96.

 \bibitem{Sindler} F. Sindler;
\emph{Nonlinear elliptic and parabolic systems of second order with upper
and lower solutions}, Dissertation, University of Utah, 179 (4) (1982), 437-451.

 \bibitem{Struwe} M. Struwe;
\emph{On the H\"{o}lder continuity of bounded weak solutions of quasilinear
 parabolic systems}. Manuscripta Math., 35 (1981) 125-145.

 \bibitem{Temam} R. Temam;
\emph{Navier-Stokes Equations Theory and Numerical Analysis}, North-Holland,
 Amsterdam (1977).

 \bibitem{Wid} K. O. Widman;
\emph{H\"{o}lder continuity of solutions of elliptic systems}, Manuscripta
 Math., 5 (1971), 299-308.

 \bibitem{Wie} M. Wiegner;
\emph{Ein optimaler Regul\"{a}rittssatz f\"{u}r schwache L\"{o}sungen
gewisser elliptischer Systeme}. Math. Z., 147 (1976),  21-28.

 \bibitem{Wie2} M. Wiegner;
\emph{das Existenz und Regularit\"{a}tsproblem bei Systemen nichtlinearer
 elliptischer Differentialgleichungen. phdthesis}, University of Bochum, 1977.

 \bibitem{JZ} Z. Wu, J. Yin, C. Wang;
\emph{Elliptic and Parabolic Equation}, World Scientific Publishing CO.,
Singapore, 2006.

\end{thebibliography}
\end{document}