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

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2012 (2012), No. 169, pp. 1--18.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2012 Texas State University - San Marcos.}
\vspace{9mm}}

\begin{document}
\title[\hfilneg EJDE-2012/169\hfil Existence of solutions]
{Existence of solutions to quasilinear elliptic systems
 with combined critical Sobolev-Hardy terms}

\author[N. Nyamoradi, M. Javidi  \hfil EJDE-2012/169\hfilneg]
{Nemat Nyamoradi, Mohamad Javidi}  

\address{Nemat Nyamoradi \newline
Department of Mathematics, Faculty of Sciences\\
Razi University, 67149 Kermanshah, Iran}
\email{nyamoradi@razi.ac.ir, neamat80@yahoo.com}

\address{Mohamad Javidi \newline
Department of Mathematics, Faculty of Sciences\\
Razi University, 67149 Kermanshah, Iran}
\email{mo\_javidi@yahoo.com}

\thanks{Submitted May 17, 2012. Published October 4, 2012}
\subjclass[2000]{35B33, 35J60, 35J65}
\keywords{Ekeland variational principle;
 critical Hardy-Sobolev exponent;  \hfill\break\indent
concentration-compactness principle}

\begin{abstract}
 This article is devoted to the study of multiple positive solutions
 to a singular elliptic system where the nonlinearity involves a
 combination of concave and convex terms. Using the effect  of the
 coefficient of the critical nonlinearity, and a variational
 method, we establish the main result which is based on a compactness
 argument.
\end{abstract}

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

%\newtheorem{remark}[theorem]{Remark}

\allowdisplaybreaks 

 \section{Introduction}

The aim of this paper is to establish the existence
 of nontrivial solution to the  elliptic system
\begin{equation}\label{1}
      \begin{gathered}
    - \Delta_p u  - \mu \frac{|u|^{p - 2} u}{|x|^p} =
    \frac{|u|^{p^* (s_1) - 2} u}{|x|^{s_1}} +   \frac{\alpha}{\alpha
  + \beta} Q (x ) \frac{ |u|^{\alpha - 2} |v|^\beta u}{|x - x_0|^t}
     + \lambda h (x) \frac{ |u|^{q - 2} u}{|x|^s}, \\
%    \quad x \in   \Omega,\\
   - \Delta_p v  - \mu \frac{|v|^{p - 2} v}{|x|^p} =
    \frac{|v|^{p^* (s_2) - 2} v}{|x|^{s_2}} +   \frac{\beta}{\alpha + \beta} Q (x )
 \frac{ |u|^\alpha  |v|^{\beta -2} v}{|x- x_0|^t} +
    \lambda h (x) \frac{ |v|^{q - 2} v}{|x|^s},\\
     x \in   \Omega,\\
        u = v = 0,
      \quad x \in  \partial \Omega
        \end{gathered}
\end{equation}
where $\Delta_p u = \operatorname{div} (|\nabla u|^{p - 2} \nabla u)$,
$0 \in \Omega$ is a bounded domain in $\mathbb{R}^N$ ($N \geq 3$)
with smooth boundary $\partial \Omega$, $\lambda > 0$ is a
parameter, $1 \leq q < p$, $1 < p < N$, $0 \leq \mu <
\overline{\mu} \triangleq \big(\frac{N - p}{p} \big)^p$; $Q (x)$
is nonnegative and continuous on $\overline{\Omega}$ satisfying
some additional conditions which will be given later,
$Q (x_0) =
\|Q\|_\infty$ for $0 \ne x_0 \ne \Omega$, $h (x) \in C (\overline{\Omega})$;
 $\alpha, \beta > 1$,
$\alpha + \beta = p^*(t) \triangleq \frac{p (N - t)}{N - p}$,
$p^*(s) \triangleq \frac{p (N - s)}{N - p}$
($0 < s, s_1, s_2 \leq t < p $) are critical Sobolev-Hardy exponents.
 Note that $ p^*(0) = p^* : =  \frac{N p}{N - p}$ is the critical Sobolev exponent.

We denote by $W_0^{1,p} (\Omega)$ the completion of
 $C_0^\infty (\Omega)$ with respect to the norm
 $\big(\int_\Omega |\nabla \cdot|^p dx \big)^{1/p}$.

Problem \eqref{1} is related to the well known Caffarelli-Kohn-Nirenberg 
inequality in \cite{L. Caffarelli}:
\begin{equation} \label{ine1}
 \Big(\int_\Omega \frac{|u|^r}{|x|^t } dx \Big)^{p/r}
\leq C_{r,t, p} \int_\Omega  |\nabla  u|^p
 dx, \quad \text{for all } u \in W_0^{1,p} (\Omega),
\end{equation}
where $p \leq r < p^*(t)$. When $t = r = p$, the above inequality
becomes the well known Hardy inequality \cite{L. Caffarelli,J.
Garcia Azorero,N. Ghoussoub}:
\begin{equation} \label{ine11}
 \int_\Omega \frac{|u|^p}{|x|^p } dx
\leq \frac{1}{\overline{\mu}} \int_\Omega  |\nabla  u|^p
 dx, \quad \text{for all }  u \in W_0^{1,p} (\Omega).
\end{equation}
In the space $W_0^{1,p} (\Omega)$ we use the  norm
\begin{align*}
 \|u\|_\mu = \|u\|_{D^{1,p} (\Omega)} : = \Big(\int_\Omega \Big(|\nabla
u|^p - \mu \frac{|u|^p}{|x|^p } \Big) dx \Big)^{1/p},
\quad \mu \in [0, \overline{\mu}).
\end{align*}
By using the Hardy inequality \eqref{ine11} this norm is
equivalent to the usual norm $\big(\int_\Omega |\nabla u|^p dx
\big)^{1/p}$. The elliptic operator $L := \big(|\nabla
\cdot|^{p - 2} \nabla \cdot - \mu \frac{|\cdot|^{p - 2}}{|x|^p }
\big)$ is positive in $W_0^{1,p} (\Omega)$ if $0 \leq \mu <
\overline{\mu}$.

Now, we define the space $W =  W_0^{1,p}(\Omega) \times
W_0^{1,p}(\Omega)$ with the norm
\[
\|(u, v)\|^p = \|u\|_\mu^p + \|v\|_\mu^p.
\]
Also, by Hardy inequality and Hardy-Sobolev inequality, for
 $0 \leq \mu < \overline{\mu}$, $0 \leq t < p$ and
$p \leq r \leq p^*(t)$ we can define the best Hardy-Sobolev constant:
\begin{equation} \label{hardy}
 A_{\mu, t, r} (\Omega) = \inf_{u \in W_0^{1,p}(\Omega) \setminus
\{0\}} \frac{\int_\Omega  \Big(|\nabla u|^p - \mu
\frac{|u|^p}{|x|^p } \Big) dx}{\Big(\int_\Omega
\frac{|u|^r}{|x|^t} dx \Big)^{p/r}}.
\end{equation}
In the important case when $r = p^*(t)$, we simply denote $A_{\mu,
t, p^*(t)}$ as $A_{\mu, t}$.

For any $0 \leq \mu < \overline{\mu}$, $\alpha, \beta > 1$ and
$\alpha + \beta = p^* (t)$, by \eqref{ine1}, \eqref{ine11}, $0 <
s_1, s_2 \leq t < p$, Set
\begin{gather}
A_{\mu, s} := \inf_{u \in W_0^{1,p} (\Omega) \setminus \{0\}}
\frac{\int_\Omega \big(|\nabla u|^p - \mu \frac{|u|^p}{|x|^p }
\big) dx}{\big(\int_\Omega \frac{|u|^{p^* (s)}}{|x|^s} dx
\big)^{\frac{p}{p^* (s)}}}\,, \label{san}
\\
S_{s,\alpha, \beta} := \inf_{(u, v) \in W \setminus \{(0, 0)\}}
\frac{\int_\Omega \big(|\nabla u|^p + |\nabla v|^p - \mu
\frac{|u|^p  + |v|^p}{|x|^p } \big) dx}{\big(\int_\Omega
\frac{|u|^\alpha |v|^\beta}{|x|^s} dx \big)^{\frac{p}{\alpha +
\beta}}}\,. \label{hardy2}
\end{gather}
Then we have the following equality (whose proof is the same as that of Theorem 5 in
\cite{C. Alves})
\begin{align*}
S_{s,\alpha, \beta} (\mu) = \Big(\big (\frac{\alpha}{\beta}
\big)^\frac{\beta}{\alpha + \beta} + \big (\frac{\beta}{\alpha}
\big)^\frac{\alpha}{\alpha + \beta} \Big)A_{\mu, s}\,.
\end{align*}

Throughout this paper, let $R_0$ be the positive constant such
that $\Omega \subset B(0; R_0)$, where 
$B(0; R_0) = \{x \in \mathbb{R}^N : |x| < R_0 \}$.
 By Holder and Sobolev-Hardy inequalities, for all $u \in W^{1,p}_ 0 (\Omega)$,
 we obtain
\begin{equation} \label{hold}
\begin{split}
 \int_\Omega  \frac{ |u|^q}{|x|^s}
& \leq  \Big(\int_{B(0; R_0)} |x|^{-s}
 \Big)^\frac{p^* (s) - q}{p^* (s)} \Big(\int_\Omega \frac{ |u|^{p^*(s)}}{|x|^s}
 \Big)^\frac{q}{p^* (s)} \\
&\leq  \Big(\int_0^ {R_0} r^{N -s + 1} dr
 \Big)^\frac{p^* (s) - q}{p^* (s)} A_{\mu, s}^{-\frac{q}{p}}
 \|u\|^q \\
&\leq  \Big(\frac{N \omega_N R_0^{N - s}}{N - s}
 \Big)^\frac{p^* (s) - q}{p^* (s)} A_{\mu, s}^{-\frac{q}{p}}
 \|u\|^q,
\end{split}
\end{equation}
where
$\omega_N = \frac{2 \pi^{N/2}}{N \Gamma(N/2)}$
is the volume of the unit ball in $\mathbb{R}^N$ .

Existence of nontrivial non-negative solutions for elliptic
equations with singular potentials were recently studied by
several authors, but, essentially, only with a solely critical
exponent. We refer, e.g., in bounded domains and for $p = 2$ to
\cite{D. Cao, N. Ghoussoub,N. Ghoussoub1, Han1, Han2, Han4, Han3}, 
and for general $p > 1$ to
\cite{ D. Cao1, Degiovanni, F. Demengel, R. Filippucci, N. Ghoussoub2,
 P. Han, Kang, S. Peng, Xuan1} 
and the references therein. For
example, Han and Liu \cite{Han3} studied the  problem
\begin{equation}\label{ne2}
      \begin{gathered}
    - \Delta u  - \mu \frac{ u}{|x|^2} =
    \lambda u +    Q (x ) |u|^{2^* - 2} u,    \quad x \in   \Omega,\\
        u =  0, \quad x \in  \partial \Omega
        \end{gathered}
\end{equation}
where  $0 \in \Omega$ is a smooth bounded domain in $\mathbb{R}^N$
($N \geq 5$), $\lambda > 0$, $0 \leq \mu < \overline{\mu}
\triangleq \big(\frac{N - 2}{2} \big)^2$,
 $2^* = \frac{2 N}{N -2}$ and $Q (x)$ is nonnegative and continuous on
$\overline{\Omega}$ satisfying some suitable conditions. using
critical point theory, the authors proved the existence of
nontrivial solutions to problem \eqref{ne2}. Also, by
investigating the effect of the coefficient $Q$, Han \cite{Han2}
studied problem \eqref{ne2} and proved that there exists
$\lambda_0 > 0$ such that \eqref{ne2} has at least $k$ positive 
solutions for $\lambda \in (0, \lambda_0)$.

Kang in \cite{Kang} studied the following elliptic equation via
the generalized Mountain-Pass theorem \cite{P. Rabinowitz},
\begin{equation}\label{ne1}
      \begin{gathered}
    - \Delta_p u  - \mu \frac{|u|^{p - 2} u}{|x|^p} =
    \frac{|u|^{p^* (t) - 2} u}{|x|^t} +   \lambda \frac{ |u|^{p - 2} u}{|x|^s},
    \quad x \in \Omega,\\
           u =  0,      \quad x \in  \partial \Omega
        \end{gathered}
\end{equation}
where $\Omega \subset \mathbb{R}^N$ is a bounded domain, 
$1 < p < N$, $0 \leq  s, t < p$ and $0 \leq \mu < \overline{\mu} \triangleq
\big(\frac{N - p}{p} \big)^p$. Degiovanni and Lancelotti
\cite{Degiovanni} studied problem \eqref{ne1} with $\mu = s = t =
0$ and proved that \eqref{ne1} has at least one positive solutions
for $\lambda \geq \lambda_1 : = A_{0,0}$ ($A_{0,0}$ is defined in
\eqref{san}). Indeed, in \cite{Degiovanni} the much more difficult
case $\lambda \geq \lambda_1$ is treated.

The authors  in \cite{R. Filippucci}, via the Mountain-Pass
Theorem of Ambrosetti and Rabinowitz \cite{A. Ambrosetti}, proved
that
\begin{equation*}
        - \Delta_p u  - \mu \frac{u^{p - 1} }{|x|^p} = |u|^{p^* - 1} +
    \frac{u^{p^* (s) - 1}}{|x|^s} ,
   \quad  \text{in  } \mathbb{R}^N,
\end{equation*}
 admits a positive solution in $\mathbb{R}^N$,
whenever $\mu < \overline{\mu} \triangleq \big(\frac{N - p}{p}
\big)^p$ and $0 < s < p$.

Recently, in \cite{Xuan1} the author studied the following
equation via the Mountain-Pass theorem,
\begin{equation*}
        - \operatorname{div} \Big(\frac{|D u|^{p - 2} D u}{|x|^{a p}} \Big)  - \mu \frac{|u|^{p - 2} u}{|x|^{(a + 1)p}} = \frac{|u|^{p^* (b) - 2} u}{|x|^{b p^*}} +
    \frac{|u|^{p^* (c) - 2} u}{|x|^{c p^*}},
   \quad  \text{in  } \mathbb{R}^N
\end{equation*}
where $1 < p < N$, $0 \leq \mu < \overline{\mu} \triangleq
\big(\frac{N - (a + 1) p}{p} \big)^p$, $0 \leq a < \frac{N -
p}{p}$, $a \leq b, c < a + 1$, $p^* (b) = \frac{N p}{N - (a + 1 -
b) p}$ and $p^* (c) = \frac{N p}{N - (a + 1 - c) p}$.

Zhang and Wei \cite{J. Zhang} studied the existence of multiple
positive solutions for \eqref{1} with $t = s = 0$, $Q(x) = f (x)$
and $h (x) = 1$. Set $s_1 = s_2 = t$, $s = t$, $x_0 = 0$ and $Q(x)
= h (x) \equiv 1$, then problem \eqref{1} reduces to the
quasilinear elliptic system
\begin{equation}\label{1nia}
      \begin{gathered}
    - \Delta_p u  - \mu \frac{|u|^{p - 2} u}{|x|^p} = \frac{|u|^{p^* (t) - 2}
    u}{|x|^t}+
    \frac{\eta \alpha}{\alpha + \beta} \frac{ |u|^{\alpha - 2} |v|^\beta u}{|x|^t}
     + \lambda \frac{ |u|^{q - 2} u}{|x|^s},\\
   - \Delta_p v  - \mu \frac{|v|^{p - 2} v}{|x|^p} = \frac{|v|^{p^* (t) - 2} v}{|x|^t}
     + \frac{\eta \beta}{\alpha + \beta} \frac{ |u|^\alpha  |v|^{\beta -2} v}{|x|^t} +
    \theta \frac{ |v|^{q - 2} v}{|x|^s},\\
 x \in   \Omega,\\
        u = v = 0,
      \quad  x \in  \partial \Omega
        \end{gathered}
\end{equation}
where $\lambda > 0$, $\theta > 0$, $0 < \eta < \infty$, $1 < p <
N$, $0 \leq \mu < \overline{\mu} \triangleq \big(\frac{N - p}{p}
\big)^p$, $0 \leq s, t < p$, $1 \leq q < p$, $\alpha + \beta =
p^*(t) \triangleq \frac{p (N - t)}{N - p}$ is the Hardy- Sobolev
critical exponent. The author \cite{Nia1} have studied
\eqref{1nia} via the Nehari manifold. In \cite{Y. Li}, Li et al.
studied the following quasilinear elliptic problem
\begin{equation}\label{li}
      \begin{gathered}
    - \Delta_p u  - \mu \frac{|u|^{p - 2} u}{|x|^p} =
    K (x ) \frac{|u|^{p^* (s) - 2} u}{|x|^{s}} +    Q (x ) \frac{ |u|^{p^* (t) - 2}}{|x - x_0|^t}
     + \lambda f (x, u),    \quad  x \in   \Omega,\\
           u = 0,      \quad  x \in  \partial \Omega
        \end{gathered}
\end{equation}
where  $1 < p < N$, $K (x), Q (x)$ are nonnegative continuous
functions on $\overline{\Omega}$, $f$ satisfying some suitable
conditions and obtained the existence of solutions via variational
methods. For $p = 2$, $x_0 = 0$, $ K (x ) \equiv 1$ and $ Q (x )
\equiv 0$, the problem \eqref{li} has been studied.

Motivated  by the above works we study  problem \eqref{1} by using the
 Mountain-Pass Theorem of Ambrosetti and Rabinowitz. We shall show 
that  system \eqref{1} has at least two positive weak solutions.

In this article, we assume that $0 < s_1, s_2 \leq t < p$,
 $\alpha, \beta > 1$ and $\alpha + \beta = p^*(t)$.
 For $0 \leq \mu < \overline{\mu} $, we set
\begin{gather*}
\theta (\mu, s) : = \frac{p - s}{p (N - s)} A_{\mu, s}^\frac{N -s}{p - s},
\\
\theta^* : = \big\{\theta (\mu, s_1), \theta (\mu, s_2), \frac{p
- t}{p(N - t)} \frac{1}{\|Q\|_\infty^\frac{N - p}{p- t}} S_{t,
\alpha, \beta}^\frac{N - t}{p - t} \big\}.
\end{gather*}
Moreover, we assume that $Q (x)$ satisfies some of the following
assumptions:
\begin{itemize}
\item[(H1)] $Q \in C(\overline{\Omega})$, $Q (x) \geq 0$ and
$\text{meas}(\{x \in \Omega, \; h (x) > 0 \}) > 0$.

\item[(H2)] There exist $\vartheta > 0$ such that
 $Q (x_0) = \|Q\|_\infty > 0$ and $Q (x) = Q (x_0) + O(|x - x_0|^\varrho)$, as 
$x \to x_0$.

\item[(H3)] There exist $\beta_0$ and $\rho > 0$ such that $B_{2 \rho_0}
(x_0) \subset \Omega$ and $h (x) \geq \beta_0$ for all $x \in B_{2
\rho_0} (x_0)$.
\end{itemize}

Set
\begin{align*}
h_+ : = \max \{h, 0 \}, \quad  h_- : = \max \{-h, 0 \}.
\end{align*}
The main results of this article are stated in the following
two theorems.

\begin{theorem}\label{thmA}
Assume that $N \geq 3, \mu \in [0,\overline{\mu}),\;1<q<p$ and
{\rm (H1)}. Then there exists $\Lambda_{11}^* > 0$, such that for $0 <
\lambda < \Lambda_{11}^*$ problem \eqref{1} has at lest one
positive solutions.
\end{theorem}

\begin{theorem}\label{the1}
Assume that $N \geq p^2$, $0 \leq \mu < \overline{\mu}$,
 $\theta^* = \frac{p - t}{p(N - t)} \frac{1}{\|Q\|_\infty^\frac{N - p}{p- t}}
S_{t, \alpha,\beta}^\frac{N - t}{p - t}$, (H1)-(H3), $Q (0) = 0$,
$\varrho > b (\mu) p + p - N + t$ and $\frac{N -s}{b (\mu)} < q <
p$ hold, and $b(\mu)$ is the constant defined as in Lemma \ref{lem
kang}. Then there exists $\Lambda^{**}> 0$, such that for
 $0 < \lambda < \Lambda^{**}$, problem \eqref{1} has at least two 
positive solutions.
\end{theorem}

This article is divided into three sections, organized as follows.
In Section 2, we establish some elementary results. In Section 3,
we prove our main results (Theorems \ref{thmA} and \ref{the1}).

 \section{Preliminary lemmas}

The corresponding energy functional of problem
 \eqref{1} is defined by
\begin{align*}
J (u, v)  
&= \frac{1}{p} \int_\Omega \Big(|\nabla u|^p - \mu
\frac{|u|^p}{|x|^p } + |\nabla v|^p - \mu \frac{|v|^p}{|x|^p }
\Big) d x - \frac{\lambda}{q} \int_\Omega h (x)
(\frac{|u|^q}{|x|^s} +
 \frac{|v|^q}{|x|^s})  dx\\
&\quad - \frac{1}{p^*(s_1)} \int_\Omega \frac{|u|^{p^*(s_1)}
}{|x|^{s_1}} dx - \frac{1}{p^*(s_2)} \int_\Omega
\frac{|v|^{p^*(s_2)} }{|x|^{s_2}} dx \\
&\quad - \frac{1}{\alpha + \beta}
\int_\Omega Q (x) \frac{|u|^\alpha |v|^\beta }{|x - x_0|^t} dx,
\end{align*}
for each $(u, v) \in W$. Then $J \in C^1 (W, \mathbb{R})$.

\begin{lemma}\label{ps0}
Assume that $N \geq 3$,  $0 \leq \mu < \overline{\mu}$, (H1), 
$h_+\ne 0$ and $(u, v)$ is a weak solution of problem \eqref{1}. Then
there exists a positive constant $d$ depending on 
$N, |\Omega|, |h_+|_\infty, A_{\mu, s}, s_1, s_2$ and $q$ such that
\[
J(u, v) \geq - d \lambda^\frac{p}{p - q}.
\]
\end{lemma}

\begin{proof} 
Since $(u, v)$ is a weak solution of 
\eqref{1}, then, Note that $\langle J'(u, v), (u, v) \rangle = 0$,
we have
\begin{equation} \label{j0}
\begin{split}
&\langle J'(u, v), (u, v) \rangle\\
& = \int_\Omega \Big(|\nabla
u|^p - \mu \frac{|u|^p}{|x|^p } + |\nabla v|^p - \mu
\frac{|v|^p}{|x|^p } \Big) d x - \lambda \int_\Omega h (x)
(\frac{|u|^q}{|x|^s} + \frac{|v|^q}{|x|^s})  dx \\
&\quad - \int_\Omega \frac{|u|^{p^*(s_1)} }{|x|^{s_1}} dx -
\int_\Omega \frac{|v|^{p^*(s_2)} }{|x|^{s_2}} dx -  \int_\Omega Q
(x) \frac{|u|^\alpha |v|^\beta }{|x - x_0|^t} dx = 0.
\end{split}
\end{equation}
Now, by using $h_+ \ne 0$, \eqref{j0}, \eqref{hold}, the
H\"{o}lder inequality and the Sobolev-Hardy inequality, we have
 \begin{align*}
J (u, v)
& \geq  \frac{1}{p} \int_\Omega \Big(|\nabla u|^p - \mu
\frac{|u|^p}{|x|^p } + |\nabla v|^p - \mu \frac{|v|^p}{|x|^p }
\Big) d x - \frac{\lambda}{q} \int_\Omega h (x)
(\frac{|u|^q}{|x|^s} +
 \frac{|v|^q}{|x|^s})  dx\\
&\quad - \frac{1}{p^*(t)} \Big[ \int_\Omega \frac{|u|^{p^*(s_1)}
}{|x|^{s_1}} dx -  \int_\Omega \frac{|v|^{p^*(s_2)} }{|x|^{s_2}}
dx - \int_\Omega Q (x) \frac{|u|^\alpha |v|^\beta }{|x - x_0|^t}
dx \Big]
\\
& =  \frac{1}{p} \int_\Omega \Big(|\nabla u|^p - \mu
\frac{|u|^p}{|x|^p } + |\nabla v|^p - \mu \frac{|v|^p}{|x|^p }
\Big) d x - \frac{\lambda}{q} \int_\Omega h (x)
(\frac{|u|^q}{|x|^s} +
 \frac{|v|^q}{|x|^s})  dx\\
&\quad - \frac{1}{p^*(t)} \Big[ \int_\Omega \Big(|\nabla u|^p - \mu
\frac{|u|^p}{|x|^p } + |\nabla v|^p - \mu \frac{|v|^p}{|x|^p }
\Big) d x \\
&\quad - \lambda \int_\Omega h (x) (\frac{|u|^q}{|x|^s} +
 \frac{|v|^q}{|x|^s})  dx \Big]\\
& \geq  \Big (\frac{1}{p} - \frac{1}{p^*(t)} \Big) \int_\Omega \Big(|\nabla u|^p - \mu
 \frac{|u|^p}{|x|^p } + |\nabla v|^p - \mu \frac{|v|^p}{|x|^p }
 \Big) d x \\
&\quad - \lambda \Big (\frac{1}{q} - \frac{1}{p^*(t)} \Big )
\int_\Omega h (x) (\frac{|u|^q}{|x|^s} +  \frac{|v|^q}{|x|^s})  dx
\\
& \geq  \Big (\frac{1}{p} - \frac{1}{p^*(t)} \Big) (\|u\|_\mu^p + \|v\|_\mu^p) \\
&\quad -   \lambda \Big (\frac{1}{q} - \frac{1}{p^*(t)} \Big)
 \Big(\frac{N \omega_N R_0^{N - s}}{N - s}
 \Big)^\frac{p^* (s) - q}{p^* (s)} A_{\mu, s}^{-\frac{q}{p}}
 |h_+|_\infty(\|u\|_\mu^q + \|v\|_\mu^q)
\\
& \geq  2 \inf_{t \geq 0} \Big[\Big (\frac{1}{p} - \frac{1}{p^*(s)} \Big) t^p
    - \lambda \Big (\frac{1}{q} - \frac{1}{p^*(s)} \Big)
    \Big(\frac{N \omega_N R_0^{N - s}}{N - s}
 \Big)^\frac{p^* (s) - q}{p^* (s)} A_{\mu, s}^{-\frac{q}{p}}
 |h_+|_\infty t^q \Big]\\
& \geq  - d \lambda^\frac{p}{p - q}.
\end{align*}
Here $d_\Omega : = \sup_{x, y \in \Omega} |x - y|$ is the diameter
of $\Omega$ and $d$ is a positive constant depending on $N,
|\Omega|, |h_+|_\infty, A_{\mu, s}, s_1, s_2$ and $q$.
\end{proof}

Recall that a sequence $(u_n, v_n)_{n \in \mathbb{N}}$ is a
$(PS)_c$ sequence for the functional $J$ if $J (u_n, v_n)
\to c$ and $J' (u_n, v_n) \to 0$. If any $(PS)_c$
sequence $(u_n, v_n)_{n \in \mathbb{N}}$ has a convergent
subsequence, we say that $J$ satisfies the $(PS)_c$ condition.

\begin{lemma}\label{ps2}
Assume that $N \geq 3$,  $0 \leq \mu < \overline{\mu}$, {\rm (H1)}, 
$h_+ \ne 0$ and $Q (0) = 0$. Then $J(u, v)$ satisfies the $(PS)_c$
condition with $c$ satisfying
\begin{equation} \label{c}
\begin{split}
c < c_* : = \min \Big\{&\frac{p - s_1}{p (N - s_1)} A_{\mu,
s_1}^\frac{N - s_1}{p - s_1}, \frac{p - s_2}{p (N - s_2)} A_{\mu,
s_2}^\frac{N - s_2}{p - s_2}, \\
& \frac{p - t}{p(N - t)} \frac{1}{\|Q\|_\infty^\frac{N - p}{p- t}} S_{t, \alpha,
\beta}^\frac{N - t}{p - t} \Big\} 
 - d \lambda^\frac{p}{p - q}.
 \end{split}
\end{equation}
\end{lemma}

\begin{proof} 
It is easy to see that the $(PS)_c$ sequence $(u_n,
v_n)_{n \in \mathbb{N}}$ of $J (u, v)$ is bounded in $W$. Then
$(u_n, v_n) \rightharpoonup (u, v)$ weakly in $W$ as $n
\to \infty$, which implies $u_n \rightharpoonup u$ weakly
and $v_n \rightharpoonup v$ weakly in $W^{1,p}_0 (\Omega)$ as $n
\to \infty$. Passing to a subsequence we may assume that
\begin{gather*}
|\nabla u_n|^p dx \rightharpoonup \overline{\alpha}, \quad
|\nabla v_n|^p dx \rightharpoonup \widetilde{\alpha},
\\
\frac{|u_n|^p}{|x|^p} dx \rightharpoonup \overline{\beta},
\quad \frac{|v_n|^p}{|x|^p} dx \rightharpoonup \widetilde{\beta},
\\
\frac{|u_n|^{p^*(s_1)}}{|x|^{s_1}} dx \rightharpoonup \overline{\gamma}, \quad 
\frac{|v_n|^{p^*(s_2)}}{|x|^{s_2}} dx
\rightharpoonup \widetilde{\gamma},
\\
Q(x) \frac{|u_n|^\alpha |v_n|^\beta }{|x - x_0|^t} dx
\rightharpoonup \nu
\end{gather*}
weakly in the sense of measures. Using the
concentration-compactness principle in \cite{P.L. Lions}, there
exist an at most countable set $I$, a set of points 
$\{x_i \}_{i \in I} \in \Omega \setminus \{ 0 \}$, real numbers
$\overline{a}_{x_i} , \widetilde{a}_{x_i} , d_{x_i}$, $i \in I$,
$\overline{a}_0, \widetilde{a}_0, \overline{b}_0, \widetilde{b}_0,
\overline{c}_0, \widetilde{c}_0$ and $d_0$, such that
\begin{gather}
 \overline{\alpha} \geq |\nabla u|^p dx +\sum_{i
\in I}\overline{a}_{x_i}\delta_{x_i}+\overline{a}_0\delta_0, \label{n1}\\
 \widetilde{\alpha} \geq |\nabla u|^p dx +\sum_{i
\in I}\widetilde{a}_{x_i}\delta_{x_i}+\widetilde{a}_0\delta_0, \label{n2}\\
\overline{\beta} = \frac{|u|^p}
{|x|^p} dx + \overline{b}_0\delta_0, \label{n3}\\
 \widetilde{\beta} = \frac{|v|^p}
{|x|^p} dx + \widetilde{b}_0\delta_0, \label{n4}\\
 \overline{\gamma}=
\frac{|u|^{p^*(s_1)}}{|x|^{s_1}}+\overline{c}_0\delta_0, \label{n5}\\
 \widetilde{\gamma} = \frac{|v|^{p^*(s_2)}}{|x|^{s_2}} dx
+\widetilde{c}_0\delta_0, \label{n6}
\\
\nu = Q(x) \frac{|u|^\alpha |v|^\beta }{|x - x_0|^t} dx +\sum_{i
\in I}Q (x_i) d_{x_i}\delta_{x_i}+ Q (0) d_0\delta_0,\label{n7}
\end{gather}
where $\delta_x$ is the Dirac-mass of mass 1 concentrated at the
point $x$.

 First, we consider the possibility of the concentration at
$\{x_i\}_{i \in I}\in\Omega \setminus \{0\}$.
Let $\epsilon>0$ be
small enough, take $\eta_{x_i} \in C_c^\infty (B_{2 \varepsilon}
(x_i))$, such that $\eta_{x_i}|_{B_{\varepsilon} (x_i)} = 1$, $0
\leq \eta_{x_i} \leq 1$ and $|\nabla \eta_{x_i} (x)| \leq
\frac{C}{\varepsilon}$. Then
\begin{align*}
o (1)&= 
\langle J' (u_n, v_n), (\eta_{x_i}^p u_n, \eta_{x_i}^p v_n) \rangle\\
&= \int_{\Omega}\left(|\nabla u_n|^{p-2}\nabla u_n
\nabla(\eta_{x_i}^p u_n)+ |\nabla v_n|^{p-2}\nabla v_n
\nabla(\eta_{x_i}^p v_n)\right)dx\\
&\quad - \int_{\Omega} Q(x) \frac{|u_n|^\alpha |v_n|^\beta}{|x - x_0|^t}
\eta_{x_i}^p dx
 -\mu \int_{\Omega}\Big(\frac{|u_n|^{p}}{|x|^p} \eta_{x_i}^p +
\frac{|v_n|^{p}}{|x|^p} \eta_{x_i}^p \Big)dx\\
 &\quad - \lambda\int_{\Omega} h (x) \Big(
\frac{|u_n|^q}{|x|^s} \eta_{x_i}^p + \frac{|v_n|^q}{|x|^s} \eta_{x_i}^p \Big)dx\\
&\quad -\int_{\Omega}\frac{|u_n|^{p^*(s_1)}}{|x|^{s_1}} \eta_{x_i}^p
dx - \int_{\Omega} \frac{|v_n|^{p^*(s_2)}}{|x|^{s_2}}\eta_{x_i}^p dx.
\end{align*}
From \eqref{n3}-\eqref{n7}, one can obtain
\begin{gather}
\lim_{\varepsilon \to 0} \lim_{n \to \infty}
\int_{\Omega}\Big(\frac{|u_n|^{p}} {|x|^p} \eta_{x_i}^p
+\frac{|v_n|^{p}}{|x|^p} \eta_{x_i}^p \Big)dx =
\lim_{\varepsilon \to 0}\Big(\int_{\Omega} \eta_{x_i}^p d
\overline{\beta} + \int_{\Omega}
\eta_{x_i}^p d\widetilde{\beta} \Big) = 0, \label{s3}\\
 \lim_{\varepsilon \to 0}
\lim_{n\to\infty}\int_{\Omega}
\Big(\frac{|u_n|^{p^*(s_1)}} {|x|^{s_1}} \eta_{x_i}^p
+\frac{|v_n|^{p^*(s_2)}}{|x|^{s_2}} \eta_{x_i}^p \Big)dx =
\lim_{\varepsilon \to 0}\Big(\int_{\Omega}\eta_{x_i}^p d
\overline{\gamma} + \int_{\Omega}
\eta_{x_i}^p d \widetilde{\gamma}\Big) = 0, 
\nonumber\\
 \lim_{\varepsilon \to 0}
\lim_{n\to\infty}\int_{\Omega} h(x)
\Big(\frac{|u_n|^q}{|x|^s} \eta_{x_i}^p + \frac{|v_n|^q}{|x|^s}
\eta_{x_i}^p \Big) dx = 0,
\nonumber\\
\lim_{\varepsilon \to 0}\lim_{n \to \infty}
\int_{\Omega} Q(x) \frac{|u_n|^\alpha |v_n|^\beta}{|x - x_0|^t}
\eta_{x_i}^p dx =\lim_{\varepsilon \to
0}\int_{\Omega}\eta_{x_i}^p d \nu = Q (x_i) d x_i. \nonumber
\end{gather}
Thus,
\begin{equation}\label{e2.3}
\begin{split}
0&=\lim_{\varepsilon \to 0}
\lim_{n\to\infty}\int_{\Omega} \Big(|\nabla
u_n|^{p-2}\nabla u_n\nabla(\eta_{x_i}^p u_n)\\
&\quad + |\nabla
v_n|^{p-2}\nabla v_n\nabla(\eta_{x_i}^p v_n)\Big) dx - Q (x_i) d
x_i.
\end{split}
\end{equation}
Moreover, we have
\begin{equation} \label{e2.4}
\begin{split}
&\lim_{\varepsilon\to0}\lim_{n\to\infty}\big|\int_{\Omega}u_n
 |\nabla u_n|^{p-2}\nabla u_n \nabla \eta_{x_i}^p
dx\big|\\
&\leq  \lim_{\varepsilon\to0}
\lim_{n\to\infty}\Big(\int_{\Omega}|u_n|^p
 |\nabla \eta_{x_i}^p|^p dx \Big)^{1/p}
\Big(\int_{\Omega}|\nabla u_n|^p dx \Big)^{\frac{p-1}{p}}
\\
&\leq  C \lim_{\varepsilon\to0}\lim_{n\to\infty}
\Big(\int_{\Omega}|u_n|^p
|\nabla \eta_{x_i}^p|^p dx \Big)^{1/p}
\\
& =  C \lim_{\varepsilon\to0} \Big(\int_{\Omega}|u|^p
|\nabla \eta_{x_i}^p|^p dx \Big)^{1/p}
\\
&\leq  C\lim_{\varepsilon \to 0} \Big(\int_{B_\varepsilon (x_i)}
|\nabla \eta_{x_i}^p |^N dx\Big)^{1/N}
\Big(\int_{B_\varepsilon (x_i)}|u|^{p^*} dx\Big)^{1/P^*}
\\
&\leq  C\lim_{\varepsilon \to 0}  \Big(\int_{B_\varepsilon (x_i)}
|u|^{p^*} dx\Big)^{1/P^*} = 0.\\
\end{split}
\end{equation}
Similarly,
\begin{equation}\label{e2.5}
\lim_{\varepsilon
\to0}\lim_{n\to\infty}\big|\int_{\Omega} v_n
|\nabla v_n|^{p-2}\nabla v_n \nabla \eta_{x_i}^p dx\big|=0.
\end{equation}
Combining \eqref{e2.3}-\eqref{e2.5}, there holds
\begin{equation}\label{e2.6}
\begin{split}
0 & =  \lim_{\varepsilon \to0}\lim_{n\to\infty}
\int_{\Omega}(|\eta_{x_i} \nabla u_n|^p + |\eta_{x_i} \nabla
v_n|^p) dx - Q (x_i) d_{x_i} \\
& =  \lim_{\varepsilon \to 0}
\int_{\Omega}\left(\eta_{x_i} ^p d\overline{\alpha}+ \eta_{x_i} d
\widetilde{\alpha}\right)- Q (x_i) d_{x_i}.
\end{split}
\end{equation}
On the other hand, \eqref{hardy2} implies
\begin{equation} \label{s4}
\begin{split}
& \frac{1}{\|Q\|_\infty^\frac{p}{p^* (t)}} S_{t,\alpha, \beta}
 \Big(\int_\Omega Q (x) \frac{|\eta_{x_i} u_n|^\alpha
|\eta_{x_i}
v_n|^\beta}{|x - x_0|^t} dx \Big)^{\frac{p}{p^* (t)}} \\
& \leq \int_\Omega \Big(|\nabla (\eta_{x_i} u_n)|^p
+ |\nabla (\eta_{x_i} v_n)|^p - \mu \frac{|\eta_{x_i} u_n|^p  +
|\eta_{x_i} v_n|^p}{|x|^p } \Big) dx.
\end{split}
\end{equation}
Note that
\[
\lim_{\varepsilon \to0}\lim_{n\to\infty}
\int_{\Omega} |\nabla \eta_{x_i}|^p |u_n|^p dx  =
\lim_{\varepsilon \to0}\lim_{n\to\infty}
\int_{\Omega} |\nabla \eta_{x_i}|^p |v_n|^p dx = 0.
\]
From this equality, \eqref{e2.4} and \eqref{e2.5}, we obtain
\begin{gather}\label{s1}
\lim_{\varepsilon \to0}\lim_{n\to\infty}
\int_{\Omega} |\eta_{x_i} \nabla u_n|^p dx  = \lim_{\varepsilon
\to0}\lim_{n\to\infty} \int_{\Omega} |\nabla(
\eta_{x_i} u_n)|^p dx,
\\ \label{s2}
\lim_{\varepsilon \to0}\lim_{n\to\infty}
\int_{\Omega} |\eta_{x_i} \nabla v_n|^p dx  = \lim_{\varepsilon
\to0}\lim_{n\to\infty} \int_{\Omega} |\nabla(
\eta_{x_i} v_n)|^p dx.
\end{gather}
Relations \eqref{n7}, \eqref{s3} and \eqref{s4}-\eqref{s2}
imply
\begin{equation} \label{s5}
 \frac{1}{\|Q\|_\infty^\frac{p}{p^* (t)}} S_{t,\alpha, \beta}
 \big(Q (x_i) d_{x_i} \big)^{\frac{p}{p^* (t)}} \leq
\lim_{\varepsilon \to0} \int_{\Omega} |\eta_{x_i}|^p d
\overline{\alpha} + \lim_{\varepsilon \to0} \int_{\Omega}
|\eta_{x_i}|^p d \widetilde{\alpha}.
\end{equation}
Combining \eqref{e2.6} and \eqref{s5},
\begin{equation}\label{e2.7}
\frac{1}{\|Q\|_\infty^\frac{p}{p^* (t)}} S_{t,\alpha, \beta}
\big(Q (x_i) d_{x_i} \big)^{\frac{p}{p^* (t)}} \leq Q (x_i)
d_{x_i},
\end{equation}
which implies that either
\begin{equation}\label{e2.9}
 Q (x_i) d_{x_i} = 0,\quad\text{or}\quad
 Q (x_i) d_{x_i} \geq \frac{1}{\|Q\|_\infty^\frac{N - p}{p - t}}
S^{\frac{N - t}{p - t}}_{t,\alpha, \beta} .
\end{equation}
Now, we consider the possibility of the concentration at 0.
For $\epsilon>0$ be small enough, take
 $\eta_0 \in C_c^\infty (B_{2 \varepsilon} (0))$, such that
 $\eta_0|_{B_{\varepsilon} (0)} = 1$, $0 \leq \eta_0 \leq 1$ and
$|\nabla \eta_0 (x)| \leq \frac{C}{\varepsilon}$. Then
\begin{align*}
o (1)&= \langle J' (u_n, v_n), (\eta_0^p u_n, 0) \rangle\\
&= \int_{\Omega} |\nabla u_n|^{p-2}\nabla u_n \nabla(\eta_0^p
u_n)dx
 -\mu \int_{\Omega}\frac{|u_n|^{p}}{|x|^p} \eta_0^p dx - \lambda \int_{\Omega} h (x)
\frac{|u_n|^q}{|x|^s} \eta_0^p dx \\
&\quad -\int_{\Omega}\frac{|u_n|^{p^*(s_1)}}{|x|^{s_1}} \eta_0^p dx -
\frac{\alpha}{\alpha + \beta} \int_{\Omega} Q(x)
\frac{|u_n|^\alpha |v_n|^\beta}{|x - x_0|^t} \eta_0^p dx.
\end{align*}
From \eqref{n3}, \eqref{n5}, \eqref{n7} and $Q (0) = 0$, we obtain
\begin{gather*}
 \lim_{\varepsilon \to 0} \lim_{n \to \infty}
\int_{\Omega} \frac{|u_n|^{p}} {|x|^p} \eta_0^p dx =
\overline{b}_0, \;\;\;\;\; \lim_{\varepsilon \to 0}
\lim_{n \to \infty} \int_{\Omega} \frac{|u_n|^{p^* (s_1)}}
{|x|^{s_1}} \eta_0^p dx = \overline{c}_0,
\\
\lim_{\varepsilon \to 0}\lim_{n \to \infty}
\int_{\Omega} Q(x) \frac{|u_n|^\alpha |v_n|^\beta}{|x - x_0|^t}
\eta_0^p dx = \lim_{\varepsilon \to 0} \lim_{n \to
\infty} \int_{\Omega} h (x) \frac{|u_n|^q}{|x|^s} \eta_0^p dx = 0.
\end{gather*}
Thus,
\begin{equation}\label{s6}
0=\lim_{\varepsilon \to 0}
\lim_{n\to\infty}\int_{\Omega} |\nabla u_n|^{p-2}\nabla
u_n\nabla(\eta_0^p u_n) dx - \mu \overline{b}_0 - \overline{c}_0.
\end{equation}
Note that
\begin{equation*}
\lim_{\varepsilon \to0}\lim_{n\to\infty}
\int_{\Omega} u_n |\nabla u_n|^{p-2}\nabla u_n \nabla \eta_0^p dx
=0.
\end{equation*}
This equality and \eqref{s6} yield
\begin{equation}\label{s7}
\lim_{\varepsilon \to 0} \int_{\Omega} \eta_0^p
d\overline{\alpha} - \mu \overline{b}_0 = \overline{c}_0.
\end{equation}
On the other hand, \eqref{san} implies
\[
A_{\mu, s_1} \Big(\int_\Omega \frac{|\eta_0 u_n|^{p^*
(s_1)}}{|x|^{s_1}} dx \Big)^{\frac{p}{p^* (s_1)}} \leq \int_\Omega
\Big(| \nabla (\eta_0 u_n)|^p - \mu \frac{|\eta_0 u_n|^p}{|x|^p }
\Big) dx.
\]
Thus
\begin{equation} \label{s11}
A_{\mu, s_1} \overline{c}_0^{\frac{p}{p^* (s_1)}} \leq
\lim_{\varepsilon \to0}\lim_{n\to\infty}
 \int_\Omega
| \nabla (\eta_0 u_n)|^p d x - \mu \overline{b}_0.
\end{equation}
Note that
\[
\lim_{\varepsilon \to0}\lim_{n\to\infty}
\int_{\Omega} |\eta_0 \nabla u_n|^p dx  = \lim_{\varepsilon
\to0}\lim_{n\to\infty} \int_{\Omega} |\nabla(
\eta_0 u_n)|^p dx,
\]
which together with \eqref{s11} imply
\begin{equation} \label{s12}
A_{\mu, s_1} \overline{c}_0^{\frac{p}{p^* (s_1)}} \leq
\lim_{\varepsilon \to0} \int_\Omega |\eta_0|^p d
\overline{\alpha} - \mu \overline{b}_0.
\end{equation}
Therefore, from \eqref{s7} and \eqref{s12},
\begin{equation} \label{s13}
A_{\mu, s_1} \overline{c}_0^{\frac{p}{p^* (s_1)}} \leq
\overline{c}_0,
\end{equation}
which implies that either
\begin{equation}\label{s21}
 \overline{c}_0 = 0,\quad\text{or}\quad
\overline{c}_0 \geq A_{\mu, s_1} ^{\frac{N - s_1}{p - s_1}}.
\end{equation}
Similarly, either
\begin{equation}\label{s22}
 \overline{c}_0 = 0,\quad\text{or}\quad
\overline{c}_0 \geq A_{\mu, s_2} ^{\frac{N - s_2}{p - s_2}}.
\end{equation}
Recall that $u_n \rightharpoonup u$ weakly and
 $v_n \rightharpoonup v$ weakly in $W^{1,p}_0 (\Omega)$, we have
\begin{equation} \label{s14}
\begin{split}
&c + o (1) \\
&=  J (u_n, v_n)  \\
&=  \frac{1}{p} \int_\Omega
\Big(|\nabla u_n - \nabla u|^p - \mu \frac{|u_n - u|^p}{|x|^p } +
|\nabla v_n - \nabla v|^p - \mu \frac{|v_n - v|^p}{|x|^p } \Big) d
x  \\
&\quad - \frac{1}{p^*(s_1)} \int_\Omega \frac{|u_n - u|^{p^*(s_1)}
}{|x|^{s_1}} dx - \frac{1}{p^*(s_2)} \int_\Omega \frac{|v_n -
v|^{p^*(s_2)} }{|x|^{s_2}} dx  \\
&\quad - \frac{1}{p^* (t)} \int_\Omega Q (x) \frac{|u_n - u|^\alpha
|v_n - v|^\beta }{|x - x_0|^t} dx + J (u, v).
\end{split}
\end{equation}
On the other hand, from $o(1) = J' (u_n, v_n)$, we obtain
\begin{equation} \label{s15}
J' (u_n, v_n)  = 0.
\end{equation}
Thus, $0 = \langle J '(u, v), (u, v) \rangle$, which together with
$o(1) = \langle J '(u_n, v_n), (u_n, v_n) \rangle$ imply
\begin{equation} \label{s16}
\begin{split}
o (1) 
&=   \int_\Omega \Big(|\nabla u_n - \nabla u|^p - \mu
\frac{|u_n - u|^p}{|x|^p } + |\nabla v_n - \nabla v|^p - \mu
\frac{|v_n - v|^p}{|x|^p } \Big) dx  \\
&\quad -  \int_\Omega \frac{|u_n - u|^{p^*(s_1)} }{|x|^{s_1}} dx -
\int_\Omega \frac{|v_n - v|^{p^*(s_2)} }{|x|^{s_2}} dx  \\
&\quad -\int_\Omega Q (x) \frac{|u_n - u|^\alpha |v_n - v|^\beta }{|x -
x_0|^t} dx.
\end{split}
\end{equation}
From \eqref{s14}-\eqref{s16} and Lemma \ref{ps0},
\begin{equation} \label{s17}
\begin{split}
c + o (1)
& \geq   \frac{p - s_1}{p (N - s_1)} \int_\Omega
\frac{|u_n - u|^{p^*(s_1)} }{|x|^{s_1}} dx + \frac{p - s_2}{p (N -
s_2)} \int_\Omega \frac{|v_n - v|^{p^*(s_2)} }{|x|^{s_2}} dx  \\
&\quad + \frac{p - t}{p(N - t)} \int_\Omega Q (x) \frac{|u_n -
u|^\alpha |v_n - v|^\beta }{|x - x_0|^t} dx - d \lambda^\frac{p}{p- q}.
\end{split}
\end{equation}
Passing to the limit in \eqref{s17} as $n \to \infty$, we
have
\begin{equation} \label{s18}
c   \geq   \frac{p - s_1}{2 (N - s_1)} \overline{c}_0 + \frac{p
- s_2}{p (N - s_2)} \widetilde{c}_0 + \frac{p - t}{p(N - t)}
\sum_{i \in I} Q (x_i) d_{x_i} - d \lambda^\frac{p}{p - q}.
\end{equation}
By the assumption $c < c_*$ and in view of \eqref{e2.9},
\eqref{s21} and \eqref{s22}, there holds
$\overline{c}_0 = \widetilde{c}_0 = 0$, $Q (x_i) d_{x_i} = 0$, $i \in I$.
 Up to a subsequence, $(u_n, v_n) \to (u, v)$ strongly in $W$ as $n
\to \infty$.
\end{proof}

When the restriction $Q (0) = 0$ is removed, we establish the
following version of Lemma \ref{ps2}.


\begin{lemma}\label{ps3}
Assume that $N \geq 3$,  $0 \leq \mu < \overline{\mu}$, (H1) and
$h_+ \ne 0$. Then $J(u, v)$ satisfies the $(PS)_c$ condition with
$c$ satisfying
\begin{equation} \label{c2}
\begin{split}
c < c_0 : = \min \Big\{&\frac{p - s_1}{p (N - s_1)}
\Big (\frac{1}{p} A_{\mu, s_1} \Big)^\frac{N - s_1}{p - s_1},
\frac{p -s_2}{p (N - s_2)} \Big ( \frac{1}{p} A_{\mu, s_2} \Big)^\frac{N -
s_2}{p - s_2},  \\
& \frac{p - t}{p(N - t)} \frac{\Big ( \frac{1}{p}
S_{t, \alpha, \beta} \Big )^\frac{N - t}{p -
t}}{\|Q\|_\infty^\frac{N - p}{p- t}}  \Big\} - d
\lambda^\frac{p}{p - q}.
 \end{split}
\end{equation}
\end{lemma}

The proof of the above lemma is similar to Lemma \ref{ps2} and is
omitted. 

\begin{lemma}[\cite{Kang}] \label{lem kang}
Assume that $1 < p < N$, $0 \leq t < p$ and $0 \leq \mu <
\overline{\mu}$. Then the limiting problem
\begin{gather*}
 - \Delta_p u  - \mu \frac{|u|^{p - 1} }{|x|^p} =
    \frac{|u|^{p^* (t) - 1} }{|x|^t},
    \quad  \text{in  } \mathbb{R}^N \setminus \{0 \},\\
   u \in D^{1, p} (\mathbb{R}^N), \quad 
 u > 0,    \quad  \text{in  } \mathbb{R}^N \setminus \{0 \},
\end{gather*}
has positive radial ground states
\begin{equation} \label{ve}
V_\epsilon (x) \triangleq \epsilon^\frac{p - N}{p} U_{p,
\mu}(\frac{x}{\epsilon}) = \epsilon^\frac{p - N}{p} U_{p, \mu}
(\frac{|x|}{\epsilon}), \quad \forall \epsilon > 0,
\end{equation}
that satisfy
\[
\int_\Omega \Big(|\nabla V_\epsilon (x)|^p - \mu \frac{|V_\epsilon
(x)|^p}{|x|^p } \Big) dx = \int_\Omega \frac{|V_\epsilon
(x)|^{p^*(t)}}{|x|^t } dx = (A_{\mu, t})^\frac{N - t}{p - t},
\]
where $U_{p, \mu} (x) = U_{p, \mu} (|x|)$ is the unique radial
solution of the limiting problem with
\[
U_{p, \mu}(1) = \Big ( \frac{(N - t) (\overline{\mu} - \mu)}{N -
p} \Big)^\frac{1}{p^*(t) - p}.
\]
Furthermore, $U_{p, \mu}$ has the following properties:
\begin{gather*}
 \lim_{r \to 0} r^{a(\mu)} U_{p, \mu}(r) = C_1 > 0,\quad
 \lim_{r \to + \infty} r^{b(\mu)} U_{p, \mu}(r) = C_2 > 0,\\
 \lim_{r \to 0} r^{a(\mu) + 1} |U'_{p, \mu}(r)| = C_1 a(\mu) \geq 0,\\
 \lim_{r \to + \infty} r^{b(\mu) + 1} |U'_{p, \mu}(r)| = C_2 b (\mu) > 0,
\end{gather*}
where $C_i$ $(i = 1, 2)$ are positive constants and $a(\mu)$ and
$b(\mu)$ are zeros of the function
\[
f(\zeta) = (p - 1) \zeta^p - (N - p) \zeta^{p - 1} + \mu, \quad
\zeta \geq 0, \; 0 \leq \mu < \overline{\mu},
\]
that satisfy
\[
0 \leq a(\mu) < \frac{N - p}{p} < b(\mu) \leq \frac{N - p}{p - 1}.
\]
\end{lemma}

\begin{lemma}\label{lem11}
Under the assumptions of Theorem \ref{the1}, there exists $(u_1,
v_1) \in W \setminus \{(0, 0)\}$ and $\Lambda_1 > 0$, such that
for $0 < \lambda < \Lambda_1$, there holds
\begin{equation} \label{3.11}
\sup_{t \geq 0} J (t u_1, t v_1) < \frac{p - t}{p(N - t)}
\frac{1}{\|Q\|_\infty^\frac{N - p}{p - t}} S^\frac{N - t}{p -
t}_{t,\alpha, \beta} (\mu) - d \lambda^\frac{p}{p - q}.
\end{equation}
\end{lemma}

\begin{proof}
 First, we give some estimates on the extremal
function $V_\epsilon (x)$ defined in \eqref{ve}. For $m \in
\mathbb{N}$ large, choose $\varphi (x) \in C_0^\infty
(\mathbb{R}^N)$, $0 \leq \varphi (x) \leq 1$, $\varphi (x) = 1$
for $|x| \leq \frac{1}{2m}$, $\varphi (x) = 0$ for $|x| \geq
\frac{1}{m}$, $\|\nabla \varphi (x)\|_{L^p (\Omega)} \leq 4 m$,
set $u_\epsilon (x) = \varphi (x) V_\epsilon (x)$. For $\epsilon
\to 0$, the behavior of $u_\epsilon$ has to be the same as
that of $V_\epsilon$, but we need precise estimates of the error
terms. For $1 < p < N$, $0 \leq t < p$ and $1 < q < p^*(s)$, we
have the following estimates \cite{D. Cao}:
\begin{gather}
\int_\Omega \Big(|\nabla u_\epsilon|^p - \mu
\frac{|u_\epsilon|^p}{|x|^p } \Big) dx  
 = (A_{\mu, t})^\frac{N-t}{p - t} + O (\epsilon ^{b(\mu) p + p - N}), \label{e8}
\\
\int_\Omega \frac{|u_\epsilon|^{p^*(t)}}{|x|^t } dx  
=  (A_{\mu, t})^\frac{N - t}{p - t} + O (\epsilon ^{b(\mu) p^*(t) - N +
t}), \label{e9}
\\
\int_\Omega \frac{|u_\epsilon|^q}{|x|^s } dx 
 \geq  \begin{cases}
   C \epsilon^{N - s + (1 - \frac{N}{p})q},
    \ & q > \frac{N - s}{b (\mu)},\\
    C \epsilon^{N - s + (1 - \frac{N}{p})q} |\ln \epsilon|,
     & q = \frac{N - s}{b (\mu)},\\
     C \epsilon^{q (b(\mu) + 1 - \frac{N}{p})q},
     & q < \frac{N - s}{b (\mu)}.
        \end{cases}\label{e10}
\end{gather}

Now, we consider the functional $I : W \to \mathbb{R}$
defined by
\[
I (u, v)  =  \frac{1}{p} \int_\Omega \Big(|\nabla u|^p - \mu
\frac{|u|^p}{|x|^p } + |\nabla v|^p - \mu \frac{|v|^p}{|x|^p }\Big) d x  
  - \frac{1}{p^* (t)} \int_\Omega Q (x) \frac{|u|^\alpha
|v|^\beta }{|x - x_0|^t} dx.
\]
 Let $u_1 = \alpha^\frac{1}{p}
u_\epsilon$, $v_1 = \beta^\frac{1}{p} u_\epsilon$ and define the
function $g_1 (t) : = J (t u_1,
 t v_1)$, $t \geq 0$.  Note that $\lim_{t \to +
\infty} g_1(t) = - \infty$ and $g_1(t) > 0$ as $t$ is close to
$0$. Thus $\sup_{t \geq 0} g_1(t)$ is attained at some finite
$t_\epsilon > 0$ with $g'_1(t_\epsilon) = 0$. Furthermore, 
$C' < t_\epsilon < C''$; where $C'$ and $C''$ are the positive 
constants independent of $\epsilon$. We have
\begin{equation}\label{sa3}
I(t u_1, t v_1) =   y (t u_1, t v_1) -  \frac{t^{p^*(t)}}{p^*(t)}
(\alpha^\frac{\alpha}{p} \beta^\frac{\beta}{p}) \int_\Omega (Q (x)
- Q (x_0)) \frac{|u_\epsilon|^{p^*(t)}}{|x - x_0|^t } dx.
\end{equation}
where
\begin{align*}
&y (t u_1, t v_1)\\
&:=
\Big [\frac{t^p}{p} (\alpha + \beta)
\int_\Omega \Big(|\nabla u_\epsilon|^p - \mu
\frac{|u_\epsilon|^p}{|x|^p } \Big) dx - \frac{t^{p^*(t)}}{p^*(t)}
(\alpha^\frac{\alpha}{p} \beta^\frac{\beta}{p}) \int_\Omega Q
(y_0) \frac{|u_\epsilon|^{p^*(t)}}{|x - x_0|^t } dx \Big ]\,.
\end{align*}
Note that
\begin{equation} \label{ina}
\sup_{t \geq 0} \big(\frac{t^p}{p} A - \frac{t^{p^*(t)}}{p^*(t)} B \big) 
=\big(\frac{1}{p} - \frac{1}{p^*(t)} \big)
\big(\frac{A}{B^\frac{p}{p^*(t)}}\big)^{\frac{p^*(t)}{p^*(t) -
p}}, \quad A, B > 0.
\end{equation}
From (H2), \eqref{e8}, \eqref{e9} and \eqref{ina} it follows that
\begin{align}
&\sup_{t \geq 0} y(t u_1, t v_1) \nonumber\\
& =  y (t_\epsilon u_1, t_\epsilon v_1) \nonumber\\
&\leq  \frac{p - t}{p(N - t)} \Big( \frac{(\alpha +
\beta) \int_\Omega \Big(|\nabla u_\epsilon|^p - \mu
\frac{|u_\epsilon|^p}{|x|^p } \Big) dx }{((\alpha^\frac{\alpha}{p}
\beta^\frac{\beta}{p}) \|Q\|_\infty \int_\Omega
\frac{|u_\epsilon|^{p^*(t)}}{|x - x_0|^t } dx)^\frac{N - p}{N -t}}
\Big)^\frac{N - t}{p - t}
\nonumber\\
&\leq  \frac{p - t}{p(N - t)} \frac{1}{\|Q\|_\infty^\frac{N -
p}{p - t}} \Big[ \Big(\Big (\frac{\alpha}{\beta}
\Big)^\frac{\beta}{\alpha + \beta} + \Big (\frac{\beta}{\alpha}
\Big)^\frac{\alpha}{\alpha + \beta} \Big) \frac{(A_{\mu,
t})^\frac{N - t}{p - t} + O (\epsilon ^{b(\mu) p + p -
N})}{(A_{\mu,t})^\frac{N - t}{p - t} + O (\epsilon ^{b(\mu) p^*(t)
- N + t})} \Big]^{\frac{N - t}{p - t}}
\nonumber \\
&\leq  \frac{p - t}{p(N - t)} \frac{1}{\|Q\|_\infty^\frac{N -
p}{p - t}} \Big[ \Big(\Big (\frac{\alpha}{\beta}
\Big)^\frac{\beta}{\alpha + \beta} + \Big (\frac{\beta}{\alpha}
\Big)^\frac{\alpha}{\alpha + \beta} \Big) (A_{\mu, t})^\frac{N -
t}{p - t} \Big]^{\frac{N - t}{p - t}} + O (\epsilon
^{b(\mu) p + p - N})
\nonumber\\
&\leq  \frac{p - t}{p(N - t)} \frac{1}{\|Q\|_\infty^\frac{N -
p}{p - t}} \Big(S_{t, \alpha, \beta}  \Big)^\frac{N - t}{p - t}  +
O (\epsilon ^{b(\mu) p + p - N}). \label{sa4}
\end{align}
On the other hand, (H2) implies that there exists $r_1 < r$, such
that for $x \in B_{r_1} (y_0)$, $|Q (x) - Q (x_0)| \leq C |x -
x_0|^{\vartheta}$. Thus
\begin{equation} \label{sa5}
\begin{split}
\big| \int_\Omega (Q (x) - Q (x_0))
\frac{|u_\epsilon|^{p^*(t)}}{|x - x_0|^t } dx \big|
&\leq  C \int_\Omega |Q (x) - Q (x_0)| \frac{|u_\epsilon|^{p^*(t)}}{|x -
x_0|^t } dx \\
& =  C \int_{B_{2 r} (x_0)}  \frac{|x - x_0|^\vartheta
|u_\epsilon |^{p^*(t)}}{|x -
x_0|^t } dx
=  O ( \epsilon^{\vartheta - t})
\end{split}
\end{equation}
From \eqref{sa3}, \eqref{sa4} and \eqref{sa5}, we conclude that
\begin{equation} \label{sa6}
\begin{split}
\sup_{t \geq 0} I(t u_1, t v_1)
&= I (t_\epsilon u_1, t_\epsilon v_1) \\
&\leq \frac{p - t}{p(N - t)} \frac{1}{\|Q\|_\infty^\frac{N -
p}{p - t}} \Big(S_{t, \alpha, \beta}  \Big)^\frac{N - t}{p - t}  +
O (\epsilon ^{b(\mu) p + p - N}).
\end{split}
\end{equation}
Observe that there exists $\Lambda_1^* > 0$, such that for
 $0 < \lambda < \Lambda_1^*$ and
\[
\frac{p - t}{p(N - t)}
\frac{1}{\|Q\|_\infty^\frac{N - p}{p - t}} \Big(S_{t, \alpha,
\beta}  \Big)^\frac{N - t}{p - t}  - d \lambda^\frac{p}{p - q} >0.
\]
Then for $0 < \lambda < \Lambda_1^*$,
 there exists $t_1 \in (0, 1)$, such that
\begin{equation} \label{sa1}
\begin{split}
\sup_{0 \leq t \leq t_1} J (t u_1, t v_1)
&\leq  \sup_{0 \leq t \leq t_1} \frac{1}{p} t^p \int_\Omega \Big(|\nabla u_1|^p +
|\nabla v_1|^p  \Big) dx \\
& <  \frac{p - t}{p(N - t)} \frac{1}{\|Q\|_\infty^\frac{N - p}{p
- t}} \Big(S_{t, \alpha, \beta}  \Big)^\frac{N - t}{p - t}  - d
\lambda^\frac{p}{p - q}.
\end{split}
\end{equation}
On the other hand,
\begin{equation} \label{sa2}
\begin{split}
&\sup_{t \geq t_1} J (t u_1, t v_1)\\
&\leq  \sup_{t \geq t_1} \Big[I (t u_1, t v_1)
 - \frac{\lambda}{q} t^q \int_\Omega h (x)\frac{|u_1|^q}{|x|^s}  dx
 - \frac{1}{p^*(s_1)} t^{p^*(s_1)} \int_\Omega
\frac{|u_1|^{p^*(s_1)} }{|x|^{s_1}} dx \Big ] \\
&\leq  \sup_{t \geq t_1} I (t u_1, t v_1) - \frac{\lambda}{q}
t_1^q \int_\Omega h (x) |u_1|^q  dx - \frac{1}{p^*(s_1)}
t_1^{p^*(s_1)} \int_\Omega \frac{|u_1|^{p^*(s_1)} }{|x|^{s_1}} dx\Big ]
\\
&\leq  \frac{p - t}{p(N - t)} \frac{1}{\|Q\|_\infty^\frac{N -
p}{p - t}} \Big(S_{t, \alpha, \beta}  \Big)^\frac{N - t}{p - t}  +
O (\epsilon ^{b(\mu) p + p - N})  \\
&\quad - C \int_\Omega \frac{|u_\epsilon|^{p^*(s_1)} }{|x|^{s_1}} dx -
\lambda C \int_\Omega h (x) \frac{|u_\epsilon|^q}{|x|^s}  dx
\end{split}
\end{equation}
From \eqref{e9},
\begin{equation} \label{san1}
\int_\Omega \frac{|u_\epsilon|^{p^*(s_1)}}{|x|^{s_1} } dx \geq O
(\epsilon ^{b(\mu) p^*(s_1) - N + s_1}).
\end{equation}
Also, from \eqref{e10}, it follows that
\begin{equation} \label{san2}
\int_\Omega h (x) \frac{|u_\epsilon|^q}{|x|^s} dx \geq \beta_0
\int_\Omega \frac{|u_\epsilon|^q}{|x|^s} dx \geq
\begin{cases}
   C \epsilon^{N - s + (1 - \frac{N}{p})q},
      & q > \frac{N - s}{b (\mu)},\\
    C \epsilon^{N - s + (1 - \frac{N}{p})q} |\ln \epsilon|,
      & q = \frac{N - s}{b (\mu)},\\
     C \epsilon^{q (b(\mu) + 1 - \frac{N}{p})q},
      & q < \frac{N - s}{b (\mu)}.
        \end{cases}
\end{equation}
Since $q \geq \frac{N - s}{b (\mu)}$, by \eqref{sa2}-\eqref{san2}
we have
\begin{align*}
\sup_{t \geq t_1} J (t u_1, t v_1)
&\leq  \frac{p - t}{p(N - t)}
\frac{1}{\|Q\|_\infty^\frac{N - p}{p - t}} 
\big(S_{t, \alpha,\beta}\big)^\frac{N - t}{p - t}  +
O (\epsilon ^{b(\mu) p + p - N})  \\
&\quad + O (\epsilon ^{b(\mu) p^*(s_1) - N + s_1}) - \lambda
\begin{cases}
   C \epsilon^{N - s + (1 - \frac{N}{p})q},
      & q > \frac{N - s}{b (\mu)},\\
    C \epsilon^{N - s + (1 - \frac{N}{p})q} |\ln \epsilon|,
      & q = \frac{N - s}{b (\mu)}.
        \end{cases}
\end{align*}
Note that $b(\mu) p + p - N < b(\mu) p^*(s_1) - N + s_1$, then we
have
\begin{equation} \label{san3}
\begin{split}
\sup_{t \geq t_1} J (t u_1, t v_1)
&\leq   \frac{p - t}{p(N -t)} \frac{1}{\|Q\|_\infty^\frac{N -
p}{p - t}} \big(S_{t, \alpha, \beta}  \big)^\frac{N - t}{p - t}    \\
&\quad + O (\epsilon ^{b(\mu) p^*(s_1) - N + s_1}) - \lambda
\begin{cases}
   C \epsilon^{N - s + (1 - \frac{N}{p})q},
      & q > \frac{N - s}{b (\mu)},\\
    C \epsilon^{N - s + (1 - \frac{N}{p})q} |\ln \epsilon|,
     & q = \frac{N - s}{b (\mu)}.
        \end{cases}
\end{split}
\end{equation}
Note that $N > p^2$, $b (\mu) \geq \frac{N - s}{q}$. Thus
\[
\big[N -s + (1 - \frac{N}{p})q \big] \frac{p - q}{q} < b(\mu) p + p - N -
\big[N  -s + (1 - \frac{N}{p})q \big].
\]
 Choose $\lambda = \epsilon^\tau$, where $\big[N - s + (1 - \frac{N}{p})q \big]
\frac{p - q}{q} < \tau < b(\mu) p + p - N - \big[N - s + (1 -
\frac{N}{p})q \big]$. Then
\[
\lambda O (\epsilon^{N - s + (1 - \frac{N}{p})q}) = O (\epsilon^{\tau
+ N - s + (1 - \frac{N}{p})q}), \quad
d \lambda^\frac{p}{p - q} = O (\epsilon^\frac{p \tau}{p - q}).
\]
Since $\tau + N - s + (1 - \frac{N}{p})q < \frac{p \tau}{p - q}$,
$\tau + N - s + (1 - \frac{N}{p})q < b(\mu) p + p - N$, taking
$\epsilon$ small enough we deduce that there exists $\delta > 0$,
such that
\begin{equation} \label{san4}
O (\epsilon ^{b(\mu) p^*(s_1) - N + s_1}) - \lambda O (\epsilon^{N
- s + (1 - \frac{N}{p})q}) < - d \lambda^\frac{p}{p - q}, \quad
\forall  \lambda : 0 < \lambda^\frac{p}{p - q} < \delta.
 \end{equation}
Choose $\Lambda_1 = \min \{\Lambda_1^*, \frac{p - q}{p}
 \delta\}$. Then for all $\lambda \in (0, \Lambda_1)$ we have
\[
\sup_{t \geq t_1} J (t u_1, t v_1)
\leq   \frac{p - t}{p(N -t)} \frac{1}{\|Q\|_\infty^\frac{N - p}{p - t}}
\big(S_{t, \alpha,\beta}  \big)^\frac{N - t}{p - t}  - d \lambda^\frac{p}{p - q}.
\]
Together with \eqref{sa1}, we get the conclusion of Lemma \ref{lem11}.
\end{proof}

\section{Proof of the main results}

\begin{proof}[Proof of Theorem \ref{thmA}]
Let 
\begin{gather*}
r:  =  \|(u,v)\|\,,
\\
f(r):  =  \frac{1}{p}r^p-\frac{1}{p^*(s_1)}
 A_{\mu,s_1}^{-\frac{p^*(s_1)}{p}}r^{p^*(s_1)}
-\frac{1}{p^*(s_2)}A_{\mu,s_2}^{-\frac{p^*(s_2)}{p}}r^{p^*(s_2)} -
\frac{1}{p^* (t)}
S_{t, \alpha, \beta}^{-\frac{p^* (t)}{p}} \|Q\|_\infty\,,
\\
h(r):  =  \frac{\lambda}{q} \Big(\frac{N \omega_N R_0^{N - s}}{N- s}
 \Big)^\frac{p^* (s) - q}{p^* (s)} A_{\mu, s}^{-\frac{q}{p}} r^q\,.
\end{gather*}
From \eqref{san}, \eqref{hardy2} and \eqref{hold},
\[
J (u,v)  \geq f(r)-h(r).
\]
Note that $p < p^*(s_1)$, $p^*(s_2)$, $p^*(t)$, it is easy to
see that there exists $\varrho > 0$ such that $f(r)$ achieves its
maximum at $\varrho$ and $f(\varrho)>0$. Therefore, there exists
$\Lambda_{11}>0$, such that for $0 < \lambda <\Lambda_{11}$,
\begin{equation}\label{e3.1}
\inf_{\|(u,v)\|= \varrho} I(u,v) \geq f(\varrho) - h(\varrho) > 0.
\end{equation}
On the other hand, set $B_\varrho=\{(u,v);\;\|(u,v)\|\leq
\varrho \}$. For $(u,v) \neq  (0,0)$, we can choose $d > 0$
small enough, such that $(d u,d v)\in B_\varrho$ and
\begin{equation}\label{e3.2}
I(d u,d v) < 0.
\end{equation}
Thus,
\begin{equation}\label{e3.3}
-\infty < \inf_{(u,v)\in B_\varrho} I(u,v) < 0.
\end{equation}
Now we can apply the Ekeland variational principle in \cite{MW}
and obtain $\{(u_n,v_n)\}\subset B_\varrho$, such that
\begin{gather}\label{e3.4}
I(u_n,v_n)\leq\inf_{(u,v)\in B_\varrho}I(u,v)+\frac{1}{n},\\
\label{e3.5}
I(u_n,v_n)\leq I(u,v)+\frac{1}{n}\|(u_n-u,v_n-v)\|,
\end{gather}
for all $(u,v)\in B_R$. Define
\begin{equation}\label{e3.6}
J_1(u,v):= J (u,v)+\frac{1}{n}\|(u_n-u,v_n-v)\|.
\end{equation}
By \eqref{e3.5}, we have $(u_n,v_n)$ is the minimizer of
$J_1(u,v)$ on $B_\varrho$. \eqref{e3.1}, \eqref{e3.3} and
\eqref{e3.4} imply that there exists $\epsilon>0$ and $k\in N$,
such that for $n\geq k$, 
$\{(u,v),\:\|(u,v)\|\leq \varrho -\epsilon\}$. Therefore, for $n\geq k$ and
$(\phi,\varphi)\in W$, we can choose $t>0$ small enough, such that
$(u_n+t\phi,v_n+t\varphi)\in B_\varrho$ and
\[
\frac{J_1 (u_n+t\phi,v_n+t\varphi)- J_1 (u_n,v_n)}{t}\geq 0.
\]
That is,
\begin{equation}\label{e3.7}
\frac{J (u_n+t\phi,v_n+t\varphi)- J (u_n,v_n)}{t} +
\frac{1}{n}\|(\phi,\varphi)\|\geq 0.
\end{equation}
Passing to the limit in \eqref{e3.7} as $n\to0$, one can obtain
\[
\langle J'(u_n,v_n),(\phi,\varphi)\rangle\geq-\frac{1}{n}\|(\phi,\varphi)\|,
\]
which implies 
\begin{equation}\label{e3.8}
\|J'(u_n,v_n)\|\leq\frac{1}{n}.
\end{equation}
Combining \eqref{e3.4} and \eqref{e3.8}, there holds
\begin{gather}\label{e3.9}
\lim_{n\to\infty} J (u_n,v_n)= \inf_{(u,v)\in B_\varrho}J(u,v) < 0,\\
\label{e3.10}
\lim_{n\to\infty}J'(u_n,v_n)=0.
\end{gather}
We note that there exists $\Lambda_{11}^* \in (0, \Lambda_{11})$,
such that for $0<\lambda < \Lambda_{11}^*$, and
$c_0>\inf_{(u,v)\in B_\varrho } I(u,v)$, where $c_0$ is defined in
Lemma \ref{ps3}. Thus, \eqref{e3.9} and \eqref{e3.10} and Lemma
Lemma \ref{ps3} imply that for $0<\lambda <
\Lambda_{11}^*,\;(u_n,v_n)\to(u,v)$ strongly in $W$.
Therefore, $(u,v)$ is a nontrivial solution of problem \eqref{1}
satisfying $J(u,v)=\inf_{(u,v)\in B_\varrho} J (u,v) < 0$. Note
that $J (u,v)=J (|u|,|v|)$ and
$$
(|u|,|v|)\in \{(u,v),\:\|(u,v)\|\leq \varrho - \epsilon \},
$$
we have $I(|u|,|v|)=\inf_{(u,v)\in B_\varrho} J(u,v)$ and
$J'(|u|,|v|)=0$. Then problem \eqref{1} has a nontrivial
nonnegative solution. By the strongly maximum principle, we get
the conclusion of Theorem \ref{thmA}. 
\end{proof}

\begin{proof}[Proof of Theorem \ref{the1}]
 In view of the proof of
Theorem \ref{thmA}, we know that for $0<\lambda< \Lambda_{11}$,
there exists $\varrho > 0$, such that $\inf_{\|(u,v)\|=\varrho}
I(u,v) \geq \vartheta^* > 0$. Moreover, \eqref{e3.9} and
\eqref{e3.10} hold. We note that there exists $\Lambda_{12} \in
(0,\Lambda_{11})$, such that for $0< \lambda < \Lambda_{12},\;c_*>
\inf_{(u,v) \in B_\varrho} J(u,v)$, where $c_*$ is defined in
Lemma \ref{ps2}. Thus \eqref{e3.9} and \eqref{e3.10} and Lemma
\ref{ps2} imply that $(u_n,v_n) \to (u,v)$ strongly in
$W$. Standard argument shows that for $0 < \lambda< \Lambda_{12}$,
problem \eqref{1} has at least one positive solution satisfying
 $J(u,v)<0$ and $J'(u,v)=0$.

Now we prove a second positive solution. It is easy to see $J
(0,0)=0$. Set $\Lambda^{**}= \min \{\Lambda_{12}, \Lambda_1\}$,
where $\Lambda_1$ is given in Lemma \ref{lem11}. Then it follows
from Lemma \ref{lem11} there exists $(u',v')\in W \setminus \{0
\}$, such that for $0< \lambda < \Lambda^{**}$,
$$
\sup_{t \geq 0} J (tu',tv')< c_*.
$$
On the other hand we obtain that $\lim_{l\to\infty} J(lu',lv')=-\infty$.
Thus there exists $l'>0$ such that $\|(l'u',l'v')\|> \varrho$ 
and $J(l'u',l'v')<0$. Let
$$
c:=\inf_{\gamma\in\Gamma}\sup_{t\in[0,1]}J(\gamma(t)),
$$
where
$$
\Gamma:=\{\gamma\in C^0([0,1],W):\gamma(0)=(0,0),\,\gamma(1)=(l'u',l'v')\}.
$$
Thus, it follows from the mountain pass lemma in 
\cite{A. Ambrosetti} that there exists a sequence $(u_n,v_n) \in W$ such
that
\begin{gather*}
\lim_{n\to\infty}J(u_n,v_n)=c,\\
\lim_{n\to\infty} J'(u_n,v_n)=0.
\end{gather*}
Moreover, $c\in(0,c_*)$. From Lemma \ref{ps2},
$(u_n,v_n)\to(\overline{u},\overline{v})$ strongly in $W$,
 which implies that $J (\overline{u},\overline{v})=c$ and 
$J'(\overline{u},\overline{v})=0$, Therefore,
$(\overline{u},\overline{v})$ is a second nontrivial solution of \eqref{1}.
Set $u^+=\max\{u,0\},v^+=\max\{v,0\}$. Replacing
\[
\int_{\Omega} \frac{|u|^q}{|x|^s} dx,\quad
\int_{\Omega} \frac{|v|^q}{|x|^s}  dx,\quad
\int_{\Omega}\frac{|u|^{p^*(s_1)}}{|x|^{s_1}}dx,\quad
\int_{\Omega}\frac{|v|^{p^*(s_2)}}{|x|^{s_2}}dx, \quad
\int_{\Omega} Q(x) \frac{|u|^\alpha |v|^\beta}{|x - x_0|^t}dx
\]
by
\begin{gather*}
\int_{\Omega} \frac{(u^+)^q}{|x|^s} dx,\quad
\int_{\Omega} \frac{(v^+)^q}{|x|^s} dx,\quad
\int_{\Omega}\frac{(u^+)^{p^*(s_1)}}{|x|^{s_1}}dx,\\
\int_{\Omega}\frac{(v^+)^{p^*(s_2)}}{|x|^{s_2}}dx, \quad
\int_{\Omega} Q(x) \frac{(u^+)^\alpha (v^+)^\beta}{|x - x_0|^t}dx
\end{gather*}
and repeating the above process, we have a nonnegative solution 
$(\widetilde{u},\widetilde{v})$ of problem \eqref{1} satisfying 
$J (\widetilde{u},\widetilde{v}) > 0$.
 Then by the strongly maximum principle, we have a second positive solution.
\end{proof}

\subsection*{Acknowledgments}
 The authors would like to thank the
anonymous referees for his/her valuable suggestions and comments.

\begin{thebibliography}{99}

\bibitem{C. Alves} C. Alves, D. Filho, M.\ Souto;
\emph{On systems of elliptic equations involving subcritical or
critical Sobolev exponents}, Nonlinear Anal. 42 (2000) 771-787.

\bibitem{A. Ambrosetti} A. Ambrosetti, P. H. Rabinowitz;
\emph{Dual variational methods in critical point theory and
applications}, J. Funct. Anal. 14 (1973) 349-381.

\bibitem{L. Caffarelli}  L. Caffarelli, R. Kohn, L. Nirenberg;
\emph{First order interpolation inequality with weights}, Compos.
Math. 53 (1984) 259-275.

\bibitem{D. Cao} D. Cao, P. Han;
\emph{Solutions for semilinear elliptic equations with critical
exponents and Hardy potential}, J. Differential Equations 205
(2004) 521-537.

\bibitem{D. Cao1} D. Cao, D. Kang;
\emph{Solutions of quasilinear elliptic problems involving a
Sobolev exponent and multiple Hardy-type terms}, J. Math. Anal.
Appl. 333 (2007) 889-903.

\bibitem{Degiovanni} M. Degiovanni, S. Lancelotti;
\emph{Linking solutions for p-Laplace equations with nonlinearity
at critical growth}, J. Funct. Anal. 256 (2009) 3643-3659.

\bibitem{F. Demengel} F. Demengel, E. Hebey;
\emph{On some nonlinear equations involving the p-Laplacian with
critical Sobolev growth}, Adv. Differential Equations 3 (1998)
533-574.

\bibitem{R. Filippucci} R. Filippucci, P. Pucci, F. Robert;
\emph{On a p-Laplace equation with multiple critical
nonlinearities}, J. Math. Pures Appl. 91 (2009) 156-177.

\bibitem{J. Garcia Azorero} J. Garcia Azorero, I. Peral;
\emph{Hardy inequalities and some critical elliptic and parabolic
problems}, J. Differential Equations 144 (1998) 441-476.

\bibitem{N. Ghoussoub} N. Ghoussoub, F. Robert;
\emph{The effect of curvature on the best constant in the
Hardy–Sobolev inequality}, Geom. Funct. Anal. 16 (2006) 897-908.

\bibitem{N. Ghoussoub1} N. Ghoussoub, C. Yuan;
\emph{Multiple solutions for quasi-linear PDEs involving the
critical Sobolev and Hardy exponents}, Amer. Math. Soc. 352 (2000)
5703-5743.

\bibitem{N. Ghoussoub2} N. Ghoussoub, C. Yuan;
\emph{Multiple solutions for quasi-linear PDEs involving the
critical Sobolev and Hardy exponents}, Trans. Amer. Math. Soc. 352
(2000) 5703-5743.

\bibitem{Han1} P. Han;
\emph{Multiple positive solutions for a critical growth problem
with Hardy potential}, Proc. Edinb. Math. Soc. 49 (2006) 53-69.

\bibitem{Han2} P. Han;
\emph{Multiple solutions to singular critical elliptic equations},
Israel J. Math. 156 (2006) 359-380.

\bibitem{Han4} P. Han;
\emph{Many solutions for elliptic equations with critical
exponents}, Israel J. Math. 164 (2008) 125-152.


\bibitem{P. Han} P. Han;
\emph{Quasilinear elliptic problems with critical exponents and
Hardy terms}, Nonlinear Anal. 61 (2005) 735-758.

\bibitem{Han3} P. Han, L. Zhaoxia;
\emph{Solutions for a singular critical growth problem with a
weight}, J. Math. Anal. Appl. 327 (2007) 1075-1085.

\bibitem{Kang} D. Kang;
\emph{On the quasilinear elliptic problem with a critical
Hardy–Sobolev exponent and a Hardy term}, Nonlinear Anal. 69
(2008) 2432-2444.

\bibitem{S. Peng} D. Kang, S. Peng;
\emph{Solutions for semilinear elliptic problems with critical
Sobolev–Hardy exponents and Hardy potential}, Appl. Math. Lett. 18
(2005) 1094-1100.

\bibitem{Y. Li} Y. Li, Q. Guo, P. Niu;
\emph{The existence of solutions for quasilinear elliptic problems
with combined critical Sobolev-Hardy terms}, J. Math. Anal. Appl.
388 (2012) 525-538.

\bibitem{P.L. Lions} P. L. Lions;
\emph{The concentration-compactness principle in the calculus of
variations. the limit case, part 1}, Rev. Mat. Iberoam. 1 (1985)
145-201.

\bibitem{MW} J. Mawhin, M. Willem;
\emph{Critical Point Theory and Hamiltonian Systems}, Appl. Math.
Sci. 74, Springer, New York (1989).

\bibitem{Nia1} N. Nyamoradi;
\emph{Existence and multiplicity of solutions to a singular
elliptic system with critical Sobolev-Hardy exponents and
concave-convex nonlinearities}, J. Math. Anal. Appl. 396 (2012)
280-293.

\bibitem{P. Rabinowitz} P. Rabinowitz;
\emph{Minimax methods in critical points theory with applications
to differential equations}, in: CBMS Regional Conference Series in
Mathematics, vol. 65, Providence, RI, 1986.

\bibitem{Shang} Y. Y. Shang;
\emph{Existence and multiplicity of positive solutions for some
Hardy-Sobolev critical elliptic equation with boundary
singularities}, Nonlinear Anal. 75 (2012) 2724-2734.

\bibitem{Xuan1} B. Xuan, J. Wang'
\emph{Existence of a nontrivial weak solution to quasilinear
elliptic equations with singular weights and multiple critical
exponents}, Nonlinear Analysis 72 (2010) 3649-3658.

\bibitem{J. Zhang} J. Zhang, Z. Wei;
\emph{Existence of multiple positive solutions to singular
elliptic systems involving critical exponents}, Nonlinear Anal. 75
(2012) 559-573.

\end{thebibliography}

\end{document}