In this paper, we consider the following elliptic problem with singular coefficient

$\begin{equation}\label{1.1} \left\{ \begin{array}{ll} -\mbox{div}(|x|^{-2a}\nabla u)-\mu \frac{u}{|x|^{2(1+a)}}=\frac{1}{|x|^{bp}}(|u|^{p-2}u+\frac{\eta \alpha}{\alpha+\beta}|u|^{\alpha-2}u|v|^\beta)+\frac{a_1u+a_2v}{|x|^{dD}}, &x\in \Omega, \\ -\mbox{div}(|x|^{-2a}\nabla v)-\mu \frac{v}{|x|^{2(1+a)}}=\frac{1}{|x|^{bp}}(|v|^{p-2}v+\frac{\eta \beta}{\alpha+\beta}|u|^{\alpha}|v|^{\beta-2}v)+\frac{a_2u+a_3v}{|x|^{dD}}, &x\in\Omega, \\ u =v =0, &x\in \partial \Omega, \\ \end{array} \right. \end{equation}$

(1.1)

where $\Omega $ is a smooth bounded domain in ${{\mathbb{R}}^{N}}(N\ge 3),0\in \Omega ,\eta \ge 0,$${{a}_{i}}\in \mathbb{R},i=1,2,3,0\le \mu ＜{{(\sqrt{{\bar{\mu }}}-a)}^{2}},$$\bar{\mu }=:{{(\frac{N-2}{2})}^{2}},0\le a＜\sqrt{{\bar{\mu }}},a\le b＜a+1,$$a\le d＜a+1,\alpha ,\beta 1,\alpha +\beta =p=p(a,b)$$a\le d＜a+1,\alpha ,\beta＞1,\alpha +\beta =p=p(a,b)$$=:\frac{2N}{N-2(1+a-b)},D=D(a,d)=:\frac{2N}{N-2(1+a-d)}$. For problem (1.1), we are interested in the existence and non-existence of a nontrivial solution $(u, v)$, that is to say that $u\not\equiv0$ and $v\not\equiv0$. Moreover, we call a solution $(u, v)$ semi-trivial if $(u, v)$ is type of $(u, 0)$ or $(0, v)$.

Problem (1.1) can be seen as a counterpart of the following elliptic equation

$\label{1.2} \left\{ \begin{array}{ll} { -\mbox{div}(|x|^{-2a}\nabla u)-\mu \frac{u}{|x|^{2(1+a)}}=\frac{|u|^{p-2}u}{|x|^{bp}}+\lambda \frac{u}{|x|^{dD}}}, &x\in \Omega, \\ u = 0 , &x\in \partial \Omega. \end{array} \right. $

(1.2)

In particular, when $a=b=d=\mu=0$, problem (1.2) reduces to the Brezis-Nirenberg problem

$\begin{equation}\label{1.3} \left\{ \begin{array}{ll}{ -\Delta u=|u|^{2^*-2} u+\lambda u}, &x\in \Omega, \\ \hspace{0.5em}u = 0 , &x\in \partial \Omega. \end{array} \right. \end{equation}$

(1.3)

In the well-known literature [5], Brezis and Nirenberg proved the existence of positive solutions to (1.3), when $0<\lambda<\lambda_1(\Omega)$, $ N\geq 4$ and $\lambda^*<\lambda<\lambda_1(\Omega)$, $N=3$, where $\lambda_1(\Omega)$ is the first eigenvalue of $-\Delta$ on $\Omega$ with Dirichlet boundary condition and $\lambda^*\in (0, \lambda_1(\Omega))$. Moreover, in [11, 13, 28, 29], sign-changing solutions to (1.3) were obtained. For (1.2), when $a=b=d=0, \mu\neq 0$, i.e.,

$\begin{equation}\label{1.4} \left\{ \begin{array}{ll}{ -\Delta u-\mu \frac{u}{|x|^2}=|u|^{2^*-2} u+\lambda u}, &x\in \Omega, \\ \hspace{0.5em}u = 0 , &x\in \partial \Omega. \end{array} \right. \end{equation}$

(1.4)

For (1.4), Enrico Jannelli in [20] studied the role of space dimension on the existence of solutions, on one hand, the existence of positive solutions was obtained when $\mu\leq \bar{\mu}-1, 0<\lambda<\overline{\lambda}_1(\mu)$ and $\bar{\mu}-1<\mu<\bar{\mu}$, $\lambda_*(\mu)<\lambda<\overline{\lambda}_1(\mu)$; on the other hand, the non-existence of positive solutions was also proved in the case $\bar{\mu}-1<\mu<\bar{\mu}, \lambda\leq\lambda_*(\mu)$ and $\Omega=B_R(0)$, where $\gamma=\sqrt{\bar{\mu}}+\sqrt{\bar{\mu}-\mu}$,

$ \begin{array}{ll} &\overline{\lambda}_1(\mu)=\inf\limits_{u\in H_0^1(\Omega)\setminus\{0\}}\int_\Omega\bigg(|\nabla u |^2-\mu\frac{u^2}{|x|^{2}}\bigg)\Bigg/\int_\Omega u^2, \\ &\lambda_*(\mu)=\inf\limits_{v\in H_0^1(\Omega, |x|^{-2a})\setminus\{0\} }\int_\Omega\frac{|\nabla v |^2}{|x|^{2\gamma}}\Bigg/\int_\Omega\frac{v^2}{|x|^{2\gamma}}. \end{array} $

Meanwhile, in [10, 25], sign-changing solutions were proved to exist when $N\ge 7$, $\mu \in [0, \ \bar \mu-4)$ and $\lambda \in (0, \overline{\lambda}_1(\mu))$. While for the nonexistence result, it was proved in \cite{DW} that (1.4) has no radial sign-changing solutions for $\lambda\in (0, \lambda(N))$ when $3\le N\le 6$, $\Omega=B_1(0)$, where $\lambda(N)>0$ depending on $N$.

For (1.2), it is clear that singularity occurs, the singularity of potential $\frac{\mu}{|x|^{2(1+a)}}$ is critical both from the mathematical and the physical point of view. As it does not belong to the Kato's class, it cannot be regarded as a lower order perturbation of the laplacian but strongly influences the properties of the associated elliptic operator. To be mentioned, singular potentials arise in many fields, such as quantum mechanics, nuclear physics, molecular physics, and quantum cosmology, we refer to [18] for further discussion and motivation.

Mathematically, (1.2) is related to the following well-known Caffarelli-Kohn-Nirenberg inequality (see [9])

$\begin{equation}\label{1.5} \Bigl(\int_{\mathbb{R}^N}|x|^{-bp} |u|^{p}\Bigl)^{\frac{2}{p}}\leq C_{a, b}\int_{\mathbb{R}^N}|x|^{-2a}|\nabla u|^2,\quad \forall u\in C_0^\infty(\mathbb{R}^N), \end{equation}$

(1.5)

where $0\le a<\sqrt{\bar{\mu }}, a \le b<a+1, p=p(a, b)=:\frac{2N}{N-2(1+a-b)}$. $p=p(a, b)$ is called the critical Sobolev-Hardy exponent, since (1.5) is classical Sobolev and Hardy inequality respectively in the case $a=b=0$ and $a=0, b=1$. (1.5) plays an important role in many applications by virtue of the complete knowledge about the best constant $C_{a, b}$ and the extremal functions (see [9, 12]). Concerning (1.2), the existence and non-existence of sign-changing solutions were studied in [23] and [24]. For some other related results, we refer to [1-4, 7, 8, 10, 12, 15-17, 20, 25, 27] and the references therein.

Based on these results, a nature problem is: can we obtain the existence of positive solution and sign-changing solution for system (1.1)? In this paper, we will investigate the above problems and we obtain an affirmative answer.

To state our main results, we need to introduce some notations.

For $\mu \in [0, (\sqrt{\bar{\mu }}-a)^2), $ define $H:=H_0^1(\Omega, |x|^{-2a}) $ to be the completion of $C_0^\infty(\Omega )$ with respect to the norm

$\begin{equation}\label{1.51} \|u\|=\|u\|_{H}=:\bigg(\int_\Omega \Bigr(|x|^{-2a}|\nabla u|^2 - \mu \frac{|u|^2}{|x|^{2(1+a)}}\Bigr)\bigg)^\frac{1}{2}. \end{equation}$

(1.6)

Set $b=a+1$ in (1.5), we have the following weighted Hardy inequality (see [7, 12])

$\begin{equation}\label{1.6} \int_{\mathbb{R}^N}\frac{|u|^2}{|x|^{2(1+a)}} \le \frac{1}{(\sqrt{\bar{\mu }}-a)^2}\int_{\mathbb{R}^N}|x|^{-2a}|\nabla u|^2, \hspace{.5em}\forall u\in C_0^\infty (\mathbb{R}^N). \end{equation}$

(1.7)

Hence norm (1.6) is well defined and equivalent to the usual norm $(\displaystyle\int_\Omega |x|^{-2a}|\nabla u|^2 dx)^\frac{1}{2}$.

Denote $W:=H\times H$ to be the completion of $C_0^\infty(\Omega)$ with respect to the norm $\|(u, v)\|^2:=\|u\|^2+\|v\|^2.$

Define the energy functional corresponding to problem (1.1)

$\begin{equation}\label{1.7} J(u, v)=\frac{1}{2}\|(u, v)\|^2-\frac12\int_\Omega\frac{A(u, v)}{|x|^{dD}}-\frac{1}{p}\int_{\Omega} \frac{F(u, v)}{|x|^{bp}}, \forall (u, v) \in W, \end{equation}$

(1.8)

where $A(u, v):=a_1u^2+2a_2uv+a_3v^2, F(u, v):=|u|^p+|v|^p+\eta|u|^\alpha|v|^\beta$. Then $J\in C^1(W, \mathbb{R})$. The duality product between $W$ and its dual space $W^{-1}$ is defined as

$\begin{eqnarray*} \langle J'(u, v), (\varphi, \psi)\rangle&=& \int_\Omega(\frac{\nabla u\nabla\varphi+\nabla v\nabla\psi}{|x|^{2a}}-\mu\frac{u\varphi+v\psi}{|x|^{2(1+a)}})\\ &&- \int_{\Omega}\frac{a_1u\varphi+a_2v\varphi+a_2u\psi+a_3v\psi}{|x|^{dD}}\\ &&- \int_\Omega\frac{|u|^{p-2}u\varphi+|v|^{p-2}v\psi+\frac{\eta \alpha}{\alpha+\beta}|u|^{\alpha-2}|v|^\beta u\varphi+\frac{\eta \beta}{\alpha+\beta}|u|^\alpha|v|^{\beta-2} v\psi}{|x|^{bp}}, \end{eqnarray*}$

where $u, v, \varphi, \psi\in H$. A pair of functions $(u, v)\in W$ is said to be a solution of problem (1.1) if

$ \langle J'(u, v), (\varphi, \psi)\rangle=0, \forall (\varphi, \psi)\in W. $

Define $\beta^*=\sqrt{(\sqrt{\bar{\mu}}-a)^2-\mu}$, $\nu=\sqrt{\bar{\mu}}-a-\beta^*$, $\gamma=\sqrt{\bar{\mu}}-a+\beta^*$. Set operator $L(\cdot)= -\mbox{div}(|x|^{-2a}\nabla\cdot)-\frac{\mu}{|x|^{2(1+a)}}\cdot$. Define $\lambda^*_1(\mu)$ to be the first eigenvalue of problem

$\begin{equation}\label{1.8} L(u)=\lambda\frac{u}{|x|^{dD}}, u\in H \end{equation}$

(1.9)

and $\lambda_1(\mu)$ the first eigenvalue of problem

$\begin{equation}\label{1.9} -\mbox{div}(|x|^{-2(a+\gamma)}\nabla \varphi)=\lambda \frac{\varphi}{|x|^{2\gamma+dD}}, \varphi\in H_0^1(\Omega, |x|^{-2(a+\gamma)}), \end{equation}$

(1.10)

that is

$ \begin{array}{ll} &\lambda^*_1(\mu)=\inf\limits_{u\in H\setminus\{0\}}\int_\Omega\bigg(\frac{|\nabla u |^2}{|x|^{2a}}-\mu\frac{u^2}{|x|^{2(1+a)}}\bigg)\Bigg/\int_\Omega\frac{u^2}{|x|^{dD}}, \\ &\lambda_1(\mu)=\inf\limits_{v\in H^1_0(\Omega, |x|^{-2(a+\gamma)})\setminus\{0\} }\int_\Omega\frac{|\nabla v |^2}{|x|^{2(a+\gamma)}}\Bigg/\int_\Omega\frac{v^2}{|x|^{2\gamma+dD}}. \end{array} $

By Sobolev inequality and Young inequality, the following best constants are well defined

$S(\mu):=\inf_{u\in H\setminus \{0\}}\frac{\|u\|^2}{\Bigr(\int_{\mathbb{R}^N}\frac{|u|^p}{|x|^{bp}}\Bigr)^{\frac{2}{p}}}, $

(1.11)

${S_{\eta ,\alpha ,\beta }}(\mu ): = \mathop {\inf }\limits_{(u,v) \in W \setminus \{ (0,0)\} } \frac{{||(u,v)||{^2}}}{{{{(\int_{{^N}} {\frac{{|u{|^p} + |v{|^p} + \eta |u{|^\alpha }|v{|^\beta }}}{{|x{|^{bp}}}}} )}^{\frac{2}{p}}}}}.$

(1.12)

Set

$f(\tau ): = \frac{{1 + {\tau ^2}}}{{{{(1 + \eta {\tau ^\beta } + {\tau ^p})}^{\frac{2}{p}}}}},f({\tau _{\min }}): = \mathop {\min }\limits_{\tau \ge 0} f(\tau ) \le 1.$

Throughout this paper, we always assume that the following conditions:

(H1)$0\le a<\sqrt{\bar{\mu }}, 0 \le \mu ＜ (\sqrt{\bar{\mu }}-a)^2, \eta\ge 0, \alpha, \beta>1, \alpha+\beta=p.$

(H2)$a_i\geq0, i=1, 2, 3, a_1a_3-a_2^2>0, 0<\Lambda_1\leq \Lambda_2<\lambda_1^*(\mu), $ where $\Lambda_1$ and $\Lambda_2$ are the eigenvalues of the matrix

$ A:=\left ( \begin{array}{ll} a_1 a_2\\ a_2 a_3 \end{array} \right ). $

Our main results are as follows:

Theorem 1.1  Suppose $N \ge 4+4a-dD$ and (H1), (H2).

(ⅰ) If $0 \le \mu \le (\sqrt{\bar{\mu}}-a)^2-(1+a-\frac{dD}{2})^2$, then $(1.1)$ has a positive solution in $W$ when $ 0 < \Lambda_1\leq\Lambda_2 <\lambda^*_1(\mu)$.

(ⅱ) If $(\sqrt{\bar{\mu}}-a)^2-(1+a-\frac{dD}{2})^2<\mu<(\sqrt{\bar{\mu}}-a)^2$, then $(1.1)$ has a positive solution in $W$ when $\lambda_1(\mu)<\Lambda_1\leq \Lambda_2<\lambda^*_1(\mu)$.

Theorem 1.2  Suppose (H1), (H2), $\eta=0$, $N\geq \max\{6(1+a)-2dD, 4+2a\}, \ 0\leq \mu <(\sqrt{\bar{\mu}}-a)^2-(\max\{2(1+a)-dD, 1\})^2$ and $\Lambda_1, \Lambda_2\in(0, \ \lambda^*_{1}(\mu))$, then $(1.1)$ has a pair of sign-changing solutions.

Remark 1.3  Theorem 1.2 says that when $a=b=d=0$, $(1.1)$ has a pair of sign-changing solutions. This result generalizes the results of Theorem 1.3 (ⅰ) in [19].

To verify Theorem 1.1, we mainly employ the framework in [5, 20]. However, the singularity of the solutions and the non-uniform ellipticity of the operator $-\mbox{div}(|x|^{-2a}\nabla\cdot)$ bring us more difficulties, so we need to find new arguments. On one hand, to obtain positive solutions, a new maximum principle should be established; on the other hand, we need to estimate the asymptotic behavior (near the origin) of (1.2). Moreover, whether or not $\lambda_1(\mu)$ can be attained is not clear and we also need to estimate $\lambda_1(\mu)$ and $\lambda^*_1(\mu)$.

To obtain Theorem 1.2, our methods are inspired by the work of [19]. However, comparing with [19], since the generality of (1.1), more complex calculation will be needed.

This paper is organized as follows. In Section 2, we will give some important preliminaries. A positive solution will be obtained in Section 3 by using the mountain pass lemma. In the last section, we will discuss the existence of sign-changing solutions. In this paper, for simplicity, we denote $C$ (may be different in different places) positive constants, $B_r(x):=\{y\in \mathbb{R}^N:|y-x|<r\}$ and we omit $dx$ in the integral.

In this section, we shall give some preliminaries and a non-existence result.

Lemma 2.1  Suppose $0\le a<\sqrt{\bar{\mu }},a\le b<a+1$ and $0 \le \mu < (\sqrt{\bar{\mu }}-a)^2.$ Then

(ⅰ) $S(\mu)$ is independent of $\Omega$.

(ⅱ) When $\Omega = \mathbb{R}^N$, $S(\mu)$ can be achieved by the functions

$ U_\varepsilon (x) =\Bigr(2\varepsilon p\beta ^2 \Bigr)^{\frac{1}{p-2}}\Bigr/\Bigr(|x|^\nu \bigr(\varepsilon + |x|^{(p-2)\beta } \bigr)^{\frac{2}{p-2}} \Bigr) $

for all $\varepsilon >0$. The functions $U_\varepsilon(x) $ solve the equation

$\begin{equation}\label{2.1} -\mbox{div}(|x|^{-2a}\nabla u) -\mu \frac{u}{|x|^{2(1+a)}} = \frac{|u|^ {p-2}u}{|x|^{bp}}, x\in \mathbb{R}^N \setminus \{0\} \end{equation}$

(2.1)

satisfying

$ \int_{\mathbb{R}^N}\Bigr( |x|^{-2a}|\nabla U_\varepsilon |^ 2-\mu \frac{U_\varepsilon^ 2}{|x|^{2(1+a)}} \Bigr) = \int_{\mathbb{R}^N}\frac{U_\varepsilon^{p}}{|x|^{bp}} = S(\mu)^{\frac{p}{p-2}}. $

Proof  The result was proved in [7, 12].

Lemma 2.2  Suppose (H1) and (H2), then

(ⅰ) $S_{\eta, \alpha, \beta}(\mu)=f(\tau_{\min})S(\mu).$

(ⅱ) $S_{\eta, \alpha, \beta}(\mu)$ has the minimizers $(U_\varepsilon (x), \tau_{\min}U_\varepsilon (x)), \forall \varepsilon>0$, where $f(\tau):=\frac{1+\tau^2}{\Bigr(1+\eta\tau^\beta+\tau^p\Bigr)^{\frac{2}{p}}}$, $f(\tau_{\min}):=\min_{\tau\geq 0}f(\tau)\leq 1$ and $\tau_{\min}$ satisfies

$ p+\eta\alpha \tau^\beta-\eta\beta \tau^{\beta-2}-p\tau^{p-2}=0. $

Proof  The proof is similar to Theorem 1.1 in [19]. Here we omit it.

Lemma 2.3  Let $\tau>2-N$. Suppose that $u\in C^2(\Omega\setminus\{0\}), u\geq0, u\not\equiv 0$ satisfies $ -\mbox{div}(|x|^\tau\nabla u)\geq 0, $ then $u>0$ in $\Omega\setminus\{0\}$.

Proof   The proof is similar to [6] or [7]. Here we omit it.

Lemma 2.4  Suppose that (H1), (H2) and $(u_1(x), v_1(x))\in W$ is a positive solution of (1.1), then

(ⅰ) if $0 \le \mu <(\sqrt{\bar{\mu }}-a)^2$, then for any $B_\rho(0)\subset\Omega$, there exist $0<C_1<C_2<\infty$ such that

$\begin{equation}\label{2.2} c_1|x|^{-\nu}\leq u_1(x), v_1(x)\leq c_2|x|^{-\nu}, \forall x\in B_\rho(0)\setminus\{0\}. \end{equation}$

(2.2)

(ⅱ) $0\leq\lambda_1(\mu)<\lambda^*_1(\mu)$.

Proof  The proof is similar to [6] and [20]. Here we omit it.

To complete this section, we give a nonexistence result of solutions for (1.1).

Lemma 2.5  If $\Omega$ is star-shaped with respect to the origin and $\Lambda_2\leq 0$, then (1.1) has no solution in $W$.

Proof  The proof is based on a Pohozaev's type identity which can be verified by the similar method as [7].

$\begin{equation}\label{2.3} \frac{1}{2}\int_{\partial\Omega}|x|^{-2a}(|\nabla u|^2+|\nabla v|^2) (\vec{n}\cdot x)=\Bigl(\frac{N-dD}{2}-\frac{N-2-2a}{2}\Bigl)\int_{\Omega}\frac{a_1u^2+2a_2uv+a_3v^2}{|x|^{dD}}. \end{equation}$

(2.3)

But $\frac{N-dD}{2}-\frac{N-2-2a}{2}>0$ by our assumptions, hence (2.3) is impossible in the case $\Lambda_2\leq 0$ since the left hand side of (2.3) is positive. So we complete our proof.

From Lemma 2.5, to obtain positive solution of (1.1), we impose the condition $\Lambda_1, \Lambda_2>0$.

In this section, we will prove Theorem 1.1. Since $J\in C^2(W, \mathbb{R})$, we see that critical points of functional $J$ correspond to the weak solution of (1.1).

Define $ \mathcal{D}=\{\varphi\in C_0^\infty(\Omega): \varphi\equiv 1 \hbox{in a neightbourhood of } x=0\}. $ Let $\varphi(x)\in \mathcal{D}$ and set $u_\varepsilon(x)$ = $\varphi (x) U_\varepsilon (x)$, then $u_\varepsilon(x)\in H$. By direct calculation, we have the following estimates: as $\varepsilon\rightarrow 0$,

$\int_\Omega (\varphi^2-1)\frac{U_\varepsilon^p}{|x|^{bp}}=O(\varepsilon ^\frac{p}{p-2}),$

(3.1)

$\int_\Omega (\varphi^p-1)\frac{U_\varepsilon^p}{|x|^{bp}}=O(\varepsilon ^\frac{p}{p-2}), $

(3.2)

$\int_{\mathbb{R}^N\setminus\Omega}\frac{U_\varepsilon^p}{|x|^{bp}}=O(\varepsilon^\frac{p}{p-2}), $

(3.3)

$\int_{\Omega }{\frac{|{{u}_{\varepsilon }}{{|}^{2}}}{|x{{|}^{dD}}}}=\left\{ \begin{align} & O({{\varepsilon }^{\frac{2(1+a)-dD}{(p-2){{\beta }^{*}}}}}),{{\beta }^{*}}>1+a-dD/2, \\ & O({{\varepsilon }^{\frac{2}{p-2}}}|\ln \varepsilon |),{{\beta }^{*}}=1+a-dD/2, \\ & O({{\varepsilon }^{\frac{2}{p-2}}}){{\beta }^{*}}<1+a-dD/2, \\ \end{align} \right.$

(3.4)

$\int_{\Omega }{\frac{|{{u}_{\varepsilon }}{{|}^{2}}}{|x{{|}^{2a}}}}=\left\{ \begin{align} & O({{\varepsilon }^{\frac{2}{(p-2){{\beta }^{*}}}}}),{{\beta }^{*}}>1, \\ & O({{\varepsilon }^{\frac{2}{p-2}}}|\ln \varepsilon |),{{\beta }^{*}}=1, \\ & O({{\varepsilon }^{\frac{2}{p-2}}}),{{\beta }^{*}}<1. \\ \end{align} \right.$

(3.5)

Lemma 3.1  Suppose (H1) and (H2) hold. Then $J(u, v)$ satisfies the $(\mbox{PS})_c$ condition for $c<c^*:=(\frac12-\frac{1}{p})S_{\eta, \alpha, \beta}(\mu).$

Proof  The proof is standard (see [5] for example) and we omit it.

Set

$\tilde{c}=\inf\limits_{\psi\in \Psi}\sup\limits_{t\in [0, 1]}J(\gamma(t)),$

where $\Psi:=\{\gamma \in C([0, 1], W), \gamma(0)=(0, 0), J(\gamma(1))<0\}$. It is easy to prove that $J(u, v)$ satisfies the geometry structure of mountain pass lemma. Therefore by mountain pass lemma, $J$ admits a P.S. sequence at level $\tilde{c}$.

Set

$\begin{equation}\label{3.6} c^*=\inf\limits_{(u, v)\in W}\bigg\{\sup\limits_{t\geq 0}J(tu, tv); (u, v)\not\equiv (0, 0)\bigg\}. \end{equation}$

(3.6)

Then $ \tilde{c}=c^*$ (see [21]). Hence from Lemma 3.1, to find solutions for (1.1), we only need to verify

$\begin{equation}\label{3.7} c^*\leq \sup\limits_{t\geq 0}J(tu_\varepsilon, t\tau_{\min}u_\varepsilon)<(\frac{1}{2}-\frac{1}{p})S_{\eta, \alpha, \beta}(\mu)^{\frac{p}{p-2}}. \end{equation}$

(3.7)

Proof of Theorem 1.1 Under assumption (H2), we have

$ \Lambda_1(u^2+v^2)\leq a_1u^2+2a_2uv+a_3v^2 \leq \Lambda_2(u^2+v^2), \forall u, v\in H. $

Meanwhile, for any $v\in D^{1, 2}(\mathbb{R}^N, |x|^{-2a})$ and $\varphi\in \mathcal{D}$, we see

$\begin{eqnarray*} \int_\Omega \varphi v L(\varphi v)=\int_\Omega \varphi^2 v L(v)+\int_\Omega|x|^{-2a} v^2|\nabla\varphi|^2. \end{eqnarray*}$

Taking $v=U_\varepsilon$, we obtain

$\begin{eqnarray*} \int_\Omega \varphi U_\varepsilon L(\varphi U_\varepsilon)&=&\int_\Omega \varphi^2 U_\varepsilon L(U_\varepsilon)+\int_\Omega|x|^{-2a} U_\varepsilon^2|\nabla\varphi|^2\\ &=&\int_\Omega \varphi^2 \frac{U^p_\varepsilon}{|x|^{bp}}+\int_\Omega|x|^{-2a} U_\varepsilon^2|\nabla\varphi|^2. \end{eqnarray*}$

So for $t>0$,

$\begin{eqnarray*} J(tu_\varepsilon,t\tau_{\min}u_\varepsilon)&=&\frac{t^2(1+\tau_{\min}^2)}{2}\int_\Omega \Bigl(|x|^{-2a}|\nabla u_\varepsilon|^2-\mu\frac{ u_\varepsilon^2}{|x|^{2(1+a)}}\Bigl)\\ &&-\frac{t^2}{2}\int_\Omega\frac{(a_1+2a_2\tau_{\min}+a_3\tau_{\min}^2) u_\varepsilon^2}{|x|^{dD}}-\frac{t^p}{p}\int_\Omega \frac{(1+\eta\tau_{\min}^\beta+\tau_{\min}^p)|u_\varepsilon|^p}{|x|^{bp}}\\ &=&\frac{t^2(1+\tau_{\min}^2)}{2}\int_\Omega \varphi U_\varepsilon L(\varphi U_\varepsilon)-\frac{(a_1+2a_2\tau_{\min}+a_3\tau_{\min}^2) t^2}{2}\int_\Omega \frac{U^2_\varepsilon\varphi^2}{|x|^{dD}}\\ &&-\frac{t^p(1+\eta\tau_{\min}^\beta+\tau_{\min}^p)}{p}\int_\Omega \frac{\varphi^{p}U_\varepsilon^p}{|x|^{bp}}\\ &\leq&\frac{t^2(1+\tau_{\min}^2)}{2}\int_\Omega \varphi^2 \frac{U^p_\varepsilon}{|x|^{bp}}+\frac{t^2(1+\tau_{\min}^2)}{2}\int_\Omega\Bigl(|x|^{-2a} |\nabla\varphi|^2 -\Lambda_1\frac{\varphi^2}{|x|^{dD}}\Bigl)U_\varepsilon^2\\ &&-\frac{t^p(1+\eta\tau_{\min}^\beta+\tau_{\min}^p)}{p}\int_\Omega \frac{\varphi^{p}U_\varepsilon^p}{|x|^{bp}}\\ &=&\Bigr(\frac{t^2(1+\tau_{\min}^2)}{2}-\frac{t^p(1+\eta\tau_{\min}^\beta+\tau_{\min}^p)}{p}\Bigr)\int_{{\mathbb R}^N} \frac{U_{\varepsilon}^{p}}{|x|^{p}}\\ &&+\frac{t^2(1+\tau_{\min}^2)}{2}\int_\Omega\Bigl(|x|^{-2a} |\nabla\varphi|^2 -\Lambda_1 \frac{\varphi^2}{|x|^{dD}}\Bigl)U_\varepsilon^2\\ &&+\Bigl\{\frac{t^2(1+\tau_{\min}^2)}{2}\int_\Omega (\varphi^2-1)\frac{U_\varepsilon^p}{|x|^{bp}}-\frac{t^p(1+\eta\tau_{\min}^\beta+\tau_{\min}^p)}{p}\int_\Omega (\varphi^p-1)\frac{U_\varepsilon^p}{|x|^{bp}}\\ &&-\Bigr(\frac{t^2(1+\tau_{\min}^2)}{2}-\frac{t^p(1+\eta\tau_{\min}^\beta+\tau_{\min}^p)}{p}\Bigr)\int_{\mathbb{R}^N\setminus\Omega}\frac{U_\varepsilon^p}{|x|^{bp}}\Bigl\}\\ &=:& \Bigr(\frac{t^2(1+\tau_{\min}^2)}{2}-\frac{t^p(1+\eta\tau_{\min}^\beta+\tau_{\min}^p)}{p}\Bigr)\int_{{\mathbb R}^N} \frac{U_{\varepsilon}^{p}}{|x|^{p}}\\ &&+\frac{t^2(1+\tau_{\min}^2)}{2}\int_\Omega\Bigl(|x|^{-2a} |\nabla\varphi|^2-\Lambda_1\frac{\varphi^2}{|x|^{dD}}\Bigl)U_\varepsilon^2+H(t,\varphi,\varepsilon). \end{eqnarray*}$

From (3.1)-(3.5), we see that for $\varepsilon$ sufficiently small, there exists bounded $t_{\varepsilon}$ such that $\sup\limits_{t\geq 0}J(tu_\varepsilon, t\tau_{\min}u_\varepsilon)=J(t_{\varepsilon}u_\varepsilon, t_{\varepsilon}\tau_{\min}u_\varepsilon)$. Hence $H(t, \varphi, \varepsilon)=O(\varepsilon^{\frac{p}{p-2}})$ and

$\begin{eqnarray*} \underset{t\ge 0}{\mathop{\sup }}\,J(tu_\varepsilon, t\tau_{\min}u_\varepsilon)&\leq& \max\limits_{t\geq 0}\Bigr(\frac{t^2(1+\tau_{\min}^2)}{2}-\frac{t^p(1+\eta\tau_{\min}^\beta+\tau_{\min}^p)}{p}\Bigr)\int_{\mathbb{R}^N} \frac{U_{\varepsilon}^{p}}{|x|^{p}}\\ &&+C\int_\Omega\Bigl(|x|^{-2a} |\nabla\varphi|^2 -\Lambda_1\frac{\varphi^2}{|x|^{dD}}\Bigl)U_\varepsilon^2+O(\varepsilon^{\frac{p}{p-2}})\\ &=&(\frac{1}{2}-\frac{1}{p})S_{\eta, \alpha, \beta}(\mu)^{\frac{p}{p-2}}+C\int_\Omega\Bigl(|x|^{-2a} |\nabla\varphi|^2 -\Lambda_1\frac{\varphi^2}{|x|^{dD}}\Bigl)U_\varepsilon^2+O(\varepsilon^{\frac{p}{p-2}}). \end{eqnarray*}$

Denote $G(0, x)=\frac{1}{|x|^\gamma}$, which is the Green function of operator $L$.

(ⅰ) If $0 \le \mu \le (\sqrt{\bar{\mu}}-a)^2-(1+a-\frac{dD}{2})^2$, then $2\gamma+dD\geq N$.

On the other hand, as $\varepsilon\to 0$,

$\begin{equation}\label{3.8} \frac{\varepsilon^{-\frac{2}{p-2}}U_\varepsilon^2}{|x|^{2a}}\to \frac{CG^2(0, x)}{|x|^{2a}}=\frac{C}{|x|^{2a+2\gamma}}, \frac{\varepsilon^{-\frac{2}{p-2}}U_\varepsilon^2}{|x|^{dD}}\to \frac{CG^2(0, x)}{|x|^{dD}}=\frac{C}{|x|^{2\gamma+dD}}, \end{equation}$

(3.8)

thus

$\begin{equation}\label{3.9} \varepsilon^{-\frac{2}{p-2}}\int_\Omega\frac{U_\varepsilon^2}{|x|^{dD}}\to+\infty. \end{equation}$

(3.9)

Hence, considering $2a<dD$, we see that for a fixed $\varphi\in\mathcal{D}$, and any $\Lambda_1\in(0, \lambda^*_1(\mu))$, we can choose $\varepsilon$ sufficiently small such that

$\begin{equation}\label{3.10} \varepsilon^{\frac{p}{p-2}}=o\Bigl(\int_\Omega\Bigl(|x|^{-2a} |\nabla\varphi|^2 -\Lambda_1\frac{\varphi^2}{|x|^{dD}}\Bigl)U_\varepsilon^2\Bigl), \int_\Omega\Bigl(|x|^{-2a} |\nabla\varphi|^2 -\Lambda_1\frac{\varphi^2}{|x|^{dD}}\Bigl)U_\varepsilon^2<0. \end{equation}$

(3.10)

Therefore, for $\varepsilon$ small enough,

$ \sup\limits_{t\geq 0}J(tu_\varepsilon, t\tau_{\min}u_\varepsilon)<(\frac{1}{2}-\frac{1}{p})S_{\eta, \alpha, \beta}(\mu)^{\frac{p}{p-2}}, $

which is exactly (3.7).

(ⅱ) If $(\sqrt{\bar{\mu}}-a)^2-(1+a-\frac{dD}{2})^2<\mu<(\sqrt{\bar{\mu}}-a)^2$, then $2\gamma+dD< N$. So when $\varepsilon\to 0$,

$\begin{equation}\label{3.11} \varepsilon^{-\frac{2}{p-2}}\int_\Omega\Bigl(|x|^{-2a} |\nabla\varphi|^2 -\Lambda_1\frac{\varphi^2}{|x|^{dD}}\Bigl)U_\varepsilon^2 \to C\int_\Omega\Bigl(\frac{|\nabla\varphi|^2}{|x|^{2a+2\gamma}} -\Lambda_1\frac{\varphi^2}{|x|^{dD+2\gamma}}\Bigl)<\infty. \end{equation}$

(3.11)

By Lemma 2.4 and density arguments, for any $\Lambda_1\in (\lambda_1(\mu), \lambda^*_1(\mu))$, there exists $\varphi\in \mathcal{D}$, such that (3.10) holds for $\varepsilon$ sufficiently small. Hence we also obtain (3.7).

Let $(u_0, v_0)$ be the positive solution of (1.1) obtained in Theorem 1.1 and set $c_0:=J(u_0, v_0)$. From [26], we can infer that $c_0$ can be characterized by ${c_0} = \mathop {\min }\limits_{(u,v) \in \Lambda } J(u,v)$, where

$\begin{eqnarray*} \Lambda:&=&\{(u, v)\in W, (u, v)\geq 0, \langle J'(u, v), (u, v)\rangle=0\}\\ &=&\left\{ \begin{array}{ll} (u, v)|(u, v)\in W\setminus \{(0, 0)\}, (u, v)\geq0, \\ \frac{\displaystyle\int_\Omega(\frac{|u|^p+|v|^p+\eta|u|^\alpha |v|^\beta}{|x|^{bp}})}{\displaystyle\int_\Omega(\frac{|\nabla u|^2+|\nabla v|^2}{|x|^{2a}}-\mu\frac{u^2+v^2}{|x|^{2(1+a)}}-\frac{a_1u^2+2a_2uv+a_3v^2}{|x|^{dD}})}=1 \end{array} \right\}. \end{eqnarray*}$

Let $g(u, v)$ be the functional defined in $W$ by

$\begin{eqnarray*} g(u, v):&=&\left\{ \begin{array}{ll} 0, (u, v)=(0, 0), \\ \frac{\displaystyle\int_\Omega(\frac{|u|^p+|v|^p+\eta|u|^\alpha |v|^\beta}{|x|^{bp}})}{\displaystyle\int_\Omega(\frac{|\nabla u|^2+|\nabla v|^2}{|x|^{2a}}-\mu\frac{u^2+v^2}{|x|^{2(1+a)}}-\frac{a_1u^2+2a_2uv+a_3v^2}{|x|^{dD}})}, (u, v)\neq (0, 0). \end{array} \right. \end{eqnarray*}$

Set $u^+=\max\{u, 0\}, u^-=\max\{-u, 0\}$. Define

$\begin{eqnarray*} &&E:=\{(u, v)\in W|g(u^+, v^+)=g(u^-, v^-)=1\}, \\ &&F:=\{(u, v)\in W||g(u^+, v^+)-1|<\frac12, |g(u^-, v^-)-1|<\frac12\}, \end{eqnarray*}$

then $E\neq\emptyset$. Arguing as in [11], let $\Im$ denote the cone of nonnegative functions in $W$ and $\Sigma$ be the set of maps $\sigma$ such that

(ⅰ) $\sigma\in \mathcal{C}(\mathcal{D}, W), \mathcal{D}=[0,1]\times [0,1].$

(ⅱ) $\sigma(s, 0)\in \Im, \sigma(0, s)\in -\Im, \sigma(1, s)\in -\Im, \forall s\in [0,1].$

(ⅲ) $(J\cdot\sigma)(s, 1)\leq 0, (g\cdot\sigma)(s, 1)\geq2, \forall s\in [0,1].$

We claim that $\Sigma\neq \emptyset$. In fact, for any $(u, v)\in W$ with $(u^+, v^+)\neq (0, 0), (u^-, v^-)\neq (0, 0)$. Set

$ \sigma=\sigma(s_1, s_2)=ks_2(1-s_1)(u^+, v^+)-ks_1s_2(u^-, v^-), (s_1, s_2)\in \mathcal{D}, $

then $\sigma \in \Sigma$ for $k>0$ large enough.

Let $\bar{F}$ be the closure of $F$. Then we have the following result.

Lemma 4.1  There exists a sequence $\{(u_n, v_n)\}\subset \bar{F}$ such that

$ J(u_n, v_n)\rightarrow c_1, J'(u_n, v_n)\rightarrow 0, n\rightarrow \infty. $

Furthermore,

$\mathop {\inf }\limits_{\sigma \in \Sigma } \mathop {\sup }\limits_{(u,v) \in E} J(u,v) = \mathop {\inf }\limits_{(u,v) \in E} J(u,v).$

Proof  The proof is similar to that of [25]. Here we omit it.

Lemma 4.2  Suppose that (H1)--(H3) hold. If $c_1<c_0+c^*$ and $\{(u_n, v_n)\}\subset \bar{F}$ satisfies

$ J(u_n, v_n)\rightarrow c_1, J'(u_n, v_n)\rightarrow 0, n\rightarrow \infty, $

then $\{(u_n, v_n)\}$ is relatively compact in $W$.

Proof  According to Lemma 2.1 and following the same lines as in [25], we can obtain the result. Here we omit it.

Lemma 4.3  Suppose that (H1), (H2), $\eta=0$ and $\beta^*>\max\{2(1+a)-dD, 1\}$, then $c_1<c_0+c^*.$

Proof  By the proof of Theorem 1.1, we infer that $\tau_{\min}=0$ and $S_{0, \alpha, \beta}=S(\mu)$. In this case, $c^*=\frac{1}{N}S(\mu)^{\frac{N}{2}}.$ By Lemma 4.1, it suffices to show that

$\mathop {\sup }\limits_{{s_1},{s_2} \in \mathbb{R}} {\mathop{\sup }}\,J(s_1(u_0, v_0)+s_2(u_\varepsilon, 0))<c_0+\frac{1}{N}S(\mu)^{\frac{N}{2}}. $

Since

$\begin{eqnarray*} &&\lim_{|s_1|+|s_2|\rightarrow 0}J(s_1(u_0, v_0)+s_2(u_\varepsilon, 0))=0, \\ &&\lim_{|s_1|+|s_2|\rightarrow \infty}J(s_1(u_0, v_0)+s_2(u_\varepsilon, 0))=-\infty, \end{eqnarray*}$

we may assume that there exist constants $0<C_1<C_2$ such that $ C_1\leq |s_i|\leq C_2, i=1, 2. $ Note that the following elementary inequality holds: $\forall q\in [1, +\infty)$, there exists a constant $C=C(q)>0$ such that

$ ||a+b|^q-|a|^q-|b|^q|\leq C(|a|^{q-1}|b|+|a||b|^{q-1}), \forall a, b \in \mathbb{R}. $

Since $(u_0, v_0)$ is a positive solution of (1.1), we have that $ \langle J'(u_0, v_0), (\varphi, \psi)=0, $ i.e.,

$\begin{eqnarray*} \langle J'(u_0, v_0), (\varphi, \psi)\rangle&=&\int_\Omega(\frac{\nabla u_0\nabla\varphi+\nabla v_0\nabla\psi}{|x|^{2a}}-\mu\frac{u_0\varphi+v_0\psi}{|x|^{2(1+a)}})\\ &&-\int_{\Omega}\frac{a_1u_0\varphi+a_2v_0\varphi+a_2u_0\psi+a_3v_0\psi}{|x|^{dD}}\\ &&-\int_\Omega\frac{|u_0|^{p-2}u_0\varphi+|v_0|^{p-2}v_0\psi}{|x|^{bp}}. \end{eqnarray*}$

In particular, $ \langle J'(u_0, v_0), (u_\varepsilon, 0)\rangle =0. $ Consequently,

$\begin{eqnarray*} && J(s_1(u_0, v_0)+s_2(u_\varepsilon, 0))\\ &=&\frac{1}{2}\int_\Omega \frac{|\nabla (s_1u_0+s_2u_\varepsilon)|^2+|\nabla (s_1v_0)|^2|}{|x|^{2a}}-\mu\frac{ |s_1u_0+s_2u_\varepsilon|^2+|s_1v_0|^2}{|x|^{2(1+a)}}\Bigl)\\ &&- \frac{1}{2}\int_\Omega\frac{a_1(s_1u_0+s_2u_\varepsilon)^2+2a_2(s_1u_0+s_2u_\varepsilon)(s_1v_0)+a_3(s_1v_0)^2}{|x|^{dD}}\\ &&- \frac{1}{p}\int_\Omega \frac{|s_1u_0+s_2u_\varepsilon|^p)+|s_1v_0|^p}{|x|^{bp}}\\ &\leq&J(s_1u_0, s_1v_0)+J(s_2u_\varepsilon, 0)+C\int_\Omega \frac{u_0u_\varepsilon^{p-1}+u_0^{p-1}u_\varepsilon}{|x|^{bp}}. \end{eqnarray*}$

From Lemma 2.4, it follows that

$\begin{equation}\label{4.1} \begin{array}{ll} &\int_\Omega \frac{u_0u_\varepsilon^{p-1}}{|x|^{bp}} \leq C \int_{B_\rho(0)}\frac{u_\varepsilon^{p-1}}{|x|^{bp+\nu}}\\ =&C \int_{B_\rho(0)}|x|^{-\nu-bp}\frac{\varepsilon^{\frac{p-1}{p-2}}}{|x|^{(p-1)\nu}(\varepsilon +|x|^{(p-2)\beta})^{2(p-1)/p-2}}dx\\ =&C \int_0^{\rho \varepsilon^{-1/(p-2)\beta}}\varepsilon^{\frac{1}{(p-2)\beta}}\frac{r^{N-1-p(b+\nu)}}{(1+r^{(p-2)\beta})^{2(p-1)/p-2}}dr\\ =&C \int_0^{\rho \varepsilon^{-1/(p-2)\beta}}\varepsilon^{p\beta-1}\frac{1}{(1+r^{(p-2)\beta})^{2(p-1)/p-2}}dr\\ =&C \varepsilon^{\frac{1}{p-2}}. \end{array} \end{equation}$

(4.1)

Similarly,

$\begin{equation}\label{4.2} \int_\Omega \frac{u_0^{p-1}u_\varepsilon}{|x|^{bp}} \leq C\varepsilon^{\frac{1}{p-2}}. \end{equation}$

(4.2)

Arguing as the proof of Theorem 1.1 and by (3.1)-(3.4), (4.1)-(4.2), we have

$\begin{eqnarray*} \begin{array}{ll} \mathop {\sup }\limits_{{s_1},{s_2} \in \mathbb{R}} J(s_1(u_0, v_0)+s_2(u_\varepsilon, 0)) &\leq c_0+\frac{1}{N}S(\mu)^{\frac{N}{2}}+O(\varepsilon^{\frac{1}{p-2}})+O(\varepsilon^{\frac{2}{(p-2)\beta}})-C\varepsilon^{\frac{2(1+a)-dD}{(p-2)\beta}}\\ &<c_0+\frac{1}{N}S(\mu)^{\frac{N}{2}}, \end{array} \end{eqnarray*}$

where we use the fact that $\beta^*>\max\{2(1+a)-dD, 1\}$.

Proof of Theorem 1.2  By Lemma 4.1--Lemma 4.3, there exists a sequence $\{(u_n, v_n)\}\subset \bar{F}$ such that

$ J(u_n, v_n)\rightarrow c_1<c_0+c^*, J'(u_n, v_n)\rightarrow 0, n\rightarrow \infty. $

Passing to a subsequence if necessary, $(u_n, v_n)\rightarrow (u, v)$ in $W$ as $n\rightarrow \infty$. Therefore $(u, v)$ is a critical point of $J$ and solves (1.1). Since $(u_n, v_n)\in \bar{F}$, we infer that $(u, v)\in \bar{F}$. Moreover, we have $u\not\equiv 0, v\not\equiv 0$. It follows from the Hölder and Young inequality that there exists a constant $\delta>0$ such that

$ \|(u^+, v^+)\|\geq \delta, \|(u^-, v^-)\|\geq \delta. $

Therefore $(u, v)$ is a sign-changing solution of (1.1) and $(-u, -v)$ is also a solution. So far, the proof of Theorem 1.2 is completed.