# Existence and uniqueness results to a quasilinear singular Lane-Emden
Neumann system

Nouredine Medjoudj  Nouredine Medjoudj
Applied Mathematics Laboratory, Faculty of Exact Sciences,
A. Mira Bejaia University, Targa Ouzemour, 06000 Bejaia, Algeria
<EMAIL_ADDRESS>and  Abdelkrim Moussaoui  Abdelkrim
Moussaoui
Applied Mathematics Laboratory (LMA), Faculty of Exact Sciences
and Department of Physico-Chemical Biology (BPC), Faculty of Natural and Life
Sciences,
A. Mira Bejaia University, Algeria<EMAIL_ADDRESS>
###### Abstract.

We establish the existence and uniqueness of solutions for quasilinear
singular Lane-Emden type systems subjected to Neumann boundary conditions. The
approach is chiefly based on sub-supersolutions method.

###### Key words and phrases:

sub-supersolutions, Neumann boundary conditions.

###### Key words and phrases:

Singularity, Lane-Emden system, Neumann boundary conditions, sub-
supersolutions, uniqueness

###### 2010 Mathematics Subject Classification:

35J75, 35J62, 35J92

##  1\. Introduction

Let $\Omega$ be a bounded domain in $\mathbb{R}^{N}$ ($N\geq 2)$ having a
smooth boundary $\partial\Omega.$ Given $1<p_{i}<N$ for $i=1,2,$ we consider
the following quasilinear Lane-Emden type system

$\left(\mathrm{P}\right)\qquad\left\\{\begin{array}[]{ll}-\Delta_{p_{1}}u+|u|^{p_{1}-2}u=u^{\alpha_{1}}+v^{\beta_{1}}&\text{in}\;\Omega,\\\
-\Delta_{p_{2}}v+|v|^{p_{2}-2}v=u^{\alpha_{2}}+v^{\beta_{2}}&\text{in}\;\Omega,\\\
u,v>0&\text{in }\;\Omega,\\\ \frac{\partial u}{\partial\eta}=\frac{\partial
v}{\partial\eta}=0&\text{on}\;\partial\Omega,\end{array}\right.$

where $\eta$ is the unit outer normal to $\partial\Omega$, while
$\Delta_{p_{i}}$ denotes the $p_{i}$-Laplace operator, namely
$\Delta_{p_{i}}:=\mathrm{div}(|\nabla w|^{p_{i}-2}\nabla w)$, $\forall\,w\in
W^{1,p_{i}}(\Omega).$ We consider system $\left(\mathrm{P}\right)$ in a
singular case assuming that the exponents verify the condition

(1.1) $-1<\alpha_{1},\beta_{2}<0,\text{ \ }-1<\beta_{1}<p_{1}-1\text{ \ and \
}-1<\alpha_{2}<p_{2}-1.$

A solution $(u,v)\in W^{1,p_{1}}(\Omega)\times W^{1,p_{2}}(\Omega)$ of
$\left(\mathrm{P}\right)$ is understood in the weak sense, that is,

(1.2) $\left\\{\begin{array}[]{l}\int_{\Omega}|\nabla u|^{p_{1}-2}\nabla
u\nabla\varphi\,\mathrm{d}x+\int_{\Omega}|u|^{p_{1}-2}u\varphi\,\mathrm{d}x=\int_{\Omega}(u^{\alpha_{1}}+v^{\beta_{1}})\varphi\,\mathrm{d}x\\\
\int_{\Omega}|\nabla v|^{p_{2}-2}\nabla
v\nabla\psi\,\mathrm{d}x+\int_{\Omega}|v|^{p_{2}-2}v\,\psi\,\mathrm{d}x=\int_{\Omega}(u^{\alpha_{2}}+v^{\beta_{2}})\,\psi\,\mathrm{d}x\end{array}\right.$

for all $(\varphi,\psi)\in W^{1,p_{1}}(\Omega)\times W^{1,p_{2}}(\Omega)$.

System $\left(\mathrm{P}\right)$ is a natural extension and generalization of
the celebrated Lane-Emden equation

(1.3) $\Delta w+w^{\gamma}=0\text{ in }\Omega,$

subjected to Neumann boundary conditions. It is introduced by Homer Lane [21],
who was interested in computing both the temperature and the density of mass
on the surface of the sun. The Lane-Emden equation (1.3) has been the focus of
a huge number of works which we shall not discuss merely mentioning some of
them [2, 3, 12]. It is involved in wide range of phenomena in mathematical
physics and chemistry, specifically in the areas of conformal geometry,
thermal explosion, isothermal gas spheres and thermionic currents [31]. In
astrophysics, it describes the behavior of the density of a gas sphere in
hydrostatic equilibrium. The index $\gamma,$ called the polytropic index, is
related to the ratio $r$ of the specific heats of the gas throught
$\gamma=\frac{1}{r-1}.$ It is used to determine the structure of the interior
of polytropic stars, which are subject to the influence of their own
gravitational field [5]. The case $\gamma<0$ in (1.3) is highly challenging as
it brings out a singularity at the origin. It has attracted considerable
interest in recent years. We mention [8] where it is shown that problem (1.3)
subject to Dirichlet boundary conditions admits a unique positive solution $u$
in $\mathcal{C}^{2}(\Omega)\cap\mathcal{C}(\overline{\Omega})$. In addition,
this solution $u$ belongs to
$\mathcal{C}^{2}(\Omega)\cap\mathcal{C}^{1}(\overline{\Omega})$ once
$-1<\gamma<0.$ When the quasilinear $p$-Laplacian operator is involved in
singular problem (1.3), [15] provides a unique solution in
$\mathcal{C}^{1,\tau}(\overline{\Omega})$, $\tau\in(0,1),$ for $-1<\gamma<0$.

Singularities are an important feature in the study of the problem
$\left(\mathrm{P}\right)$. They occur near the origin under assumption (1.1)
on nonlinearities. This fact represents a serious difficulty to overcome
especially since a very marked singularity character is considered for
$\left(\mathrm{P}\right)$, resulting from (1.1) when all exponents are
negative. Actually, singularities are involved in a wide range of important
elliptic problems which have been studied extensively in recent years, see for
instance [1, 8, 10, 11, 13, 14, 15, 17, 19, 24, 25, 26, 28, 27] and there
references. The semilinear case, that is when $p_{1}=p_{2}=2$, has been widely
investigated, especially in the context of the Gierer Meinhardt system (see,
e.g., [27] and the references therein). However, as far as we know, [16] is
the only paper where Neumann boundary conditions is considered for singular
quasilinear systems. We emphasize that system $\left(\mathrm{P}\right)$ cannot
be incorporated neither in Gierer-Meinhardt system addressed in [27] nor in
the convective system studied in [16] even if gradient terms involved are
canceled.

In the aforementioned papers, two complementary structures are separately
discussed. In [1, 10, 11, 14, 19, 26], the systems examined are cooperative
while in [13, 28, 24, 25], they are competitive. For system
$\left(\mathrm{P}\right)$, these both structures are closely related to the
sign of exponents $\alpha_{2}$ and $\beta_{1}$. Namely,
$\left(\mathrm{P}\right)$ is cooperative if $\min\\{\alpha_{2},\beta_{1}\\}>0$
whereas for opposite strict inequality, $\left(\mathrm{P}\right)$ is
competitive. It should be noted that, in the present paper, these both
complementary structures for the system $\left(\mathrm{P}\right)$ are handled
simultaneously without referring to them despite their important structural
disparity that makes the right hand side in $\left(\mathrm{P}\right)$ behaving
in a drastically different way.

Motivated by the aforementioned facts, our main concern is the question of
existence of solutions for system $\left(\mathrm{P}\right)$. By adequate
truncation and owing to Schauder’s fixed point theorem, we first develop two
new sub-supersolutions Theorems for a general singular Neumann Systems
involving $p$-Laplacian operator (cf. Theorems 2.1 and 2.2, section 2). These
results address respectively the case of bounded and unbounded nonlinearities,
restricted to the rectangle formed by the sub-supersolution pair. Thence, the
nature and properties of the sub-supersolutions are meaningfully impacted
which, henceforth, can be constructed on a wider choice of functions than
specified in the literature. Namely, this enable to consider subsolutions
subject to homogeneous Dirichlet boundary conditions, what was not possible in
[16] under assumption $($H$)$ stated therein. It is worth noting that, unlike
to what has been stated in [27], sub-supersolutions cannot have both zero
traces on the boundary $\partial\Omega$. This would lead to a solution for
problem $\left(\mathrm{P}\right)$ with zero trace condition on
$\partial\Omega$ and therefore, according to [17, Lemma 3.1], its normal
derivative would be nontrivial, which is absurd. Another important issue being
addressed by Theorems 2.1 and 2.2 concern the property sign of normal
derivative on $\partial\Omega$ of sub-supersolutions. It is established that
the condition (2.3) is crucial to handle quasilinear Neumann problems via sub-
supersolutions method and therefore, in no case could be ignored. Furthermore,
in the aforementioned Theorems, we mention that no sign condition is required
on the right-hand side nonlinearities. Hence, they can be used for large
classes of quasilinear singular problems, including those discussed in [16,
27].

In the prospect of applying Theorems 2.1 and 2.2, pairs of sub-supersolutions
for system $\left(\mathrm{P}\right)$ are constructed by exploiting spectral
properties of the $p$-Laplacian operator as well as properties of torsion
problems, both subject to Dirichlet or Neumann boundary conditions (cf.
section 4). A suitable adjustment of constants is also required. Moreover, by
making some specific and necessary adjustments, it is quite possible to
construct a pair of sub-supersolution in the spirit of [16] which verifies
assumption (2.3). Nevertheless, although their new form is appropriate, it no
longer leads to an infinite sequence of solutions, as it is stated in [16].
This seems to be unfeasible under assumption (2.3).

The sub-supersolutions pairs establish a location of solutions of
$\left(\mathrm{P}\right)$ provided by Theorems 3.1 and 3.2 (cf. section 3). We
show that their regularity property depend on the behavior of the subsolution
near the boundary $\partial\Omega$. Precisely, when the subsolution has a zero
trace condition on $\partial\Omega,$ solutions $(u,v)$ of
$\left(\mathrm{P}\right)$ are bounded in $L^{\infty}(\Omega)\times
L^{\infty}(\Omega)$, whereas in the complementary case, $(u,v)$ are bounded in
$\mathcal{C}^{1}(\overline{\Omega})\times\mathcal{C}^{1}(\overline{\Omega})$.
Note that the absence of $\mathcal{C}^{1}$-bound for solutions in the previous
case is rather due to the fact that we have not been able to find in the
literature an equivalent of the regularity theorem for singular Dirichlet
problems [17, Lemma 3.1] to the case of Neumann boundary conditions.
Regularity results for singular Neumann problems remains an open question.

A second main objective of our work is to provide a criterion ensuring the
uniqueness of a solution to problem $\left(\mathrm{P}\right)$. In this
direction, the uniqueness result is established for $\mathcal{C}^{1}$-bound
solutions. The proof is based on the adaption of argument by Krasnoselskii
[18] where properties of sub-supersolutions are crucial. In this part of our
work, we assume all exponents in (1.1) negative, subject to some additional
restrictions as described in our result (cf. section 3). It is worth noting
that the previous argument does not apply to $L^{\infty}$\- bound solutions
simply because the obtained supersolutions are not comparable to the distance
function $d(x)$ in $\Omega$.

The rest of this article is organized as follows. Section 2 contains the
existence theorems involving sub-supersolutions. Section 3 presents abstract
existence and uniqueness results. Section 4 deals with the existence and
uniqueness of solutions for problem $\left(\mathrm{P}\right)$.

##  2\. Sub-super-solution theorems

The Sobolev spaces $W^{1,p_{i}}(\Omega)$ and $W_{0}^{1,p_{i}}(\Omega)$ will be
equipped with the norm

$\|w\|_{1,p_{i}}:=\left(\|w\|_{p_{i}}^{p_{i}}+\|\nabla
w\|_{p_{i}}^{p_{i}}\right)^{\frac{1}{p_{i}}},\quad w\in W^{1,p_{i}}(\Omega),$
$\|w\|_{1,p_{i}}:=\|\nabla w\|_{p_{i}},\quad w\in W_{0}^{1,p_{i}}(\Omega),$

where, as usual,

$\|w\|_{p_{i}}:=\left\\{\begin{array}[]{ll}\left(\int_{\Omega}|w(x)|^{p_{i}}\mathrm{d}x\right)^{\frac{1}{p_{i}}}&\text{
if }p_{i}<+\infty,\\\ &\\\ \underset{x\in\Omega}{ess\sup}\,|w(x)|&\text{
otherwise.}\end{array}\right.$

We denote by $(W^{1,p_{i}}(\Omega))^{\ast}$ the topological dual space of
$W^{1,p_{i}}(\Omega)$. Moreover,

$W_{+}^{1,p_{i}}(\Omega):=\\{w\in W^{1,p_{i}}(\Omega):w\geq 0\\},\quad
W_{b}^{1,p_{i}}(\Omega):=W^{1,p_{i}}(\Omega)\cap L^{\infty}(\Omega),$

and

$W_{0,b}^{1,p_{i}}(\Omega):=W_{0}^{1,p_{i}}(\Omega)\cap L^{\infty}(\Omega)$

We also utilize

$\displaystyle\mathcal{C}_{+}^{1,\tau}(\overline{\Omega})$ $\displaystyle:$
$\displaystyle=\\{w\in\mathcal{C}^{1,\tau}(\overline{\Omega}):w\geq 0\text{
for all }x\in\overline{\Omega}\\},$ $\displaystyle
int\mathcal{C}_{+}^{1,\tau}(\overline{\Omega})$ $\displaystyle:$
$\displaystyle=\\{w\in\mathcal{C}_{+}^{1,\tau}(\overline{\Omega}):w>0\text{
for all }x\in\overline{\Omega}\\},\text{ for }\tau\in(0,1).$

In what follows, we set $r^{\pm}:=\max\\{\pm r,0\\}$ and we denote by
$\gamma_{0}$ the unique continuous linear map
$\gamma_{0}:W^{1,p_{i}}(\Omega)\rightarrow L^{p_{i}}(\partial\Omega)$ known as
the trace map such that $\gamma_{0}(u)=u_{/\partial\Omega},$ for all $u\in
W^{1,p}(\Omega)$ and verifies the property

(2.1) $\gamma_{0}(u^{+})=\gamma_{0}(u)^{+}\text{ \ for all }u\in
W^{1,p_{i}}(\Omega),$

(see, e.g., [23]).

Hereafter, we denote by $d(x)$ the distance from a point
$x\in\overline{\Omega}$ to the boundary $\partial\Omega$, where
$\overline{\Omega}=\Omega\cup\partial\Omega$ is the closure of
$\Omega\subset\mathbb{R}^{N}$. For $1<r<N$ and $-r<s\leq 0,$ it is known that

$\left(\int_{\Omega}d(x)^{s}|u(x)|^{r}\mathrm{d}x\right)^{\frac{1}{r}}\leq
C\|u\|_{1,r}\quad\forall\,u\in W^{1,r}(\Omega),$

with suitable $C>0$; see [29, Theorem 19.9, case (19.29)]. Accordingly, by
Hölder’s inequality, if $-1<\beta\leq 0$ then

(2.2)
$\int_{\Omega}|d(x)^{\beta}u(x)|\,\mathrm{d}x\leq|\Omega|^{\frac{1}{r^{\prime}}}\left(\int_{\Omega}d(x)^{\beta
r}|u(x)|^{r}\mathrm{d}x\right)^{\frac{1}{r}}\leq
C|\Omega|^{\frac{1}{r^{\prime}}}\|u\|_{1,r},\;\;u\in W^{1,r}(\Omega).$

Finally, we say that $j:\Omega\times\mathbb{R}^{2}\rightarrow\mathbb{R}$ is a
Carathéodory function provided

  * $\bullet$

$x\mapsto j(x,s,t)$ is measurable for every $(s,t)\in\mathbb{R}^{2}$, and

  * $\bullet$

$(s,t)\mapsto j(x,s,t)$ is continuous for almost all $x\in\Omega$.

This section investigates the existence of solutions to system

$\left(\mathrm{P}_{f_{1},f_{2}}\right)\qquad\left\\{\begin{array}[]{ll}-\Delta_{p_{1}}u+|u|^{p_{1}-2}u=f_{1}(x,u,v)&\text{in}\;\Omega,\\\
-\Delta_{p_{2}}v+|v|^{p_{2}-2}v=f_{2}(x,u,v)&\text{in}\;\Omega,\\\
\frac{\partial u}{\partial\eta}=\frac{\partial
v}{\partial\eta}=0&\text{on}\;\partial\Omega,\end{array}\right.$

where $f_{i}:\Omega\times\mathbb{R}^{2}\rightarrow\mathbb{R}$ satisfy
Carathéodory’s conditions. The following assumptions will be posited.

  * $(\mathrm{H}_{1})$

With appropriate $(\underline{u},\underline{v}),(\overline{u},\overline{v})\in
W_{b}^{1,p_{1}}(\Omega)\times W_{b}^{1,p_{2}}(\Omega)$ such that

(2.3)
$\max\\{\frac{\partial\underline{u}}{\partial\eta},\frac{\partial\underline{v}}{\partial\eta}\\}\leq
0\leq\min\\{\frac{\partial\overline{u}}{\partial\eta},\frac{\partial\overline{v}}{\partial\eta}\\}\text{
on }\partial\Omega,$

one has $\underline{u}\leq\overline{u}$, $\underline{v}\leq\overline{v}$, as
well as

(2.4)
$\left\\{\begin{array}[]{l}\int_{\Omega}|\nabla\underline{u}|^{p_{1}-2}\nabla\underline{u}\nabla\varphi\,\mathrm{d}x+\int_{\Omega}|\underline{u}|^{p_{1}-2}\underline{u}\varphi\,\mathrm{d}x\\\
-\int_{\partial\Omega}|\nabla\underline{u}|^{p_{1}-2}\frac{\partial\underline{u}}{\partial\eta}\gamma_{0}(\varphi)\text{
}\mathrm{d}s\leq\int_{\Omega}f_{1}(\cdot,\underline{u},v)\varphi\,\mathrm{d}x,\\\
\\\
\int_{\Omega}|\nabla\underline{v}|^{p_{2}-2}\nabla\underline{v}\,\nabla\psi\,\mathrm{d}x+\int_{\Omega}|\underline{v}|^{p_{2}-2}\underline{v}\,\psi\,\mathrm{d}x\\\
-\int_{\partial\Omega}|\nabla\underline{v}|^{p_{2}-2}\frac{\partial\underline{v}}{\partial\eta}\gamma_{0}(\psi)\text{
}\mathrm{d}s\leq\int_{\Omega}f_{2}(\cdot,u,\underline{v})\psi\,\mathrm{d}x,\end{array}\right.$
(2.5)
$\left\\{\begin{array}[]{l}\int_{\Omega}|\nabla\overline{u}|^{p_{1}-2}\nabla\overline{u}\,\nabla\varphi\,\mathrm{d}x+\int_{\Omega}|\overline{u}|^{p_{1}-2}\overline{u}\,\varphi\,\mathrm{d}x\\\
-\int_{\partial\Omega}|\nabla\overline{u}|^{p_{1}-2}\frac{\partial\overline{u}}{\partial\eta}\gamma_{0}(\varphi)\text{
}\mathrm{d}s\geq\int_{\Omega}f_{1}(\cdot,\overline{u},v)\varphi\,\mathrm{d}x,\\\
\\\
\int_{\Omega}|\nabla\overline{v}|^{p_{2}-2}\nabla\overline{v}\,\nabla\psi\,\mathrm{d}x+\int_{\Omega}|\overline{v}|^{p_{2}-2}\overline{v}\,\psi\,\mathrm{d}x\\\
-\int_{\partial\Omega}|\nabla\overline{v}|^{p_{2}-2}\frac{\partial\overline{v}}{\partial\eta}\gamma_{0}(\psi)\text{
}\mathrm{d}s\geq\int_{\Omega}f_{2}(\cdot,u,\overline{v})\psi\,\mathrm{d}x\end{array}\right.$

for all $(\varphi,\psi)\in W_{+}^{1,p_{1}}(\Omega)\times
W_{+}^{1,p_{2}}(\Omega)$, $(u,v)\in W^{1,p_{1}}(\Omega)\times
W^{1,p_{2}}(\Omega)$ such that
$(u,v)\in[\underline{u},\overline{u}]\times[\underline{v},\overline{v}]$.

  * $(\mathrm{H}_{2})$

For appropriate $C>0$ and $\gamma\in(-1,0)$ one has

$\max|f_{i}(x,s,t)|\leq
Cd(x)^{\gamma}\quad\text{in}\quad\Omega\times[\underline{u},\overline{u}]\times[\underline{v},\overline{v}].$

Note that under $(\mathrm{H}_{2})$ and by virtue of the Hardy-Sobolev
inequality type in (2.2), the integrals involving $f_{1}$ and $f_{2}$ in (2.4)
and (2.5) take sense .

######  Theorem 2.1.

Suppose $(\mathrm{H}_{1})$–$(\mathrm{H}_{2})$ hold true. Then, problem
$(\mathrm{P}_{f_{1},f_{2}})$ possesses a solution $(u,v)\in
W_{b}^{1,p_{1}}(\Omega)\times W_{b}^{1,p_{2}}(\Omega)$ such that

(2.6) $\underline{u}\leq
u\leq\overline{u}\quad\text{and}\quad\underline{v}\leq v\leq\overline{v}.$

Moreover, $\frac{\partial u}{\partial\eta}=\frac{\partial v}{\partial\eta}=0$
on $\partial\Omega$.

###### Proof.

Given $(z_{1},z_{2})\in L^{p_{1}}(\Omega)\times L^{p_{2}}(\Omega)$, we define

(2.7)
$\mathrm{T}_{1}(z_{1}):=\left\\{\begin{array}[]{ll}\underline{u}&\text{when
}z_{1}\leq\underline{u},\\\ z_{1}&\text{if }\underline{u}\leq
z_{1}\leq\overline{u},\\\
\overline{u}&\text{otherwise,}\end{array}\right.\quad\mathrm{T}_{2}(z_{2}):=\left\\{\begin{array}[]{ll}\underline{v}&\text{when
}z_{2}\leq\underline{v},\\\ z_{2}&\text{if }\underline{v}\leq
z_{2}\leq\overline{v},\\\ \overline{v}&\text{otherwise.}\end{array}\right.$

Assumption $(\mathrm{H}_{2})$ together with Hardy-Sobolev type inequality
(2.2) we infer that

$f_{i}(x,\mathrm{T}_{1}(z_{1}),\mathrm{T}_{2}(z_{2}))\in\left(W^{1,p_{i}}(\Omega)\right)^{*},\text{
for }i=1,2.$

Then, from Minty-Browder Theorem (see, e.g., [4, Theorem V.15]), it follows
that the auxiliary problem

(2.8)
$\left\\{\begin{array}[]{ll}-\Delta_{p_{1}}u+|u|^{p_{1}-2}u=f_{1}(x,\mathrm{T}_{1}(z_{1}),\mathrm{T}_{2}(z_{2}))&\text{in
}\Omega,\\\
-\Delta_{p_{2}}v+|v|^{p_{2}-2}v=f_{2}(x,\mathrm{T}_{1}(z_{1}),\mathrm{T}_{2}(z_{2}))&\text{in
}\Omega,\\\ \frac{\partial u}{\partial\eta}=\frac{\partial
v}{\partial\eta}=0&\text{on }\partial\Omega,\end{array}\right.$

admits a unique solution $(u,v)\in W^{1,p_{1}}(\Omega)\times
W^{1,p_{2}}(\Omega).$

Let us introduce the operator

$\begin{array}[]{cccc}\mathcal{T}:&L^{p_{1}}(\Omega)\times
L^{p_{2}}(\Omega)&\rightarrow&W^{1,p_{1}}(\Omega)\times W^{1,p_{2}}(\Omega)\\\
&(z_{1},z_{2})&\mapsto&(u,v).\end{array}$

We note from (2.8) that any fixed point of $\mathcal{T}$ within
$[\underline{u},\overline{u}]\times[\underline{v},\overline{v}]$ coincides
with the weak solution of $(\mathrm{P}_{f_{1},f_{2}})$.

Let us show that $\mathcal{T}$ is continuous. Let
$(z_{1,n},z_{2,n})\rightarrow(z_{1},z_{2})$ in $L^{p_{1}}(\Omega)\times
L^{p_{2}}(\Omega)$ for all $n$. Denote
$\left(u_{n},v_{n}\right)=\mathcal{T}(z_{1,n},z_{2,n})$, which reads as

(2.9) $\begin{array}[]{c}\int_{\Omega}|\nabla u_{n}|^{p_{1}-2}\nabla
u_{n}\nabla\varphi_{1}\,\mathrm{d}x+\int_{\Omega}|u_{n}|^{p_{1}-2}u_{n}\varphi_{1}\,\mathrm{d}x\\\
=\int_{\Omega}f_{1}(x,\mathrm{T}_{1}(z_{1,n}),\mathrm{T}_{2}(z_{2,n}))\varphi_{1}\,\mathrm{d}x\end{array}$

and

(2.10) $\begin{array}[]{c}\int_{\Omega}|\nabla v_{n}|^{p_{2}-2}\nabla
v_{n}\nabla\varphi_{2}\,\mathrm{d}x+\int_{\Omega}|v_{n}|^{p_{2}-2}v_{n}\varphi_{2}\,\mathrm{d}x\\\
=\int_{\Omega}f_{2}(x,\mathrm{T}_{1}(z_{1,n}),\mathrm{T}_{2}(z_{2,n}))\varphi_{2}\,\mathrm{d}x\end{array}$

for all $\varphi_{i}\in W^{1,p_{i}}(\Omega),\ i=1,2$. Inserting
$(\varphi_{1},\varphi_{2})=(u_{n},v_{n})$ in (2.9) and (2.10), using
$(\mathrm{H}_{2})$ and (2.2), we get

(2.11)
$\left\|u_{n}\right\|_{1,p_{1}}^{p_{1}}=\int_{\Omega}f_{1}(x,\mathrm{T}_{1}(z_{1,n}),\mathrm{T}_{2}(z_{2,n}))u_{n}\,\mathrm{d}x\leq
C\int_{\Omega}d(x)^{\gamma}u_{n}\,\mathrm{d}x$

and

(2.12)
$\left\|v_{n}\right\|_{1,p_{2}}^{p_{2}}=\int_{\Omega}f_{2}(x,\mathrm{T}_{1}(z_{1,n}),\mathrm{T}_{2}(z_{2,n}))v_{n}\,\mathrm{d}x\leq
C\int_{\Omega}d(x)^{\gamma}v_{n}\,\mathrm{d}x.$

Since $-1<\gamma<0$, by virtue of the Hardy-Sobolev inequality type in (2.2),
the last integrals in (2.11) and (2.12) are finite and it holds

(2.13) $\left\|u_{n}\right\|_{1,p_{1}}^{p_{1}}\leq
C_{1}\left\|u_{n}\right\|_{1,p_{1}}$

and

(2.14) $\left\|v_{n}\right\|_{1,p_{2}}^{p_{2}}\leq
C_{1}\left\|v_{n}\right\|_{1,p_{2}},$

for a certain $C_{1}>0$ independent of $n$. Thus, $\\{u_{n}\\}$ and
$\\{v_{n}\\}$ are bounded in $W^{1,p_{1}}(\Omega)$ and $W^{1,p_{2}}(\Omega)$,
respectively. So, passing to relabeled subsequences, we get

(2.15) $\begin{array}[]{c}(u_{n},v_{n})\rightharpoonup\left(u,v\right)\text{
in }W^{1,p_{1}}(\Omega)\times W^{1,p_{2}}(\Omega),\end{array}$

for certain $\left(u,v\right)$ in $W^{1,p_{1}}(\Omega)\times
W^{1,p_{2}}(\Omega)$. Setting $\varphi_{1}=u_{n}-u$ in (2.9) we find that

$\begin{array}[]{l}\int_{\Omega}|\nabla u_{n}|^{p_{1}-2}\nabla
u_{n}\nabla(u_{n}-u)\,\mathrm{d}x+\int_{\Omega}|u_{n}|^{p_{1}-2}u_{n}(u_{n}-u)\,\mathrm{d}x\\\
=\int_{\Omega}f_{1}(x,\mathrm{T}_{1}(z_{1,n}),\mathrm{T}_{2}(z_{2,n}))(u_{n}-u)\
\mathrm{d}x\end{array}$

Note that $(\mathrm{H}_{2})$ as well as (2.2) ensure that

(2.16) $f_{1}(x,\mathrm{T}_{1}(z_{1,n}),\mathrm{T}_{2}(z_{2,n}))(u_{n}-u)\in
L^{1}(\Omega).$

Thus, Fatou’s Lemma implies

$\begin{array}[]{l}\underset{n\rightarrow\infty}{\lim}\sup\int_{\Omega}f_{1}(x,\mathrm{T}_{1}(z_{1,n}),\mathrm{T}_{2}(z_{2,n}))(u_{n}-u)\
dx\\\
\leq\int_{\Omega}\underset{n\rightarrow\infty}{\lim}\sup\left(f_{1}(x,\mathrm{T}_{1}(z_{1,n}),\mathrm{T}_{2}(z_{2,n}))(u_{n}-u)\right)\
dx\rightarrow 0,\end{array}$

showing that

$\underset{n\rightarrow\infty}{\lim}\sup\left\langle-\Delta_{p_{1}}u_{n}+|u_{n}|^{p_{1}-2}u_{n},u_{n}-u\right\rangle\leq
0.$

Likewise, we prove that

$\underset{n\rightarrow\infty}{\lim}\sup\left\langle-\Delta_{p_{2}}v_{n}+|v_{n}|^{p_{2}-2}v_{n},v_{n}-v\right\rangle\leq
0.$

Then, the $S_{+}$-property of $-\Delta_{p_{i}}$ on $W^{1,p_{i}}(\Omega)$ (see,
e.g., [23, Proposition 2.72]) along with (2.15) implies

(2.17) $\begin{array}[]{c}(u_{n},v_{n})\rightarrow(u,v)\text{ in
}W^{1,p_{1}}(\Omega)\times W^{1,p_{2}}(\Omega).\end{array}$

Through (2.9), (2.10) and the invariance of $L^{p_{1}}(\Omega)\times
L^{p_{2}}(\Omega)$ by $\mathcal{T}$, we infer that
$\left(u,v\right)=\mathcal{T}(z_{1},z_{2})$, as desired.

Let us verify that $\mathcal{T}(L^{p_{1}}(\Omega)\times L^{p_{2}}(\Omega))$ is
a relatively compact subset. If $(u_{n},v_{n}):=\mathcal{T}(y_{1,n},y_{2,n}),$
$n\in\mathbb{N}$, (2.9) and (2.10) can be written. Hence, the previous
argument yields a pair $(u,v)\in W^{1,p_{1}}(\Omega)\times
W^{1,p_{2}}(\Omega)$ fulfilling (2.17), possibly along a subsequence.

We are thus in a position to apply Schauder’s fixed point theorem to the map
$\mathcal{T}$, which establishes the existence of $(u,v)\in
W^{1,p_{1}}(\Omega)\times W^{1,p_{2}}(\Omega)$ satisfying
$(u,v)=\mathcal{T}(u,v).$ Due to [7, Theorem 3], one has

$\frac{\partial u}{\partial\eta}=\frac{\partial v}{\partial\eta}=0\;\;\text{on
}\partial\Omega.$

Hence, $(u,v)$ is a solution of $(\mathrm{P}_{f_{1},f_{2}})$. Let us show that
(2.6) is fulfilled. Put $\zeta=(\underline{u}-u)^{+}$ and suppose $\zeta\neq
0$. Then, from $(\mathrm{H}_{1})$, (2.7) and (2.8) , we get

$\begin{array}[]{l}\int_{\\{u<\underline{u}\\}}|\nabla u|^{p_{1}-2}\nabla
u\nabla\zeta\ \mathrm{d}x+\int_{\\{u<\underline{u}\\}}|u|^{p_{1}-2}u\zeta\
\mathrm{d}x\\\ =\int_{\Omega}|\nabla u|^{p_{1}-2}\nabla u\nabla\zeta\
\mathrm{d}x+\int_{\Omega}|u|^{p_{1}-2}u\zeta\
\mathrm{d}x=\int_{\Omega}f_{1}(x,\mathrm{T}_{1}(u),\mathrm{T}_{2}(v))\zeta\
\mathrm{d}x\\\
=\int_{\\{u<\underline{u}\\}}f_{1}(x,\mathrm{T}_{1}(u),\mathrm{T}_{2}(v))\zeta\
\mathrm{d}x=\int_{\\{u<\underline{u}\\}}f_{1}(x,\underline{u},\mathrm{T}_{2}(v))\zeta\
\mathrm{d}x\\\
\geq\int_{\\{u<\underline{u}\\}}|\nabla\underline{u}|^{p_{1}-2}\nabla\underline{u}\nabla\zeta\
\mathrm{d}x+\int_{\\{u<\underline{u}\\}}|\underline{u}|^{p_{1}-2}\underline{u}\zeta\
\mathrm{d}x-\int_{\partial\Omega}|\nabla\underline{u}|^{p_{1}-2}\frac{\partial\underline{u}}{\partial\eta}\gamma_{0}(\zeta)\text{
}\mathrm{d}s.\end{array}$

which by $(\mathrm{H}_{1})$ and (2.1), implies that

$\begin{array}[]{c}\int_{\\{u<\underline{u}\\}}(|\nabla\underline{u}|^{p_{1}-2}\nabla\underline{u}-|\nabla
u|^{p_{1}-2}\nabla u)\nabla\zeta\
\mathrm{d}x+\int_{\\{u<\underline{u}\\}}(|\underline{u}|^{p_{1}-2}\underline{u}-|u|^{p_{1}-2}u)\zeta\
\mathrm{d}x\\\
\leq\int_{\partial\Omega}|\nabla\underline{u}|^{p_{1}-2}\frac{\partial\underline{u}}{\partial\eta}\gamma_{0}(\zeta)\text{
}\mathrm{d}s\leq 0,\end{array}$

a contradiction. Hence $u\geq\underline{u}$ in $\Omega$. Arguing similarly,
set $\hat{\zeta}=(u-\overline{u})^{+}$ and assume that $\hat{\zeta}\neq 0.$
Then

$\begin{array}[]{l}\int_{\\{u>\overline{u}\\}}|\nabla u|^{p_{1}-2}\nabla
u\nabla\hat{\zeta}\
\mathrm{d}x+\int_{\\{u>\overline{u}\\}}|u|^{p_{1}-2}u\hat{\zeta}\
\mathrm{d}x\\\ =\int_{\Omega}|\nabla u|^{p_{1}-2}\nabla u\nabla\hat{\zeta}\
\mathrm{d}x+\int_{\Omega}|u|^{p_{1}-2}u\hat{\zeta}\
\mathrm{d}x=\int_{\Omega}f_{1}(x,\mathrm{T}_{1}(u),\mathrm{T}_{2}(v))\hat{\zeta}\
\mathrm{d}x\\\
=\int_{\\{u>\overline{u}\\}}f_{1}(x,\mathrm{T}_{1}(u),\mathrm{T}_{2}(v))\hat{\zeta}\
\mathrm{d}x=\int_{\\{u>\overline{u}\\}}f_{1}(x,\overline{u},\mathrm{T}_{2}(v))\hat{\zeta}\
\mathrm{d}x\\\
\leq\int_{\\{u>\overline{u}\\}}|\nabla\overline{u}|^{p_{1}-2}\nabla\overline{u}\nabla\hat{\zeta}\
\mathrm{d}x+\int_{\\{u>\overline{u}\\}}|\overline{u}|^{p_{1}-2}\overline{u}\hat{\zeta}\
\mathrm{d}x-\int_{\partial\Omega}|\nabla\overline{u}|^{p_{1}-2}\frac{\partial\overline{u}}{\partial\eta}\gamma_{0}(\hat{\zeta})\text{
}\mathrm{d}s\end{array}$

which leads to

$\begin{array}[]{c}\int_{\\{u>\overline{u}\\}}(|\nabla u|^{p_{1}-2}\nabla
u-|\nabla\overline{u}|^{p_{1}-2}\nabla\overline{u})\nabla\hat{\zeta}\
\mathrm{d}x+\int_{\\{u>\overline{u}\\}}(|u|^{p_{1}-2}u-|\overline{u}|^{p_{1}-2}\overline{u})\hat{\zeta}\
\mathrm{d}x\\\
\leq-\int_{\partial\Omega}|\nabla\overline{u}|^{p_{1}-2}\frac{\partial\overline{u}}{\partial\eta}\gamma_{0}(\hat{\zeta})\text{
}\mathrm{d}s\leq 0,\end{array}$

a contradiction. Thus, we have $u\leq\overline{u}$ in $\Omega$. A quite
similar argument provides that $\underline{v}\leq v\leq\overline{v}$ in
$\Omega$. Therefore, since
$(\underline{u},\underline{v}),(\overline{u},\overline{v})\in
L^{\infty}(\Omega)\times L^{\infty}(\Omega)$, we conclude that $(u,v)\in
W_{b}^{1,p_{1}}(\Omega)\times W_{b}^{1,p_{2}}(\Omega)$. This completes the
proof. ∎

Instead of $(\mathrm{H}_{2}),$ if we assume that $f_{1}$ and $f_{2}$ are
bounded in
$\Omega\times[\underline{u},\overline{u}]\times[\underline{v},\overline{v}]$,
the conclusion of Theorem 2.1 remains true but assigning more regularity to
the solution.

  * $(\mathrm{H}_{3})$

There exists a constant $M>0$ such that

$\left|f_{i}(x,u,v)\right|\leq M\text{ in
}\Omega\times[\underline{u},\overline{u}]\times[\underline{v},\overline{v}]\text{,
for }i=1,2.$

######  Theorem 2.2.

Suppose $(\mathrm{H}_{1})$ and $(\mathrm{H}_{3})$ hold true. Then, problem
$(\mathrm{P}_{f_{1},f_{2}})$ possesses a solution
$(u,v)\in\mathcal{C}^{1,\tau}(\overline{\Omega})\times\mathcal{C}^{1,\tau}(\overline{\Omega})$
with suitable $\tau\in(0,1)$ such that

(2.18) $\underline{u}\leq
u\leq\overline{u}\quad\text{and}\quad\underline{v}\leq v\leq\overline{v}.$

Moreover, $\frac{\partial u}{\partial\eta}=\frac{\partial v}{\partial\eta}=0$
on $\partial\Omega$.

###### Proof.

The proof is similar in the spirit to the one of Theorem 2.1. Given
$(z_{1},z_{2})\in\mathcal{C}(\overline{\Omega})\times\mathcal{C}(\overline{\Omega})$,
we introduce the auxiliary system

(2.19)
$\left\\{\begin{array}[]{ll}-\Delta_{p_{1}}u+|u|^{p_{1}-2}u=f_{1}(x,\mathrm{T}_{1}(z_{1}),\mathrm{T}_{2}(z_{2}))&\text{in
}\Omega,\\\
-\Delta_{p_{2}}v+|v|^{p_{2}-2}v=f_{2}(x,\mathrm{T}_{1}(z_{1}),\mathrm{T}_{2}(z_{2}))&\text{in
}\Omega,\\\ \frac{\partial u}{\partial\eta}=\frac{\partial
v}{\partial\eta}=0&\text{on }\partial\Omega,\end{array}\right.$

where the operators $\mathrm{T}_{1}$ and $\mathrm{T}_{2}$ are defined by
(2.7). Notice that $(\mathrm{H}_{3})$ together with Minty-Browder Theorem
(see, e.g., [4, Theorem V.15]) ensure that (2.19) admits a unique solution
$(u,v)\in W^{1,p_{1}}(\Omega)\times W^{1,p_{2}}(\Omega).$ Let us introduce the
operator

$\begin{array}[]{cccc}\mathcal{T}:&\mathcal{C}(\overline{\Omega})\times\mathcal{C}(\overline{\Omega})&\rightarrow&\mathcal{C}^{1}(\overline{\Omega})\times\mathcal{C}^{1}(\overline{\Omega})\\\
&(z_{1},z_{2})&\mapsto&(u,v).\end{array}$

Observe from (2.19) that any fixed point of $\mathcal{T}$ within
$[\underline{u},\overline{u}]\times[\underline{v},\overline{v}]$ coincides
with the weak solution of $(\mathrm{P}_{f_{1},f_{2}})$.

By $(\mathrm{H}_{3})$ and according to the regularity result [22], it follows
that
$(u,v)\in\mathcal{C}^{1,\tau}(\overline{\Omega})\times\mathcal{C}^{1,\tau}(\overline{\Omega})$
and

(2.20)
$\left\|u\right\|_{\mathcal{C}^{1,\tau}(\overline{\Omega})},\left\|v\right\|_{\mathcal{C}^{1,\tau}(\overline{\Omega})}\leq
L_{0},$

for some constant $L_{0}>0$ independent of $u$ and $v$. Then, the compactness
of the embedding
$\mathcal{C}^{1,\tau}(\overline{\Omega})\subset\mathcal{C}^{1}(\overline{\Omega})$
implies that
$\mathcal{T}(\mathcal{C}(\overline{\Omega})\times\mathcal{C}(\overline{\Omega}))$
is a relatively compact subset of
$\mathcal{C}^{1}(\overline{\Omega})\times\mathcal{C}^{1}(\overline{\Omega})$.
This proves that $\mathcal{T}$ is compact.

Next, we show the continuity of $\mathcal{T}$ with respect to the topology of
$\mathcal{C}(\overline{\Omega})\times\mathcal{C}(\overline{\Omega})$. Let
$(z_{1,n},z_{2,n})\rightarrow(z_{1},z_{2})$ in
$\mathcal{C}(\overline{\Omega})\times\mathcal{C}(\overline{\Omega})$ for all
$n$ and denote $\left(u_{n},v_{n}\right)=\mathcal{T}(z_{1,n},z_{2,n})$.
Repeating the previous argument in the proof of Theorem 2.1 we obtain

$\begin{array}[]{c}(u_{n},v_{n})\rightarrow(u,v)\text{ in
}W^{1,p_{1}}(\Omega)\times W^{1,p_{2}}(\Omega),\end{array}$

for certain $\left(u,v\right)$ in $W^{1,p_{1}}(\Omega)\times
W^{1,p_{2}}(\Omega)$ and $\left(u,v\right)=\mathcal{T}(z_{1},z_{2})$. On the
basis of [22], the sequence $\\{\left(u_{n},v_{n}\right)\\}$ is bounded in
$\mathcal{C}^{1,\tau}(\overline{\Omega})\times\mathcal{C}^{1,\tau}(\overline{\Omega})$
for certain $\tau\in(0,1)$. Then, through the compact embedding
$\mathcal{C}^{1,\tau}(\overline{\Omega})\subset\mathcal{C}^{1}(\overline{\Omega})$,
along a relabeled subsequence, there holds $(u_{n},v_{n})\rightarrow(u,v)$ in
$\mathcal{C}^{1}(\overline{\Omega})\times\mathcal{C}^{1}(\overline{\Omega})$,
showing that $\mathcal{T}$ is continuous.

We are thus in a position to apply Schauder’s fixed point theorem to the map
$\mathcal{T}$, which establishes the existence of
$(u,v)\in\mathcal{C}^{1}(\overline{\Omega})\times\mathcal{C}^{1}(\overline{\Omega})$
satisfying $(u,v)=\mathcal{T}(u,v).$ By $(\mathrm{H}_{3})$ together with [22],
we infer that
$(u,v)\in\mathcal{C}^{1,\tau}(\overline{\Omega})\times\mathcal{C}^{1,\tau}(\overline{\Omega}),$
$\tau\in(0,1)$. The rest of the proof runs as the one of Theorem 2.1. ∎

######  Remark 2.3.

Theorems 2.1 and 2.2 remain valid if we replace Neumann boundary conditions
with Dirichlet ones.

##  3\. Abstract results

######  Theorem 3.1.

For a constant $\Lambda>0$, let
$(\underline{u}_{\Lambda},\underline{v}_{\Lambda})$ and
$(\overline{u}_{\Lambda},\overline{v}_{\Lambda})$ be sub-supersolutions pairs
of problem $(\mathrm{P})$. Assume (1.1) holds and suppose there exists a
constant $\rho>0$ such that

(3.1) $\underline{u}_{\Lambda},\underline{v}_{\Lambda}>\rho\text{ a.e. in
}\overline{\Omega}.$

Then, problem $(\mathrm{P})$ admits a solution $(u,v)\in
int\mathcal{C}_{+}^{1,\tau}(\overline{\Omega})\times
int\mathcal{C}_{+}^{1,\tau}(\overline{\Omega}),$ for certain $\tau\in(0,1)$,
within
$[\underline{u}_{\Lambda},\overline{u}_{\Lambda}]\times[\underline{v}_{\Lambda},\overline{v}_{\Lambda}]$.
Moreover, if $\alpha_{2},\beta_{1}\in(-1,0)$ and

(3.2)
$\max\\{\frac{-\gamma_{1}\alpha_{i}}{p_{1}-1},\frac{-\gamma_{2}\beta_{i}}{p_{2}-1}\\}<p_{i}-1,\text{
}i=1,2,$

with

$\gamma_{i}=\max\\{-\alpha_{i},-\beta_{i}\\},$

then, $(u,v)$ is unique.

###### Proof.

By (1.1), (3.1) and for all
$(u,v)\in[\underline{u}_{\Lambda},\overline{u}_{\Lambda}]\times[\underline{v}_{\Lambda},\overline{v}_{\Lambda}],$
we have

$u^{\alpha_{1}}+v^{\beta_{1}}\leq\left\\{\begin{array}[]{ll}\underline{u}_{\Lambda}^{\alpha_{1}}+\overline{v}_{\Lambda}^{\beta_{1}}&\text{
if }\beta_{1}>0\\\
\underline{u}_{\Lambda}^{\alpha_{1}}+\underline{v}_{\Lambda}^{\beta_{1}}&\text{
if
}\beta_{1}<0\end{array}\right.\leq\left\\{\begin{array}[]{ll}\rho^{\alpha_{1}}+\left\|\overline{v}_{\Lambda}\right\|_{\infty}^{\beta_{1}}&\text{
if }\beta_{1}>0\\\ \rho^{\alpha_{1}}+\rho^{\beta_{1}}&\text{ if
}\beta_{1}<0,\end{array}\right.$

as well as

$u^{\alpha_{2}}+v^{\beta_{2}}\leq\left\\{\begin{array}[]{ll}\overline{u}_{\Lambda}^{\alpha_{2}}+\underline{v}_{\Lambda}^{\beta_{2}}&\text{
if }\alpha_{2}>0\\\
\underline{u}_{\Lambda}^{\alpha_{2}}+\underline{v}_{\Lambda}^{\beta_{2}}&\text{
if
}\alpha_{2}<0\end{array}\right.\leq\left\\{\begin{array}[]{ll}\left\|\overline{u}_{\Lambda}\right\|_{\infty}^{\alpha_{2}}+\rho^{\beta_{2}}&\text{
if }\alpha_{2}>0\\\ \rho^{\alpha_{2}}+\rho^{\beta_{2}}&\text{ if
}\alpha_{2}<0.\end{array}\right.$

Then, thanks to Theorem 2.2, we conclude that $(\mathrm{P})$ possesses a
solution
$(u,v)\in\mathcal{C}^{1,\tau}(\overline{\Omega})\times\mathcal{C}^{1,\tau}(\overline{\Omega}),$
with suitable $\tau\in(0,1),$ within
$[\underline{u}_{\Lambda},\overline{u}_{\Lambda}]\times[\underline{v}_{\Lambda},\overline{v}_{\Lambda}]$.
Due to (3.1), we deduce that $(u,v)\in
int\mathcal{C}_{+}^{1,\tau}(\overline{\Omega})\times
int\mathcal{C}_{+}^{1,\tau}(\overline{\Omega})$.

We proceed to show that $(u,v)$ is a unique solution in
$[\underline{u}_{\Lambda},\overline{u}_{\Lambda}]\times[\underline{v}_{\Lambda},\overline{v}_{\Lambda}]$
for $\alpha_{2},\beta_{1}\in(-1,0)$. To this end, let $(u_{1},v_{1})$ and
$(u_{2},v_{2})$ be two distinct positive solutions of $(\mathrm{P})$ within
$[\underline{u}_{\Lambda},\overline{u}_{\Lambda}]\times[\underline{v}_{\Lambda},\overline{v}_{\Lambda}]$.
Set

$\tau=\sup\\{c\in\mathbb{R}_{+}\text{, \ }cu_{2}\leq u_{1}\text{ \ and
}cv_{2}\leq v_{1}\text{ \ in\ }\Omega\\}.$

Then $0<\tau<\infty$ because

(3.3)
$\min\\{\frac{\rho}{||u_{2}||_{\infty}},\frac{\rho}{||v_{2}||_{\infty}}\\}\leq\tau\leq\max\\{\frac{||u_{1}||_{\infty}}{\rho},\frac{||v_{1}||_{\infty}}{\rho}\\}$

If we manage to show that $\tau\geq 1$, we are done as this entails $u_{2}\leq
u_{1}$ and $v_{2}\leq v_{1}$ in $\Omega$ and thus, by interchanging the roles
of $(u_{1},v_{1})$ and $(u_{2},v_{2})$ we get $u_{2}\geq u_{1}$ and $v_{2}\geq
v_{1}$ in $\Omega$. By contradiction, suppose that $0<\tau<1.$ Then, by (1.1)
and since $\alpha_{2},\beta_{1}<0$, we infer that

$\min\\{\tau^{-\alpha_{i}},\tau^{-\beta_{i}}\\}\geq\tau^{\gamma_{i}}\text{ \
and \
}\min\\{\tau^{-\frac{\gamma_{1}\alpha_{i}}{p_{1}-1}},\tau^{-\frac{\gamma_{2}\beta_{i}}{p_{2}-1}}\\}\geq\tau^{\hat{\gamma}_{i}}\text{,}$

with

$\hat{\gamma}_{i}=\max\\{\frac{-\gamma_{1}\alpha_{i}}{p_{1}-1},\frac{-\gamma_{2}\beta_{i}}{p_{2}-1}\\}>0,\text{for
}i=1,2.$

A direct computations show that

$\begin{array}[]{l}-\Delta_{p_{1}}u_{2}+|u_{2}|^{p_{1}-2}u_{2}=u_{2}^{\alpha_{1}}+v_{2}^{\beta_{1}}=\frac{(\tau
u_{2})^{\alpha_{1}}}{\tau^{\alpha_{1}}}+\frac{(\tau
v_{2})^{\beta_{1}}}{\tau^{\beta_{1}}}\\\ \geq\tau^{\gamma_{1}}\left((\tau
u_{2})^{\alpha_{1}}+(\tau
v_{2})^{\beta_{1}}\right)\geq\tau^{\gamma_{1}}(u_{1}^{\alpha_{1}}+v_{1}^{\beta_{1}})\\\
=\tau^{\gamma_{1}}\left(-\Delta_{p_{1}}u_{1}+|u_{1}|^{p_{1}-2}u_{1}\right)\\\
=-\Delta_{p_{1}}(\tau^{\frac{\gamma_{1}}{p_{1}-1}}u_{1})+|\tau^{\frac{\gamma_{1}}{p_{1}-1}}u_{1}|^{p_{1}-2}(\tau^{\frac{\gamma_{1}}{p_{1}-1}}u_{1})\end{array}$

and similarly

$\begin{array}[]{l}-\Delta_{p_{2}}v_{2}+|v_{2}|^{p_{2}-2}v_{2}=u_{2}^{\alpha_{2}}+v_{2}^{\beta_{2}}=\frac{(\tau
u_{2})^{\alpha_{2}}}{\tau^{\alpha_{2}}}+\frac{(\tau
v_{2})^{\beta_{2}}}{\tau^{\beta_{2}}}\\\ \geq\tau^{\gamma_{2}}\left((\tau
u_{2})^{\alpha_{2}}+(\tau
v_{2})^{\beta_{2}}\right)\geq\tau^{\gamma_{2}}(u_{1}^{\alpha_{2}}+v_{1}^{\beta_{2}})\\\
=-\Delta_{p_{2}}(\tau^{\frac{\gamma_{2}}{p_{2}-1}}v_{1})+|\tau^{\frac{\gamma_{2}}{p_{2}-1}}v_{1}|^{p_{2}-2}(\tau^{\frac{\gamma_{2}}{p_{2}-1}}v_{1}).\end{array}$

The weak comparison principle (see [30, Lemma 3.2]) yields

(3.4) $u_{2}\geq\tau^{\frac{\gamma_{1}}{p_{1}-1}}u_{1}\text{ \ and \
}v_{2}\geq\tau^{\frac{\gamma_{2}}{p_{2}-1}}v_{1}\text{ \ in }\Omega\text{.}$

Using (3.4) in the equations for $(u_{1},v_{1})$, we get

$\begin{array}[]{l}-\Delta_{p_{1}}u_{1}+|u_{1}|^{p_{1}-2}u_{1}=u_{1}^{\alpha_{1}}+v_{1}^{\beta_{1}}=\frac{(\tau^{\frac{\gamma_{1}}{p_{1}-1}}u_{1})^{\alpha_{1}}}{\tau^{\frac{\gamma_{1}\alpha_{1}}{p_{1}-1}}}+\frac{(\tau^{\frac{\gamma_{2}}{p_{2}-1}}v_{1})^{\beta_{1}}}{\tau^{\frac{\gamma_{2}\beta_{1}}{p_{2}-1}}}\\\
\geq\tau^{\hat{\gamma}_{1}}(u_{2}^{\alpha_{1}}+v_{2}^{\beta_{1}})=\tau^{\hat{\gamma}_{1}}\left(-\Delta_{p_{1}}u_{2}+|u_{2}|^{p_{1}-2}u_{2}\right)\\\
=-\Delta_{p_{1}}(\tau^{\frac{\hat{\gamma}_{1}}{p_{1}-1}}u_{2})+|\tau^{\frac{\hat{\gamma}_{1}}{p_{1}-1}}u_{2}|^{p_{1}-2}(\tau^{\frac{\hat{\gamma}_{1}}{p_{1}-1}}u_{2}).\end{array}$

and

$\begin{array}[]{l}-\Delta_{p_{2}}v_{1}+|v_{1}|^{p_{2}-2}v_{1}=u_{1}^{\alpha_{2}}+v_{1}^{\beta_{2}}=\frac{(\tau^{\frac{\gamma_{1}}{p_{1}-1}}u_{1})^{\alpha_{2}}}{\tau^{\frac{\gamma_{1}\alpha_{2}}{p_{1}-1}}}+\frac{(\tau^{\frac{\gamma_{2}}{p_{2}-1}}v_{1})^{\beta_{2}}}{\tau^{\frac{\gamma_{2}\beta_{2}}{p_{2}-1}}}\\\
\geq\tau^{\hat{\gamma}_{2}}(u_{2}^{\alpha_{2}}+v_{2}^{\beta_{2}})=-\Delta_{p_{2}}(\tau^{\frac{\hat{\gamma}_{2}}{p_{2}-1}}v_{2})+|\tau^{\frac{\hat{\gamma}_{2}}{p_{2}-1}}v_{2}|^{p_{2}-2}(\tau^{\frac{\hat{\gamma}_{2}}{p_{2}-1}}v_{2})\text{
in }\Omega.\end{array}$

Owing to [30, Lemma 3.2] we derive

$u_{1}\geq\tau^{\frac{\hat{\gamma}_{1}}{p_{1}-1}}u_{2}\text{ \ and \
}v_{1}\geq\tau^{\frac{\hat{\gamma}_{2}}{p_{2}-1}}v_{2}\text{ \ in }\Omega.$

In view of (3.2), we deduce that

$1>\tau^{\frac{\hat{\gamma}_{1}}{p_{1}-1}}>\tau\text{ \ and \
}1>\tau^{\frac{\hat{\gamma}_{2}}{p_{2}-1}}>\tau,$

a contradiction with the definition of $\tau$. Thus, $\tau\geq 1$ and
therefore, $(u_{1},v_{1})=(u_{2},v_{2})$, completing the proof of the theorem.
∎

The existence result in the case of zero trace subsolutions on the boundary
$\partial\Omega$ is formulated as follows.

######  Theorem 3.2.

For a constant $\Lambda>0$, let
$(\underline{u}_{\Lambda},\underline{v}_{\Lambda})$ and
$(\overline{u}_{\Lambda},\overline{v}_{\Lambda})$ be sub-supersolutions pairs
of problem $(\mathrm{P})$ and assume that assumption (1.1) holds. If
$(\underline{u}_{\Lambda},\underline{v}_{\Lambda})\in
W_{0,b}^{1,p_{1}}(\Omega)\times W_{0,b}^{1,p_{2}}(\Omega)$ and there exists a
constant $c>0$ such that

(3.5) $\underline{u}_{\Lambda},\underline{v}_{\Lambda}\geq cd(x)\text{ a.e. in
}\Omega,$

then, problem $(\mathrm{P})$ admits a solution $(u,v)\in
W_{b}^{1,p_{1}}(\Omega)\times W_{b}^{1,p_{2}}(\Omega)$ within
$[\underline{u}_{\Lambda},\overline{u}_{\Lambda}]\times[\underline{v}_{\Lambda},\overline{v}_{\Lambda}]$.

###### Proof.

By (1.1) and (3.5), for all
$(u,v)\in[\underline{u}_{\Lambda},\overline{u}_{\Lambda}]\times[\underline{v}_{\Lambda},\overline{v}_{\Lambda}],$
we have

$\displaystyle u^{\alpha_{1}}+v^{\beta_{1}}$ $\displaystyle\leq$
$\displaystyle\left\\{\begin{array}[]{ll}\underline{u}_{\Lambda}^{\alpha_{1}}+\overline{v}_{\Lambda}^{\beta_{1}}&\text{
if }\beta_{1}>0\\\
\underline{u}_{\Lambda}^{\alpha_{1}}+\underline{v}_{\Lambda}^{\beta_{1}}&\text{
if }\beta_{1}<0\end{array}\right.$ $\displaystyle\leq$
$\displaystyle\left\\{\begin{array}[]{ll}(cd(x))^{\alpha_{1}}+\left\|\overline{v}_{\Lambda}\right\|_{\infty}^{\beta_{1}}&\text{
if }\beta_{1}>0\\\ (cd(x))^{\alpha_{1}}+(cd(x))^{\beta_{1}}&\text{ if
}\beta_{1}<0\end{array}\right.$ $\displaystyle\leq$
$\displaystyle\left\\{\begin{array}[]{ll}(c^{\alpha_{1}}+d(x)^{-\alpha_{1}}\left\|\overline{v}_{\Lambda}\right\|_{\infty}^{\beta_{1}})d(x)^{\alpha_{1}}&\text{
if }\beta_{1}>0\\\
(c^{\alpha_{1}}+c^{\beta_{1}})\max\\{d(x))^{\alpha_{1}},d(x))^{\beta_{1}}\\}&\text{
if }\beta_{1}<0\end{array}\right.$ $\displaystyle\leq$ $\displaystyle
C_{0}\left\\{\begin{array}[]{ll}d(x)^{\alpha_{1}}&\text{ if }\beta_{1}>0\\\
\max\\{d(x)^{\alpha_{1}},d(x)^{\beta_{1}}\\}&\text{ if
}\beta_{1}<0\end{array}\right.\text{, in }\Omega$

and similarly

$\displaystyle u^{\alpha_{2}}+v^{\beta_{2}}$ $\displaystyle\leq$
$\displaystyle\left\\{\begin{array}[]{ll}\overline{u}_{\Lambda}^{\alpha_{2}}+\underline{v}_{\Lambda}^{\beta_{2}}&\text{
if }\alpha_{2}>0\\\
\underline{u}_{\Lambda}^{\alpha_{2}}+\underline{v}_{\Lambda}^{\beta_{2}}&\text{
if
}\alpha_{2}<0\end{array}\right.\leq\left\\{\begin{array}[]{ll}\left\|\overline{u}_{\Lambda}\right\|_{\infty}^{\alpha_{2}}+(cd(x))^{\beta_{2}}&\text{
if }\alpha_{2}>0\\\ (cd(x))^{\alpha_{2}}+(cd(x))^{\beta_{2}}&\text{ if
}\alpha_{2}<0.\end{array}\right.$ $\displaystyle\leq$
$\displaystyle\tilde{C}_{0}\left\\{\begin{array}[]{ll}d(x)^{\beta_{2}}&\text{
if }\alpha_{2}>0\\\ \max\\{d(x)^{\alpha_{2}},d(x)^{\beta_{2}}\\}&\text{ if
}\alpha_{2}<0\end{array}\right.\text{, in }\Omega,$

for certain constants $C_{0},\tilde{C}_{0}>0$. Then, thanks to Theorem 2.1, we
conclude that there exists a solution $(u,v)\in W_{b}^{1,p_{1}}(\Omega)\times
W_{b}^{1,p_{2}}(\Omega)$ of problem $(\mathrm{P})$ within
$[\underline{u}_{\Lambda},\overline{u}_{\Lambda}]\times[\underline{v}_{\Lambda},\overline{v}_{\Lambda}]$.
This proves the theorem. ∎

######  Remark 3.3.

It is worth noting that the supersolution
$(\overline{u}_{\Lambda},\overline{v}_{\Lambda})$ in Theorem 3.2 should not
belong to $W_{0}^{1,p_{1}}(\Omega)\times W_{0}^{1,p_{2}}(\Omega)$ because if
so, the solution $(u,v)$ would be of zero trace on the boundary
$\partial\Omega.$ Hence, [17, Lemma 3.1] ensures that
$(u,v)\in\mathcal{C}^{1,\tau}(\overline{\Omega})\times\mathcal{C}^{1,\tau}(\overline{\Omega})$
which, by Hopf’s Lemma (see, e.g., [1]), implies that
$\frac{\partial\overline{u}}{\partial\eta},\frac{\partial\overline{v}}{\partial\eta}<0$
on $\partial\Omega$, absurd.

######  Remark 3.4.

The argument used to show the uniqueness result of Theorem 3.1 is not
applicable in the context of Theorem 3.2. Being impossible to construct
supersolutions for problem $(\mathrm{P})$ behaving like the distance function
$d(x)$ (see Remark 3.3), then it would not be possible to get an estimate of
type (3.3).

##  4\. Existence and uniqueness results

Our goal is to construct sub- and super-solution pairs of $(\mathrm{P})$. With
this aim, consider the following nonlinear Dirichlet and Neumann eigenvalue
problems

(4.1)
$\left\\{\begin{array}[]{l}-\Delta_{p_{i}}\phi_{1,p_{i}}+|\phi_{1,p_{i}}|^{p_{i}-2}\phi_{1,p_{i}}=\lambda_{1,p_{i}}|\phi_{1,p_{i}}|^{p_{i}-2}\phi_{1,p_{i}}\text{
in }\Omega\\\ \phi_{1,p_{i}}=0\text{ on }\partial\Omega,\text{
}i=1,2,\end{array}\right.$ (4.2)
$\left\\{\begin{array}[]{l}-\Delta_{p_{i}}\hat{\phi}_{1,p_{i}}+|\hat{\phi}_{1,p_{i}}|^{p_{i}-2}\hat{\phi}_{1,p_{i}}=\hat{\lambda}_{1,p_{i}}|\hat{\phi}_{1,p_{i}}|^{p_{i}-2}\hat{\phi}_{1,p_{i}}\text{
in }\Omega\\\ \frac{\partial\hat{\phi}_{1,p_{i}}}{\partial\eta}=0\text{ on
}\partial\Omega,\text{ }i=1,2,\end{array}\right.$

where $\phi_{1,p_{i}}\in\mathcal{C}_{+}^{1}(\overline{\Omega})$ and
$\hat{\phi}_{1,p_{i}}\in int\mathcal{C}_{+}^{1}(\overline{\Omega})$ are the
eigenfunctions corresponding to the first eigenvalues
$\lambda_{1,p_{i}},\hat{\lambda}_{1,p_{i}}>0,$ respectively (see [23]). Recall
that

(4.3) $\phi_{1,p_{i}}\geq c_{0}d(x)\text{ in }\Omega,\text{ \ and \
}\frac{\partial\phi_{1,p_{i}}}{\partial\eta}<0\;\text{on}\;\partial\Omega,$

for a constant $c_{0}>0,$ and there exists a constant $\mu>0$ such that

(4.4) $\hat{\phi}_{1,p_{i}}(x)>\mu\text{ for all }x\in\overline{\Omega}.$

Consider the homogeneous Dirichlet and Neumann problems

(4.5)
$\left\\{\begin{array}[]{l}-\Delta_{p_{i}}y_{i}+|y_{i}|^{p_{i}-2}y_{i}=1\text{
in }\Omega,\\\ y_{i}=0\text{ on }\partial\Omega,\text{
}i=1,2\end{array}\right.$ (4.6)
$\left\\{\begin{array}[]{l}-\Delta_{p_{i}}\hat{y}_{i}+|\hat{y}_{i}|^{p_{i}-2}\hat{y}_{i}=1\text{
in }\Omega\\\ \frac{\partial\hat{y}_{i}}{\partial\eta}=0\text{ on
}\partial\Omega,\text{ }i=1,2,\end{array}\right.$

which admit unique solutions $y_{i}\in\mathcal{C}^{1,\tau}(\overline{\Omega})$
and $\hat{y}_{i}\in int\mathcal{C}_{+}^{1}(\overline{\Omega})$ satisfying

(4.7) $\frac{d}{c}\leq y_{i}\leq cd\;\text{in}\;\Omega,\;\;\frac{\partial
y_{i}}{\partial\eta}<0\;\text{on}\;\partial\Omega,$ (4.8)
$\frac{\hat{\phi}_{1,p_{i}}}{\hat{c}}\leq\hat{y}_{i}\leq\hat{c}\hat{\phi}_{1,p_{i}}\;\text{in}\;\Omega,$

for some constants $c,\hat{c}>1$ and $L>0$ (see [9, 30, 32] ).

###  4.1. $C^{1}$\- bound solutions

######  Theorem 4.1.

Assume (1.1) is satisfied. Then, if $\Lambda>0$ is big enough, problem
$(\mathrm{P})$ admits a solution $(u,v)\in
int\mathcal{C}_{+}^{1,\tau}(\overline{\Omega})\times
int\mathcal{C}_{+}^{1,\tau}(\overline{\Omega}),$ for certain $\tau\in(0,1)$,
such that

(4.9)
$(u,v)\in[\Lambda^{-1}\hat{\phi}_{1,p_{1}},\Lambda\hat{y}_{1}]\times[\Lambda^{-1}\hat{\phi}_{1,p_{2}},\Lambda\hat{y}_{2}].$

which is unique once $\alpha_{2},\beta_{1}\in(-1,0)$ and assumption (3.2) is
fulfilled .

###### Proof.

For a constant $\Lambda>0$ which will be specified later, let us show that
$\Lambda(\hat{y}_{1},\hat{y}_{2})$ satisfy (2.5). With this aim, pick
$(u,v)\in W^{1,p_{1}}(\Omega)\times W^{1,p_{2}}(\Omega)$ within
$[\Lambda^{-1}\hat{\phi}_{1,p_{1}},\Lambda\hat{y}_{1}]\times[\Lambda^{-1}\hat{\phi}_{1,p_{2}},\Lambda\hat{y}_{2}]$.
From (4.8) and (1.1), if $\beta_{1},\alpha_{2}>0$, we get

(4.10)
$\begin{array}[]{l}(\Lambda\hat{y}_{1})^{\alpha_{1}}+(\Lambda\hat{y}_{2})^{\beta_{1}}\leq(\Lambda\frac{\hat{\phi}_{1,p_{1}}}{\hat{c}})^{\alpha_{1}}+(\Lambda\hat{c}\hat{\phi}_{1,p_{2}})^{\beta_{1}}\\\
\leq\hat{c}^{-\alpha_{1}}+(\Lambda\hat{c}||\hat{\phi}_{1,p_{2}}||_{\infty})^{\beta_{1}}\leq\Lambda^{\beta_{1}}(1+(\hat{c}||\hat{\phi}_{1,p_{2}}||_{\infty})^{\beta_{1}})\text{
\ in }\Omega,\end{array}$ (4.11)
$\begin{array}[]{l}(\Lambda\hat{y}_{1})^{\alpha_{2}}+(\Lambda\hat{y}_{2})^{\beta_{2}}\leq(\Lambda\hat{c}\hat{\phi}_{1,p_{1}})^{\alpha_{2}}+(\Lambda\frac{\hat{\phi}_{1,p_{2}}}{\hat{c}})^{\beta_{2}}\\\
\leq(\Lambda\hat{c}||\hat{\phi}_{1,p_{1}}||_{\infty})^{\alpha_{2}}+\hat{c}^{-\beta_{2}}\leq\Lambda^{\alpha_{2}}(1+(\hat{c}||\hat{\phi}_{1,p_{1}}||_{\infty})^{\alpha_{2}})\text{
\ in }\Omega,\end{array}$

while if $\beta_{1},\alpha_{2}<0$, we obtain

(4.12)
$\begin{array}[]{l}(\Lambda\hat{y}_{1})^{\alpha_{1}}+(\Lambda^{-1}\hat{\phi}_{1,p_{2}})^{\beta_{1}}\leq(\Lambda\frac{\hat{\phi}_{1,p_{1}}}{\hat{c}})^{\alpha_{1}}+(\Lambda^{-1}\mu)^{\beta_{1}}\\\
\leq\hat{c}^{-\alpha_{1}}+(\Lambda^{-1}\mu)^{\beta_{1}}\leq\Lambda^{-\beta_{1}}(1+\mu^{\beta_{1}})\text{
\ in }\Omega,\end{array}$ (4.13)
$\begin{array}[]{l}(\Lambda^{-1}\hat{\phi}_{1,p_{1}})^{\alpha_{2}}+(\Lambda\hat{y}_{2})^{\beta_{2}}\leq(\Lambda^{-1}\mu)^{\alpha_{2}}+(\Lambda\frac{\hat{\phi}_{1,p_{2}}}{\hat{c}})^{\beta_{2}}\\\
\leq(\Lambda^{-1}\mu)^{\alpha_{2}}+\hat{c}^{-\beta_{2}}\leq\Lambda^{-\alpha_{2}}(\mu^{\alpha_{2}}+1)\text{
\ in }\Omega,\end{array}$

provided $\Lambda>0$ is sufficiently large. By (4.6) one has

(4.14)
$\begin{array}[]{l}-\Delta_{p_{i}}(\Lambda\hat{y}_{i})+|\Lambda\hat{y}_{i}|^{p_{i}-2}(\Lambda\hat{y}_{i})=\Lambda^{p_{i}-1}\text{\
in }\Omega,\text{ for }i=1,2.\end{array}$

Then, in view of (1.1), gathering (4.10)-(4.14) together imply

(4.15)
$\displaystyle-\Delta_{p_{1}}(\Lambda\hat{y}_{1})+|\Lambda\hat{y}_{1}|^{p_{1}-2}(\Lambda\hat{y}_{1})$
$\displaystyle\geq$
$\displaystyle\left\\{\begin{array}[]{ll}(\Lambda\hat{y}_{1})^{\alpha_{1}}+(\Lambda\hat{y}_{2})^{\beta_{1}}&\text{if
}\beta_{1}>0\\\
(\Lambda\hat{y}_{1})^{\alpha_{1}}+(\Lambda^{-1}\hat{\phi}_{1,p_{2}})^{\beta_{1}}&\text{if
}\beta_{1}<0\end{array}\right.$ $\displaystyle\geq$
$\displaystyle(\Lambda\hat{y}_{1})^{\alpha_{1}}+v^{\beta_{1}}\text{ \ in
}\Omega$

and

(4.19)
$\displaystyle-\Delta_{p_{2}}(\Lambda\hat{y}_{2})+|\Lambda\hat{y}_{2}|^{p_{2}-2}(\Lambda\hat{y}_{2})$
$\displaystyle\geq$
$\displaystyle\left\\{\begin{array}[]{ll}(\Lambda\hat{y}_{1})^{\alpha_{2}}+(\Lambda\hat{y}_{2})^{\beta_{2}}&\text{if
}\alpha_{2}>0\\\
(\Lambda^{-1}\hat{\phi}_{1,p_{1}})^{\alpha_{2}}+(\Lambda\hat{y}_{2})^{\beta_{2}}&\text{if
}\alpha_{2}<0\end{array}\right.$ $\displaystyle\geq$ $\displaystyle
u^{\alpha_{1}}+(\Lambda\hat{y}_{2})^{\beta_{2}}\text{ \ in }\Omega,$

for $\Lambda>1$ large enough. According to (4.6), (2.3) is fulfilled for
$(\overline{u},\overline{v}):=\Lambda(\hat{y}_{1},\hat{y}_{2})$ and therefore,
in view of (4.15)-(4.19), (2.5) too.

We claim that (2.4) holds for
$(\underline{u},\underline{v}):=\Lambda^{-1}(\hat{\phi}_{1,p_{1}},\hat{\phi}_{1,p_{2}})$.
From (4.2), (4.8) and (1.1), we have

(4.23)
$(\Lambda^{-1}\hat{\phi}_{1,p_{1}})^{\alpha_{1}}+v^{\beta_{1}}\geq(\Lambda^{-1}\hat{\phi}_{1,p_{1}})^{\alpha_{1}}\geq(\Lambda^{-1}||\hat{\phi}_{1,p_{1}}||_{\infty})^{\alpha_{1}}\text{\
\ in }\Omega,$ (4.24)
$u^{\alpha_{2}}+(\Lambda^{-1}\hat{\phi}_{1,p_{2}})^{\beta_{2}}\geq(\Lambda^{-1}\hat{\phi}_{1,p_{2}})^{\beta_{2}}\geq(\Lambda^{-1}||\hat{\phi}_{1,p_{2}}||_{\infty})^{\beta_{2}}\text{\
\ in }\Omega,$

as well as

(4.25)
$\begin{array}[]{l}-\Delta_{p_{i}}(\Lambda^{-1}\hat{\phi}_{1,p_{i}})+|\Lambda^{-1}\hat{\phi}_{1,p_{i}}|^{p_{i}-2}(\Lambda^{-1}\hat{\phi}_{1,p_{i}})=\hat{\lambda}_{1,p_{i}}(\Lambda^{-1}\hat{\phi}_{1,p_{i}})^{p_{i}-1}\\\
\leq\hat{\lambda}_{1,p_{i}}(\Lambda^{-1}||\hat{\phi}_{1,p_{1}}||_{\infty})^{p_{i}-1}\text{\
in }\Omega,\text{ for }i=1,2.\end{array}$

Then, gathering (4.23)-(4.25) together, since (2.3) is fulfilled (in view of
(4.2)), we deduce that (2.4) is achieved for $\Lambda>0$ large enough.
Consequently, considering that functions $\Lambda^{-1}\hat{\phi}_{1,p_{1}}$
and $\Lambda^{-1}\hat{\phi}_{1,p_{2}}$ fulfill (3.1), Theorem 3.1 ensures the
existence of a positive solution $(u,v)\in
int\mathcal{C}_{+}^{1,\tau}(\overline{\Omega})\times
int\mathcal{C}_{+}^{1,\tau}(\overline{\Omega}),$ verifying (4.9). Moreover,
$(u,v)$ is unique once $\alpha_{2},\beta_{1}\in(-1,0)$ and fulfill assumption
(3.2) in Theorem 3.1. This ends the proof. ∎

By making some specific and necessary adjustments, it is quite possible to
construct a pair of sub-supersolution in the spirit of [16]. However, it
should be noted that the sub-supersolution produced in [16] does not fulfill
the required assumptions (2.4)-(2.5), in particular the condition (2.3) which
in no case could be ignored.

Let $\Lambda>0$ be a constant so large such that

(4.26)
$\Lambda>\max_{i=1,2}\\{2(1+\frac{3}{(3^{p_{i}-1}-2^{p_{i}-1})^{\frac{1}{p_{i}-1}}}+||\hat{\phi}_{1,p_{i}}||_{\infty}+\|y_{i}\|_{\infty})\\}.$

Define

(4.27)
$\underline{u}_{\Lambda}:=\Lambda^{-1}(\Lambda-\hat{\phi}_{1,p_{1}}),\quad\underline{v}_{\Lambda}:=\Lambda^{-1}(\Lambda-\hat{\phi}_{1,p_{2}}),$
(4.28) $\overline{u}_{\Lambda}:=\Lambda(\Lambda-
y_{1}),\quad\overline{v}_{\Lambda}:=\Lambda(\Lambda-y_{2}).$

Obviously, $\underline{u}_{\Lambda}\leq\overline{u}_{\Lambda}$ and
$\underline{v}_{\Lambda}\leq\overline{v}_{\Lambda}$. Moreover, via (4.5) and
(4.2) we have

(4.29)
$\frac{\partial\underline{u}_{\Lambda}}{\partial\eta}=-\Lambda\frac{\partial\hat{\phi}_{1,p_{1}}}{\partial\eta}=0,\text{
\
}\frac{\partial\underline{v}_{\Lambda}}{\partial\eta}=-\Lambda\frac{\partial\hat{\phi}_{1,p_{1}}}{\partial\eta}=0\text{
\ on \ }\partial\Omega$

and

(4.30)
$\frac{\partial\overline{u}_{\Lambda}}{\partial\eta}=-\Lambda\frac{\partial
y_{1}}{\partial\eta}>0,\text{ \ }\
\frac{\partial\overline{v}_{\Lambda}}{\partial\eta}=-\Lambda\frac{\partial
y_{2}}{\partial\eta}>0\text{ \ \ on \ }\partial\Omega,$

showing that $(\underline{u}_{\Lambda},\underline{v}_{\Lambda})$ and
$(\overline{u}_{\Lambda},\overline{v}_{\Lambda})$ fulfill assumption (2.3).

######  Theorem 4.2.

Assume (1.1) holds. Then, for $\Lambda>0$ is big enough, problem
$(\mathrm{P})$ admits a solution $(u,v)\in
int\mathcal{C}_{+}^{1,\tau}(\overline{\Omega})\times
int\mathcal{C}_{+}^{1,\tau}(\overline{\Omega}),$ $\tau\in(0,1)$, within
$[\underline{u}_{\Lambda},\overline{u}_{\Lambda}]\times[\underline{v}_{\Lambda},\overline{u}_{\Lambda}]$,
which is unique once $\alpha_{2},\beta_{1}\in(-1,0)$ and assumption (3.2) is
fulfilled .

###### Proof.

Let $(u,v)\in W^{1,p_{1}}(\Omega)\times W^{1,p_{2}}(\Omega)$ within
$[\underline{u}_{\Lambda},\overline{u}_{\Lambda}]\times[\underline{v}_{\Lambda},\overline{v}_{\Lambda}]$.
From (4.26) observe that

(4.31) $\Lambda-||\hat{\phi}_{1,p_{i}}||_{\infty}>\frac{\Lambda}{2}\text{ and
}\Lambda-\|y_{i}\|_{\infty}>\frac{\Lambda}{2},\text{ for }i=1,2.$

Moreover, on the assumption that

(4.32) $\Lambda>\frac{6}{(3^{p_{i}-1}-2^{p_{i}-1})^{\frac{1}{p_{i}-1}}}$

we derive that

(4.33) $(\frac{\Lambda}{2})^{p_{i}-1}-(\frac{\Lambda}{3})^{p_{i}-1}>1,\text{
for }i=1,2.$

By (4.31), (4.28) and (4.33), a direct computations give

(4.34)
$\begin{array}[]{l}-\Delta_{p_{1}}\overline{u}_{\Lambda}+|\overline{u}_{\Lambda}|^{p_{1}-2}\overline{u}_{\Lambda}=\Delta_{p_{1}}(\Lambda
y_{1})+\overline{u}_{\Lambda}^{p_{1}-1}\\\
=\Lambda^{p_{1}-1}(-1+y_{1}^{p_{1}-1}+(\Lambda-y_{1})^{p_{1}-1})\\\
\geq\Lambda^{p_{1}-1}(-1+(\Lambda-\left\|y_{1}\right\|_{\infty})^{p_{1}-1})\\\
\geq\Lambda^{p_{1}-1}(-1+(\frac{\Lambda}{2})^{p_{1}-1})\\\
\geq\Lambda^{p_{1}-1}(\frac{\Lambda}{3})^{p_{1}-1}\geq(\frac{\Lambda^{2}}{3})^{p_{1}-1}\text{
in }\Omega.\end{array}$

If $\beta_{1}>0$, it follows from (4.31), (4.28) and (1.1) that

(4.35)
$\begin{array}[]{l}\overline{u}_{\Lambda}^{\alpha_{1}}+v^{\beta_{1}}\leq\overline{u}_{\Lambda}^{\alpha_{1}}+\overline{v}_{\Lambda}^{\beta_{1}}=\Lambda^{\alpha_{1}}(\Lambda-
y_{1})^{\alpha_{1}}+\Lambda^{\beta_{1}}(\Lambda-y_{2})^{\beta_{1}}\\\
\leq\Lambda^{\alpha_{1}}(\Lambda-\left\|y_{1}\right\|_{\infty})^{\alpha_{1}}+\Lambda^{2\beta_{1}}\leq
1+\Lambda^{2\beta_{1}}\text{ in }\Omega,\end{array}$

while, if $\beta_{1}<0$, we get

(4.36)
$\begin{array}[]{l}\overline{u}_{\Lambda}^{\alpha_{1}}+v^{\beta_{1}}\leq\overline{u}_{\Lambda}^{\alpha_{1}}+\underline{v}_{\Lambda}^{\beta_{1}}=\Lambda^{\alpha_{1}}(\Lambda-
y_{1})^{\alpha_{1}}+\Lambda^{-\beta_{1}}(\Lambda-\hat{\phi}_{1,p_{2}})^{\beta_{1}}\\\
\leq\Lambda^{\alpha_{1}}(\Lambda-\left\|y_{1}\right\|_{\infty})^{\alpha_{1}}+\Lambda^{-\beta_{1}}(\Lambda-||\hat{\phi}_{1,p_{2}}||_{\infty})^{\beta_{1}}\\\
\leq
1+2^{-\beta_{1}}\Lambda^{-\beta_{1}}\Lambda^{\beta_{1}}=1+2^{-\beta_{1}}\text{
in }\Omega.\end{array}$

Then, for $\Lambda>0$ large enough, (4.34)-(4.36) result in

$-\Delta_{p_{1}}\overline{u}_{\Lambda}+|\overline{u}_{\Lambda}|^{p_{1}-2}\overline{u}_{\Lambda}\geq\overline{u}_{\Lambda}^{\alpha_{1}}+v^{\beta_{1}}\text{
in }\Omega,$

for all $v\in[\underline{v}_{\Lambda},\overline{v}_{\Lambda}]$. Similar
arguments apply to the second equation in $(\mathrm{P})$ lead to

$-\Delta_{p_{2}}\overline{v}_{\Lambda}+|\overline{v}_{\Lambda}|^{p_{2}-2}\overline{v}_{\Lambda}\geq
u^{\alpha_{2}}+\overline{v}_{\Lambda}^{\beta_{1}}\text{ in }\Omega,$

for all $u\in[\underline{u}_{\Lambda},\overline{u}_{\Lambda}]$. Bearing in
mind (4.30), it follows that $(\overline{u}_{\Lambda},\overline{v}_{\Lambda})$
in (4.28) verifies assumption (2.5). Consequently,
$(\overline{u}_{\Lambda},\overline{v}_{\Lambda})$ is a supersolution of
$(\mathrm{P})$.

Next, we prove that $(\underline{u}_{\Lambda},\underline{v}_{\Lambda})$ in
(4.27) is a subsolution of $(\mathrm{P}).$ It is important to mention that the
eigenvalues in (4.2) verify
$\hat{\lambda}_{1,p_{1}},\hat{\lambda}_{1,p_{2}}>1$. Hence, using (1.1), (4.2)
and (4.31), we obtain

(4.37)
$\begin{array}[]{l}-\Delta_{p_{1}}\underline{u}_{\Lambda}+|\underline{u}_{\Lambda}|^{p_{1}-2}\underline{u}_{\Lambda}=\Delta_{p_{1}}(\Lambda^{-1}\hat{\phi}_{1,p_{1}})+\underline{u}_{\Lambda}^{p_{1}-1}\\\
=\Lambda^{-(p_{1}-1)}(1-\hat{\lambda}_{1,p_{1}})\hat{\phi}_{1,p_{1}}^{p_{1}-1}+\Lambda^{-(p_{1}-1)}(\Lambda-\hat{\phi}_{1,p_{1}})^{p_{1}-1}\\\
\leq\Lambda^{-(p_{1}-1)}(\Lambda-\hat{\phi}_{1,p_{1}})^{p_{1}-1}\leq\Lambda^{-(p_{1}-1)}(\Lambda-||\hat{\phi}_{1,p_{1}}||_{\infty})^{p_{1}-1}\\\
\leq(\frac{\Lambda-||\hat{\phi}_{1,p_{1}}||_{\infty}}{\Lambda})^{p_{1}-1}\leq
1\text{ in }\Omega\end{array}$

and

(4.38)
$\begin{array}[]{l}-\Delta_{p_{2}}\underline{v}_{\Lambda}+|\underline{v}_{\Lambda}|^{p_{2}-2}\underline{v}_{\Lambda}=\Delta_{p_{2}}(\Lambda^{-1}\hat{\phi}_{1,p_{2}})+\underline{v}_{\Lambda}^{p_{2}-1}\\\
=\Lambda^{-(p_{2}-1)}(1-\hat{\lambda}_{1,p_{2}})\hat{\phi}_{1,p_{2}}^{p_{2}-1}+\Lambda^{-(p_{2}-1)}(\Lambda-\hat{\phi}_{1,p_{2}})^{p_{2}-1}\\\
\leq\Lambda^{-(p_{2}-1)}(\Lambda-\hat{\phi}_{1,p_{2}})^{p_{2}-1}\leq\Lambda^{-(p_{2}-1)}(\Lambda-||\hat{\phi}_{1,p_{2}}||_{\infty})^{p_{2}-1}\\\
\leq(\frac{\Lambda-||\hat{\phi}_{1,p_{2}}||_{\infty}}{\Lambda})^{p_{2}-1}\leq
1\text{ in }\Omega\end{array}$

as well as

(4.39)
$\underline{u}_{\Lambda}^{\alpha_{1}}+v^{\beta_{1}}\geq\underline{u}_{\Lambda}^{\alpha_{1}}=\Lambda^{-\alpha_{1}}(\Lambda-\hat{\phi}_{1,p_{1}})^{\alpha_{1}}\geq\Lambda^{-\alpha_{1}}\Lambda^{\alpha_{1}}=1\text{
in }\Omega,$

and

(4.40)
$u_{\Lambda}^{\alpha_{2}}+\underline{v}_{\Lambda}^{\beta_{2}}\geq\underline{v}_{\Lambda}^{\beta_{2}}=\Lambda^{-\beta_{2}}(\Lambda-\hat{\phi}_{1,p_{2}})^{\beta_{2}}\geq\Lambda^{-\beta_{2}}\Lambda^{\beta_{2}}=1\text{
in }\Omega,$

for all
$(u,v)\in[\underline{u}_{\Lambda},\overline{u}_{\Lambda}]\times[\underline{v}_{\Lambda},\overline{v}_{\Lambda}]$.
Then, in view of (4.29), we deduce from (4.37)-(4.40) that
$(\underline{u}_{\Lambda},\underline{v}_{\Lambda})$ in (4.27) fulfill (2.4)
and therefore, $(\underline{u}_{\Lambda},\underline{v}_{\Lambda})$ is a
subsolution of $(\mathrm{P}).$ Using (4.31), we derive that

$\begin{array}[]{l}\underline{u}_{\Lambda}\geq\Lambda^{-1}(\Lambda-\hat{\phi}_{1,p_{1}})\geq\Lambda^{-1}(\Lambda-||\hat{\phi}_{1,p_{1}}||_{\infty})\geq\Lambda^{-1}(\frac{\Lambda}{2})=\frac{1}{2}\end{array}$

and

$\begin{array}[]{l}\underline{v}_{\Lambda}\geq\Lambda^{-1}(\Lambda-\hat{\phi}_{1,p_{2}})\geq\Lambda^{-1}(\Lambda-||\hat{\phi}_{1,p_{2}}||_{\infty})\geq\Lambda^{-1}(\frac{\Lambda}{2})=\frac{1}{2}.\end{array}$

Consequently, (3.1) is verified and thence, Theorem 3.1 guarantees the
existence of a positive solution $(u,v)\in
int\mathcal{C}_{+}^{1,\tau}(\overline{\Omega})\times
int\mathcal{C}_{+}^{1,\tau}(\overline{\Omega}),$ verifying (4.9), which is
unique once $\alpha_{2},\beta_{1}\in(-1,0)$ and assumption (3.2) is fulfilled.
This completes the proof. ∎

###  4.2. $L^{\infty}$\- bound solutions

In this part, we focus on subsolutions with zero trace condition on the
boundary $\partial\Omega$. Combined with (1.1), this type of subsolutions
bring to the forefront the singularity at the origin, generated by negative
exponents, which, until now, has been rather discreet due to the strict
positivity of the sub-supersolutions constructed. Obviously, this will create
additional difficulties in getting sub-supersolutions, especially when the
exponents are assumed all negative. Moreover, singularities combined with zero
trace subsolutions on the boundary $\partial\Omega$ constitute a significant
barrier to get regular solutions. At this point, note that we were not able to
find in the literature the Neumann counterpart of the regularity result for
Dirichlet singular problems in [17, Lemma 3.1]. This question remains open.

######  Theorem 4.3.

Assume (1.1) is fulfilled with $\alpha_{2},\beta_{1}>0$. Then, problem
$(\mathrm{P})$ admits a solution $(u,v)\in W_{b}^{1,p_{1}}(\Omega)\times
W_{b}^{1,p_{2}}(\Omega)$ such that

(4.41)
$(u,v)\in[\Lambda^{-1}\phi_{1,p_{1}},\Lambda\hat{y}_{1}]\times[\Lambda^{-1}\phi_{1,p_{2}},\Lambda\hat{y}_{2}],$

provided $\Lambda>1$ is big enough.

###### Proof.

First, putting
$(\underline{u},\underline{v}):=(\Lambda^{-1}\phi_{1,p_{1}},\Lambda^{-1}\phi_{1,p_{2}})$
note that (4.3) implies that (2.3) is satisfied. Let $(u,v)\in
W^{1,p_{1}}(\Omega)\times W^{1,p_{2}}(\Omega)$ within
$[\Lambda^{-1}\phi_{1,p_{1}},\Lambda\hat{y}_{1}]\times[\Lambda^{-1}\phi_{1,p_{2}},\Lambda\hat{y}_{2}]$.
In view of (4.1) and (4.6), we get

(4.42)
$\begin{array}[]{l}-\Delta_{p_{i}}\left(\Lambda^{-1}\phi_{1,p_{i}}\right)+|\Lambda^{-1}\phi_{1,p_{i}}|^{p_{i}-2}\left(\Lambda^{-1}\phi_{1,p_{i}}\right)=\lambda_{1,p_{i}}(\Lambda^{-1}\phi_{1,p_{i}})^{p_{i}-1}\\\
\leq\lambda_{1,p_{i}}(\Lambda^{-1}||\phi_{1,p_{1}}||_{\infty})^{p_{i}-1}\text{
\ in }\Omega\end{array}$

and

(4.43)
$-\Delta_{p_{i}}\left(\Lambda\hat{y}_{i}\right)+|\Lambda\hat{y}_{i}|^{p_{i}-2}\left(\Lambda\hat{y}_{i}\right)=\Lambda^{p_{i}-1}\text{\
\ in }\Omega,\text{ for }i=1,2.$

By (1.1) one has

(4.44)
$(\Lambda^{-1}\phi_{1,p_{1}})^{\alpha_{1}}+v^{\beta_{1}}\geq(\Lambda^{-1}\phi_{1,p_{1}})^{\alpha_{1}}\geq(\Lambda^{-1}||\phi_{1,p_{1}}||_{\infty})^{\alpha_{1}}\text{
\ in }\Omega,$ (4.45)
$u^{\alpha_{2}}+(\Lambda^{-1}\phi_{1,p_{2}})^{\beta_{2}}\geq(\Lambda^{-1}\phi_{1,p_{2}})^{\beta_{2}}\geq(\Lambda^{-1}||\phi_{1,p_{2}}||_{\infty})^{\beta_{2}}\text{
\ in }\Omega,$

while, from (4.8) and (4.4), it hold

(4.46)
$\begin{array}[]{l}(\Lambda\hat{y}_{1})^{\alpha_{1}}+v^{\beta_{1}}\leq(\Lambda\hat{y}_{1})^{\alpha_{1}}+(\Lambda\hat{y}_{2})^{\beta_{1}}\leq(\Lambda\frac{\hat{\phi}_{1,p_{1}}}{\hat{c}})^{\alpha_{1}}+(\Lambda\hat{c}\hat{\phi}_{1,p_{2}})^{\beta_{1}}\\\
\leq(\Lambda\frac{\mu}{\hat{c}})^{\alpha_{1}}+(\Lambda\hat{c}||\hat{\phi}_{1,p_{2}}||_{\infty})^{\beta_{1}}\leq\Lambda^{\beta_{1}}((\frac{\mu}{\hat{c}})^{\alpha_{1}}+(\hat{c}||\hat{\phi}_{1,p_{2}}||_{\infty})^{\beta_{1}})\text{
\ in }\Omega\end{array}$

and

(4.47)
$\begin{array}[]{l}u^{\alpha_{2}}+(\Lambda\hat{y}_{2})^{\beta_{2}}\leq(\Lambda\hat{y}_{1})^{\alpha_{2}}+(\Lambda\hat{y}_{2})^{\beta_{2}}\leq(\Lambda\hat{c}\hat{\phi}_{1,p_{1}})^{\alpha_{2}}+(\Lambda\frac{\hat{\phi}_{1,p_{2}}}{\hat{c}})^{\beta_{2}}\\\
\leq(\Lambda\hat{c}||\hat{\phi}_{1,p_{1}}||_{\infty})^{\alpha_{2}}+(\Lambda\frac{\mu}{\hat{c}})^{\beta_{2}}\leq\Lambda^{\alpha_{2}}((\hat{c}||\hat{\phi}_{1,p_{1}}||_{\infty})^{\alpha_{2}}+(\frac{\mu}{\hat{c}})^{\beta_{2}})\text{
\ in }\Omega.\end{array}$

Then, gathering (4.42)-(4.47) together, we conclude that (2.4) and (2.5) are
achieved for $\Lambda>0$ large enough. Hence,
$\Lambda^{-1}(\phi_{1,p_{1}},\phi_{1,p_{2}})$ and
$\Lambda(\hat{y}_{1},\hat{y}_{2})$ are a sub-supersolutions of problem
$(\mathrm{P})$. By (4.3), assumption (3.5) is fulfilled. Consequently, owing
to Theorem 3.2, there exists a solution $(u,v)\in
W_{b}^{1,p_{1}}(\Omega)\times W_{b}^{1,p_{2}}(\Omega)$ of system
$(\mathrm{P})$ verifying (4.41). ∎

To deal with the case when all exponents are negative in (1.1), we slightly
modify the location area of potential solutions by moving the upper limits of
the rectangle formed through sub-supersolutions. To do so, assume
$-1<\alpha_{2},\beta_{1}<0$ and let $\hat{z}_{1},\hat{z}_{2}\in
W_{b}^{1,p_{i}}(\Omega)$ be the unique solutions of Neumann Problems

(4.48)
$-\Delta_{p_{1}}\hat{z}_{1}+|\hat{z}_{1}|^{p_{1}-2}\hat{z}_{1}=d(x)^{\beta_{1}}\text{
in }\Omega,\text{ }\frac{\partial\hat{z}_{1}}{\partial\eta}=0\text{ on
}\partial\Omega,$ (4.49)
$-\Delta_{p_{2}}\hat{z}_{2}+|\hat{z}_{2}|^{p_{2}-2}\hat{z}_{2}=d(x)^{\alpha_{2}}\text{
in }\Omega,\text{ \ }\frac{\partial\hat{z}_{2}}{\partial\eta}=0\text{ on
}\partial\Omega.$

Note that, the Hardy-Sobolev type inequality (2.2) guarantees that the right-
hand sides of (4.48) and (4.49) belong to $W^{-1,p_{1}}(\Omega)^{*}$ and
$W^{-1,p_{2}}(\Omega)^{*}$, respectively. Consequently, Minty-Browder theorem
(see, e.g., [4, Theorem V.15]) implies the existence and uniqueness of
$\hat{z}_{i}\in W^{1,p_{i}}(\Omega)$ in (4.48) and (4.49), $i=1,2$. Moreover,
by weak comparison principle (see [30]), it is readily seen that there is a
constant $c_{1}>0$ such that

(4.50) $\hat{z}_{i}\geq c_{1}\hat{\phi}_{1,p_{i}}\text{ in }\Omega,\text{ for
}i=1,2.$

Next, we provide the $L^{\infty}$-bound of $\hat{z}_{i}.$

######  Lemma 4.4.

Under assumption

(4.51) $0>\beta_{1},\alpha_{2}>\frac{-1}{N},$

solutions $\hat{z}_{1}\in W^{1,p_{1}}(\Omega)$ and $\hat{z}_{2}\in
W^{1,p_{2}}(\Omega)$ of problems (4.48) and (4.49) are bounded in
$L^{\infty}(\Omega)$.

###### Proof.

We only show the $L^{\infty}$-bound of $\hat{z}_{1}$ in (4.48) because that of
$\hat{z}_{2}$ in (4.49) can be justified similarly. Inspired by [AM, Lemma 2],
for each $k\in\mathbb{N}$, set

$A_{k}=\left\\{x\in\Omega:\hat{z}_{1}(x)>k\right\\}.$

It is readily seen that $|A_{k}|\rightarrow 0,$ as well as,

(4.52) $\|d(x)^{\beta_{1}}\|_{L^{N}(A_{k})}\rightarrow 0\text{ \ as
}k\rightarrow\infty,$

because $\hat{z}_{1}\in L^{1}(\Omega)$ and $d(x)^{\beta_{1}}\in L^{N}(\Omega)$
(due to (4.51)).

Test with $(\hat{z}_{1}-k)^{+}$ in (4.48), it results in

$\begin{array}[]{l}\int_{\Omega}|\nabla\hat{z}_{1}|^{p_{1}-2}\nabla\hat{z}_{1}\nabla(\hat{z}_{1}-k)^{+}\mathrm{d}x+\int_{\Omega}\hat{z}_{1}^{p_{1}-1}(\hat{z}_{1}-k)^{+}\,\mathrm{d}x=\int_{\Omega}d(x)^{\beta_{1}}(\hat{z}_{1}-k)^{+}\mathrm{d}x.\end{array}$

Thus

(4.53)
$\begin{array}[]{l}\int_{\Omega}|\nabla\hat{z}_{1}|^{p_{1}-2}\nabla\hat{z}_{1}\nabla(\hat{z}_{1}-k)^{+}\mathrm{d}x+\int_{\Omega}\hat{z}_{1}^{p_{1}-1}(\hat{z}_{1}-k)^{+}\,\mathrm{d}x\\\
\geq\int_{\Omega}|\nabla\hat{z}_{1}|^{p_{1}-2}\nabla\hat{z}_{1}\nabla(\hat{z}_{1}-k)^{+}\mathrm{d}x+\int_{\Omega}((\hat{z}_{1}-k)^{+})^{p_{1}-1}(\hat{z}_{1}-k)^{+}\,\mathrm{d}x\\\
\geq\int_{A_{k}}|\nabla\hat{z}_{1}|^{p_{1}}\mathrm{d}x+\int_{A_{k}}(\hat{z}_{1}-k)^{p_{1}}\,\mathrm{d}x\end{array}$

while, by Hölder inequality together with Sobolev embedding
$W^{1,1}(\Omega)\hookrightarrow L^{\frac{N}{N-1}}(\Omega)$, one has

(4.54)
$\begin{array}[]{l}\int_{\Omega}d(x)^{\beta_{1}}(\hat{z}_{1}-k)^{+}\mathrm{d}x=\int_{A_{k}}d(x)^{\beta_{1}}(\hat{z}_{1}-k)^{+}\mathrm{d}x\\\
\leq\|d^{\beta_{1}}\|_{L^{N}(A_{k})}\|(\hat{z}_{1}-k)^{+}\|_{L^{\frac{N}{N-1}}(A_{k})}\leq\|d^{\beta_{1}}\|_{L^{N}(A_{k})}\|(\hat{z}_{1}-k)^{+}\|_{L^{\frac{N}{N-1}}(\Omega)}\\\
\leq
c_{p_{1}}\|d^{\beta_{1}}\|_{L^{N}(A_{k})}\|(\hat{z}_{1}-k)^{+}\|_{W^{1,1}(\Omega)}=c_{p_{1}}\|d^{\beta_{1}}\|_{L^{N}(A_{k})}\|\hat{z}_{1}-k\|_{W^{1,1}(A_{k})}.\end{array}$

Therefore, from (4.53)-(4.54), it turn out that

(4.55)
$\int_{A_{k}}|\nabla\hat{z}_{1}|^{p_{1}}\mathrm{d}x+\int_{A_{k}}(\hat{z}_{1}-k)^{p_{1}}\,\mathrm{d}x\leq
c_{p_{1}}\|d^{\beta_{1}}\|_{L^{N}(A_{k})}\int_{A_{k}}(|\nabla\hat{z}_{1}|+(\hat{z}_{1}-k))\mathrm{d}x.$

Young’s inequality yields

(4.56)
$\int_{A_{k}}(|\nabla\hat{z}_{1}|+(\hat{z}_{1}-k))\mathrm{d}x\leq\int_{A_{k}}|\nabla\hat{z}_{1}|^{p_{1}}\mathrm{d}x+\int_{A_{k}}(\hat{z}_{1}-k)^{p_{1}}\mathrm{d}x+2|A_{k}|,$

which, combined with (4.55), immediately leads to

$\begin{array}[]{l}\int_{A_{k}}|\nabla\hat{z}_{1}|^{p_{1}}\mathrm{d}x+\int_{A_{k}}(\hat{z}_{1}-k)^{p_{1}}\,\mathrm{d}x\\\
\leq
c_{p_{1}}\|d^{\beta_{1}}\|_{L^{N}(A_{k})}\left(\int_{A_{k}}|\nabla\hat{z}_{1}|^{p_{1}}\mathrm{d}x+\int_{A_{k}}(\hat{z}_{1}-k)^{p_{1}}\mathrm{d}x+2|A_{k}|\right).\end{array}$

Thereby, for $k$ large enough, the limit (4.52) implies

$\begin{array}[]{l}\int_{A_{k}}|\nabla\hat{z}_{1}|^{p_{1}}\mathrm{d}x+\int_{A_{k}}(\hat{z}_{1}-k)^{p_{1}}\,\mathrm{d}x\leq
C_{1}|A_{k}|.\end{array}$

Replacing the last inequality in (4.56), we achieve

$\begin{array}[]{l}\int_{A_{k}}(|\nabla\hat{z}_{1}|+(\hat{z}_{1}-k))\mathrm{d}x\leq
C_{2}|A_{k}|.\end{array}$

Now, once again applying Hölder inequality and the above Sobolev embedding, it
follows that

$\begin{array}[]{l}\int_{A_{k}}(\hat{z}_{1}-k)\mathrm{d}x\leq|A_{k}|^{\frac{1}{N}}\|(\hat{z}_{1}-k)\|_{L^{\frac{N}{N-1}}(A_{k})}\leq
C|A_{k}|^{\frac{1}{N}}\int_{A_{k}}(|\nabla\hat{z}_{1}|+(\hat{z}_{1}-k))\mathrm{d}x\end{array}$

and so,

$\int_{A_{k}}(\hat{z}_{1}-k)\mathrm{d}x\leq C_{3}|A_{k}|^{1+\frac{1}{N}}.$

Then, owing to [20, Lemma 5.1 , Chaper 2], we conclude that there is
$\mathrm{K}>0$, independent of $\hat{z}_{1}$, such that

$\begin{array}[]{l}\hat{z}_{1}(x)\leq\mathrm{K}\text{ a.e. in
}\Omega,\end{array}$

showing that $\hat{z}_{1}\in L^{\infty}(\Omega)$. Similar argument produces
that $\hat{z}_{2}\in L^{\infty}(\Omega)$ in (4.49). This ends the proof. ∎

######  Theorem 4.5.

Assume (1.1) is fulfilled with

(4.57) $\max\\{-1,-(p_{2}-1)\\}<\alpha_{2}<0\text{ \ and \
}\max\\{-1,-(p_{1}-1)\\}<\beta_{1}<0.$

Then, problem $(\mathrm{P})$ admits a solution $(u,v)\in
W^{1,p_{1}}(\Omega)\times W^{1,p_{2}}(\Omega)$ verifying

(4.58)
$(u,v)\in[\Lambda^{-1}\phi_{1,p_{1}},\Lambda\hat{z}_{1}]\times[\Lambda^{-1}\phi_{1,p_{2}},\Lambda\hat{z}_{2}],$

provided $\Lambda>0$ is big enough. Moreover, if (4.51) is fullfield, then
$(u,v)\in W_{b}^{1,p_{1}}(\Omega)\times W_{b}^{1,p_{2}}(\Omega)$.

###### Proof.

It is readily seen that inequalities in (4.42), (4.44) and (4.45) hold true
even for $-1<\alpha_{2},\beta_{1}<0$. Then,
$\Lambda^{-1}(\phi_{1,p_{1}},\phi_{1,p_{2}})$ is a subsolution of system
$(\mathrm{P})$ satisfying (3.5). Let us show that
$\Lambda(\hat{z}_{1},\hat{z}_{2})$ is a supersolution of $(\mathrm{P})$. By
(4.50), (4.4) and (4.3), we get

(4.59)
$\begin{array}[]{l}(\Lambda\hat{z}_{1})^{\alpha_{1}}+v^{\beta_{1}}\leq(\Lambda\hat{z}_{1})^{\alpha_{1}}+(\Lambda^{-1}\phi_{1,p_{2}})^{\beta_{1}}\\\
\leq(\Lambda
c_{1}\hat{\phi}_{1,p_{1}})^{\alpha_{1}}+(\Lambda^{-1}\phi_{1,p_{2}})^{\beta_{1}}\\\
\leq 1+(\Lambda^{-1}c_{0}d(x))^{\beta_{1}}\\\
\leq((\Lambda^{-1}d(x))^{-\beta_{1}}+c_{0}^{\beta_{1}})(\Lambda^{-1}d(x))^{\beta_{1}}\\\
\leq(1+c_{0}^{\beta_{1}})(\Lambda^{-1}d(x))^{\beta_{1}}\text{\ in
}\Omega,\end{array}$

and similarly

(4.60)
$\begin{array}[]{l}u^{\alpha_{2}}+(\Lambda\hat{z}_{2})^{\beta_{2}}\leq(\Lambda^{-1}\phi_{1,p_{1}})^{\alpha_{2}}+(\Lambda\hat{z}_{2})^{\beta_{2}}\\\
\leq(\Lambda^{-1}\phi_{1,p_{1}})^{\alpha_{2}}+(\Lambda
c_{1}\hat{\phi}_{1,p_{2}})^{\beta_{2}}\\\
\leq(\Lambda^{-1}c_{0}d(x))^{\alpha_{2}}+1\\\
\leq(c_{0}^{\alpha_{2}}+(\Lambda^{-1}d(x))^{-\alpha_{2}})(\Lambda^{-1}d(x))^{\alpha_{2}}\\\
\leq(c_{0}^{\alpha_{2}}+1)(\Lambda^{-1}d(x))^{\alpha_{2}}\text{\ in
}\Omega,\end{array}$

for all
$(u,v)\in[\Lambda^{-1}\phi_{1,p_{1}},\Lambda\hat{z}_{1}]\times[\Lambda^{-1}\phi_{1,p_{2}},\Lambda\hat{z}_{2}]$,
provided that $\Lambda>0$ is large. On the other hand, (4.48) and (4.49) imply

(4.61)
$-\Delta_{p_{1}}(\Lambda\hat{z}_{1})+|\Lambda\hat{z}_{1}|^{p_{1}-2}(\Lambda\hat{z}_{1})=\Lambda^{p_{1}-1}d(x)^{\beta_{1}}\text{
in }\Omega$

and

(4.62)
$-\Delta_{p_{2}}(\Lambda\hat{z}_{2})+|\Lambda\hat{z}_{2}|^{p_{2}-2}(\Lambda\hat{z}_{2})=\Lambda^{p_{2}-1}d(x)^{\alpha_{2}}\text{
in }\Omega.$

Then, it turns out from (4.59)-(4.62) and (4.57) that (2.5) is fulfilled for
$\Lambda>0$ large enough. Hence, $\Lambda(\hat{z}_{1},\hat{z}_{2})$ is a
supersolution of problem $(\mathrm{P})$. Consequently, owing to Theorem 3.2,
there exists a solution $(u,v)\in W^{1,p_{1}}(\Omega)\times
W^{1,p_{2}}(\Omega)$ of system $(\mathrm{P})$ verifying (4.58). Moreover,
according to Lemma 4.4 and (4.58), we infer that $(u,v)\in
W_{b}^{1,p_{1}}(\Omega)\times W_{b}^{1,p_{2}}(\Omega)$ once (4.51) is
fulfilled. This ends the proof. ∎

## References

  * [1] C.O. Alves & F.J.S.A. Correa, _On the existence of positive solution for a class of singular systems involving quasilinear operators_ , Appl. Math. Comput. 185 (2007), 727-736.
  * [2] C.O. Alves, V.D. Radulescu, _The Lane-Emden equation with variable double-phase and multiple regime_ , Proc. Amer. Math. Soc. 148 (2020), 2937-2952.
  * [3] R.D. Benguria, _The Lane-Emden equation revisited_ , Contemporary Math. 327 (2003), 11-19.
  * [4] H. Brézis, _Analyse fonctionnelle theorie et applications_ , Masson, Paris, 1983.
  * [5] S. Chandrasekhar, _An Introduction to the Study of Stellar Structure_ , Dover Publications, New York, 1957.
  * [6] S. Carl, V. K. Le & D. Motreanu, _Nonsmooth Variational Problems and Their Inequalities_ , Springer Monogr. Math., Springer, New York, 2007.
  * [7] E. Casas and L.A. Fernandez, _A Green’s formula for quasilinear elliptic operators_ , J. Math. Anal. Appl., 142 (1989), 62-73.
  * [8] M.G. Crandall, P.H. Rabinowitz, L. Tartar, _On a Dirichlet problem with a singular nonlinearity_ , Comm. Partial Diff. Eqts. 2 (1977), 193-222.
  * [9] H. Dellouche, A. Moussaoui, _Singular quasilinear elliptic systems with gradient dependence_ , Positivity 26 (2022), doi:10.1007/s11117-022-00868-3.
  * [10] H. Didi, A. Moussaoui, _Multiple positive solutions for a class of quasilinear singular elliptic systems_ , Rend. Circ. Mat. Palermo, II. Ser 69 (2020), 977-994.
  * [11] H. Didi, B. Khodja, A. Moussaoui, _Singular Quasilinear Elliptic Systems With (super-) Homogeneous Condition_ , J. Sibe. Fede. Univ. Math. Phys. 13(2) (2020), 1-9.
  * [12] A. Farina, _On the classification of solutions of the Lane-Emden equation on unbounded domains of_ $\mathbb{R}^{N}$, J. Math. Pures et App. 87 (5) (2007), 537-561.
  * [13] M. Ghergu, _Lane-Emden systems with negative exponents_ , J. Funct. Anal. 258 (2010), 3295-3318.
  * [14] J. Giacomoni, J. Hernandez, A. Moussaoui, _Quasilinear and singular systems: the cooperative case_ , Contemporary Math. 540 (2011).
  * [15] J. Giacomoni, I. Schindler & P. Takac, _Sobolev versus Hölder local minimizers and existence of multiple solutions for a singular quasilinear equation_ , A. Sc. N. Sup. Pisa (5) 6 (2007), 117-158.
  * [16] U. Guarnotta & S.A. Marano, _Infinitely many solutions to singular convective Neumann systems with arbitrarily growing reactions_ , J. Diff. Eqts. 271 (2021), 849-863.
  * [17] D. D. Hai, _On a class of singular p-Laplacian boundary value problems_ , J. Math. Anal. Appl. 383 (2011), 619-626.
  * [18] M. A. Krasnoselskii, _Topological Methods in the Theory of Nonlinear Integral Equations_ , Pergamon Press, Oxford-London-Paris, 1964. (Translated from the Russian by A. H. Armstro)
  * [19] B. Khodja & A. Moussaoui, _Positive solutions for infinite semipositone_ $/$_positone quasilinear elliptic systems with singular and superlinear terms_ , Diff. Eqts. App. 8(4) (2016), 535-546.
  * [20] O. A. Ladyzenskaja & N. N. Ural’tzeva, _Linear and Quasilinear Elliptic Equations_ , Academic Press, New York, 1968.
  * [21] J. H. Lane, _On the Theoretical Temperature of the Sun under the Hypothesis of a Gaseous Mass Maintaining Its Volume by Its Internal Heat and Depending on the Laws of Gases Known to Terrestrial Experiment_ , American J. Sci. Arts, 50 (1870), 57-74.
  * [22] G. M. Lieberman, _Boundary regularity for solutions of degenerate elliptic equations_ , Nonl. Anal. 12 (1988), 1203-1219.
  * [23] D. Motreanu, V.V. Motreanu & N.S. Papageorgiou, _Topological and Variational Methods with Applications to Nonlinear Boundary Value Problems_ , Springer, New York (2014).
  * [24] D. Motreanu, A. Moussaoui, _An existence result for a class of quasilinear singular competitive elliptic systems_ , Appl. Math. Lett.__ 38 (2014), 33-37.
  * [25] D. Motreanu, A. Moussaoui, A quasilinear singular elliptic system without cooperative structure, Acta Math. Sci.__ 34 (B) (2014), 905-916.
  * [26] D. Motreanu, A. Moussaoui, _Existence and boundedness of solutions for a singular cooperative quasilinear elliptic system_ , Complex Var. Elliptic Equ. 59 (2014), 285-296.
  * [27] A. Moussaoui, _Multiple solutions to Gierer-Meinhardt systems_ , Disc. Continuous Dyn. Systs. 43 (7) (2023), 2835-2851.
  * [28] A. Moussaoui, B. Khodja, S. Tas, _A singular Gierer-Meinhardt system of elliptic equations in_ $\mathbb{R}^{N}$, Nonl. Anal. 71 (2009), 708-716.
  * [29] B. Opic & A. Kufner, _Hardy-type inequalities_ , Pitman Res. Notes Math., Longman, Harlow, 1990.
  * [30] K. Sreenadh & S. Tiwari, _Global multiplicity results for_ $p(x)$_-Laplacian equation with nonlinear Neumann boundary condition_ ,
  * [31] M. Struwe, _Variational Methods: Applications to nonlinear partial differential equations and Hamiltonian systems_ , Springer-Verlag, Berlin, 1990.
  * [32] J.L. Vázquez, _A strong maximum principle for some quasilinear elliptic equation_ s, Appl. Math. Optimization 12 (1984), 191-202.