\documentclass[twoside]{article}
\usepackage{amssymb} % font used for R in Real numbers
\pagestyle{myheadings}
\markboth{\hfil Solutions to a nonlinear drift-diffusion model \hfil EJDE--1999/15}
{EJDE--1999/15\hfil Weifu Fang \& Kazufumi Ito \hfil}
\begin{document}
\title{\vspace{-1in}\parbox{\linewidth}{\footnotesize\noindent
{\sc  Electronic Journal of Differential Equations},
Vol. {\bf 1999}(1999), No.~15, pp. 1--38. \newline
ISSN: 1072-6691. URL: http://ejde.math.swt.edu or http://ejde.math.unt.edu
\newline ftp  ejde.math.swt.edu \quad ftp ejde.math.unt.edu (login: ftp)}
 \vspace{\bigskipamount} \\
%
 Solutions to a nonlinear drift-diffusion model for semiconductors 
\thanks{ {\em 1991 Mathematics Subject Classifications:} 35K57, 35K55, 35J60, 78A35.
\hfil\break\indent
{\em Key words and phrases:} Drift-diffusion model, semiconductors, 
  nonlinear diffusion, \hfill\break\indent
degenerated parabolic and elliptic equations, attractors.
\hfil\break\indent
\copyright 1999 Southwest Texas State University  and University of
North Texas. \hfil\break\indent
Submitted May 28, 1998. Published May 10, 1999. \hfill\break\indent
(WF) partially supported by WVU Senate Research Grant R-97-051 \hfill\break\indent
and by Army Research Office grant DAAG55-98-1-0261. \hfill\break\indent     
(IK) paritally supported by Air Force Office gant AFOSR-F49620-95-1-0447. 
} }
\date{}
%
\author{Weifu Fang \& Kazufumi Ito}
\maketitle

\begin{abstract} 
A nonlinear drift-diffusion model for semiconductors is analyzed to show the
existence of non-vacuum global solutions and stationary solutions.  The long
time behavior of the solutions is studied by establishing the existence of an
absorbing set and a compact attractor of the dynamical system.  Parallel results
on vacuum solutions are also obtained under weaker conditions on model
parameters.  
\end{abstract}

\renewcommand{\theequation}{\thesection.\arabic{equation}}    
\newtheorem{theorem}{Theorem}[section]    
\newtheorem{lemma}[theorem]{Lemma}    
\newtheorem{corollary}[theorem]{Corollary}    
\newtheorem{remark}[theorem]{Remark} 
\newcommand{\ds}{\displaystyle}    
\newcommand{\ep}{\varepsilon}    
   
\section{Introduction}     
\setcounter{equation}{0}    

In this paper we analyze a nonlinear drift-diffusion model for semiconductors.
Let $n=n(t,x)$ and $p=p(t,x)$ be the densities of electrons and holes, 
respectively, in a region $\Omega\subset {mathbb R}^l$ ($l\le 3$) occupied 
by a semiconductor. Consider the flow of the
charge carriers in $\Omega$. Conservation of mass yields
\begin{equation}\label{np}
\frac{\partial n}{\partial t}=G_n+\mbox{div }J_n \quad\mbox{and}\quad
\frac{\partial p}{\partial t}=G_p-\mbox{div }J_p
\end{equation}
where $J_n$ and $J_p$ are the current densities, and $G_n$ and $G_p$ are 
generation-recombination rates for electrons and holes, respectively. 

In the standard drift-diffusion model, temperature is assumed to be constant $T_0$ 
(the lattice temperature) and thus the current densities are given by
\begin{equation} \label{current0}
J_n=  \mu_n (\nabla (k_b n T_0)- n\nabla u)\;\;\;\mbox{and}\;\;\;
J_p= - \mu_p(\nabla (k_b p T_0)+ p\nabla u)
\end{equation}
where $\mu_n$, $\mu_p$ are the mobilities, $k_b$ is the Boltzmann constant, 
and the Einstein relations that relate the diffusion coefficients to the mobilities 
are already assumed. The electric potential $u=u(t,x)$ is governed by the Poisson's equation
\begin{equation}\label{u}
\nabla\cdot (\epsilon_0\nabla u)=n-p-N 
\end{equation}
with electrical permitivity $\epsilon_0$ and impurity doping profile $N(x)$.
Noting that $P_n=k_b n T_0$ and $P_p=k_b p T_0$ are the pressures, if
the electron and hole temperatures, $T_n$ and $T_p$, are variables, then
the pressures should be expressed as
$$P_n=k_b n T_n \quad\mbox{and}\quad P_p=k_b p T_p.$$
Using an isentropic assumption analogous to the case of gas dynamics, we further 
assume that $P_n$ and $P_p$ are functions of $n$ and $p$ only, i.e., 
\begin{equation}\label{isen}
k_b T_n n=P_n(n)\quad\mbox{and}\quad k_b T_p p=P_p(p)
\end{equation}
where $P_n$ and $P_p$ are given functions. For example, in the case of polytropic gas, 
$P_n=k n^\gamma$ $(\gamma>1)$. 
For simplicity we will assume
$$P_n(\cdot)=P_p(\cdot)=P(\cdot).$$
Thus the current relations in (\ref{current0}) should be modified to become
\begin{equation}\label{current}
J_n=\mu_n( \nabla P(n)- n\nabla u)\;\;\;\mbox{and}\;\;\;
J_p=-\mu_p( \nabla P(p)+ p\nabla u).
\end{equation}

Here we obtain this model by using the analogy to gas dynamics,
but such current expressions are also known as results of the so-called 
Fermi-Dirac statistics where $P(n)$ is related to a Fermi integral and 
approximately equal to $k n^{5/3}$ when $n$ is large. 
The standard model (\ref{current0}) is the special case of (\ref{current})
when $P(\cdot)$ is linear, which corresponds to the Boltzmann statistics
(see, e.g., \cite{MRS}).

For the generation-recombination rates $G_n$ and $G_p$ for the carriers, it is 
often assumed that 
$$G_n=G_p=G$$
since electrons and holes are generated or combined in pairs. 
This rate $G$ is determined by certain deviation of the density 
levels from the equilibrium level. The common form used in the literature
for the standard model (\ref{current0}) is
\begin{equation}\label{G0}
G=g-Q(n,p)(np-1)
\end{equation}
with different models for the coefficients $Q=Q(n,p)$. For example, in the Shockley-Read-Hall
model, $Q_{\rm SRH}\sim (1+n+p)^{-1}$, and in the Auger model, $Q_{\rm AU}\sim n+p$ (see \cite{MRS,J0}).
The nonnegative function $g$ in (\ref{G0}) represents a distributed generation rate applied to 
the system (see \cite{BFI} for an application). 
For the nonlinear model (\ref{current}), the equilibrium state is given by
$$H(n)+H(p)=h_i$$
where $H(\cdot)$ is the enthalpy function
\begin{equation}\label{enthalpy}
H(s)=\int_1^s \frac{P'(s)}{s} ds
\end{equation}
associated with the pressure $P(\cdot)$, and $h_i$ is a material constant. 
This equilibrium relation is derived from the zero quasi-Fermi potential 
conditions (see \cite{MU}), 
and includes the standard case $np=1$ since $H(s)=\ln s$ when $P$ is linear. 
Therefore we assume that the generation-recombination rate has the form 
\begin{equation}\label{G}
G=g-Q(n,p)(S(H(n)+H(p))-1)
\end{equation}
where the scaling function $S$ is nonnegative, nondecreasing with $S(h_i)=1$. 
When $P(\cdot)$ is linear and $S$ is exponential,
the above model becomes the common model (\ref{G0}).
But when $P$ is nonlinear, (\ref{G0}) cannot be included as a special case of (\ref{G}).
We will include both models (\ref{G0}) and (\ref{G}) in our study here. 
For (\ref{G}), the two limits of $H(\cdot)$ also play an important role in our analysis:
\begin{equation}\label{h0-infty}
h_0=-H(0^+)=\int_0^1\frac{P'(s)}{s}ds \quad\mbox{and}\quad
h_\infty=H(\infty)=\int_1^\infty\frac{P'(s)}{s}ds.
\end{equation}

Thus our model equations read
\begin{equation}\label{n}
\frac{\partial n}{\partial t}-\mbox{div }\mu_n(\nabla P(n)- n\nabla u)=G(n,p)
\end{equation}
\begin{equation}\label{p}
\frac{\partial p}{\partial t}-\mbox{div }\mu_p(\nabla P(p)+p\nabla u)=G(n,p)
\end{equation}
where $u$ satisfies (\ref{u}). 

The boundary conditions are the same as for standard drift-diffusion model
(see, e.g., \cite{MRS,J0}). Let $\partial\Omega=\Sigma_N\cup\Sigma_D$, 
and $\nu$ the outward normal on $\partial\Omega$. $\Sigma_N$ represents the insulated portion 
of the boundary, and $\Sigma_D$ the contact portion. Thus we have the
following boundary conditions:
\begin{equation}\label{bc}
\begin{array}{cl}
 n=\bar n(x),\;\; p=\bar p(x),\;\;u=\bar u(x) & \mbox{on}\;\;\Sigma_D\\[5pt]
{\ds\frac{\partial n}{\partial\nu}
=\frac{\partial p}{\partial\nu}
=\frac{\partial u}{\partial\nu}=0} & \mbox{on}\;\;\Sigma_N.
\end{array}
\end{equation}
Initial conditions $(n^0,p^0)$ are prescribed for $(n,p)$.
For compatibility, we also assume that $u^0$ is given by (\ref{u}) at $t=0$.

There have been some mathematical studies devoted to the analysis of 
the standard drift-diffusion model since the 1970's. 
For example, existence of stationary solutions has been shown 
(e.g., \cite{J,BFI}) by using decoupling mappings and fixed point theorems,
and local and global solutions for the evolutionary problem have been 
established (e.g., \cite{GG,ST,JDE1,JDE2,JDE3}).
We refer the interested reader to the books \cite{MRS,J0} for more complete
references. As for this nonlinear model under consideration, very few studies
can be found in the literature. Gajewski and Gr\"{o}ger \cite{GG2}
study the time-dependent model based on the Fermi-Dirac statistics with the 
assumptions that $S$ is the exponential function and $h_0=h_\infty=\infty$.
In the common example $P(s)=s^\gamma$ with $\gamma>1$,
$h_0=\gamma/(\gamma-1)$ is finite. In \cite{MU} Markowich and Untererreiter study 
the stationary problem for the case $h_0$ is finite and $G=0$, 
and show the existence of positive solutions under some smallness assumptions on the data,
and in general the existence of nonnegative solutions (vacuum solutions).  
Similar results on nonnegative solutions for the time-dependent problem are obtained by 
J\"{u}ngel \cite{Ju} for a more general model, but it does not include the generation-recombination 
rates like (\ref{G}). It is also worth noting that this model can be derived from
the hydrodynamic model when the equations of momentum conservation are reduced to 
the current relations (\ref{current}) by neglecting the convective terms (see \cite{RO}).
In \cite{MN}, it is shown rigorously in the case of one space dimension that such a model 
is the limit of the hydrodynamic model with the isentropic assumption as the relaxation time
tends to zero. Models of this type are known in the literature
as electro-diffusion systems (\cite{R}), and appear in other applications such as 
electrophoresis in electro-chemistry (\cite{CL}) and ion flow through
membrane channels in physiology (\cite{B,PJ}).

In this paper, we study the model (\ref{n})-(\ref{p}) with (\ref{u}) and 
boundary conditions (\ref{bc}) for generation-recombination rates given by either (\ref{G0}) 
or (\ref{G}). In particular, we will establish the existence of global non-vacuum solutions for the 
time-dependent problem and non-vacuum stationary solutions. 
Under some general assumptions (Assumption (A) below), we show 
the existence of non-vacuum global solutions when $G$ is from the common model (\ref{G0}).
For $G$ given by (\ref{G}), we need to require a relation (see (\ref{new}) below) among 
the constants $h_0$, $h_\infty$ and $h_i$ for the positive lower bound 
for $(n,p)$; this relation is always satisfied if $h_0=\infty$.
If we allow vacuum solutions (nonnegative solutions), then this restriction on the data
can be removed. We will also study the long-time behavior of the system, showing the
existence of an absorbing set and a compact global attractor. We use the following examples
to illustrate some of our results in this paper.

\vspace{.05in}

\noindent{\bf Example 1.} Consider
$$P(s)=s^{\gamma} \quad\mbox{and}\quad G(n,p)=g(x)-\bar{q}\frac{np-1}{n+p+1}$$
where $\gamma>1$ and $\bar q>0$ (the Shockley-Read-Hall recombination model).
Then, from any positive initial data $(n^0, p^0)$, 
the system (\ref{n})-(\ref{p}) with (\ref{u}) has a global positive (non-vacuum) 
solution $(n,p)$ in $H^1\cap L^\infty$ (Theorem \ref{thm:exist}), 
and there is a positive stationary solution (Corollary \ref{cor:stat}).
As a dynamical system, 
(\ref{n})-(\ref{p}) has an absorbing set 
$${\cal B}=\{ (n, p)\in L^\infty(\Omega)^2: 0<\underline{c}\le n(x), p(x)\le\overline{c}\}$$
(Theorem \ref{thm:absorb}). In some cases (in particular, the one-dimensional 
$\Omega$ case) the global solution is unique (Theorem \ref{thm:unique}), 
and the dynamical system has a compact attractor $\cal A$ that attracts 
all sets that are bounded in terms of the $L^\infty$-norm 
(Theorem \ref{thm:attract}).

\vspace{.05in}

\noindent{\bf Example 2.} For the case
$$P(s)=s^{\gamma} \quad\mbox{and}\quad G(n,p)=g(x)-\bar{q}(\exp( n^{\gamma-1}+p^{\gamma-1})-1)$$
with $\gamma>1$ and $\bar q>0$, beginning with any non-negative initial 
data $(n^0, p^0)$, the system has a global solution and a stationary 
solution, all in $L^\infty$ with $n$, $p\ge 0$, 
and $n^\gamma$, $p^\gamma$, $u\in H^1(\Omega)$, and
there is an absorbing set as in Example 1 above but with $\underline{c}=0$
(Theorem \ref{thm:vac}).

\vspace{.05in}

Throughout this paper we make the following general assumptions on the model parameters:

\begin{description}\itemsep=-3pt
\item[{Assumption (A):}] { }
\item[{\rm (i)}] {\em $P(\cdot)$ is from $[0,\infty)$ to $[0,\infty)$, $P(0)=0$,
 and $P'(\cdot)>0$ is locally Lipschitz continuous in $(0,\infty)$.} 
\item[{\rm (ii)}] {\em $\epsilon_0$, $\mu_n$ and $\mu_p$ are positive constants.}
\item[{\rm (iii)}] {\em $\Omega$ is an open, connected, bounded domain in ${mathbb R}^l$ 
 ($l=1, 2$ or $3$) of class $C^{0,1}$, and $\partial\Omega=\Sigma_D\cup\Sigma_N$ with
$\Sigma_D$ nonempty and relatively closed. } 
\item[{\rm (iv)}] {\em $N$, $g\in L^\infty(\Omega)$ with $g\ge 0$ $a.e.$ in $\Omega$. }
\item[{\rm (v)}] {\em $Q=Q(n,p)$ is locally Lipschitz continuous and 
$0\le \frac{Q(n,p)}{1+n+p} \le \overline Q$ (constant) for all $n$, $p\ge 0$. }
\item[{\rm (vi)}] {\em $S$ is locally Lipschitz continuous, nonnegative, nondecreasing on its domain 
 Dom$(S)\supseteq (-2h_0, 2h_\infty)$, and $S(h_i)=1$. }
\item[{\rm (vii)}] {\em $\bar n$, $\bar p$ and $\bar u\in H^1(\Omega)\cap L^\infty(\Sigma_D)$
 with $\inf_{\Sigma_D}\bar n>0$,  $\inf_{\Sigma_D}\bar p>0$. }
\item[{\rm (viii)}] {\em $n^0$, $p^0\in L^\infty(\Omega)$ with $\inf_{\Omega}n^0>0$, $\inf_{\Omega}p^0>0$. }
\end{description}

For the model (\ref{G0}), the condition (vi) on $S$ is not needed.
Under (v) on $Q$ above, (\ref{G0}) includes 
the most commonly used Shockley-Read-Hall model 
$Q_{\rm SRH}\sim (1+n+p)^{-1}$ and the Auger model $Q_{\rm AU}\sim n+p$. 
See e.g., \cite{MRS,J0}.

Define 
$$V=\{\phi\in H^1(\Omega):\;\phi=0\;\;\mbox{on}\;\;\Sigma_D\},$$
which is the space known as $H^1_0(\Omega\cup\Sigma_N)$ (see, e.g., \cite{Tr}).
Denote the usual $L^2$ product as $\langle\phi,\psi\rangle$. $V$ is
equipped with the $H^1$ equivalent norm $|\nabla\cdot|_2$. 
By the Poincar\'{e} inequality, for $2\le r\le 6$, 
there exists constant $\alpha_r>0$ such that
\begin{equation}\label{poincare}
|\phi|_r\le\alpha_r |\nabla\phi|_2\;\;\;\mbox{for all }\;\phi\in V.
\end{equation}
Moreover, the embedding of $V$ into each $L^r(\Omega)$ is compact.

Given $T>0$, a triplet $(n, p, u)$ is said to be a weak solution to the 
model equations (\ref{n})--(\ref{p}), (\ref{u}) with (\ref{bc}) if 
$$(n, p, u)\in (\bar n, \bar p, \bar u)+L^2(0,T; V)^3\cap H^1(0,T; V^\ast)^3$$
with $(n,p)=(n^0,p^0)$ at $t=0$, $n>0$ and $p>0$,  
and the triplet $(n,p,u)$ satisfies
\begin{equation}\label{w-n}
\langle\frac{\partial n}{\partial t},\phi\rangle_{V^\ast\times V}
+\mu_n\langle \nabla P(n)-n\nabla u,\nabla\phi\rangle =\langle G(n,p), \phi\rangle
\end{equation}
\begin{equation}\label{w-p}
\langle\frac{\partial p}{\partial t},\psi\rangle_{V^\ast\times V}
 +\mu_p\langle \nabla P(p)+p\nabla u,\nabla\psi\rangle = \langle G(n,p),\psi\rangle
\end{equation}
\begin{equation}\label{w-u}
\epsilon_0\langle\nabla u,\nabla\chi\rangle+\langle n-p-N, \chi\rangle=0
\end{equation}
for all $(\phi, \psi, \chi)\in V^3$, and a.e. in $(0,T)$. 

As in the case of the standard drift-diffusion model, it is
more convenient to use a set of new variables, the so-called
quasi-Fermi potentials, in place of $(n,p)$ when showing existence of solutions. 
For positive $(n,p)$, define new variables $(v,w)$ as follows:
\begin{equation}\label{quasi}
v=H(n)-u \quad\mbox{and}\quad w=H(p)+u
\end{equation}
where the enthalpy function $H$ is as defined in (\ref{enthalpy}).
Note that $H(\cdot)$ is strictly increasing from $(0,\infty)$ onto $(-h_0, h_\infty)$.
If we denote 
$$\Gamma=H^{-1}: \;\;\;s=\Gamma(\tau)\;\;\mbox{iff}\;\;\tau=H(s),$$ 
then $\Gamma(\cdot)$ is strictly increasing from $(-h_0, h_\infty)$ 
onto $(0,\infty)$ with $\Gamma(0)=1$. 
The relations (\ref{quasi}) can then be expressed as
\begin{equation}\label{def-np}
n=\Gamma(v+u)\quad\mbox{and}\quad p=\Gamma(w-u),
\end{equation}
and the current expressions in (\ref{current}) can be simplified as
$$J_n=\mu_n(\nabla P(n)-n\nabla u)=\mu_n\Gamma(v+u)\nabla v,$$
$$J_p=-\mu_p(\nabla P(p)+p\nabla u)=-\mu_p\Gamma(w-u)\nabla w.$$
The boundary conditions for $(v,w,u)$ are the same as in (\ref{bc})
with $\Gamma(\bar v+\bar u)=\bar n$ and $\Gamma(\bar w-\bar u)=\bar p$.
Hence we can solve for $(v,w,u)$ instead of $(n,p,u)$.
When showing the existence of non-vacuum solutions, we will establish positive lower
and upper bounds for $(n,p)$ or $(\Gamma(v+u),\Gamma(w-u))$ and
hence $v+u$ and $w-u$ will be always inside the domain of $\Gamma$. 
With such bounds established, the two sets of variables, $(v,w,u)$ and $(n,p,u)$,
become equivalent. In studying vacuum solutions, we will not use the $(v,w)$
variables since $(n,p)$ may become zero.

The rest of our presentation is organized as follows.
In section 2 we analyze the time-dependent model (\ref{w-n})--(\ref{w-u}),
and show the existence of global
positive $H^1\cap L^\infty$ solutions under Assumption (A). We carry out the proof by
first showing existence of the time-discretized systems, and then passing to 
the limit to obtain solutions to the continuous system. In this process, 
it is necessary for us to establish proper lower and upper a.e. bounds
for the solutions. When the generation-recombination term $G$ is given
by (\ref{G}), we need to assume $h_0$ sufficiently large to obtain
a positive lower bound. With an additional regularity assumption, 
we can also establish the uniqueness of the global solution. 
We use the $H^{-1}$ topology to deal with the nonlinear diffusion 
in $P(\cdot)$ for the uniqueness.
In section 3, we investigate the system as an infinite-dimensional dynamical
system acting on a set of admissible initial data $(n^0,p^0)$. 
We first sharpen the $L^\infty$ bounds for the solutions, so that an absorbing
set can be constructed. Again, for the lower bound, more restrictions on the 
model parameters are needed. We then establish the existence of a global 
attractor of the dynamical system. The semigroup defined by the system
is continuous with respect to the $H^{-1}$ topology, and we develop 
these results by using both the $H^{-1}$ and $L^\infty$ topologies.
In section 4 we prove the existence of solutions to the corresponding 
stationary system. All the discussions above are for the non-vacuum solutions,
and the extra conditions (e.g. $h_0$ is sufficiently large) 
we impose to obtain the positive lower bounds for $(n,p)$ are not needed for the upper bounds.
In section 5, we will discuss parallel results about vacuum solutions
by removing these conditions. In general vacuum solutions $(n,p)$ may have less regularity 
than the non-vacuum solutions, 
and we prove that in this case the solutions $(n, p, u)$ are in $L^\infty$ 
with $P(n)$, $P(p)$ and $u$ all in $H^1$. 

\section{Existence of Global Solutions}\nopagebreak
\setcounter{equation}{0} 
In this section we establish the existence of global solutions $(n,p,u)$ 
to the model equations (\ref{n})--(\ref{p}) with (\ref{u}) and boundary
conditions (\ref{bc}).

\subsection{The Time-Discretized System}\label{sec:dis}\nopagebreak
Given an integer $m>0$, let $h=T/m$. 
Consider the time-discretized version of the system (\ref{w-n})--(\ref{w-u}):
$(n^k, p^k, u^k)\in (\bar n,\bar p,\bar u)+V^3$ satisfying
\begin{equation}\label{nk}
\langle\frac{n^k- n^{k-1}}{h},\phi\rangle
+\mu_n\langle \nabla P(n^k)-n^k\nabla u^k,\nabla\phi\rangle 
=\langle G(n^k,p^k), \phi\rangle
\end{equation}
\begin{equation}\label{pk}
\langle\frac{p^k-p^{k-1}}{h},\psi\rangle
 +\mu_p\langle\nabla P(p^k)+p^k\nabla u^k,\nabla\psi\rangle 
 =\langle G(n^k,p^k),\psi\rangle
\end{equation}
\begin{equation}\label{uk}
\epsilon_0\langle\nabla u^k,\nabla\chi\rangle+\langle n^k-p^k-N, \chi\rangle=0
\end{equation}
for all $(\phi, \psi, \chi)\in V^3$. For the $(v,w)$ variables, we replace 
(\ref{nk})--(\ref{pk}) by
\begin{equation}\label{vk}
\langle\frac{n^k- n^{k-1}}{h},\phi\rangle
+\mu_n\langle n^k \nabla v^k,\nabla\phi\rangle = \langle G(n^k,p^k), \phi\rangle
\end{equation}
\begin{equation}\label{wk}
\langle\frac{p^k-p^{k-1}}{h},\psi\rangle
 +\mu_p\langle p^k \nabla w^k,\nabla\psi\rangle = \langle G(n^k,p^k), \psi\rangle
\end{equation}
where $(n^k, p^k)=(\Gamma(v^k+u^k), \Gamma(w^k-u^k))$, 
$k=1,2,\cdots,m$.  The two sets of equations are the same if one
can establish positive lower and upper bounds for $(n,p)=(\Gamma(v+u), \Gamma(w-u))$.

To show existence, we will consider the $(v,w)$-system.
We first make certain truncations of the function $\Gamma$
to ensure the properness of the nonlinearity in our setting, and then
``remove'' these truncations by establishing proper
$L^\infty$ bounds for the solutions and choosing the truncation constants
accordingly.

For $d_k\in (0, h_\infty)$ and $c_k\in (0, h_0)$ to be determined 
later (see (\ref{dk}), (\ref{ck}) below), define
\begin{equation}\label{Gammak}
\Gamma_k(s)=\left\{\begin{array}{ll}
     \Gamma(-c_k)\;\;\; &\mbox{when }\; s<-c_k\\[3pt]
     \Gamma(s) &\mbox{when }\;-c_k\le s\le d_k\\[3pt]
     \Gamma(d_k) &\mbox{when }\; s> d_k
  \end{array}\right.
\end{equation}
($k=1,2,\cdots, m$). Note that $\Gamma_k$'s are defined on 
$(-\infty,\infty)$. Then we consider (\ref{vk})--(\ref{wk}) 
and (\ref{uk}) with the replacement of 
\begin{equation}\label{npk}
n^k=\Gamma_k(v^k+u^k)\quad\mbox{and}\quad p^k=\Gamma_k(w^k-u^k),\;\;k=1,2,\cdots, m.
\end{equation}

Given a triplet $(v^{k-1}, w^{k-1}, u^{k-1})$, we will show the
existence of a solution to (\ref{vk})--(\ref{wk}) and (\ref{uk}) 
with (\ref{npk}).
To this end, we construct the solution map $\cal T$ from 
$(\tilde v,\tilde w,\tilde u)\in L^2(\Omega)^3$ to
the solution $(v, w, u)\in (\bar v,\bar w, \bar u)+V^3$ of the system
\begin{equation}\label{vk-T}
\langle\frac{\Gamma_k(v+u)- n^{k-1}}{h},\phi\rangle
+\mu_n\langle \tilde n \nabla v,\nabla\phi\rangle
=\langle G(\Gamma_k(v+\tilde u), \tilde p), \phi\rangle
\end{equation}
\begin{equation}\label{wk-T}
\langle\frac{\Gamma_k(w-u)-p^{k-1}}{h},\psi\rangle
 +\mu_p\langle \tilde p \nabla w,\nabla\psi\rangle 
 =\langle G(\tilde n, \Gamma_k(w-\tilde u) ),\psi\rangle
\end{equation}
\begin{equation}\label{uk-T}
\epsilon_0\langle\nabla u,\nabla\chi\rangle+\langle \Gamma_k(v+u)-\Gamma_k(w-u)-N,\chi\rangle=0
\end{equation}
for all $(\phi, \psi, \chi)\in V^3$, where
$n^{k-1}=\Gamma_{k-1}(v^{k-1}+u^{k-1})$,
 $p^{k-1}=\Gamma_{k-1}(w^{k-1}-u^{k-1})$,
$\tilde n=\Gamma_k(\tilde v+\tilde u)$, and
$\tilde p=\Gamma_k(\tilde w-\tilde u)$.

\paragraph{Well-definedness of $\cal T$.}  It is straightforward to verify that 
(\ref{vk-T})--(\ref{uk-T}) is the gradient system of the following 
convex lower-semicontinuous functional
$$\begin{array}{rcl}
I(v,w,u)&=& \ds\frac{\mu_n}{2}\langle\tilde n\nabla v,\nabla v\rangle
      +\frac{\mu_p}{2}\langle\tilde p\nabla w,\nabla w\rangle
      +\frac{\epsilon_0}{2h}\langle\nabla u,\nabla u\rangle\\[5pt]
& & -\ds\int_\Omega \left( \int_0^{v+\tilde u}G(\Gamma_k(s),\tilde p) ds
                   +\int_0^{w-\tilde u}G(\tilde n, \Gamma_k(s))ds \right) dx \\[5pt]
& & +\ds\frac{1}{h}\int_\Omega
 \left(\int_0^{v+u}+\int_0^{w-u}\right)\Gamma_k(s)ds dx\\[5pt]
& & -\ds\frac{1}{h}\left(
 \langle n^{k-1}, v+u\rangle +\langle p^{k-1}, w-u\rangle+\langle N,u\rangle\right).
\end{array}$$
Therefore the operator defined by (\ref{vk-T})--(\ref{uk-T}) is maximal
monotone, and hence standard result (see, e.g., \cite{Bz}) yields the 
existence of a solution $(v,w,u)$ to (\ref{vk-T})--(\ref{uk-T}).

\paragraph{Compactness of $\cal T$.} Choose $\phi=v-\bar v$, $\psi=w-\bar w$
and $\chi=\frac{1}{h}(u-\bar u)$ in (\ref{vk-T})--(\ref{uk-T}) and
sum up the three to obtain
\begin{eqnarray*}
\lefteqn{ \mu_n\langle\tilde n\nabla v,\nabla(v-\bar v)\rangle
+\mu_p\langle\tilde p\nabla w,\nabla(w-\bar w)\rangle
+\frac{\epsilon_0}{h}\langle\nabla u,\nabla(u-\bar u)\rangle }\\
&=& \langle G(\Gamma_k(v+\tilde u), \tilde p),v-\bar v\rangle
  +\langle G(\tilde n, \Gamma_k(w-\tilde u)),w-\bar w\rangle  \\
&& -\frac{1}{h}\langle \Gamma_k(v+u), v-\bar v+u-\bar u \rangle
  -\frac{1}{h}\langle \Gamma_k(w-u),  w-\bar w-u+\bar u \rangle \\
&&  +\frac{1}{h}\langle N, u-\bar u\rangle \,.
\end{eqnarray*}
Then using the boundedness and monotonicity of $\Gamma_k$ we can show
\begin{eqnarray*}
\lefteqn{ \frac{\mu_n}{2}\Gamma(-c_k)|\nabla (v-\bar v)|_2^2
+\frac{\mu_p}{2}\Gamma(-c_k)|\nabla (w-\bar w)|_2^2
+\frac{\epsilon_0}{2h}|\nabla (u-\bar u)|_2^2   }\\
&\le & \frac{\mu_n}{2}\Gamma(d_k)|\nabla \bar v|_2^2
+\frac{\mu_p}{2}\Gamma(d_k)|\nabla \bar w|_2^2
+\frac{\epsilon_0}{2h}|\nabla \bar u|_2^2 \\
&&+ \overline{G}|\Omega|^{1/2}(|v-\bar v|_2+|w-\bar w|_2) \\
&& + \frac{1}{h}|\Gamma_k(\bar v+\bar u)|_2 |v-\bar v|_2
 + \frac{1}{h}|\Gamma_k(\bar w-\bar u)|_2 |w-\bar w|_2 \\
&& +\frac{1}{h}(|\Gamma_k(\bar v+\bar u)|_2+|\Gamma_k(\bar w-\bar u)|_2
     +|N|_2)|u-\bar u|_2
\end{eqnarray*}
where $\overline{G}=\max_{0\le n,p\le \Gamma(d_k)} G(n,p)$.
Hence, by the Poincar\'{e} inequality, we obtain
\begin{equation}\label{Vbound}
|\nabla(v-\bar v)|_2+|\nabla (w-\bar w)|_2+|\nabla(u-\bar u)|_2\le\mbox{const.}
\end{equation}
where the constant depends on $h$, $c_k$ and $d_k$. This immediately
implies that the range of map $\cal T$ is a bounded set in $H^1(\Omega)^3$, 
and hence $\cal T$ is compact from $L^2(\Omega)^3$ to $L^2(\Omega)^3$
due to the compact embedding of $H^1(\Omega)$ into $L^2$.

\paragraph{Continuity of $\cal T$.} 
Suppose $({\tilde v}_j, {\tilde w}_j, {\tilde u}_j)
\to (\tilde v, \tilde w,\tilde u)$ in $L^2(\Omega)^3$ as $j\to\infty$, 
and denote their images under $\cal T$ by $(v_j,w_j,u_j)$ and $(v,w,u)$, 
respectively. Then from (\ref{vk-T})--(\ref{uk-T}), by using the boundedness
and monotonicity of $\Gamma_k$, we can obtain
\begin{eqnarray*}
\lefteqn{ \mu_n\Gamma(-c_k)|\nabla (v_j-v)|_2^2
 +\mu_p\Gamma(-c_k)|\nabla (w_j-w)|_2^2
 +\frac{\epsilon_0}{h}|\nabla (u_j-u)|_2^2 } \\
&\le & -\mu_n \langle (\tilde n_j-\tilde n)\nabla v,\nabla (v_j-v)\rangle
  -\mu_p \langle (\tilde p_j-\tilde p)\nabla w,\nabla (w_j-w)\rangle\\
&& +\langle G(\Gamma_k(v_j+\tilde u_j), \tilde p_j)
        -G(\Gamma_k(v+\tilde u),\tilde p), v_j-v\rangle\\
&& +\langle G(\tilde n_j, \Gamma_k(w_j-\tilde u_j))
        -G(\tilde n, \Gamma_k(w-\tilde u)), w_j-w\rangle\\
&\le & \mu_n |(\tilde n_j-\tilde n)\nabla v|_2 |\nabla (v_j-v)|_2
  +\mu_p |(\tilde p_j-\tilde p)\nabla w|_2 |\nabla (w_j-w)|_2\\
&& +|G(\Gamma_k(v_j+\tilde u_j), \tilde p_j)
-G(\Gamma_k(v+\tilde u),\tilde p)|_2 |v_j-v|_2\\
&& +|G(\tilde n_j, \Gamma_k(w_j-\tilde u_j))
-G(\tilde n, \Gamma_k(w-\tilde u))|_2 |w_j-w|_2.
\end{eqnarray*}
To show $\cal T$ is continuous, we argue that every subsequence of 
$(v_j,w_j,u_j)$ has a further subsequence converging to $(v,w,u)$. 
This can be shown from the above inequality, the Poincar\'{e} inequality,
and using the following facts:
(i) every $L^2$ convergent sequence has a pointwise convergent subsequence;
(ii) $G$ and $\Gamma_k$ are Lipschitz continuous;
(iii) the $V$-norm of the solution $(v,w,u)$ is bounded by (\ref{Vbound});
(iv) the dominated convergence theorem. (For details, see, e.g., \cite{J}.)

Therefore, by Schauder's fixed point theorem, $\cal T$ has at least one
fixed point $(v,w,u)$, which obviously is a solution to the system 
(\ref{vk})--(\ref{wk}) and (\ref{uk}) with the substitution of (\ref{npk}). 

By choosing proper positive constants $d_k$ and $c_k$, 
we will show in the next subsection that such a solution 
$(v^k,w^k,u^k)$ has the following $L^\infty$ bounds:
\begin{equation}\label{Linfty}
-c_k\le v^k(x)+u^k(x)\le d_k \quad\mbox{and}\quad -c_k\le w^k(x)-u^k(x)\le d_k.
\end{equation}
Thus the truncations in $\Gamma$ become unnecessary, and the solution 
$(v^k, w^k, u^k)$ becomes a solution to (\ref{vk})--(\ref{wk}), (\ref{uk}) 
with the original $\Gamma$.
Furthermore, we see that $n^k=\Gamma(v^k+u^k)$, 
$p^k=\Gamma(w^k-u^k)$ and $u^k$ form a solution to (\ref{nk})--(\ref{uk}).

\subsection{The $L^\infty$ Bounds for Solutions}\label{sec:Linfty}\nopagebreak

Now we establish independently the $L^\infty$ bounds (\ref{Linfty}) for 
the solution $(v,w,u)$ to the truncated version of (\ref{vk})--(\ref{wk}), 
(\ref{uk}). Hence when we choose the truncation constants $d_k$ and $c_k$
accordingly, we can replace $\Gamma_k$ by $\Gamma$ and $(v,w,u)$ becomes
a solution to (\ref{vk})--(\ref{wk}). 

{\bf The Upper Bounds.} Set 
\begin{equation}\label{DD}
D_0=\max\{\sup_\Omega n^0, \; \sup_\Omega p^0,\;
\sup_{\Sigma_D} \bar n,\; \sup_{\Sigma_D} \bar p,\; 1\}
\end{equation}
For $B>0$ to be determined below (in (\ref{B})), set 
$d_k\in (0, h_\infty)$ such that
$$\Gamma(d_k)=D_0(1+B h)^k$$
and let $\phi_d=(v+u-d_k)_+$.
We also require $h<1/B$. Clearly $\phi_d\in V$. Choose this $\phi_d$ 
as the test function in both (\ref{vk}) and (\ref{uk}) to obtain
\begin{eqnarray*}
\lefteqn{ \frac{1}{h}\langle\Gamma(d_k)-\Gamma_{k-1}(v^{k-1}+u^{k-1}),\phi_d\rangle
 +\mu_n\Gamma(d_k)\langle\nabla v,\nabla\phi_d\rangle }\\
 &=&\langle G(\Gamma(d_k), \Gamma_k(w-u) ),\phi_d\rangle , \hspace{6cm}
 \end{eqnarray*}
 and
$$\epsilon_0\langle\nabla u,\nabla\phi_d\rangle
 +\langle\Gamma(d_k)-\Gamma_k(w-u)-N,\phi_d\rangle=0\,.$$
Since $\phi_d\ge 0$ and $\nabla\phi_d=\nabla (v+u)$ when $v+u>d_k$, 
combining the two inequalitites, we have
\begin{eqnarray}\label{phid} 
\lefteqn{ \mu_n\Gamma(d_k)|\nabla\phi_d|_2^2 }\\
&\le&  \langle G(\Gamma(d_k), \Gamma_k(w-u)),\phi_d\rangle
 +\langle\frac{1}{h}(\Gamma(d_{k-1})-\Gamma(d_k)),\phi_d\rangle
 +\frac{\mu_n}{\epsilon_0}\Gamma(d_k)\langle N,\phi_d\rangle . \nonumber
\end{eqnarray}
Note that
$$\frac{1}{h}(\Gamma(d_{k-1})-\Gamma(d_k))=-\Gamma(d_k)\frac{B}{1+Bh}
\le-\Gamma(d_k)\frac{B}{2},$$
and for both models (\ref{G0}) and (\ref{G}) for $G$, we have
$$G(n,p)\le g(x)+Q(n,p)\le |g|_\infty+\overline Q (1+n+p).$$
Then (\ref{phid}) becomes
$$\mu_n|\nabla\phi_d|_2^2
\le\langle \frac{|g|_\infty+\overline{Q}(1+2D_0)}{D_0}
-\frac{B}{2}+\frac{\mu_n}{\epsilon_0}N,\phi_d\rangle\le 0$$
if we choose
\begin{equation}\label{B}
B=2(\frac{|g|_\infty+\overline{Q}(1+2D_0)}{D_0}
+\frac{\max(\mu_n,\mu_p)}{\epsilon_0}|N|_\infty).
\end{equation}
Therefore $\phi_d=0$ a.e. and $v+u\le d_k$ a.e. in $\Omega$.

The upper bounds above can be improved to be independent of $k$ and $h$
(i.e. time-independent) if the number $h_\infty$ is sufficiently large.
In the following we give the proof of this claim.
For $\tau$ such that $H(D_0)\le\tau<h_\infty$, let 
$d_k=\tau$ for all $k$. Set $\phi_\tau=(v+u-\tau)_+$ as the test function 
in both (\ref{vk}) and (\ref{uk}) and then combine the two to obtain 
$$\mu_n\Gamma(\tau)|\nabla\phi_\tau|_2^2
\le \langle G(\Gamma(\tau),\Gamma_k(w-u)),\phi_\tau\rangle
+\frac{\mu_n}{\epsilon_0}\Gamma(\tau)\langle N,\phi_\tau\rangle .$$
Hence
$$|\nabla\phi_\tau|_2^2\le 
\langle\frac{g+\overline{Q}(1+2D_0)}{D_0\min(\mu_n,\mu_p)}
    +\frac{N}{\epsilon_0},\phi_\tau\rangle$$
Then by applying Lemma \ref{lm:Linfty} in the Appendix, we obtain
$v(x)+u(x)\le d_0^\ast$ a.e. in $\Omega$ with
\begin{equation}\label{d0*}
d_0^\ast=H(D_0)+4\alpha_6^2|\Omega|^{1/6} 
 \left(\frac{|g|_2+\overline{Q}(1+2D_0)|\Omega|^{1/2}}{D_0\min(\mu_n,\mu_p)}
    +\frac{|N|_2}{\epsilon_0}\right)
\end{equation}
provided that $h_\infty> d_0^\ast$. Thus we obtain the time-independent bound 
$d_0^\ast$.

The same upper bound for $w(x)-u(x)$ can be established similarly.

In summary, we have established the upper bound $\Gamma(d_k)$ for both $n=\Gamma(v+u)$ and
$p=\Gamma(w-u)$ as
\begin{equation}\label{dk}
\Gamma(d_k)=\left\{\begin{array}{ll}
            D_0 (1+Bh)^k \;\;\;& \mbox{in general}\\
            \Gamma(d_0^\ast) & \mbox{if}\;\;h_\infty> d_0^\ast\end{array}\right.
\end{equation}
where $D_0$ is given in (\ref{DD}), $B$ in (\ref{B}) and $d_0^\ast$ in (\ref{d0*}).

{\bf The Lower Bounds.} Set 
\begin{equation}\label{CC}
C_0=\min\{\inf_\Omega n^0,\; \inf_\Omega p^0\,\;
\inf_{\Sigma_D} \bar n,\; \inf_{\Sigma_D} \bar p\}.
\end{equation}
For $0<\tilde{C}_0\le C_0$ to be chosen below, set $c_k\in (0,h_0)$ such that
$$\Gamma(-c_k)=\tilde{C}_0 (1+B h)^{-k} $$
and let $\phi_c=(v+u+c_k)_-$.
Then $\phi_c\le 0$ and $\phi_c\in V$. 
Choose this $\phi_c$ as the test function in both (\ref{vk}) and (\ref{uk}), 
and combine the two to obtain
$$\begin{array}{rl}
& \mu_n\Gamma(-c_k)|\nabla\phi_c|^2_2\\[5pt]
= &\frac{1}{h}\langle \Gamma(-c_k)-\Gamma_{k-1}(v^{k-1}+u^{k-1}),-\phi_c\rangle
   -\langle G(\Gamma(-c_k), \Gamma_k(w-u) ), -\phi_c\rangle\\[5pt]
&  +\frac{\mu_n}{\epsilon_0}\Gamma(-c_k)
   \langle \Gamma(-c_k)-\Gamma_k(w-u)-N,-\phi_c\rangle\\[5pt]
\le &\langle\frac{1}{h}(\Gamma(-c_k)-\Gamma(-c_{k-1}))
    -G(\Gamma(-c_k), \Gamma_k(w-u)) -\frac{\mu_n}{\epsilon_0}\Gamma(-c_k) N, -\phi_c\rangle\\[5pt]
\le & -\Gamma(-c_k)\langle B +\frac{\mu_n}{\epsilon_0} N, -\phi_c\rangle
  \;\;\;\mbox{(if $G(\Gamma(-c_k), \Gamma(w-u))\ge 0$)}\\[5pt]
\le & 0 \;\;\;\mbox{(since $B$ is from (\ref{B}))}
\end{array}$$
Then we have $\phi_c\equiv 0$, i.e., $v(x)+u(x)\ge -c_k$ a.e. in $\Omega$.

Now we verify that $G(\Gamma(-c_k), \Gamma(w-u))\ge 0$ for the two models (\ref{G0}) and (\ref{G})
separately, by choosing different constants $\tilde{C}_0$. 
For $G$ given by (\ref{G0}), choose $\tilde{C}_0=\min(C_0, 1/D_0)$. Then
$$G(\Gamma(-c_k), \Gamma(w-u))\ge g-Q(\Gamma(-c_k)\Gamma(d_k)-1)\ge Q(1-\tilde{C}_0 D_0)\ge 0.$$
For $G$ given by (\ref{G}), we need to assume
\begin{equation}\label{new}
\min(d_0^\ast, h_\infty)<h_i+ h_0.
\end{equation}
From (\ref{dk}) we see $H(\Gamma_k(w-u))\le d_k$ is no larger than either 
$d_0^\ast$ or $h_\infty$.  
By (\ref{new}), $h_i-\min(d_0^\ast,h_\infty)>-h_0$, hence we can 
choose $\tilde{C}_0=\min(C_0, \Gamma(h_i-\min(d_0^\ast,h_\infty))$. Then
$$H(\Gamma(-c_k))+H(\Gamma(d_k))\le H(\tilde{C}_0)+\min(d_0^\ast,h_\infty)\le h_i$$
and hence
$$G(\Gamma(-c_k), \Gamma(w-u))\ge g-Q( S(H(\Gamma(-c_k))+H(\Gamma(d_k))-1)\ge 0$$
since $S$ is nondecreasing and $S(h_i)=1$.

Hence we have the positive lower bound $\Gamma(-c_k)=\tilde{C}_0(1+Bh)^{-k}$ 
for $n$ and $p$, where
\begin{equation}\label{tC0}
\tilde{C}_0=\left\{\begin{array}{ll}
     \min\{C_0, 1/D_0\} &
          \mbox{for $G$ in (\ref{G0})}\\
     \min\{C_0, \Gamma(h_i-\min(d_0^\ast, h_\infty) ) \}\;\;&
          \mbox{for $G$ in (\ref{G}) with } \\
          & \min(d_0^\ast, h_\infty)<h_i+h_0
\end{array}\right.
\end{equation}

In the case $h_\infty>d_0^\ast$, $\Gamma(d_0^\ast)$ is a time-independent
upper bound for both $n$ and $p$, then we can also establish a time-independent
positive lower bound for $n$ and $p$ provided that $h_0>c_0^\ast$ (to be 
given below). To prove this claim, 
let $c_k=\tau$ for all $k$ with $-H(C_0)\le\tau<h_0$. 
Set $\phi_\tau=(v+u+\tau)_-=-(-(v+u)-\tau)_+\in V$. Following the same
calculations as above leads to 
$$\mu_n\Gamma(-\tau)|\nabla\phi_\tau|_2^2
\le \langle G(\Gamma(-\tau), \Gamma_k(w-u)), \phi_\tau\rangle
 +\frac{\mu_n}{\epsilon_0}\Gamma(-\tau)\langle N, \phi_\tau\rangle .$$
As in choosing $\tilde{C}_0$ above, we can see that
$G(\Gamma(-\tau), \Gamma_k(w-u))\ge 0$
if we choose $\tau$ such that $\Gamma(-\tau)\le 1/\Gamma(d_0^\ast)$
when $G$ is given by (\ref{G0}), and choose $\tau\ge d_0^\ast-h_i$ when 
$G$ is by (\ref{G}). Thus we have
$$|\nabla\phi_\tau|_2^2\le\langle\frac{N}{\epsilon_0},\phi_\tau\rangle .$$
From Lemma \ref{lm:Linfty} we conclude that 
$v+u\ge -c_0^\ast$ a.e. in $\Omega$ for
\begin{equation}\label{c0*}
c_0^\ast=\tilde{c}+4\frac{\alpha_6^2}{\epsilon_0}|\Omega|^{1/6}|N|_2
\end{equation}
where
\begin{equation}\label{tc0}
\tilde{c}=\left\{\begin{array}{ll}
       \max\{-H(C_0), -H(1/\Gamma(d_0^\ast))\}\;\;&\mbox{for model (\ref{G0})}\\
       \max\{-H(C_0), d_0^\ast-h_i\} & \mbox{for model (\ref{G})}\end{array}\right.
\end{equation}
provided that $h_0>c_0^\ast$. 

Similarly we can show the same lower bounds for $w(x)-u(x)$ in all situations.

In summary, we have established the positive lower bound $\Gamma(-c_k)$ for both 
$n=\Gamma(v+u)$ and $p=\Gamma(w-u)$ as
\begin{equation}\label{ck}
\Gamma(-c_k)=\left\{\begin{array}{ll} 
  \tilde{C}_0(1+Bh)^{-k}\;\;\; & \mbox{in general}\\
  \Gamma(-c_0^\ast)   &\mbox{if $h_\infty>d_0^\ast$ and $h_0>c_0^\ast$}
\end{array}\right.
\end{equation}
where $\tilde{C}_0$ is from (\ref{tC0}) and $c_0^\ast$ from (\ref{c0*}).
    
{\bf Bounds for $u$.} To establish the $L^\infty$ bounds for $u$ in (\ref{uk}), 
we use the same $L^\infty$ estimate technique as in showing the time-independent
upper bound for $v+u$ above.
Let $\bar\tau=\sup_{\Sigma_D}|\bar u|$.
For $\tau\ge \bar\tau$, choose $u_\tau=(u-\tau)^+$ as the test function 
in (\ref{uk}) to obtain
$$|\nabla u_\tau|_2^2\le\langle \frac{1}{\epsilon_0}(n-p-N),u_\tau\rangle .$$
By Lemma \ref{lm:Linfty} we obtain
$u\le \bar\tau+4\frac{\alpha_6^2}{\epsilon_0}|\Omega|^{1/6} |n-p-N|_2$.
Applying the same argument for $-u$ we can obtain the lower bound estimate. 
Hence we have the following $L^\infty$ estimates for $u$ in (\ref{uk}):
\begin{equation}\label{Linfty-u}
|u(t,x)|\le \sup_{\Sigma_D} |\bar u(x)|
 +4\frac{\alpha_6^2}{\epsilon_0} |\Omega|^{1/6} 
  |n(t,\cdot)-p(t,\cdot)-N(\cdot)|_2\;\;\;
\mbox{a.e. in}\;\;(0,T)\times\Omega.
\end{equation}

\subsection{Energy Estimates}\nopagebreak
Now we establish some energy estimates for solutions to
the $(n,p)$-system (\ref{nk})--(\ref{uk}), which will be needed 
in the next subsection for convergence argument.

First we note that 
\begin{equation}\label{npk-Linfty}
C_T\equiv \tilde{C}_0 e^{-BT}\le \Gamma(-c_k)\le
n^k,\; p^k\le \Gamma(d_k)\le  D_0 e^{BT}\equiv D_T 
\end{equation}
for all $k=1, 2, \cdots, m$. Set
$$\kappa_1=\min P'(s)\quad\mbox{and}\quad \kappa_2=\max P'(s)
\;\;\;\mbox{over}\;\;s\in [C_T, D_T].$$

Choose $\chi=u^k-\bar u\in V$ in (\ref{uk}) to obtain
$$\begin{array}{rl}
& \frac{\epsilon_0}{2}(|\nabla u^k|_2^2-|\nabla\bar u|_2^2
      +|\nabla(u^k-\bar u)|_2^2) \\[5pt]
\le & |n^k-p^k-N|_2 |u^k-\bar u|_2
\le \alpha_2|n^k-p^k-N|_2 |\nabla(u^k-\bar u)|_2,
\end{array}$$
and hence
$$|\nabla u^k|_2^2
\le |\nabla\bar u|_2^2+\frac{\alpha_2^2}{\epsilon_0^2}|n^k-p^k-N|_2^2
\le |\nabla\bar u|_2^2+\frac{\alpha_2^2}{\epsilon_0^2}
    (2 D_T |\Omega|^{1/2}+|N|_2)^2\equiv M_1. $$
Then choose $\phi=n^k-\bar n\in V$ in (\ref{nk}) to yield
$$\begin{array}{rl}
& \frac{1}{2h}(|n^k-\bar n|^2_2-|n^{k-1}-\bar n|_2^2+|n^k-n^{k-1}|_2^2)
    +\mu_n\kappa_1 |\nabla(n^k-\bar n)|_2^2\\[5pt]
\le & \kappa_2|\nabla\bar n|_2 |\nabla (n^k-\bar n)|_2
      + D_T|\nabla u^k|_2 |\nabla (n^k-\bar n)|_2
      +\overline{G}_T|\Omega|^{1/2}|n^k-\bar n|_2,
\end{array}$$
where $\overline{G}_T=\max G(n,p)$ over $n, p\in [C_T, D_T]$,
which leads to
$$\begin{array}{rl}
& \frac{1}{2h}(|n^k-\bar n|^2_2-|n^{k-1}-\bar n|_2^2+|n^k-n^{k-1}|_2^2)
    +\frac{\mu_n\kappa_1}{4} |\nabla(n^k-\bar n)|_2^2\\[5pt]
\le & \frac{1}{\mu_n\kappa_1}\left\{\kappa_2^2|\nabla\bar n|_2^2 
      + D_T^2 M_1 
      +\alpha_2^2\overline{G}_T^2|\Omega|\right\}
 \equiv \frac{1}{\mu_n} M_2.
\end{array}$$
Sum up the above for $1\le k\le m$ to obtain
\begin{equation}\label{energy-n}
 \sum_{k=1}^m |n^k-n^{k-1}|_2^2
 + h\frac{\mu_n\kappa_1}{2}\sum_{k=1}^m |\nabla(n^k-\bar n)|_2^2
\le |n^0-\bar n|_2^2 + \frac{2}{\mu_n} M_2 T.
\end{equation}
Similarly we can obtain for $p^k$ the estimate
\begin{equation}\label{energy-p}
 \sum_{k=1}^m |p^k-p^{k-1}|_2^2
 + h\frac{\mu_p\kappa_1}{2}\sum_{k=1}^m |\nabla(p^k-\bar p)|_2^2
\le |p^0-\bar p|_2^2 + \frac{2}{\mu_p} M_2 T.
\end{equation}

\subsection{Convergence to the Continuous System}\nopagebreak
In this subsection we use the solutions $\{(n^k, p^k, u^k)\}_{k=0}^m$
to construct solutions of continuous time to (\ref{w-n})--(\ref{w-u}).

Let $t_k= k h$, and define the following interpolation (in time) functions
$n^{(1)}_t(t,x)$ and $n^{(2)}_h(t,x)$ as
$$n^{(1)}_h(t,x)=n^k(x)+\frac{t-t_k}{h}(n^k(x)-n^{k-1}(x))$$
and 
$n^{(2)}_h(t,x)=n^k(x)$ on $(t_{k-1},t_k]$,  
($k=1,2,\cdots,m$). From this definition, we see that
$$ \int_0^T |n^{(1)}_h(t)- n^{(2)}_h(t)|_2^2 dt 
 \le \frac{h}{3} \sum_{k=1}^{m}|n^k-n^{k-1}|_2^2,$$
$$\int_0^T |\nabla (n^{(1)}_h(t)-\bar n)|_2^2 dt 
\le h \sum_{k=0}^{m}|\nabla (n^k-\bar n)|_2^2,$$
and 
$$\int_0^T |\nabla (n^{(2)}_h(t)-\bar n)|_2^2 dt 
\le h \sum_{k=1}^{m}|\nabla (n^k-\bar n)|_2^2. $$
Similarly we define such functions for the other two variables $p$ and $u$,
and the same estimates hold. 
From (\ref{energy-n}) and (\ref{energy-p}) we then have
\begin{equation}\label{n-diff}
\int_0^T |n^{(1)}_h(t)- n^{(2)}_h(t)|_2^2 dt 
 \le \mbox{const.} h,
\end{equation}
\begin{equation}\label{n-V}
\int_0^T |\nabla (n^{(1)}_h(t)-\bar n)|_2^2 dt,\;\;\;
\int_0^T |\nabla (n^{(2)}_h(t)-\bar n)|_2^2 dt \le\mbox{const.}
\end{equation}
and the same hold for $p$.

Using these interpolation functions, we can write the equations 
(\ref{nk})--(\ref{uk}) as follows
\begin{eqnarray}\label{w-n-i}
\lefteqn{ \langle  n^{(1)}_h(t)- n^0,\phi\rangle }\\
&&+\int_0^t\!\! \left\{\mu_n\langle P'(n^{(2)}_h(s)) \nabla n^{(2)}_h(s)
           -n^{(2)}_h(s)\nabla u^{(2)}_h(s),\nabla\phi\rangle
- \langle G_h, \phi\rangle\right\} ds =0 \nonumber
\end{eqnarray}
\begin{eqnarray}\label{w-p-i}
\lefteqn{ \langle  p^{(1)}_h(t)- p^0,\psi\rangle }\\
&&+\int_0^t\!\! \left\{\mu_p\langle P'(p^{(2)}_h(s)) \nabla p^{(2)}_h(s)
           +p^{(2)}_h(s)\nabla u^{(2)}_h(s),\nabla\psi\rangle
 -\langle G_h , \psi\rangle\right\} ds =0 \nonumber 
\end{eqnarray}
\begin{equation}\label{w-u-i}
\epsilon_0\langle\nabla u^{(2)}_h(t),\nabla\chi\rangle
+\langle n^{(2)}_h(t)-p^{(2)}_h(t)-N, \chi\rangle=0
\end{equation}
for all $(\phi, \psi, \chi)\in V^3$, $t\in [0,T]$. Here we use
$G_h=G(n^{(2)}_h(s), p^{(2)}_h(s))$. By the $L^\infty$ bounds 
(\ref{npk-Linfty}), we can see easily from (\ref{w-n-i})--(\ref{w-p-i})
that 
\begin{equation}\label{V*-bound}
\int_0^T\left|\frac{\partial n^{(1)}_h(t)}{\partial t}\right|_{V^\ast}^2 dt,\;\;\;
\int_0^T\left|\frac{\partial p^{(1)}_h(t)}{\partial t}\right|_{V^\ast}^2 dt
\le \mbox{const.}
\end{equation} 
Therefore (\ref{n-V}) and (\ref{V*-bound}) together yield that
\begin{equation}\label{Aubin}
\{n^{(1)}_h(t)-\bar n\}_{h>0},\;\;\;
\{p^{(1)}_h(t)-\bar p\}_{h>0}\;\;\;
\mbox{are bounded in $L^2(0,T; V)\cap H^1(0,T;V^\ast)$} 
\end{equation}
(uniformly in $h>0$), and hence Aubin's Lemma implies that they are precompact
in $L^2(0,T;L^2(\Omega))$. Therefore we can conclude that there exist 
a sequence of positive numbers $\{h\}$ and a pair of functions $(n(t),p(t))$ 
in $L^2(0,T;L^2(\Omega))$  such that, as $h\to 0^+$,
$$(n^{(1)}_h(t), p^{(1)}_h(t))\to (n(t), p(t))\;\;\;
\mbox{strongly in $L^2(0,T; L^2(\Omega))$} $$
and, from (\ref{n-V}), by taking further subsequences if necessary, 
$$(n^{(1)}_h(t), p^{(1)}_h(t))\to (n(t), p(t))\;\;\;
\mbox{weakly in $L^2(0,T;H^1(\Omega))$}.$$ 
Hence $(n(t)-\bar n, p(t)-\bar p)$ are in $L^2(0,T; V)$. In view of (\ref{n-diff}), the above
convergence results also hold for $(n^{(2)}_h(t), p^{(2)}_h(t) )$. Since
every $L^2$-convergent sequence has an a.e. convergent subsequence, we
can further assume that
$$(n^{(2)}_h(t,x), p^{(2)}_h(t,x) )\to (n(t,x), p(t,x))\;\;\;
\mbox{a.e. in}\;\;\;(0,T)\times\Omega.$$
Hence we see from (\ref{npk-Linfty}) that $(n,p)$ are in 
$L^\infty((0,T)\times \Omega)$, and the expressions for $d_k$ 
and $c_k$ in (\ref{dk}) and (\ref{ck}) clearly lead to
\begin{equation}\label{npt}
\tilde{C}_0 e^{-Bt}\le n(t,x),\;\; p(t,x)\le D_0 e^{Bt} 
\;\;\;\mbox{a.e. in }\;\;(0,T)\times\Omega
\end{equation}
(or with the upper bound replaced by $\Gamma(d^\ast_0)$ when $h_\infty>h^\ast$, and the lower
bound by $\Gamma(-c^\ast_0)$ if additionally $h_0>c^\ast_0$.)
For the $u$-component, it is easy to see from (\ref{w-u-i}) that there exists
$u(t)\in L^2(0,T;H^1(\Omega))$ satisfying
$$\epsilon_0\langle\nabla u(t),\nabla\chi\rangle
+\langle n(t)-p(t)-N, \chi\rangle=0$$
for all $\chi\in V$ such that
$$u(t)-\bar u\in V \quad\mbox{and}\quad  u^{(2)}_h(t)\to u(t) \;\;\;
\mbox{strongly in $H^1(\Omega)$, a.e. in $(0,T)$.} $$
Finally, we can pass to the limit $h\to 0^+$ in (\ref{w-n-i})--(\ref{w-p-i}),
by the dominated convergence theorem, to obtain
$$\displaylines{
 \langle  n(t)- n^0,\phi\rangle 
+\int_0^t \left\{\mu_n\langle \nabla P(n(s))-n(s)\nabla u(s),\nabla\phi\rangle
  -\langle G(n(s), p(s)), \phi\rangle\right\} ds =0 \cr
\langle  p(t)- p^0,\psi\rangle 
+\int_0^t \left\{\mu_p\langle \nabla P(p(s))+p(s)\nabla u(s),\nabla\psi\rangle
 -\langle G(n(s), p(s)), \psi\rangle\right\} ds =0
}$$
a.e. in $(0,T)$, which imply that $n(t)$ and $p(t)$ are in $H^1(0,T;V^\ast)$,
and the triplet $(n,p,u)$ is a solution to (\ref{w-n})--(\ref{w-u}).

\vspace{.1in}
 
We summarize our discussions of \S2.1-\S2.4 in the following theorem
about the existence of global solutions.
Recall the various constants $D_0$ defined in (\ref{DD}), 
$\tilde{C}_0$ in (\ref{ck}), $B$ in (\ref{B}), $d_0^\ast$ in (\ref{d0*}),
and $c_0^\ast$ in (\ref{c0*}). They are all independent of $T>0$.

\begin{theorem}\label{thm:exist} {\rm (Existence of global non-vacuum solutions.)}
Suppose Assumption (A) holds, and the generation-recombination rate $G$ is given either by
(\ref{G0}) or by (\ref{G}) under the condition $\min(d_0^\ast, h_\infty)<h_i+h_0$.
Then for any given $T>0$, there exists at least one 
solution $(n, p, u)$ to (\ref{w-n})--(\ref{w-u}) that has the following estimates: 
$$\displaylines{ 
\tilde{C}_0 e^{-Bt}\le n(t,x),\;\;p(t,x)\le  D_0 e^{Bt},\cr
 |u(t,x)|\le \sup_{\Sigma_D} |\bar u| 
+4\frac{\alpha_6^2}{\epsilon_0}|\Omega|^{1/6}
 |n(t,\cdot)-p(t,\cdot)-N(\cdot)|_2 
 } $$
a.e. in $(0,T)\times\Omega$. In the case $h_\infty>d_0^\ast$, the upper bound for $(n,p)$
can be replaced by the constant $\Gamma(d_0^\ast)$. If additionally $h_0>c_0^\ast$,
the lower bound for $(n,p)$ can be also replaced by the constant $\Gamma(-c_0^\ast)$.
\end{theorem}

\begin{remark}\label{rmk:exist}
The condition $\min(d_0^\ast, h_\infty)<h_i+h_0$ is used only to obtain
the lower bound constant $\tilde{C}_0$ when $G$ is given by (\ref{G}),
and it is true if $h_0=\infty$. 
\end{remark}

\subsection{Uniqueness}\label{sec:unique}\nopagebreak
In this subsection we discuss the uniqueness of the global solutions
for the system (\ref{w-n})--(\ref{w-p}) with (\ref{w-u}).
It turns out that we need to impose a
regularity assumption on the Laplacian operator with the 
mixed boundary conditions as in (\ref{bc}), and to use the weaker $H^{-1}$ 
topology to deal with the nonlinear $P(\cdot)$.
These extra cares are not necessary in the case of linear diffusion ($P(\cdot)$ linear),
as shown in \cite{JDE2,JDE3}. 

We make the following additional regularity assumption for the uniqueness:


\paragraph{Assumption (A1):} {\em  For $f\in L^2(\Omega)$, the weak solution 
 $z\in V$ to
$$\langle\nabla z,\nabla\chi\rangle=\langle f,\chi\rangle\;\;\;\mbox{for all}\;\;\chi\in V$$
is in $H^2(\Omega)$, and there exists a constant (depending on $\Omega$ and 
$\Sigma_D$, $\Sigma_N$) such that
$$|z|_{H^2(\Omega)}\le \mbox{\rm const.} |f|_2.$$  }


Such a regularity assumption is in fact an assumption on the domain $\Omega$
and the partition of its boundary $\partial\Omega=\Sigma_D\cup\Sigma_N$.
In the case of one-dimensional domain, it is always satisfied.
When $\Omega$ is a rectangle in ${mathbb R}^2$ and $\Sigma_D$ consists of the two opposite
sides of the rectangle, we know that $u$ is in $H^2$. In general, such regularity
for the solution is true if $\Omega$ is of the class $C^{1,1}$ and
$\Sigma_N$ is both open and closed (see \cite[Theorem 2.24]{Tr}), which includes
cases such as $\Sigma_N=\emptyset$ (the pure Dirichlet problem) or when $\Sigma_N$ and
$\Sigma_D$ are disjoint portions of $\partial\Omega$ (annulus domains).

We also need to use the weaker $H^{-1}$ topology to 
deal with the nonlinear $P(\cdot)$. Let $\Delta_0$ be the
Laplacian operator defined on $H^2(\Omega)\cap V$. Then $-\Delta_0$ is
a positive self-adjoint operator on $L^2(\Omega)$, 
and $-\Delta_0\in{\cal L}(V,V^\ast)$ with
$$\langle-\Delta_0\phi,\psi\rangle_{V^\ast\times V}=\langle\nabla\phi,\nabla\psi\rangle
\;\;\;\mbox{for }\;\;\phi,\;\psi\in V.$$
Let $(-\Delta_0)^{-1}\in {\cal L}(V^\ast, V)$ be the Riesz map of $V^\ast$
equipped with the norm
$$ |\phi|^2_{V^\ast}=\langle\phi, (-\Delta_0)^{-1}\phi\rangle_{V^\ast\times V}
\;\;\;\mbox{for }\;\;\phi\in V^\ast.$$
That is, $|\phi|_{V^\ast}=|(-\Delta_0)^{-1}\phi|_{V}$ for $\phi\in V^\ast$.
We also see that
$$\langle\nabla\phi,\nabla(-\Delta_0)^{-1}\xi\rangle=\langle\phi,\xi\rangle\;\;\;
\mbox{for }\;\phi\in V,\;\;\xi\in L^2(\Omega),$$
and, by Assumption (A1), 
$$|(-\Delta_0)^{-1}\xi|_{H^2(\Omega)}\le\mbox{const.}|\xi|_2
\;\;\;\mbox{for}\;\;\xi\in L^2(\Omega).$$

\begin{theorem}\label{thm:unique} {\rm (Uniqueness.)}
Under Assumptions (A) and (A1), the $H^1\cap L^\infty$ weak solution
$(n,p,u)$ to (\ref{w-n})--(\ref{w-u}) is unique.
\end{theorem}

{\em Proof.} Let $(n,p,u)$ and $(\hat n,\hat p,\hat u)$ be $H^1\cap L^\infty$
weak solutions to (\ref{w-n})--(\ref{w-u}) corresponding to initial data 
$(n^0,p^0)$ and $({\hat n}^0, {\hat p}^0)$.
From (\ref{w-n}) we have the equation for $n-\hat n$:
\begin{eqnarray}
\lefteqn{ \langle \ds\frac{\partial (n-\hat n)}{\partial t},\phi\rangle_{V^\ast\times V}
 +\mu_n\langle\nabla(P(n)-P(\hat n)),\nabla\phi\rangle
 -\mu_n\langle n\nabla u-\hat n\nabla\hat u,\nabla\phi\rangle  }\nonumber\\
 &=& \langle G(n, p)- G(\hat n, \hat p),\phi\rangle  \hspace{75mm}
\label{n-hat}
\end{eqnarray}
for all $\phi\in V$ and a similar equation for $p-\hat p$ from (\ref{w-p}). 
For $u-\hat u$ there holds
$$\epsilon_0\langle\nabla (u-\hat u),\nabla\chi\rangle+\langle n-\hat n- (p-\hat p),\chi\rangle=0
$$
for all $\chi\in V$, and hence
\begin{equation}\label{u-hat}
|\nabla(u-\hat u)|_2\le \frac{\alpha_2}{\epsilon_0}(|n-\hat n|_2+|p-\hat p|_2).
\end{equation}
Now choose $\phi=(-\Delta_0)^{-1}(n-\hat n)\in V$ in
(\ref{n-hat}), and we estimate the terms as follows.
The first and the second terms in (\ref{n-hat}) become:
$$\langle\frac{\partial (n-\hat n)}{\partial t},(-\Delta_0)^{-1}(n-\hat n)\rangle_{V^\ast\times V}
=\frac{d}{dt}\left(\frac{1}{2}|n-\hat n|^2_{V^\ast}\right)$$
and
$$\langle\nabla (P(n)-P(\hat n)),\nabla(-\Delta_0)^{-1}(n-\hat n)\rangle
 =\langle P(n)-P(\hat n), n-\hat n\rangle\ge \kappa_1|n-\hat n|^2_2.$$
The third term is
\begin{equation}\label{reg}
\begin{array}{rl}
& |\langle n\nabla u-\hat n\nabla\hat u, \nabla(-\Delta_0)^{-1}(n-\hat n)\rangle| \\[5pt]
\le & |n-\hat n|_2 |\nabla u|_4 |\nabla(-\Delta_0)^{-1}(n-\hat n)|_4+
      |\hat n|_\infty |\nabla(u-\hat u)|_2 |\nabla(-\Delta_0)^{-1}(n-\hat n)|_2 \\[5pt]
\le & |n-\hat n|_2 |\nabla u|_4 \mbox{const.} |n-\hat n|_2^{3/4}|n-\hat n|_{V^\ast}^{1/4}
                      \;\;\;\;\;\;\mbox{(see (\ref{3/4}) below)}\\[5pt]
  &   + |\hat n|_\infty \frac{\alpha_2}{\epsilon_0}(|n-\hat n|_2+|p-\hat p|_2) 
       |n-\hat n|_{V^\ast}\;\;\;\;\;\;\mbox{(from (\ref{u-hat}))}\\[5pt]
\le & \mbox{const.}\left(|n-\hat n|_2^{7/4} |n-\hat n|_{V^\ast}^{1/4}+
  (|n-\hat n|_2+|p-\hat p|_2) |n-\hat n|_{V^\ast}\right)
\end{array}
\end{equation}
where the constant depends on $|\nabla u|_4$ and $|\hat n|_\infty$, 
and we have used the following inequality:
\begin{equation}\label{3/4}
\begin{array}{rl}
|\nabla(-\Delta_0)^{-1}(n-\hat n)|_4 
& \le |\nabla(-\Delta_0)^{-1}(n-\hat n)|_6^{3/4} |\nabla(-\Delta_0)^{-1}(n-\hat n)|_2^{1/4}\\[5pt]
& \le \mbox{const.} |\nabla(-\Delta_0)^{-1}(n-\hat n)|_{V}^{3/4} |n-\hat n|_{V^\ast}^{1/4}\\[5pt]
 &\hspace{1in}\mbox{(by the embedding of $V$ into $L^6$)}\\[5pt]
& \le \mbox{const.}|(-\Delta_0)^{-1}(n-\hat n)|_{H^2}^{3/4} |n-\hat n|_{V^\ast}^{1/4}\\[5pt]
& \le  \mbox{const.} |n-\hat n|_2^{3/4} |n-\hat n|_{V^\ast}^{1/4}.
\end{array}
\end{equation}
The right side of (\ref{n-hat}) becomes
\begin{eqnarray*}
\lefteqn{ |\langle G(n, p)-G(\hat n, \hat p), (-\Delta_0)^{-1}(n-\hat n)
\rangle|   }\\
& \le& \mbox{const.} (|n-\hat n|_2+|p-\hat p|_2)|(-\Delta_0)^{-1}(n-\hat n)|_2\\
& =& \mbox{const.} (|n-\hat n|_2+|p-\hat p|_2)|n-\hat n|_{V^\ast}
\end{eqnarray*}
where the constant depends on the Lipschitz continuity of $G$
and $|(n,p)|_\infty$, $|(\hat n,\hat p)|_\infty$. 

Hence we obtain
\begin{eqnarray*}
\lefteqn{ \ds\frac{d}{dt}\left(\frac{1}{2}|n-\hat n|^2_{V^\ast}\right)
+\mu_n\kappa_1 |n-\hat n|^2_2  }\\
&\le & \mbox{const.}\left(|n-\hat n|_2^{7/4} |n-\hat n|_{V^\ast}^{1/4}
 +(|n-\hat n|_2+|p-\hat p|_2) |n-\hat n|_{V^\ast}\right).
 \end{eqnarray*}
A similar inequality can be obtained for $p-\hat p$. 
Combining the two yields
\begin{eqnarray*}
\lefteqn{ \frac{d}{dt}(|n-\hat n|^2_{V^\ast}+|p-\hat p|^2_{V^\ast})
+2\kappa_1\min\{\mu_n,\mu_p\}(|n-\hat n|_2^2+|p-\hat p|_2^2) }\\
&\le &  \mbox{const.}\left(|n-\hat n|_2^{7/4} |n-\hat n|_{V^\ast}^{1/4}+
 |p-\hat p|_2^{7/4} |p-\hat p|_{V^\ast}^{1/4} \right.\\
&&\hspace{15mm}\left. +(|n-\hat n|_2+|p-\hat p|_2)(|n-\hat n|_{V^\ast}+|p-\hat p|_{V^\ast} )\right)
\end{eqnarray*} 
and hence, by Young's inequality ($ab\le a^p/p+b^q/q$),
$$\frac{d}{dt}(|n-\hat n|^2_{V^\ast}+|p-\hat p|^2_{V^\ast})
\le K (|n-\hat n|^2_{V^\ast}+|p-\hat p|^2_{V^\ast})$$
where the constant $K$ depends only on the $L^\infty$ bounds of the solutions
and $|\nabla u|_4$. Thus Gronwall's inequality leads to
\begin{equation}\label{gron}
|n-\hat n|^2_{V^\ast}+|p-\hat p|^2_{V^\ast}
\le e^{Kt} (|n^0-{\hat n}^0|^2_{V^\ast}+|p^0-{\hat p}^0|^2_{V^\ast}).
\end{equation}
Hence (\ref{gron}) and (\ref{u-hat}) imply the uniqueness of $H^1\cap L^\infty$ solutions. $\Box$

\begin{remark}
The $H^2$-regularity of the Laplacian operator in Assumption (A1) can be replaced by 
$W^{1,\infty}$-regularity. 
\end{remark}
Indeed, in the proof above, the only place we use (A1) is in the third term estimate (\ref{reg}), 
and when $u$ is in $W^{1,\infty}$, this term can be easily estimated as
\begin{eqnarray*}
\lefteqn{ |\langle n\nabla u-\hat n\nabla\hat u, \nabla(-\Delta_0)^{-1}
(n-\hat n)\rangle|  }\\
&\le& \left(|n-\hat n|_2 |\nabla u|_\infty+|\hat n|_\infty|\nabla(u-\hat u)|_2\right)
   |n-\hat n|_{V^\ast}
\end{eqnarray*}
and then proceed as above to obtain (\ref{gron}) with $K$ depending on $|\nabla u|_\infty$.

%%%%%%%%%%%%%%%%%%%%

\section{Long-Time Behavior of Solutions}\nopagebreak
\setcounter{equation}{0}

In this section we study the system (\ref{w-n})--(\ref{w-p}) for 
$(n,p)$ as an infinite-dimensional dynamical system acting on an 
admissible set of initial data $(n^0,p^0)$ for its long-time behavior. 
In particular, we will establish sharper $L^\infty$
bounds for any $L^\infty\cap H^1$ global solutions so that an absorbing 
set can be constructed, and then show the existence of a global compact 
attractor of the dynamical system.

\subsection{Construction of an Absorbing Set}\label{sec:absorb}\nopagebreak
The existence of an absorbing set amounts to
careful analysis of the dependence of bounds for solutions on the initial
conditions. To this end, we will work on the continuous system 
(\ref{w-n})--(\ref{w-p}) for sharper $L^\infty$ bounds of the solutions.
Suppose $(n,p)$ is a $H^1\cap L^\infty$ global solution to 
(\ref{w-n})--(\ref{w-p}) with initial data $(n^0,p^0)$. 
We will construct two positive functions $C(t)$ and $D(t)$ such that 
$$C(t)\le n(t,x),\;p(t,x)\le D(t)\;\;\;\mbox{a.e. in }\;(0,\infty)\times\Omega$$
where $C(\infty)$ and $D(\infty)$ are independent
of the initial conditions $(n^0,p^0)$. With these bounds we can easily
construct an absorbing set.

We make the following assumptions in addition to the general Assumption (A). 

\begin{description}
\item[{Assumption (A2):}] {\em  $\ds\inf_{s\ge 1} P'(s)\ge r_0>0$ and 
 $0\le Q(n,p)\le\overline{Q}_0$. }
\end{description}
The condition on $P$ yields $h_\infty=\infty$ and thus we have the time-independent upper
bound $\Gamma(d_0^\ast)$ from Theorem \ref{thm:exist} (but it depends on initial conditions).
It is satisfied for $P(s)=k s^{\gamma}$ when $\gamma\ge 1$.
We also remark that it is possible to relax the condition on $Q$ to the case 
of linear growth $Q(n,p)\le \overline{Q}_0+ \alpha(n+p)$ with 
sufficiently small $\alpha>0$. This can be easily seen in the following
upper bound estimates.

{\bf Upper Bounds.} Let 
\begin{equation}\label{dbd}
d_0=(\max\{\sup_\Omega n^0, \; \sup_\Omega p^0\}-1)^+
\quad\mbox{and}\quad
\bar d=\max\{\sup_{\Sigma_D}\bar n, \; \sup_{\Sigma_D}\bar p,\;1 \}.
\end{equation}
For $d\ge \bar d$ and $\lambda>0$ (to be determined below), set
$$D(t)=d(1+d_0 e^{-\lambda t})\quad\mbox{and}\quad
\phi_d=(n-D(t))^+,\;\;\;\psi_d=(p-D(t))^+.$$
Clearly $\phi_d$, $\psi_d\in V$ and $\phi_d(0)=\psi_d(0)=0$. 
Since $(n, p, u)$ is in $L^\infty((0,T)\times\Omega)$ for any given $T>0$, we see that 
$\phi_d$ and $\psi_d$ are in $L^\infty((0,T)\times\Omega)$.
Use this $\phi_d$ as a test function in (\ref{w-n}). Since $n=\phi_d+D(t)$ and
$\nabla n=\nabla\phi_d$ when $\phi_d\ge 0$ and 
$\langle\frac{\partial}{\partial t}\phi_d,\phi_d\rangle=\frac{d}{dt}|\phi_d|_2^2$,  we have
\begin{equation}\label{nd1}
\frac{d}{dt}(\frac{1}{2}|\phi_d|_2^2)+\langle D'(t),\phi_d\rangle
+\mu_n r_0 |\nabla\phi_d|_2^2\le \mu_n \langle n\nabla u,\nabla\phi_d\rangle
  +\langle G(n, p),\phi_d\rangle .
\end{equation}
Since $n\nabla\phi_d=(\phi_d+D(t))\nabla\phi_d
=\nabla(\frac{1}{2}\phi_d^2+D(t)\phi_d)$, we use 
$\frac{1}{2}\phi_d^2+D(t)\phi_d\in V$ as a test function in (\ref{w-u}) to obtain
$$\langle n\nabla u,\nabla\phi_d\rangle
=\langle\nabla u,\nabla (\frac{1}{2}\phi_d^2+D(t)\phi_d)\rangle
=-\frac{1}{\epsilon_0}\langle n-p-N, \frac{1}{2}\phi_d^2+D(t)\phi_d\rangle .$$
For $G$ from either (\ref{G0}) or (\ref{G}), we have 
$G(n,p)\le g+Q(n,p)\le g+\overline{Q}_0$ by Assumption (A2).
Thus (\ref{nd1}) becomes
\begin{equation}\label{nd2}
\frac{d}{dt}(\frac{1}{2\mu_n}|\phi_d|_2^2)+ r_0 |\nabla\phi_d|_2^2
\le -\frac{1}{\epsilon_0}\langle n-p-N, \frac{1}{2}\phi_d^2+D(t)\phi_d\rangle
  +\frac{1}{\mu_n}\langle Q+g-D'(t),\phi_d\rangle .
\end{equation}
For the $p$-component we obtain similarly
\begin{equation}\label{pd2}
\frac{d}{dt}(\frac{1}{2\mu_p}|\psi_d|_2^2)+r_0 |\nabla\psi_d|_2^2
\le \frac{1}{\epsilon_0}\langle n-p-N, \frac{1}{2}\psi_d^2+D(t)\psi_d\rangle
  +\frac{1}{\mu_p}\langle Q+g-D'(t),\psi_d\rangle .
\end{equation}
It is straightforward to verify that
$$(n-p)(\phi_d-\psi_d)\ge 0,\;\;\;
(n-p-N)(\phi_d^2-\psi_d^2)\ge -\frac{N^2}{4}(\phi_d+\psi_d).$$
Hence combining (\ref{nd2})--(\ref{pd2}) yields
\begin{equation}\label{nd3}
\begin{array}{rl}
&\ds\frac{d}{dt}(\frac{1}{2\mu_n}|\phi_d|_2^2+\frac{1}{2\mu_p}|\psi_d|_2^2)
+r_0 (|\nabla\phi_d|_2^2 +|\nabla\psi_d|_2^2)\\[5pt]
\le &\ds \frac{1}{8\epsilon_0}\langle N^2, \phi_d+\psi_d \rangle
  +\frac{D(t)}{\epsilon_0}\langle N, \phi_d-\psi_d\rangle
  +\langle Q+g-D'(t),\frac{1}{\mu_n}\phi_d+\frac{1}{\mu_p}\psi_d\rangle .
\end{array}
\end{equation}
The following is similar to the proof of Lemma \ref{lm:Linfty}, 
but we need to use Lemma \ref{lm:stamp2} instead of Lemma \ref{lm:stamp} 
because the upper bound $D(t)$ is a function of $t$.
Denote 
$$\Omega_d^{(n)}(t)=\{x\in\Omega:\;n(t,x)>D(t)\},\;\;\;
\Omega_d^{(p)}(t)=\{x\in\Omega:\;p(t,x)>D(t)\}$$ 
and
$$\varphi(d)=\sup_{t\ge 0}\left(|\Omega_d^{(n)}(t)|+|\Omega_d^{(p)}(t)|\right).$$
Clearly $\varphi$ is nonnegative,  non-increasing in $d$.
Note that, for $f\in L^\infty$,
\begin{eqnarray*}
\langle f,\phi_d\rangle&\le &|f|_\infty |\Omega_d^{(n)}|^{3/4}|\phi_d|_4 \\
&\le& |f|_\infty  |\Omega_d^{(n)}|^{3/4} \alpha_4 |\nabla\phi_d|_2
\le \frac{r_0}{2} |\nabla\phi_d|_2^2 
 +\frac{\alpha_4^2}{2 r_0}|f|_\infty^2|\Omega_d^{(n)}|^{3/2}
\end{eqnarray*}
and similarly for $\psi_d$. Hence by properly identifying $f$ in (\ref{nd3}) 
we obtain
\begin{eqnarray*}
\lefteqn{ \frac{d}{dt}(\frac{1}{\mu_n}|\phi_d|_2^2+\frac{1}{\mu_p}|\psi_d|_2^2)
+r_0 (|\nabla\phi_d|_2^2 +|\nabla\psi_d|_2^2) }\\
&\le & \frac{\alpha_4^2}{r_0}\left[\frac{|N|^2_\infty}{8\epsilon_0}
       +\frac{\overline{Q}_0+|g|_\infty}{\min(\mu_n,\mu_p)}
+ \left(\frac{|N|_\infty}{\epsilon_0}
    +\frac{\lambda}{\min(\mu_n,\mu_p)}\right) D(t) \right]^2 \times \\
&& \left( |\Omega_d^{(n)}|^{3/2}+ |\Omega_d^{(p)}|^{3/2} \right)\,,
\end{eqnarray*}
where we have also used the fact that $-D'(t)\le \lambda D(t)$. Using the
Poincar\'{e} inequality again and the fact that $D(t)\ge 1$, we reduce 
the above further into
\begin{equation}\label{nd4}
\frac{d}{dt}\left( |\phi_d|_2^2+|\psi_d|_2^2\right)
+\frac{r_0\max(\mu_n,\mu_p)}{\alpha_2^2}\left(|\phi_d|_2^2+|\psi_d|_2^2\right)
\le K_1 D^2(t)\varphi(d)^{3/2}
\end{equation}
where we denote
\begin{equation}\label{K1}
K_1=\frac{\alpha_4^2\max(\mu_n,\mu_p)}{r_0}
 \left(\frac{|N|^2_\infty+8|N|_\infty}{8\epsilon_0}
  +\frac{\overline{Q}_0+|g|_\infty+\lambda}{\min(\mu_n,\mu_p)}\right)^2.
\end{equation}
Choose 
\begin{equation}\label{lambda}
\lambda=\frac{r_0\max(\mu_n,\mu_p)}{3\alpha_2^2}
\end{equation}
(independent of initial data), and multiply (\ref{nd4}) with $e^{3\lambda t}$
and then integrate to yield
$$|\phi_d|_2^2+|\psi_d|_2^2
\le K_1 \varphi(d)^{3/2} e^{-3\lambda t}\int_0^t e^{3\lambda s} D^2(s) ds
\le \frac{K_1}{\lambda} \varphi(d)^{3/2} D^2(t).$$
Now since for $\hat d> d\ge\bar d$,
$$|\phi_d|_2^2\ge\int_{\Omega_{\hat d}^{(n)}}|\phi_d|^2 dx
\ge (\hat d-d)^2(1+d_0e^{-\lambda t})^2 |\Omega_{\hat d}^{(n)}|
 =d^{-2}(\hat{d}-d)^2|\Omega_{\hat d}^{(n)}| D^2(t)$$
and similarly for $|\psi_d|_2^2$, we obtain from the above that
$$ |\Omega_{\hat d}^{(n)}(t)|+|\Omega_{\hat d}^{(p)}(t)|
 \le \frac{K_1}{\lambda} d^2 (\hat d- d)^{-2} \varphi(d)^{3/2}.$$
Hence,
$$\varphi(\hat d)
\le\frac{K_1}{\lambda} d^2 (\hat d- d)^{-2} \varphi(d)^{3/2}$$
for $\hat d>d\ge\bar d$. By applying Lemma \ref{lm:stamp2} to $\varphi(d)$,
we conclude that $\varphi(d)=0$ for
\begin{equation}\label{d}
d=2\bar d(1+2^8 K_1^{3/2}\lambda^{-3/2}|\Omega|^{3/4}). 
\end{equation}
Hence $\phi_d=\psi_d=0$ a.e. for $t\ge 0$ and $x\in\Omega$. That is,
$$n(t,x),\;\;p(t,x)\le d(1+d_0 e^{-\lambda t})$$
where $d$ in (\ref{d}) and $\lambda$ in (\ref{lambda}) are independent of the
initial data.

{\bf Lower Bounds.} Let
\begin{equation}\label{c0bc}
c_0=\min\{1/d_0,\;\inf_{\Omega} n^0,\;\inf_{\Omega}p^0\},\;\;\;
\bar c=\min\{\inf_{\Sigma_D}\bar n,\;\inf_{\Sigma_D}\bar p,\;1\}.
\end{equation}
For $C(t)$ to be determined below such that $C(t)\le\bar c$ and $C(0)\le c_0$,
let
$$ \phi_c=(n-C(t))^-,\;\;\;\psi_c=(p-C(t))^-\in V$$
Choose $\phi_c$ and $\psi_c$ as test functions in (\ref{w-n})--(\ref{w-p}) 
and follow the same calculations as in (\ref{nd1}) to (\ref{pd2}) above to 
obtain (note that here $\phi_c$ and $\psi_c$ are non-positive)
$$\frac{d}{dt}(\frac{1}{2\mu_n}|\phi_c|_2^2)
\le -\frac{1}{\epsilon_0}\langle n-p-N, \frac{1}{2}\phi_c^2+C(t)\phi_c\rangle
  +\frac{1}{\mu_n}\langle -G(n, p)+C'(t),-\phi_c\rangle$$
and
$$\frac{d}{dt}(\frac{1}{2\mu_p}|\psi_c|_2^2)
\le \frac{1}{\epsilon_0}\langle n-p-N, \frac{1}{2}\psi_c^2+C(t)\psi_c\rangle
 +\frac{1}{\mu_p}\langle -G(n, p)+C'(t),-\psi_d\rangle .$$
Using the fact that $(n-p)(\phi_c-\psi_c)\ge 0$ and $(n-p)(\phi_c^2-\psi_c^2)\ge 0$,
we combine the two above to yield
\begin{eqnarray} 
\frac{d}{dt}(\frac{1}{2\mu_n}|\phi_c|_2^2+\frac{1}{2\mu_p}|\psi_c|_2^2)
&\le & \frac{1}{2\epsilon_0}\langle N, \phi_c^2-\psi_c^2 \rangle
  +\langle \frac{C(t)}{\epsilon_0} N, \phi_c-\psi_c\rangle  \nonumber\\
&&  +\langle -G(n, p)+C'(t),-\frac{1}{\mu_n}\phi_c-\frac{1}{\mu_p}\psi_c\rangle .
\label{nc1}\end{eqnarray}

Our choice of $C(t)$ will depend on the form of $G(n,p)$, and we need
to make additional assumptions. Hence we discuss
the two models (\ref{G0}) and (\ref{G}) separately below.

\underline{\em Case 1: $G$ is by (\ref{G0})}. In this case we assume
further that
\begin{equation}\label{SRH}
\inf_{n, p\ge 0}(1+n+p)Q(n,p)\ge\bar{q}>0.
\end{equation}
This is true for the Shockley-Read-Hall model $Q_{\rm SRH}\sim (1+n+p)^{-1}$, 
and it also includes the case when $Q$ has a positive lower bound.
We choose
\begin{equation}\label{C-G0}
C(t)=\frac{c}{1+e^{-\omega t}/c_0}
\;\;\;\mbox{with}\;\;\omega=\min\{c_0,\lambda\}
\end{equation}
for $0<c\le\bar c$, where $c_0$ and $\bar c$ are as in (\ref{c0bc}) and
$\lambda$ in (\ref{lambda}). Thus $C'(t)\le C(t)e^{-\omega t}\omega/c_0\le C(t)$.
By (\ref{SRH}) on $Q$,
$$Q(n,p)\ge \frac{\bar{q}}{1+n+p}\ge \frac{\bar{q}}{2D(t)}
\ge \frac{\bar{q}}{2d(1+e^{-\omega t}/c_0)}=\frac{\bar{q}}{2cd}C(t)$$
since $d_0\le 1/c_0$ and $\omega\le\lambda$.
When $n<C(t)$, since $p\le D(t)$ from the upper bound, we have
$$np\le C(t)D(t)\le c d\le 1/2\;\;\;\mbox{if we choose}\;c\le 1/2d.$$
Hence,
$$\langle -G(n, p), -\phi_c\rangle=\langle -g+Q(np-1),-\phi_c\rangle\le\langle -\frac{1}{2}Q,-\phi_c\rangle
\le \langle-\frac{\bar{q}}{4cd} C(t), -\phi_c\rangle .$$
Similarly we can obtain
$$\langle -G(n, p), -\psi_c\rangle\le \langle -\frac{\bar{q}}{4cd} C(t), -\psi_c\rangle .$$
Substituting these estimates in to (\ref{nc1}) we obtain
\begin{equation}\label{nc2}
\frac{d}{dt}(\frac{1}{2\mu_n}|\phi_c|_2^2+\frac{1}{2\mu_p}|\psi_c|_2^2)
\le  \frac{|N|_\infty}{2\epsilon_0}(|\phi_c|_2^2+|\psi_c|_2^2 )
  +C(t)\langle F, -\frac{1}{\mu_n}\phi_c-\frac{1}{\mu_p}\psi_c\rangle
\end{equation}
where
$$F=\frac{\max(\mu_n,\mu_p)}{\epsilon_0}|N|_\infty
 +1-\frac{\bar{q}}{4cd}.$$
Therefore when we choose
\begin{equation}\label{c}
c=\min\left\{\bar c,\; 1/2d,\;
 \frac{\bar{q}}{4d(\max(\mu_n,\mu_p)|N|_\infty/\epsilon_0+1)}\right\},
\end{equation}
we have $F\le 0$. Hence 
\begin{equation}\label{nc3}
\frac{d}{dt}(\frac{1}{2\mu_n}|\phi_c|_2^2+\frac{1}{2\mu_p}|\psi_c|_2^2)
\le \frac{\max(\mu_n,\mu_p)}{\epsilon_0}|N|_\infty
 (\frac{1}{2\mu_n}|\phi_c|_2^2+\frac{1}{2\mu_p}|\psi_c|_2^2)
\end{equation}
and Gronwall's inequality leads to $\phi_c=\psi_c=0$ a.e. for $t\ge 0$ and
$x\in\Omega$. Therefore $n(t,x)$, $p(t,x)\ge C(t)$ 
when $C(t)$ is given by (\ref{C-G0}) with $c$ in (\ref{c}).

\underline{\em Case 2: $G$ is by (\ref{G})}. Here we need to assume that
\begin{equation}\label{new2}
\inf_{0< s\le 1}P'(s)=s_0>0,\;\;\inf_{n, p\ge 0}Q(n,p)\ge \bar{q}_0>0\quad\mbox{and}\quad
S(-\infty)=0.
\end{equation}
The first condition implies that $h_0=\infty$, and when $S$ is an exponential
function, $S(-\infty)=0$. We need to choose a different form of $C(t)$:
\begin{equation}\label{C-G}
C(t)=\Gamma(-\theta-\theta_0 e^{-\hat\omega t})
\end{equation}
where $\theta_0$ and $\hat\omega$ are given by
\begin{equation}\label{theta0}
\theta_0=\max\{-H(c_0),\; K_0 d d_0\},\;\;\;
\hat\omega=\min\{\frac{s_0}{\theta_0},\;\lambda\}
\;\;\;\mbox{with}\;\;\;K_0=\max_{d\le\xi\le d(1+d_0)}H'(\xi),
\end{equation}
and for $\theta\ge-H(\bar{c})$. Then $C(t)\le \Gamma(-\theta)\le\bar c$,
$C(0)\le\Gamma(-\theta_0)\le c_0$, and 
$$C'(t)=\Gamma'(-\theta-\theta_0 e^{-\hat\omega t}) 
   e^{-\hat\omega t}\theta_0\hat\omega
\le C(t)\frac{\theta_0}{s_0}\hat\omega\le C(t)$$
since $\Gamma$ is the inverse of $H$ and hence 
$\Gamma'(\cdot)=\Gamma(\cdot)/P'(\Gamma(\cdot))\le \Gamma(\cdot)/s_0$
when $\Gamma\le 1$ by (\ref{new2}). We also have
$$\displaylines{
H(D(t))=H(d+dd_0e^{-\lambda t})=H(d)+H'(\xi) d d_0 e^{-\lambda t}  \cr
\le H(d)+K_0 d d_0 e^{-\lambda t}\le H(d)+\theta_0 e^{-\hat\omega t},
}$$
and hence $H(C(t))+H(D(t))\le-\theta+H(d)$.
From (\ref{new2}), there exists $\bar k>0$ such that
$$S(-k)\le 1/2\;\;\;\mbox{for all}\;\;k\ge\bar k.$$
Hence, when $n< C(t)$, since $p\le D(t)$, 
$$ S(H(n)+H(p))\le S(-\theta+H(d))\le 1/2
\;\;\;\mbox{if we choose}\;\;\theta\ge \bar k + H(d).$$
Therefore,
\begin{eqnarray*}
\langle -G(n,p),-\phi_c\rangle&=&\langle -g+Q(S(H(n)+H(p))-1),-\phi_c\rangle\\
&\le& \langle -Q/2,-\phi_c\rangle\le \langle-\bar{q}_0/2,-\phi_c\rangle
\end{eqnarray*}
and similarly for $\langle-G,-\psi_c\rangle$. Substituting these estimates in (\ref{nc1})
and following a similar calculation as in the previous case from (\ref{nc2})
to (\ref{nc3}), we can obtain
$n(t,x)$, $p(t,x)\ge C(t)$ for $C(t)$ given in (\ref{C-G}) with 
\begin{equation}\label{theta}
\theta=\max\{-H(c), \bar{k}+H(d), 
-H(\frac{\bar{q}_0}{2(\max(\mu_n,\mu_p)|N|_\infty/\epsilon_0+1)})\}.
\end{equation}

In summary, we have established for global solutions $(n,p)$ the
upper bound $D(t)$ given by
\begin{equation}\label{Dt}
D(t)=d(1+d_0 e^{-\lambda t})
\end{equation}
where $d$ is independent of the initial data and given in (\ref{d}),
$d_0$ in (\ref{dbd}), and $\lambda$ in (\ref{lambda}), and the lower 
bound $C(t)$ given by
\begin{equation}\label{Ct}
C(t)=\left\{\begin{array}{ll}
    \ds\frac{c}{1+e^{-\omega t}/c_0}\;\;\;&
        \mbox{when $G$ is by (\ref{G0}) with (\ref{SRH})}\\[5pt]
    \ds \Gamma(-\theta-\theta_0 e^{-\hat\omega t})\;\;
        &\mbox{when $G$ is by (\ref{G}) with (\ref{new2})}
    \end{array}\right.
\end{equation}
where $c$ and $\theta$ are independent of the initial data and given in
(\ref{c}) and (\ref{theta}) respectively, and $c_0$ in (\ref{c0bc}),
$\omega$ in (\ref{C-G0}), and $\theta_0$ and $\hat\omega$ in (\ref{theta0}).

\begin{theorem}\label{thm:exist1} {\rm (Absorbing $L^\infty$ bounds.)}
Suppose Assumptions (A) and (A2) hold, and $G$ is given either by
(\ref{G0}) with the condition (\ref{SRH}) or by (\ref{G}) with (\ref{new2}). 
Then for any global $L^\infty\cap H^\infty$ solutions $(n,p,u)$ to 
(\ref{w-n})--(\ref{w-u}) from initial data $(n^0, p^0)$, there 
hold the following bounds:
\begin{equation}\label{Linfty-s}
C(t)\le n(t,x),\; p(t,x)\le D(t)
\end{equation}
a.e. for $t\ge 0$ and $x\in\Omega$, where $D(t)$ and $C(t)$ are
given in (\ref{Dt}) and (\ref{Ct}), respectively. In particular 
the positive constants $C(\infty)$ and $D(\infty)$ are independent of the 
initial data $(n^0,p^0)$. 
\end{theorem}

\begin{remark}\label{rmk:exist1} 
The extra conditions (\ref{SRH}) or (\ref{new2})
are only needed for establishing the lower bound $C(t)$.
\end{remark}

\vspace{.1in}

Now we study the system (\ref{w-n})--(\ref{w-p}) as a dynamical system acting
on the set of admissible initial data
$$H_{\rm ad}=\{(n^0,p^0)\in L^\infty(\Omega)^2: \;\inf_{\Omega} n^0>0,\;\;
\inf_{\Omega} p^0>0\}.$$
We will follow the framework of analyzing infinite dimensional dynamical systems 
given in \cite{Tm}. 
Let ${\cal S}(t)$ ($t\ge 0$) denote the semigroup of solution maps defined on 
$H_{\rm ad}$, i.e., ${\cal S}(t)(n^0,p^0)=(n(t), p(t))$ where $(n(t),p(t))$ is 
a solution to (\ref{w-n})--(\ref{w-p}) with initial data $(n^0,p^0)$. 
From Theorem \ref{thm:exist}, we see for each $t>0$, ${\cal S}(t)$ is well-defined 
from $H_{\rm ad}$ into itself. In general it may be set-valued, and under 
Assumption (A1), it becomes single-valued. 

We first introduce a few definitions. A set $B\subset H_{\rm ad}$ is said to be
$L^\infty$-bounded in $H_{\rm ad}$ if there exist positive constants $M_1$ and 
$M_2$ such that 
$$M_1\le n^0(x),\;p^0(x)\le M_2\;\;\;\mbox{for all}\;\;\;(n^0,p^0)\in B.$$
A set ${\cal B}\subset H_{\rm ad}$ is said to be an absorbing set for ${\cal S}(t)$ 
if it ``absorbs'' all $L^\infty$-bounded sets in $H_{\rm ad}$; that is, 
for any $L^\infty$-bounded set $B\subset H_{\rm ad}$, there exists $t_B>0$ 
(depending on $B$) such that ${\cal S}(t) B\subset {\cal B}$ for all $t\ge t_B$.

Given any $L^\infty$-bounded set $B\subset H_{\rm ad}$, suppose $M_1$ and
$M_2$ are the two bounds as defined above. 
From the $L^\infty$ estimates established for the solutions $(n,p)$ in 
Theorem~\ref{thm:exist1}, and the way $d_0$ and $c_0$ or $\theta_0$ are 
chosen in (\ref{dbd}) and (\ref{c0bc}) or (\ref{theta0}), 
we see that, for each $t>0$,
\begin{equation}\label{SB}
 \frac{c}{1+e^{-\omega t}/m_1}\;\;\mbox{or}\;\;
 \Gamma(-\theta-\hat{m}_1 e^{-\hat\omega t})\le {\cal S}(t)(n^0,p^0)\le d(1+m_2 e^{-\lambda t})
\;\;\;\mbox{a.e. in}\;\;\Omega
\end{equation}
for all $(n^0,p^0)\in B$, where 
$$m_2=(M_2-1)^+,\;\;
m_1=\min\{1/m_2, M_1\}\;\;\mbox{or}\;\;\hat{m}_1=\max\{-H(m_1), K d M_2\}$$
(with $K=\max_{d\le\xi\le d M_2} H'(\xi)$.) Set 
$$t_B=\max\{\ln m_2/\lambda, -\ln m_1/\omega\}\;\;\mbox{or}\;\;
\max\{\ln m_2/\lambda, \ln\hat{m}_1/\hat{\omega}\},$$
then for all $t\ge t_B$, there holds
$$\frac{c}{2} \;\mbox{or}\;\Gamma(-\theta-1)\le {\cal S}(t)(n^0,p^0)\le 2d.$$
Thus we have an absorbing set as stated in the following theorem.

\begin{theorem}\label{thm:absorb} {\rm (Existence of an absorbing set.)}
Under the assumptions of Theorem \ref{thm:exist1}, 
the dynamical system $\{{\cal S}(t)\}_{t\ge 0}$ given by
(\ref{w-n})--(\ref{w-p}) has an absorbing set $\cal B$ defined by
\begin{equation}\label{absorb}
{\cal B}=\{(n,p)\in H_{\rm ad}:\; c_\ast \le n(x),\;p(x)\le 2 d\}
\end{equation}
where $c_\ast=c/2$ if $G$ is by (\ref{G0}) and $c_\ast=\Gamma(-\theta-1)$
if $G$ is by (\ref{G}), and $d$, $c$, or $\theta$ are as given in 
(\ref{d}), (\ref{c}) or (\ref{theta}).
\end{theorem}

\subsection{Existence of a Global Compact Attractor}

In this subsection we study the existence of an attractor for the dynamical 
system under the assumptions in Theorem \ref{thm:exist1} and (A1). 
With these assumptions the semigroup $\{{\cal S}(t)\}$ is of single-valued, and there is an 
absorbing set as constructed in (\ref{absorb}). 

We also equip $H_{\rm ad}$ with the $V^\ast$-topology; i.e., with respect
to the norm 
$$\|(n^0,p^0)\|_\ast^2
=\langle n^0,(-\Delta_0)^{-1} n^0\rangle_{V^\ast\times V}+
\langle p^0,(-\Delta_0)^{-1} p^0\rangle_{V^\ast\times V}.$$
Then from (\ref{gron}) in the proof of Theorem~\ref{thm:unique}, 
we have the Lipschitz continuity of ${\cal S}(t)$ on $H_{\rm ad}$ with 
respect this $V^\ast$-norm. 

\begin{lemma}\label{lm:cont}
For each $t\ge 0$, ${\cal S}(t)$ is Lipschitz continuous from $H_{\rm ad}$ to 
$H_{\rm ad}$ with respect to the $V^\ast$-norm.
\end{lemma}

\begin{lemma}\label{lm:compact}
For each $t>0$ and each $L^\infty$-bounded set $B\in H_{\rm ad}$, 
$\overline{{\cal S}(t)B}$ is an $L^\infty$-bounded set in $H_{\rm ad}$, 
where the closure is taken with the $V^\ast$-topology in $H_{\rm ad}$.
Moreover $\overline{{\cal S}(t)B}$ is compact in $H_{\rm ad}$.
\end{lemma}
{\em Proof.} Suppose $B$ is an $L^\infty$-bounded set in $H_{\rm ad}$.
Then from (\ref{SB}) we see that ${\cal S}(t)B$ is $L^\infty$-bounded in $H_{\rm ad}$. 
To show that its closure in the $V^\ast$-topology is also $L^\infty$-bounded,
we assume $z=\lim_{i\to\infty} z_i$ where each $z_i\in {\cal S}(t)B$.
Since $\{z_i\}$ is bounded in $L^\infty$, there exists 
${\tilde z}\in L^\infty$ and a subsequence $\{z_{i_k}\}$ such that 
$z_{i_k}$ converges weakly-$\ast$ to ${\tilde z}$.  Clearly the same
bounds as in (\ref{SB}) hold for $\tilde z$. For any $\phi\in V^2$,
$\langle z_{i_k},\phi\rangle\to \langle{\tilde z},\phi\rangle$. Hence
$\langle{\tilde z}-z,\phi\rangle_{V^\ast\times V}=0$ for all $\phi\in V^2$.
Thus $z={\tilde z}$. Therefore $z\in L^\infty$ and the same bounds
as in (\ref{SB}) hold for $z$. That is, $\overline{{\cal S}(t)B}$ is also
an $L^\infty$-bounded set.
Since $H_{\rm ad}$ is equipped with the $V^\ast$-topology, and $L^2$ is 
compactly embedded in $V^\ast$, from the bounds in (\ref{SB})
we see immediately that ${\cal S}(t)B$ is precompact in $H_{\rm ad}$. $\Box$

Define the $\omega$-limit set of the absorbing set $\cal B$:
\begin{equation}\label{AA}
{\cal A}=\omega({\cal B})
=\bigcap_{s\ge 0}\;\overline{\bigcup_{t\ge s} {\cal S}(t){\cal B}}
\end{equation}
where the closure is taken with the $V^\ast$-topology in $H_{\rm ad}$.
Then from Lemma \ref{lm:compact}, we see that $\cal A$ is $L^\infty$-bounded
and compact in $H_{\rm ad}$. As an $\omega$-limit set of a nonempty set, 
$\cal A$ is nonempty and invariant: 
${\cal S}(t){\cal A}={\cal A}$ for all $t\ge 0$ (see e.g. Lemma 1.1 in \cite{Tm}). 

We claim that $\cal A$ attracts every $L^\infty$-bounded set in $H_{\rm ad}$;
that is, for any $L^\infty$-bounded set $B\subset H_{\rm ad}$, 
$${\rm dist}({\cal S}(t)B, {\cal A})\to 0\;\;\;\mbox{as}\;\;t\to\infty$$
where the distance is by the $V^\ast$-norm:
$${\rm dist}(X, Y)=\sup_{x\in X}\inf_{y\in Y} \|x-y\|_\ast.$$
We argue by contradiction. If otherwise, there is an $L^\infty$-bounded 
set $B_0\subset H_{\rm ad}$ such that ${\rm dist}({\cal S}(t)B_0,{\cal A})$ does not 
tend to 0 as $t\to\infty$.  That is, there exist $\delta>0$, $t_i\to\infty$ 
and $z_i\in B_0$ such that
\begin{equation}\label{zi}
{\rm dist}({\cal S}(t_i)z_i,{\cal A})\ge\delta>0
\end{equation}
for all $i=1,2,\cdots$. By (\ref{SB}), there
is $t_0$ such that ${\cal S}(t)B_0\subset {\cal B}$ when $t\ge t_0$. 
Without loss of generality, we assume $t_i\ge 2 t_0$ for all $i=1,2,\cdots$. 
Then the sequence ${\cal S}(t_i-t_0)z_i$ is contained in ${\cal B}$, and 
${\cal S}(t_i)z_i={\cal S}(t_0){\cal S}(t_i-t_0)z_i$ is precompact in $H_{\rm ad}$ 
by Lemma \ref{lm:compact}. Thus ${\cal S}(t_i)z_i$ has a convergent 
subsequence ${\cal S}(t_{i_k})z_{i_k}\to z_\ast\in H_{\rm ad}$ as $k\to\infty$. Note
that $z_\ast=\lim_{k\to\infty} {\cal S}(t_{i_k}-t_0)y_k$ with 
$y_k={\cal S}(t_0) z_{i_k}\in {\cal B}$. Hence $z_\ast\in{\cal A}$ by the definition of
$\omega$-limit set. This contradicts (\ref{zi}).
Therefore $\cal A$ attracts sets that are $L^\infty$-bounded in $H_{\rm ad}$.

$\cal A$ is also maximal among such attractors. 
Suppose there is another $L^\infty$-bounded attractor ${\cal A}_1\supset {\cal A}$. Since ${\cal B}$ absorbs $L^\infty$-bounded set, ${\cal S}(t){\cal A}_1\subset {\cal B}$ 
for $t$ large enough. Hence ${\cal A}_1={\cal S}(t){\cal A}_1\subset {\cal B}$, and 
${\cal A}_1=\omega({\cal A}_1)\subset\omega({\cal B})={\cal A}$. 
Therefore, ${\cal A}_1={\cal A}$. 

Finally, we show $\cal A$ is connected. 
Suppose otherwise. Then there are open sets $O_1$ and $O_2$
in $H_{\rm ad}$ such that $O_1\cap O_2=\emptyset$ and 
${\cal A}\subset O_1\cup O_2$.
Let $K=\mbox{co}{\cal A}$ be the convex hull of $\cal A$. Then $K$ is
connected. $K\subset H_{\rm ad}$ because $H_{\rm ad}$ is convex, and $K$ is 
$L^\infty$-bounded since $\cal A$ is. Note that
${\cal A}={\cal S}(t){\cal A}\subset {\cal S}(t)K$ and ${\cal S}(t)K$ is also connected since 
${\cal S}(t)$ is continuous for each $t>0$ (Lemma \ref{lm:cont}). 
Therefore $O_1\cup O_2$ does not cover ${\cal S}(t)K$ for each $t>0$.
That is, there are $z_j\in K$ such that $S(j)z_j\not\in O_1\cup O_2$ for all
$j=1,2,\cdots$. Since it is easy to see that $S(j)z_j$ is precompact in 
$H_{\rm ad}$ for $j$ large enough (by Lemma \ref{lm:compact}), 
there exists a subsequence of $S(j)z_j$ that
converges to some $y\in H_{\rm ad}$. Clearly $y\not\in O_1\cup O_2$.
On the other hand, since $\cal A$ attracts $K$, $y$ must be in $\cal A$.
This is a contradiction. Hence $\cal A$ is connected.

Note that for $(n^0,p^0)\in H_{\rm ad}$ and $t>0$, 
${\cal S}(t)(n^0,p^0)\in V^2$ from our existence results in Theorem~\ref{thm:exist}. 
Hence ${\cal A}={\cal S}(t){\cal A}$ is contained in $V^2$.

Thus we have proved the following theorem:

\begin{theorem}\label{thm:attract} {\rm (Existence of an attractor.)}
In addition to the assumptions of Theorem \ref{thm:exist1}, suppose 
Assumption (A1) holds. Then the dynamical system $\{{\cal S}(t)\}$
defined by (\ref{w-n})--(\ref{w-p}) has a global attractor $\cal A$
given in (\ref{AA}) that attracts all $L^\infty$-bounded sets in
$H_{\rm ad}$. Moreover, $\cal A$ is maximal among such attractors,
and $\cal A$ is $L^\infty$-bounded, compact, connected, and contained in $V^2$.
\end{theorem}

We make two remarks on our results above.
(1) Our proof here is similar to that of \cite[Theorem 1.1]{Tm},
which cannot be applied directly because the set ${\cal B}$ absorbs 
only sets that are $L^\infty$-bounded, not in terms of the $V^\ast$-topology.
(2) Using the $L^\infty$ estimates we established for the solutions and
the continuity of the solution map ${\cal S}(t)$ with respect to the $V^\ast$-topology,
we can improve the continuity of ${\cal S}(t)$ to an interpolation space between
$H^{-1}$ and $L^2$, and the discussion above can be easily carried over to that
setting. 

%%%%%%%%%%%%%%%%    

\section{Existence of Stationary Solutions}\label{sec:stat}\nopagebreak
\setcounter{equation}{0}    

In this section we study the stationary version of the model equations 
(\ref{n})--(\ref{p}) and (\ref{u}) with boundary conditions (\ref{bc}). 
We are looking for a triplet $(n,p,u)\in (\bar n,\bar p,\bar u)+V^3$ 
with $n>0$, $p>0$ satisfying
\begin{equation}\label{n0}
\mu_n\langle \nabla P(n)-n\nabla u,\nabla\phi\rangle=\langle G(n,p), \phi\rangle
\end{equation}
\begin{equation}\label{w0}
\mu_p\langle \nabla P(p)+p\nabla u,\nabla\psi\rangle=\langle G(n,p),\psi\rangle
\end{equation}
\begin{equation}\label{u0}
\epsilon_0\langle\nabla u,\nabla\chi\rangle+\langle n-p-N, \chi\rangle=0
\end{equation}
for all $(\phi, \psi, \chi)\in V^3$, where $G(n,p)$ is given by either (\ref{G0}) or (\ref{G}).

The proof for existence consists of similar arguments as in \S\ref{sec:dis}
and \S\ref{sec:Linfty}. Consider the solution map ${\cal T}_0$ defined
from $(\tilde v,\tilde w,\tilde u)\in L^2(\Omega)^3$ to 
$(v,w,u)\in (\bar v,\bar w,\bar u)+V^3$ in
$$\mu_n\langle \tilde n \nabla v,\nabla\phi\rangle
=\langle G(\Gamma_0(v+\tilde u), \tilde p), \phi\rangle$$
$$\mu_p\langle \tilde p \nabla w,\nabla\psi\rangle
 =\langle G(\tilde n, \Gamma_0(w-\tilde u)),\psi\rangle$$
$$\epsilon_0\langle\nabla u,\nabla\chi\rangle
+\langle \Gamma_0(\tilde v+u)-\Gamma_0(\tilde w-u)-N,\chi\rangle=0 $$
for all $(\phi, \psi, \chi)\in V^3$, where
$$\Gamma_0(s)=\left\{\begin{array}{ll}
\Gamma(-c_0)\;\; &\mbox{when}\;\; s<-c_0\\[5pt]
\Gamma(s) & \mbox{when}\;\; -c_0\le s\le d_0\\[5pt]
\Gamma(d_0) & \mbox{when}\;\; s>d_0 \end{array}\right.$$
($c_0\in (0,h_0)$ and $d_0\in (0, h_\infty)$ to be determined)
and $\tilde n=\Gamma_0(\tilde v+\tilde u)$, 
$\tilde p=\Gamma_0(\tilde w-\tilde u)$.

Since the equations above are all decoupled, and each defines a maximum
monotone operator, the well-definedness of ${\cal T}_0$ is obvious.
The compactness and continuity of ${\cal T}_0$ can be shown similarly as 
for $\cal T$ in \S\ref{sec:dis}. Therefore ${\cal T}_0$ has a fixed point 
$(v,w,u)\in (\bar v, \bar w, \bar u)+V^3$ satisfying
\begin{equation}\label{v-00}
\mu_n\langle \Gamma_0(v+u) \nabla v,\nabla\phi\rangle
=\langle G(\Gamma_0(v+u), \Gamma_0(w-u)), \phi\rangle
\end{equation}
\begin{equation}\label{w-00}
\mu_p\langle \Gamma_0(w-u) \nabla w,\nabla\psi\rangle
 =\langle G(\Gamma_0(v+u), \Gamma_0(w-u)),\psi\rangle
\end{equation}
\begin{equation}\label{u-00}
\epsilon_0\langle\nabla u,\nabla\chi\rangle
+\langle \Gamma_0(v+u)-\Gamma_0(w-u)-N,\chi\rangle=0
\end{equation}
for all $(\phi, \psi, \chi)\in V^3$.

Now we establish the $L^\infty$ bounds for $v+u$ and $w-u$, 
and hence choose the truncation constants $d_0$ and $c_0$ accordingly. 

{\bf The Upper Bounds.} Here we need to assume $h_\infty$ is sufficiently large
($h_\infty>d_0$ given in (\ref{d0}) below.) 
Set $\hat d\in (0, h_\infty)$ such that
$$\Gamma(\hat d)=\max\{\sup_{\Sigma_D}\bar n,\;\sup_{\Sigma_D}\bar p,\; 1 \}$$
and for $d\in (\hat d, h_\infty)$ set $\phi_d=(u+v-d)^+\in V$ 
as the test function in both (\ref{v-00}) and (\ref{u-00}) to obtain
$$
\mu_n\Gamma(d) |\nabla\phi_d|_2^2 \le
\langle Q+g+\frac{\mu_n}{\epsilon_0}\Gamma(d) N, \phi_d\rangle .$$
Following the same argument as for the bounds of $v+u$ in 
\S\ref{sec:Linfty}, we can conclude that
$u(x)+v(x)\le d_0$ a.e. in $\Omega$ where
\begin{equation}\label{d0}
d_0=\hat d +4\alpha_6^2|\Omega|^{1/6}
\left(\frac{\overline{Q}(1+2\Gamma(\hat{d}))|\Omega|^{1/2}+|g|_2}{\Gamma(\hat d)\min(\mu_n,\mu_p)} 
+2\frac{|N|_2}{\epsilon_0}\right)
\end{equation}
provided that $d_0<h_\infty$. The same bound can be shown for $w(x)-u(x)$.

{\bf The Lower Bounds.} Here we also need to assume that $h_0>c_0$
where $c_0$ is to be given below in (\ref{c0}). 
Set $\hat{c}\in (0, h_0)$ such that
$$\Gamma(-\hat{c})=\min\{\inf_{\Sigma_D}\bar{n}, \inf_{\Sigma_D}\bar{p}\}.$$
Following the same argument as in \S\ref{sec:Linfty} for the lower bound, 
we can obtain $v+u\ge -c_0$ when
\begin{equation}\label{c0}
c_0=\hat{c}_0+4\frac{\alpha_6^2}{\epsilon_0}|\Omega|^{1/6}|N|_2
\end{equation}
provided that $h_0>c_0$, where
$$\hat{c}_0=\left\{\begin{array}{ll}
         \max\{\hat c, H(1/\Gamma(d_0))\}\;\;&\mbox{for $G$ in (\ref{G0})}\\
         \max\{\hat c, d_0-h_i\} &\mbox{for $G$ in (\ref{G})}
         \end{array}\right.$$

The same lower bound can be shown for $w(x)-u(x)$.

The $L^\infty$ bound estimates for $u$ can be obtained using the same argument
as in \S\ref{sec:Linfty}.

Hence when we choose the $d_0$ and $c_0$ as the truncation constants
in $\Gamma_0$, the fixed point $(v,w,u)$ of ${\cal T}_0$
is also a solution to (\ref{v-00})--(\ref{u-00}) with $\Gamma_0(\cdot)$ replaced
by $\Gamma(\cdot)$, and hence $n=\Gamma(v+u)$, $p=\Gamma(w-u)$ and $u$
form a solution to the stationary system (\ref{n0})--(\ref{u0}).

The above discussions can be summarized in the following theorem.

\begin{theorem}\label{thm:stat} {\rm (Existence of non-vacuum stationary solutions.)}
In addition to Assumption (A), suppose $h_\infty>d_0$ and $h_0>c_0$.
Then the stationary system (\ref{n0})--(\ref{u0}) has at least one solution
$(n, p, u)\in (\bar n, \bar p, \bar u)+V^3$ satisfying
$$\begin{array}{c}
\Gamma(-c_0)\le n(x),\;p(x)\le \Gamma(d_0),\\[5pt]
|u(x)|\le\sup_{\Sigma_D}\bar u + 
\ds 4\frac{\alpha_6^2}{\epsilon_0} (2\Gamma(d_0)|\Omega|^{1/2}+|N|_2)|\Omega|^{1/6}
\end{array}$$
a.e. in $\Omega$, where the positive constants $d_0$ and $c_0$ are as given in 
(\ref{d0})--(\ref{c0}). 
\end{theorem} 

\begin{remark}\label{rmk:stat}
The condition $h_0>c_0$ is needed only for the lower bounds of $(n,p)$. 
\end{remark}

In the case $G$ is given by (\ref{G0}), it is also possible to replace
the condition $h_0>c_0$ by a condition on the recombination coefficient
$Q$: $Q(n,p)>0$. For example, the Shockley-Read-Hall model 
$Q_{\rm SRH}\sim (1+n+p)^{-1}$ satisfies this condition. 
Then 
$$q_0=\min Q(n,p)\;\;\;\mbox{over}\;\;n, p\in [0,\Gamma(d_0)]$$ 
is positive, where $d_0$ is the upper bound 
in (\ref{d0}). For $c_0'\in (0,h_0)$ such that
\begin{equation}\label{c0'}
\Gamma(-c_0')=\min\{\Gamma(-\hat c),\; 
\frac{\epsilon_0 q_0}{\epsilon_0 q_0 \Gamma(d_0)+|N|_\infty \max(\mu_n,\mu_p)}\},
\end{equation}
set $\chi_c=(v+u+c_0')_-\in V$ as the test function in both
(\ref{v-00}) and (\ref{u-00}) to obtain
\begin{eqnarray*}
\mu_n\Gamma(-c_0')|\nabla\chi_c|_2^2 
&=& \langle Q (\Gamma(-c_0')\Gamma_0(w-u)-1)-g,-\chi_c\rangle  \\
&&+\frac{\mu_n}{\epsilon_0}\Gamma_0(-c_0')
  \langle \Gamma(-c_0')-\Gamma_0(w-u)-N,-\chi_c\rangle\\
&\le & \langle Q (\Gamma(-c_0')\Gamma(d_0)-1) 
  - \frac{\mu_n}{\epsilon_0}\Gamma(-c_0') N, -\chi_c\rangle\\
&\le &  \langle -\frac{q_0 |N|_\infty\max(\mu_n,\mu_p)}
 {\epsilon_0 q_0\Gamma(d_0)+|N|_\infty\max(\mu_n,\mu_p)}  \\
&& + \frac{\mu_n}{\epsilon_0}\Gamma(-c_0') |N|_\infty, 
   -\chi_c\rangle 
\;\le\;  0 
\end{eqnarray*}
Therefore, $\chi_c\equiv 0$ and hence $v(x)+ u(x)\ge -c_0'$ a.e. 
in $\Omega$. 

\begin{corollary}\label{cor:stat} 
Under Assumption (A), suppose $h_\infty>d_0$ and $G$ is given by (\ref{G0}) with
$Q(n,p)>0$. Then the stationary system (\ref{n0})--(\ref{u0}) has at least one solution
$(n, p, u)\in (\bar n, \bar p, \bar u)+V^3$ satisfying
$$\Gamma(-c_0')\le n(x),\;\;p(x)\le \Gamma(d_0)$$
a.e. in $\Omega$, where $d_0$ is from (\ref{d0}) and $c_0'$ from (\ref{c0'}). 
\end{corollary} 

Although there is no rigorous proof, the standard stationary drift-diffusion model 
is believed to have multiple solutions (see, e.g., \cite{MRS,J0}).
Hence we do not expect the uniqueness of solutions to be true for 
this nonlinear model.

%%%%%%%%%%%%%%%%%%%

\section{Vacuum Solutions}\label{sec:vacuum}\nopagebreak
\setcounter{equation}{0}
So far we have considered only non-vacuum solutions $(n,p)$ (i.e. $n$, $p>0$) for the 
system (\ref{w-n})--(\ref{w-u}). Due to the nonlinear diffusion $P(\cdot)$ 
(e.g., $P(n)=n^\gamma$), such a system in general is degenerated if there are 
vacuum solutions (i.e. nonnegative solutions). We have shown in the previous 
sections that, under the conditions stated, the system (\ref{w-n})--(\ref{w-u})
does not allow vacuum solutions to develop. In this section, we discuss
how to remove some of the conditions we imposed so to obtain 
parallel results for the case of vacuum solutions.

In the section we use the general Assumption (A) but with nonnegative 
boundary and initial data in (A)(vii)-(viii).

Suppose $P$ is as given by Assumption (A) and $H$ is the associated enthalpy defined in 
(\ref{enthalpy}). For $0<\ep<1$, define
\begin{equation}\label{Pe}
P_\ep(s)=\left\{\begin{array}{ll}
         \frac{P(\ep)}{\ep} s\;\;&\mbox{if }\;0<s\le\ep,\\[5pt]
          P(s) &\mbox{if }\;s>\ep \end{array}\right.
\end{equation}
if $P'(0^+)=0$ (degeneracy). 
When $P'(0^+)>0$, we set $P_\ep=P$ for all $\ep$.
Let $H_\ep$ be the enthalpy associated with $P_\ep$:
\begin{equation}\label{He}
H_\ep(s)=\int_1^s \frac{P_\ep'(s)}{s} ds
   =\left\{\begin{array}{ll}
      \frac{P(\ep)}{\ep}\ln\frac{s}{\ep} + H(\ep)\;\;&\mbox{if}\;\;0<s\le\ep,\\[5pt]
     H(s) &\mbox{if}\;\;s>\ep.\end{array}\right.
\end{equation}
Clearly $\inf_{0 < s\le 1} P_\ep'(s)>0$ and hence $H_\ep(0^+)=\infty$ for each $\ep$. 
We will also perturb the generation-recombination rate $G$ to $G_\ep$ by replacing
the coefficient $Q$ in either (\ref{G0}) or (\ref{G}) with $Q_\ep=Q+\ep$ so that 
$\inf Q_\ep\ge \ep >0$. 
For given nonnegative $(n^0, p^0)$ and $(\bar n,\bar p)$, 
we use positive initial data $(n^0_\ep, p^0_\ep)=(n^0+\ep, p^0+\ep)$
and positive boundary data $({\bar n}_\ep, {\bar p}_\ep)=(\bar n+\ep, \bar p+\ep)$ 
for the perturbed system
\begin{equation}\label{ne}
\langle\frac{\partial n_\ep}{\partial t},\phi\rangle_{V^\ast\times V}
+\mu_n\langle \nabla P_\ep(n_\ep)-n_\ep\nabla u_\ep,\nabla\phi\rangle=\langle G_\ep(n_\ep,p_\ep), \phi\rangle
\end{equation}
\begin{equation}\label{pe}
\langle\frac{\partial p_\ep}{\partial t},\psi\rangle_{V^\ast\times V}
 +\mu_p\langle \nabla P_\ep(p_\ep)+p_\ep\nabla u_\ep,\nabla\psi\rangle= \langle G_\ep(n_\ep,p_\ep),\psi\rangle
\end{equation}
\begin{equation}\label{ue}
\epsilon_0\langle\nabla u_\ep,\nabla\chi\rangle+\langle n_\ep-p_\ep-N, \chi\rangle=0
\end{equation}
for all $(\phi, \psi, \chi)\in V^3$. 
By Theorem \ref{thm:exist}, there exists a solution 
$(n_\ep, p_\ep, u_\ep)\in H^1\cap L^\infty$ such that 
$(n_\ep-{\bar n}_\ep, p_\ep-{\bar p}_\ep, u_\ep-\bar u)\in V^3$ and
\begin{equation}\label{np-bounds}
0< n_\ep(t,x),\;p_\ep(t,x)\le (D_0+1) e^{Bt}
\end{equation}
a.e. in $\Omega$ and $t>0$, where $D_0$ is as in (\ref{DD}) and $B$ in (\ref{B}) 
with $\overline{Q}$ replaced by $\overline{Q}+1$. 
Note also that (\ref{ne})--(\ref{pe}) can be equivalently written as
\begin{eqnarray}\label{ne-i}
\lefteqn{ \langle n_\ep(t)- n^0_\ep, \phi\rangle }\\
&&+\int_0^t\left\{\mu_n \langle \nabla P_\ep(n_\ep(s))-n_\ep(s)\nabla u_\ep(s),\nabla\phi\rangle
- \langle G_\ep(n_\ep(s),p_\ep(s)), \phi\rangle\right\} ds=0 \nonumber
\end{eqnarray}
\begin{eqnarray}\label{pe-i}
\lefteqn{ \langle p_\ep(t)- p^0_\ep, \psi\rangle  }\\
&&+\int_0^t\left\{\mu_p \langle \nabla P_\ep(p_\ep(s))+p_\ep(s)\nabla u_\ep(s),\nabla\psi\rangle
-\langle G_\ep(n_\ep(s),p_\ep(s)), \psi\rangle\right\}ds=0 \nonumber
\end{eqnarray}
for all $\phi$, $\psi\in V$.

From (\ref{np-bounds}), for each $T>0$, 
there exist $n(t,x)$ and $p(t,x)\in L^\infty((0,T)\times\Omega)$ such that
\begin{equation}\label{w*}
(n_\ep, p_\ep)\to (n,p)\;\;\;\mbox{as $\ep\to 0$ in $L^\infty$ weakly-$\ast$}
\end{equation}
(for a sequence of $\{\ep\}$, but still denoted as the same.)
Let $u$ be such that $u-\bar u\in L^2(0,T;V)$ and solve
\begin{equation}\label{uu}
\epsilon_0\langle\nabla u,\nabla\chi\rangle+\langle n-p-N,\chi\rangle=0
\end{equation}
for all $\chi\in V$. We show below that this triplet $(n,p,u)$ is 
a non-negative solution to the original system (\ref{w-n})-(\ref{w-u}). 

From (\ref{ue}) above, we see that $\nabla u_\ep$
is uniformly bounded in $L^2((0,T)\times\Omega)$, and hence, from (\ref{ne-i})--(\ref{pe-i})
and (\ref{np-bounds}) we have
\begin{equation}\label{P-bounds}
\int_0^T|\nabla P_\ep(n_\ep)|_2^2 dt,\;\;\int_0^T|\nabla P_\ep(p_\ep)|_2^2 dt \le \mbox{const.}
\end{equation}
Consequently,
\begin{equation}\label{dnp-bounds}
\int_0^T\left|\frac{\partial n_\ep}{\partial t}\right|_{V^\ast}^2 dt,\;\;\;
\int_0^T\left|\frac{\partial p_\ep}{\partial t}\right|_{V^\ast}^2 dt\le\mbox{const.}
\end{equation}
From (\ref{dnp-bounds}) and (\ref{np-bounds}), we see that $\{n_\ep, p_\ep\}$ is uniformly
bounded in 
$L^2(0,T; L^2(\Omega))\cap H^1(0,T; V^\ast)$, hence compact in $L^2(0,T;V^\ast)$.
Therefore there exist a subsequence of $\{\ep\}$ (still denoted by $\{\ep\}$) and 
a set $E\subset [0,T]$ with $m E=0$ such that 
$$(n_\ep(t), p_\ep(t))\to (n(t), p(t))\;\;\;\mbox{in}\;\;V^\ast\;\;\;
\mbox{for all}\;\;t\in [0,T]\setminus E$$
as $\ep\to 0$. From (\ref{ue}), 
$$u_\ep=\tilde u+\frac{1}{\epsilon_0}(-\Delta_0)^{-1}(n_\ep-p_\ep-N)$$
where $\tilde u\in\bar u+ V$ satisfies $\Delta\tilde u=0$.
Since $(-\Delta_0)^{-1}$ is an isometry from $V^\ast$ to $V$, we have
$$u_\ep(t) \to u(t)\;\;\;\mbox{in }\;H^1(\Omega)\;\;\;\mbox{a.e.}\;\;t\in (0,T)$$
and hence, together with (\ref{w*}),
$$ n_\ep\nabla u_\ep \to n\nabla u
\quad\mbox{and}\quad p_\ep\nabla u_\ep \to p\nabla u\;\;\;\mbox{weakly in}\;\;L^2(\Omega)
\;\;\mbox{a.e.}\;\;t\in (0,T). $$

From (\ref{P-bounds}), there exists $\eta\in L^2(0,T; V)$ such that
$$P_\ep(n_\ep)-P_\ep(\bar n)\to \eta\;\;\;\mbox{weakly in}\;\;
L^2(0,T; V)$$
(a further subsequence if necessary.) Since $\frac{\partial P_\ep(n_\ep)}{\partial t}
=P_\ep'(n_\ep)\frac{\partial n_\ep}{\partial t}$ and $P_\ep'(s)\le \sup P'(s)$
on bounded intervals, we see from (\ref{P-bounds})-(\ref{dnp-bounds}) that
$$\{P_\ep(n_\ep)-P_\ep(\bar n)\}\;\;
\mbox{is uniformly bounded in }\;\; 
L^2(0,T; V)\cap H^1(0,T; V^\ast).$$
Hence by Aubin's Lemma, we can also assume
$$P_\ep(n_\ep)-P_\ep(\bar n)\to \eta
\;\;\;\mbox{in}\;\;L^2((0,T)\times\Omega)$$
(a further subsequence if necessary.) Note that 
$$|P_\ep(\phi)-P(\phi)|_2^2
= \int_0^T\int_{\phi<\ep}(P_\ep(\phi)-P(\phi))^2dxdt 
\le 2P(\ep)^2 T|\Omega|\to 0$$
for any $\phi\in L^\infty((0,T)\times\Omega)$. Hence we have
$$P(n_\ep)-P(\bar n)\to\eta\;\;\;\mbox{in}\;\;L^2((0,T)\times\Omega).$$
Thus
$$P(n_\ep)\to P(\bar n)+\eta\;\;\;\mbox{or}\;\;\;
n_\ep\to P^{-1}(P(\bar n)+\eta)\;\;\;
\mbox{a.e. in }\;(0,T)\times\Omega.$$
Together with (\ref{w*}) we see $n=P^{-1}(P(\bar n)+\eta)$ and
$n_\ep\to n$ a.e. in $(0, T)\times\Omega$.
The same can be shown for $p_\ep$:
$$P_\ep(p_\ep)\to P(p)\;\;\;\mbox{weakly in $L^2(0,T;H^1(\Omega))$
and }\;\; p_\ep\to p\;\;\;\mbox{a.e. in}\;\;(0,T)\times\Omega.$$
Therefore, $G_\ep(n_\ep, p_\ep)\to G(n,p)$ a.e. in $(0,T)\times\Omega$.
Passing to the limit of $\ep\to 0$ in (\ref{ne-i})--(\ref{pe-i}) we obtain
\begin{eqnarray*}
\lefteqn{ \langle n(t)- n^0,\phi\rangle  }\\
&&+\int_0^t\left\{\mu_n \langle \nabla P(n(s))-n(s)\nabla u(s),\nabla\phi\rangle
- \langle G(n(s),p(s)), \phi\rangle\right\} ds=0,
\end{eqnarray*}
\begin{eqnarray*}
\lefteqn{ \langle p(t)- p^0,\psi\rangle   }\\
&&+\int_0^t\left\{\mu_p \langle \nabla P(p(s))+p(s)\nabla u(s),\nabla\psi\rangle
-\langle G(n(s),p(s)), \psi\rangle\right\}ds=0
\end{eqnarray*}
and $u$ as in (\ref{uu}). That is, $(n,p,u)$ is a solution that
satisfies (\ref{w-n})--(\ref{w-u}) with (\ref{u}) in the sense that
\begin{equation}\label{vac}
\begin{array}{l}
(n,p, u)\in L^\infty((0,T)\times\Omega)^3\cap H^1(0,T; V^\ast)^3\;\;\;
\mbox{with}\;\;\;n, \;p \ge 0, \\[5pt]
(P(n), P(p), u)\in (P(\bar n), P(\bar p),\bar u) + L^2(0,T; V)^3.
\end{array}
\end{equation}

Although we cannot obtain $H^1$-regularity for $(n,p)$ in general, 
it can be shown that 
\begin{equation}\label{almost}
\mbox{for every}\;\delta>0:\;\;(n-\delta)^+,\;(p-\delta)^+\in H^1(\Omega). 
\end{equation}
Indeed, since $P(n)\in H^1(\Omega)$, so is $(P(n)-P(\delta))^+$, and
$$|\nabla (P(n)-P(\delta))^+|_2^2=
\int_{n>\delta} |P'(n)\nabla n|^2 dx
\ge\kappa_0\int_{n>\delta}|\nabla n|^2 dx
=\kappa_0|\nabla (n-\delta)^+|_2^2$$
where $\kappa_0=\inf_{\delta\le s\le |n|_\infty}P'(s)>0$.
Thus, under Assumption (A2), the argument in \S\ref{sec:absorb} for 
the upper bound can be applied to obtain the absorbing $L^\infty$ upper bound 
$D(t)$ for $n$ and $p$:
$$0\le n(t,x),\;p(t,x)\le d(1+d_0 e^{-\lambda t})\;\;\;\mbox{a.e. }\;t>0,\;x\in\Omega$$
where the constants $d$, $d_0$ and $\lambda$ are as in \S\ref{sec:absorb},
and $d$ is independent of the initial data. 
The set of admissible initial data is now extended to
$$\tilde{H}_{\rm ad}
=\{(n, p)\in L^\infty(\Omega)^2:\;\;n(x)\ge 0,\;p(x)\ge 0\;\;\mbox{a.e.}\;\}.$$

Our uniqueness argument in \S\ref{sec:unique} requires positive lower bounds of
solutions, and this restriction cannot be removed by the above perturbation technique.
Hence we do not have uniqueness result for vacuum solutions, and thus there is no
result on global attractors.

As for the stationary system, the existence of positive solutions $\{n_\ep, p_\ep\}$
to the perturbed systems is by Theorem \ref{thm:stat} if $h_\infty>d_0$ ($d_0$ given by
(\ref{d0})). The convergence argument of
$\{n_\ep, p_\ep\}$ to a non-negative solution $(n,p)$ 
is similar (but simpler), and we skip the details.

These parallel results can be summarized in the following theorem.

\begin{theorem}\label{thm:vac}
{\rm (Vacuum solutions.)} Suppose Assumption (A) holds and $G$ is
given by either (\ref{G0}) or (\ref{G}). 
\begin{itemize}\itemsep=-2pt
\item[{\rm (i)}]{\rm (Global solutions.)} For each $T>0$, there exists
at least one solution $(n, p, u)$ to (\ref{w-n})--(\ref{w-u}) 
with the regularities stated in (\ref{vac}). 
\item[{\rm (ii)}] {\rm (Absorbing set.)} If Assumption (A2) also holds, then any global solution
$(n,p,u)$ from the initial data $(n^0, p^0)\in\tilde{H}_{\rm ad}$ satisfies
$$0\le n(t,x),\;p(t,x)\le d(1+d_0 e^{-\lambda t})$$
where $d$, $d_0$ and $\lambda$ are as in \S\ref{sec:absorb},
and $d$ is independent of the initial data $(n^0, p^0)$. Therefore,
the dynamical system has an absorbing set
$${\cal B}=\{(n, p)\in\tilde{H}_{\rm ad}:\;0\le n(x), p(x)\le 2d\}.$$
\item[{\rm (iii)}] {\rm (Stationary solutions.)} If $h_\infty>d_0$ ($d_0$ as in 
(\ref{d0})), then the stationary system (\ref{n0})--(\ref{u0}) has
at least a solution $(n, p, u)$ such that
$(n, p, u)\in L^\infty(\Omega)^3$ with $n(x)$, $p(x)\ge 0$ and 
$(P(n), P(p), u)\in (P(\bar n), P(\bar p),\bar u) + V^3$.
\end{itemize}
\end{theorem}


%%%%%%%%%%%%%%%%%%%%

\renewcommand{\thesection}{\Alph{section}}\nopagebreak
\setcounter{section}{0}

\section{Appendix}
\setcounter{equation}{0}
%%\setcounter{theorem}{0}

In this appendix we give the proofs of the lemmas we have used in this
paper for establishing the $L^\infty$ estimates for solutions of partial
differential equations. This technique is an extension of the classical 
maximum principle for parabolic and elliptic equations.

\begin{lemma}\label{lm:Linfty}    
Suppose $\phi\in H^1(\Omega)$, and $\phi_\tau=(\phi-\tau)^+\in V$ for 
$\tau\ge\bar\tau$. If  $\phi_\tau$ satisfies 
$$\langle \nabla\phi_\tau,\nabla\phi_\tau\rangle\le \langle f,\phi_\tau\rangle$$
for some $f\in L^2(\Omega)$. Then     
$$\phi(x)\le \bar\tau +4\alpha_6 |\Omega|^{1/6} |f|_2 .$$    
\end{lemma}    
{\it Proof.} Let $\Omega_\tau=\{ x\in\Omega:\; \phi(x)>\tau\}$. 
By the Poincar\'{e} inequality,
$$\alpha_6^{-2}|\phi_\tau|_6^2\le |\nabla \phi_\tau|_2^2
\le \langle f,\phi_\tau\rangle\le |f|_2 |\phi_\tau|_6 |\Omega_\tau|^{1/3},$$
that is,
$$|\phi_\tau|_6\le\alpha_6^2|f|_2|\Omega_\tau|^{1/3}.$$
For $\hat\tau>\tau\ge \bar\tau$,
$$|\phi_\tau|_6^6\ge\int_{\Omega_{\hat\tau}}|\phi_\tau|^6 dx
\ge (\hat\tau-\tau)^6 |\Omega_{\hat\tau}|.$$
Therefore
$$|\Omega_{\hat\tau}|
\le (\alpha_6^2 |f|_2)^6 (\hat\tau-\tau)^{-6}|\Omega_\tau|^2.$$
Then by applying Lemma \ref{lm:stamp} below to $|\Omega_\tau|$, 
a non-increasing function in $\tau$, we obtain
$$|\Omega_\tau|=0\;\;\;\mbox{for}\;\;\; \tau\ge \tau^\ast$$
where
$$\tau^\ast=\bar\tau+4\alpha_6^2|\Omega|^{1/6} |f|_2. $$
Therefore we have $\phi_{\tau^\ast}=0$ a.e. 
in $\Omega$. Hence $\phi\le \tau^\ast$ as stated.  $\Box$

For the sake of completeness, we now state and prove the elementary lemma used 
above. See \cite[Lemma 2.9]{Tr}. 

\begin{lemma}\label{lm:stamp}
Suppose $\varphi(\tau)$ is a nonnegative, non-increasing function on 
$[\tau_0,\tau_1]$, and there are positive constants $M$, $\gamma$ and $\delta$
such that
\begin{equation}\label{stamp}
\varphi(\hat\tau)\le M(\hat\tau-\tau)^{-\gamma}\varphi(\tau)^{1+\delta}
\;\;\;\mbox{for all}\;\;\hat\tau, \tau\in [\tau_0,\tau_1]\;\;\;
\mbox{with}\;\;\hat\tau> \tau.
\end{equation}
Then
$$\varphi(\tau^*)=0\;\;\;\mbox{provided that}\;\;\;
\tau^*\equiv \tau_0+M^{1/\gamma} 2^{(1+\delta)/\delta}
\varphi(\tau_0)^{\delta/\gamma}\le \tau_1.$$
\end{lemma}
{\em Proof.} Let 
$\sigma=M^{1/\gamma} 2^{(1+\delta)/\delta}\varphi(\tau_0)^{\delta/\gamma}$
and $x_k=\tau_0+\sigma (1- 2^{-k})$ for $k=0,1,\cdots$. 
By induction and (\ref{stamp}), it is straightforward to show that
$$\varphi(x_k)\le\varphi(x_0) 2^{-k\gamma/\delta}\;\;\;k=0,1,\cdots.$$
Therefore, when $\tau^*=\tau_0+\sigma\le \tau_1$, 
$\varphi(\tau^*)\le\varphi(x_k)\le\varphi(x_0) 2^{-k\gamma/\delta}\to 0$
as $k\to\infty$.  $\Box$

We have also used the following extension to the above lemma in \S\ref{sec:absorb} 
(see (\ref{d})). Other variations can be found in \cite{Lady,JDE2}.

\begin{lemma}\label{lm:stamp2}
Suppose $\varphi(\tau)$ is a nonnegative, non-increasing function on 
$[\tau_0,\infty)$, and there are positive constants $M$, $\gamma$ and $\delta$
such that 
\begin{equation}\label{stamp2}
\varphi(\hat\tau)\le 
M\tau^\gamma(\hat\tau-\tau)^{-\gamma}\varphi(\tau)^{1+\delta}
\;\;\;\mbox{for all}\;\;\hat\tau>\tau\ge\tau_0
\end{equation}
Then
$$\varphi(\tau^*)=0\;\;\;\mbox{for}\;\;\;
\tau^*\equiv 2\tau_0(1+2^{(1+2\delta)/\delta^2}M^{(1+\delta)/\delta\gamma}
\varphi(\tau_0)^{(1+\delta)/\gamma}).$$
\end{lemma}
{\em Proof.} For any $m>1$, from (\ref{stamp2}) with $\tau=\tau_0$ and 
$\hat\tau=m\tau_0$ we have
$$\varphi(m\tau_0)\le M (m-1)^{-\gamma}\varphi(\tau_0)^{1+\delta}.$$
Now apply Lemma \ref{lm:stamp} to $\varphi$ on the interval 
$[m\tau_0, 2m\tau_0]$ to obtain
$$\varphi(\tau^\ast)=0\;\;\;\mbox{for}\;\;
\tau^\ast\ge\tau^\ast_0\equiv m\tau_0
+M^{1/\gamma} 2m\tau_0 2^{(1+\delta)/\delta}\varphi(m\tau_0)^{\delta/\gamma}$$
provided that $\tau^\ast\le 2m\tau_0$. From the above estimate 
for $\varphi(m\tau_0)$,
$$\frac{\tau^\ast_0}{m\tau_0}
\le 1+M^{(1+\delta)/\gamma} 2^{(1+2\delta)/\delta}
\varphi(\tau_0)^{(1+\delta)\delta/\gamma}(m-1)^{-\delta}=2$$
if we choose
$$m=1+M^{(1+\delta)/\delta\gamma} 2^{(1+2\delta)/\delta^2}
\varphi(\tau_0)^{(1+\delta)/\gamma}.$$
Hence for this chosen $m$ we set $\tau^\ast=2m\tau_0\ge\tau^\ast_0$ and thus
$\varphi(\tau^\ast)=0$. $\Box$



%%%%%%%%%%%%%%%%%
\begin{thebibliography}{99}    
\bibitem{Bz} H.~Brezis, Monotonicity methods in Hilbert spaces and some
  applications to nonlinear partial differential equations, in ``Contributions
  to Nonlinear Functional Analysis,'' E.~Zarantonello ed., Academic Press, 
 (1979), pp.~101--156.  
\bibitem{B} V.~Barcilon, D.-P.~Chen, and R.~S.~Eisenberg, Ion flow through narrow membrane channels: Part II, SIAM J. Appl. Math., 52(1992), pp.~1405--1425.
\bibitem{BFI} S.~Busenberg, W.~Fang, and K.~Ito, Modeling and analysis    
  of laser-beam-induced current images in semiconductors, SIAM J. Appl. Math.,  
  53(1993), pp.~187--204. 
\bibitem{CL} Y.~S.~Choi and R. Lui, Multi-dimensional electro-chemistry model,
     Arch. Rational Mech. Anal.,  130(1995), pp.~315--342.    
\bibitem{JDE1} W.~Fang and K.~Ito, On the time-dependent     
 drift-diffusion model for semiconductors,  J. Differential Equations,     
 117(1995), pp.~245--280.    
\bibitem{JDE2} W.~Fang and K.~Ito,  Global solutions of     
 the time-dependent drift-diffusion semiconductor equations,     
  J. Differential Equations, 123(1995), pp.~523--566. 
\bibitem{JDE3} W.~Fang and K.~Ito,  Asymptotic behavior of     
 the drift-diffusion semiconductor equations,     
  J. Differential Equations, 123(1995), pp.~567--587.
\bibitem{GG} H.~Gajewski and K.~Gr\"{o}ger, On the basic equations for    
  carrier transport in semiconductors,  J. Math. Anal. Appl.,     
  113(1986), pp.~12--35. 
\bibitem{GG2} H.~Gajewski and K.~Gr\"{o}ger, Semiconductor equations for variable mobilities
  based on Boltzmann statistics or Fermi-Dirac statistics, Math. Nachr., 140(1989), pp.~7--36.   
\bibitem{J} J.~W.~Jerome, Consistency of semiconductor modeling:     
 On existence/stability analysis for the stationary Van Roosbroeck system,    
 SIAM J. Appl. Math., 45(1985), pp.~565--590.     
%%\bibitem{JS} J.~W.~Jerome and C.-W.~Shu, Energy models for one-carrier     
%% transport in semiconductor devices,  in  ``Semiconductors II,'' 
%% W.~M.~Coughhran~Jr., J.~Cole, P.~Lloyd, and J.~K.~White eds,  IMA Vol.~59, 
%% pp.~185--207,  Springer-Verlag,  New York, 1994.
\bibitem{J0} J.~W.~Jerome, {\it Analysis of Charge Transport}, Springer-Verlag, 
 Berlin-Heidelberg, 1996.
\bibitem{Ju} A.~J\"{u}ngel, A nonlinear drift-diffusion system with electric convection
 arising in electrophoretic and semiconductor modeling, Math. Nachr., 185(1997), pp.~85--110.
\bibitem{Lady} O.~A.~Ladyzenskaja, V.~A.~Solonnikov, and N.~N.~Uralceva,
 ``Linear and Quasilinear Equations of Parabolic Type,'' 
 American Mathematical Society, Providence, RI, 1968.
\bibitem{MN} P.~Marcati and R.~Natalini, Weak solutions to a hydrodynamic model
  for semiconductors and relaxation to the drift-diffusion equation,
 Arch. Rational Mech. Anal.,  129(1995), pp.~129--145.
%\bibitem{M} P.~A.~Markowich, ``The Stationary Semiconductor Device     
%     Equations,'' Springer-Verlag, Wein-New York, 1986.      
\bibitem{MRS} P.~A.~Markowich, C.~A.~Ringhofer, and C.~Schmeiser,    
  {\it Semiconductor Equations,} Springer-Verlag, Wein-New York, 1990.
\bibitem{MU} P.~A.~Markowich and A.~Unterreiter, Vacuum solutions of a stationary drift-diffusion
 model, Ann. Sc. Norm. Sup. di Pisa, 20(1993), pp.~371--386. 
\bibitem{PJ} J.-H. Park and J. W. Jerome,  Qualitative properties of steady-state
 PNP systems: mathematical study,  SIAM J. Appl. Math., 57(1997),  pp.~609--630. 
\bibitem{R} I.~Rubinstein, ``Electro-Diffusion of Ions,'' SIAM Studies in Applied Math., SIAM,
 Philadelphia, PA, 1990.
\bibitem{RO} M.~Rudan and F.~Odeh, Multi-dimensional discretization     
 scheme for the hydrodynamic model of semiconductor devices,      
 COMPEL, 5(1986), pp.~149--183.    
\bibitem{ST} T.~I.~Seidman and G.~M.~Troianiello, Time-dependent solutions of 
a nonlinear system arising in semiconductor theory, Nonlinear Anal., 9(1985), pp.~1137--1157.   
\bibitem{Tm} R.~Temam, ``Infinite-Dimensional Dynamical Systems in Mechanics
  and Physics,'' Springer-Verlag, New York, 1987.
\bibitem{Tr} G.~M.~Troianiello, ``Elliptic Differential Equations and   
 Obstacle Problems,'' Plenum Press, New York, 1987.    
\end{thebibliography} \bigskip

\noindent{\sc Weifu Fang}\\
Department of Mathematics, West Virginia University\\
Morgantown, WV 26506. USA \\
Email address: wfang@math.wvu.edu \medskip
  
\noindent{\sc Kazufumi Ito}\\
Department of Mathematics, North Carolina State University\\
Raleigh, NC 27695. USA \\
Email address: kito@eos.ncsu.edu 



\subsection*{Addendum: June 21, 1999.}

 
On page 33, there is an error in the argument of
the strong convergence of $P(n_\ep)\to P(\bar n)+\eta$
in $L^2((0,T)\times\Omega)$. A correct proof of the
same result is given below.

First we note that, on the top of page 33, the convergence of
$(n_\ep, p_\ep)\to (n,p)$ is strong in $L^2(0,T;V^\ast)$
and weak in $L^2(0,T;H)$.
Then the paragraph from line 13 to line 26 on page 33 should
be replaced by the following:

From (\ref{P-bounds}), there exists $\eta\in L^2(0,T; H^1(\Omega))$ 
such that 
$$P_\ep(n_\ep)\to \eta\;\;\;\mbox{weakly in}\;\;L^2(0,T;H^1(\Omega))$$
(a further subsequence if necessary.) 
Since $n_\ep$, $p_\ep$ are in $L^\infty((0,T)\times\Omega)$, 
there exists $G_0$ such that 
$$G_\ep(n_\ep,p_\ep)\to G_0 \;\;\;\mbox{weakly in}\;\;L^2((0,T)\times\Omega).$$
Hence, passing to the limit in (\ref{ne-i}) we obtain
$$\langle n(t)- n^0, \phi\rangle
+\int_0^t\mu_n\langle \nabla\eta,\nabla\phi\rangle ds
=\int_0^t\left\{\mu_n\langle n\nabla u,\nabla\phi\rangle
+\langle G_0, \phi\rangle\right\} ds.$$
By subtracting this equation from (\ref{ne-i}) and 
setting $\phi=(-\Delta_0)^{-1}(n_\ep-n)$, we can obtain
$$\frac{1}{2}|n_\ep(t)-n(t)|_{V^\ast}^2
 +\mu_n\int_0^t \langle P_\ep(n_\ep)-\eta, n_\ep-n\rangle ds 
=O(|n_\ep-n|_{L^2(0,T;V^\ast)}).$$
Note that $n_\ep\to n$ strongly in $L^2(0,T;V^\ast)$ and, for any 
$\phi\in L^\infty((0,T)\times\Omega)$, there holds
$$|P_\ep(\phi)-P(\phi)|_2^2
= \int_0^T\int_{\phi<\ep}(P_\ep(\phi)-P(\phi))^2dxdt 
\le 2P(\ep)^2 T|\Omega|\to 0.$$
Hence we have
$$\int_0^T\langle P(n_\ep)-\eta, n_\ep-n\rangle ds\to 0
  \eqno{\mbox{(B.1)}}$$
which is equivalent to
$$\int_0^T\langle P(n_\ep), n_\ep\rangle
\to\int_0^T\langle \eta, n\rangle ds.$$
Since $P$ is maximum monotone on $L^2(0,T;H)$, we can see that 
the above gives $P(n)=\eta$. 
Therefore the limit in (B.1) leads to $P(n_\ep)\to P(n)$ and
$n_\ep\to n$ strongly in $L^2((0,T)\times\Omega)$, and thus 
$n_\ep\to n$ a.e. in $(0,T)\times\Omega$. \hfill$\diamondsuit$
 

\end{document}