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

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

\begin{document}
\title[\hfilneg EJDE-2018/200\hfil Isentropic quantum drift-diffusion model]
{Solution to a multi-dimensional isentropic
 quantum drift-diffusion model for bipolar semiconductors}

\author[J. Ri, S. Ra \hfil EJDE-2018/200\hfilneg]
{Jinmyong Ri, Sungjin Ra}

\address{Jinmyong Ri (corresponding author)\newline
Institute of Mathematics,
State Academy of Sciences,
Pyongyang, Korea}
\email{jmri2015@163.com}

\address{Sungjin Ra \newline
Department of Mathematics,
University of Science,
Pyongyang, Korea}
\email{math.inst@star-co.net.kp}


\dedicatory{Communicated by Jesus Ildefonso Diaz}

\thanks{Submitted January 12, 2016. Published December 21, 2018.}
\subjclass[2010]{35A01, 35D30, 35J25, 35K35}
\keywords{Quantum drift-diffusion;  bipolar semiconductor; time-discretization; 
\hfill\break\indent mixed boundary value problem; semiclassical limit}

\begin{abstract}
 We study the existence of weak solution and semiclassical limit for mixed
 Dirichlet-Neumann boundary value problem of
 1,2,3-dimensional isentropic transient quantum
 drift-diffusion models for  bipolar semiconductors.
 A time-discrete approximate scheme for the model constructed employing
 the quantum quasi-Fermi potential is composed of non-degenerate
 elliptic systems, and the system in each time step has a solution in
 which the components of carrier's densities are strictly positive.
 Some stability estimates guarantee convergence of the approximate
 solutions and performance of the semiclassical limit.
\end{abstract}

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

\section{Introduction}

In this article we consider the bipolar isentropic quantum
drift-diffusion model
\begin{equation}\label{1.1}
\begin{gathered}
\frac{\partial n}{\partial t}
=\operatorname{div}\Big(-\varepsilon^2n\nabla\Big(\frac{\Delta\sqrt{n}}{\sqrt{n}}\Big)
+\theta\nabla n^{r}-n\nabla V\Big), \\
\frac{\partial p}{\partial t}
=\operatorname{div}\Big(-\varepsilon^2p\nabla\Big(\frac{\Delta\sqrt{p}}{\sqrt{p}}\Big)
+\theta\nabla p^{r}+p\nabla V\Big), \\
\lambda^2\Delta V=n-p-f    \quad  \text{in } \Omega\times(0,T),
\end{gathered}
\end{equation}
where $\Omega$ is a bounded domain of $R^{d}$ $(d=1,2,3)$ occupied by
semiconductor,  $n, p$ are the electron and hole's densities,
 $V$ is the electrostatic potential and $f(x)$ is
doping profile which describes the fixed background charges. The
parameters $\varepsilon>0$, $\lambda>0$, $\theta>0$ are the scaled
Planck constant, permittivity and temperature, respectively. $r>1$
is a constant. The model is the relaxation time limit version of the
quantum hydrodynamic model which is derived from the mixed state
Schr\"{o}dinger-Poisson system or, equivalently, the Wigner-Poisson
system.

The system is supplemented with the following initial and mixed
Dirichlet-Neumann boundary conditions which are physically motivated
and commonly employed in the quantum semiconductor modeling
\begin{equation}\label{1.2}
\begin{gathered}
(n, p, V)=(n_D,p_D,V_D),\quad
\varepsilon^2\Delta\sqrt{n}=\varepsilon^2\Delta\sqrt{p}=0\quad
\text{on }\Gamma_D,\\
\nabla n^{r}\cdot\nu=\nabla p^{r}\cdot\nu=\nabla V\cdot\nu=
\varepsilon^2\nabla(\frac{\Delta
\sqrt{n}}{\sqrt{n}})\cdot\nu=
\varepsilon^2\nabla(\frac{\Delta\sqrt{p}}{\sqrt{p}})\cdot\nu=0
\\ \text{on } \Gamma_{N}, \\
n(x,0)=n_0(x),\quad  p(x,0)=p_0(x) \quad \text{in }\Omega.
\end{gathered}
\end{equation}
The boundary $\partial\Omega\in C^{0,1}$ is piecewise regular  and
splits into two disjoint parts $\Gamma_D$(Ohmic contacts),
$\Gamma_{N}$(insulating parts).   $\nu$ denotes the unit outward
normal vector on $\partial\Omega$.  Putting $\varepsilon=0$ formally
in \eqref{1.1}, \eqref{1.2}, we obtain mixed boundary value problem
of classical drift-diffusion model. The geometrical figure and
boundary conditions affect greatly the flow of carriers because the
device becomes more and more smaller nowadays. Hence, the
investigation of the multi-dimensional mixed boundary value problem
is very important for engineers, especially.

From the view point of mathematics the model is a fourth-order
parabolic system for carrier's densities $n, p$ coupled with the
Poisson equation for electrostatic potential   $V$ and the main
difficulty in the investigation is that maximum principle is not
available in general for fourth-order parabolic equation to ensure
the positivity of the carrier's densities. Another difficulty is
that one could not expect, in general, smooth solution due to the
mixed boundary condition in \eqref{1.2} and singular behavior of the
solution may occur near
$\overline{\Gamma}_D\cap\overline{\Gamma}_{N}$ no matter how
smooth the known data are.


 For the 1-dimensional problem, many  results  were
obtained. In \cite{cj,ccj,ccj1,jp2,jv} the existence of weak solution
and semiclassical or quasineutral limit were studied for various boundary
value problems.
In \cite{nss} the existence of a unique strong solution near the
stationary solution and classical limit were studied for the
Dirichlet boundary value problem when the scaled Planck constant and
disturbance of boundary data are small.  In \cite{lzzh} it was
proved the unique existence of strong solution of initial value
problem which has the character of self-similarity in large time
when the doping profile is zero.


Also, for the multi-dimensional problem and related models, some results
were obtained only in the case of single boundary value problems.
The existence of strong solution for initial value problem near the
stationary solution was proved by the relaxation time limit argument
of quantum hydrodynamic model in \cite{jlm} and in \cite{c,cc} Neumann
or periodic boundary value problem was studied.
Concerning the zero-electric field and zero-temperature
approximation of the model we refer \cite{gst,gjt,jm}.

The investigation of stationary problem with both single and mixed
boundary value problems is very active and we refer
\cite{au,BCD-12,D-13,DW-15,j1,zlz}. However,
there is no any result for the multi-dimensional transient quantum
drift-diffusion model \eqref{1.1} with mixed boundary conditions
\eqref{1.2} which we are interested in.


In this article we prove the existence of  weak solution
and semiclassical limit for \eqref{1.1}, \eqref{1.2}. For this, we
construct a time-discrete approximate system for the original
transient system using the quantum quasi-Fermi potential. The
approximate system in each time step is non-degenerate elliptic
 and the approximate solution of the system exists.
Furthermore, the approximate carrier's densities in each time step
are strictly positive and some stability estimates needed for
convergence of the approximate solutions and semiclassical limit
hold.

 Throughout this paper we assume the following:
$\Gamma_D$ is nonempty open subset of $\partial\Omega$,
$\Gamma_{N}=\partial\Omega\backslash\overline{\Gamma}_D$.
$f\in L_{\infty}(\Omega)$ and $r>1$. $V_D$ is a trace of some function
$\widetilde{V}_D\in W^{1}_{\infty}(\Omega)$.
$n_0, p_0\in L_{\infty}(\Omega)\cap  H^{1}(\Omega)$,
$\exists m_0>0; n_0, p_0\geq{m_0},
 n_0|_{\Gamma_D}=n_D, p_0|_{\Gamma_D}=p_D$.
$n_D,  p_D>0$ are constants. With
 out loss of generality, we assume $n_D, p_D<1$.
In fact, if  $ n_D, p_D<k$, $k>1$, then the new functions
$\zeta=n/k$,  $\xi=p/k$ satisfy
  \eqref{1.1}, \eqref{1.2} with $\lambda^2/k$, $\theta k^{r-1}$, $f/k$
instead of $\lambda^2$,$\theta$, $f$,
 and with $n_D/k$,  $p_D/k$,  $n_0/k$,  $p_0/k$
instead of  $n_D$,  $p_D$,  $n_0$,  $p_0$.

The quantum quasi-Fermi potentials  $F,  G$ for
 \eqref{1.1} are defined as
 \begin{equation}
F=-\varepsilon^2\frac{\Delta\sqrt{n}}{\sqrt{n}}+\theta h(n)-V, \quad
G=-\varepsilon^2\frac{\Delta\sqrt{p}}{\sqrt{p}}+\theta h(p)+V.
\label{1.3}
\end{equation}
where $h(x)=\frac{r}{r-1}(x^{r-1}-1)$.
 For the time-discretization of \eqref{1.1}, we divide the time interval $(0,T]$
into $N$  subintervals $(t_{i-1}, t_i]$,  $i=1, 2,\dots, N$
 with mesh size $\tau=t_i-t_{i-1}=\frac{T}{N}$ and $t_0=0$.
 Given $\rho_{i-1}$,  $\eta_{i-1}$,  $i=1, 2,\dots , N$, we solve the
 following elliptic system recursively
 \begin{equation}
 \begin{gathered}
 \frac{1}{\tau}((\rho_i+a(\tau))^2-(\rho_{i-1}+a(\tau))^2)
 =\operatorname{div}((\rho_i+a(\tau))^2\nabla F_i),\\
\varepsilon^2\Delta \rho_i=\rho_i(\theta h(\rho^2_i)
 +\tau\ln \rho^2_i-F_i-V_i),\\
 \frac{1}{\tau}((\eta_i+a(\tau))^2-(\eta_{i-1}+a(\tau))^2)
 =\operatorname{div}((\eta_i+a(\tau))^2\nabla G_i),\\
\varepsilon^2\Delta \eta_i=\eta_i(\theta h(\eta^2_i)
 +\tau\ln \eta^2_i-G_i+V_i),\\
\lambda^2\Delta V_i=(\rho_i+a(\tau))^2-(\eta_i+a(\tau))^2-f\quad
\text{in } \Omega,\\
(\rho_i, \eta_i, F_i, G_i, V_i)=(\rho_D, \eta_D, F_D, G_D, V_D)\quad
\text{on } \Gamma_D,\\
\nabla\rho_i\cdot\nu=\nabla\eta_i\cdot\nu=\nabla
F_i\cdot\nu=\nabla G_i\cdot\nu=\nabla V_i\cdot\nu=0 \quad
\text{on } \Gamma_{N}.
\end{gathered}\label{1.4}
\end{equation}
where $(\rho_0,\eta_0)=(\sqrt{n_0},\sqrt{p_0})$,
$(\rho_D,\eta_D)=(\sqrt{n_D},\sqrt{p_D})$,
$F_D=\theta h(n_D)+\tau\ln n_D-V_D$,
$G_D=\theta h(p_D)+\tau\ln p_D+V_D$ and continuous function
$a(\tau)$ satisfies $a(\tau)>0$, there exists
$c>0$ such that $\frac{a(\tau)}{\tau^2}\leq c$,for all $\tau>0$ and
$a(0)=0$. We define the approximate solutions for \eqref{1.1},
\eqref{1.2} as follows
\begin{equation}
\big(\rho^{(N)},
\eta^{(N)}, F^{(N)}, G^{(N)}, V^{(N)}\big)(x,t)
=(\rho_i,\eta_i,  F_i, G_i, V_i), \quad  t\in (t_{i-1},  t_i].
\label{1.5}
\end{equation}
where $(\rho_i,  \eta_i,  F_i, G_i,  V_i)$ is the
solution to \eqref{1.4}.  Unlike the previous works
(see \cite{cj,jp2,jv}) where the embedding
$H^{1}(\Omega)\hookrightarrow L_{\infty}(\Omega)$ in 1-dimensional
case and exponential transformation are used essentially, we
introduce a new ``relaxation parameter'' $a(\tau)$ in the
semi-discretization to ensure non-degeneracy of the first and third
equations. The appearance of
$\tau\ln \rho_i^2, \tau\ln \eta_i^2$ in the
second and forth equations gives the positive lower bounds of
$\rho_i, \eta_i$ in each steps. So, we can prove the existence
of the weak solution to \eqref{1.4}  using the Stampacchia's
truncation method and Leray-Schauder fixed point theorem.
(Theorem \ref{thm1.1})

\begin{theorem}\label{thm1.1}
Let $N=N_0, N_0+1,\dots$ be the integers such
that
\begin{equation}\label{1.6}
\frac{T}{2\theta N}\|\nabla
V_D\|^2_{L_{\infty}(\Omega)}\leq\min \{-h(n_D), -h(p_D)\}.
\end{equation}
Then there exist weak solutions
$(\rho_i,\eta_i,F_i,G_i,V_i)\in(L_{\infty}(\Omega)\cap
H^{1}(\Omega))^{5}$, $i=1,2,\dots N$ to the recurrent elliptic system
\eqref{1.4} satisfying
\[
m_{i,N}\leq\rho_i, \eta_i\leq M_{i,N}
\]
for some $m_{i,N}, M_{i,N}>0$.
\end{theorem}

The entropy inequality in the previous works (see \cite{jp1,jp2})
 which show boundedness of the first-order derivatives of
the approximate solutions for the case of 1-dimensional model also
holds for our case (Lemma \ref{lem3.2}). This is enough for the upper bound of
the approximate carrier's densities independent of the mesh size in
1-dimensional case because of the embedding
$H^{1}(\Omega)\hookrightarrow L_{\infty}(\Omega)$. Furthermore, the
upper bound gives  boundedness
 of their second-order derivatives needed for  convergence
 of the scheme. Hence, for 1-dimensional problems the embedding
 $H^{1}(\Omega)\hookrightarrow L_{\infty}(\Omega)$ plays a
crucial role in the convergence of the approximate solution as well
as in the existence of the approximate solution. However, such
embedding does not hold in multi-dimensional case.
 So, we employ the  functions
\[
\frac{\rho_i-\rho_D}{\rho_i+a(\tau)}, \frac{\eta_i-\eta_D}{\eta_i+a(\tau)}\in
 H^{1}_0(\Omega\cup\Gamma_{N}):=\{u\in H^{1}(\Omega);
u=0 \text{ on } \Gamma_D\}
\]
as  test functions of the first  and third equation in
 \eqref{1.4} and, with careful calculation, get the boundedness of
 $ \{(\varepsilon\frac{\Delta\rho^{(N)}}{\sqrt{\rho^{(N)}}},
 \varepsilon\frac{\Delta\eta^{(N)}}{\sqrt{\eta^{(N)}}})\}$ and
$\{(\varepsilon\nabla(\rho^{(N)})^{2r}, \varepsilon\nabla(\eta^{(N)})^{2r})\}$.
 Using these facts and employing the Stampacchia's truncation method, we obtain
  boundedness  of $\{(\varepsilon\Delta\rho^{(N)}, \varepsilon\Delta\eta^{(N)})\}$.
Through the obtained estimates and the compactness result for piecewise constant
functions in time(see \cite{dj}) we prove for any fixed $\varepsilon\in(0,1)$
the compactness of
 $\{(\rho^{(N)}, \eta^{(N)})\}$ in $L_{p}(0, T;  H^{1}(\Omega))$ for all
$p\in  (1, \infty)$.  Furthermore, when $r\geq 9/5$, some estimates of the
 solutions independent of $\varepsilon\in(0,1)$ as well as $N$ are
 obtained.
  See the details of these stability estimates in section 3.

 \begin{theorem} \label{thm1.2}
For each fixed $\varepsilon\in(0, 1)$  there exist $ (\rho,  \eta, V )$
 and  a subsequence of approximate solutions obtained in Theorem \ref{thm1.1}
(again denoted by
$(\rho^{(N)},\eta^{(N)},V^{(N)})$) such that
\begin{equation}\label{1.8}
\begin{gathered}
\rho^{(N)}\to \rho,  \quad \eta^{(N)}\to \eta\quad
\ast\text{-weakly in } X \cap L_{4/3}(0, T;H^{3/2-\delta}(\Omega)), \quad \delta>0, \\
\rho^{(N)}\to \rho, \quad \eta^{(N)}\to \eta\quad
\text{in } L_{p}(0,T;H^{1}(\Omega)), \quad \forall p\in(1, \infty),\\
\nabla(\rho^{(N)})^{2r}\to \nabla{\rho}^{2r},  \quad
\nabla(\eta^{(N)})^{2r}\to \nabla{\eta}^{2r}\quad
\text{weakly in }L_{1}(0,T; L_{6/5}(\Omega)),\\
V^{(N)}\to V\quad \ast\text{-weakly in }L_{\infty}(0, T;H^{1}(\Omega))
\end{gathered}
\end{equation}
as $N\to\infty$ where $X:=\{u\in L_{\infty}(0, T;H^{1}(\Omega))$
$\Delta u\in L_{4/3}(0,T;L_2(\Omega))$,
$\nabla u\cdot\nu=0$ on $\Gamma_{N}\}$. Furthermore,
$(\rho^2, \eta^2, V)$ is a solution of
\eqref{1.1}, \eqref{1.2}  in the sense of
\begin{equation}\label{1.9}
\begin{gathered}
\rho, \eta\geq 0, \quad
\frac{\partial\rho^2}{\partial t}, \frac{\partial\eta^2}{\partial t} \in
L_2(0,T;W^{-1}_{r+1}(\Omega\cup\Gamma_{N})), \\
\rho-\rho_D, \eta-\eta_D,  V-V_D\in L_{\infty}(0,T; H^{1}_0(\Omega\cup\Gamma_{N})), \\
\begin{aligned}
\int^{T}_0\langle \frac{\partial\rho^2}{\partial t},  \phi\rangle\,dt
&=-2\varepsilon^2\int_{Q}\Delta \rho\nabla\rho\cdot\nabla\phi \,dx\,dt
 -\varepsilon^2\int_{Q}\Delta \rho\rho\Delta\phi \,dx\,dt \\
&\quad -\theta\int_{Q}\nabla\rho^{2r}\cdot\nabla\phi \,dx\,dt
+\int_{Q}\rho^2\nabla V\cdot\nabla\phi \,dx\,dt,
\end{aligned} \\
\begin{aligned}
\int^{T}_0\langle \frac{\partial\eta^2}{\partial t},  \phi\rangle \,dt
&=-2\varepsilon^2\int_{Q}\Delta \eta\nabla\eta\cdot\nabla\phi \,dx\,dt
 -\varepsilon^2\int_{Q}\Delta \eta\eta\Delta\phi \,dx\,dt \\
&\quad -\theta\int_{Q}\nabla\eta^{2r}\cdot\nabla\phi \,dx\,dt
-\int_{Q}\eta^2\nabla V\cdot\nabla\phi \,dx\,dt,
\end{aligned}\\
-\lambda^2\int_{Q}\nabla V\cdot \nabla\phi
\,dx\,dt=\int_{Q}(\rho^2-\eta^2-f)\phi \,dx\,dt,  \; \forall\phi\in
C_0^{\infty}(Q)
\end{gathered}
\end{equation}
where $W^{-1}_{r+1}(\Omega\cup\Gamma_{N})$ is dual space of
$\{u\in W^{1}_{(r+1)/r}(\Omega);  u=0  \text{ on } \Gamma_D\}$.
\end{theorem}

\begin{remark} \label{rmk1.3} \rm
We note that in the case of unipolar model one can
also obtain such kind of existence result by the same method.
However,  we discuss here only  the bipolar model.
\end{remark}

\begin{theorem} \label{thm1.3}
Let $(n^{(\varepsilon)}, p^{(\varepsilon)},
V^{(\varepsilon)}),  \varepsilon\in(0,1)$ be the solution of
\eqref{1.1}, \eqref{1.2} with $r \geq 9/5$ obtained in the
Theorem \ref{thm1.2}.  Then there exist some $n,  p, V$
 and sequence $\varepsilon\to0$ such that
\begin{equation}\label{1.10}
\begin{gathered}
n^{(\varepsilon)}\to n, \quad p^{(\varepsilon)}\to p\quad
\ast\text{-weakly in } L_{\infty}(0,T; L_{r}(\Omega)),\\
\nabla(n^{(\varepsilon)})^{r}\to \nabla n^{r},  \quad
\nabla(p^{(\varepsilon)})^{r}\to \nabla p^{r}\quad
\text{weakly in }L_{1}(Q),\\
V^{(\varepsilon)}\to V\;  \quad  \ast\text{-weakly in }
 L_{\infty}(0, T; H^{1}(\Omega)).
\end{gathered}
\end{equation}
Furthermore, $(n, p,  V)$ is a weak solution to the mixed boundary
value problem of classical drift-diffusion model in the sense of
\begin{equation}\label{1.12}
\begin{gathered}
n, p\geq 0, \quad \frac{\partial n}{\partial t}, \frac{\partial
p}{\partial t} \in L_2(0,T; W^{-1}_{r+1}(\Omega)),\\
\int^{T}_0\langle \frac{\partial n}{\partial t},  \phi\rangle \,dt
=-\theta\int_{Q}\nabla n^{r}\cdot\nabla\phi \,dx\,dt
+\int_{Q}n\nabla V\cdot\nabla\phi \,dx\,dt, \\
\int^{T}_0\langle \frac{\partial p}{\partial t},  \phi\rangle\,dt
=-\theta\int_{Q}\nabla p^{r}\cdot\nabla\phi \,dx\,dt
-\int_{Q}p\nabla V\cdot\nabla\phi \,dx\,dt,\\
-\lambda^2\int_{Q}\nabla V\cdot \nabla\phi
\,dx\,dt=\int_{Q}(n-p-f)\phi \,dx\,dt,  \quad \forall\phi\in
C_0^{\infty}(Q).
\end{gathered}
\end{equation}
\end{theorem}

The article is organized as follows. In section 2 we prove the
Theorem \ref{thm1.1}. Some stability estimates of the approximate solutions
are presented in section 3. Section 4 is devoted to the proof of the
convergence of the approximate solutions and semiclassical limit.


\section{Existence of approximate solutions}

First of all, we introduce Stampacchia's lemma
which will be used later.

 \begin{lemma}[{\cite[lemma 5.2.4]{j1}}]  \label{lem2.1}
 Let $\varphi:(a,b)\to R^{1}$ be a nonnegative,
 nonincreasing function where $a<b\leq+\infty$.
 Suppose that there exist constants $K>0,  r>0, \alpha>1$ such that
 \[
\varphi(\xi)\leq K(\xi-\zeta)^{-r}\varphi(\zeta)^{\alpha}, \quad  a<\zeta<\xi<b.
\]
If the number  $\xi^{\ast}=K^{\frac{1}{r}}2^{\frac{\alpha}{\alpha-1}}
\varphi(a)^{\frac{\alpha-1}{r}}$
 is such that $a+\xi^{\ast}<b,$ then $\varphi(a+\xi^{\ast})=0$.
\end{lemma}

\begin{lemma}\label{lem2.2}
The weak solution $u\in H^{1}(\Omega)$ to
\begin{equation}
\begin{gathered}
\operatorname{div}(a(x)\nabla u)=f^2_{1}-f^2_2-f_0, \;\text{in}\; \Omega,\\
u=u_D \text{ on }\Gamma_D, \quad
\frac{\partial u}{\partial\nu}=0 \text{ on } \Gamma_{N}.
\end{gathered}\label{2.1}
\end{equation}
with  $f_i\in L_{4}(\Omega)$,  $i=0, 1, 2$,
$u_D\in W^{1}_{\infty}(\Omega)$ and $a(\cdot)\in L_{\infty}(\Omega)$
satisfying $a(x)\geq a_0$ for some constant $a_0>0$ satisfies
for some constant $c(\Omega)>0$ depending only on $\Omega$
\begin{equation}\label{2.2}
\begin{gathered}
 u(x)\leq u^{\ast}:=\|u_D\|_{L_{\infty}(\Gamma_D)}
+c(\Omega)a^{-1}_0(\|f_2\|^2_{L_{4}(\Omega)}
+\|f_0\|_{L_2(\Omega)}),\\
 u(x)\geq -u_{\ast}:=-\Big(\|u_D\|_{L_{\infty}(\Gamma_D)}
 +c(\Omega)a^{-1}_0(\|f_{1}\|^2_{L_{4}(\Omega)}
+\|f_0\|_{L_2(\Omega)})\Big),
\end{gathered}
\end{equation}
a.e.\ in $\Omega$.
\end{lemma}

The above lemma can be proved as in \cite[Lemma 5.2.5]{j1},
using Lemma \ref{lem2.1}.

For $\rho_{i-1}, \eta_{i-1}\in H^{1}(\Omega)\cap
L_{\infty}(\Omega)$ satisfying
$$
\rho_{i-1}|_{\Gamma_D}=\rho_D, \quad
\eta_{i-1}|_{\Gamma_D}=\eta_D, \quad
\exists m_{i-1}>0; \rho_{i-1}, \eta_{i-1}\geq m_{i-1},
$$
we consider the auxiliary boundary-value problem
\begin{equation}\label{2.3}
 \begin{gathered}
\operatorname{div}((S_{M}(\rho)+a(\tau))^2\nabla F)
 =\frac{\delta}{\tau}((\rho+a(\tau))^2-(\xi +a(\tau))^2),\\
\varepsilon^2\Delta \rho=\rho_{+}(\theta\delta h(S^2_{M}(\rho))
+\tau\delta \ln \rho^2-F-\delta V),\\
\operatorname{div}((S_{M}(\eta)+a(\tau))^2\nabla G)
 =\frac{\delta}{\tau}((\eta+a(\tau))^2-(\zeta +a(\tau))^2),\\
\varepsilon^2\Delta \eta=\eta_{+}(\theta\delta h(S^2_{M}(\eta))
+\tau\delta \ln \eta^2-G+\delta V),\\
\lambda^2\Delta V=(\rho+a(\tau))^2-(\eta+a(\tau))^2-f\quad
\text{in } \Omega,
\end{gathered}
\end{equation}
with
\begin{equation}\label{2.4}
\begin{gathered}
(\rho, \eta, F, G, V)=(\rho_{D\delta}, \eta_{D\delta}, F_{D\delta},
G_{D\delta}, V_D)\quad \text{on } \Gamma_D,\\
\nabla\rho\cdot\nu=\nabla\eta\cdot\nu=\nabla F\cdot\nu=\nabla
G\cdot\nu=\nabla V\cdot\nu=0\quad \text{on } \Gamma_{N}
\end{gathered}
\end{equation}
where $\delta\in (0,1]$, $M\geq1$ are constants,
$S_{M}(\varphi)=\min \{M, \max \{\varphi,0\}\}$,
$\varphi_{+}=\max \{\varphi,0\}$, and
\begin{gather*}
\xi=\delta \rho_{i-1}, \quad \zeta=\delta\eta_{i-1}, \quad
\rho_{D\delta}=\delta\rho_D, \quad \eta_{D\delta}=\delta\eta_D,\\
F_{D\delta}= \delta(\theta h(\delta^2\rho^2_D)+\tau\ln (\delta\rho_D)^2-V_D),\\
G_{D\delta}= \delta (\theta h(\delta^2\eta^2_D)+\tau\ln (\delta\eta_D)^2+ V_D).
\end{gather*}

\begin{lemma} \label{lem2.3}
The weak solution
$(\rho, \eta, F, G, V)\in (H^{1}(\Omega))^{5}$ to
\eqref{2.3}, \eqref{2.4} satisfies
\begin{gather}\label{2.5}
\rho(x), \eta(x)\geq m, \quad \text{a. e. in } \Omega, \\
\label{2.6}
\|(\rho, \eta)\|_{(L_{\infty}(\Omega))^2}\leq c(
\|(\rho, \eta)\|_{(L_{4}(\Omega))^2})
\end{gather}
for a constant $m>0$, and
$c( \|(\rho, \eta)\|_{(L_{4}(\Omega))^2})$ which  is also
bounded if $\|(\rho, \eta)\|_{(L_{4}(\Omega))^2}$ is bounded.
\end{lemma}

We note that the $L_{\infty}$-bound depends on  $M$ in the
definition of the function $S_{M}(\cdot)$.

\begin{proof}
Taking $\rho_{-}=\min \{\rho, 0\}$,
$\eta_{-}=\min \{\eta, 0\}\in H_0^{1}(\Omega\cup \Gamma_{N})$
as test functions of the second and fourth equation of \eqref{2.3},
respectively, we can easily
verify the nonnegativity of $\rho, \eta$. Also, by Lemma \ref{lem2.2} we
obtain the lower and upper bounds for $F, G, V$
\begin{equation} \label{2.7}
\begin{gathered}
V(x)\leq V^{\ast}:=\|V_D\|_{L_{\infty}(\Gamma_D)}
+c(\Omega, \lambda)\Big(\|\eta+a(\tau)\|_{L_{4}(\Omega)}^2
+\|f\|_{L_{\infty}(\Omega)}\Big),\\
V(x)\geq -V_{\ast}:=-\Big(\|V_D\|_{L_{\infty}(\Gamma_D)}
+c(\Omega, \lambda)(\|\rho+a(\tau)\|_{L_{4}(\Omega)}^2
 +\|f\|_{L_{\infty}(\Omega)})\Big),\\
F(x)\leq F^{\ast}:=\|F_{D\delta}\|_{L_{\infty}(\Gamma_D)}
+c(\Omega, \tau, \delta)\|\xi+a(\tau)\|_{L_{4}(\Omega)}^2,\\
F(x)\geq -F_{\ast}:=-\Big(\|F_{D\delta}\|_{L_{\infty}(\Gamma_D)}
+c(\Omega, \tau, \delta)\|\rho+a(\tau)\|_{L_{4}(\Omega)}^2\Big),\\
G(x)\leq G^{\ast}:=\|G_{D\delta}\|_{L_{\infty}(\Gamma_D)}
 +c(\Omega, \tau, \delta)\|\zeta+a(\tau)\|_{L_{4}(\Omega)}^2,\\
 G(x)\geq -G_{\ast}:=-\Big(\|G_{D\delta}\|_{L_{\infty}(\Gamma_D)}
+c(\Omega, \tau, \delta)\|\eta+a(\tau)\|_{L_{4}(\Omega)}^2\Big)
\end{gathered}
\end{equation}
for some constants
$c(\Omega, \lambda), c(\Omega, \tau, \delta)>0$. Hence, we have
\begin{equation}\label{2.8}
\rho(x)\leq
K_{\rho}:=\max \{1, \exp (\frac{1}{2\tau\delta}(F^{\ast}+V^{\ast}))\}
=c(\|\xi\|_{L_{4}(\Omega)}, \|\eta\|_{L_{4}(\Omega)}),
\end{equation}
a. e.\ in $\Omega$.
The second equation in \eqref{2.3} yields
\begin{align*}
\varepsilon^2\int_{\Omega}|\nabla (\rho-K_{\rho})_{+}|^2dx
&=-\int_{\Omega}\rho_{+}(\theta\delta
h(S^2_{M}(\rho))+\tau\delta\ln \rho^2-F
-\delta V)(\rho-K_{\rho})_{+}dx\\
&\leq \int_{\Omega}\rho_{+}(F^{\ast}+V^{\ast}-\tau\delta\ln K_{\rho}^2)
(\rho-K_{\rho})_{+}dx\leq0.\\
\end{align*}
In the same way we obtain the upper bound for $\eta$ as
\begin{equation}\label{2.9}
\eta(x)\leq
K_{\eta}:=\max \{1, \exp (\frac{1}{2\tau\delta}(G^{\ast}
+V_{\ast}))\}=c(\|\zeta\|_{L_{4}(\Omega)}, \|\rho\|_{L_{4}(\Omega)}),
\end{equation}
a. e.\ in $\Omega$.

From \eqref{2.8}, \eqref{2.9} we obtain \eqref{2.6}.
To obtain the lower bounds
\begin{equation} \label{2.10}
\begin{gathered}
\rho(x)\geq
m_{\rho}:=\min \{\delta\rho_D,
\exp (-\frac{1}{2\tau\delta}(F_{\ast}+V_{\ast}))\},\quad
\text{a.e.\ in }\Omega,\\
\eta(x)\geq
m_{\eta}:=\min \{\delta\eta_D, \exp (-\frac{1}{2\tau\delta}(G_{\ast}+V^{\ast}))\},
\quad\text{a.e.\ in }\Omega,
\end{gathered}
\end{equation}
we take $(\rho-m_{\rho})_{-}$, $(\eta-m_{\eta})_{-}\in
H_0^{1}(\Omega\cup \Gamma_{N})$ as test functions of the second
and fourth equation of \eqref{2.3} respectively and use \eqref{2.7}.
\end{proof}

\begin{lemma} \label{lem2.4}
Assume that $\tau=N/T$ satisfies \eqref{1.6}. Then the weak solution\\
$(\rho, \eta, F, G, V)\in (H^{1}(\Omega))^{5}$  to
\eqref{2.3}, \eqref{2.4} satisfies
\begin{equation}\label{2.12}
\|(\rho, \eta)\|_{(H^{1}(\Omega))^2}\leq c
\end{equation}
for some constant $c>0$ independent of the solution,
$\delta\in (0, 1]$ and the choice of $M\geq1$.
\end{lemma}

\begin{proof}
By Lemma \ref{lem2.3},
\[
F=-\varepsilon^2\frac{\Delta\rho}{\rho}
 +\theta\delta h(S^2_{M}(\rho))+\tau\delta\ln \rho^2-\delta V, \quad
G=-\varepsilon^2\frac{\Delta\eta}{\eta}+\theta\delta
h(S^2_{M}(\eta))+\tau\delta\ln \eta^2+\delta V.
\]
Now, we take $(F-F_{D\delta})\in H^{1}_0(\Omega\cap\Gamma_{N})$ as test
function of the first equation in \eqref{2.3} to obtain
\begin{equation} \label{2.13}
\begin{aligned}
&\frac{1}{2}\int_{\Omega}(S_{M}(\rho)+a(\tau))^2|\nabla
F_{D\delta}|^2dx \\
&\geq\frac{1}{2}\int_{\Omega}(S_{M}(\rho)+a(\tau))^2|\nabla F|^2dx\\
&\quad +\frac{\delta}{\tau}\int_{\Omega}(\rho^2-\xi^2+2a(\tau)(\rho-\xi))
\Big(-\varepsilon^2\frac{\Delta \rho}{\rho}\Big)\,dx\\
&\quad +\delta^2\int_{\Omega}(\rho^2-\xi^2+2a(\tau)(\rho-\xi))\ln \rho^2dx\\
&\quad +\frac{\delta^2}{\tau}\int_{\Omega}((\rho+a(\tau))^2
 -(\xi+a(\tau))^2)(-V+V_D)\,dx\\
&\quad +\frac{\theta\delta^2}{\tau}\int_{\Omega}(\rho^2-\xi^2
 +2a(\tau)(\rho-\xi))h(S^2_{M}(\rho))\,dx\\
&\quad +\frac{\delta}{\tau}\int_{\Omega}((\rho+a(\tau))^2
 -(\xi+a(\tau))^2)(-F_{D\delta}-\delta
V_D)\,dx \\
&=\sum_{j=1}^{6}R_{j}.
\end{aligned}
\end{equation}
We estimate term by term. Integration by parts and Young's inequality
yield
\begin{equation}\label{2.14}
\begin{aligned}
R_2&=\frac{\delta\varepsilon^2}{\tau}\int_{\Omega}
 \nabla\Big(\frac{\rho^2-\xi^2}{\rho}\Big)
\cdot\nabla\rho dx
+2\frac{\delta\varepsilon^2a(\tau)}{\tau}\int_{\Omega}
 \nabla\Big(\frac{\rho-\xi}{\rho}\Big) \cdot\nabla\rho dx\\
&=\frac{\delta\varepsilon^2}{\tau}\Big(\int_{\Omega}|\nabla\rho|^2dx
-\int_{\Omega}|\nabla\xi|^2dx
+\int_{\Omega}|\nabla\xi-\frac{\xi}{\rho}\nabla\rho|^2dx\Big)\\
&\quad +2\frac{\delta\varepsilon^2a(\tau)}{\tau}\int_{\Omega}
\Big(\frac{\xi}{\rho^2}|\nabla\rho|^2
-\frac{\nabla\rho\cdot\nabla\xi}{\rho}\Big)\,dx\\
&\geq\frac{\delta\varepsilon^2}{\tau}\int_{\Omega}|\nabla\rho|^2dx
-\frac{\delta\varepsilon^2}{\tau}\int_{\Omega}|\nabla\xi|^2dx
-\frac{\delta\varepsilon^2a(\tau)}{\tau}\int_{\Omega}\frac{|\nabla\xi|^2}{\xi}dx.
\end{aligned}
\end{equation}
The estimate for  $R_{3}$ is
\begin{equation}\label{2.15}
\begin{aligned}
R_{3}&=\delta^2\int_{\Omega}(\rho^2-\xi^2)\ln \rho^2dx
+4\delta^2a(\tau)\int_{\Omega}(\rho-\xi)\ln \rho dx\\
&\geq\delta^2\int_{\Omega}(H(\rho^2)-H(\xi^2))\,dx
+4\delta^2a(\tau)\int_{\Omega}(H(\rho) -H(\xi))\,dx,
\end{aligned}
\end{equation}
where $H(\alpha):=\alpha(\ln \alpha-1)+1$ and $\alpha>0$,
which is well-known (see \cite[Lemma 2.2]{jp2}).
 We rewrite $R_5$ as
\begin{align*}
R_{5}&=\frac{r\theta\delta^2}{(r-1)\tau}\int_{\Omega(\rho>1, \rho>\xi)}
 ((\rho^2-\xi^2)+2a(\tau)(\rho-\xi))(S^{2(r-1)}_{M}(\rho)-1)\,dx\\
&\quad +\frac{r\theta\delta^2}{(r-1)\tau}\int_{\Omega(\rho>1, \rho\leq\xi)}
 ((\rho^2-\xi^2)+2a(\tau)(\rho-\xi))(S^{2(r-1)}_{M}(\rho)-1)\,dx\\
&\quad +\frac{r\theta\delta^2}{(r-1)\tau}\int_{\Omega(\rho\leq1)}
 ((\rho^2-\xi^2)+2a(\tau)(\rho-\xi))(S^{2(r-1)}_{M}(\rho)-1)\,dx\\
&=R_{5,1}+R_{5,2}+R_{5,3}.
\end{align*}
Since $R_{5,1}\geq0$,
\begin{align*}
R_{5,2}&\geq-\frac{r\theta\delta^2}{(r-1)\tau}\int_{\Omega(\rho>1,
\rho\leq\xi)}(\xi^2+2a(\tau)\xi)\Big(S^{2(r-1)}_{M}(\rho)-1\Big)\,dx\\
&\geq-\frac{r\theta\delta^2}{(r-1)\tau}\int_{\Omega(\rho>1, \rho\leq\xi)}
(\xi^{2r}+2a(\tau)\xi^{2r-1})\,dx,
\end{align*}
and
$$
 R_{5,3}\geq-\frac{r\theta\delta^2}{(r-1)\tau}\int_{\Omega(\rho\leq1)}
(\rho^2+2a(\tau)\rho)\,dx\geq-\frac{r\theta\delta^2(1+2a(\tau))}{(r-1)
\tau}\operatorname{meas}(\Omega),
$$
we obtain
\begin{equation}\label{2.16}
R_{5}\geq-\frac{r\theta\delta^2}{(r-1)\tau}
\Big(\int_{\Omega}(\xi^{2r}+2a(\tau)\xi^{2r-1})\,dx+(1+2a(\tau))
\operatorname{meas}(\Omega)\Big).
\end{equation}
Using the above inequalities we have
\begin{align*}
&\varepsilon^2\int_{\Omega}|\nabla\rho|^2dx
+\tau\delta\int_{\Omega}H(\rho^2)\,dx
 +4\tau\delta a(\tau)\int_{\Omega}H(\rho)\,dx\\
&+\frac{\tau}{2\delta}\int_{\Omega}(S_{M}(\rho)+a(\tau))^2|\nabla
F|^2dx
-\delta\int_{\Omega}((\rho+a(\tau))^2-(\xi+a(\tau))^2)(V-V_D)\,dx\\
&\leq\varepsilon^2\int_{\Omega}|\nabla\xi|^2dx
+\tau\delta\int_{\Omega}H(\xi^2)\,dx+4\tau\delta a(\tau)\int_{\Omega}H(\xi)\,dx\\
&\quad +\varepsilon^2a(\tau)\int_{\Omega}\frac{|\nabla\xi|^2}{\xi}dx
+\frac{\tau}{2\delta}\int_{\Omega}(\rho+a(\tau))^2|\nabla F_{D\delta}|^2dx\\
&\quad +\frac{r\theta\delta}{(r-1)}(\int_{\Omega}(\xi^{2r}
 +2a(\tau)\xi^{2r-1})\,dx+(1+2a(\tau))\operatorname{meas}(\Omega))\\
&\quad +\int_{\Omega}((\rho+a(\tau))^2-(\xi+a(\tau))^2)(F_{D\delta}+\delta
V_D)\,dx.
\end{align*}
We obtain a similar inequality for $\eta$ by taking
$(G-G_{D\delta})\in H^{1}_0(\Omega\cap\Gamma_{N})$ as test
function of the third equation in \eqref{2.3}. Adding the resulting
two inequalities and taking into account that
\begin{gather*}
\begin{aligned}
&\frac{\tau}{2\delta}\Big(\int_{\Omega}(\rho+a(\tau))^2
|\nabla F_{D\delta}|^2dx+\int_{\Omega}(\eta+a(\tau))^2|\nabla G_{D\delta}|^2dx\Big)\\
&\leq\frac{\delta\tau}{2}\|\nabla V_D\|_{L_{\infty}(\Omega)}^2
\int_{\Omega}((\rho+a(\tau))^2+(\eta+a(\tau))^2)\,dx,
\end{aligned}\\
\begin{aligned}
&\int_{\Omega}(\rho+a(\tau))^2(\delta V_D+F_{D\delta})\,dx
 +\int_{\Omega}(\eta+a(\tau))^2(G_{D\delta}-\delta V_D)\,dx\\
&\leq\theta\delta h(\rho^2_D)\int_{\Omega}(\rho+a(\tau))^2dx
 +\theta\delta h(\eta^2_D)\int_{\Omega}(\eta+a(\tau))^2dx,
\end{aligned}\\
\begin{aligned}
&-\delta\int_{\Omega}(((\rho+a(\tau))^2
-(\eta+a(\tau))^2)-((\xi+a(\tau))^2-(\zeta+a(\tau))^2))(V-V_D)\,dx\\
&=\delta\lambda^2\int_{\Omega}\nabla(V-V_D)\cdot\nabla((V-V_D)-(V'-V_D))\,dx\\
&\geq\frac{\delta\lambda^2}{2}\Big(\int_{\Omega}|\nabla(V-V_D)|^2dx
-\int_{\Omega}|\nabla(V'-V_D)|^2dx\Big)
\end{aligned}
\end{gather*}
we have
\begin{align*}
&\varepsilon^2\int_{\Omega}(|\nabla\rho|^2+|\nabla\eta|^2)\,dx
-\theta\delta\Big(h(\rho^2_D)\int_{\Omega}(\rho+a(\tau))^2dx
+h(\eta^2_D)\int_{\Omega}(\eta+a(\tau))^2dx\Big)\\
&\leq\varepsilon^2\int_{\Omega}(|\nabla\xi|^2+|\nabla\zeta|^2)\,dx
+\tau\int_{\Omega}(H(\xi^2)+H(\zeta^2))\,dx\\
&\quad +4\tau a(\tau)\int_{\Omega}(H(\xi)+H(\zeta))\,dx
+\frac{\lambda^2}{2}\int_{\Omega}|\nabla(V'-V_D)|^2dx\\
&\quad +\frac{\delta\tau}{2}\|\nabla
V_D\|_{L_{\infty}(\Omega)}^2\int_{\Omega}((\rho+a(\tau))^2+(\eta+a(\tau))^2)\,dx\\
&\quad +\varepsilon^2a(\tau)\int_{\Omega}\Big(\frac{|\nabla\xi|^2}{\xi}
+\frac{|\nabla\zeta|^2}{\zeta}\Big)\,dx
+c\int_{\Omega}((\xi+a(\tau))^2
+(\zeta+a(\tau))^2)\,dx\\
&\quad +\frac{r\theta}{(r-1)}\Big(\int_{\Omega}(\xi^{2r}+\zeta^{2r}
+2a(\tau)(\xi^{2r-1}+\zeta^{2r-1}))\,dx+2(1+2a(\tau))\operatorname{meas}(\Omega)\Big)
\end{align*}
for some $c>0$ independent of the solution, $\delta\in (0, 1]$ and
the choice of $M\geq1$. Here $V'$ is the solution to
\begin{gather*}
\lambda^2\triangle V=(\xi+a(\tau))^2-(\zeta+a(\tau))^2-f \quad \text{in } \Omega,\\
V=V_D  \text{ on }\Gamma_D, \quad
\frac{\partial V}{\partial\nu}=0 \text{ on } \Gamma_{N},
\end{gather*}
This completes the proof with the help of
the choice of $\tau$ and $\xi=\delta \rho_{i-1},
\zeta=\delta\eta_{i-1}$.
\end{proof}


\begin{proof}[Proof of Theorem \ref{thm1.1}]
 We define a fixed point mapping
$R:\;(L_{4}(\Omega))^2\times[0,1]\to(L_{4}(\Omega))^2$
as follows. Let $\delta\in [0,1]$, $u, w\in(L_{4}(\Omega))^2$
be given. Then, by the theory of elliptic equation, there exists a
unique weak solution $(\rho, \eta, F, G, V)\in
(H^{1}(\Omega))^{5}$ to
\begin{equation}\label{2.16b}
 \begin{gathered}
 \operatorname{div}((S_{M}(u)+a(\tau))^2\nabla F)
 =\frac{\delta}{\tau}((u+a(\tau))^2-(\xi+a(\tau))^2),\\
\varepsilon^2\Delta \rho=u_{+}(\theta\delta h(S^2_{M}(u))
+\tau\delta\ln u^2-F-\delta V),\\
\operatorname{div}((S_{M}(w)+a(\tau))^2\nabla G)
 =\frac{\delta}{\tau}((w+a(\tau))^2-(\zeta+a(\tau))^2),\\
\varepsilon^2\Delta \eta=w_{+}(\theta\delta h(S^2_{M}(w))
+\tau\delta \ln w^2-G+\delta V),\\
\lambda^2\Delta V=(u+a(\tau))^2-(w+a(\tau))^2-f\quad
\text{in } \Omega
\end{gathered}
\end{equation}
with the boundary condition \eqref{2.4} and
$\xi=\delta\rho_{i-1}$, $\zeta=\delta\eta_{i-1}$. Hence, the mapping
given by $R((u, w), \delta)=(\rho, \eta)$ is well defined.
Moreover, $R((u, w), 0)=0$ for any $u, w\in L_{4}(\Omega)$.
We can easily verify that $R$ is continuous and compact by standard
argument. Lemma \ref{lem2.4} shows that there is a constant $c>0$ such
that for all $\rho, \eta\in L_{4}(\Omega), \delta\in[0, 1]$
satisfying $R((\rho, \eta), \delta)=(\rho, \eta)$ it holds
$\|(\rho, \eta)\|_{(L_{4}(\Omega))^2}\leq c$. Hence, by the
Leray-Schauder fixed point theorem there exists  $(\rho, \eta)$
satisfying $R((\rho, \eta), 1)=(\rho, \eta)$ which is the
solution to \eqref{2.3}, \eqref{2.4} with $\delta=1$ depending on
the choice of $M$. However, taking into account of the estimate in
Lemmas \ref{lem2.3} and \ref{lem2.4}, and the embedding
$H^{1}(\Omega)\hookrightarrow L_{4}(\Omega)$, the solution satisfies
$\rho(x), \eta(x)\leq M$ for some $M$ large enough. This completes
the proof.
\end{proof}

\section{Stability estimates}

Let $N=N_0,N_0+1,\dots,$
$(\rho_i, \eta_i, F_i, G_i, V_i)\in(L_{\infty}(\Omega)\cap
H^{1}(\Omega))^{5}$,  $i=1,2,\dots N$ be the recursively
defined solutions to \eqref{1.4} and $(\rho^{(N)},
\eta^{(N)}, F^{(N)},  G^{(N)}, V^{(N)})$ be the approximate
solutions defined by \eqref{1.5}.

\begin{lemma}\label{lem3.1}
There exist  constant $c>0$ and integer
$N^{\ast}\geq N_0$ such that for all
$N=N^{\ast}, N^{\ast}+1,\dots$ and
$\varepsilon\in(0,1)$ it holds
\begin{equation}\label{3.1}
\begin{aligned}
& \|(\rho^{(N)},\eta^{(N)})\|_{(L_{\infty}(0,T;L_2(\Omega)))^2}^2+
\|(\varepsilon\frac{\Delta\rho^{(N)}}{\sqrt{\rho^{(N)}}},
 \varepsilon\frac{\Delta\eta^{(N)}}{\sqrt{\eta^{(N)}}})\|_{(L_2(Q))^2}^2\\
&+\tau\|(\frac{\nabla\rho^{(N)}}{\sqrt{\rho^{(N)}}},
 \frac{\nabla\eta^{(N)}}{\sqrt{\eta^{(N)}}})\|_{(L_2(Q))^2}^2
+\sum^{N}_{i=1}\int_{\Omega}(\frac{(\rho_i-\rho_{i-1})^2}
{\rho_i+a(\tau)}+\frac{(\eta_i-\eta_{i-1})^2}{\eta_i+a(\tau)})\,dx\\
&+\|(\nabla(\rho^{(N)})^{r-1/2},\nabla(\eta^{(N)})^{r-1/2})\|_{(L_2(Q))^2}^2
\leq c.
\end{aligned}
\end{equation}
\end{lemma}


\begin{proof}
We take $\phi=\frac{\rho_i-\rho_D}{\rho_i+a(\tau)}\in
H^{1}_0(\Omega\cup\Gamma_{N})$ as test function in the first
equation of \eqref{1.4} to obtain
\begin{equation}\label{3.2}
\begin{aligned}
& -\int_{\Omega}(\rho_i+a(\tau))^2\nabla
F_i \cdot\nabla\frac{\rho_i-\rho_D}{\rho_i+a(\tau)} dx\\
&=\frac{1}{\tau}\int_{\Omega}((\rho_i+a(\tau))^2-(\rho_{i-1}+a(\tau))^2)
\frac{\rho_i-\rho_D}{\rho_i+a(\tau)}dx.
\end{aligned}
\end{equation}
Using the  equality
$\alpha(\alpha-\beta)=\frac{1}{2}(\alpha^2-\beta^2+(\alpha-\beta)^2),$
we rewrite the  right hand side of \eqref{3.2} as
\begin{align*}
R&=\frac{1}{\tau}\int_{\Omega}\Big(2(\rho_i+a(\tau))(\rho_i-\rho_{i-1})-
(\rho_i-\rho_{i-1})^2\Big)\frac{\rho_i-\rho_D}{\rho_i+a(\tau)}dx\\
&=\frac{1}{\tau}\Big(\int_{\Omega}(\rho_i-\rho_D)^2dx+
\int_{\Omega}\frac{\rho_D+a(\tau)}{\rho_i+a(\tau)}(\rho_i-\rho_{i-1})^2dx
-\int_{\Omega}(\rho_{i-1}-\rho_D)^2dx\Big).
\end{align*}
The left-hand side of \eqref{3.2} can be written as
\begin{align*}
L&=-\int_{\Omega}(\rho_i+a(\tau))^2\nabla\Big(-\varepsilon^2\frac{\Delta\rho_i}{\rho_i}
+\tau\ln \rho^2_i+\theta
h(\rho^2_i)-V_i\Big)
\cdot\nabla\frac{\rho_i-\rho_D}{\rho_i+a(\tau)}dx\\
& =-(\rho_D+a(\tau))\int_{\Omega}\frac{\varepsilon^2
(\Delta\rho_i)^2+2\tau|\nabla\rho_i|^2}{\rho_i}dx\\
&\quad -(\rho_D+a(\tau))\Big(2r\theta\int_{\Omega}\rho^{2r-3}_i|\nabla\rho_i|^2dx
-\int_{\Omega}\nabla
V_i\cdot\nabla(\rho_i-\rho_D)\,dx\Big).
\end{align*}
Hence, from \eqref{3.2} we have
\begin{equation}\label{3.3}
\begin{aligned}
& \frac{1}{\rho_D+a(\tau)}\int_{\Omega}(\rho_i-\rho_D)^2dx
+\int_{\Omega}\frac{(\rho_i-\rho_{i-1})^2}{\rho_i+a(\tau)}dx\\
&+\tau\int_{\Omega}\frac{\varepsilon^2(\Delta\rho_i)^2+2\tau|\nabla\rho_i|^2}{\rho_i}dx
+2r\theta\tau\int_{\Omega}\rho^{2r-3}_i|\nabla\rho_i|^2dx\\
&=\frac{1}{\rho_D+a(\tau)}\int_{\Omega}(\rho_{i-1}-\rho_D)^2dx
+\tau\int_{\Omega} \nabla V_i\cdot\nabla(\rho_i-\rho_D)\,dx.
\end{aligned}
\end{equation}
Similar inequality can be obtained by taking
$\eta=\frac{\eta_i-\eta_D}{\eta_i+a(\tau)}$ as test function
in the third equation of \eqref{1.4}. We observe that
\begin{align*}
& \tau\int_{\Omega}\nabla
V_i\cdot\nabla((\rho_i-\rho_D)-(\eta_i-\eta_D))\,dx\\
&=\tau\lambda^{-2}\int_{\Omega}(f-(\rho_i+a(\tau))^2+
(\eta_i+a(\tau))^2)\Big((\rho_i+a(\tau)) -(\eta_i+a(\tau))\\
&\quad -(\rho_D-\eta_D)\Big)\,dx\\
&\leq c\tau\int_{\Omega} ((\rho_i-\rho_D)^2
+(\eta_i-\eta_D)^2)\,dx+c\tau
\end{align*}
for some constant $c>0$.  Here and in the following we denote the
constants independent of $N=N_0,N_0+1,\dots$, $i=1,2,\dots,N$,
$\varepsilon\in(0,1)$ and solution as $c$. Add the resulting
inequalities for $\rho, \eta$  with respect $i=1, 2,\dots, k\leq
N$  to have
\begin{align*}
& \Big(\frac{1}{\rho_D+a(\tau)}-c\tau\Big)
\int_{\Omega}(\rho_{k}-\rho_D)^2dx
+\Big(\frac{1}{\eta_D+a(\tau)}-c\tau\Big)
\int_{\Omega}(\eta_{k}-\eta_D)^2dx\\
&+\sum^{k}_{i=1}\int_{\Omega}\Big(\frac{(\rho_i-\rho_{i-1})^2}
{\rho_i+a(\tau)}+\frac{(\eta_i-\eta_{i-1})^2}{\eta_i+a(\tau)}\Big)\,dx\\
&+\tau\sum^{k}_{i=1}\int_{\Omega}\Big(\frac{\varepsilon^2(\Delta\rho_i)^2
+2\tau|\nabla\rho_i|^2}{\rho_i}+\int_{\Omega}\frac{\varepsilon^2(\Delta\eta_i)^2
+2\tau|\nabla\eta_i|^2}{\eta_i}\Big)\,dx\\
&+\frac{2r\theta\tau}{(r+\frac{1}{2})^2}\sum^{k}_{i=1}
\Big(\int_{\Omega}|\nabla\rho_i^{r-1/2}|^2dx
+\int_{\Omega}|\nabla\eta_i^{r-1/2}|^2dx\Big)\\
&\leq\frac{1}{\rho_D+a(\tau)}\int_{\Omega}(\rho_0-\rho_D)^2dx
+\frac{1}{\eta_D+a(\tau)}\int_{\Omega}(\eta_0-\eta_D)^2dx\\
&\quad +c\tau\sum^{k-1}_{i=1}\int_{\Omega}((\rho_i-\rho_D)^2
 +(\eta_i-\eta_D)^2)\,dx+cT.
\end{align*}
This proves  \eqref{3.1} with the help of the discrete Gronwall
lemma \cite[Theorem 2.18]{br}.
\end{proof}

\begin{lemma} \label{lem3.2}
 There exist a constants $c>0$ independent
$N=N^{*}, N^{*}+1,\dots$ and $\varepsilon\in(0,1)$ such that
\begin{equation}\label{3.5}
\begin{aligned}
& \|(\varepsilon\rho^{(N)},\varepsilon\eta^{(N)},V^{(N)})
 \|_{(L_{\infty}(0,T;H^{1}(\Omega)))^{3}}^2
+\|(\rho^{(N)},\eta^{(N)})\|_{(L_{\infty}(0,T;L_{2r}(\Omega)))^2}^{2r}\\
&+\|((\rho^{(N)}+a(\tau))\nabla F^{(N)},(\eta^{(N)}+a(\tau))\nabla
G^{(N)})\|_{(L_2(Q))^2}^2\leq c
\end{aligned}
\end{equation}
\end{lemma}

\begin{proof}  Using
\begin{gather*}
F_i-F_D=-\varepsilon^2\frac{\Delta\rho_i}{\rho_i}+\theta
h(\rho_i^2)+\tau\ln \rho_i^2- V-F_D,\\
G_i-G_D=-\varepsilon^2\frac{\Delta\eta_i}{\eta_i}+\theta
h(\eta_i^2)+\tau\ln \eta_i^2+ V-G_D.
\end{gather*}
as test functions in the first and third equations of \eqref{1.4}
similarly to the proof of Lemma \ref{lem2.4}, we  obtain that for any
$i=1, 2,\dots, N$,
\begin{align*}
& \varepsilon^2\int_{\Omega}|\nabla\rho_i|^2dx
+\tau\int_{\Omega}H(\rho_i^2)\,dx+4\tau a(\tau)\int_{\Omega}H(\rho_i)\,dx\\
& +\frac{\tau}{2}\int_{\Omega}(\rho_i+a(\tau))^2|\nabla
F_i|^2dx-\int_{\Omega}((\rho_i+a(\tau))^2-(\rho_{i-1}+a(\tau))^2)(V_i-V_D)\,dx\\
& +\frac{\theta}{r-1}\int_{\Omega}\Big((\rho_i^{2r}-\rho_{i-1}^{2r})
+\frac{2a(\tau)r}{2r-1}(\rho_i^{2r-1}-\rho_{i-1}^{2r-1})\Big)\,dx\\
& -\frac{\theta r}{r-1}\int_{\Omega}((\rho_i^2-\rho_{i-1}^2)+2a(\tau)(\rho_i-\rho_{i-1}))\,dx\\
&\leq\varepsilon^2\int_{\Omega}|\nabla\rho_{i-1}|^2dx
+\tau\int_{\Omega}H(\rho_{i-1}^2)\,dx+4\tau a(\tau)\int_{\Omega}H(\rho_{i-1})\,dx\\
&\quad +\varepsilon^2a(\tau)\int_{\Omega}\frac{|\nabla\rho_{i-1}|^2}{\rho_{i-1}}dx
+\frac{\tau}{2}\int_{\Omega}(\rho_i+a(\tau))^2|\nabla F_D|^2dx\\
&\quad +\int_{\Omega}((\rho_i+a(\tau))^2-(\rho_{i-1}+a(\tau))^2)(F_D+V_D)\,dx
\end{align*}
and
\begin{align*}
& \varepsilon^2\int_{\Omega}|\nabla\eta_i|^2dx
+\tau\int_{\Omega}H(\eta_i^2)\,dx+4\tau a(\tau)\int_{\Omega}H(\eta_i)\,dx\\
& +\frac{\tau}{2}\int_{\Omega}(\eta_i+a(\tau))^2|\nabla
G_i|^2dx+\int_{\Omega}((\eta_i+a(\tau))^2-(\eta_{i-1}+a(\tau))^2)(V_i-V_D)\,dx\\
& +\frac{\theta}{r-1}\int_{\Omega}((\eta_i^{2r}-\eta_{i-1}^{2r})+\frac{2a(\tau)r}{2r-1}(\eta_i^{2r-1}-\eta_{i-1}^{2r-1}))\,dx\\
& -\frac{\theta r}{r-1}\int_{\Omega}((\eta_i^2-\eta_{i-1}^2)+2a(\tau)(\eta_i-\eta_{i-1}))\,dx\\
&\leq\varepsilon^2\int_{\Omega}|\nabla\eta_{i-1}|^2dx
+\tau\int_{\Omega}H(\eta_{i-1}^2)\,dx+4\tau a(\tau)\int_{\Omega}H(\eta_{i-1})\,dx\\
&\quad +\varepsilon^2a(\tau)\int_{\Omega}\frac{|\nabla\eta_{i-1}|^2}{\eta_{i-1}}dx
+\frac{\tau}{2}\int_{\Omega}(\eta_i+a(\tau))^2|\nabla G_D|^2dx\\
&\quad +\int_{\Omega}((\eta_i+a(\tau))^2-(\eta_{i-1}+a(\tau))^2)(G_D-
V_D)\,dx.
\end{align*}
Here we use Young's inequality
 $\alpha\beta\leq\frac{1}{p}\alpha^{p}+\frac{1}{q}\beta^{q}, \frac{1}{p}+\frac{1}{q}=1$
 to estimate the terms of the type
 $((\alpha+a(\tau))^2-(\beta+a(\tau))^2)h(\alpha^2)$.
 Adding the above two inequalities similarly to the proof of Lemma \ref{lem2.4} and
 again adding the resulting inequalities for $i=1, 2,\dots,k\leq N$, we obtain the conclusion 
with the help of Lemma \ref{lem3.1}.
\end{proof}

\begin{lemma} \label{lem3.3}
There exist constants $c>0$ independent
$N=N^{*}, N^{*}+1,\dots$ and $\varepsilon\in(0,1)$ such that
\begin{equation} \label{3.6}
\begin{gathered}
\|(\varepsilon\nabla(\rho^{(N)})^{2r}, \varepsilon\nabla(\eta^{(N)})^{2r})\|_
{(L_{1}(0,T; L_{6/5}(\Omega)))^2}\leq c\,,\\
\|(\nabla(\rho^{(N)})^{2r-1}, \nabla(\eta^{(N)})^{2r-1})\|_
{(L_{1}(0,T; L_{3/2}(\Omega)))^2}\leq c.
\end{gathered}
\end{equation}
Especially, if $r\geq 9/5$, then
\begin{equation} \label{3.7}
\begin{gathered}
\|(|\varepsilon^{3/2}\Delta\rho^{(N)}\nabla\rho^{(N)}|, |\varepsilon^{3/2}\Delta\eta^{(N)}\nabla\eta^{(N)}|)\|_
{(L_{1}(Q))^2}\leq c,\\
\|(\varepsilon\Delta\rho^{(N)}\rho^{(N)}, \varepsilon\Delta\eta^{(N)}\eta^{(N)})\|_
{(L_{1}(Q))^2}\leq c,\\
\|(\nabla(\rho^{(N)})^{2r}, \nabla(\eta^{(N)})^{2r})\|_
{(L_{1}(Q))^2}\leq c.
\end{gathered}
\end{equation}
\end{lemma}

\begin{proof}
We estimate only $\rho$ because the same estimates
hold also for $\eta$. By the Sobolev's imbedding theorem and
Lemma \ref{lem3.1}, it holds
\[
\|\rho^{(N)}\|_{L_{2r-1}(0,T;L_{6r-3}(\Omega))}^{2r-1}
=\|(\rho^{(N)})^{r-1/2}\|_{L_2(0,T;L_{6}(\Omega))}^2\leq c.
\]
We apply H\"{o}lder's inequality to obtain
\begin{gather*}
\begin{aligned}
& \|\varepsilon\nabla(\rho^{(N)})^{2r}\|_{L_{1}(0,T; L_{6/5}(\Omega))}
\leq 2r\|(\rho^{(N)})^{2r-1}\|_{L_{1}(0,T; L_{3}(\Omega))}
\|\varepsilon\nabla\rho^{(N)}\|_{L_{\infty}(0,T; L_2(\Omega))}\\
&=2r\|\rho^{(N)}\|_{L_{2r-1}(0,T;L_{6r-3}(\Omega))}^{2r-1}
\|\varepsilon\nabla\rho^{(N)}\|_{L_{\infty}(0,T; L_2(\Omega))}\leq c,
\end{aligned}\\
\begin{aligned}
&\|\nabla(\rho^{(N)})^{2r-1}\|_{L_{1}(0,T; L_{3/2}(\Omega))} \\
&\leq 2\|\nabla(\rho^{(N)})^{r-1/2}\|_
{L_2(Q)}\|(\rho^{(N)})^{r-1/2}\|_ {L_2(0,T; L_{6}(\Omega))}
\leq c.
\end{aligned}
\end{gather*}
In the  $r\geq 9/5$, using the interpolation
$\alpha=\frac{4r}{10r-3}$,
\[
[{L_{\infty}(0,T;L_{2r}(\Omega))}, {L_{2r-1}(0,T;L_{6r-3}(\Omega))}]_{\alpha}=L_{\frac{10r-3}{3}}(Q)\subseteq
L_{5}(Q)
\]
and Lemma \ref{lem3.1}, we have
\begin{align*}
 \||\sqrt{\varepsilon\rho^{(N)}}\nabla\rho^{(N)}|\|_{L_2(Q)}^2
&=\frac{1}{2}\int_0^{T}\int_{\Omega}\varepsilon\nabla\rho^{(N)}
 \cdot\nabla((\rho^{(N)})^2-\rho_D^2)\,dx\,dt\\
&=-\frac{1}{2}\int_0^{T}\int_{\Omega}\varepsilon
\frac{\Delta\rho^{(N)}}{\sqrt{\rho^{(N)}}}\sqrt{\rho^{(N)}}
((\rho^{(N)})^2-\rho_D^2)\,dx\,dt\leq c
\end{align*}
which yields
\[
\|\varepsilon^{3/2}\Delta\rho^{(N)}\nabla\rho^{(N)}\|_{L_{1}(Q)}
=\|\varepsilon\frac{\Delta\rho^{(N)}}{\sqrt{\rho^{(N)}}}
\sqrt{\varepsilon\rho^{(N)}}\nabla\rho^{(N)}\|_{L_{1}(Q)}
\leq c.
\] By  Lemmas \ref{lem3.1} and \ref{lem3.2} it follows the second estimate
of \eqref{3.7},
\[
\|\varepsilon\Delta\rho^{(N)}\rho^{(N)}\|_
{L_{1}(Q)}=\|\varepsilon\frac{\Delta\rho^{(N)}}{\sqrt{\rho^{(N)}}}
(\rho^{(N)})^{3/2}\|_{L_{1}(Q)}\leq c.
\]
 The third estimate of \eqref{3.7} is immediate
consequence of Lemma \ref{lem3.1} and
\begin{gather*}
\nabla(\rho^{(N)})^{2r}=\frac{4r}{2r-1}(\rho^{(N)})^{r+1/2}
\nabla(\rho^{(N)})^{r-1/2},\\
\|(\rho^{(N)})^{r+1/2}\|_{L_{\frac{20r-6}{6r+3}}(Q)}
=\|\rho^{(N)}\|_{L_{\frac{10r-3}{3}}(Q)}^{\frac{2r+1}{2}}\leq
c, \\
\frac{20r-6}{6r+3}\geq 2.
\end{gather*}
\end{proof}

In the following we denote all constants dependent on
$\varepsilon\in(0, 1)$ as $c(\varepsilon)$.

\begin{lemma}\label{lem3.4}
There exists a constant $c(\varepsilon)>0$
independent $N=N^{*}, N^{*}+1,\dots$ such that
\begin{equation}\label{3.8}
\begin{gathered}
\|(\rho^{(N)}, \eta^{(N)})\|
_{(L_2(0,T; L_{\infty}(\Omega))\cap
L_{4/3}(0,T; H^{3/2-\delta}(\Omega)))^2}
\leq c(\varepsilon),\\
\|(\Delta\rho^{(N)}, \Delta\eta^{(N)})\|_{(L_{4/3}(0,T; L_2(\Omega)))^2}
\leq c(\varepsilon).
\end{gathered}
\end{equation}
where $\delta\in(0,  3/2)$.
\end{lemma}

\begin{proof}
By Lemma \ref{lem3.1} and \ref{lem3.2}, and Sobolev's embedding
theorem, we have
\begin{equation}\label{3.9}
\begin{aligned}
\|\Delta\rho^{(N)}\|_{L_2(0,T; L_{12/7}(\Omega))}^2
&\leq\int^{T}_0\Big(\int_{\Omega}|\frac{\Delta\rho^{(N)}}{\sqrt{\rho^{(N)}}}|^2dx\Big)
\Big(\int_{\Omega}|\rho^{(N)}|^{6}dx\Big)^{1/6}dt\\
&\leq c(\varepsilon).
\end{aligned}
\end{equation}
Here we use the embedding $H^{1}(\Omega)\hookrightarrow
L_{6}(\Omega)$ for the $d(\leq 3)$-dimensional domain $\Omega$ and
$\|\rho^{(N)}\|_{L_{\infty}(0,T;H^{1}(\Omega))}\leq c(\varepsilon)$.

Using this fact, let us prove that the set
$\{\rho^{(N)}:N=N^{*}, N^{*}+1,\dots\}$ is bounded in
$L_2(0,T; L_{\infty}(\Omega))$. We set
$Z_i:=\rho_i(\theta
h(\rho^2_i)+\tau\ln \rho^2_i-F_i-V_i)$
and take $\phi=(\rho_i-K)_{+}\in
H^{1}_0(\Omega\cup\Gamma_{N}), \;K\geq\rho_D$ as test function
of the second equation of \eqref{1.4}. Then we have
\[
\varepsilon^2\int_{\Omega}|\nabla(\rho_i-K)_{+}|^2dx
\leq \|Z_i\|_{L_{12/7}(\Omega)}\|(\rho_i-K)_{+}\|_{L_{6}(\Omega)}
|\Omega(\rho_i>K)|^{1/4}.
\]
Using the Sobolev's embedding theorem,  from the above
inequality we obtain
\begin{align*}
 (K'-K)|\Omega(\rho_i>K')|^{1/6}
&\leq\|(\rho_i-K)_{+}\|_{L_{6}(\Omega)}\\
&\leq c\varepsilon^{-2}\|Z_i\|_{L_{12/7}(\Omega)}
|\Omega(\rho_i>K)|^{1/4},  \quad K'>K
\end{align*}
which allows us to use Lemma \ref{lem2.1}. Hence it follows
\begin{equation}\label{3.10}
\begin{gathered}
\rho_i\leq\rho_D+c\varepsilon^{-2}\|Z_i\|_{L_{12/7}(\Omega)}
=\rho_D+c\|\Delta\rho_i\|_{L_{12/7}(\Omega)},\\
\|\rho^{(N)}\|_{L_2(0,T; L_{\infty}(\Omega))}^2
\leq\tau\sum^{N}_{i=1}\Big(\rho_D+c\|\Delta\rho_i\|_{L_{12/7}(\Omega)}\Big)^2
\leq c(\varepsilon).
\end{gathered}
\end{equation}
Also, by \eqref{3.10} and Lemma \ref{lem3.1} we obtain
\begin{equation}\label{3.11}
\begin{aligned}
\|\Delta\rho^{(N)}\|_{L_{4/3}(0,T; L_2(\Omega))}^{4/3}
&\leq\int^{T}_0\Big(\int_{\Omega}|
 \frac{\Delta\rho^{(N)}}{\sqrt{\rho^{(N)}}}|^2dx\Big)^{2/3}
\|\rho^{(N)}(t)\|^{2/3}_{L_{\infty}(\Omega)}dt \\
&\leq c(\varepsilon).
\end{aligned}
\end{equation}
Furthermore, by  \cite[Theorem 1]{s1} it holds
\[
\|\rho_i\|_{H^{3/2-\delta}(\Omega)}\leq c(\|\rho_i\|_{H^{1}(\Omega)}
+\|Z_i\|_{L_2(\Omega)}+\|\rho_D\|_{L_2(\Omega)}).
\]
which with \eqref{3.11} implies
\begin{equation}\label{3.12}
\|\rho^{(N)}\|_{L_{4/3}(0,T;H^{3/2-\delta}(\Omega))}\leq
c(\varepsilon).
\end{equation}
We can obtain similar estimates to \eqref{3.10}-\eqref{3.12} for
$\eta$.
\end{proof}

\section{Convergence}


\begin{lemma}\label{lem4.1}
 For any fixed $\varepsilon\in(0, 1)$ the set
$\{(\rho^{(N)},\eta^{(N)}); N=N^{*},N^{*}+1,\dots\}$ is precompact
in $(L_{p}(0,T; H^{1}(\Omega)))^2$ for all $p\in(1,\infty)$.
\end{lemma}

\begin{proof} 
 Dividing the first equation in \eqref{1.4} by
$\rho_i+a(\tau)$, we have
\begin{equation}\label{4.1}
\widehat{\rho^{(N)}}=\frac{1}{2(\rho^{(N)}+a(\tau))}
\operatorname{div}\Big((\rho^{(N)}+a(\tau))^2\nabla
F^{(N)}\Big)+g^{(N)}
\end{equation}
where 
\[
\widehat{\rho^{(N)}}(x,t)=\frac{\rho_i-\rho_{i-1}}{\tau}, \quad
g^{(N)}(x,t)=\frac{(\rho_i-\rho_{i-1})^2}{2\tau(\rho_i+a(\tau))}, \quad
t\in(t_{i-1}, t_i].
\]
By Lemma \ref{lem3.1},
\begin{equation}\label{4.2} 
\|g^{(N)}\|_{L_{1}(Q)}\leq c.
\end{equation}
We estimate the first term of the right-hand side in \eqref{4.1} as
\begin{align*}
& \frac{1}{2(\rho^{(N)}+a(\tau))}
\operatorname{div}\Big((\rho^{(N)}+a(\tau))^2\nabla F^{(N)}\Big)\\
&=\frac{1}{2}\operatorname{div}((\rho^{(N)}+a(\tau))\nabla F^{(N)})
+\tau\frac{|\nabla\rho^{(N)}|^2}{\rho^{(N)}}-\frac{1}{2}\nabla\rho^{(N)}\cdot\nabla V^{(N)}\\
&\quad -\frac{\varepsilon^2}{2}\nabla\rho^{(N)}
\cdot\nabla\Big(\frac{\Delta\rho^{(N)}}{\rho^{(N)}}\Big)
+r\theta(\rho^{(N)})^{2r-3}|\nabla\rho^{(N)}|^2\\
&=\sum^{5}_{i=1}R_i^{(N)}.
\end{align*}
By Lemmas \ref{lem3.1} and \ref{lem3.2},  for some $c(\varepsilon)>0$ it holds
\[
\|R_{1}^{(N)}\|_{L_2(0,T;H^{-1}(\Omega))}, \;
\|R_2^{(N)}\|_{L_{1}(Q)},
\|R_{3}^{(N)}\|_{L_{\infty}(0,T;L_{1}(\Omega))}, \;
\|R_{5}^{(N)}\|_{L_{1}(Q)}\leq c(\varepsilon).
\] 
Now, let us estimate $R_{4}$. Lemma\ref{lem3.1} and \ref{lem3.2}, and the equality
\[
\rho^{(N)}\nabla F^{(N)}
=-\varepsilon^2\rho^{(N)}\nabla(\frac{\Delta\rho^{(N)}}{\rho^{(N)}})
+2r\theta(\rho^{(N)})^{2(r-1)}\nabla\rho^{(N)}-\rho^{(N)}\nabla
V^{(N)}
\] 
yield the estimate:  for some $c(\varepsilon)>0$,
\[
\|\rho^{(N)}\nabla(\frac{\Delta\rho^{(N)}}{\rho^{(N)}})\|_{L_{1}(Q)}\leq
c(\varepsilon).
\] 
Hence  Lemma \ref{lem3.2} and \ref{lem3.4}, and the equality
\begin{align*}
\nabla\rho^{(N)}\cdot\nabla(\frac{\Delta\rho^{(N)}}{\rho^{(N)}})
&=\operatorname{div}(\nabla\rho^{(N)}\frac{\Delta\rho^{(N)}}{\rho^{(N)}})
-\frac{(\Delta\rho^{(N)})^2}{\rho^{(N)}}\\
&=\operatorname{div}(\nabla(\Delta\rho^{(N)})-\rho^{(N)}\nabla(\frac{\Delta\rho^{(N)}}{\rho^{(N)}}))
-\frac{(\Delta\rho^{(N)})^2}{\rho^{(N)}}
\end{align*}
give the estimate: for some $c(\varepsilon)>0$,
\[ 
\|R_{4}^{(N)}\|_{L_{1}(0,T; Z)}\leq c(\varepsilon),
\]
where $Z=H^{-2}(\Omega)$. Thus we obtain the estimate
\begin{equation}\label{4.3}
\frac{1}{\tau}\|\rho^{(N)}(\cdot)-\rho^{(N)}(\cdot-\tau)\|_{L_{1}(0,T;Z)}
=\|\widehat{\rho^{(N)}}\|_{L_{1}(0,T;Z)}
\leq c(\varepsilon).
\end{equation}
This fact and Lemma \ref{lem3.4} allows us to use the compactness theorem for
piecewise constant functions \cite[Theorem 1]{dj} to derive the
compactness of $\{\rho^{(N)}; N=N^{*}, N^{*}+1,\dots\}$ in
$L_{p}(0,T; H^{1}(\Omega)), \;\forall p\in(1,\infty)$.
We can also prove the compactness of $\{\eta^{(N)};
N=N^{*}, N^{*}+1,\dots\}$ similarly.
\end{proof}


\begin{proof}[Proof of Theorem \ref{thm1.2}]
 Using Lemmas \ref{lem3.2} and \ref{lem3.4}, and \eqref{3.6} of 
Lemmas \ref{lem3.3} and  \ref{lem4.1},  we can easily verify the
convergence estimate \eqref{1.8} for some $\rho, \eta, V$.

It remains only to prove that $\rho, \eta, V$ satisfy
\eqref{1.9}. From Theorem \ref{thm1.1} and Lemma \ref{lem3.2} we have
\begin{equation}\label{4.4}
\begin{gathered}
 \frac{1}{\tau}(n^{(N)}-\sigma_{N}n^{(N)})
=\operatorname{div}((\rho^{(N)}+a(\tau))^2\nabla F^{(N)}),\\
\begin{aligned}
& \|\frac{1}{\tau}(n^{(N)}-\sigma_{N}n^{(N)})\|_{L_2(0,T; W^{-1}_{r+1}
(\Omega\cup\Gamma_{N}))}\\
&=\|(\rho^{(N)}+a(\tau))^2\nabla
F^{(N)}\|_{L_2(0,T; L_{(r+1)/r}(\Omega))}\leq c
\end{aligned}
\end{gathered}
\end{equation}
for some constant $c>0$ independent of $N=N^{*}, N^{*}+1, \dots$
and $\varepsilon\in(0,1)$ where $n^{(N)}$ is defined as
\[
n^{(N)}(x,t)=(\rho_i+a(\tau))^2,t\in(t_{i-1}, t_i]
\] 
and  $\sigma_{N}$ is the shift operator
\[
\sigma_{N}n^{(N)}(x,t)=(\rho_{i-1}+a(\tau))^2,
t\in(t_{i-1}, t_i].
\] 
Since it holds by the strong convergence of
$\rho^{(N)}$ to $\rho$
\[
n^{(N)}\to\rho^2\quad  \text{in }L_2(Q)
\]  
as $N\to\infty$, we can easily verify in
the sense of distribution $[0,T]\to W^{-1}_{r+1}
(\Omega\cup\Gamma_{N})$,
\[
\frac{1}{\tau}\Big(n^{(N)}-\sigma_{N}n^{(N)}\Big)\to
 \frac{\partial\rho^2}{\partial t}.
\] 
This and \eqref{4.4} gives us, up to a subsequence,
\begin{equation}\label{4.5}
\frac{1}{\tau}\Big(n^{(N)}-\sigma_{N}n^{(N)}\Big)\rightharpoonup
\frac{\partial\rho^2}{\partial t} \quad  \text{weakly in } 
L_2(0,T; W^{-1}_{r+1} (\Omega\cup\Gamma_{N})).
\end{equation}

We take  $\phi\in C^{\infty}_0(Q)$ as test function of the first
equation in \eqref{4.4} to obtain
\begin{align}
&\int^{T}_0 \langle \frac{1}{\tau}(n^{(N)}-\sigma_{N}n^{(N)}),\phi\rangle\,dt \nonumber\\
&=-\varepsilon^2\int_{Q}\frac{\Delta\rho^{(N)}}{\rho^{(N)}}
\operatorname{div}((\rho^{(N)}+a(\tau))^2\nabla\phi)\,dx\,dt \nonumber\\
&\quad -\theta\int_{Q}(\rho^{(N)}+a(\tau))^2
\nabla h((\rho^{(N)})^2)\cdot\nabla\phi \,dx\,dt \nonumber\\
&\quad -2\tau\int_{Q}(\rho^{(N)}+a(\tau))^2\frac{\nabla\rho^{(N)}}{\rho^{(N)}}
\cdot\nabla\phi \,dx\,dt
+\int_{Q}(\rho^{(N)}+a(\tau))^2\nabla V^{(N)}\cdot\nabla\phi \,dx\,dt \nonumber\\
&=-2\varepsilon^2\int_{Q}\Delta\rho^{(N)}\nabla\rho^{(N)}\cdot\nabla\phi
\,dx\,dt -\varepsilon^2\int_{Q}\Delta\rho^{(N)}\rho^{(N)}\Delta\phi \,dx\,dt 
\nonumber\\
&\quad -\theta\int_{Q}\nabla(\rho^{(N)})^{2r}\cdot\nabla\phi \,dx\,dt
 +\int_{Q}(\rho^{(N)})^2\nabla V^{(N)}\cdot\nabla\phi \,dx\,dt \nonumber\\
&\quad -2\tau\int_{Q}\rho^{(N)}\nabla\rho^{(N)}\cdot\nabla\phi \,dx\,dt
 -a(\tau)^2\int_{Q}\nabla F^{(N)}\cdot \nabla\phi \,dx\,dt \nonumber\\
&\quad -2a(\tau)\int_{Q}\Big(\varepsilon^2\frac{\Delta\rho^{(N)}}{\rho^{(N)}}
 \nabla\rho^{(N)}\cdot\nabla\phi +\varepsilon^2\Delta\rho^{(N)}
 \Delta\phi\Big)\,dx\,dt\nonumber\\
&\quad -2a(\tau)\int_{Q}\Big(\frac{2\theta r}{2r-1}\nabla(\rho^{(N)})^{2r-1}
+2\tau\nabla\rho^{(N)}-\rho^{(N)}\nabla
V^{(N)}\Big)\cdot\nabla\phi \,dx\,dt \nonumber\\
&=-2\varepsilon^2\int_{Q}\Delta\rho^{(N)}\nabla\rho^{(N)}\cdot\nabla\phi
\,dx\,dt -\varepsilon^2\int_{Q}\Delta\rho^{(N)}\rho^{(N)}\Delta\phi \,dx\,dt 
\nonumber\\
& \quad -\theta\int_{Q}\nabla(\rho^{(N)})^{2r}\cdot\nabla\phi
\,dx\,dt+\int_{Q}(\rho^{(N)})^2\nabla V^{(N)}\cdot\nabla\phi \,dx\,dt \nonumber\\
&\quad -2\tau R_{1}^{N}-a(\tau)^2R_2^{N}-2a(\tau)R_{3}^{N}.  \label{4.6}
\end{align}
By Lemma \ref{lem3.2},
\begin{align*}
2\tau R_{1}^{N}=2\tau \int_{Q}\rho^{(N)}\nabla\rho^{(N)}\cdot\nabla\phi \,dx\,dt,\\
a(\tau)^2 R_2^{N}= a(\tau)^2\int_{Q}\nabla
F^{(N)}\cdot\nabla\phi \,dx\,dt
\end{align*}
approaches zero as $N\to\infty(\tau\to 0)$. 
Also, considering the property $a(\tau)\leq c\tau^2$  of $a(\tau)$ (see page 3) and
Lemmas \ref{lem3.1} and \ref{lem3.2}, \eqref{3.6} of Lemmas 
\ref{lem3.3} and  \ref{lem3.4}, it follows that
\begin{align*}
2a(\tau)R_{3}^{N}
&= 2a(\tau)\int_{Q}(\varepsilon^2\frac{\Delta\rho^{(N)}}{\sqrt{\rho^{(N)}}}
 \frac{\nabla\rho^{(N)}}{\sqrt{\rho^{(N)}}}\cdot\nabla\phi
+\varepsilon^2\Delta\rho^{(N)}\Delta\phi)\,dx\,dt\\
&\quad +2a(\tau)\int_{Q}(\frac{2\theta r}{2r-1}\nabla(\rho^{(N)})^{2r-1}
+2\tau\nabla\rho^{(N)}-\rho^{(N)}\nabla
V^{(N)})\cdot\nabla\phi \,dx\,dt \\
&\to 0 \quad\text{as }N\to \infty.
\end{align*}
Thus, using the obtained convergence  estimates \eqref{1.8} and
taking the limit $\tau\to 0$ in the above equation
\eqref{4.6}, we arrive at the first equation of \eqref{1.9}.
Similarly, we obtain the second  equation in \eqref{1.9} for $\eta$ .
The third equation of \eqref{1.9} can also be obtained by the limit
of
\[
\lambda^2\Delta V^{(N)}=(\rho^{(N)}+a(\tau))^2-(\eta^{(N)}+a(\tau))^2-f.
\]
\end{proof}


\begin{proof}[Proof of Theorem \ref{thm1.3}]
 Using \eqref{3.7} of lemmas \ref{lem3.2} and \ref{lem3.3}, 
\eqref{4.4} and \eqref{4.5},  we can easily verify for some $c>0$
\begin{equation} \label{4.7}
\begin{gathered}
\|(\frac{\partial n^{(\varepsilon)}}{\partial t}, 
\frac{\partial p^{(\varepsilon)}}{\partial t})
\|_{(L_2(0,T; W^{-1}_{r+1}(\Omega\cup\Gamma_{N})))^2}\leq c,
\\
\|(\nabla (n^{(\varepsilon)})^{r}, \nabla (p^{(\varepsilon)})^{r})
\|_{(L_{1}(Q))^2}\leq c, \\
\|V^{(\varepsilon)}\|_{L_{\infty}(0,T; H^{1}(\Omega))}\leq c, \quad
  \forall \varepsilon\in(0,1).
\end{gathered}
\end{equation}
Hence the set $\{n^{(\varepsilon)},
p^{(\varepsilon)}; \varepsilon\in(0,1)\}$ is precompact in
$L_{p}(0,T; L_{r}(\Omega))$ for all $p<\infty$ 
(cf. \cite[Theorem  3.2.2]{j1}, \cite[Theorem 6]{s}). 
This allows us to take the limit
$\varepsilon\to 0$ in \eqref{1.9} and arrive at \eqref{1.10}
and \eqref{1.12} with the help of Lemmas \ref{lem3.2} and  \ref{lem3.3}.
\end{proof}

\begin{thebibliography}{99}


\bibitem{au}
N. Ben Abdallah, A. Unterreiter;
\emph{On the stationary quantum drift-diffusion model}, 
Z. Angew. Math. Phys., 49 (1998), 251-275.

\bibitem{br} H. Brunner;
\emph{Methods for Volterra integral and related functional
equations}, Cambridge Univ., 2004. 

\bibitem{BCD-12} S. Bian, L. Chen, M. Dreher;
\emph{Boundary layer analysis in the semiclassical limit of a quantum 
drift-diffusion model}, J. Differential Equations, 253 (2012), 356-377.

\bibitem{cct} M. J. C\'{a}ceres, J. A. Carrillo, G. Toscani;
\emph{Long-time behavoir for a nonlinear fourth-order parabolic equation}, 
Trans. Amer. Math. Soc., 357 (2005), 1161-1175.

\bibitem{cj}  L. Chen, Q.-C. Ju;
\emph{Existence of weak solution and semiclassical limit for quantum 
drift-diffusion model}, Z. Angew. Math. Phys., 58 (2007), 1-15.

\bibitem{c} X. Chen;
\emph{The isentropic quantum drift-diffusion model in two or
three dimensions}, Z. Angew. Math. Phys., 60 (2009), 416-437.

\bibitem{cc} X. Chen, L. Chen;
\emph{Initial time layer problem for quantum drift-diffusion model}, 
J. Math. Anal. Appl. 343 (2008), 64-80.

\bibitem{ccj} X.-Q. Chen, L. Chen, H.-Y. Jian;
\emph{Existence, semiclassical limit and long-time behavior of weak solution 
to quantum drift-diffusion model}, Nonlinear Anal. Real World Appl., 10 (2009), 1321-1342.

\bibitem{ccj1} X.-Q. Chen, L. Chen, H.-Y. Jian;
\emph{The existence and long-time behavior of weak solution to bipolar 
quantum drift-diffusion model},
Chin. Ann. Math. B, 28 (2007), 651-664.

\bibitem{D-13} J. Dong;
\emph{Classical solutions to the one-dimensional stationary
quantum drift-diffusion model}, J. Math. Anal. Appl., 399 (2013), 594-598.

\bibitem{DW-15} J. Dong, W. Zhang;
\emph{On the stationary quantum drift-diffusion model
in one space dimension}, Adv. Math. China, 44 (2) (2015), 263-270.

\bibitem{dj} M. Dreher, A. J\"{u}ngel;
\emph{Compact families of piecewise constant functions in $L^{p}(0,T;B)$},
 Nonlinear Analysis, 75 (2012), 3072-3077.

\bibitem{gst} U. Gianazza, G. Savare, G. Toscani;
\emph{The Wasserstein gradient flow of Fisher information and the 
quantum drift-diffusion equation},
Arch. Ration. Mech. Anal., 194 (2009), 133-220.

\bibitem{gjt} M. P. Gualdani, A. J\"{u}ngel, G. Toscani;
\emph{A nonlinear fourth-order parabolic equations with non-homogeneous boundary
conditions}, SIAM. J. Math. Anal., 37 (2006), 1761-1779.

\bibitem{j1} A. J\"{u}ngel;
\emph{Quasi-hydrodynamlc semiconductor equations},
Birkh\"{a}user,  2001.

\bibitem{jlm} A. J\"{u}ngel, H.-L. Li, A. Matsumura;
\emph{The relaxation-time limit in the quantum hydrodynamic equations for
 semiconductors}, J. Diff. Equat., 225 (2006), 440-464.

\bibitem{jm} A. J\"{u}ngel,  D. Matthes;
\emph{The Derrida-Lebowitz-Speer-Spohn equation: existence, non-uniqueness, 
and decay rates of the solutions}, SIAM J. Math. Anal., 39 (2008), 1996-2015.

\bibitem{jp1} A. J\"{u}ngel, R. Pinnau;
\emph{Global non-negative solutions of a
nonlinear fourth-order parabolic equations for quantum systems},
SIAM. J. Math. Anal., 32 (2000), 760-777.

\bibitem{jp2} A. J\"{u}ngel, R. Pinnau;
\emph{A positivity-preserving numerical scheme
for a nonlinear fourth-order parabolic equations}, SIAM. J. Num.
Anal., 39 (2001), 385-406.

\bibitem{jv} A. J\"{u}ngel, I. Violet;
\emph{The quasineutral limit in the quantum drift-diffusion equations}, 
Asymptotic Anal., 53 (2007), 139-152.

\bibitem{J-09} Q.-C. Ju;
\emph{Semiclassical limit in the quantum drift-diffusion model},
Acta Mathematica Sinica, 25 (2) (2009), 253-264.

\bibitem{jc} Q.-C. Ju, L. Chen;
\emph{Semiclassical limit for bipolar quantum drift-diffusion model}, 
Acta Mathematica Scientia, 29B (2009), 285-293.

\bibitem{lzzh} H.-L. Li, G.-J. Zhang, M. Zhang, Chengchun Hao;
\emph{Long-time self-similar asymptotic of the macroscopic quantum models},
 J. Math. Phys., 49(2008), 57-86.

\bibitem{nss} S. Nishbata, N. Shigeta, M. Suzuki;
\emph{Asymptotic behaviors and classical limits of solutions to a quantum 
drift-diffusion model for semiconductors}, 
Math. Models and Mech. Appl. Sci., 20 (2010), 909-936.

\bibitem{s1} G. Savare;
\emph{Regularity and perturbation results for mixed second
order elliptic problems}, Comm. Part. Diff. Equat., 22 (1997), 860-899.

\bibitem{s} J. Simon;
\emph{Compact sets in the space}, Ann. Mat. Pura Appl., 146 (1987), 65-96.

\bibitem{x} X. Xu;
\emph{Existence and semiclassical limit for the quantum
drift-diffusion model}, Annali di Matematica, 193 (2012), 889-908.

\bibitem{zlz} G.-J. Zhang, H.-L. Li, K. Zhang;
\emph{Semiclassical and relaxation limit
of bipolar quantum hydrodynamic model for semiconductors}, 
J. Diff. Equat., 245 (2008), 1433-1453.

\end{thebibliography}

\end{document}