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

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

\begin{document}
\title[\hfilneg EJDE-2017/183\hfil Fractional-Laplace system with singular nonlinearity]
{Multiplicity results of fractional-Laplace system with sign-changing
and singular nonlinearity}

\author[S. Goyal \hfil EJDE-2017/183\hfilneg]
{Sarika Goyal}

\address{Sarika Goyal \newline
Department of Mathematics,
Bennett University, Greater Noida,
Uttar Pradesh, India}
\email{sarika1.iitd@gmail.com, sarika.goyal@bennett.edu.in}


\dedicatory{Communicated by Jesus Ildefonso Diaz}

\thanks{Submitted July 20, 2016. Published July 18, 2017.}
\subjclass[2010]{35A15, 35J75, 35R11}
\keywords{Fractional-Laplacian system; singular nonlinearity;
\hfill\break\indent   sign-changing weight function}

\begin{abstract}
 In this article, we study the following fractional-Laplacian system
 with singular nonlinearity
 \begin{gather*}
 (-\Delta)^s u =  \lambda f(x) u^{-q}
 +  \frac{\alpha}{\alpha+\beta}b(x) u^{\alpha-1} w^\beta\quad \text{in }\Omega \\
 (-\Delta)^s w =  \mu g(x) w^{-q}+ \frac{\beta}{\alpha+\beta} b(x) u^{\alpha}
 w^{\beta-1}\; \text{in } \Omega \\
 u, w>0\text{ in }\Omega, \quad u = w = 0  \text{ in } \mathbb{R}^n \setminus\Omega,
 \end{gather*}
 where $\Omega$ is a bounded domain in $\mathbb{R}^n$ with smooth
 boundary $\partial \Omega$, $n>2s$, $s\in(0,1)$, $0<q<1$, $\alpha>1$, $\beta>1$
 satisfy $2<\alpha+\beta< 2_{s}^*-1$ with $2_{s}^*=\frac{2n}{n-2s}$, the pair
 of parameters $(\lambda,\mu)\in \mathbb{R}^2\setminus\{(0,0)\}$. The weight functions
 $f,g: \Omega\subset\mathbb{R}^n \to \mathbb{R}$ such that $0<f$,
 $g\in L^{\frac{\alpha+\beta}{\alpha+\beta-1+q}}(\Omega)$, and
 $b:\Omega\subset\mathbb{R}^n \to \mathbb{R}$ is a sign-changing function such
 that $b(x)\in L^{\infty}(\Omega)$. Using variational methods, we show existence
 and multiplicity of positive solutions  with respect to the
 pair of parameters $(\lambda,\mu)$.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{corollary}[theorem]{Corollary}
\allowdisplaybreaks

\section{Introduction}


Let $\Omega\subset \mathbb{R}^n$ is a bounded domain with smooth boundary, 
$n>2s$ and $s\in (0,1)$. We consider the following fractional system with
 singular nonlinearity
 \begin{equation} \label{ePlm}
\begin{gathered}
(-\Delta)^s u =  \lambda f(x) u^{-q}+  \frac{\alpha}{\alpha+\beta}b(x)
 u^{\alpha-1}w^\beta\quad \text{in } \Omega \\
(-\Delta)^s w =  \mu g(x) w^{-q}+ \frac{\beta}{\alpha+\beta}
  b(x) u^{\alpha} w^{\beta-1}\quad \text{in }\Omega \\
 u, w>0\text{ in }\Omega, \quad u = w = 0  \text{ in } \mathbb{R}^n \setminus\Omega.
\end{gathered}
\end{equation}
Here, $(-\Delta)^{s}$ is the fractional Laplacian operator defined as
\begin{equation*}
(-\Delta)^{s} u(x)= -\frac{1}{2} \int_{\mathbb{R}^n}\frac{u(x+y)+u(x-y)
-2u(x)}{|y|^{n+2s}} dy \quad \text{for all } x\in \mathbb{R}^n.
\end{equation*}
We assume the following assumptions on $f$ and $g$:
 \begin{itemize}
  \item[(A1)] $f,g:   \Omega\subset\mathbb{R}^n \to \mathbb{R}$ such that
$0< f,g\in L^{q^*}(\Omega)$, where $q^*=\frac{\alpha+\beta}{\alpha+\beta-1+q}$.

  \item[(A2)] $b:\Omega\subset\mathbb{R}^n \to \mathbb{R}$ is a sign-changing
function such that $b^{+}=\max\{f,0\}\not\equiv 0$ and
$b(x)\in L^{\infty}(\Omega)$.
\end{itemize}
Also the pair of parameters $(\lambda,\mu)\in \mathbb{R}^2\setminus\{(0,0)\}$,
$0<q <1$ and $\alpha>1$, $\beta>1$ satisfy $2<\alpha+\beta<2^{*}_{s}-1$,
with $2^{*}_{s}= \frac{2n}{n-2s}$, known as fractional critical Sobolev exponent.

In this work, we prove the existence of multiple positive solutions for a 
system of fractional operator with singular and sign changing nonlinearity
 by studying the nature of Nehari manifold with respect to the pair of
 parameter $(\lambda,\mu)$. These same result can be easily extended to 
the following $p-$fractional Laplacian system with singular and sign-changing 
nonlinearity
 \begin{equation}   \label{ePlmp}
\begin{gathered}
(-\Delta)^s_p u =  \lambda f(x) u^{-q}
+  \frac{\alpha}{\alpha+\beta}b(x) u^{\alpha-1}w^\beta\quad \text{in }
\Omega \\
(-\Delta)^s_p w =  \mu g(x) w^{-q}+ \frac{\beta}{\alpha+\beta}
 b(x) u^{\alpha} w^{\beta-1}\quad \text{in }\Omega \\
 u, w>0\text{ in }\Omega, \quad u = w = 0 \text{ in } \mathbb{R}^n \setminus\Omega.
\end{gathered}
\end{equation}
where $0<q<1\leq p-1<\alpha+\beta <p_{s}^{*}-1$ and 
$(-\Delta)_{p}^s$ is a $p-$fractional operator which is defined as
\[
  (-\Delta)_{p}^s u(x) = 2 \lim_{\epsilon\to 0} \int_{\mathbb{R}^n
\setminus B_{\epsilon}(x)} \frac{|u(y)-u(x)|^{p-2}(u(y)-u(x))}{|x-y|^{n+ps}}dy .
\]
This definition is consistent, up to a normalization constant depending on $n$ and
 $s$, with the fractional Laplacian $(-\Delta)^s$, for the case $p=2$.

 For $u=v$, $\alpha=\beta$, $\alpha+\beta=r$, $\lambda=\mu$ and $f=g$, 
 problem \eqref{ePlmp} reduces  to the  $p$-fractional Laplace equation with
 singular nonlinearities
\begin{equation} \label{ePl}
 \begin{gathered}
(-\Delta)_{p}^s =  f(x) w^{-q}+ \lambda b(x) w^r\quad \text{in }
\Omega ,\\
 w>0\text{ in }\Omega, \quad w = 0  \text{ in } \mathbb{R}^n \setminus\Omega.
\end{gathered}
\end{equation}
In \cite{cv}, the author studied the existence and multiplicity of positive 
solutions to problem \eqref{ePl} with sign-changing and singular nonlinearity. 
In the scalar case, the problems involving the fractional operator 
 with singular nonlinearity have been studied by many authors,
 see \cite{tu, bb, fa} and references therein. Also, in \cite{s1}, 
the author used the Caffarelli and Silvestre \cite{cs} approach to obtain the 
multiplicity results for singular and sign-changing nonlinearity.

The fractional power of Laplacian is the infinitesimal generator of 
L\'{e}vy stable diffusion process and arise in anomalous diffusions in plasma, 
population dynamics, geophysical fluid dynamics, flames propagation, 
chemical reactions in liquids and American options in finance. 
For more details, one can see \cite{da, gm} and reference therein. 
Recently the fractional elliptic problem have been investigated by many 
authors for polynomial type nonlinearities see \cite{mp,bb,xy} and 
reference therein. Moreover, by Nehari manifold and fibering maps, 
the author obtained the existence of multiple solutions for $p$-fractional 
equations \cite{ssa,ssi} and reference therein.

 For $s=1$, the paper by Crandall, Robinowitz and Tartar \cite{cr} is the 
starting point on semilinear problem with singular nonlinearity. 
There is a large body of literature on singular nonlinearity see 
\cite{AJ, Am, cr, cp,dh, dhr, gr, gr1, hss, hnc, lm,lmp, J1, j2, jh} 
and reference therein. In \cite{yc}, Chen showed the existence and multiplicity
 of the problem
\begin{gather*}
-\Delta w -\frac{\lambda}{|x|^2}w=  \frac{f(x)}{w^q} + \mu g(x) w^p \quad \text{in }
\Omega\setminus\{0\} \\
  w>0\text{ in } \Omega\setminus\{0\}, \quad w = 0  \text{ in } \partial \Omega,
\end{gather*}
where $0\in \Omega$ is a bounded smooth domain of $\mathbb{R}^n$ with smooth
boundary, $0<\lambda<\frac{(n-2)^2}{4}$, $0<q<1<p<\frac{n+2}{n-2}$, $f(x)>0$ and 
$g$ is sign-changing continuous function.

 The natural space to look for solutions of the problem \eqref{ePlmp} is the 
product space
$W^{s,p}_{0}(\Omega)\times W^{s,p}_{0}(\Omega)$. 
To study \eqref{ePlm}, it is important to encode the `boundary condition'
$u=v=0$ in $\mathbb{R}^n\setminus\Omega$ in the weak formulation.
 Servadei and Valdinoci  \cite{mp}  introduced
the new function spaces  to study the variational functionals related to 
the fractional Laplacian by observing the interaction between $\Omega$ and
$\mathbb{R}^n\setminus\Omega$.

 To the best of our knowledge, there is no work related the fractional
 Laplacian system with singular and sign-changing nonlinearity. 
In this work, we studied the multiplicity results for the system of 
fractional Laplacian equation with singular nonlinearity and sign-changing 
weight function with respect to the pair of parameter $(\lambda,\mu)$.
This work is motivated by the work of Chen and Chen in \cite{yc}. 
But one can not directly extend all the results for fractional Laplacian, 
due to the non-local behavior of the operator and the bounded support of 
the test function is not preserved. Also due to the singularity of the problem, 
the associated functional is not differentiable in the sense of G\^{a}teaux.  
The results obtained here are somehow expected but we show how the results 
arise out of nature of the Nehari manifold.


This article is organized as follows: Section 2 is devoted to some preliminaries 
and notations. We also state our main results. In section 3, we study the 
decomposition of Nehari manifold and the associated energy functional is 
bounded below and coercive. Section 3 contains the existence of a nontrivial 
solutions in $\mathcal{N}_{\lambda,\mu}^{+}$ and $\mathcal{N}_{\lambda,\mu}^{-}$.

 We will use the following notation throughout this paper: 
$\|f\|_{q^*}$, $\|g\|_{q^*}$ denote the norm in 
$L^{\frac{\alpha+\beta}{\alpha+\beta-1+q}}(\Omega)$.

\section{Preliminaries}

In this section we give some definitions and functional settings.
 At the end of this section, we state our main results. 
For this we define $H^{s}(\Omega)$, the usual fractional Sobolev
space 
\[
H^{s}(\Omega):= \big\{w\in L^{2}(\Omega): \frac{(w(x)-w(y))}{|x-y|^{\frac{n}{2}+s}}
\in L^{2}(\Omega\times\Omega)\big\}
\]
 endowed with the norm
\begin{align}\label{2}
\|w\|_{H^{s}(\Omega)}=\|w\|_{L^2(\Omega)}+ \Big(\int_{\Omega\times\Omega}
\frac{|w(x)-w(y)|^{2}}{|x-y|^{n+2s}}\,dx\,dy \Big)^{1/2}.
\end{align}
for details on fractional Sobolev spaces, we refer the reader to \cite{hi}.

Because of the non-localness of the operator, we define the linear space
\begin{align*}
X_0= \Big\{& w:\mathbb{R}^n \to\mathbb{R}: w \text{ is measurable},
w|_{\Omega} \in L^p(\Omega)\\
 &  \frac{w(x)- w(y)}{|x-y|^{\frac{n+2s}{2}}}\in
L^2(Q), \;   w = 0 \text{ a.e. in } \mathbb{R}^n\setminus \Omega\Big\}
\end{align*}
where $Q=\mathbb{R}^{2n}\setminus(\mathcal C\Omega\times \mathcal C\Omega)$ and
 $\mathcal C\Omega := \mathbb{R}^n\setminus\Omega$. The space $X_0$ was firstly
introduced by Servadei and Valdinoci \cite{mp}. 
The space $X_0$ endowed with the norm
\begin{equation} \label{01}
 \|w\|=\Big(\int_{Q}\frac{|w(x)-w(y)|^{2}}{|x-y|^{n+2s}}\,dx\, dy\Big)^{1/2}
\end{equation}
is a Hilbert space.
We note that, the norms in \eqref{2} and \eqref{01} are not same because
$\Omega\times \Omega$ is strictly contained in $Q$.
Let $Y= X_0\times X_0$ be the cartesian product of two reflexive Banach spaces,
 which is also reflexive Banach space with the norm
\begin{align*}
\|(u,w)\|
&=(\|u\|_{X_0}^{2}+\|w\|_{X_0}^{2})^{1/2} \\
&= \left(\int_{Q}\frac{|u(x)-u(y)|^{2}}{|x-y|^{n+2s}}\,dx\, dy
 +\int_{Q}\frac{|w(x)- w(y)|^{2}}{|x-y|^{n+2s}}\,dx\,dy\right)^{1/2}.
\end{align*}
Now we define the space
\[
C_{Y}:= \{(u,w): u,w \in C_{c}^{\infty}(\mathbb{R}^n): u=w=0 \text{ in }
 \mathbb{R}^n\setminus \Omega\}.
\]
Then $C_{ Y}$ is a dense in the space $ Y$.

 Denote 
\begin{gather*}
S:=\inf_{u\in X_{0}}\Big\{\frac{\int_{\mathbb{R}^{2n}}|u(x)-u(y)|^{2}
 |x-y|^{-(n+2s)}\,dx\,dy}
{(\int_{\mathbb{R}}|u|^{\alpha+\beta}dx)^{\frac{2}{\alpha+\beta}}}\Big\}, \\
K_{\lambda,\mu} := \lambda\int_{\Omega}f(x)(u_+)^{1-q}dx + \mu\int_{\Omega}g(x)(w_+)^{1-q} dx.
\end{gather*}

\begin{definition} \label{def2.1} \rm
 A weak solution of  problem \eqref{ePlm} is a function $(u,w)\in  Y$, 
with $u$, $w>0$ in $\Omega$ such that for  every $(\phi,\psi)\in Y$,
\begin{align*}
&\int_{Q}\frac{(u(x)-u(y))(\phi(x)-\phi(y))}{|x-y|^{(n+2s)}}  \,dx\,dy
+\int_{Q}\frac{(w(x)-w(y))(\psi(x)-\psi(y))}{|x-y|^{(n+2s)}}  \,dx\,dy \\
&= \lambda\int_{\Omega} f(x) (u^{-q} \phi)(x)dx  +\mu\int_{\Omega} g(x) (w^{-q} \psi)(x)dx   \\
&\quad + \frac{\alpha}{\alpha+\beta}\int_{\Omega} b(x) (u^{\alpha-1}w^{\beta} \phi)(x) dx
+ \frac{\beta}{\alpha+\beta}\int_{\Omega} b(x) (u^{\alpha}w^{\beta-1} \psi)(x) dx.
\end{align*}
\end{definition}

To show the existence of positive solution of \eqref{ePlm}, we consider the  problem
  \begin{equation*} \label{ePplm}
\begin{gathered}
(-\Delta)^s u =  \lambda f(x) u_{+}^{-q} 
+  \frac{\alpha}{\alpha+\beta}b(x) u_{+}^{\alpha-1}w_{+}^\beta\quad \text{in }
\Omega \\
(-\Delta)^s w =  \mu g(x) w_{+}^{-q}+ \frac{\beta}{\alpha+\beta} 
b(x) u_{+}^{\alpha} w_{+}^{\beta-1}\quad \text{in }\Omega \\
 u, w>0\text{ in }\Omega, \quad u = w = 0  \text{ in } \mathbb{R}^n \setminus\Omega,
\end{gathered}
\end{equation*}
where  $w_+:= \max\{w,0\}$ denotes the positive part of $w$.
Then the function $(u,w)\in  Y$ with  $u, w>0$ in $\Omega\times\Omega$ is a
 weak solution of the problem \eqref{ePplm} if for  every 
$(\phi,\psi)\in Y$, we have
\begin{align*}
&\int_{Q}\frac{(u(x)-u(y))(\phi(x)-\phi(y))}{|x-y|^{(n+2s)}}  \,dx\,dy
+\int_{Q}\frac{(w(x)-w(y))(\psi(x)-\psi(y))}{|x-y|^{(n+2s)}}  \,dx\,dy \\
&= \lambda\int_{\Omega} f(x) (u_{+}^{-q} \phi)(x)dx
 +\mu\int_{\Omega} g(x) (w_{+}^{-q} \psi)(x)dx   \\
&\quad + \frac{\alpha}{\alpha+\beta}\int_{\Omega} b(x) (u_{+}^{\alpha-1}w_{+}^{\beta} \phi)(x) dx
 + \frac{\beta}{\alpha+\beta}\int_{\Omega} b(x) (u_{+}^{\alpha}w_{+}^{\beta-1} \psi)(x) dx.
\end{align*}
We note that if $(u,w)>0$ is a solution of $\eqref{ePplm}$ then one can easily see 
that $(u,w)$ is also a solution $\eqref{ePlm}$. To find the solution of \eqref{ePplm},
 we will use variational approach. So we define the associated functional 
$J_{\lambda,\mu}:  Y \to \mathbb [-\infty,\infty)$ as
\begin{align*}
J_{\lambda,\mu}(u,w)
& =  \frac{1}{2}\|(u,w)\|^2- \frac{1}{1-q}\int_{\Omega} \left(\lambda f(x)u_{+}^{1-q}
 + \mu g(x) w_{+}^{1-q}\right) dx \\
&\quad -\frac{1}{\alpha+\beta} \int_{\Omega} b(x) u_{+}^{\alpha}w_{+}^{\beta}dx.
\end{align*}
Here $J_{\lambda,\mu}$ is not bounded below on $Y$ but is bounded below on
appropriate subset $\mathcal{N}_{\lambda,\mu}$ of $Y$. Therefore in order to obtain the existence
 results, we introduce the Nehari manifold
\[
\mathcal{N}_{\lambda,\mu}= \big\{(u,w)\in Y: \langle J_{\lambda,\mu}'(u,w),(u,w)\rangle=0\big\}
=\big\{(u,w)\in Y : \phi_{u,w}'(1)=0\big\}
\]
where $\langle\cdot,\cdot \rangle$ denotes the duality between $Y$ and its dual
space. Thus $(u,w)\in \mathcal{N}_{\lambda,\mu}$ if and only if
\[
\|(u,w)\|^2-\int_{\Omega}\left(\lambda f(x)u_{+}^{1-q}+ \mu g(x) w_{+}^{1-q}\right)dx
- \int_{\Omega} b(x)u_{+}^{\alpha}w_{+}^{\beta}\,dx =0.
\]
We note that $\mathcal{N}_{\lambda,\mu}$ contains every solution of
\eqref{ePplm}. Now as we know that the Nehari manifold is closely
related to the behavior of the functions $\phi_{u,w}: \mathbb{R}^+\to \mathbb{R}$
defined as $\phi_{u,w}(t)=J_{\lambda,\mu}(tu,tw)$. Such maps are called fiber
maps and were introduced by Drabek and Pohozaev in \cite{DP}. For
$(u,w)\in Y$, we have
\begin{gather*}
\phi_{u,w}(t) 
= \frac{t^2}{2} \|(u,w)\|^2- \frac{t^{1-q}}{1-q} K_{\lambda,\mu}(u,w)
 - \frac{t^{\alpha+\beta}}{\alpha+\beta}\int_{\Omega} b(x)u_{+}^{\alpha}w_{+}^{\beta}\,dx ,\\
\phi_{u,w}'(t)
=  t\|(u,w)\|^2- { t^{-q}} K_{\lambda,\mu}(u,w)  - t^{\alpha+\beta-1} \int_{\Omega} b(x)u_{+}^{\alpha}w_{+}^{\beta}\,dx,\\
\phi_{u,w}''(t)
= \|(u,w)\|^2 +  qt^{-q-1}  K_{\lambda,\mu}(u,w) - (\alpha+\beta-1) t^{\alpha+\beta-2}\int_{\Omega} b(x)u_{+}^{\alpha}w_{+}^{\beta}\,dx.
\end{gather*}
Then it is easy to see that $(tu,tw)\in \mathcal{N}_{\lambda,\mu}$ if and only if
$\phi_{u,w}'(t)=0$ and in particular, $(u,w)\in \mathcal{N}_{\lambda,\mu}$ if
and only if $\phi_{u,w}'(1)=0$. Thus it is natural to split
$\mathcal{N}_{\lambda,\mu}$ into three parts corresponding to local minima,
local maxima and points of inflection. For this we set
\begin{gather*}
\mathcal{N}_{\lambda,\mu}^{\pm}:= \big\{(u,w)\in \mathcal{N}_{\lambda,\mu}: \phi_{u,w}''(1) \gtrless0\big\}
 =\big\{(tu,tw)\in Y : \phi_{u,w}'(t)=0,\; \phi_{u,w}''(t)\gtrless  0\big\},\\
\mathcal{N}_{\lambda,\mu}^{0}:= \big\{(u,w)\in \mathcal{N}_{\lambda,\mu}: \phi_{u,w}''(1) = 0\big\}
=\big\{(tu,tw)\in Y : \phi_{u,w}'(t)=0,\; \phi_{u,w}''(t)= 0\big\}.
\end{gather*}
We also observe that if $(u,w)\in\mathcal{N}_{\lambda,\mu}$ then
$$
\phi_{u,w}''(1)
=\begin{cases}
(1+q) \|(u,w)\|^2 - (\alpha+\beta-1+q)\int_{\Omega} 
b(x)u_{+}^{\alpha}w_{+}^{\beta}\,dx \\
\quad \text{(when we solve for $K$ in $\phi'_{u.w}(1)=0$)},\\[4pt]
(2-\alpha-\beta) \|(u,w)\|^2 +(\alpha+\beta-1+q)  K_{\lambda,\mu}(u,w) \\
\quad \text{(when we solve for $\int_\Omega$ in $\phi'_{u,w}(1)=0$)}.
\end{cases}
$$
 Inspired by \cite{yc}, we show that how variational methods can be used 
to established some existence and multiplicity results for \eqref{ePplm}. 
Our results are as follows.

\begin{theorem}\label{th1}
Suppose $(\lambda,\mu)\in \Gamma$, where
\begin{equation} \label{k1}
\begin{aligned}
\Gamma :=\Big\{&(\lambda,\mu)\in \mathbb{R}^2\setminus\{(0,0)\}: \\
&0<\Lambda:= (|\lambda|\|f\|_{q^*})^{\frac{2}{1+q}}+(|\mu|\|g\|_{q^*})^{\frac{2}{1+q}}
< C(n,\alpha,\beta,q,S)\Big\},
\end{aligned}
\end{equation}
and
\begin{align*}
C(n,\alpha,\beta,q,S)&= \Big(\frac{(1+q)}{(\alpha+\beta-1+q)}
 \Big)^{\frac{2}{\alpha+\beta-2}}
\Big(\frac{\alpha+\beta-2}{\alpha+\beta-1+q}\Big)^{\frac{2}{1+q}} \\
&\quad\times \Big(\frac{1}{\|b\|_{\infty}}\Big)^{\frac{2}{\alpha+\beta-2}}
 S^{\frac{2(\alpha+\beta-1+q)}{(1+q)(\alpha+\beta-2)}}.
\end{align*}
Then  problem \eqref{ePlm} has at least two solutions
\[
(u,w)\in\mathcal{N}_{\lambda,\mu}^{+}, (U,W)\in \mathcal{N}_{\lambda,\mu}^{-}\quad
\text{with } \|(U,W)\|>\|(u,w)\|.
\]
\end{theorem}

\section{Fibering map analysis}

In this section, we show that $\mathcal{N}_{\lambda,\mu}^{\pm}$ is nonempty and
$\mathcal{N}_{\lambda,\mu}^{0}=\{(0,0)\}$. Moreover, $J_{\lambda,\mu}$ is bounded below
and coercive.

\begin{lemma}\label{le1}
Let $(\lambda,\mu)\in\Gamma$. Then for each $(u,w)\in  Y$ with $ K_{\lambda,\mu}(u,w)>0$, we have
the following:
\begin{itemize}
\item[(i)] $ \int_{\Omega}  b(x)u_{+}^{\alpha}w_{+}^{\beta}\,dx\leq 0$, then there exists a unique $0<t_1<t_{\rm max}$ 
such that $(t_1 u,t_1 w)\in \mathcal{N}_{\lambda,\mu}^{+}$ and
$J_{\lambda,\mu}(t_1u,t_1 w)= \inf_{t> 0} J_{\lambda,\mu}(tu,tw)$,

\item[(ii)] $\int_{\Omega}  b(x)u_{+}^{\alpha}w_{+}^{\beta}\,dx > 0$, then there exists  unique $t_1$ and $t_2$ 
with $0<t_1<t_{\rm max}<t_2$ such that $(t_1u,t_1 w)\in \mathcal{N}_{\lambda,\mu}^{+} $,
$(t_2 u,t_2 w)\in \mathcal{N}_{\lambda,\mu}^{-}$ and
$J_{\lambda,\mu}(t_1 u, t_1 w)= \inf_{0\leq t\leq t_{\rm max}} J_{\lambda,\mu}(tu, tw)$,
$J_{\lambda,\mu}(t_2u, t_2 w)= \sup_{t\geq t_1} J_{\lambda,\mu}(tu, tw)$.
\end{itemize}
\end{lemma}

\begin{proof}
For $t>0$, we define
\[
\psi_{u,w}(t)= t^{2-\alpha-\beta}\|(u,w)\|^2 -t^{-\alpha-\beta+1-q}  K_{\lambda,\mu}(u,w)- \int_{\Omega} b(x)u_{+}^{\alpha}w_{+}^{\beta}\,dx.
\]
One can easily see that $\psi_{u,w}(t)\to -\infty$ as $t\to 0^{+}$. Now
\begin{gather*}
\psi_{u,w}'(t)= (2-\alpha-\beta)t^{1-\alpha-\beta}\|(u,w)\|^2 
 +(\alpha+\beta-1+q)t^{-\alpha-\beta-q}  K_{\lambda,\mu}(u,w), \\
\begin{aligned}
\psi_{u,w}''(t)
&= (2-\alpha-\beta)(1-\alpha-\beta)t^{-\alpha-\beta}\|(u,w)\|^2 \\
 &\quad -(\alpha+\beta-1+q)(\alpha+\beta+q)t^{-\alpha-\beta-q-1}  K_{\lambda,\mu}(u,w).
\end{aligned}
\end{gather*}
Then $\psi_{u,w}'(t)=0$ if and only if
\[ 
t=t_{\rm max}:= \Big[\frac{(\alpha+\beta-2)\|(u,w)\|^2}{(\alpha+\beta-1+q) K_{\lambda,\mu}(u,w)}
\Big]^{-\frac{1}{1+q}}. 
\]
Also
\begin{align*}
&\psi_{u,w}''(t_{\rm max}) \\
&= (2-\alpha-\beta)(1-\alpha-\beta)
 \Big[\frac{(\alpha+\beta-2)\|(u,w)\|^2}{(\alpha+\beta-1+q) K_{\lambda,\mu}(u,w)}\Big]
 ^{\frac{\alpha+\beta}{1+q}}\|(u,w)\|^2\\
&\quad - (\alpha+\beta-1+q)(\alpha+\beta+q)
 \Big[\frac{(\alpha+\beta-2)\|(u,w)\|^2}{(\alpha+\beta-1+q) K_{\lambda,\mu}(u,w)}
 \Big]^{\frac{\alpha+\beta+q+1}{1+q}}  K_{\lambda,\mu}(u,w) \\
&= -\|(u,w)\|^2 (\alpha+\beta-2)(1+q)
 \Big[\frac{(\alpha+\beta-2)\|(u,w)\|^2}{(\alpha+\beta-1+q) K_{\lambda,\mu}(u,w)}
 \Big]^{\frac{\alpha+\beta}{1+q}}<0.
\end{align*}
Thus $\psi_{u,w}$ achieves its maximum at $t=t_{\rm max}$. 
Now using the H\"{o}lder's inequality and fractional Sobolev inequality, 
we obtain
\begin{equation} \label{e2}
\begin{aligned}
 K_{\lambda,\mu}(u,w) &\leq |\lambda| \int_{\Omega} |f(x)| |u|^{1-q}dx +|\mu| \int_{\Omega} |g(x)| |w|^{1-q} dx \\
&\leq |\lambda|\|f\|_{q^*} \|u\|_{\alpha+\beta}^{1-q}+ |\mu|\|g\|_{q^*} \|w\|_{\alpha+\beta}^{1-q} \\
&\leq \left((|\lambda|\|f\|_{q^*})^{\frac{2}{1+q}}+(|\mu|\|g\|_{q^*})^{\frac{2}{1+q}}
 \right)^{\frac{1+q}{2}}\Big(\frac{ \|(u,w)\|}{\sqrt{S}}\Big)^{1-q} \\
&= \Lambda^{\frac{1+q}{2}}\Big(\frac{ \|(u,w)\|}{\sqrt{S}}\Big)^{1-q},
\end{aligned}
\end{equation}
where $\Lambda:= (|\lambda|\|f\|_{q^*})^{\frac{2}{1+q}}+(|\mu|\|g\|_{q^*})^{\frac{2}{1+q}}$.
\begin{equation} \label{e3}
\begin{aligned}
\int_{\Omega}  b(x)u_{+}^{\alpha}w_{+}^{\beta}\,dx
&\leq   \|b\|_{\infty}\Big(\frac{\alpha}{\alpha+\beta}\int_{\Omega} |u|^{\alpha+\beta}dx
 +\frac{\beta}{\alpha+\beta}\int_{\Omega} |v|^{\alpha+\beta} dx\Big) \\
&\leq \|b\|_{\infty} \Big(\frac{\|(u,w)\|}{\sqrt{S}}\Big)^{\alpha+\beta}.
\end{aligned}
\end{equation}
Using \eqref{e2} and \eqref{e3} we obtain,
\begin{equation} \label{e40}
\begin{aligned}
\psi_{u,w}(t_{\rm max}) 
&= \frac{(1+q)}{(\alpha+\beta-1+q)}
\Big(\frac{\alpha+\beta-2}{\alpha+\beta-1+q}\Big)^{\frac{\alpha+\beta-2}{1+q}}
 \frac{\|(u,w)\|^{\frac{2(\alpha+\beta-1+q)}{(1+q)}}}
 {[ K_{\lambda,\mu}(u,w)]^{\frac{\alpha+\beta-2}{1+q}}} \\
&\quad - \int_{\Omega} b(x)u_{+}^{\alpha}w_{+}^{\beta}\,dx  \\
&\geq \Big[\frac{(1+q)}{(\alpha+\beta-1+q)}
 \Big(\frac{\alpha+\beta-2}{\alpha+\beta-1+q}\Big)^{\frac{\alpha+\beta-2}{1+q}}
 \Big(\frac{(\sqrt{S})^{(1-q)}}{\Lambda^{\frac{1+q}{2}}}\Big)
 ^{\frac{(\alpha+\beta-2)}{(1+q)}} \\
&\quad -  \|b\|_{\infty} \big(\frac{1}{\sqrt{S}}\big)^{\alpha+\beta} \Big]
 \|(u,w)\|^{\alpha+\beta} \\
&\equiv E_{\lambda,\mu} \|(u,w)\|^{\alpha+\beta}.
\end{aligned}
\end{equation}
where
\begin{align*}
E_{\lambda,\mu} &= \Big[\frac{(1+q)}{(\alpha+\beta-1+q)}
 \Big(\frac{\alpha+\beta-2}{\alpha+\beta-1+q}\Big)^{\frac{\alpha+\beta-2}{1+q}}
 \Big(\frac{(\sqrt{S})^{(1-q)}}
{\Lambda^{\frac{1+q}{2}}}\Big)^{\frac{\alpha+\beta-2}{1+q}}  \\
&\quad - \|b\|_{\infty} \big(\frac{1}{\sqrt{S}}\big)^{\alpha+\beta} \Big]
\end{align*}
Then we see that $E_{\lambda,\mu}=0$ if and only if $\Lambda= C(n,\alpha,\beta,q,S)$, where
\begin{align*}
 C(n,\alpha,\beta,q,S)
&=  \Big(\frac{(1+q)}{(\alpha+\beta-1+q)}  \Big)^{\frac{2}{\alpha+\beta-2}}
\Big(\frac{\alpha+\beta-2}{\alpha+\beta-1+q}\Big)^{\frac{2}{1+q}} \\
&\quad\times \Big(\frac{1}{\|b\|_{\infty}}\Big)^{\frac{2}{\alpha+\beta-2}}
 S^{\frac{2(\alpha+\beta-1+q)}{(1+q)(\alpha+\beta-2)}}.
\end{align*}
 Thus for $(\lambda,\mu)\in\Gamma$, we have $E_{\lambda,\mu}>0$, and therefore it follows
from \eqref{e40} that $\psi_{u,w}(t_{\rm max})>0$.

 (i) If $\int_{\Omega} b(x)u_{+}^{\alpha}w_{+}^{\beta}\,dx\geq 0$, then $\psi_{u,w}(t)\to - \int_{\Omega}  b(x)u_{+}^{\alpha}w_{+}^{\beta}\,dx<0$ as $t\to\infty$.
Consequently, $\psi_{u,w}(t)$ has exactly two points $0<t_{1}<t_{\rm max}<t_{2}$ 
such that
\[ 
\psi_{u,w}(t_{1})=0=\psi_{u,w}(t_{2})\quad\text{and}\quad
 \psi'_{u,w}(t_{1})>0> \psi'_{u,w}(t_{2}).
\]
Now we show that if $\psi_{u,w}(t)=0$ and $\psi'_{u,w}(t)>0$, then 
$(tu,tw)\in\mathcal{N}_{\lambda,\mu}^{+}$.
 \begin{align*}
 \psi_{u,w}(t)=0  
&\Leftrightarrow \|(tu,tw)\|^2 =  K_{\lambda,\mu}(tu,tw)
 +  \int_{\Omega} b(x)(t u)_{+}^{\alpha}(t w)_{+}^{\beta} dx \\
&\Leftrightarrow (tu,tw)\in \mathcal{N}_{\lambda,\mu},
 \end{align*}
 and therefore
 \begin{align*}
&\psi'_{u,w}(t)>0  \\
&\Rightarrow (2-\alpha-\beta)t^{1-\alpha-\beta}\|(u,w)\|^2 
 -(-\alpha-\beta+1-q) t^{-\alpha-\beta-q}  K_{\lambda,\mu}(u,w) >0\\
&\Rightarrow (2-\alpha-\beta)\|(tu,tw)\|^2 +(\alpha+\beta-1+q) 
 \Big[\|(tu,tw)\|^2  \\
&\quad -  \int_{\Omega} b(x)(t u)_{+}^{\alpha}(t w)_{+}^{\beta} dx\Big]>0,\\
&\Rightarrow (1+q)\|(tu,tw)\|^2 - (\alpha+\beta-1+q)  \int_{\Omega} b(x)(t u)_{+}^{\alpha}(t w)_{+}^{\beta} dx >0 \\
&\Rightarrow (tu,tw)\in \mathcal{N}_{\lambda,\mu}^{+}.
\end{align*}
Similarly one can show that if $\psi_{u,w}(t)=0$ and $\psi'_{u,w}(t)<0$, 
then $(tu,tw)\in\mathcal{N}_{\lambda,\mu}^{-}$.

Now $\phi_{u,w}'(t)= t^{\alpha+\beta-1}\psi_{u,w}(t)$. 
Thus $\phi_{u,w}'(t)<0$ in $(0,t_1)$, $\phi_{u,w}'(t)>0$ in $(t_1,t_2)$ and 
$\phi_{u,w}'(t)<0$ in $(t_2,\infty)$. Hence 
$J_{\lambda,\mu}(t_1 u,t_1 w)= \inf_{0\leq t\leq t_{\rm max}} J_{\lambda,\mu}(tu, tw)$,
$J_{\lambda,\mu}(t_2 u,t_2 w)= \sup_{t\geq t_1} J_{\lambda,\mu}(tu,t w)$.
Moreover $(t_1u,t_{1} w)\in \mathcal{N}_{\lambda,\mu}^{+}$ and
$(t_2u,t_{2} w)\in \mathcal{N}_{\lambda,\mu}^{-}$.

(ii) If $\int_{\Omega} b(x)u_{+}^{\alpha}w_{+}^{\beta}\,dx<0$ and $\psi_{u,w}(t)\to - \int_{\Omega}  b(x)u_{+}^{\alpha}w_{+}^{\beta}\,dx>0$ as $t\to\infty$.
 Consequently, $\psi_{u,w}(t)$ has exactly one point $0<t_{1}<t_{\rm max}$ such that
\[ 
\psi_{u,w}(t_{1})=0\; \text{and}\;\psi'_{u,w}(t_{1})>0.
\]
Using $\phi_{u,w}'(t)= t^{\alpha+\beta-1}\psi_{u,w}(t)$, we have 
$\phi_{u,w}'(t)<0$ in $(0,t_1)$, $\phi_{u,w}'(t)>0$ in $(t_1,\infty)$. 
So, $J_{\lambda,\mu}(t_1u,t_1 w)= \inf_{t\geq 0} J_{\lambda,\mu}(tu, tw)$.
Hence, it follows that $(t_1u,t_{1} w)\in \mathcal{N}_{\lambda,\mu}^{+}$.
\end{proof}

\begin{corollary} \label{coro3.2}
Suppose that $(\lambda,\mu)\in\Gamma$, then $\mathcal{N}_{\lambda,\mu}^{\pm}\ne\emptyset$.
\end{corollary}

\begin{proof}
From assumptions (A1) and (A2), we can choose $(u,w) \in Y\setminus\{(0,0)\}$ 
such that $ K_{\lambda,\mu}(u,w)>0$ and $\int_{\Omega} b(x)u_{+}^{\alpha}w_{+}^{\beta}\,dx>0$.
By (ii) of Lemma \ref{le1}, there exists unique $t_1$ and $t_2$ such 
that $(t_1u,t_1 w)\in \mathcal{N}_{\lambda,\mu}^{+}$, $(t_2u, t_2 w)\in \mathcal{N}_{\lambda,\mu}^{-}$.
In conclusion, $\mathcal{N}_{\lambda,\mu}^{\pm}\ne \emptyset$.
\end{proof}

 \begin{lemma}
 For $(\lambda,\mu)\in\Gamma$, we have $\mathcal{N}_{\lambda,\mu}^{0}=\{(0,0)\}$.
 \end{lemma}

\begin{proof}
 We prove this by contradiction. Assume that there exists 
$(0,0)\not\equiv (u,w)\in \mathcal{N}_{\lambda,\mu}^{0}$. Then it follows from
$(u,w)\in \mathcal{N}_{\lambda,\mu}^{0}$ that
\[
(1+q)\|(u,w)\|^2 =  (\alpha+\beta-1+q) \int_{\Omega} b(x)u_{+}^{\alpha}w_{+}^{\beta}\,dx 
\]
and consequently
\begin{align*}
0&= \|(u,w)\|^2 - K_{\lambda,\mu}(u,w) -\int_{\Omega} b(x)u_{+}^{\alpha}w_{+}^{\beta}\,dx\\
&=\frac{(\alpha+\beta-2)}{(\alpha+\beta-1+q)}\|(u,w)\|^2 - K_{\lambda,\mu}(u,w).
\end{align*}
 Therefore, as $(\lambda,\mu)\in \Gamma$ and $(u,w)\not\equiv (0,0)$, we use similar
arguments as those in \eqref{e40} to obtain
\begin{align*}
0&< E_{\lambda,\mu} \|(u,w)\|^{\alpha+\beta}\\
&\leq \frac{(1+q)}{(\alpha+\beta-1+q)}
 \left(\frac{\alpha+\beta-2}{\alpha+\beta-1+q}\right)^{\frac{\alpha+\beta-2}{1+q}} 
 \frac{\|(u,w)\|^{\frac{2(\alpha+\beta-1+q)}{1+q}}}
 {\left[ K_{\lambda,\mu}(u,w)\right]^{\frac{\alpha+\beta-2}{1+q}}}\\
&\quad - \int_{\Omega} b(x)u_{+}^{\alpha}w_{+}^{\beta}\,dx \\
&= \frac{(1+q)}{(\alpha+\beta-1+q)} 
 \left(\frac{\alpha+\beta-2}{\alpha+\beta-1+q}\right)^{\frac{\alpha+\beta-2}{1+q}} 
 \frac{\|(u,w)\|^{\frac{2(\alpha+\beta-1+q)}{1+q}}}
 {\left(\frac{\alpha+\beta-2}{\alpha+\beta-1+q}\|(u,w)\|^2
 \right)^{\frac{\alpha+\beta-2}{1+q}}} \\
&\quad -\frac{(1+q)}{(\alpha+\beta-1+q)}\|(u,w)\|^2
=0,
\end{align*}
a contradiction. Hence $(u,w)=(0,0)$. That is, $\mathcal{N}_{\lambda,\mu}^{0}=\{(0,0)\}$.
\end{proof}

 We note that $\Gamma$ is also related to a gap structure in $\mathcal{N}_{\lambda,\mu}$.

\begin{lemma}\label{le2}
Suppose that $(\lambda,\mu) \in \Gamma$, then there exist a gap structure in
$\mathcal{N}_{\lambda,\mu}$,
\[
\|(U,W)\|> A_{0}> A_{\lambda,\mu}> \|(u,w)\|\quad \text{for all }
 (u,w)\in \mathcal{N}_{\lambda,\mu}^{+}, (U,W)\in \mathcal{N}_{\lambda,\mu}^{-},
\]
where
\begin{gather*}
A_{0}=  \Big[\frac{(1+q)}{(\alpha+\beta-1+q)\|b\|_{\infty}} 
(\sqrt{S})^{\alpha+\beta} \Big]^{\frac{1}{\alpha+\beta-2}}, \\
A_{\lambda,\mu}=\Big[\frac{(\alpha+\beta-1+q)}{(\alpha+\beta-2)}
 \big(\frac{1}{\sqrt{S}}\big)^{1-q} \Big]^{\frac{1}{1+q}}\Lambda^{1/2}.
\end{gather*}
\end{lemma}

\begin{proof}
If $w\in \mathcal{N}_{\lambda,\mu}^{+}\subset \mathcal{N}_{\lambda,\mu}$, then
\begin{align*}
0&< (1+q)\|(u,w)\|^2 - (\alpha+\beta-1+q)\int_{\Omega}  b(x)u_{+}^{\alpha}w_{+}^{\beta}\,dx \\
&= (2-\alpha-\beta) \|(u,w)\|^2 + (\alpha+\beta-1+q) K_{\lambda,\mu}(u,w).
\end{align*}
Hence it follows from \eqref{e2}
\begin{align*}
&(\alpha+\beta-2) \|(u,w)\|^2  \\
&< (\alpha+\beta-1+q) K_{\lambda,\mu}(u,w)\\
& \leq (\alpha+\beta-1+q) \left((|\lambda|\|f\|_{q^*})^{\frac{2}{1+q}}
+(|\mu|\|g\|_{q^*})^{\frac{2}{1+q}}\right)^{\frac{1+q}{2}}
\Big(\frac{ \|(u,w)\|}{\sqrt{S}}\Big)^{1-q}
\end{align*}
which yields
\begin{align*}
\|(u,w)\| &< \Big[\frac{(\alpha+\beta-1+q) }{(\alpha+\beta-2)}
  \big(\frac{1}{\sqrt{S}}\big)^{1-q} \Big]^{\frac{1}{1+q}}
 \left(\left(|\lambda|\|f\|_{q^*})^{\frac{2}{1+q}}
 +(|\mu|\|g\|_{q^*}\right)^{\frac{2}{1+q}}\right)^{1/2}\\
&\equiv A_{\lambda,\mu}.
\end{align*}
 If $(U,W)\in \mathcal{N}_{\lambda,\mu}^{-}$, then  from \eqref{e3} it follows that
\begin{align*}
(1+q)\|(U,W)\|^2 
&< (\alpha+\beta-1+q)\int_{\Omega} b(x)U_{+}^{\alpha}W_{+}^{\beta} dx \\
&\leq  (\alpha+\beta-1+q)\|b\|_{\infty}
 \Big(\frac{\|(U,W)\|}{\sqrt{S}}\Big)^{\alpha+\beta}
\end{align*}
which yields
\[
\|(U,W)\|
> \Big[\frac{(1+q)}{(\alpha+\beta-1+q)\|b\|_{\infty}}
  (\sqrt{S})^{\alpha+\beta} \Big]^{\frac{1}{\alpha+\beta-2}}
\equiv A_0.
\]
Now we show that $A_{\lambda,\mu}=A_0$ if and only if $\Lambda= C(n,\alpha,\beta,q,S)$.
\begin{align*}
\Lambda 
&= C(n,\alpha,\beta,q,S) \\
&= \left(\frac{(1+q)}{\|b\|_{\infty}(\alpha+\beta-1+q)}\right)^{\frac{2}{\alpha+\beta-2}} \left(\frac{\alpha+\beta-2}{\alpha+\beta-1+q}\right)^{\frac{2}{1+q}}  S^{\frac{2(\alpha+\beta-1+q)}{(1+q)(\alpha+\beta-2)}}.\\
&\Leftrightarrow A_{\lambda,\mu} =\Lambda^{1/2}
\Big[\frac{(\alpha+\beta-1+q) }{(\alpha+\beta-2)} \Big(\frac{1}{\sqrt{S}}
 \Big)^{1-q} \Big]^{\frac{1}{1+q}}\\
&= \Big(\frac{(1+q)}{\|b\|_{\infty}(\alpha+\beta-1+q)}
 \Big)^{\frac{1}{\alpha+\beta-2}} 
 \Big(\frac{\alpha+\beta-2}{\alpha+\beta-1+q}\Big)^{\frac{1}{1+q}} 
 S^{\frac{\alpha+\beta-1+q}{(1+q)(\alpha+\beta-2)}} \\
&\quad\times \Big[\frac{(\alpha+\beta-1+q) }{(\alpha+\beta-2)} 
\big(\frac{1}{\sqrt{S}}\big)^{1-q} \Big]^{\frac{1}{1+q}}
 \equiv A_0.
\end{align*}
Thus for all $(\lambda,\mu)\in \Gamma$, we can conclude that
\[
\|(U,W)\|> A_{0}> A_{\lambda,\mu}> \|(u,w)\|\quad
 \text{for all }  (u,w)\in \mathcal{N}_{\lambda,\mu}^{+},\;
(U,W)\in \mathcal{N}_{\lambda,\mu}^{-}.
\]
This completes the proof.
\end{proof}

\begin{lemma}\label{le3}
Suppose that $(\lambda,\mu)\in\Gamma$, then $\mathcal{N}_{\lambda,\mu}^{-}$ is a closed set in
$ Y$-topology.
\end{lemma}

\begin{proof}
Let $\{(U_k,W_k)\}$ be a sequence in $\mathcal{N}_{\lambda,\mu}^{-}$ with
$(U_k,W_k) \to (U,W)$ in $ Y$. Then we have
\begin{align*}
&\|(U_k,W_k)\|^2 \\
&=\lim_{k\to\infty}\|(U_k,W_k)\|^2\\
&= \lim_{k\to\infty} \Big[\int_{\Omega} (\lambda f(x)(U_k)_{+}^{1-q}
 +\mu g(x)(W_k)_{+}^{1-q}) dx + \int_{\Omega} b(x)(U_k)_{+}^{\alpha}(W_k)_{+}^{\beta}
  dx\Big]\\
&= \int_{\Omega} (\lambda f(x)U_{+}^{1-q}+\mu g(x)W_{+}^{1-q}) dx 
 + \int_{\Omega} b(x)U_{+}^{\alpha}W_{+}^{\beta} dx
\end{align*}
and
\begin{align*}
&(1+q) \|(U,W)\| -  (\alpha+\beta-1+q) \int_{\Omega} b(x)U_{+}^{\alpha}W_{+}^{\beta} dx \\
& =\lim_{k\to\infty} \Big[(1+q) \|(U_k,W_k)\|^2 
- (\alpha+\beta-1+q) \int_{\Omega} b(x)(U_k)_{+}^{\alpha}(W_k)_{+}^{\beta} dx\Big]\leq 0,
\end{align*}
i.e. $(U,W)\in  \mathcal{N}_{\lambda,\mu}^{-}\cap \mathcal{N}_{\lambda,\mu}^{0}$.
Since $\{(U_k,W_k)\}\subset \mathcal{N}_{\lambda,\mu}^{-}$, from Lemma \ref{le2} we have
\[
\|(U,W)\|= \lim_{k\to\infty} \|(U_k,W_k)\|\geq A_{\lambda,\mu}>0;
\]
that is, $(U,W)\not\equiv (0,0)$. It follows from Lemma \ref{le1}, that 
$(U,W)\not\in \mathcal{N}_{\lambda,\mu}^{0}$ for any $(\lambda,\mu)\in\Gamma$.
Thus $(U,W)\in \mathcal{N}_{\lambda,\mu}^{-}$.
That is, $\mathcal{N}_{\lambda,\mu}^{-}$ is a closed set in $Y$- topology for any
$(\lambda,\mu)\in\Gamma$.
\end{proof}

\begin{lemma}\label{le4}
 Let $(u,w)\in \mathcal{N}_{\lambda,\mu}^{\pm}$, then for any $\Phi=(\phi,\psi)\in C_{Y}$,
there exists a number $\epsilon>0$ and a continuous function 
$f:B_{\epsilon}(0):=\{v=(v_1,v_2)\in Y : \|v\|<\epsilon\}\to \mathbb{R}^{+}$ such that
\[
f(v_1,v_2)>0, f(0,0)=1\quad\text{and}\quad
 f(v_1,v_2)(u+v_1\phi, w+v_2\psi)\in \mathcal{N}_{\lambda,\mu}^{\pm}\]
for all $ v\in B_{\epsilon}(0)$.
\end{lemma}

\begin{proof}
We give the proof only for the case $(u,w)\in \mathcal{N}_{\lambda,\mu}^{+}$,
the case $\mathcal{N}_{\lambda,\mu}^{-}$ may be proved the same way.
For any $\Phi=(\phi,\psi)\in C_{ Y}$, we define 
$F: Y\times \mathbb{R}^{+}\to \mathbb{R}$ as follows:
\begin{align*}
F(v,t):=&F((v_1,v_2),t)\\
=& t^{1+q}\|(u+v_1 \phi, w+v_2\psi)\|^2 
- t^{\alpha+\beta-1+q}\int_{\Omega} b(x) (u+v_1\phi)_{+}^{\alpha}(w+v_2\psi)_{+}^{\beta} dx\\
&-  K_{\lambda,\mu}(u+v_1\phi,w+v_2\psi)
\end{align*}
Since $w\in \mathcal{N}_{\lambda,\mu}^{+}(\subset \mathcal{N}_{\lambda,\mu})$, we have
\[
F((0,0),1)= \|(u,w)\|^2-  K_{\lambda,\mu}(u,w) - \int_{\Omega}  b(x)u_{+}^{\alpha}w_{+}^{\beta}\,dx =0,
\]
and
\[
\frac{\partial F}{\partial t}((0,0),1)= (1+q)\|(u,w)\|^2 - (\alpha+\beta-1+q) \int_{\Omega}  b(x)u_{+}^{\alpha}w_{+}^{\beta}\,dx>0  .
\]
Applying the implicit function theorem at the point $((0,0),1),$ we have 
that there exists $\bar{\epsilon}>0$ such that for $\|v\|<\bar{\epsilon}$, $v\in  Y$, 
the equation $F((v_1,v_2),t)=0$ has a unique continuous solution 
$t=f(v_1,v_2)>0$. It follows from $F((0,0),1)=0$ that $f(0,0)=1$ and from 
$F((v_1,v_2),f(v_1,v_2))=0$ for $\|v\|<\bar\epsilon$, $v\in Y$ that
\begin{align*}
 0& = f^{1+q}(v) \|(u+v_1\phi,w+v_2\psi)\|^2 - K_{\lambda,\mu}(u+v_1\phi,w+v_2\psi)\\
  &\quad - f^{\alpha+\beta-1+q}(v) \int_{\Omega} b(x) (u+v_1\phi)_{+}^{\alpha}(x)
  (w+v_2\psi)_{+}^{\beta}(x) dx\\
 & = \frac{\|f(v)(u+v_1\phi,w+v_2\psi)\|^2-K_{\lambda,\mu}(f(v)(u+v_1\phi),
 f(v)(w+v_2\psi)) }{f^{1-q}(v)}\\
 &\quad -\frac{ \int_{\Omega} b(x)(f(v)(u+v_1\phi))_{+}^{\alpha}(x) 
 (f(v)(w+v_2\psi))_{+}^{\beta}(x)dx}{f^{1-q}(v)};
\end{align*}
that is,
\[
f(v_1,v_2)(u+v_1\phi,w+v_2\psi)\in \mathcal{N}_{\lambda,\mu}\quad \text{for all }
 v\in  Y, \|v\|<\tilde{\epsilon}.
\]
Since $\frac{\partial F}{\partial t}((0,0),1)>0$ and
\begin{align*}
 &\frac{\partial F}{\partial t}((v_1,v_2),f(v_1,v_2))\\
 & = (1+q)f^{q}(v) \|(u+v_1\phi,w+v_2\psi\|^2  \\
&\quad - (\alpha+\beta-1+q) 
 f^{\alpha+\beta-1+q-1}(v) \int_{\Omega} b(x)(u+v_1\phi)_{+}^{\alpha} (w+v_2\psi)_{+}^{\beta}\\
 & = \frac{(1+q)\|(f(v)(u+v_1\phi),f(v)(w+v_2\psi))\|^2}{f^{2-q}(v)}\\
 &\quad-\frac{(\alpha+\beta-1+q) \int_{\Omega} b(x)(f(v)(u+v_1\phi))_{+}^{\alpha} 
(f(v)(w+v_2\psi))_{+}^{\beta} dx}{f^{2-q}(v)}
\end{align*}
we can take $\epsilon>0$ possibly smaller $(\epsilon<\bar\epsilon)$ such that for any 
$v=(v_1,v_2)\in Y$, $\|v\|<\epsilon$,
\begin{align*}
& (1+q)\|(f(v)(u+v_1\phi),f(v)(w+v_2\psi))\|^2 \\
&- (\alpha+\beta-1+q)\int_{\Omega} b(x)(f(v)(u+v_1\phi))_{+}^{\alpha} 
(f(v)(w+v_2\psi))_{+}^{\beta} dx >0; 
\end{align*}
that is,
\[
f(v_1,v_2)(u+v_1\phi, w+v_2\psi) \in \mathcal{N}_{\lambda,\mu}^{+}\quad\text{for all }
 v=(v_1,v_2)\in B_{\epsilon}(0).
\]
This completes the proof. 
\end{proof}

\begin{lemma}\label{lem3.7}
$J_\lambda$ is bounded below and coercive on $\mathcal{N}_{\lambda,\mu}$.
\end{lemma}

\begin{proof}
For $(u,w)\in \mathcal{N}_{\lambda,\mu}$,  from \eqref{e2} we obtain
\begin{equation} \label{s1}
\begin{aligned}
&J_{\lambda,\mu}(u,w) \\
&= \Big(\frac{1}{2} -\frac{1}{\alpha+\beta}\Big)\|(u,w)\|^2
  - \Big(\frac{1}{1-q} -\frac{1}{\alpha+\beta}\Big) K_{\lambda,\mu}(u,w) \\
&\geq \Big(\frac{1}{2} -\frac{1}{\alpha+\beta}\Big)\|(u,w)\|^2
 - \Big(\frac{1}{1-q} -\frac{1}{\alpha+\beta}\Big)
\Lambda^{\frac{1+q}{2}} \Big(\frac{\|(u,w)\|}{\mathrm{\sqrt{S}}}\Big)^{1-q},
\end{aligned}
\end{equation}
where $\Lambda$ is given in \eqref{k1}. Now consider the function
$\rho: \mathbb{R}^{+}\to \mathbb{R}$ as $\rho(t)= c t^2 - d t^{1-q}$,
where $c,$ $d$ are both positive constants.
One can easily show that $\rho$ is convex($\rho''(t)>0$ for all $t>0$)
 with $\rho(t)\to 0$ as $t\to 0$ and $\rho(t)\to \infty$ as $t\to\infty$.
$\rho$ achieves its minimum at $t_{min}=[\frac{d(1-q)}{2c}]^{\frac{1}{1+q}}$ and
\[
\rho(t_{min})
= c\Big[\frac{d(1-q)}{2c}\Big]^{\frac{2}{1+q}}
- d\Big[\frac{d(1-q)}{2c}\Big]^{\frac{1-q}{1+q}}
= -\frac{(1+q)}{2} d^{\frac{2}{1+q}}\Big(\frac{1-q}{2c}\Big)^{\frac{1-q}{1+q}}.
\]
Applying $\rho(t)$ with $c=\big(\frac{1}{2} -\frac{1}{\alpha+\beta}\big)$,
$d=\big(\frac{1}{1-q} -\frac{1}{\alpha+\beta}\big)\Lambda^{\frac{1+q}{2}}
\big(\frac{1}{\sqrt{S}}\big)^{1-q}$ and
$t=\|(u,w)\|$, $(u,w)\in \mathcal{N}_{\lambda,\mu}$, we obtain from \eqref{s1} that
\[
\lim_{\|(u,w)\|\to\infty} J_{\lambda,\mu}(u,w)
\geq \lim_{t\to\infty} \rho(t) = \infty,
\]
since $0<q<1$. That is $J_{\lambda,\mu}$ is coercive on $\mathcal{N}_{\lambda,\mu}$.
Moreover it follows from \eqref{s1} that
\begin{equation}\label{a1}
J_{\lambda,\mu}(u,w)\geq \rho(t)\geq \rho(t_{min})\quad (\text{a constant}),
\end{equation}
i.e.
\begin{align*}
J_{\lambda,\mu}(u,w)
&\geq -\frac{(1+q)}{2} d^{\frac{2}{1+q}}\left(\frac{1-q}{2c}\right)^{\frac{1-q}{1+q}} \\
&= -\frac{(1+q)(\alpha+\beta-2)}{(1-q)(\alpha+\beta)}
\left(\frac{\alpha+\beta-1+q}{2(\alpha+\beta-2)}\right)^{\frac{2}{1+q}}
\Lambda \big(\frac{1}{\sqrt{S}}\big)^{\frac{2(1-q)}{1+q}}.
\end{align*}
Thus $J_{\lambda,\mu}$ is bounded below on $\mathcal{N}_{\lambda,\mu}$.
\end{proof}

\section{Existence of Solutions in $\mathcal{N}_{\lambda,\mu}^{\pm}$}

Now from Lemma \ref{le3}, $\mathcal{N}_{\lambda,\mu}^{+}\cup \mathcal{N}_{\lambda,\mu}^{0}$ and
$\mathcal{N}_{\lambda,\mu}^{-}$ are two closed sets in $ Y$ provided $(\lambda,\mu)\in\Gamma$.
Consequently, the Ekeland variational principle can be applied to the problem 
of finding the infimum of $J_{\lambda,\mu}$ on both
$\mathcal{N}_{\lambda,\mu}^{+}\cup \mathcal{N}_{\lambda,\mu}^{0}$ and
$\mathcal{N}_{\lambda,\mu}^{-}$.
First, consider $\{(u_k,w_k)\}\subset\mathcal{N}_{\lambda,\mu}^{+}
\cup \mathcal{N}_{\lambda,\mu}^{0}$  with the following properties:
\begin{gather}
J_{\lambda,\mu}(u_k,w_k)
< \inf_{(u,w)\in\mathcal{N}_{\lambda,\mu}^{+}\cup \mathcal{N}_{\lambda,\mu}^{0}}
J_{\lambda,\mu}(u,w) +\frac{1}{k},\label{c1} \\
\begin{aligned}
J_{\lambda,\mu}(u,w) & \geq J_{\lambda,\mu}(u_k,w_k) -\frac{1}{k}
\|(u-u_k,w-w_k)\|  \\
&\quad \text{for all } (u,w)\in \mathcal{N}_{\lambda,\mu}^{+} \cup \mathcal{N}_{\lambda,\mu}^{0}
\end{aligned} \label{c2}.
\end{gather}

\begin{lemma}
The sequence $\{(u_k,w_k)\}$ is bounded in $\mathcal{N}_{\lambda,\mu}$.
Moreover, there exists $0\not\equiv (u,w) \in  Y$ such that 
$(u_k,w_k)\rightharpoonup (u,w)$ weakly in $ Y$.
\end{lemma}

\begin{proof}
From equations \eqref{a1} and \eqref{c1}, we have
\[
c t^2 - d t^{1-q} = \rho(t)\leq J_{\lambda,\mu}(u,w)
<\inf_{(u,w)\in\mathcal{N}_{\lambda,\mu}^{+}
\cup \mathcal{N}_{\lambda,\mu}^{0}} J_{\lambda,\mu}(u,w)
+\frac{1}{k}\leq C_5,
\]
for sufficiently large $k$ and a suitable positive constant.
Hence, putting $t= \|(u_k,w_k)\|$ in the above equation, we obtain the 
sequence $\{(u_k,w_k)\}$ is bounded.

Let $\{(u_k,w_k)\}$ is bounded sequence in $ Y$. Then, there exists a subsequence 
of $\{(u_k,w_k)\}_k$, still denoted by $\{(u_k,w_k)\}_k$ and $(u,w)\in  Y$ 
such that $(u_k,w_k)\rightharpoonup (u,w) $  weakly in $ Y$, 
$(u_k,w_k)(\cdot)\to (u,w)(\cdot)$ strongly in $(L^{r}(\Omega))^2$ for
 $1\leq  r< 2^{*}_s$ and $u_k(\cdot)\to u(\cdot)$, $w_k(\cdot)\to w(\cdot)$
a.e. in $\Omega$.

 For any $(u,w)\in \mathcal{N}_{\lambda,\mu}^{+}$,  from $0<q<1$, $2<\alpha+\beta<2^*_s$ we have
\begin{align*}
J_{\lambda,\mu}(u,w)
&= \Big(\frac{1}{2} -\frac{1}{1-q}\Big)\|(u,w)\|^2 
 + \Big(\frac{1}{1-q} -\frac{1}{\alpha+\beta}\Big)\int_{\Omega} b(x)u_{+}^{\alpha}w_{+}^{\beta}\,dx\\
&< \Big(\frac{1}{2} -\frac{1}{1-q}\Big)\|(u,w)\|^2 
 + \Big(\frac{1}{1-q} -\frac{1}{\alpha+\beta}\Big)\frac{1+q}{\alpha+\beta-1+q}\|(u,w)\|^2\\
&= \Big(\frac{1}{\alpha+\beta} -\frac{1}{2}\Big)\frac{(1+q)}{(1-q)}\|(u,w)\|^2<0,
\end{align*}
which means that $\inf_{\mathcal{N}_{\lambda,\mu}^{+}} J_{\lambda,\mu}<0$.
Now for $(\lambda,\mu)\in\Gamma$, we know from Lemma \ref{le1}, that
$\mathcal{N}_{\lambda,\mu}^{0}=\{(0,0)\}$.
Together, these imply that $(u_k,w_k) \in \mathcal{N}_{\lambda,\mu}^{+}$
for $k$ large and
\[
\inf_{(u,w)\in\mathcal{N}_{\lambda,\mu}^{+}
\cup\mathcal{N}_{\lambda,\mu}^{0}}J_{\lambda,\mu}(u,w)
= \inf_{(u,w)\in\mathcal{N}_{\lambda,\mu}^{+}} J_{\lambda,\mu}(u,w)<0.
\]
Therefore, by weak lower semi-continuity of the norm,
\[
J_{\lambda,\mu}(u,w) \leq \liminf_{k\to\infty} J_{\lambda,\mu}(u_k,w_k)
= \inf_{\mathcal{N}_{\lambda,\mu}^{+}\cup\mathcal{N}_{\lambda,\mu}^{0}}J_{\lambda,\mu}<0,
\]
that is, $(u,w)\not\equiv 0$ and $(u,w)\in Y$.
\end{proof}

\begin{lemma}\label{le6}
Suppose $(u_k,w_k)\in \mathcal{N}_{\lambda,\mu}^{+}$ such that
$(u_k,w_k)\rightharpoonup (u,w)$ weakly in $ Y$. Then for $(\lambda,\mu)\in \Gamma$,
\begin{align}\label{s2}
(1+q) \int_{\Omega} (\lambda f(x)u_{+}^{1-q}+\mu g(x)w_{+}^{1-q}) dx
 - (\alpha+\beta-2) \int_{\Omega} b(x)u_{+}^{\alpha}w_{+}^{\beta} dx >0.
\end{align}
Moreover, there exists a constant $C_2>0$ such that
\begin{equation} \label{s5}
(1+q) \|(u_k,w_k)\|^2 - (\alpha+\beta-1+q) \int_{\Omega}  b(x)(u_k)_{+}^{\alpha}(w_k)_{+}^{\beta} \geq C_2 > 0.
\end{equation}
\end{lemma}

\begin{proof}
For $\{(u_k,w_k)\}\subset \mathcal{N}_{\lambda,\mu}^{+}(\subset \mathcal{N}_{\lambda,\mu})$,
 we have
\begin{align*}
&(1+q)  K_{\lambda,\mu}(u,w)  - (\alpha+\beta-2) \int_{\Omega}  b(x)u_{+}^{\alpha}w_{+}^{\beta}\,dx\\
&= \lim_{k\to\infty} \Big[(1+q) K_{\lambda,\mu}(u_k,w_k) - (\alpha+\beta-2) \int_{\Omega}  b(x)(u_k)_{+}^{\alpha}(w_k)_{+}^{\beta} \Big]\\
&=\lim_{k\to\infty}\Big[(1+q)\|(u_k,w_k)\|^2 - (\alpha+\beta-1+q) \int_{\Omega} b(x)(u_k)_{+}^{\alpha}(w_k)_{+}^{\beta}\Big]\geq 0.
\end{align*}
Now, we can argue by a contradiction and assume that
\begin{align}\label{s3}
(1+q)  K_{\lambda,\mu}(u,w) - (\alpha+\beta-2) \int_{\Omega}  b(x)u_{+}^{\alpha}w_{+}^{\beta}\,dx = 0.
\end{align}
Using $(u_k,w_k)\in \mathcal{N}_{\lambda,\mu}$, the weak lower semi
continuity of norm and \eqref{s3} we have that
\begin{align*}
0&= \lim_{k\to\infty} \left[\|(u_k,w_k)\|^2-  K_{\lambda,\mu}(u_k,w_k) - \int_{\Omega}  b(x)(u_k)_{+}^{\alpha}(w_k)_{+}^{\beta} \right]\\
&\geq  \|(u,w)\|^2-  K_{\lambda,\mu}(u,w)
 -  \int_{\Omega}  b(x)u_{+}^{\alpha}w_{+}^{\beta}\, dx\\
&=\begin{cases}
\|(u,w)\|^2-  \frac{\alpha+\beta-1+q}{1+q}\int_{\Omega} 
 b(x)u_{+}^{\alpha}w_{+}^{\beta}\,dx 
& \text{(solving for $K$ in \eqref{e3})},\\[4pt]
\|(u,w)\|^2-  \frac{\alpha+\beta-1+q}{\alpha+\beta-2} K_{\lambda,\mu}(u,w)
& \text{(solving for $\int_\Omega$ in \eqref{e3})}.
\end{cases}
\end{align*}
Thus for any $(\lambda,\mu)\in \Gamma$ and $(u,w)\not\equiv 0$, by similar arguments as those in \eqref{e40} we have that
\begin{align*}
0&< E_{\lambda,\mu}  \|(u,w)\|^{\alpha+\beta}\\
&\leq \frac{(1+q)}{(\alpha+\beta-1+q)} 
\Big(\frac{\alpha+\beta-2}{\alpha+\beta-1+q}\Big)^{\frac{\alpha+\beta-2}{1+q}}
 \frac{\|(u,w)\|^{\frac{2(\alpha+\beta-1+q)}{1+q}}}
{\left[ K_{\lambda,\mu}(u,w)\right]^{\frac{\alpha+\beta-2}{1+q}}} 
 - \int_{\Omega}  b(x)u_{+}^{\alpha}w_{+}^{\beta}\,dx \\
&= \frac{(1+q)}{(\alpha+\beta-1+q)} 
 \Big(\frac{\alpha+\beta-2}{\alpha+\beta-1+q}\Big)^{\frac{\alpha+\beta-2}{1+q}}
 \frac{\|(u,w)\|^{ \frac{2(\alpha+\beta-1+q)}{1+q}}}
 {\left(\frac{\alpha+\beta-2}{\alpha+\beta-1+q}\|(u,w)\|^2\right)
 ^{\frac{\alpha+\beta-2}{1+q}}}  \\
&\quad -\frac{(1+q)}{(\alpha+\beta-1+q)}\|(u,w)\|^2
=0,
\end{align*}
which is clearly impossible. Now by \eqref{s2}, we have that
\begin{align}\label{s4}
(1+q)   K_{\lambda,\mu}(u_k,w_k) - (\alpha+\beta-2) \int_{\Omega}  b(x)(u_k)_{+}^{\alpha}(w_k)_{+}^{\beta} \geq C_2
\end{align}
for sufficiently large $k$ and a suitable positive constant $C_2$. 
Then \eqref{s4}, together with the fact that $(u_k,w_k)\in \mathcal{N}_{\lambda,\mu}$
we obtain equation \eqref{s5}.
\end{proof}

 Fix $(\phi,\psi)\in C_{ Y}$ with $\phi,\psi\geq 0$. 
Then we apply Lemma \ref{le4} with $(u_k,w_k)\in \mathcal{N}_{\lambda,\mu}^{+}$
($k$ large enough such that $\frac{(1-q)C_1}{k}<C_2$), we obtain a 
sequence of functions $f_k: B{\epsilon_k}(0)\subset Y\to \mathbb{R}$ such that
  $f_k(0,0)=1$ and
$f_k(s_1,s_2)(u_k+s_1\phi,w_k+ s_2\psi)\in \mathcal{N}_{\lambda,\mu}^{+}$ for all
$s=(s_1,s_2)\in B_{\epsilon_k}(0)$. 
It follows from $(u_k,w_k)\in \mathcal{N}_{\lambda,\mu}$ and
$f_k(s_1,s_2)(u_k+s_1\phi,w_k+ s_2\psi)\in \mathcal{N}_{\lambda,\mu}$ that
\begin{align}\label{f1}
\|(u_k,w_k)\|^2 -  K_{\lambda,\mu}(u_k,w_k) - \int_{\Omega} b(x)(u_k)_{+}^{\alpha}(w_k)_{+}^{\beta} dx=0
\end{align}
and
\begin{equation} \label{f2}
\begin{aligned}
& f_k^{2}(s_1,s_2)\|(u_k+s_1\phi,w_k+s_2\psi)\|^2
 - f_k^{1-q}(s_1,s_2) K(u_k+s_1\phi,w_k+s_2\psi)  \\
&- f_k^{\alpha+\beta}(s_1,s_2)
 \int_{\Omega} b(x)(u_k+s_1\phi)_{+}^{\alpha}(w_k+s_2\psi)_{+}^{\beta} dx =0.
\end{aligned}
\end{equation}
 Choose $0<\rho <\epsilon_k$, and $(s_1,s_2)=(\rho v_1,\rho v_2)$  with $\|v\|<1$
then we find $f_k(v_1,v_2)$ such that $f_k(0,0)=1$ and
$f_k(v_1,v_2)(u_k+v_1\phi,w_k+ v_2\psi)\in \mathcal{N}_{\lambda,\mu}^{+}$
for all $v\in B_{\rho}(0)$.

\begin{lemma}\label{le7}
For $(\lambda,\mu)\in \Gamma$ we have $|\langle f'_k(0,0), (v_1,v_2) \rangle|$
is finite for every $0\leq v=(v_1,v_2)\in C_{ Y}$ with $\|v\|\leq1$.
\end{lemma}

\begin{proof}
From \eqref{f1} and \eqref{f2} we have
\begin{align*}
0&=[f_k^2(\rho v_1,\rho v_2) -1]\|(u_k+\rho v_1\phi,w_k+ \rho v_2\psi)\|^2 \\
&\quad +\|(u_k+\rho v_1\phi,w_k+\rho v_2\psi)\|^2-\|(u_k,w_k)\|^2\\
&\quad - \big[f_k^{1-q}(\rho v_1,\rho v_2) -1\big]
 \int_{\Omega} \left(\lambda f(x)(u_k+\rho v_1\phi)_{+}^{1-q} 
 +\mu g(x)(w_k+\rho v_2\psi)_{+}^{1-q}\right)dx \\
&\quad- \lambda\int_{\Omega} f(x)\left[(u_k+ \rho v_1\phi)_{+}^{1-q}
 - (u_k)_{+}^{1-q}\right]dx \\
&\quad - \mu\int_{\Omega} g(x)\left[(w_k+\rho v_2\psi)_{+}^{1-q}
 - (w_k)_{+}^{1-q}\right]dx\\
&\quad- \big[f_k^{\alpha+\beta}(\rho v_1,\rho v_2) -1\big]
 \int_{\Omega} b(x)(u_k+\rho v_1\phi)_{+}^{\alpha}(x)(w_k+\rho v_2\psi)_{+}^{\beta}(x) dx \\
&\quad- \int_{\Omega} b(x)[(u_k+\rho v_1\phi)^{\alpha}_{+}(w_k+\rho v_2\psi)_{+}^{\beta}(x)
 - (u_k)_{+}^{\alpha}(w_k)_{+}^{\beta}(x)]dx,\\
&\leq  \big[f_k^2(\rho v_1,\rho v_2) -1\big]\|(u_k+\rho v_1\phi,w_k
 + \rho v_2\psi)\|^2 +\|(u_k+\rho v_1\phi,w_k+\rho v_2\psi)\|^2 \\
&\quad -\|(u_k,w_k)\|^2
 - \left[f_k^{1-q}(\rho v_1,\rho v_2) -1\right]
 \int_{\Omega} \Big(\lambda f(x)(u_k+\rho v_1\phi)_{+}^{1-q}(x) \\
&\quad  +\mu g(x)(w_k+\rho v_2\psi)_{+}^{1-q}(x)\Big)dx \\
&\quad - \big[f_k^{\alpha+\beta}(\rho v_1,\rho v_2) -1\big]
 \int_{\Omega} b(x)(u_k+\rho v_1\phi)_{+}^{\alpha}(x)(w_k+\rho v_2\psi)_{+}^{\beta}(x)dx \\
&\quad - \int_{\Omega} b(x)\left[(u_k+\rho v_1\phi)_{+}^{\alpha}
 (w_k+\rho v_2\psi)_{+}^{\beta}(x)- (u_k)_{+}^{\alpha}(w_k)_{+}^{\beta}(x)\right]dx.
\end{align*}
Since
\begin{equation} \label{r1}
\begin{aligned}
&(u_k+\rho v_1\phi)_{+}^{1-q}(x)- (u_k)_{+}^{1-q}(x) \\
&=\begin{cases}
 (u_k+\rho v_1\phi)^{1-q}(x)- (u_k)^{1-q}(x)& \text{if } u_k\geq 0\\
 0 &\text{if } u_k\leq 0, u_k+ \rho v_1\phi\leq 0 \\
 (u_k+\rho v_1\phi)^{1-q}(x) &\text{if } u_k\leq 0, u_k+\rho v_1\phi\geq 0,
\end{cases}
\end{aligned}
\end{equation}
we have
\[
\int_{\Omega} f(x)[((u_k+\rho v_1\phi)_{+}^{1-q}- (u_k)_{+}^{1-q})(x)]dx\geq 0.
\]
Similarly, one can see that
\[
\int_{\Omega} g(x)[((w_k+\rho v_2\psi)_{+}^{1-q}- (w_k)_{+}^{1-q})(x)]dx\geq 0.
\]
 Now dividing by $\rho>0$ and passing to the limit $\rho\to 0$, we derive that
\begin{equation} \label{s6}
 \begin{aligned}
0&\leq\langle f_k'(0,0),(v_1,v_2) \rangle\Big[ 2 \|(u_k,w_k)\|^2 - (1-q)  K_{\lambda,\mu}(u_k,w_k) \\
&\quad -(\alpha+\beta)\int_{\Omega} b(x)(u_k)_{+}^{\alpha}(w_k)_{+}^{\beta} dx\Big] \\
&\quad + 2 \int_{Q}
\Big[(u_k(x)-u_k(y))\big((v_1\phi)(x)-(v_1\phi)(y)\big) \\
&\quad +(w_k(x)-w_k(y))
((v_2\psi)(x)-(v_2\psi)(y))\Big]
\big/ |x-y|^{n+2s} \,dx\,dy \\
&\quad -\alpha \int_{\Omega} b(x)(u_k)_{+}^{\alpha-1} (w_k)_{+}^{\beta} v_1\phi dx
 -\beta \int_{\Omega} b(x) (u_k)_{+}^{\alpha}(w_k)_{+}^{\beta-1} v_2\psi dx \\
&= \langle f_k'(0,0),(v_1,v_2) \rangle\Big[ (1+q)\|(u_k,w_k)\|^2  \\
&\quad -(\alpha+\beta-1+q)\int_{\Omega} b(x)(u_k)_{+}^{\alpha}(w_k)_{+}^{\beta}dx\Big]  \\
&\quad + 2 \int_{Q}
 \Big[(u_k(x)-u_k(y))((v_1\phi)(x)-(v_1\phi)(y)) \\
&\quad +(w_k(x)-w_k(y))((v_2\psi)(x)-(v_2\psi)(y))\Big]
\big/|x-y|^{n+2s} \,dx\,dy \\
&\quad -\alpha \int_{\Omega} b(x)(u_k)_{+}^{\alpha-1} (w_k)_{+}^{\beta} v_1\phi dx
 -\beta \int_{\Omega} b(x) (u_k)_{+}^{\alpha}(w_k)_{+}^{\beta-1} v_2\psi dx.
\end{aligned}
\end{equation}
From \eqref{s5} and \eqref{s6}, we know immediately that $\langle f_k'(0,0),(v_1,v_2) \rangle\ne -\infty$. Now we show that $\langle f_k'(0,0),(v_1,v_2) \rangle\ne +\infty$. Arguing by contradiction, we assume that $\langle f_k'(0,0),(v_1,v_2) \rangle=+\infty$. Since
\begin{align}\label{s7}
|f_k(\rho v_1,\rho v_2)-1|&\|(u_k, w_k)\| + \rho f_k(\rho v_1,\rho v_2) \|(v_1 \phi, v_2\psi)\| \\
& \geq \|[f_k(\rho v_1,\rho v_2)-1]\|(u_k,w_k\| + f_k(\rho v_1,\rho v_2)\|(\rho v_1 \phi, \rho v_2\psi)\| \\
&= \|f_k(\rho v_1,\rho v_2) (u_k+\rho v_1\phi, w_k +\rho v_2\psi)- (u_k,w_k)\|
\end{align}
and
\[
f_k(\rho v_1,\rho v_2)> f_k(0,0) =1
\]
for sufficiently large $k$. From the definition of derivative
$\langle f_k'(0,0),(v_1,v_2) \rangle$, applying equation \eqref{c2}
with $(u,w)= f_k(\rho v_1,\rho v_2)(u_k+\rho v_1\phi, w_k+\rho v_2\psi)
\in \mathcal{N}_{\lambda,\mu}^{+}$, we clearly have
\begin{align*}
&[f_k(\rho v_1,\rho v_2)-1]\frac{\|(u_k,w_k)\|}{k}
 + \rho f_k(\rho v_1,\rho v_2)\frac{ \|(v_1 \phi, v_2\psi)\|}{k}\\
& \geq \frac{1}{k}\|f_k(\rho v_1,\rho v_2)(u_k +\rho v_1 \phi, w_k
 +\rho v_2 \psi)-(u_k,w_k)\|\\
&\geq J_{\lambda,\mu}(u_k,w_k) - J_{\lambda,\mu}(f_k(\rho v_1,\rho v_2)
 (u_k +\rho v_1 \phi, w_k+ \rho v_2\psi))\\
&= \Big(\frac{1}{2}-\frac{1}{1-q}\Big) \|(u_k,w_k)\|^2
 + \Big(\frac{1}{1-q}-\frac{1}{2}\Big) f_k^2(\rho v_1,\rho v_2)
  \|(u_k+\rho v_1\phi, w_k+\rho v_2\psi)\|^2 \\
& \quad+ \Big(\frac{1}{1-q}-\frac{1}{\alpha+\beta}\Big)
 \Big(\int_{\Omega} b(x) (u_k)_{+}^{\alpha}(w_k)_{+}^{\beta} \\
&\quad -  f_k^{\alpha+\beta}(\rho v_1,\rho v_2)\int_{\Omega} b(x)
 (u_k+\rho v_1\phi)_{+}^{\alpha}(w_k+\rho v_2\psi)_{+}^{\beta} dx\Big)\\
&= \frac{1+q}{1-q}\Big(\|(u_k+\rho v_1\phi, w_k+\rho v_2\psi)\|^2
  -\|(u_k,w_k)\|^2 + [f_k^2(\rho v_1,\rho v_2)-1] \\
&\quad\times   \|(u_k+\rho v_1\phi, w_k+\rho v_2\psi)\|^2\Big)\\
&\quad - \Big(\frac{1}{1-q}-\frac{1}{\alpha+\beta}\Big)
  f_k^{\alpha+\beta}(\rho v_1,\rho v_2) \int_{\Omega} b(x)
 \Big[(u_k+\rho v_1\phi)_{+}^{\alpha}(w_k +\rho v_2\psi)_{+}^{\beta} \\
&\quad  - (u_k)_{+}^{\alpha}(w_k)_{+}^{\beta}\Big] dx\\
&\quad - \Big(\frac{1}{1-q}-\frac{1}{\alpha+\beta}\Big)
 \big[f_k^{\alpha+\beta}(\rho v_1,\rho v_2)-1\big]
  \int_{\Omega} b(x)(u_k)_{+}^{\alpha}(w_k)_{+}^{\beta} dx.
\end{align*}
Dividing by $\rho>0$ and passing to the limit as $\rho\to 0$, we  obtain
\begin{align*}
&\langle f_k'(0,0),(v_1,v_2) \rangle\frac{\|(u_k,w_k)\|}{k}
  + \frac{\|(v_1\phi, v_2\psi)\|}{k}\\
&\geq \Big(\frac{1+q}{1-q}\Big)\langle f_k'(0,0),(v_1,v_2) \rangle\|(u_k,w_k)\|^2\\
&\quad - \Big(\frac{\alpha+\beta-1+q}{1-q}\Big)
 \langle f_k'(0,0),(v_1,v_2) \rangle
 \int_{\Omega} b(x) (u_k)_{+}^{\alpha} (w_k)_{+}^{\beta}\\
&+ \Big(\frac{1+q}{1-q}\Big)\int_{Q}
 \Big(u_k(x)-u_k(y))((v_1\phi)(x)-(v_1\phi)(y))\\
&\quad +(w_k(x)-w_k(y)) ((v_2\psi)(x)-(v_2\psi)(y))\Big]
 \big/ |x-y|^{n+2s} \,dx\,dy\\
& - \Big(\frac{\alpha+\beta-1+q}{1-q}\Big)
 \Big(\frac{\alpha}{\alpha+\beta} \int_{\Omega} b(x)(u_k)_{+}^{\alpha-1}
 (w_k)_{+}^{\beta} v_1\phi  dx  \\
&\quad  + \frac{\beta}{\alpha+\beta}\int_{\Omega} b(x)(u_k)_{+}^{\alpha}
  (w_k)_{+}^{\beta-1}  v_2\psi dx \Big) \\
&= \frac{\langle f_k'(0,0),(v_1,v_2) \rangle}{1-q}
 \Big[ (1+q)\|(u_k,w_k)\|^2 - (\alpha+\beta-1+q) \int_{\Omega}  b(x)(u_k)_{+}^{\alpha}(w_k)_{+}^{\beta} \Big]\\
&\quad +\Big(\frac{1+q}{1-q}\Big)\int_{Q}
 \Big[(u_k(x)-u_k(y))((v_1\phi)(x)-(v_1\phi)(y)) \\
&\quad  +(w_k(x)-w_k(y))((v_2\psi)(x)-(v_2\psi)(y))\Big]
\big/ |x-y|^{n+2s} \,dx\,dy \\
&- \Big(\frac{\alpha+\beta-1+q}{1-q}\Big)
 \Big[\frac{\alpha}{\alpha+\beta} \int_{\Omega} b(x)
 (u_k)_{+}^{\alpha-1} (w_k)_{+}^{\beta}v_1\phi dx\\
&\quad  + \frac{\beta}{\alpha+\beta}\int_{\Omega} b(x) (u_k)_{+}^{\alpha}
 (w_k)_{+}^{\beta-1}v_2\psi dx\Big];
\end{align*}
that is,
 \begin{align}
&\frac{\|(v_1\phi, v_2\psi)\|}{k} \nonumber \\
&\geq  \frac{\langle f_k'(0,0),(v_1,v_2) \rangle}{1-q}
\Big[ (1+q)\|(u_k,w_k)\|^2  \nonumber \\
&\quad - (\alpha+\beta-1+q) \int_{\Omega}  b(x)(u_k)_{+}^{\alpha}(w_k)_{+}^{\beta}  
 - \frac{(1-q)\|(u_k,w_k)\|}{k} \Big]  \nonumber \\
& +\Big(\frac{1+q}{1-q}\Big) \int_{Q}
\Big[(u_k(x)-u_k(y))((v_1\phi)(x)-(v_1\phi)(y))  \nonumber \\
&\quad +(w_k(x)-w_k(y))((v_2\psi)(x)-(v_2\psi)(y))\Big]
\big/ |x-y|^{n+2s} \,dx\,dy  \label{s8}\\
&- \Big(\frac{\alpha+\beta-1+q}{1-q}\Big)
 \Big[ \frac{\alpha}{\alpha+\beta}\int_{\Omega} b(x) (u_k)_{+}^{\alpha-1}
 (w_k)_{+}^{\beta}v_1\phi dx \nonumber \\
&\quad + \frac{\beta}{\alpha+\beta}\int_{\Omega} b(x) (u_k)_{+}^{\alpha}
 (w_k)_{+}^{\beta-1}v_2\psi\Big]   \nonumber
\end{align}
which is impossible because $\langle f_k'(0,0),(v_1,v_2) \rangle=+\infty$ and
\begin{align*}
&(1+q)\|(u_k,w_k)\|^2 -(\alpha+\beta-1+q)\int_{\Omega}  b(x)(u_k)_{+}^{\alpha}(w_k)_{+}^{\beta}
 -\frac{(1-q)\|(u_k,w_k)\|}{k}  \\
&\geq C_2 - \frac{(1-q)C_1}{k}>0.
\end{align*}
 In conclusion, $|\langle f_k'(0,0),(v_1,v_2) \rangle|<+\infty$.
 Furthermore \eqref{s5} with $\|(u_k,w_k)\|\leq C_1$ and two inequalities
\eqref{s6} and \eqref{s8} also imply that
\[
|\langle f_k'(0,0),(v_1,v_2) \rangle|\leq C_3
\]
for $k$ sufficiently large and a suitable constant $C_3$.
\end{proof}

\begin{lemma}\label{le8}
For each $0\leq (\phi,\psi)\in C_{ Y}$ and for every $0\leq v=(v_1,v_2)\in Y$ 
with $\|v\|\leq 1$, we have 
$\lambda f(x)u_{+}^{-q} v_1\phi+\mu g(x) w_{+}^{-q} v_2\psi \in L^{1}(\Omega)$ 
and
\begin{equation} \label{s9}
 \begin{aligned}
&\int_{Q}\frac{(u(x)-u(y))((v_1\phi)(x)-(v_1\phi)(y))}{|x-y|^{n+2s}} \,dx\,dy \\
&+ \int_{Q}\frac{(w(x)-w(y))((v_2\psi)(x)-(v_2\psi)(y))}{|x-y|^{n+2s}} \,dx\,dy \\
&- \int_{\Omega} \left(\lambda f(x)u_{+}^{-q} v_1\phi +\mu g(x) w_{+}^{-q} v_2\psi\right) dx
 - \int_{\Omega} b(x) u_{+}^{\alpha-1}v_{+}^{\beta}v_1\phi dx \\
&-\int_{\Omega} b(x)u_{+}^{\alpha} w_{+}^{\beta-1}v_2\psi dx\geq 0.
\end{aligned}
\end{equation}
\end{lemma}

\begin{proof}
Applying \eqref{s7} and $\eqref{c2}$ again, we have that
\begin{align*}
&[f_k(\rho v_1,\rho v_2)-1]\frac{\|(u_k,w_k)\|}{k}
 + \rho f_k(\rho v_1,\rho v_2)\frac{ \|(v_1 \phi, v_2\psi)\|}{k}\\
& \geq \frac{1}{k}\|f_k(\rho v_1,\rho v_2)(u_k +\rho v_1 \phi, w_k
 +\rho v_2 \psi)-(u_k,w_k)\|\\
&\geq J_{\lambda,\mu}(u_k,w_k) - J_{\lambda,\mu}(f_k(\rho v_1,\rho v_2)
 (u_k+\rho v_1 \phi, w_k+ \rho v_2\psi))\\
&= \frac{1}{2}\|(u_k,w_k)\|^2 -\frac{1}{2} \|f_k(\rho v_1,\rho v_2)
 (u_k+\rho v_1 \phi,w_k+\rho v_2\psi)\|^2 dx \\
&\quad -\frac{1}{1-q}\int_{\Omega} \left(\lambda f (u_k)_{+}^{1-q}
 +\mu g (w_k)_{+}^{1-q}\right)\\
&\quad +\frac{1}{1-q}\int_{\Omega} \Big(\lambda f(x)(f_k(\rho v_1,\rho v_2)
 (u_k+\rho v_1\phi))_{+}^{1-q} \\
&\quad +\mu g(x)(f_k(\rho v_1,\rho v_2)(w_k+\rho v_2\psi))_{+}^{1-q} \Big)dx \\
&\quad + \frac{1}{\alpha+\beta} \int_{\Omega} b(x)
 \Big[(f_k(\rho v_1,\rho v_2)(u_k+\rho v_1\phi))_{+}^{\alpha}(f_k
 (\rho v_1,\rho v_2)(w_k+\rho v_2\psi))_{+}^{\beta}  \\
&\quad -(u_k)_{+}^{\alpha}(w_k)_{+}^{\beta}\Big]dx\\
&=- \frac{f^{2}_k(\rho v_1,\rho v_2)-1}{2} \|(u_k,w_k)\|^2
 - \frac{f_k^2(\rho v_1,\rho v_2)}{2}(\|(u_k+\rho v_1\phi,
  w_k+\rho v_2\psi)\|^2 \\
&\quad - \|(u_k,w_k)\|^2)\\
&\quad +\frac{f_k^{1-q}(\rho v_1,\rho v_2)-1}{1-q}
 \int_{\Omega} \left(\lambda f(x) (u_k+\rho v_1\phi)_{+}^{1-q}
 +\mu g(x) (w_k+\rho v_2\psi)_{+}^{1-q}\right)dx\\
&\quad+\frac{1}{1-q} \int_{\Omega} \Big(\lambda f(x)
 \left((u_k+\rho v_1\phi)_{+}^{1-q} - (u_k)_{+}^{1-q}\right)(x) \\
&\quad +\mu g(x) \left((w_k+\rho v_2\psi)_{+}^{1-q}
 - (w_k)_{+}^{1-q}\right)(x)\Big) dx\\
&\quad + \frac{f_k^{\alpha+\beta}(\rho v_1,\rho v_2)-1}{\alpha+\beta}
  \int_{\Omega} b(x) (u_k+\rho v_1\phi)_{+}^{\alpha}(x)(w_k+\rho v_2\psi)_{+}^{\beta}(x) dx\\
 &\quad+\frac{1}{\alpha+\beta} \int_{\Omega} b(x)
 \left[\left((u_k+\rho v_1\phi)_{+}^{\alpha}(w_k+\rho v_2\psi)_{+}^{\beta}
 - (u_k)_{+}^{\alpha}(w_k)_{+}^{\beta}\right)(x)\right]dx.
\end{align*}
 Dividing by $\rho>0$ and passing to the limit $\rho\to 0^+$, we obtain
 \begin{align*}
&|\langle f'_k(0,0),(v_1,v_2) \rangle|\frac{\|(u_k,w_k)\|}{k}
 + \frac{\|(v_1\phi,v_2\psi)\|}{k}\\
&\geq- \langle f'_k(0,0), (v_1,v_2)\rangle
 \Big[\|(u_k,w_k)\|^2 - K_{\lambda,\mu}(u_k,w_k) -\int_{\Omega} b(x)(u_k)_{+}^{\alpha}(w_k)_{+}^{\beta} dx\Big]\\
&\quad -\int_{Q}
\Big(u_k(x)-u_k(y))((v_1\phi)(x)-(v_1\phi)(y)) \\
&\quad +(w_k(x)-w_k(y))((v_2\psi)(x)-(v_2\psi)(y))\Big]
\big/ |x-y|^{n+2s}\,dx\,dy \\
&\quad+\frac{\alpha}{\alpha+\beta}\int_{\Omega} b(x) (u_k)_{+}^{\alpha-1}
 (w_k)_{+}^{\beta} v_1\phi dx + \frac{\beta}{\alpha+\beta}
 \int_{\Omega} b(x) (u_k)_{+}^{\alpha}(w_k)_{+}^{\beta-1} v_2\psi dx\\
&\quad + \frac{1}{1-q}\liminf_{\rho\to 0^{+}}
\Big[ \int_{\Omega}\frac{ \lambda f(x)((u_k+\rho v_1\phi)_{+}^{1-q}
- (u_k)_{+}^{1-q})}{\rho} \\
&\quad + \int_{\Omega}\frac{ \mu g(x)((w_k+\rho v_2\psi)_{+}^{1-q}
- (w_k)_{+}^{1-q})}{\rho}\Big]\\
&= -\int_{Q}
\Big[(u_k(x)-u_k(y))((v_1\phi)(x)-(v_1\phi)(y))\\
&\quad +(w_k(x)-w_k(y)) ((v_2\psi)(x)-(v_2\psi)(y))\Big]
\big/|x-y|^{n+2s}\,dx\,dy\\
&\quad +\frac{\alpha}{\alpha+\beta}\int_{\Omega} b(x) (u_k)_{+}^{\alpha-1}
 (w_k)_{+}^{\beta} v_1\phi dx
+ \frac{\beta}{\alpha+\beta} \int_{\Omega} b(x) (u_k)_{+}^{\alpha}(w_k)_{+}^{\beta-1} v_2\psi dx\\
&+ \frac{1}{1-q}\liminf_{\rho\to 0^{+}}
\Big[\int_{\Omega}\frac{ \lambda f(x)((u_k+\rho v_1\phi)_{+}^{1-q}
- (u_k)_{+}^{1-q})}{\rho} \\
&\quad + \int_{\Omega}\frac{ \mu g(x)((w_k+\rho v_2\psi)_{+}^{1-q}
- (w_k)_{+}^{1-q})}{\rho}\Big].
\end{align*}
Then by the above inequality, one can see that
\begin{align*}
&\liminf_{\rho\to 0^{+}}\Big[ \int_{\Omega}\frac{ \lambda f(x)((u_k+\rho v_1\phi)_{+}^{1-q}
- (u_k)_{+}^{1-q})}{\rho}dx \\
&+ \int_{\Omega}\frac{ \mu g(x)((w_k+\rho v_2\psi)_{+}^{1-q}
 - (w_k)_{+}^{1-q})}{\rho}dx\Big]
\end{align*}
is finite. Now, using \eqref{r1}, we have
\[
f(x)\left((u_k+ \rho v_1\phi)_{+}^{1-q} - (u_k)_{+}^{1-q}\right)\geq0.
\]
Similarly we have
\[
g(x)\left((w_k+\rho v_2\psi)_{+}^{1-q} - (w_k)_{+}^{1-q}\right)\geq 0,\quad
\text{for all }  x\in \Omega,\text{ for all} \;t>0.
\]
Then by Fatou's Lemma, we have
\begin{align*}
&\int_{\Omega} \left(\lambda f(x)(u_k)_{+}^{-q}v_1\phi
 +\mu g(x)(w_k)_{+}^{-q}v_2\psi\right)  dx\\
&\leq \frac{1}{1-q} \liminf_{\rho\to 0^{+}}
 \Big[\lambda\int_{\Omega}\frac{ f(x)((u_k+\rho v_1\phi)_{+}^{1-q} - (u_k)_{+}^{1-q})}{\rho}\\
&\quad + \mu \int_{\Omega}\frac{ g(x) ((w_k+\rho v_2\psi)_{+}^{1-q}
 - (w_k)_{+}^{1-q})}{\rho}dx\Big]\\
&\leq \frac{|\langle f'_k(0,0),(v_1,v_2) \rangle|\|(u_k,w_k)\|
 + \|(v_1\phi, v_2\psi)\|}{k}\\
&\quad+\int_{Q}
\Big[(u_k(x)-u_k(y))((v_1\phi)(x)-(v_1\phi)(y)) \\
&\quad +(w_k(x)-w_k(y))((v_2\psi)(x)-(v_2\psi)(y))\Big]
\big/ |x-y|^{n+2s} \,dx\,dy\\
&\quad- \frac{\alpha}{\alpha+\beta} \int_{\Omega} b(x)(u_k)_{+}^{\alpha-1}
 (w_k)_{+}^{\beta} v_1\phi dx
-\frac{\beta}{\alpha+\beta} \int_{\Omega} b(x)(u_k)_{+}^{\alpha} (w_k)_{+}^{\beta-1}
  v_2\psi dx\\
&\leq \frac{C_1C_3 \|(v_1,v_2)\|+ \|(v_1\phi, v_2\psi)\|}{k}
 - \frac{\alpha}{\alpha+\beta} \int_{\Omega} b(u_k)_{+}^{\alpha-1} (w_k)_{+}^{\beta} v_1\phi\\
&\quad -\frac{\beta}{\alpha+\beta} \int_{\Omega} b(x)(u_k)_{+}^{\alpha}
 (w_k)_{+}^{\beta-1} v_2\psi\\
&\quad+\int_{Q}
\Big[(u_k(x)-u_k(y))((v_1\phi)(x)-(v_1\phi)(y)) \\
&\quad +(w_k(x)-w_k(y))((v_2\psi)(x)-(v_2\psi)(y))\Big]
\big/ |x-y|^{n+2s}\,dx\,dy
\end{align*}
 Again using the Fatou's Lemma and the above relation, we have
\begin{align*}
 &\lambda \int_{\Omega} f(x) u_{+}^{-q} v_1\phi dx
 + \mu \int_{\Omega} g(x) w_{+}^{-q}v_2\psi dx \\
& \leq \int_{\Omega} \left[\liminf_{k\to\infty}\left(\lambda f(x) u_{+}^{-q}v_1\phi
 + \mu  g(x) w_{+}^{-q}v_2\psi\right)  \right] dx\\
&\leq  \liminf_{k\to\infty}\int_{\Omega} \left(\lambda f(x)(u_k)_{+}^{-q}v_1\phi
 +\mu g(x)(w_k)_{+}^{-q}v_2\psi\right) dx \\
&\leq \frac{C_1C_3 \|(v_1,v_2)\|+ \|(v_1\phi, v_2\psi)\|}{k}
 - \frac{\alpha}{\alpha+\beta} \int_{\Omega} b(x)(u_k)_{+}^{\alpha-1}
 (w_k)_{+}^{\beta} v_1\phi \\
&\quad -\frac{\beta}{\alpha+\beta} \int_{\Omega} b(u_k)_{+}^{\alpha} (w_k)_{+}^{\beta-1}
  v_2\psi\\
&\quad+\int_{Q}\Big[(u_k(x)-u_k(y))((v_1\phi)(x)-(v_1\phi)(y)) \\
&\quad +(w_k(x)-w_k(y))((v_2\psi)(x)-(v_2\psi)(y))\Big]
\big/|x-y|^{n+2s}\,dx\,dy
\end{align*}
which completes the proof.
\end{proof}

\begin{corollary}\label{cc1}
For every $0\leq (\phi,\psi)\in  Y$, we have 
$\left(\lambda f(x) u_{+}^{-q}\phi+\mu g(x) w_{+}^{-q}\psi\right) \in L^{1}(\Omega)$,
 $u_{+}$, $w_{+}>0$ in ${\Omega}$ and
\begin{equation} \label{s91}
\begin{aligned}
&\int_{Q}\frac{(u(x)-u(y))(\phi(x)-\phi(y))}{|x-y|^{n+2s}} \,dx\,dy\\
&+\int_{Q}\frac{(w(x)-w(y))(\psi(x)-\psi(y))}{|x-y|^{n+2s}} \,dx\,dy 
-\lambda \int_{\Omega} f(x) u_{+}^{-q}\phi dx  \\
& - \mu \int_{\Omega} g(x) w_{+}^{-q}\psi dx
 -\frac{\alpha}{\alpha+\beta} \int_{\Omega} b(x) u_{+}^{\alpha-1}w_{+}^{\beta}\phi dx \\
& -\frac{\beta}{\alpha+\beta} \int_{\Omega} b(x) u_{+}^{\alpha}w_{+}^{\beta-1}\psi dx\geq 0.
\end{aligned}
\end{equation}
\end{corollary}

\begin{proof}
Choosing $v=(v_1,v_2)\in Y$ such that $v\geq 0$, $v\equiv l$ in the neighborhood 
of support of $\phi$ and $\|v\|\leq 1$, for some $l>0$ is a constant. 
Then we note that $ \lambda\int_{\Omega} f(x) u_{+}^{-q}\phi dx 
+ \mu \int_{\Omega} g(x) w_{+}^{-q}\psi dx<\infty$, for every $0\leq (\phi,\psi)\in C_{Y}$
 which guarantees that $u_{+}$, $w_{+}>0$  a.e in $\Omega$. Putting this choice
of $v$ in \eqref{s9},  for every $0\leq (\phi,\psi)\in C_{Y}$ we have
\begin{align*}
&\int_{Q}\frac{(u(x)-u(y))(\phi(x)-\phi(y))}{|x-y|^{n+2s}} \,dx\,dy
 +\int_{Q}\frac{(w(x)-w(y))(\psi(x)-\psi(y))}{|x-y|^{n+2s}} \,dx\,dy\\
&-\lambda \int_{\Omega} f(x) u_{+}^{-q}\phi dx  
 - \mu \int_{\Omega} g(x) w_{+}^{-q}\psi dx
 -\frac{\alpha}{\alpha+\beta} \int_{\Omega} b(x) u_{+}^{\alpha-1}w_{+}^{\beta}\phi dx \\
&-\frac{\beta}{\alpha+\beta} \int_{\Omega} b(x) u_{+}^{\alpha}w_{+}^{\beta-1}\psi dx\geq 0.
\end{align*}
Hence by density argument, \eqref{s91} holds for every $0\leq (\phi,\psi)\in  Y$,
 which completes the proof.
\end{proof}

\begin{lemma}\label{le9}
We have that $u>0$, $w>0$ and $(u,w)\in \mathcal{N}_{\lambda,\mu}^{+}$.
\end{lemma}

\begin{proof}
Using \eqref{s91} with $\phi=u^-$, $\psi=w^-$, we obtain that
\begin{align*}
0&\leq \int_{Q}\frac{(u(x)-u(y))(u^-(x)- u^-(y))}{|x-y|^{n+2s}} \,dx\,dy \\
&\quad + \int_{Q}\frac{(w(x)-w(y))(w^-(x)- w^-(y))}{|x-y|^{n+2s}} \,dx\,dy \\
&\leq - \|u^-\|^2 - \|w^-\|^2 - 2\int_{Q}
 \frac{u^-(x) u^+(y)+w^-(x) w^+(y)}{|x-y|^{n+2s}}\,dx\,dy \\
& \leq - \|u^-\|^2 - \|w^-\|^2 \leq 0.
\end{align*}
i.e, $u^{-}=w^-=0$ a.e. So, $u=u^+ >0$, $w=w^+ >0$ a.e.
 by Corollary \ref{cc1}. Hence $u$, $w>0$ in $\Omega$.
Now using \eqref{s91} with $\phi=u$, $\psi=w$, we obtain that
\begin{align*}
\|(u,w)\|^2\geq  \int_{\Omega} \left(\lambda f(x) u_{+}^{1-q}+\mu g(x) w_{+}^{1-q}\right) dx 
+ \int_{\Omega}  b(x)u_{+}^{\alpha}w_{+}^{\beta}\,dx.
\end{align*}
On the other hand, by the weak lower semi-continuity of the norm, we have 
\begin{align*}
\|(u,w)\|^2
&\leq \liminf_{k\to\infty}\|(u_k,w_k)\|^2
 \leq \limsup_{k\to\infty}\|(u_k,w_k)\|^2\\
&= \int_{\Omega} \left(\lambda f(x) u_{+}^{1-q}+\mu g(x) w_{+}^{1-q}\right) dx 
+ \int_{\Omega}  b(x)u_{+}^{\alpha}w_{+}^{\beta}\,dx.
\end{align*}
Thus
\begin{equation} \label{s10}
\|(u,w)\|^2
= \int_{\Omega} \left(\lambda f(x) u_{+}^{1-q}+\mu g(x) w_{+}^{1-q}\right)dx
+  \int_{\Omega}  b(x)u_{+}^{\alpha}w_{+}^{\beta}\,dx.
\end{equation}
Consequently, $(u_k,w_k)\to (u,w)$ in $ Y$ and $(u,w)\in \mathcal{N}_{\lambda,\mu}$.
Now from \eqref{s2} it follows that
\begin{align*}
&(1+q)\|(u,w)\|^2 -(\alpha+\beta-1+q)\int_{\Omega}  b(x)u_{+}^{\alpha}w_{+}^{\beta}\,dx\\
=& (1+q)\int_{\Omega} \left(\lambda f(x) u_{+}^{1-q}+\mu g(x) w_{+}^{1-q}\right) dx
- (\alpha+\beta-2) \int_{\Omega}  b(x)u_{+}^{\alpha}w_{+}^{\beta}\,dx>0,
\end{align*}
that is, $(u,w)\in \mathcal{N}_{\lambda,\mu}^{+}$.
\end{proof}

\begin{lemma}\label{le10}
the pair $(u,w)$ is  a positive weak solution of problem \eqref{ePlm}.
\end{lemma}

\begin{proof}
Let $(u,w)=(u_1,u_2)$, $(\phi_1,\phi_2)\in Y$ and $\epsilon>0$, then we define
 \[
\Psi(x)=(\Psi_1,\Psi_2)= ((u_1+\epsilon\phi_1)_{+},(u_2+\epsilon\phi_2)_{+})
\]
For $i=1,2$, let $\Omega =\Omega_i\times\Gamma_i$  with
\[
\Omega_i:= \{x\in \Omega: u_i(x)+ \epsilon \phi_i(x)>0\}, \quad
\Gamma_i:= \{x\in\Omega: u_i(x)+ \epsilon \phi_i(x)\leq0\}.
\]
Then $\Psi_i|_{\Omega_i}(x)= (u_i+\epsilon\phi_i)_{+}(x)$, and 
$\Psi_i|_{\Gamma_i}(x)= 0$. Decompose
\begin{align*}
Q :=& (\Omega_{i}\times \Omega^c)\cup (\Gamma_i\times \Omega^c)\cup(\Omega^c\times \Omega_i) 
  \cup(\Omega^c\times \Gamma_i)\\
&\cup(\Gamma_i\times \Omega_i)\cup(\Omega_i\times \Gamma_i)
 \cup(\Omega_i\times \Omega_i)\cup(\Gamma_i\times \Gamma_i).
\end{align*}
Let $ M_i(x,y)= u_i(x,y)((u_i+\epsilon\phi_i)^{-}(x)- (u_i+\epsilon\phi_i)^{-}(y))K(x,y)$,
 where $u_i(x,y)=(u_i(x)-u_i(y))$ and $K(x,y)=\frac{1}{|x-y|^{n+2s}}$. Then we have
\begin{itemize}
\item[(1)] $\int_{\Omega_i\times \Omega^c} M_i(x,y)\,dx\,dy 
 = \int_{\Omega^c\times \Omega_i} M_i(x,y)\,dx\,dy= 0$,
\item[(2)] $\int_{\Gamma_i\times \Omega^c} M_i(x,y)\,dx\,dy 
 = -\int_{\Gamma_i\times \Omega^c} u_i(x)(u_i+\epsilon\phi_i)(x)K(x,y)\,dx\,dy$,
\item[(3)] $\int_{\Omega^c\times \Gamma_i} M_i(x,y)\,dx\,dy
  = -\int_{\Omega^c\times \Gamma_i} u_i(x)(u_i+\epsilon\phi_i)(x)K(x,y)\,dx\,dy$,
\item[(4)] $\int_{\Gamma_i\times \Omega_i} M_i(x,y)\,dx\,dy 
 = - \int_{\Gamma_i\times \Omega_i}u_i(x,y)(u_i+\epsilon\phi_i)(x)K(x,y)\,dx\,dy $,
\item[(5)] $\int_{\Omega_i\times \Gamma_i} M_i(x,y)\,dx\,dy 
= - \int_{\Omega_i\times \Gamma_i}u_i(x,y)(u_i+\epsilon\phi_i)(x)K(x,y)\,dx\,dy$,
\item[(6)] $\int_{\Omega_i\times \Omega_i} M_i(x,y)\,dx\,dy = 0$,
\item[(7)] $\int_{\Gamma_i\times \Gamma_i} M_i(x,y)\,dx\,dy 
= - \int_{\Gamma_i\times \Gamma_i}u_i(x,y)((u_i+\epsilon\phi_i)(x) \\
- (u_i+\epsilon\phi_i)(y))K(x,y)\,dx\,dy$.
\end{itemize}
 Now relabeling $(\Psi_1,\Psi_2)=(\Phi,\Psi)$, $(u_1,u_2)=(u,w)$ and 
$(\phi_1,\phi_2)=(\phi,\psi)$. Then putting $(\Phi,\Psi)$ into \eqref{s9} 
and using \eqref{s10}, we see that
\begin{align*}
0&\leq  \int_{Q}\frac{u(x,y)(\Phi(x)-\Phi(y))+w(x,y)(\Psi(x) 
 -\Psi(y))}{|x-y|^{n+2s}}\,dx\,dy \\
&\quad - \int_{\Omega} \left(\lambda f(x) u_{+}^{-q}\Phi+ \mu g(x) w_{+}^{-q}\Psi\right) dx\\
&- \frac{\alpha}{\alpha+\beta}\int_{\Omega} b(x) u_{+}^{\alpha-1}w_{+}^{\beta}\Phi dx
 -\frac{\beta}{\alpha+\beta}\int_{\Omega} b(x) u_{+}^{\alpha}w_{+}^{\beta-1}\Psi dx\\
&= \int_{Q}\Big[u(x,y)((u+\epsilon\phi)(x)- (u+\epsilon\phi)(y)) \\
&\quad  + w(x,y)((w+\epsilon\psi)(x)- (w+\epsilon\psi)(y))\Big]
\big/ |x-y|^{n+2s}\,dx\,dy\\
&\quad+\int_{Q}\Big[u(x,y)((u+\epsilon\phi)^{-}(x)- (u+\epsilon\phi)^{-}(y))\\
&\quad + w(x,y)((w+\epsilon\psi)^{-}(x)- (w+\epsilon\psi)^{-}(y))\Big]
\big/ |x-y|^{n+2s}\,dx\,dy\\
&\quad - \int_{\Omega}(\lambda f(x) u_{+}^{-q}(u+\epsilon\phi)+\mu g(x) w_{+}^{-q}(w+\epsilon\psi))\\
&\quad - \int_{\Omega}(\lambda f(x) u_{+}^{-q}(u+\epsilon\phi)^{-}+ \mu g(x) 
 w_{+}^{-q}(w+\epsilon\psi)^{-})\\
&\quad - \frac{\alpha}{\alpha+\beta} 
 \int_{\Omega} b(x) u_{+}^{\alpha-1}w_{+}^{\beta}(u+\epsilon\phi) dx
 -\frac{\beta}{\alpha+\beta}\int_{\Omega} b(x) u_{+}^{\alpha}w_{+}^{\beta-1}(w+\epsilon\psi) dx\\
 &- \frac{\alpha}{\alpha+\beta}\int_{\Omega} b(x) u_{+}^{\alpha-1}w_{+}^{\beta}(u+\epsilon\phi)^{-} dx
-\frac{\beta}{\alpha+\beta}\int_{\Omega} b(x) u_{+}^{\alpha}w_{+}^{\beta-1}(w+\epsilon\psi)^{-} dx\\
&=\epsilon\Big(\int_{Q}\frac{u(x,y)(\phi(x)- \phi(y))+w(x,y)(\psi(x)
 - \psi(y))}{|x-y|^{n+2s}} \,dx\,dy \\
&\quad - \int_{\Omega} (\lambda f(x) u_{+}^{-q}\phi + \mu g(x) w_{+}^{-q}\psi)dx 
 - \frac{\alpha}{\alpha+\beta}\int_{\Omega} b(x) u_{+}^{\alpha-1}w_{+}^{\beta}\phi \\
&\quad -\frac{\beta}{\alpha+\beta}\int_{\Omega} b(x)
u_{+}^{\alpha}w_{+}^{\beta-1}\psi\Big) \\
&\quad +\int_{Q}\frac{|u(x)-u(y)|^{2}+|w(x)-w(y)|^{2}}{|x-y|^{n+2s}} \,dx\,dy\\
 &\quad +\int_{Q}
 \Big[u(x,y)((u+\epsilon\phi)^{-}(x)- (u+\epsilon\phi)^{-}(y)) \\
&\quad +w(x,y)((w+\epsilon\psi)^{-}(x) - (w+\epsilon\psi)^{-}(y))\Big]
\big/ |x-y|^{n+2s}\,dx\,dy \\
& - \int_{\Omega} (\lambda f(x) u_{+}^{1-q}+ \mu g(x) w_{+}^{1-q}) dx
 +\lambda \int_{\Gamma_1} f(x) u_{+}^{-q}(u+\epsilon\phi)dx \\
&\quad + \mu \int_{\Gamma_2}g(x) w_{+}^{-q}(w+\epsilon\psi) dx
+\frac{\alpha}{\alpha+\beta} \int_{\Gamma_1} 
 b(x) u_{+}^{\alpha-1}w_{+}^{\beta}(u+\epsilon\phi)\, dx \\
&\quad +\frac{\beta}{\alpha+\beta} \int_{\Gamma_2} b(x) u_{+}^{\alpha}w_{+}^{\beta-1}
 (w+\epsilon\psi)\, dx
 -  \int_{\Omega} b(x) u_{+}^{\alpha}w_{+}^{\beta} dx \\
&=\epsilon\Big(\int_{Q}\frac{u(x,y)(\phi(x)- \phi(y))+w(x,y)(\psi(x)
 - \psi(y))}{|x-y|^{n+2s}} \,dx\,dy \\
&\quad - \int_{\Omega} (\lambda f(x) u_{+}^{-q}\phi + \mu g(x) w_{+}^{-q}\psi)dx  
 - \frac{\alpha}{\alpha+\beta}\int_{\Omega} b(x) u_{+}^{\alpha-1}w_{+}^{\beta}\phi \\
&\quad -\frac{\beta}{\alpha+\beta}\int_{\Omega} b(x) u_{+}^{\alpha}w_{+}^{\beta-1}\psi \Big)
 -2\int_{\Gamma_1\times \Omega^c}\frac{u(x)(u+\epsilon\phi)(x)}{|x-y|^{n+2s}}\\
&\quad -2\int_{\Gamma_2\times \Omega^c}\frac{w(x)(w+\epsilon\psi)(x)}{|x-y|^{n+2s}}\,dx\,dy
 - 2\int_{\Gamma_1\times \Omega_1}\frac{u(x,y)(u+\epsilon\phi)(x)}{|x-y|^{n+2s}}\,dx\,dy\\
&\quad - 2\int_{\Gamma_2\times \Omega_2}\frac{w(x,y)(w+\epsilon\psi)(x)}{|x-y|^{n+2s}}\,dx\,dy \\
&\quad -2\int_{\Gamma_1\times\Gamma_1}\frac{u(x,y)((u+\epsilon\phi)(x)-
(u+\epsilon\phi)(y))}{|x-y|^{n+2s}}\\
&\quad -2\int_{\Gamma_2\times\Gamma_2}\frac{w(x,y)((w+\epsilon\psi)(x)
 - (w+\epsilon\psi)(y))}{|x-y|^{n+2s}}\,dx\,dy \\
&\quad + \lambda\int_{\Gamma_1} f(x) u_{+}^{-q}(u+\epsilon\phi)dx
 + \mu \int_{\Gamma_2}g(x) w_{+}^{-q}(w+\epsilon\psi) \\
&\quad + \frac{\alpha}{\alpha+\beta}\int_{\Gamma_1} b(x) u_{+}^{\alpha-1}
 w_{+}^{\beta}(u+\epsilon\phi) 
 + \frac{\beta}{\alpha+\beta}\int_{\Gamma_2} b(x) u_{+}^{\alpha}w_{+}^{\beta-1}
 (w+\epsilon\psi) \\
&=\epsilon\Big(\int_{Q}\frac{u(x,y)(\phi(x)- \phi(y))
 +w(x,y)(\psi(x)- \psi(y))}{|x-y|^{n+2s}} \,dx\,dy \\
&\quad - \int_{\Omega} (\lambda f(x) u_{+}^{-q}\phi + \mu g(x) w_{+}^{-q}\psi)dx
  - \frac{\alpha}{\alpha+\beta}\int_{\Omega} b(x) u_{+}^{\alpha-1}w_{+}^{\beta}\phi dx \\
&\quad -\frac{\beta}{\alpha+\beta}\int_{\Omega} b(x) u_{+}^{\alpha}w_{+}^{\beta-1}\psi dx\Big) 
-2\int_{\Gamma_1\times \Omega^c}\frac{|u(x)|^{2}}{|x-y|^{n+2s}}\,dx\,dy\\
&\quad - 2\int_{\Gamma_2\times \Omega^c}\frac{|w(x)|^{2}}{|x-y|^{n+2s}}\,dx\,dy
 -2\int_{\Gamma_1\times \Omega_1}\frac{u(x,y)u(x)}{|x-y|^{n+2s}}\,dx\,dy \\
&\quad -2\int_{\Gamma_2\times \Omega_2}\frac{w(x,y)w(x)}{|x-y|^{n+2s}}\,dx\,dy
-2\int_{\Gamma_1\times \Gamma_1}\frac{|u(x)-u(y)|^{2}}{|x-y|^{n+2s}}\,dx\,dy  \\
&\quad -2\int_{\Gamma_2\times \Gamma_2}\frac{|w(x)-w(y)|^{2}}{|x-y|^{n+2s}}\,dx\,dy
 + \lambda\int_{\Gamma_1} f(x) u_{+}^{-q}(u+\epsilon\phi) dx\\
&\quad +\mu \int_{\Gamma_2} g(x) w_{+}^{-q}(w+\epsilon\psi) 
 + \frac{\alpha}{\alpha+\beta}\int_{\Gamma_1} b(x) u_{+}^{\alpha-1}
 w_{+}^{\beta}(u+\epsilon\phi) \\
&\quad + \frac{\beta}{\alpha+\beta}\int_{\Gamma_2} b(x) u_{+}^{\alpha}w_{+}^{\beta-1}
 (w+\epsilon\psi)  
 -2\epsilon\Big(\int_{\Gamma_1\times \Omega^c}\frac{ u(x)\phi(x)}{|x-y|^{n+2s}}\,dx\,dy \\
&\quad  +\int_{\Gamma_2\times \Omega^c}\frac{w(x)\psi(x)}{|x-y|^{n+2s}}\,dx\,dy
  + \int_{\Gamma_1\times \Omega_1}\frac{u(x,y)\phi(x)}{|x-y|^{n+2s}}\,dx\,dy \\
&\quad +  \int_{\Gamma_2\times \Omega_2}\frac{ w(x,y)\psi(x)}{|x-y|^{n+2s}}\,dx\,dy
  +\int_{\Gamma_1\times \Gamma_1}\frac{u(x,y)(\phi(x)- \phi(y))}{|x-y|^{n+2s}} \\
&\quad + \int_{\Gamma_2\times \Gamma_2}\frac{w(x,y)(\psi(x)- \psi(y))}{|x-y|^{n+2s}}
 \,dx\,dy\Big)\\
&\leq \epsilon\Big(\int_{Q}\frac{u(x,y)(\phi(x)- \phi(y))+w(x,y)(\psi(x)
  - \psi(y))}{|x-y|^{n+2s}} \,dx\,dy \\
&\quad - \int_{\Omega}\left(\lambda f(x) u_{+}^{-q}\phi + \mu g(x) w_{+}^{-q}\psi\right)dx
 - \frac{\alpha}{\alpha+\beta}\int_{\Omega} b(x) u_{+}^{\alpha-1}w_{+}^{\beta}\phi dx \\
&\quad -\frac{\beta}{\alpha+\beta}\int_{\Omega} b(x) u_{+}^{\alpha}w_{+}^{\beta-1}\psi\Big)
 -2\int_{\Gamma_1\times \Omega_1}\frac{u(x,y)u(x)}{|x-y|^{n+2s}}\,dx\,dy \\
&\quad -2\int_{\Gamma_2\times \Omega_2}\frac{w(x,y)w(x)}{|x-y|^{n+2s}}\,dx\,dy
 + \lambda\int_{\Gamma_1} f(x) u_{+}^{-q}(u+\epsilon\phi) \\
&\quad +\mu \int_{\Gamma_2} g(x) w_{+}^{-q}(w+\epsilon\psi)) dx
 + \frac{\alpha}{\alpha+\beta}\int_{\Gamma_1} b(x) u_{+}^{\alpha-1}
 w_{+}^{\beta}(u+\epsilon\phi) dx \\
&\quad + \frac{\beta}{\alpha+\beta}\int_{\Gamma_2} b(x) 
u_{+}^{\alpha}w_{+}^{\beta-1}(w+\epsilon\psi) dx\\
&\quad -2\epsilon\Big(\int_{\Gamma_1\times \Omega^c}\frac{ u(x)\phi(x)}{|x-y|^{n+2s}}\,dx\,dy
 +\int_{\Gamma_2\times \Omega^c}\frac{w(x)\psi(x)}{|x-y|^{n+2s}}\,dx\,dy \\
&\quad +\int_{\Gamma_1\times \Omega_1}\frac{u(x,y)\phi(x)}{|x-y|^{n+2s}}\,dx\,dy 
 +  \int_{\Gamma_2\times \Omega_2}\frac{ w(x,y)\psi(x)}{|x-y|^{n+2s}} \\
&\quad + \int_{\Gamma_1\times \Gamma_1}\frac{u(x,y)(\phi(x)
 - \phi(y))}{|x-y|^{n+2s}}\,dx\,dy \\
&\quad  +\int_{\Gamma_2\times \Gamma_2}\frac{w(x,y)(\psi(x)
 - \psi(y))}{|x-y|^{n+2s}}\,dx\,dy\Big)\\
&\leq \epsilon\Big(\int_{Q}\frac{(u(x)-u(y))(\phi(x)- \phi(y))+(w(x)-w(y))
 (\psi(x)- \psi(y))}{|x-y|^{n+2s}} \,dx\,dy \\
&\quad - \int_{\Omega}(\lambda f(x) u_{+}^{-q}\phi+ \mu g(x) w_{+}^{-q}\psi) dx 
  - \frac{\alpha}{\alpha+\beta}\int_{\Omega} b(x) u_{+}^{\alpha-1}w_{+}^{\beta}\phi dx \\
&\quad -\frac{\beta}{\alpha+\beta}\int_{\Omega} b(x) u_{+}^{\alpha}w_{+}^{\beta-1}\psi\Big)\\
&\quad +2\epsilon\Big(\int_{\Gamma_1\times \Omega_1}\frac{|u(x)-u(y)|^{2}}{|x-y|^{n+2s}}\,dx\,dy
 \Big)^{1/2}
 \Big(\int_{\Gamma_1\times \Omega_1}\frac{|\phi(x)|^{2}}{|x-y|^{n+2s}}\,dx\,dy\Big)^{1/2}\\
&\quad+2\epsilon\Big(\int_{\Gamma_2\times \Omega_2}\frac{|w(x)-w(y)|^{2}}{|x-y|^{n+2s}}\,dx\,dy
 \Big)^{1/2}
 \Big(\int_{\Gamma_2\times \Omega_2}\frac{|\psi(x)|^{2}}{|x-y|^{n+2s}}\,dx\,dy\Big)^{1/2}\\
&\quad +2\epsilon\Big[\Big(\int_{\Gamma_1\times \Omega^c}\frac{|u(x)|^{2}}{|x-y|^{n+2s}}
 \,dx\,dy\Big)^{1/2}
 \Big(\int_{\Gamma_1\times \Omega^c}\frac{|\phi(x)|^{2}}{|x-y|^{n+2s}}\,dx\,dy\Big)^{1/2}\\
&\quad +\Big(\int_{\Gamma_2\times \Omega^c}\frac{|w(x)|^{2}}{|x-y|^{n+2s}}\,dx\,dy
 \Big)^{1/2}
 \Big(\int_{\Gamma_2\times \Omega^c}\frac{|\psi(x)|^{2}}{|x-y|^{n+2s}}\,dx\,dy\Big)^{1/2}\\
&\quad + \Big(\int_{\Gamma_1\times \Omega_1}\frac{|u(x)-u(y)|^{2}}{|x-y|^{n+2s}}\,dx\,dy
 \Big)^{1/2}
 \Big(\int_{\Gamma_1\times \Omega_1}\frac{|\phi(x)|^{2}}{|x-y|^{n+2s}}\,dx\,dy\Big)^{1/2}\\
&\quad + \Big(\int_{\Gamma_2\times \Omega_2}\frac{|w(x)-w(y)|^{2}}{|x-y|^{n+2s}}\,dx\,dy
 \Big)^{1/2}
 \Big(\int_{\Gamma_2\times \Omega_2}\frac{|\psi(x)|^{2}}{|x-y|^{n+2s}}\,dx\,dy\Big)^{1/2}\\
&\quad+\left({\int_{\Gamma_1\times \Gamma_1}}\frac{|u(x)-u(y)|^{2}}{|x-y|^{n+2s}}\,dx\,dy\right)^{1/2}\left(\int_{\Gamma_1\times \Gamma_1}\frac{|\phi(x)-\phi(y)|^{2}}{|x-y|^{n+2s}}\,dx\,dy\right)^{1/2}\\
&\quad+\Big({\int_{\Gamma_2\times \Gamma_2}}\frac{|w(x)-w(y)|^{2}}{|x-y|^{n+2s}}\,dx\,dy
 \Big)^{1/2}
 \Big(\int_{\Gamma_2\times \Gamma_2}\frac{|\psi(x)-\psi(y)|^{2}}{|x-y|^{n+2s}}\,dx\,dy
 \Big)^{1/2}\Big]\\
&\quad + \epsilon \epsilon^{\alpha} \|b\|_{\infty} 
 \Big(\int_{\Gamma_1}|\psi|^{\alpha+\beta}dx\Big)^{\frac{\beta}{\alpha+\beta}}
 \Big(\int_{\Gamma_1}(w_+)^{\alpha+\beta}dx\Big)^{\frac{\beta}{\alpha+\beta}}\\
&\quad +\epsilon \epsilon^{\beta} \|b\|_{\infty}
 \Big(\int_{\Gamma_2}(u_+)^{\alpha+\beta}dx\Big)^{\frac{\alpha}{\alpha+\beta}}
 \Big(\int_{\Gamma_2}|\phi|^{\alpha+\beta}dx\Big)^{\frac{\beta}{\alpha+\beta}}\\
&\quad +\frac{\epsilon\alpha}{\alpha+\beta} 
 \int_{\Gamma_1} b(x) u_{+}^{\alpha-1}w_{+}^{\beta}\phi dx
 +\frac{\epsilon\beta}{\alpha+\beta} \int_{\Gamma_2} b(x)u_{+}^{\alpha} 
 w_{+}^{\beta-1}\psi dx.
\end{align*}
Since the measure of $\Gamma_i =\{x\in \Omega| (u_i +\epsilon \phi_i)(x)\leq 0\}$ tend to zero 
as $\epsilon\to 0$, it follows that
\[
\int_{\Gamma_i\times \Omega_i}\frac{|\phi_i(x)|^{2}}{|x-y|^{n+2s}}\,dx\,dy\to 0,\quad
 \text{as } \epsilon\to 0,
\] 
and similarly
\begin{gather*}
\int_{\Gamma_i\times \Omega^c}\frac{|\phi_i(x)|^{2}}{|x-y|^{n+2s}}\,dx\,dy, \quad
\int_{\Gamma_i\times \Gamma_i}\frac{|\phi_i(x)-\phi_i(y)|^{2}}{|x-y|^{n+2s}}\,dx\,dy,\\
\int_{\Gamma_1} b(x) u_{+}^{\alpha-1}w_{+}^{\beta}\phi dx, \quad
\int_{\Gamma_2} b(x)u_{+}^{\alpha} w_{+}^{\beta-1}\psi dx,
\end{gather*}  
all are tend to $0$ as $\epsilon\to 0$. Dividing by $\epsilon$ and letting $\epsilon\to 0$, we obtain
 \begin{align*}
&\int_{Q} \frac{(u(x)-u(y)) (\phi(x)-\phi(y))+(w(x)-w(y)) 
 (\psi(x)-\psi(y))}{|x-y|^{n+2s}}\,dx\,dy\\
&- \int_{\Omega} (\lambda f(x) u_{+}^{-q}\phi+\mu g(x) w_{+}^{-q}\psi) dx 
 - \frac{\alpha}{\alpha+\beta} \int_{\Omega} b(x) u_{+}^{\alpha-1}w_{+}^{\beta}\phi dx \\
&- \frac{\beta}{\alpha+\beta} \int_{\Omega} b(x) u_{+}^{\alpha}w_{+}^{\beta-1}\psi dx \geq 0.
\end{align*}
Since this holds equally well for $(-\phi,-\psi)$, it follows that $(u,w)$ 
is indeed a positive weak solution of problem \eqref{ePplm}
and hence a positive solution of \eqref{ePlm}.
\end{proof}

\begin{lemma}\label{lm1}
 There exists a minimizing sequence $\{(U_k,W_k)\}$ in $\mathcal{N}_{\lambda,\mu}^{-}$
such that $(U_k,W_k)\to (U,W)$ strongly in $\mathcal{N}_{\lambda,\mu}^{-}$.
Moreover $(U,W)$ is a positive weak solution of $\eqref{ePlm}$.
\end{lemma}

\begin{proof}
Using the Ekeland variational principle again, we may find a minimizing sequence 
$\{(U_k,W_k)\}\subset \mathcal{N}_{\lambda,\mu}^{-}$ for the minimizing problem
 $\inf_{\mathcal{N}_{\lambda,\mu}^{-}}J_{\lambda,\mu}$ such that for
$(U_k,W_k)\rightharpoonup (U,W)$ weakly in $ Y$ and pointwise a.e. in $Q$. 
We can repeat the argument used in Lemma \ref{le6} to derive that when 
$(\lambda,\mu)\in \Gamma$
\begin{equation} \label{s11}
\begin{aligned}
&(1+q)\int_{\Omega} \left(\lambda f(x) U_{+}^{1-q}+ \mu g(x)W_{+}^{1-q}\right)dx \\
&- (\alpha+\beta-2)\int_{\Omega} b(x)U_{+}^{\alpha} W_{+}^{\beta} dx<0
\end{aligned}
\end{equation}
which yields
\begin{align*}
&(1+q)\int_{\Omega} \left(\lambda f(x) (U_k)_{+}^{1-q}+ \mu g(x)(W_k)_{+}^{1-q}\right)dx \\
&-(\alpha+\beta-2)\int_{\Omega} b(x) (U_k)_{+}^{\alpha}(W_k)_{+}^{\beta}dx
\leq -C_4
\end{align*}
for $k$ sufficiently large and a suitable positive constant $C_4$.
At this point we may proceed exactly as in
Lemmas \ref{le7}, \ref{le8}, \ref{le9}, \ref{le10} and corollary \ref{cc1}
to conclude that $U$, $W>0$ and $(U,W)$ is the required positive weak solution
of problem \eqref{ePplm}.
In particular, $(U,W)\in \mathcal{N}_{\lambda,\mu}$. Moreover from \eqref{s11}
it follows that
\begin{align*}
&(1+q)\|(U,W)\|^2 -(\alpha+\beta-1+q)\int_{\Omega} b(x)U_{+}^{\alpha} W_{+}^{\alpha} dx\\
&=(1+q)\Big[ K_{\lambda,\mu}(U,W) + \int_{\Omega} b(x)U_{+}^{\alpha} W_{+}^{\beta} dx\Big]
- (\alpha+\beta-1+q)\int_{\Omega} b(x)U_{+}^{\alpha} W_{+}^{\beta} dx\\
&=(1+q)K_{\lambda,\mu}(U,W) -(\alpha+\beta-2)\int_{\Omega} b(x) U_{+}^{\alpha}W_{+}^{\beta} dx<0,
\end{align*}
that is $(U,W)\in\mathcal{N}_{\lambda,\mu}^{-}$.
\end{proof}

\begin{proof}[Proof of the Theorem \ref{th1}]
 From Lemmas \ref{le10}, \ref{lm1} and \ref{le2}, we can conclude that the 
problem \eqref{ePlm} has at least two positive weak solutions 
$(u,w)\in\mathcal{N}_{\lambda,\mu}^{+}$, $(U,W)\in\mathcal{N}_{\lambda,\mu}^{-}$
with $\|(U,W)\|>\|(u,w)\|$ for any $(\lambda,\mu)\in\Gamma$.
\end{proof}

\begin{thebibliography}{00}

\bibitem{AJ} Adimurthi, J. Giacomoni; 
\emph{Multiplicity of positive solutions for a singular and critical elliptic 
problem in $\mathbb R^2$}, Comm. in Contemporary Mathematics, 8 (5) (2006), 621-656.

\bibitem{Am} Ahmed Mohammed; 
\emph{Positive solutions of the $p$-Laplace equation with singular nonlinearity},
 J. Math.  Anal. Appl., 352 (2009), 234-245.

\bibitem{da} D. Applebaum; 
\emph{L\`{e}vy process-from probability to finance and quantum groups}, 
Notices Amer. Math. Soc., 51 (2004), 1336-1347.

\bibitem{bb} B. Barrios, I. D. Bonis, M. Medina,  I. Peral; 
\emph{Semilinear problems for the fractional Laplacian with a singular nonlinearity},
 Open Math.,  13 (2015), 390-407.

\bibitem{cs} L. Caffarelli, L. Silvestre; 
\emph{An extension problem related to the fractional Laplacian}, 
Comm. in Partial Differential Equations, 32 (2007), 1245-1260.

\bibitem{yc} Y. Chen, J. Chen; 
\emph{Existence of multiple positive weak solutions and estimates for extremal 
values to a class of elliptic problems with Hardy term and singular nonlinearity}, 
J. Math. Anal. Appl., 429 (2015), 873-900.

\bibitem{cr} M. G. Crandall, P. H. Rabinowitz, L. Tartar;
 \emph{On a Dirichlet problem with a singular nonlinearity},
 Comm. Partial Differential Equations, 2 (1977), 193-222.

\bibitem{cp} M. M. Coclite, G. Palmieri; 
\emph{On a singular nonlinear Dirichlet problem}, 
Comm. Partial Differential Equations,  14 (10) (1989), 1315-1327.

\bibitem{dh} J. I. D\'{i}az, J. Hern\'{a}ndez, F. Mancebo; 
\emph{Branches of positive and free boundary solutions for some singular 
quasilinear elliptic problems}, J. Math. Anal. Appl., 352 (2009), 449-474.

 \bibitem{dhr} J. I. D\'{i}az, J. Hern\'{a}ndez, J. M. Rakotoson; 
\emph{On very weak positive solutions to some semilinear elliptic problems 
with simultaneous singular nonlinear and spatial dependence terms}, 
Milan J. Math., 79 (2011), 233-245.

\bibitem{DP} P. Drabek, S. I. Pohozaev; 
\emph{Positive solutions for the $p$-Laplacian: application of the fibering method},
 Proceedings of Royal Society of Edinburgh Section A, 127 (1997), 703-726.

\bibitem{gr} M. Ghergu, V. R\v{a}dulescu; 
\emph{Sublinear singular elliptic problems with two parameters}, 
J. Differential Equations, 195 (2003), 520-536.

\bibitem{gr1} M. Ghergu, V. R\v{a}dulescu; 
\emph{Bifurcation for a class of singular elliptic problems with quadratic 
convection term}, C. R. Acad. Sci. Paris Ser. I, 338 (2004), 831-836.

\bibitem{fa} Y. Fang; 
\emph{Existence, uniqueness of positive solution to a fractional laplacians 
with singular non linearity}, arXiv:1403.3149v1,  2014.

\bibitem{gm} A. Garroni, S. M\"{u}ller; 
\emph{G-limit of a phase-field model of dislocations}, 
SIMA J. Math. Anal., 36 (2005), 1943-1964.

\bibitem{J1} J. Giacomoni, H. M\^{a}agli, P. Sauvy; 
\emph{Existence of compact support solutions for a quasilinear and singular problem},
     Differential Integral Equations, 25 (2012), 629-656.

\bibitem{j2} J. Giacomoni, I. Schindler, P. Tak\'{a}\u{c}; 
\emph{Sobolev versus H\"{o}lder local minimizers and existence of multiple solutions 
for  a singular quasilinear equation}, Ann. Sc. Norm. Super. Pisa Cl. Sci., 5 (6)
 (2007), 117-158.

\bibitem{jh} J. Giacomoni, J. Hern\'{a}ndez, A. Mouassaoui; 
\emph{Quasilinear and singular systems: the cooperative case}, 
Contemp. Math., 540  (2011), 79-94.

\bibitem{cv} Sarika Goyal; 
\emph{Multiplicity results of fractional $p$-Laplace equations with 
sign-changing and singular nonlinearity}, Complex variables and Elliptic equations,
 62 (2016), 158-183.

\bibitem{s1} Sarika Goyal, K. Sreenadh;
\emph{Fractional elliptic equations with sign-changing and singular nonlinearity}, 
 Electronic Journal of Differential Equations, 2016 (2016) No.1 45, 1-23.

\bibitem{ssa} S. Goyal, K. Sreenadh; 
\emph{Existence of multiple solutions of $p$-fractional Laplace operator with 
sign-changing weight function},
 Adv. Nonlinear Anal., 4 (1) (2015), 37-58.

\bibitem{ssi} S. Goyal, K. Sreenadh; 
\emph{The Nehari manifold for non-local elliptic operator with concave-convex 
nonlinearities and sign-changing weight functions}, 
Proc. Indian Acad. Sci. Math. Sci., 125 (4) (2015), 545-558.	

\bibitem{hss} N. Hirano, C. Saccon and N. Shioji; 
\emph{Existence of multiple positive solutions for singular elliptic problems 
with concave and convex nonlinearities}, Adv. Differential Equations, 9 (2) (2004),
 197-220.

\bibitem{hnc} N. Hirano, C. Saccon, N. Shioji; 
\emph{Brezis-Nirenberg type theorems and multiplicity of positive solutions
 for a singular elliptic problem}, {J. Differential Equations}, 245 (2008), 1997-2037.

\bibitem{lm} A. C. Lazer, P. J. McKenna; 
\emph{On a Singular nonlinear elliptic boundary-value problem}, 
Proc. of Amer. Math. Soc., 111 (3) (1991), 721-730.

\bibitem{lmp} A. C. Lazer, P. J. McKenna; 
\emph{On Singular boundary value problems for the Monge-Ampère operator}, 
J. Math. Anal. Appl., 197 (1996), 341-362.

\bibitem{hi} E. Di Nezza, G. Palatucci, E. Valdinoci; 
\emph{Hitchhiker's guide to the fractional Sobolev spaces},
 Bull. Sci. Math., 136 (2012), 225-236.

\bibitem{mp} R. Servadei, E. Valdinoci; 
\emph{Mountain pass solutions for non-local elliptic operators,} 
J. Math. Anal. Appl., 389 (2012), 887-898.

\bibitem{tu} Tuhina Mukherjee, K. Sreenadh; 
\emph{Critical growth fractional elliptic problem with singular nonlinearities},
 arXiv:1602.07886 [math.AP].

\bibitem{xy} Xiaohui Yu; 
\emph{The Nehari manifold for elliptic equation involving the square root 
of the Laplacian}, J. Differential Equations, 252 (2012), 1283-1308.

\end{thebibliography}

\end{document}