=

# Combined effects in mixed local-nonlocal stationary problems

Rakesh Arora11footnotemark: 1<EMAIL_ADDRESS><EMAIL_ADDRESS>Department of Mathematics and Statistics, Masaryk University, Building 08,
Kotlářská 2, Brno, 611 37, Czech Republic  Vicenţiu D. Rădulescu
<EMAIL_ADDRESS>Faculty of Applied Mathematics, AGH University of Science
and Technology, al. Mickiewicza 30, 30-059 Kraków, Poland  Department of
Mathematics, University of Craiova, Street A.I. Cuza 13, 200585 Craiova,
Romania

###### Abstract

In this work, we study an elliptic problem involving an operator of mixed
order with both local and nonlocal aspects, and in either the presence or the
absence of a singular nonlinearity. We investigate existence or non-existence
properties, power and exponential type Sobolev regularity results, and the
boundary behavior of the weak solution, in the light of the interplay between
the summability of the datum and the power exponent in singular
nonlinearities.

###### keywords:

existence results , power and exponential type Sobolev regularity , local-
nonlocal operator , boundary behavior , singular nonlinearity , Green’s
function estimates.

###### MSC:

[2010]35A01 , 35R11 , 47G20 , 35S15 , 35B65 , 35J75.

###### Contents

  1. 1 Introduction
    1. 1.1 Understanding the notion of a solution
    2. 1.2 Previous work
    3. 1.3 Description of main results
  2. 2 Statement of main results
    1. 2.1 Lebesgue weights
    2. 2.2 Singular weights
  3. 3 Preliminary lemmas
  4. 4 Uniform apriori estimates
    1. 4.1 Lebesgue weights: Sobolev regularity estimates
    2. 4.2 Green’s function estimates
    3. 4.3 Singular weights
      1. 4.3.1 Boundary behavior
      2. 4.3.2 Sobolev regularity estimates
  5. 5 Existence, summability and Sobolev regularity results
    1. 5.1 Proof of Theorem 2.1
    2. 5.2 Proof of Theorem 2.2
    3. 5.3 Proof of Theorem 2.3
    4. 5.4 Proof of Theorem 2.4
    5. 5.5 Proof of Theorem 2.5
    6. 5.6 Proof of Theorem 2.6:
    7. 5.7 Proof of Theorem 2.7
    8. 5.8 Proof of Theorem 2.8
    9. 5.9 Proof of Theorem 2.9
  6. 6 Appendix

##  1 Introduction

In this article, we study the fine properties of the weak solution to an
elliptic problem involving a mixed type operator $\mathcal{L}$, given by

$\mathcal{L}:=(-\Delta)+(-\Delta)^{s}\quad\text{for}\ 0<s<1.$ (1.1)

Here, the word “mixed” refers to the type of the operator combining both local
and nonlocal features, and to the differential order of the operator. The
operator $\mathcal{L}$ is obtained by superposition of the classical Laplacian
$(-\Delta)$ and the fractional Laplacian $(-\Delta)^{s}$, for a fixed
parameter $s\in(0,1)$, defined as

$(-\Delta)^{s}u=C(N,s)P.V.\int_{\mathbb{R}^{N}}\frac{u(x)-u(y)}{|x-y|^{N+2s}}\leavevmode\nobreak\
dy.$

The term $``P.V."$ stands for Cauchy’s principal value and $C(N,s)$ is a
normalizing constant, whose explicit expression is given by

$C(N,s)=\left(\int_{\mathbb{R}^{N}}\frac{1-cos(\xi_{1})}{|\xi|^{N+2s}}d\xi\right)^{-1}.$

The above choice of the constant $C(N,s)$ arises from the equivalent
definition of $(-\Delta)^{s}$ to view it as a pseudo-differential operator of
symbol $|\xi|^{2s}.$ Without going into the details of the appearance of such
type of nonlocal operator in real-world phenomenons and motivation behind
studying problems involving such nonlocal operators, we refer the reader to
review the famous Hitchhiker’s guide [45] and reference within.

The mixed operators of the form $\mathcal{L}$ in (1.1) appears naturally in
applied sciences, to study the role of the impact caused by a local and a non-
local change in a physical phenomenon. More precisely, it can be understood
through the following biological scenario where the population with density
$u$ can possibly alternate both short and long-range random walks (namely, a
classical random walk and a Lévy flight), and this could be driven, for
example, by a superposition between local exploration of the environment and
hunting strategies (for a through discussion see, [36] for mixed dispersal
movement strategy, [41] for nonlocal diffusion strategy, [25] for conditional
dispersal strategy). These type of operators also arises in the models
obtained from superposition of two different scaled stochastic processes. For
a detailed presentation on this, we refer the reader to [42].

Very recently a great amount of attention has been paid in studying elliptic
problems involving mixed type of operator having both local and nonlocal
behaviours. Some questions related to structural results like existence,
maximum principle and interior Sobolev and Lipschitz regularity [1, 11, 14],
symmetry results [13], Faber-Krahn type inequality [12], Neumann problems
[42], Green functions estimates [26, 27] has been answered.

The study of elliptic or integral equations involving singular terms started
in the early sixties by the works of Fulks and Maybee [31], originating from
the models of steady-state temperature distribution in an electrically
conducting medium. On the one hand, the study of such types of equations is a
challenging mathematical problem. On the other hand, they appear in a variety
of real-world models. To demonstrate an application, let us consider $\Omega$
be an electrically conducting medium in $\mathbb{R}^{3}$ where the local
voltage drop is described by the function $f$ and $u$ be the steady state
temperature distribution in the region $\Omega$. Then, if $\sigma(u)$ is the
electrical resistivity which is, in general, a function of the temperature u,
in particular $\sigma(u)=u^{\gamma}$, the rate of generation of heat at any
point $x$ in the medium is $\frac{f(x)}{u^{\gamma}}$ and the temperature
distribution in the conducting medium satisfies the local counterpart of the
equation (see below (1.2)). For interested readers, we refer to [43, 47, 30]
for applications in the pseudo-plastic fluids, Chandrasekhar equations in
radiative transfer, and in non-Newtonian fluid flows in porous media and
heterogeneous catalysts.
Motivating from the above discussion, we study the following mixed
local/nonlocal elliptic problem in the presence of weight function $f$ and
singular nonlinearities

$\mathcal{L}u=\frac{f(x)}{u^{\gamma}},\quad u>0\quad\text{ in }\Omega,$ (1.2)

subject to the homogeneous Dirichlet boundary conditions

$u=0\ \text{ in }\ \mathbb{R}^{N}\setminus\Omega,$ (1.3)

where $\Omega\subset\mathbb{R}^{N}$, $N\geq 2$, $\gamma\geq 0$. The function
$f:\Omega\to\mathbb{R}^{+}$ either belongs to the Lebesgue class of functions
$L^{r}(\Omega)$ for some $1\leq r\leq\infty$ or has a growth of negative
powers of distance function $\delta$ near the boundary i.e.
$f(x)\sim\delta^{-\zeta}(x)$ for some $\zeta\geq 0$ and $x$ lies near the
boundary $\partial\Omega$.

###  1.1 Understanding the notion of a solution

We start by understanding the meaning of a “weak” solution for the problem
(1.2)-(1.3). An elementary way to define the notion of a solution of the
problem (1.2) is given by: a function $u$ such that

  1. (i)

$u>0$ a.e. in $\Omega$,

  2. (ii)

the following weak formulation equality holds:

$\begin{split}\int_{\Omega}\nabla u\cdot\nabla\psi\leavevmode\nobreak\
dx+\frac{C(N,s)}{2}\int_{\mathbb{R}^{N}}\int_{\mathbb{R}^{N}}\frac{(u(x)-u(y))(\psi(x)-\psi(y))}{|x-y|^{N+2s}}\leavevmode\nobreak\
dx\leavevmode\nobreak\
dy=\int_{\Omega}\frac{f(x)}{u^{\gamma}}\psi\leavevmode\nobreak\
dx.\end{split}$ (1.4)

for a class of test functions $\psi\in\mathcal{T}(\Omega)$ and a function $u$
“regular” enough so that all the integrals are well defined.

The first condition $(i)$ is imposed to give a meaning to the term
$u^{-\gamma}$ known as “singular nonlinearities” and the second condition
$(ii)$ is motivated and obtained by multiplying a smooth function $\psi$ to
the equation (1.2) and using standard integration by parts formula for a
smooth function $u$. Since solutions to equations involving a fractional
Laplacian and singular nonlinearities generally are not of class $C^{2}$
therefore a solution to (1.2) has to be understood in the “weak” sense via
(1.4).

Assuming $\mathcal{T}(\Omega)$ to be a class of smooth test functions (for
e.g. $C_{c}^{\infty}(\Omega)$), a $``$natural$"$ space to look for the
solution $u$ of the problem (1.2)-(1.3), more accurately, to define the
integrals on the left hand side of (1.4), is the following:

$\mathbb{H}(\Omega):=\\{u\in H^{1}(\mathbb{R}^{N}):\ u=0\ \text{in}\
\mathbb{R}^{N}\setminus\Omega\\}.$

In the light of the boundary regularity of $\Omega$, it is well known that
that the space $\mathbb{H}(\Omega)$ can be identified with
$H_{0}^{1}(\Omega).$ Precisely, we understand the identification (see [19,
Proposition 9.18]) in the following way: the function $u\in H_{0}^{1}(\Omega)$
as a zero extension of the function $\tilde{u}:=u\cdot{\bf
1}_{\Omega}\in\mathbb{H}(\Omega)$ and

$u\in\mathbb{H}(\Omega)\Longrightarrow u|_{\Omega}\in H_{0}^{1}(\Omega).$

In view of above identification and employing the classical embedding theorem
[45, Proposition 2.2], both the integrals on left hand side of (1.4) are well
defined in $H_{0}^{1}(\Omega).$ Due to the fact that $\gamma\geq 0$, the
nonlinearity $u^{-\gamma}$ in the our problem (1.2)-(1.3) may blow up near the
boundary and this is the reason why we regard (1.2) as an equation with
“singular nonlinearities”. By taking into account the singular nature of the
nonlinearities and the regularity of the datum $f$, specifically, when the
nature of the singularity is strong i.e. $\gamma\gg 1$ (see [38]) or $f\in
L^{r}(\Omega)$ when $r$ is close to 1 (see [16]), we cannot always expect our
solution $u$ in $H_{0}^{1}(\Omega)$ and in place of that, either we have
$u^{\alpha}\in H_{0}^{1}(\Omega)$ for some $\alpha\geq 1$ or $u\in
W_{0}^{1,q}(\Omega)$ for some $q\in[1,2)$. For these reasons, as customary in
the literature, we adopt the following definition to understand the Dirichlet
datum in a generalized sense (see [21, 8]):

######  Definition 1.1.

A function $u\leq 0$ on $\partial\Omega$, if $u=0$ in
$\mathbb{R}^{N}\setminus\Omega$ and for any $\epsilon>0$, we have

$(u-\epsilon)^{+}\in H_{0}^{1}(\Omega).$

We say that $u=0$ on $\partial\Omega$, if $u\geq 0$ and $u\leq 0$ on
$\partial\Omega.$

Now, to give a meaning to the integral present on the right hand side of
(2.5), we choose a suitable class of test functions $\mathcal{T}(\Omega)$
depending upon the value of the exponent $\gamma\geq 0$ and the regularity of
the datum $f$ and the solution $u.$ Motivating from the above discussion, the
exact notion of weak solution to our main problem (1.2)-(1.3) with different
choices of a class of test functions and sufficient regularity of the solution
$u$, is detailed in Section 2 while stating the main results.

###  1.2 Previous work

One of the seminal breakthroughs in the study of singular nonlinearities was
the work of Crandall, Rabinowitz, and Tartar [29] which majorly setup this
direction of research. Afterward, a large number of publications has been
devoted to investigate a diverse spectrum of issues rotating around
local/nonlocal elliptic equations involving the singular nonlinearities (see
for example [30, 31, 47] and monographs [32, 34]). Let us recall some known
results in the literature for both local and nonlocal elliptic equations with
singular nonlinearities.

In the local case, Crandall et al. [29] studied the singular boundary value
problem (1.2)-(1.3) with Laplace operator and $f=1$. By using the classical
method of sub-supersolutions on the non-singular approximating problem, they
proved the existence and uniqueness results of the classical solution of our
original problem. In addition, by exploiting the second-order ODE techniques
and localization near the boundary, the boundary behavior of the solution is
deduced. By Stuart [48], similar results on the existence of solutions were
obtained using, this time, an approximation argument with respect to the
boundary condition. Actually, both papers [29, 48] provide results for more
general differential operators with smooth coefficients, not necessarily in
divergence form, and for non-monotone nonlinearities as well. The same model
of elliptic equations with singular nonlinearities and $f\in
C^{\alpha}(\overline{\Omega})$, was considered by Lazer and McKenna [38] in
which they simplified the proof of boundary behavior of classical by
constructing appropriate sub and super solutions. In addition to that, they
also obtained the optimal power related the the existence of finite energy
solutions. In fact, a solution $u\in H_{0}^{1}(\Omega)$ exists if and only if
$\gamma<3.$ In [49], Yijing and Zhang analyzed the threshold value $3$ when
the datum $f\in L^{1}(\Omega)$ and a positive function, and provide a
classical Lazer-Mckenna obstruction. In [16, 23], authors studied the
existence and uniqueness results when $f\in L^{r}(\Omega)$ for $r\geq 1$ and
showed how the regularity of this solution depends upon the summability of the
datum and the singular datum. In particular, Boccardo and Orsina [16] proved
the existence and regularity of distributional solution:

$\displaystyle\left\\{\begin{array}[]{ll}u\in
W_{0}^{1,\frac{Nr(1+\gamma)}{N-r(1-\gamma)}}(\Omega)&\text{ if }\ 0<\gamma<1\
\text{and}\ f\in L^{r}(\Omega)\ \text{with}\
r\in[1,(2^{\ast}/(1-\gamma))^{\prime}),\vspace{0.1cm}\\\ u\in
H_{0}^{1}(\Omega)&\text{ if }\ 0<\gamma<1\ \text{and}\ f\in L^{r}(\Omega)\
\text{with}\ r=(2^{\ast}/(1-\gamma))^{\prime},\vspace{0.1cm}\\\ u\in
H_{0}^{1}(\Omega)&\text{ if }\ \gamma=1\ \text{and}\ f\in
L^{1}(\Omega),\vspace{0.1cm}\\\ u^{\frac{1+\gamma}{2}}\in
H_{0}^{1}(\Omega)&\text{ if }\ \gamma>1\ \text{and}\ f\in
L^{1}(\Omega),\end{array}\right.$

Extending the work of [16], Arcoya and Mérida in [5], studied some particular
cases of strongly singular elliptic equations i.e. $1<\gamma<\frac{3r-1}{r+1}$
and $f\in L^{r}(\Omega)$, $r>1$ and $f$ is strictly far away from zero on
$\Omega$ and proved the power type Sobolev regularity

$u^{\alpha}\in H_{0}^{1}(\Omega)\ \text{for all}\
\alpha\in\left(\frac{(r+1)(1+\gamma)}{4r},\frac{1+\gamma}{2}\right]$

Moreover, in connection with the same problem, in [15], Boccardo and Díaz
proved the uniqueness of finite energy solution by extending the set of
admissible test functions. For similar works concerning the local elliptic or
integral equations with purely singular nonlinearities, we refer to [51, 30,
33, 47, 18] and singular nonlinearities with source terms or absorption terms,
we refer to [28, 48, 20, 44] with no intent to furnish an exhaustive list.

Turning to the nonlocal case, the singular problems have been investigated
more recently and there are few works in the literature, in particular, with
the fractional Laplacian $(-\Delta)^{s}$ and related to Lazer-Mckeena type
problem (see for instance [4, 7, 8, 10]). In [10], Barrios et al. studied the
solvability of the nonlocal problem in presence of singular nonlinearities and
weight function $f$. In particular, they proved the existence and regularity
of solution in very weak sense depending upon the regularity of the datum $f$
and the singular exponent $\gamma$,

$\displaystyle\left\\{\begin{array}[]{ll}u\in H_{0}^{s}(\Omega)&\text{ if }\
0<\gamma\leq 1\ \text{and}\ f\in L^{r}(\Omega)\ \text{with}\
r=(2_{s}^{\ast}/(1-\gamma))^{\prime},\vspace{0.1cm}\\\
u^{\frac{1+\gamma}{2}}\in H_{0}^{s}(\Omega)&\text{ if }\ \gamma>1\ \text{and}\
f\in L^{1}(\Omega),\end{array}\right.$

and very recently in [50], Youssfi and Mahmoud extended the result of [5] to
the fractional Laplacian and established the following

$u\in W_{0}^{1,\frac{Nr(1+\gamma)}{N-rs(1-\gamma)}}(\Omega)\ \text{ if }\
0<\gamma<1\ \text{and}\ f\in L^{r}(\Omega)\ \text{with}\
r\in[1,(2_{s}^{\ast}/(1-\gamma))^{\prime}),$

and for a non-negative datum $f\in L^{r}(\Omega)$, $r>1$ and $\gamma>1$

$u^{\alpha}\in H_{0}^{s}(\Omega)\ \text{for all}\
\alpha\in\left(\max\left(\frac{1}{2},\frac{sr(1+\gamma)-r+1}{2rs}\right),\frac{1+\gamma}{2}\right].$

In case of $f\sim\delta^{-\zeta}$ for some $\zeta\in[0,2s)$, Adimurthi et al.
[4] and Arora et al. [7] studied the same problem with singular nonlinearites
in case of $N>2s$ and $N=2s$ respectively, and discussed the existence and
uniqueness results of the classical solutions with respect to the singular
parameters. Moreover, using the integral representation via Green function and
maximum principle, they proved the sharp boundary behavior of the weak
solution. For further issues on nonlocal and nonlinear singular problems, the
interested reader can consult to the bibliographic references in [8, 6, 21,
23, 35].
Notations: Throughout the paper, we assume that $\Omega\subset\mathbb{R}^{N}$
($N\geq 2$) is a bounded domain with $C^{1,1}$ boundary. Set
$\delta(x):=\operatorname{dist}(x,\partial\Omega)$ and
$\mathcal{D}_{\Omega}=\operatorname{diam}(\Omega).$ For $i\in\mathbb{N}$, we
denote by $C_{i},c_{i},d_{i}$ positive constants that may vary from line to
line. If necessary, we will write $C=C(a,b)$ to emphasize the dependence of
$C$ on $a,b.$ For a number $q\in(1,\infty)$, we denote by $q^{\prime}$ the
conjugate exponent of $q$, namely $q^{\prime}=\frac{q}{q-1}$. For two
functions $f,g$, we write $f\lesssim g$ or $f\gtrsim g$ if there exists a
constant $C>0$ such that $f\leq Cg$ or $f\geq Cg$. We write $f\sim g$ if
$f\lesssim g$ and $g\gtrsim f$. For $\eta>0$, $\chi_{\Omega_{\eta}}$ denotes
the characteristic function of
$\Omega_{\eta}:=\\{x\in\Omega:\delta(x)<\eta\\}.$ Denote

$r^{\sharp}:=\left(\frac{2^{*}}{(1-\gamma)}\right)^{\prime}\ \text{for}\
0\leq\gamma<1\ \text{and}\ 2^{*}:=\frac{2N}{N-2}$
$\mathcal{P}_{r,\gamma}:=\\{(r,\gamma):\ \text{$1\leq r\leq+\infty,\
\gamma\geq 0\ \text{and}\ (r,\gamma)\neq(1,0)$}\\},$

and

$L^{r}_{c}(\Omega):=\\{f\in
L^{r}(\Omega):\operatorname{supp}(f)\Subset\Omega\\}.$

For $\gamma>0$ and $\zeta\geq 0$, we define a class of functions

$\mathcal{A}_{\zeta}(\Omega):=\left\\{f:\Omega\to\mathbb{R}^{+}\cup\\{0\\}:f\asymp\delta^{-\zeta}\right\\}$

and for $\zeta\neq 2$, denote

$\mathfrak{L}^{*}:=\frac{\gamma+\zeta-1}{2-\zeta}.$

###  1.3 Description of main results

In the present work, we derive the qualitative properties of the weak solution
to a mixed type elliptic problem for two different class of weight functions
$f$, and in both presence or absence of singular nonlinearities i.e.
$\gamma>0$ or $\gamma=0$ respectively. In this part, we give a short
description of our main results and for a detailed presentation, we refer the
reader to Section 2.
For the first class of weight function $f$, i.e. $f\in L^{r}(\Omega)$ for
$1\leq r\leq\infty$, we show:

  * •

Existence results: For this, we use classical approach of regularizing the
singular nonlinearities $u^{-\gamma}$ by $(u+\frac{1}{n})^{-\gamma}$ and
derive uniform a priori estimates for the weak solution of the regularized
problem. The crucial step here is to choose an appropriate test function in
the energy space. By taking into the account the combined interaction of the
summability of the datum $f$ and the singular exponent $\gamma$, we obtain our
existence results in two disjoint subsets of our admissible set
$\mathcal{P}_{r,\gamma}$ with different Sobolev regularity:

$u\in\left\\{\begin{array}[]{ll}W_{0}^{1,q}(\Omega)&\text{ if }\
(r,\gamma)\in\mathcal{P}_{r,\gamma}\cap\\{(r,\gamma):r\in[1,r^{\sharp}),0\leq\gamma<1\\},\vspace{0.1cm}\\\
H^{1}_{loc}(\Omega)&\text{ if }\
(r,\gamma)\in\mathcal{P}_{r,\gamma}\setminus\\{(r,\gamma):r\in[1,r^{\sharp}),0\leq\gamma<1\\},\end{array}\right.$

with $q:=\frac{Nr(1+\gamma)}{N-r(1-\gamma)}.$ Here, the notion of solution in
two disjoints subsets may differ due to different Sobolev regularity of the
solution. Moreover, we observe that there is a kind of continuity in the
summability exponents in the sense that as $r\to(r^{\sharp})^{-}$, $q\to 2$
when $0\leq\gamma<1.$ These existence results for mixed type operator also
complements the results for classical Laplacian in [16] and fractional
Laplacian in [10, 50].

  * •

Power and exponential type Sobolev regularity of the weak solution: Here the
term power type Sobolev regularity means that $u^{\alpha}\in
H_{0}^{1}(\Omega)$ for some $\alpha>0$ and the exponential type Sobolev
regularity means that there exists $\beta>0$ such that

$\exp(\beta u)-1\in H_{0}^{1}(\Omega)\ \text{when}\ 0\leq\gamma\leq 1,$

and for $\tau\geq\gamma$,

$\left(\exp(\beta u)-1\right)^{\tau}\in H_{0}^{1}(\Omega)\ \text{when}\
\gamma>1.$

For this, we prove two types of power type Sobolev regularity result,
depending upon the value of $\alpha,$ one with the help of appropriate choice
of test functions when $\alpha$ is large and the second by using the lower
boundary behavior of the approximating solution when $\alpha$ is small.
Type 1: When $1\leq r<\frac{N}{2}$, $\gamma\geq 0$, we show that

$\text{$u^{\alpha}\in H_{0}^{1}(\Omega)$ \ for any \
$\alpha\in\left[\frac{\gamma+1}{2},\frac{\mathfrak{S}_{r}+1}{2}\right]$ \
where \ $\mathfrak{S}_{r}:=\frac{N(r-1)+\gamma r(N-2)}{N-2r}$}.$

We notice that as $r\to\frac{N}{2}^{-}$, $\mathfrak{S}_{r}\to\infty$ and so it
natural to expect the exponential type Sobolev regularity when
$r=\frac{N}{2}.$ In this regard, for $r=\frac{N}{2}$, we prove the exponential
type Sobolev regularity in the sense mentioned above. The first step in
establishing such regularity results is to find an appropriate test function
in the energy space to handle the singular nonlinearities and second step is
to derive uniform a priori estimates for the approximating sequence with the
same power type Sobolev regularity and then pass to the limits. This type of
Sobolev regularity result are even new for the classical Laplacian and
fractional Laplacian singular and non-singular problems.
Type 2: In this case, to handle the singular nonlinearities, we exploit the
boundary behavior of the weak solution in deriving uniform a priori estimates.
Precisely, we show that for any $r>1$ and $\gamma>0$,

$u^{\alpha}\in H_{0}^{1}(\Omega)\ \text{for any}\
\alpha\in\left\\{\begin{array}[]{ll}\left(\frac{1}{2},\frac{\gamma+1}{2}\right]&\text{
if }\ \gamma+\frac{1}{r}<1,\vspace{0.1cm}\\\
\left(\frac{r\gamma+1}{2r},\frac{\gamma+1}{2}\right]&\text{ if }\
\gamma+\frac{1}{r}\geq 1.\end{array}\right.$

In addition to that, we also have shown that if $\alpha\leq\frac{1}{2}$, then
$u^{\alpha}\not\in H_{0}^{1}(\Omega)$ if $\gamma+\frac{1}{r}<1$, which
highlight the optimality of this result.

  * •

Continuity with respect to datum: An an application of the Type 1 power-type
Sobolev regularity results, we show that for any $1\leq r<\frac{N}{2}$ and
$\mathfrak{S}\in[\gamma,\mathfrak{S}_{r}]$ and for minimal weak solution $u$
and $v$ with respect to datum $f$ and $g$ respectively, the following
inequality holds

$\|\nabla|u-v|^{\frac{\mathfrak{S}+1}{2}}\|^{2}_{L^{2}(\Omega)}\leq
C\leavevmode\nobreak\ \|f-g\|_{L^{r}(\Omega)}^{\frac{r(N-2)}{N-2r}}.$

The above inequality also implies the continuity of the solution with respect
to given datum and comparison estimates.

For the second class of weight function $f\in\mathcal{A}_{\zeta}(\Omega)$ for
$\zeta\geq 0$, we show

  * •

Existence results: For this, we have followed the same classical method of
regularizing the singular problem, but here to prove the uniform a priori
estimates of the approximating sequence and handle singular nonlinearities, we
cannot use the same approach of exploiting Lebesgue summability of the datum
$f$, since for $\zeta\geq 1$, $f\not\in L^{r}(\Omega)$ for any $r\geq 1$. To
resolve this issue, we prove new boundary estimates of the approximating
sequence, by using the lower and upper bound estimates of the Green’s kernel
associated to the mixed operator. By considering the interplay of both
singular exponents $\gamma>0$ and $\zeta\in[0,2)$, we obtain the following
existence results with different Sobolev regularity:

$u\in\left\\{\begin{array}[]{ll}H_{0}^{1}(\Omega)&\text{ if }\
\gamma+\zeta\leq 1,\vspace{0.1cm}\\\ H^{1}_{loc}(\Omega)&\text{ if }\
\gamma+\zeta>1.\end{array}\right.$

  * •

Optimal boundary behavior: To prove this, we study the action of the Green’s
operator on the inverse of the distance function perturbed with logarithmic
nonlinearity. Using the lower and upper bound estimates of the Green’s kernel
[26] and borrowing some techniques from [2], we show

$\left\\{\begin{array}[]{ll}u\asymp\delta&\text{ if }\
\zeta+\gamma<1,\vspace{0.1cm}\\\
u\asymp\delta\ln^{\frac{1}{2-\zeta}}(\frac{\mathcal{D}_{\Omega}}{\delta})&\text{
if }\ \zeta+\gamma=1,\vspace{0.1cm}\\\
u\asymp\delta^{\frac{2-\zeta}{\gamma+1}}&\text{ if }\
\zeta+\gamma>1.\end{array}\right.$

  * •

Optimal power type Sobolev regularity and non-existence results: As an
application of the above optimal boundary behavior, we prove

$u^{\frac{\mathfrak{L}+1}{2}}\ \text{belongs to}\
H_{0}^{1}(\Omega)\quad\text{if and only
if}\quad\mathfrak{L}>\left\\{\begin{array}[]{ll}0&\text{ if }\
\zeta+\gamma\leq 1,\vspace{0.1cm}\\\ \mathfrak{L}^{*}&\text{ if }\
\zeta+\gamma>1.\end{array}\right.$

and for $\zeta\geq 2$ and $\gamma>0$, non-existence results is established.
These results for mixed type operator in this class complements the results
for fractional Laplacian in [4, 7].

The plan of the article is organized as follows: In section 2, we provide the
detailed statement of main results of this work. In section 3, we prove some
preliminary results for the approximated problem which will be used in the
rest of the work. Section 4 is devoted to deriving uniform a priori estimates,
Green’s function estimates and boundary behavior of the solution of the
approximated problem. In section 5, we provide the proof of our main results.
At the end, we give a short appendix section 6.

##  2 Statement of main results

In this section, we write a detailed statement of our main results.

###  2.1 Lebesgue weights

We start with presenting our existence results:

######  Theorem 2.1.

Assume that $f\in L^{r}(\Omega)\setminus\\{0\\}$ is non-negative and
$(r,\gamma)\in\mathcal{P}_{r,\gamma}\cap\\{(r,\gamma):r\in[1,r^{\sharp}),0\leq\gamma<1\\}$.
Then there exists a positive weak solution $u$ of the problem (1.2) in the
following sense:

  1. (i)

$u\in W_{0}^{1,q}(\Omega)$ with $q:=\frac{Nr(1+\gamma)}{N-r(1-\gamma)}.$

  2. (ii)

for every $\omega\Subset\Omega$, there exists a constant $C=C(\omega)$ such
that $0<C(\omega)\leq u.$

  3. (iii)

for every $\psi\in W_{0}^{1,q^{\prime}}(\Omega)\cap
L^{r^{\prime}}_{c}(\Omega)$

$\begin{split}\int_{\Omega}\nabla u\cdot\nabla\psi\leavevmode\nobreak\
dx+\frac{C(N,s)}{2}\int_{\mathbb{R}^{N}}\int_{\mathbb{R}^{N}}\frac{(u(x)-u(y))(\psi(x)-\psi(y))}{|x-y|^{N+2s}}\leavevmode\nobreak\
dx\leavevmode\nobreak\
dy=\int_{\Omega}\frac{f(x)}{u^{\gamma}}\psi\leavevmode\nobreak\
dx.\end{split}$ (2.1)

######  Theorem 2.2.

Assume that $f\in L^{r}(\Omega)\setminus\\{0\\}$ is non-negative and
$(r,\gamma)\in\mathcal{P}_{r,\gamma}\setminus\\{(r,\gamma):r\in[1,r^{\sharp}),0\leq\gamma<1\\}$.
Then there exists a positive weak solution $u$ of the problem (1.2) in the
following sense:

  1. (i)

$u\in H^{1}_{loc}(\Omega)$ and $u=0\ \text{on}\ \Omega\ \text{in the sense of
Definition \ref{def:boundary}}.$

  2. (ii)

for every $\omega\Subset\Omega$, there exists a constant $C=C(\omega)$ such
that $0<C(\omega)\leq u.$

  3. (iii)

for every
$\psi\in\bigcup_{\tilde{\Omega}\Subset\Omega}H^{1}_{loc}(\tilde{\Omega})\cap
L^{r^{\prime}}_{c}(\Omega)$

$\begin{split}\int_{\Omega}\nabla u\cdot\nabla\psi\leavevmode\nobreak\
dx+\frac{C(N,s)}{2}\int_{\mathbb{R}^{N}}\int_{\mathbb{R}^{N}}\frac{(u(x)-u(y))(\psi(x)-\psi(y))}{|x-y|^{N+2s}}\leavevmode\nobreak\
dx\leavevmode\nobreak\
dy=\int_{\Omega}\frac{f(x)}{u^{\gamma}}\psi\leavevmode\nobreak\
dx.\end{split}$ (2.2)

Now, we state the regularity results showing power and exponential type
Sobolev regularity depending upon the summability of the datum $f$ and the
singular exponent $\gamma\geq 0$.

######  Theorem 2.3.

Assume that $f\in L^{r}(\Omega)\setminus\\{0\\}$ is non-negative and
$(r,\gamma)\in\mathcal{P}_{r,\gamma}$. Let $u$ be a weak solution of the
problem (1.2) obtained in Theorem 2.1 and Theorem 2.2. Then,

  * (i)

(Weak case) For $r\in[1,\frac{N}{2})$ and $N>2$,

$\text{$u^{\frac{\mathfrak{S}+1}{2}}\in H_{0}^{1}(\Omega)$ \ for any \
$\mathfrak{S}\in\left[\gamma,\mathfrak{S}_{r}\right]$ \ where \
$\mathfrak{S}_{r}:=\frac{N(r-1)+\gamma r(N-2)}{N-2r}$}.$

Moreover, $u$ belongs to $L^{\sigma_{r}}(\Omega)$ with
$\sigma_{r}:=\frac{Nr(1+\gamma)}{N-2r}.$

  * (ii)

(Limit case) For $r=\frac{N}{2}$ and $N\geq 2$. Then

    1. 1.

(Weak singularity) For $0\leq\gamma\leq 1$,

$\mathfrak{H}\left(\frac{u}{2}\right)\in
H_{0}^{1}(\Omega)\quad\text{where}\quad\mathfrak{H}(t):=\exp\left(\beta
t\right)-1$

where $\beta>0$ such that

$\frac{2}{S(N)\|f\|_{\frac{N}{2}}}\geq\left\\{\begin{array}[]{ll}\beta\max\\{1,\left(\beta\gamma^{-1}\right)^{\gamma}\\}&\text{
if }\ 0<\gamma\leq 1,\vspace{0.1cm}\\\ \beta&\text{ if }\
\gamma=0.\end{array}\right.$ (2.3)

Moreover, there exist constants $C_{1},C_{2}$ depending upon
$\beta,S(N),f,|\Omega|$ such that

$\int_{\Omega}\exp\left(\frac{\beta Nu}{N-2}\right)\leq C_{1}\quad\text{when}\
N>2\quad\text{and}\quad\|u\|_{L^{\infty}(\Omega)}\leq C_{2}\quad\text{when}\
N=2.$

    2. 2.

(Strong singularity) For $\tau\geq\gamma>1$ and $\alpha>0$ such that
$\frac{\alpha}{2^{3-2\tau}}\max\\{1,\alpha^{\tau}\\}<\frac{1}{\tau
S(N)\|f\|_{\frac{N}{2}}}$

$\mathfrak{D}\left(\frac{u}{2}\right)\in
H_{0}^{1}(\Omega)\quad\text{where}\quad\mathfrak{D}(t)=\left(\exp(\alpha
t)-1\right)^{\tau}$

and there exist constants $C_{3},C_{4}$ depending upon
$\alpha,\tau,S(N),f,|\Omega|$ such that

$\int_{\Omega}\exp\left(\frac{2N\alpha\tau u}{N-2}\right)\leq
C_{3}\quad\text{when}\
N>2\quad\text{and}\quad\|u_{n}\|_{L^{\infty}(\Omega)}\leq
C_{4}\quad\text{when}\ N=2.$

In particular, $u$ belongs to $L^{s}(\Omega)$ for every $s\in[1,\infty)$ and
$N>2.$

  * (iii)

(Strong case) For $r>\frac{N}{2}$ and $N\geq 2$, $u$ belongs to
$L^{\infty}(\Omega).$

  * (iv)

(Exact Sobolev regularity) For $r\in[1,r^{\sharp})$, $0\leq\gamma<1$ and
$N>2$, $u$ belongs to $W_{0}^{1,q}(\Omega)$ with
$q:=\frac{Nr(1+\gamma)}{N-r(1-\gamma)}$ and for any $r\in[1,\infty]$,
$\gamma=1$ or $r\geq r^{\sharp}$ and $0\leq\gamma<1$, $u$ belongs to
$H_{0}^{1}(\Omega).$

######  Theorem 2.4.

Assume that $f\in L^{r}(\Omega)\setminus\\{0\\}$ is non-negative. Let $u$ be a
weak solution of the problem (1.2) obtained in Theorem 2.1 and Theorem 2.2.
Then,

$\text{$u^{\frac{\mathfrak{S}+1}{2}}$ belongs to \ $H_{0}^{1}(\Omega)$ \ for
any \ $\mathfrak{S}\in\left(0,\gamma\right]$ \ if and only if \
$\gamma+\frac{1}{r}<1$}.$

and

$\text{$u^{\frac{\mathfrak{S}+1}{2}}$ belongs to \ $H_{0}^{1}(\Omega)$ \ for
any \ $\mathfrak{S}\in\left(\gamma-1+\frac{1}{r},\gamma\right]$ \ if \
$\gamma+\frac{1}{r}\geq 1$}.$

######  Theorem 2.5.

Let $u$ and $v$ be two solutions of the problem (1.2)-(1.3) obtained in
Theorem 2.1 and Theorem 2.2 with respect to datum $f$ and $g$ in
$L^{r}(\Omega)$ for $1\leq r<\frac{N}{2}$ respectively. Then, for any
$\mathfrak{S}\in[\gamma,\mathfrak{S}_{r}]$, there exists a constant
$C=C(S(N),\mathfrak{S})$ independent of $u,v$ such that

$\begin{split}\|\nabla|u-v|^{\frac{\mathfrak{S}+1}{2}}\|^{2}_{L^{2}(\Omega)}&\leq
C\leavevmode\nobreak\
\|f-g\|_{L^{r}(\Omega)}^{\frac{r(N-2)}{N-2r}}.\end{split}$

In addition to above, we have

$\int_{\Omega}|\nabla(u-v)_{+}^{\frac{\mathfrak{S}_{r}+1}{2}}|^{2}\leavevmode\nobreak\
dx\leq\int_{\Omega}(f(x)-g(x))(u-v)_{+}^{\mathfrak{S}_{r}-\gamma}\leavevmode\nobreak\
dx.$

######  Corollary 2.1.

Under the condition of Theorem 2.5, if $f\leq g$, then $u\leq v$ a.e. in
$\Omega.$

######  Remark 2.1.

Adopting the same arguments of [22, Theorem 1.2] for the Laplacian and [21,
Theorem 1.4] for fractional Laplacian, we can obtain the uniqueness of weak
solution proved in Theorem 2.2. The crucial step is to exploit the
$H^{1}_{loc}(\Omega)$ regularity of weak solution $u$. To avoid repeating the
same steps, we skip the proof.

######  Remark 2.2.

Let $(r,\gamma)\in\mathcal{P}_{r,\gamma}$ such that $\mathfrak{S}_{r}\geq 1.$
In this case, the weak solution $u\in H_{0}^{1}(\Omega)$ therefore by
classical density arguments, the class of test function in the weak
formulation can be extended to $H_{0}^{1}(\Omega).$

###  2.2 Singular weights

In this case of weight function, we work with $N\geq 3$. First, we start with
presenting our existence results:

######  Theorem 2.6.

Let $\gamma>0$ and $\zeta\in[0,2)$ such that $\gamma+\zeta\leq 1$, and
$f\in\mathcal{A}_{\zeta}(\Omega)$. Then there exists a positive minimal
solution $u$ of the problem (1.2) in the following sense:

  1. (i)

$u\in H_{0}^{1}(\Omega).$

  2. (ii)

for every $\psi\in H_{0}^{1}(\Omega)$

$\begin{split}\int_{\Omega}\nabla u\cdot\nabla\psi\leavevmode\nobreak\
dx+\frac{C(N,s)}{2}\int_{\mathbb{R}^{N}}\int_{\mathbb{R}^{N}}\frac{(u(x)-u(y))(\psi(x)-\psi(y))}{|x-y|^{N+2s}}\leavevmode\nobreak\
dx\leavevmode\nobreak\
dy=\int_{\Omega}\frac{f(x)}{u^{\gamma}}\psi\leavevmode\nobreak\
dx.\end{split}$ (2.4)

######  Theorem 2.7.

Let $\gamma>0$ and $\zeta\in[0,2)$ such that $\gamma+\zeta>1$ and
$f\in\mathcal{A}_{\zeta}(\Omega)$. Then there exists a positive minimal
solution $u$ of the problem (1.2) in the following sense:

  1. (i)

$u\in H^{1}_{loc}(\Omega)$ and $u=0\ \text{on}\ \Omega\ \text{in the sense of
Definition \ref{def:boundary}}.$

  2. (ii)

for every $\omega\Subset\Omega$, there exists a constant $C=C(\omega)$ such
that $0<C(\omega)\leq u.$

  3. (iii)

for every $\psi\in H^{1}_{0}(\Omega)$ in case of $\mathfrak{L}^{*}\leq 1$ and
$\psi\in\bigcup_{\tilde{\Omega}\Subset\Omega}H^{1}_{loc}(\tilde{\Omega})$ with
$\operatorname{supp}(\psi)\Subset\Omega$ in case of $\mathfrak{L}^{*}>1$,

$\begin{split}\int_{\Omega}\nabla u\cdot\nabla\psi\leavevmode\nobreak\
dx+\frac{C(N,s)}{2}\int_{\mathbb{R}^{N}}\int_{\mathbb{R}^{N}}\frac{(u(x)-u(y))(\psi(x)-\psi(y))}{|x-y|^{N+2s}}\leavevmode\nobreak\
dx\leavevmode\nobreak\
dy=\int_{\Omega}\frac{f(x)}{u^{\gamma}}\psi\leavevmode\nobreak\
dx.\end{split}$ (2.5)

Now, we state the regularity results displaying optimal power type Sobolev
regularity results and optimal boundary behavior of minimal solution depending
upon the singular exponents $\zeta$ and $\gamma$.

######  Theorem 2.8.

Let $\gamma>0$, $\zeta\in[0,2)$ and $u$ be a weak solution of the problem
(1.2) obtained in Theorem 2.6 and Theorem 2.7. Then,

$u\ \text{belongs to}\ \mathcal{B}_{\gamma,\zeta}(\Omega)\ \text{where}\
\begin{aligned}
\mathcal{B}_{\gamma,\zeta}(\Omega)=\left\\{\begin{array}[]{ll}u:u\asymp\delta&\text{
if }\ \zeta+\gamma<1,\vspace{0.1cm}\\\
u:u\asymp\delta\ln^{\frac{1}{2-\zeta}}(\frac{\mathcal{D}_{\Omega}}{\delta})&\text{
if }\ \zeta+\gamma=1,\vspace{0.1cm}\\\
u:u\asymp\delta^{\frac{2-\zeta}{\gamma+1}}&\text{ if }\
\zeta+\gamma>1,\end{array}\right.\end{aligned}$

and

$u^{\frac{\mathfrak{L}+1}{2}}\ \text{belongs to}\
H_{0}^{1}(\Omega)\quad\text{if and only
if}\quad\mathfrak{L}>\left\\{\begin{array}[]{ll}0&\text{ if }\
\zeta+\gamma\leq 1,\vspace{0.1cm}\\\ \mathfrak{L}^{*}&\text{ if }\
\zeta+\gamma>1.\end{array}\right.$

As an applications of above theorem, we prove the following non-existence
result:

######  Theorem 2.9.

Let $\zeta\geq 2.$ Then, there does not exists a weak solution of the problem
(1.2) in the sense of Theorem 2.6 and Theorem 2.7.

######  Remark 2.3.

Repeating the same proof [8, Theorem 1.1], we can obtain the uniqueness of
weak solution proved in Theorem 2.6 and Theorem 2.7 when $\gamma>0$ and
$\beta\in[0,\frac{3}{2})$.

######  Remark 2.4.

Let $\zeta\in[0,2),\gamma>0$ such that $\gamma+\zeta\leq 1$, and
$\zeta+\gamma>1$ and $\mathfrak{L}^{*}<1.$ In this case, the weak solution
$u\in H_{0}^{1}(\Omega)$ therefore by classical density arguments, the class
of test function in the weak formulation can be extended to
$H_{0}^{1}(\Omega).$ The condition $\mathfrak{L}^{*}<1$ can be rewritten as
$\gamma+2\zeta<3$ which further reduces to the classical Lazer-Mckeena
obstruction in case of $\zeta=0.$

##  3 Preliminary lemmas

######  Lemma 3.1.

Let $h\in L^{\infty}(\Omega)$, $h\geq 0$ and $h\not\equiv 0$. Then the problem

$\left\\{\begin{aligned} -\Delta u+(-\Delta)^{s}u&=h,\quad u>0\quad&&\text{ in
}\Omega,\\\ u&=0&&\text{ in }\
\mathbb{R}^{N}\setminus\Omega,\end{aligned}\right.$ ($S$)

admits a unique positive weak solution $u$. Moreover, $u\in
L^{\infty}(\Omega)\cap C^{1,\beta}(\overline{\Omega})$ for every
$\beta\in(0,1).$

###### Proof.

To prove the existence result, we use the classical minimization arguments
from the Calculus of variations. For the sake of completeness, we provide a
brief sketch of the proof. For any $h\in L^{\infty}(\Omega)$, $h\geq 0$ and
$h\not\equiv 0$, we define the energy functional
$\mathcal{I}_{h}:H_{0}^{1}(\Omega)\to\mathbb{R}$ such that

$\mathcal{I}_{h}(u):=\frac{1}{2}\int_{\Omega}|\nabla
u|^{2}\leavevmode\nobreak\
dx+\frac{C(N,s)}{4}\int_{\mathbb{R}^{N}}\int_{\mathbb{R}^{N}}\frac{|u(x)-u(y)|^{2}}{|x-y|^{N+2s}}\leavevmode\nobreak\
dx\leavevmode\nobreak\ dy-\int_{\Omega}hu\leavevmode\nobreak\ dx.$

First, we find the minimizer of the above energy functional and then look for
the solution of the problem $(S)$ as a critical point of $\mathcal{I}_{h}.$ In
fact, using [45, Proposition 2.2] and Sobolev’s embedding, we observe that
$\mathcal{I}_{h}$ is well-defined and

$\mathcal{I}_{h}(u)\geq\frac{1}{2}\|\nabla
u\|_{2}^{2}-|\Omega|^{\frac{1}{2}}\|h\|_{\infty}\|u\|_{2}\geq\|\nabla
u\|_{2}\left(\frac{1}{2}\|\nabla
u\|_{2}-S(N)|\Omega|^{\frac{1}{2}}\|h\|_{\infty}\right)\to\infty\ \text{as}\
\|\nabla u\|_{2}\to\infty$

where $S(N)$ is the Sobolev constant. This implies the energy functional
$\mathcal{I}_{h}$ is coercive. Moreover, $\mathcal{I}_{h}$ is a $C^{1}$ and
convex energy functional. Thus, $\mathcal{I}_{h}$ is weakly lower semi-
continuous. Combining all the above properties of $\mathcal{I}_{h}$, there
exists a minimizer $u\in H_{0}^{1}(\Omega)$ and which is also a critical point
of $\mathcal{I}_{h}$ i.e.

$\begin{split}\int_{\Omega}\nabla u\cdot\nabla\psi\leavevmode\nobreak\
dx+\frac{C(N,s)}{2}\int_{\mathbb{R}^{N}}\int_{\mathbb{R}^{N}}\frac{(u(x)-u(y))(\psi(x)-\psi(y))}{|x-y|^{N+2s}}\leavevmode\nobreak\
dx\leavevmode\nobreak\ dy=\int_{\Omega}h\psi\leavevmode\nobreak\
dx\quad\forall\psi\in H_{0}^{1}(\Omega).\end{split}$ (3.1)

By taking $\psi=u^{-}:=\min\\{u,0\\}\in H_{0}^{1}(\Omega)$ as a test function
in (3.1) and using the fact that $h\geq 0$ and

$(u(x)-u(y))(u^{-}(x)-u^{-}(y))\geq 0\quad\text{for all}\
(x,y)\in\mathbb{R}^{N}\times\mathbb{R}^{N},$

we obtain

$\int_{\Omega}|\nabla u^{-}|^{2}\leavevmode\nobreak\ dx\leq 0.$

Hence, $u\geq 0$ a.e. in $\Omega.$ Now, we show that the problem $(S)$ has a
unique solution. Let $u_{1},u_{2}\in H_{0}^{1}(\Omega)$ be two solution of the
problem $(S)$. Therefore, for all $\psi\in H_{0}^{1}(\Omega)$, we have

$\begin{split}\int_{\Omega}\nabla u_{1}\cdot\nabla\psi\leavevmode\nobreak\
dx+\frac{C(N,s)}{2}\int_{\mathbb{R}^{N}}\int_{\mathbb{R}^{N}}\frac{(u_{1}(x)-u_{1}(y))(\psi(x)-\psi(y))}{|x-y|^{N+2s}}\leavevmode\nobreak\
dx\leavevmode\nobreak\ dy=\int_{\Omega}h\psi\leavevmode\nobreak\
dx,\end{split}$ (3.2) $\begin{split}\int_{\Omega}\nabla
u_{2}\cdot\nabla\psi\leavevmode\nobreak\
dx+\frac{C(N,s)}{2}\int_{\mathbb{R}^{N}}\int_{\mathbb{R}^{N}}\frac{(u_{2}(x)-u_{2}(y))(\psi(x)-\psi(y))}{|x-y|^{N+2s}}\leavevmode\nobreak\
dx\leavevmode\nobreak\ dy=\int_{\Omega}h\psi\leavevmode\nobreak\
dx.\end{split}$ (3.3)

By subtracting the above two equations and inserting $\psi=u_{1}-u_{2}$, we
obtain

$\begin{split}\int_{\Omega}|\nabla(u_{1}-u_{2})|^{2}\leavevmode\nobreak\
dx+\frac{C(N,s)}{2}\int_{\mathbb{R}^{N}}\int_{\mathbb{R}^{N}}\frac{|(u_{1}-u_{2})(x)-(u_{1}-u_{2})(y)|^{2}}{|x-y|^{N+2s}}\leavevmode\nobreak\
dx\leavevmode\nobreak\ dy=0.\end{split}$

which further gives $u_{1}=u_{2}$ a.e. in $\Omega.$ The boundedness and
regularity of solution follows by employing the classical method of
Stampacchia see for e.g [11, Theorem 4.7] and using [12, Theorem 2.7]. Since
$h\not\equiv 0$, we have $u\not\equiv 0$. Finally, the strong maximum
principle in [14, Theorem 3.1], implies $u>0$ in $\Omega.$ ∎

Depending upon the class of weight function $f$, we consider a sequence of
increasing function $f_{n}$ such that $f_{n}\to f$ $a.e.$ in $\Omega$. In the
first case, when $f\in L^{r}(\Omega)$, we consider
$f_{n}(x):=\min\\{f(x),n\\}$ and in second case, when
$f\in\mathcal{A}_{\zeta}$, we consider

$f_{n}(x):=\left\\{\begin{array}[]{ll}\left(f^{\frac{-1}{\zeta}}(x)+\left(\frac{1}{n}\right)^{\frac{\gamma+1}{2-\zeta}}\right)^{-\zeta}&\text{
if }\ x\in\Omega,\vspace{0.1cm}\\\ 0&\text{else},\end{array}\right.$

and there exist positive constants $\mathcal{G}_{1},\mathcal{G}_{2}>0$ such
that, for any $x\in\Omega$

$\frac{\mathcal{G}_{1}}{\left(d(x)+\left(\frac{1}{n}\right)^{\frac{\gamma+1}{2-\zeta}}\right)^{\zeta}}\leq
f_{n}(x)\leq\frac{\mathcal{G}_{2}}{\left(d(x)+\left(\frac{1}{n}\right)^{\frac{\gamma+1}{2-\zeta}}\right)^{\zeta}}.$
(3.4)

By considering the above choice of $f_{n}$, we study the following
approximated singular problem

$\left\\{\begin{aligned} -\Delta
u+(-\Delta)^{s}u&=\frac{f_{n}}{(u+\frac{1}{n})^{\gamma}},\quad
u>0\quad&&\text{ in }\Omega,\\\ u&=0&&\text{ in }\
\mathbb{R}^{N}\setminus\Omega.\end{aligned}\right.$ ($P_{n}$)

where $\\{f_{n}\\}_{n\in\mathbb{N}}$ is a bounded increasing sequence such
that $f_{n}\to f$ in $L^{r}(\Omega)$ for $r\in[1,\infty)$ and $f_{n}=f$ when
$r=+\infty.$

######  Lemma 3.2.

For any $n\in\mathbb{N}$ and $\gamma\geq 0$, there exists a unique non-
negative weak solution $u_{n}\in H_{0}^{1}(\Omega)$ of the problem $(P_{n})$
in the sense that

$\begin{split}\int_{\Omega}\nabla u_{n}\cdot\nabla\psi\leavevmode\nobreak\
dx+\frac{C(N,s)}{2}\int_{\mathbb{R}^{N}}\int_{\mathbb{R}^{N}}&\frac{(u_{n}(x)-u_{n}(y))(\psi(x)-\psi(y))}{|x-y|^{N+2s}}\leavevmode\nobreak\
dx\leavevmode\nobreak\ dy\\\
&=\int_{\Omega}\frac{f_{n}(x)}{(u_{n}+\frac{1}{n})^{\gamma}}\psi\leavevmode\nobreak\
dx\quad\forall\psi\in H_{0}^{1}(\Omega).\end{split}$ (3.5)

Moreover,

  * (i)

The solution $u_{n}\in C^{1,\beta}(\overline{\Omega})\cap L^{\infty}(\Omega)$
for every $\beta\in(0,1)$ and $u_{n}>0$ in $\Omega.$ Moreover, when $f_{n}\in
C^{\alpha}(\overline{\Omega})$, then $u_{n}\in
C^{2,\delta}(\overline{\Omega})\cap C(\overline{\Omega})$ for some
$\delta\in(0,1).$

  * (ii)

The sequence $\\{u_{n}\\}_{n\in\mathbb{N}}$ is monotonically increasing in the
sense that $u_{n+1}\geq u_{n}$ for all $n\in\mathbb{N}.$

  * (iii)

For every compact set $K\Subset\Omega$ and $n\in\mathbb{N}$, there exists a
constant $C$ depending upon $K$ and independent of $n$ such that $u_{n}\geq
C>0.$

###### Proof.

Given $g\in L^{2}(\Omega)$ and $n\in\mathbb{N}$ fixed, set

$h:=\frac{f_{n}}{\left(g^{+}+\frac{1}{n}\right)^{\gamma}}.$

Then, in view of Lemma 3.1, there exists a unique positive bounded solution
$w\in H_{0}^{1}(\Omega)$ for the problem $(S)$ (see statement of Lemma 3.1)
with $h$ defined above. Therefore, we define a operator $T:L^{2}(\Omega)\to
L^{2}(\Omega)$ such that $T(g)=w$ where $w$ satisfies

$\begin{split}\int_{\Omega}\nabla w\cdot\nabla\psi\leavevmode\nobreak\
dx+\frac{C(N,s)}{2}\int_{\mathbb{R}^{N}}\int_{\mathbb{R}^{N}}\frac{(w(x)-w(y))(\psi(x)-\psi(y))}{|x-y|^{N+2s}}\leavevmode\nobreak\
dx\leavevmode\nobreak\
dy=\int_{\Omega}\frac{f_{n}(x)}{(g^{+}+\frac{1}{n})^{\gamma}}\psi\leavevmode\nobreak\
dx\quad\forall\psi\in H_{0}^{1}(\Omega).\end{split}$ (3.6)

Using $w$ as a test function in (3.6) and using Sobolev’s embedding, we obtain

$\int_{\Omega}|\nabla w|^{2}\leavevmode\nobreak\
dx\leq\int_{\Omega}\frac{f_{n}(x)w}{(g^{+}+\frac{1}{n})^{\gamma}}\leavevmode\nobreak\
dx\leq R\|\nabla w\|_{2}\Longrightarrow\|\nabla w\|_{2}\leq R$ (3.7)

where $R:=S(N)\|f_{n}\|_{L^{\infty}(\Omega)}n^{\gamma}|\Omega|^{\frac{1}{2}}$.
This implies that the ball $B(0,R)$ of radius $R$ in $H_{0}^{1}(\Omega)$ is
invariant under the action of the map $T$. Now, we prove the continuity and
compactness of the map $T:H_{0}^{1}(\Omega)\to H_{0}^{1}(\Omega)$ in order to
apply the Schauder’s fixed point theorem. For continuity, we claim that
$w_{k}\to w$ in $H_{0}^{1}(\Omega)$ when $h_{k}\to h$ in $H_{0}^{1}(\Omega)$
with $w_{k}=T(h_{k})$ and $w=T(h).$ Considering the corresponding the sequence
$\\{w\\}_{k\in\mathbb{N}}$ and choosing $w_{k}-w$ as a test function, we get

$\begin{split}\int_{\Omega}|\nabla w_{k}-w|^{2}\leavevmode\nobreak\
dx&=\int_{\Omega}\left(\frac{f_{n}(x)}{(h_{k}^{+}+\frac{1}{n})^{\gamma}}-\frac{f_{n}(x)}{(h^{+}+\frac{1}{n})^{\gamma}}\right)(w_{k}-w)\leavevmode\nobreak\
dx\\\ &\leq
S(N)\left(\int_{\Omega}\left(\frac{f_{n}(x)}{(h_{k}^{+}+\frac{1}{n})^{\gamma}}-\frac{f_{n}(x)}{(h^{+}+\frac{1}{n})^{\gamma}}\right)^{\frac{2N}{N+2}}\leavevmode\nobreak\
dx\right)^{\frac{N+2}{2N}}\|\nabla w_{k}-\nabla w\|_{2}.\end{split}$ (3.8)

Notice that, the integrand in the first term is dominated by
$2\|f_{n}\|_{\infty}n^{\gamma}$ and converge to $0$ a.e. in $\Omega$, then by
applying dominated convergence theorem we obtain our claim. For compactness,
let $h_{k}$ be a bounded sequence in $H_{0}^{1}(\Omega)$ and for
$w_{k}=T(h_{k})$, we claim that $w_{k}\to w$ in $H_{0}^{1}(\Omega)$ up to a
subsequence for some $w\in H_{0}^{1}(\Omega).$ In light of (3.7), both $w_{k}$
and $h_{k}$ are bounded in $H_{0}^{1}(\Omega)$, then up to a subsequence we
have

$w_{k}\rightharpoonup w\ \text{in}\ H_{0}^{1}(\Omega),\ w_{k}\to w\ \text{
in}\ L^{p}(\Omega)\ \text{ for any}\ 1\leq p<\frac{2N}{N-2}\ \text{and}\
h_{k}\to h\ \text{a.e. in}\ \Omega.$

We also know that, $w_{k}$ satisfies: for any $\psi\in H_{0}^{1}(\Omega)$

$\begin{split}\int_{\Omega}\nabla w_{k}\cdot\nabla\psi\leavevmode\nobreak\
dx+\frac{C(N,s)}{2}\int_{\mathbb{R}^{N}}\int_{\mathbb{R}^{N}}\frac{(w_{k}(x)-w_{k}(y))(\psi(x)-\psi(y))}{|x-y|^{N+2s}}\leavevmode\nobreak\
dx\leavevmode\nobreak\
dy=\int_{\Omega}\frac{f_{n}(x)}{(h_{k}^{+}+\frac{1}{n})^{\gamma}}\psi\leavevmode\nobreak\
dx.\end{split}$ (3.9)

Now to pass limits $k\to\infty$, we observe that

$\frac{(w_{k}(x)-w_{k}(y))(\psi(x)-\psi(y))}{|x-y|^{\frac{N+2s}{2}}}\ \text{in
uniformly bounded in}\ L^{2}(\mathbb{R}^{N}\times\mathbb{R}^{N})$

because of Sobolev’s embedding and by the pointwise convergence of $w_{k}$ to
$w$, we have

$\frac{(w_{k}(x)-w_{k}(y))(\psi(x)-\psi(y))}{|x-y|^{\frac{N+2s}{2}}}\to\frac{(w(x)-w(y))(\psi(x)-\psi(y))}{|x-y|^{\frac{N+2s}{2}}}\
\text{a.e. in}\ \Omega.$

Then, since

$\frac{(\psi(x)-\psi(y))(\psi(x)-\psi(y))}{|x-y|^{\frac{N+2s}{2}}}\in
L^{2}(\mathbb{R}^{N}\times\mathbb{R}^{N}),$

it is easy to see the passage of limit in (3.9) to

$\int_{\Omega}\nabla w\cdot\nabla\psi\leavevmode\nobreak\
dx+\frac{C(N,s)}{2}\int_{\mathbb{R}^{N}}\int_{\mathbb{R}^{N}}\frac{(w(x)-w(y))(\psi(x)-\psi(y))}{|x-y|^{N+2s}}\leavevmode\nobreak\
dx\leavevmode\nobreak\
dy=\int_{\Omega}\frac{f_{n}(x)}{(h^{+}+\frac{1}{n})^{\gamma}}\psi\leavevmode\nobreak\
dx.$

Now, by arguing as in (3.8), we get $w_{k}\to w$ in $H_{0}^{1}(\Omega)$ and
Schauder’s fixed point theorem implies the existence of a fixed points $u_{n}$
such that $T(u_{n})=u_{n}$ for all $n\in\mathbb{N}$ i.e.

$\int_{\Omega}\nabla u_{n}\cdot\nabla\psi\leavevmode\nobreak\
dx+\frac{C(N,s)}{2}\int_{\mathbb{R}^{N}}\int_{\mathbb{R}^{N}}\frac{(u_{n}(x)-u_{n}(y))(\psi(x)-\psi(y))}{|x-y|^{N+2s}}\leavevmode\nobreak\
dx\leavevmode\nobreak\
dy=\int_{\Omega}\frac{f_{n}(x)}{(u_{n}^{+}+\frac{1}{n})^{\gamma}}\psi\leavevmode\nobreak\
dx.$

The boundedness and regularity of solution follows by employing the classical
method of Stampacchia see for e.g [11, Theorem 4.7] and using [12, Theorem
2.7] i.e. $u_{n}\in L^{\infty}(\Omega)\cap C^{1,\beta}(\overline{\Omega})$.
Since $\frac{f_{n}}{\left(u^{+}_{n}+\frac{1}{n}\right)^{\gamma}}\not\equiv 0$,
we have $u_{n}\not\equiv 0$. Finally, the strong maximum principle in [14,
Theorem 3.1], implies $u_{n}>0$ in $\Omega$ and $u_{1}\geq C(K)>0$ for every
$K$ compact subset of $\Omega.$ To prove the monotonicity, let $u_{n}$ and
$u_{n+1}$ are positive solutions of the problem $(P_{n})$ and $(P_{n+1})$
respectively i.e. for any $\psi\in H_{0}^{1}(\Omega)$, we have

$\int_{\Omega}\nabla u_{n}\cdot\nabla\psi\leavevmode\nobreak\
dx+\frac{C(N,s)}{2}\int_{\mathbb{R}^{N}}\int_{\mathbb{R}^{N}}\frac{(u_{n}(x)-u_{n}(y))(\psi(x)-\psi(y))}{|x-y|^{N+2s}}\leavevmode\nobreak\
dx\leavevmode\nobreak\
dy=\int_{\Omega}\frac{f_{n}(x)}{(u_{n}+\frac{1}{n})^{\gamma}}\psi\leavevmode\nobreak\
dx$

and

$\int_{\Omega}\nabla u_{n+1}\cdot\nabla\psi\leavevmode\nobreak\
dx+\frac{C(N,s)}{2}\int_{\mathbb{R}^{N}}\int_{\mathbb{R}^{N}}\frac{(u_{n+1}(x)-u_{n+1}(y))(\psi(x)-\psi(y))}{|x-y|^{N+2s}}\leavevmode\nobreak\
dx\leavevmode\nobreak\
dy=\int_{\Omega}\frac{f_{n+1}(x)}{(u_{n+1}+\frac{1}{n+1})^{\gamma}}\psi\leavevmode\nobreak\
dx.$

Subtracting the above equalities by taking the test function
$\psi=(u_{n}-u_{n+1})^{+}$ and using the following inequality ([40, Lemma 9]):
$\ \text{for a.e.}\ (x,y)\in\mathbb{R}^{N}\times\mathbb{R}^{N}$

$\left(((u_{n}-u_{n+1})(x)-(u_{n}-u_{n+1})(y))((u_{n}-u_{n+1})^{+}(x)-(u_{n}-u_{n+1})^{+}(y))\right)\geq
0,$

we get

$\begin{split}\int_{u_{n}\geq
u_{n+1}}|\nabla(u_{n}-u_{n+1})|^{2}\leavevmode\nobreak\
dx&\leq\int_{\Omega}\left(\frac{f_{n}(x)}{(u_{n}+\frac{1}{n})^{\gamma}}-\frac{f_{n+1}(x)}{\left(u_{n+1}+\frac{1}{n+1}\right)^{\gamma}}\right)(u_{n}-u_{n+1})^{+}\leavevmode\nobreak\
dx\\\
&\leq\int_{\Omega}f_{n+1}(x)\frac{(u_{n}+\frac{1}{n})^{\gamma}-(u_{n+1}+\frac{1}{n+1})^{\gamma}}{\left((u_{n}+\frac{1}{n})(u_{n+1}+\frac{1}{n+1})\right)^{\gamma}}(u_{n}-u_{n+1})^{+}\leavevmode\nobreak\
dx\leq 0\end{split}$

which in turn implies that $u_{n}\leq u_{n+1}$ in $\Omega$ and $u_{n}\geq
C(K)$ for every $n\in\mathbb{N}$ and $K\Subset\Omega.$ ∎

Now, we derive the lower boundary behavior of the approximated sequence
$u_{n}$ with the help of integral representation of the solution via Green’s
operator.

######  Theorem 3.10.

Let $u_{n}$ be the weak solution of the problem $(P_{n})$ obtained in Lemma
3.2, then there exist a constant $C_{0}>0$ independent of $n$ and $x\in\Omega$
such that $C_{0}\delta(x)\leq u_{n}(x)$.

###### Proof.

Let $u_{n}$ be the sequence of weak solution of the problem $(P_{n})$ and
$G(x,y)$ the Green function associated to the mixed operator with homogeneous
Dirichlet boundary conditions in $\Omega.$ Then, using the integral
representation of the solution, we have

$u_{n}(x):=\mathbb{G}^{\Omega}\left[\frac{f_{n}(y)}{(u_{n}+\frac{1}{n})^{\gamma}}\right](x)=\int_{\Omega}\frac{G(x,y)f_{n}(y)}{(u_{n}+\frac{1}{n})^{\gamma}}\leavevmode\nobreak\
dy\quad\text{for}\ n\in\mathbb{N}.$ (3.10)

From Lemma 3.2, we know that $\\{u_{n}\\}_{n\in\mathbb{N}}$ is an increasing
sequence such that $u_{n}\in L^{\infty}(\Omega)$. Using this fact for any
$n\in\mathbb{N}$ and $x\in\Omega$, we deduce that

$\begin{split}\delta(x)\lesssim\varphi_{1}(x)&=\lambda_{1}\mathbb{G}^{\Omega}[\varphi_{1}](x)=\lambda_{1}\mathbb{G}^{\Omega}\left[\varphi_{1}\frac{(\|u_{1}\|_{L^{\infty}(\Omega)}+1)^{\gamma}}{(\|u_{1}\|_{L^{\infty}(\Omega)}+1)^{\gamma}}\right](x)\\\
&\leq\lambda_{1}\|\varphi_{1}\|_{L^{\infty}(\Omega)}(\|u_{1}\|_{L^{\infty}(\Omega)}+1)^{\gamma}(\mathcal{D}_{\Omega}+1)^{\zeta}\mathbb{G}^{\Omega}\left[\frac{f_{1}(x)}{(\|u_{1}\|_{L^{\infty}(\Omega)}+1)^{\gamma}}\right](x)\\\
&\leq
C\mathbb{G}^{\Omega}\left[\frac{f_{1}(x)}{(u_{1}+1)^{\gamma}}\right](x)\lesssim
u_{1}(x)\leq u_{n}(x).\end{split}$ (3.11)

∎

##  4 Uniform apriori estimates

###  4.1 Lebesgue weights: Sobolev regularity estimates

######  Lemma 4.1.

Assume that $f\in L^{r}(\Omega)\setminus\\{0\\}$ is non-negative and
$(r,\gamma)\in\mathcal{P}_{r,\gamma}$. Let $u_{n}$ be the weak solution of the
problem $(P_{n})$.

  * (i)

For $r\in[1,\frac{N}{2})$ and $N>2$,

$\text{$u_{n}^{\frac{\mathfrak{S}_{r}+1}{2}}$ is uniformly bounded in
$H_{0}^{1}(\Omega)$ with \ $\mathfrak{S}_{r}:=\frac{N(r-1)+\gamma
r(N-2)}{N-2r}$}.$

Moreover, $u_{n}$ is uniformly bounded in $L^{\sigma_{r}}(\Omega)$ with
$\sigma_{r}:=\frac{Nr(1+\gamma)}{N-2r}.$

  * (ii)

For $r=\frac{N}{2}$ and $N\geq 2$. Then

    1. 1.

For $0\leq\gamma\leq 1$ and $\beta>0$ satisfying (2.3),

$\mathfrak{H}\left(\frac{u_{n}}{2}\right)\ \text{is uniformly bounded in}\
H_{0}^{1}(\Omega)\quad\text{where}\quad\mathfrak{H}(t):=\exp\left(\beta
t\right)-1$

and there exist constants $C_{1},C_{2}$ depending upon $\beta,S(N),f,|\Omega|$
but independent of $n$ such that

$\int_{\Omega}\exp\left(\frac{\beta Nu_{n}}{N-2}\right)\leq
C_{1}\quad\text{when}\
N>2\quad\text{and}\quad\|u_{n}\|_{L^{\infty}(\Omega)}\leq
C_{2}\quad\text{when}\ N=2.$

    2. 2.

For $\tau\geq\gamma>1$ and $\alpha>0$ such that
$\frac{\alpha}{2^{3-2\tau}}\max\\{1,\alpha^{\tau}\\}<\frac{1}{\tau
S(N)\|f\|_{\frac{N}{2}}}$

$\mathfrak{D}\left(\frac{u_{n}}{2}\right)\ \text{is uniformly bounded in}\
H_{0}^{1}(\Omega)\quad\text{where}\quad\mathfrak{D}(t)=\left(\exp(\alpha
t)-1\right)^{\tau}$

and there exist constants $C_{3},C_{4}$ depending upon
$\alpha,\tau,S(N),f,|\Omega|$ but independent of $n$ such that

$\int_{\Omega}\exp\left(\frac{2N\alpha\tau u_{n}}{N-2}\right)\leq
C_{3}\quad\text{when}\
N>2\quad\text{and}\quad\|u_{n}\|_{L^{\infty}(\Omega)}\leq
C_{4}\quad\text{when}\ N=2.$

In particular, $u_{n}$ is uniformly bounded in $L^{s}(\Omega)$ for every
$s\in[1,\infty)$ and $N>2.$

  * (iii)

For $r>\frac{N}{2}$ and $N\geq 2$, $u_{n}$ is uniformly bounded in
$L^{\infty}(\Omega).$

  * (iv)

For $r\in[1,r^{\sharp})$, $0\leq\gamma<1$ and $N>2$, $u_{n}$ is uniformly
bounded in $W_{0}^{1,q}(\Omega)$ with $q:=\frac{Nr(1+\gamma)}{N-r(1-\gamma)}$
and for any $r\in[1,\infty]$, $\gamma=1$ or $r\geq r^{\sharp}$ and
$0\leq\gamma<1$, $u_{n}$ is uniformly bounded in $H_{0}^{1}(\Omega).$

###### Proof.

Let $n\in\mathbb{N}$ and $u_{n}$ be the weak solution of the problem $(P_{n})$
given by Lemma 3.2. To prove the uniform estimates for the sequence
$\\{u_{n}\\}_{n\in\mathbb{N}}$, we divide the proof into four steps.
Step 1: Since, $u_{n}\in L^{\infty}(\Omega)\cap H_{0}^{1}(\Omega)$ and
positive, then for any $\epsilon>0$ and $\mathfrak{S}>0$,
$(u_{n}+\epsilon)^{\mathfrak{S}}-\epsilon^{\mathfrak{S}}$ belongs to
$H_{0}^{1}(\Omega)$, therefore, an admissible test function in (3.5). Taking
it so for $\epsilon\in(0,\frac{1}{n})$ and $\mathfrak{S}\in[\gamma,\infty)$,
it yields

$\begin{split}\int_{\Omega}\nabla
u_{n}\cdot\nabla(u_{n}+\epsilon)^{\mathfrak{S}}\leavevmode\nobreak\
dx&+\frac{C(N,s)}{2}\int_{\mathbb{R}^{N}}\int_{\mathbb{R}^{N}}\frac{(u_{n}(x)-u_{n}(y))((u_{n}+\epsilon)^{\mathfrak{S}}(x)-(u_{n}+\epsilon)^{\mathfrak{S}}(y))}{|x-y|^{N+2s}}\leavevmode\nobreak\
dx\leavevmode\nobreak\ dy\\\
&\leq\int_{\Omega}\frac{f_{n}(x)}{(u_{n}+\frac{1}{n})^{\gamma}}(u_{n}+\epsilon)^{\mathfrak{S}}\leavevmode\nobreak\
dx\leq\int_{\Omega}f_{n}(x)(u_{n}+\epsilon)^{\mathfrak{S}-\gamma}\leavevmode\nobreak\
dx.\end{split}$ (4.1)

Passing $\epsilon\to 0$ in the above estimate via Fatou’s theorem and using
Lemma 6.3 (i), we obtain

$\begin{split}\frac{4\mathfrak{S}}{(\mathfrak{S}+1)^{2}}\int_{\Omega}|\nabla
u_{n}^{\frac{\mathfrak{S}+1}{2}}|^{2}\leavevmode\nobreak\
dx&+\frac{2C(N,s)\mathfrak{S}}{(\mathfrak{S}+1)^{2}}\int_{\mathbb{R}^{N}}\int_{\mathbb{R}^{N}}\frac{((u_{n}^{\frac{\mathfrak{S}+1}{2}}(x)-u_{n}^{\frac{\mathfrak{S}+1}{2}}(y))^{2}}{|x-y|^{N+2s}}\leavevmode\nobreak\
dx\leavevmode\nobreak\ dy\\\
&\leq\int_{\Omega}f_{n}(x)u_{n}^{\mathfrak{S}-\gamma}\leavevmode\nobreak\
dx.\end{split}$ (4.2)

In order to estimate the R.H.S. term of (4.2), we choose
$\mathfrak{S}=\mathfrak{S}_{r}$ such that

$\mathfrak{S}_{1}=\gamma\quad\text{for}\ N\geq
2\quad\text{and}\quad\frac{(\mathfrak{S}_{r}-\gamma)r}{(r-1)}=\frac{(\mathfrak{S}_{r}+1)}{2}\frac{2N}{N-2}\quad\text{for}\
1<r<\frac{N}{2},\ N>2.$ (4.3)

By applying the Hölder inequality in view of the above choice of
$\mathfrak{S}_{r}$ with $r>1$ and $N>2$, we get

$\begin{split}\int_{\Omega}f_{n}(x)u_{n}^{\mathfrak{S}_{r}-\gamma}\leavevmode\nobreak\
dx&\leq\|f_{n}\|_{L^{r}(\Omega)}\left(\int_{\Omega}u_{n}^{(\mathfrak{S}_{r}-\gamma)r^{\prime}}\leavevmode\nobreak\
dx\right)^{\frac{1}{r^{\prime}}}=\|f_{n}\|_{L^{r}(\Omega)}\left(\int_{\Omega}\left(u_{n}^{\frac{(\mathfrak{S}_{r}+1)}{2}}\right)^{\frac{2N}{N-2}}\leavevmode\nobreak\
dx\right)^{\frac{1}{r^{\prime}}}\\\
&\leq\|f_{n}\|_{L^{r}(\Omega)}\left(S(N)\int_{\Omega}|\nabla
u_{n}^{\frac{\mathfrak{S}_{r}+1}{2}}|^{2}\leavevmode\nobreak\
dx\right)^{\frac{N}{r^{\prime}(N-2)}}\end{split}$ (4.4)

where the last inequality follows from Sobolev’s embeddings and $S(N)$ denotes
the best Sobolev constant. Combining (4.2) and (4.4), we get

$\int_{\Omega}|\nabla
u_{n}^{\frac{\mathfrak{S}_{r}+1}{2}}|^{2}\leavevmode\nobreak\
dx\leq\left(S(N)\right)^{\frac{N(r-1)}{N-2r}}\left(\|f_{n}\|_{L^{r}(\Omega)}\frac{(\mathfrak{S}_{r}+1)^{2}}{4\mathfrak{S}_{r}}\right)^{\frac{(N-2)r}{N-2r}}\leq
C(\|f\|_{L^{r}(\Omega)},\mathfrak{S}_{r},N).$ (4.5)

The above estimates in the case $r=1$ and $N\geq 2$ holds trivially. An easy
computation with the above choice of $\mathfrak{S}_{r}$ implies

$\frac{N(\mathfrak{S}_{r}+1)}{N-2}=\frac{Nr(1+\gamma)}{N-2r}=\sigma_{r}$ (4.6)

Then, by using the uniform estimates in (4.5) and Sobolev’s embeddings, we
obtain

$\\{u_{n}\\}_{n\in\mathbb{N}}\ \text{is uniformly bounded in}\
L^{\sigma_{r}}(\Omega)\ \text{when $1\leq r<\frac{N}{2},\ N>2$}.$

Step 2: Let $r=\frac{N}{2}$ and $N\geq 2$. In case of weak singularity i.e.
$0\leq\gamma\leq 1$, let us consider an increasing, convex and locally
Lipschitz function $\mathfrak{H}:\mathbb{R}\to\mathbb{R}$ defined as

$\mathfrak{H}(t):=\exp\left(\beta t\right)-1$

where $\beta>0$ whose exact choice will be determined later. Since $u_{n}\in
L^{\infty}(\Omega)\cap H_{0}^{1}(\Omega)$ for every $n\in\mathbb{N}$, we get
$\mathfrak{H}(u_{n})\in H_{0}^{1}(\Omega).$ A simple computations leads to

$\int_{\Omega}\nabla u_{n}\cdot\nabla\mathfrak{H}(u_{n})\leavevmode\nobreak\
dx=\int_{\Omega}\mathfrak{H}^{\prime}(u_{n})|\nabla
u_{n}|^{2}\leavevmode\nobreak\
dx=\frac{4}{\beta}\int_{\Omega}\left|\nabla\mathfrak{H}\left(\frac{u_{n}}{2}\right)\right|^{2}\leavevmode\nobreak\
dx$ (4.7)

Now, by testing the energy formulation (3.5) with $\mathfrak{H}(u_{n})$ and
using (4.7), we obtain

$\begin{split}\frac{4}{\beta}&\int_{\Omega}\left|\nabla\mathfrak{H}\left(\frac{u_{n}}{2}\right)\right|^{2}\leavevmode\nobreak\
dx\\\ &\leq\int_{\Omega}\nabla
u_{n}\cdot\nabla\mathfrak{H}(u_{n})\leavevmode\nobreak\
dx+\frac{C(N,s)}{2}\int_{\mathbb{R}^{N}}\int_{\mathbb{R}^{N}}\frac{(u_{n}(x)-u_{n}(y))(\mathfrak{H}(u_{n})(x)-\mathfrak{H}(u_{n})(y))}{|x-y|^{N+2s}}\leavevmode\nobreak\
dx\leavevmode\nobreak\ dy\\\
&\leq\int_{\Omega}\frac{f_{n}(x)}{(u_{n}+\frac{1}{n})^{\gamma}}\mathfrak{H}(u_{n})\leavevmode\nobreak\
dx\leq\int_{\Omega}\frac{f(x)\mathfrak{H}(u_{n})}{u_{n}^{\gamma}}\leavevmode\nobreak\
dx=\int_{\Omega}\frac{f(x)\left(\exp(\beta
u_{n})-1\right)}{u_{n}^{\gamma}}\leavevmode\nobreak\ dx\end{split}$ (4.8)

Case 1: $0<\gamma\leq 1.$
To estimate the last term in (4.8), we use Lemma 6.3(iii) and $\exp(t)-1\leq
t\exp(t)$ for every $t\geq 0$,

$\begin{split}\int_{\Omega}\left|\nabla\mathfrak{H}\left(\frac{u_{n}}{2}\right)\right|^{2}\leavevmode\nobreak\
dx&\leq\frac{\beta}{4}\left(\beta\gamma^{-1}\right)^{\gamma}\int_{\Omega\cap\\{u_{n}\leq
1\\}}f(x)\left(\frac{\exp\left(\beta\gamma^{-1}u_{n}\right)-1}{\beta\gamma^{-1}u_{n}}\right)^{\gamma}\leavevmode\nobreak\
dx\\\ &\quad\quad+\frac{\beta}{4}\int_{\Omega\cap\\{u_{n}\geq
1\\}}f(x)\exp\left(\beta u_{n}\right)\leavevmode\nobreak\ dx\\\ &\leq
C^{\sharp}\int_{\Omega}f(x)\left(\exp\left(\frac{\beta
u_{n}}{2}\right)\right)^{2}\leavevmode\nobreak\
dx=C^{\sharp}\int_{\Omega}f(x)\left(\mathfrak{H}\left(\frac{u_{n}}{2}\right)+1\right)^{2}\leavevmode\nobreak\
dx\end{split}$

where
$C^{\sharp}:=\frac{\beta}{4}\max\\{1,\left(\beta\gamma^{-1}\right)^{\gamma}\\}.$
Case 2: $\gamma=0.$
In this case the estimate in (4.8) takes the following form

$\begin{split}\int_{\Omega}\left|\nabla\mathfrak{H}\left(\frac{u_{n}}{2}\right)\right|^{2}\leavevmode\nobreak\
dx&\leq C^{\sharp}\int_{\Omega}f(x)\left(\exp\left(\frac{\beta
u_{n}}{2}\right)\right)^{2}\leavevmode\nobreak\
dx=C^{\sharp}\int_{\Omega}f(x)\left(\mathfrak{H}\left(\frac{u_{n}}{2}\right)+1\right)^{2}\leavevmode\nobreak\
dx\end{split}$

with $C^{\sharp}=\frac{\beta}{4}.$ Now, by using the fact that
$(\mathfrak{H}(t)+1)^{2}\leq 2(\mathfrak{H}(t))^{2}+1)$ and Hölder inequality
with exponents $\frac{N}{2}$ and $\frac{N}{N-2}$ , we further estimate

$\begin{split}\int_{\Omega}\left|\nabla\mathfrak{H}\left(\frac{u_{n}}{2}\right)\right|^{2}\leavevmode\nobreak\
dx&\leq 2C^{\sharp}\left(\int_{\Omega}f(x)\leavevmode\nobreak\
dx+\int_{\Omega}f(x)\left(\mathfrak{H}\left(\frac{u_{n}}{2}\right)\right)^{2}\leavevmode\nobreak\
dx\right)\\\ &\leq
2C^{\sharp}\left(\|f\|_{1}+\|f\|_{\frac{N}{2}}\left\|\mathfrak{H}\left(\frac{u_{n}}{2}\right)\right\|_{\frac{2N}{N-2}}^{2}\right).\end{split}$

Now, by choosing $\beta$ small enough such that

$\frac{2}{S(N)\|f\|_{\frac{N}{2}}}\geq\left\\{\begin{array}[]{ll}\beta\max\\{1,\left(\beta\gamma^{-1}\right)^{\gamma}\\}&\text{
if }\ 0<\gamma\leq 1,\vspace{0.1cm}\\\ \beta&\text{ if }\
\gamma=0,\end{array}\right.$

and using Sobolev’s embeddings and $|\Omega|<\infty$, we obtain

$\left\|\nabla\mathfrak{H}\left(\frac{u_{n}}{2}\right)\right\|_{2}^{2}\leq\frac{2C^{\sharp}\|f\|_{1}}{\left(1-2C^{\sharp}S(N)\|f\|_{\frac{N}{2}}\right)}$

and

$\int_{\Omega}\exp\left(\frac{\beta Nu_{n}}{N-2}\right)\leavevmode\nobreak\
dx\leq C$

where $C$ depends upon $\beta,S(N),\|f\|_{\frac{N}{2}},|\Omega|$ but
independent of $n.$

For the case of strong singularity i.e. $\gamma>1$, we consider the following
locally Lipschitz and increasing function
$\mathfrak{D}:\mathbb{R}\to\mathbb{R}^{+}\cup\\{0\\}$ defined as

$\mathfrak{D}(t)=\left(\exp(\alpha t)-1\right)^{\tau}$

where $\tau\geq\gamma>1$ and $\alpha>0$ whose exact choices will be
highlighted later. The regularity properties of the sequence $u_{n}$ and Lemma
6.3(iii) yields $\mathfrak{D}(u_{n})\in H_{0}^{1}(\Omega)$ and

$\begin{split}\int_{\Omega}\nabla
u_{n}\cdot\nabla\mathfrak{D}(u_{n})\leavevmode\nobreak\
dx&=\int_{\Omega}\mathfrak{D}^{\prime}(u_{n})|\nabla
u_{n}|^{2}\leavevmode\nobreak\
dx=\alpha\tau\int_{\Omega}\left|\exp\left(\frac{\alpha
u_{n}}{2}\right)\left(\exp(\alpha u_{n})-1\right)^{\frac{\tau-1}{2}}\nabla
u_{n}\right|^{2}\leavevmode\nobreak\ dx\\\
&\geq\alpha\tau\int_{\Omega}\left|\exp\left(\frac{\alpha
u_{n}}{2}\right)\left(\exp\left(\frac{\alpha
u_{n}}{2}\right)-1\right)^{\tau-1}\nabla u_{n}\right|^{2}\leavevmode\nobreak\
dx\\\
&=\frac{4}{\alpha\tau}\int_{\Omega}\left|\nabla\left(\exp\left(\frac{\alpha
u_{n}}{2}\right)-1\right)^{\tau}\right|^{2}\leavevmode\nobreak\
dx=\frac{4}{\alpha\tau}\int_{\Omega}\left|\nabla\mathfrak{D}\left(\frac{u_{n}}{2}\right)\right|^{2}\leavevmode\nobreak\
dx\end{split}$ (4.9)

By taking $\mathfrak{D}(u_{n})$ as a test function in the equation (3.5) and
using (4.9), we obtain

$\begin{split}\frac{4}{\alpha\tau}&\int_{\Omega}\left|\nabla\mathfrak{D}\left(\frac{u_{n}}{2}\right)\right|^{2}\leavevmode\nobreak\
dx\\\ &\leq\int_{\Omega}\nabla
u_{n}\cdot\nabla\mathfrak{D}(u_{n})\leavevmode\nobreak\
dx+\frac{C(N,s)}{2}\int_{\mathbb{R}^{N}}\int_{\mathbb{R}^{N}}\frac{(u_{n}(x)-u_{n}(y))(\mathfrak{D}(u_{n})(x)-\mathfrak{D}(u_{n})(y))}{|x-y|^{N+2s}}\leavevmode\nobreak\
dx\leavevmode\nobreak\ dy\\\
&\leq\int_{\Omega}\frac{f_{n}(x)}{(u_{n}+\frac{1}{n})^{\gamma}}\mathfrak{D}(u_{n})\leavevmode\nobreak\
dx\leq\int_{\Omega}\frac{f(x)\mathfrak{D}(u_{n})}{u_{n}^{\gamma}}\leavevmode\nobreak\
dx=\int_{\Omega}\frac{f(x)\left(\exp(\alpha
u_{n})-1\right)^{\tau}}{u_{n}^{\gamma}}\leavevmode\nobreak\ dx\\\
&\leq\alpha^{\tau}\int_{\Omega\cap\\{u_{n}\leq
1\\}}f(x)\left(\frac{\exp\left(\alpha u_{n}\right)-1}{\alpha
u_{n}}\right)^{\tau}\leavevmode\nobreak\ dx+\int_{\Omega\cap\\{u_{n}\geq
1\\}}f(x)\exp\left(\alpha\tau u_{n}\right)\leavevmode\nobreak\ dx\\\
&\leq\max\\{1,\alpha^{\tau}\\}\int_{\Omega}f(x)\exp\left(\alpha\tau
u_{n}\right)\leavevmode\nobreak\ dx\end{split}$

where to write the last two inequalities we have used $t^{\tau}\leq
t^{\gamma}$ for $t\in(0,1]$ and $\exp(t)-1\leq t\exp(t)$ for every $t\geq 0.$
Now, by using the fact that $\exp(2\tau
t)=\left(\mathfrak{D}^{\frac{1}{\tau}}(t)+1\right)^{2\tau}\leq
2^{2\tau-1}(\mathfrak{D}(t))^{2}+1)$ and Hölder inequality with exponents
$\frac{N}{2}$ and $\frac{N}{N-2}$, we obtain

$\begin{split}\int_{\Omega}\left|\nabla\mathfrak{D}\left(\frac{u_{n}}{2}\right)\right|^{2}\leavevmode\nobreak\
dx&\leq C^{*}\int_{\Omega}f(x)\left(\exp\left(\frac{\alpha
u_{n}}{2}\right)\right)^{2\tau}\leavevmode\nobreak\ dx\leq
2^{2\tau-1}C^{*}\int_{\Omega}f(x)\left(\mathfrak{D}^{\frac{1}{\tau}}\left(\frac{u_{n}}{2}\right)+1\right)^{2\tau}\leavevmode\nobreak\
dx\\\ &\leq 2^{2\tau-1}C^{*}\left(\int_{\Omega}f(x)\leavevmode\nobreak\
dx+\int_{\Omega}f(x)\left(\mathfrak{D}\left(\frac{u_{n}}{2}\right)\right)^{2}\leavevmode\nobreak\
dx\right)\\\ &\leq
2^{2\tau-1}C^{*}\left(\|f\|_{1}+\|f\|_{\frac{N}{2}}\left\|\mathfrak{D}\left(\frac{u_{n}}{2}\right)\right\|_{\frac{2N}{N-2}}^{2}\right)\quad\quad
C^{*}:=\frac{\alpha\tau}{4}\max\\{1,\alpha^{\tau}\\}.\end{split}$

Now, by choosing $\alpha$ small enough such that
$\alpha\max\\{1,\alpha^{\tau}\\}<\frac{2^{3-2\tau}}{\tau
S(N)\|f\|_{\frac{N}{2}}}$ and using Sobolev’s embeddings, we obtain

$\left\|\nabla\mathfrak{D}\left(\frac{u_{n}}{2}\right)\right\|_{2}^{2}\leq\frac{2^{2\tau-1}C^{*}\|f\|_{1}}{\left(1-2^{2\tau-1}C^{*}S(N)\|f\|_{\frac{N}{2}}\right)}\quad\text{and}\quad\int_{\Omega}\exp\left(\frac{2N\alpha\tau
u_{n}}{N-2}\right)\leavevmode\nobreak\ dx\leq C$

where $C$ depends upon $\alpha,\tau,S(N),\|f\|_{\frac{N}{2}},|\Omega|$ but
independent of $n.$
Step 3: To prove the boundedness result when $r>\frac{N}{2}$, we use the
classical arguments from the seminal paper of Stampacchia [46]. Choosing
$G_{k}(u_{n}):=(u_{n}-k)^{+}\in H_{0}^{1}(\Omega)$ with $k\geq 1$ as a test
function in the energy formulation (3.5). Then, by using Sobolev embeddings,
the Hölder inequality, $f_{n}\leq f$ and Lemma 3.2, we get

$\begin{split}\left(\int_{A_{k}}|G_{k}(u_{n})|^{\frac{2N}{N-2}}\leavevmode\nobreak\
dx\right)^{\frac{N-2}{N}}&\leq\int_{A_{k}}|\nabla
G_{k}(u_{n})|^{2}\leavevmode\nobreak\ dx=\int_{\Omega}\nabla u_{n}\cdot\nabla
G_{k}(u_{n})\leavevmode\nobreak\ dx\\\
&\leq\int_{\Omega}\frac{f_{n}}{\left(u_{n}+\frac{1}{n}\right)^{\gamma}}G_{k}(u_{n})\leavevmode\nobreak\
dx\leq\int_{A_{k}}f_{n}(x)G_{k}(u_{n})\leavevmode\nobreak\ dx\\\ &\leq
C\|f\|_{L^{r}(\Omega)}\left(\int_{A_{k}}(G_{k}(u_{n}))^{\frac{2N}{N-2}}\right)^{\frac{N-2}{2N}}|A_{k}|^{1-\frac{N-2}{2N}-\frac{1}{r}}\end{split}$
(4.10)

where $A_{k}:=\\{x\in\Omega:u_{n}(x)\geq k\\}.$ Let $h>k\geq 1$, then
$A_{h}\subset A_{k}$ and $G_{k}(u_{n})\geq h-k$ for $x\in A_{h}.$ Now, by
manipulating the estimate in (4.10) with above facts, we obtain

$\begin{split}|h-k||A_{h}|^{\frac{N-2}{2N}}&\leq\left(\int_{A_{h}}|G_{k}(u_{n})|^{\frac{2N}{N-2}}\leavevmode\nobreak\
dx\right)^{\frac{N-2}{2N}}\leq\left(\int_{A_{k}}|G_{k}(u_{n})|^{\frac{2N}{N-2}}\leavevmode\nobreak\
dx\right)^{\frac{N-2}{2N}}\\\ &\leq
C\|f\|_{L^{r}(\Omega)}|A_{k}|^{1-\frac{N-2}{2N}-\frac{1}{r}}\end{split}$

which further implies

$|A(h)|\leq
C\frac{\|f\|_{L^{r}(\Omega)}^{\frac{2N}{(N-2)}}|A_{k}|^{\frac{2N}{N-2}\left(1-\frac{N-2}{2N}-\frac{1}{r}\right)}}{|h-k|^{\frac{2N}{N-2}}}.$

Since, $r>\frac{N}{2}$, we have that

$\frac{2N}{N-2}\left(1-\frac{N-2}{2N}-\frac{1}{r}\right)>1.$

Hence, we apply Lemma 6.4 with the choice of $\psi(k)=|A_{k}|$, consequently
there exists $k_{0}$ such that $\psi(k)=0$ for all $k\geq k_{0}$ and thus our
claim.
Step 4: Let $r\in[1,r^{\sharp})$ and $q<2$ defined as in the statement of
theorem. Observe that for any $r\in[1,r^{\sharp})$,
$\mathfrak{S}_{r}\in[\gamma,1)$ and

${\frac{(1-\mathfrak{S}_{r})q}{2-q}}=\frac{Nr(1+\gamma)}{N-2r}=\sigma_{r}.$
(4.11)

By applying Hölder inequality with exponent $\frac{2}{q}$ and $\frac{2}{2-q}$
and using claim in Step 1, we get

$\begin{split}\int_{\Omega}|\nabla u_{n}|^{q}\leavevmode\nobreak\
dx&=\int_{\Omega}\frac{|\nabla
u_{n}|^{q}}{u_{n}^{\frac{(1-\mathfrak{S}_{r})q}{2}}}u_{n}^{\frac{(1-\mathfrak{S}_{r})q}{2}}\leavevmode\nobreak\
dx\\\ &\leq\left(\int_{\Omega}u_{n}^{\mathfrak{S}_{r}-1}|\nabla
u_{n}|^{2}\leavevmode\nobreak\
dx\right)^{\frac{q}{2}}\left(\int_{\Omega}u_{n}^{\frac{(1-\mathfrak{S}_{r})q}{2-q}}\leavevmode\nobreak\
dx\right)^{\frac{2-q}{2}}\\\
&\leq\left(\frac{4}{(\mathfrak{S}_{r}+1)^{2}}\int_{\Omega}|\nabla
u_{n}^{\frac{\mathfrak{S}_{r}+1}{2}}|^{2}\leavevmode\nobreak\
dx\right)^{\frac{q}{2}}\left(\int_{\Omega}u_{n}^{\sigma_{r}}\right)^{\frac{2-q}{2}}\leq
C\ \text{(independent of $n$)}.\end{split}$

As a special case, when $r\in[1,\infty]$, $\gamma=1$ or $r\geq r^{\sharp}$ and
$0\leq\gamma<1$, the proof of Claim $(i)$ can repeated by taking
$\mathfrak{S}_{r}=1$ and which further implies that $u_{n}$ is uniformly
bounded in $H_{0}^{1}(\Omega).$ ∎

######  Remark 4.1.

  * •

From Lemma 4.1(i), (iv) and Lemma 3.2(iii) we observe that for any $\gamma\geq
1$, $r\geq 1$, and $0\leq\gamma<1$, $r\geq r^{\sharp}$, $u_{n}$ is uniformly
bounded in $H^{1}_{loc}(\Omega)$ since $\mathfrak{S}_{r}\geq 1.$

  * •

From (4.11), we observe that for any $0\leq\gamma<1$

$q=\frac{2\sigma}{1-\mathfrak{S}_{r}+\sigma}$

which implies that there is a kind of $``$continuity$"$ in the summability
exponent $q$. Precisely, as $r\to r^{\sharp}$, $\mathfrak{S}_{r}\to 1$ and
$q\to 2.$

######  Lemma 4.2.

Assume that $f\in L^{r}(\Omega)\setminus\\{0\\}$ is non-negative and
$(r,\gamma)\in\mathcal{P}_{r,\gamma}$. Let $u_{n}$ be the weak solution of the
problem $(P_{n})$. Then,

$u_{n}^{\frac{\mathfrak{S}+1}{2}}\ \text{is uniformly bounded in}\
H_{0}^{1}(\Omega)\ \text{with}\
\mathfrak{S}\in\left\\{\begin{array}[]{ll}(0,\gamma]&\text{ if }\
\gamma+\frac{1}{r}<1,\vspace{0.1cm}\\\
\left(\gamma-1+\frac{1}{r},\gamma\right]&\text{ if }\ \gamma+\frac{1}{r}\geq
1.\end{array}\right.$

###### Proof.

Taking $(u_{n}+\epsilon)^{\mathfrak{S}}-\epsilon^{\mathfrak{S}}$ as a test
function in (3.5) for $\epsilon\in(0,\frac{1}{n})$ and
$\mathfrak{S}\in(0,\gamma]$ and using lower boundary estimates in Theorem
3.10, we obtain

$\begin{split}\int_{\Omega}\nabla
u_{n}\cdot\nabla(u_{n}+\epsilon)^{\mathfrak{S}}\leavevmode\nobreak\
dx&+\frac{C(N,s)}{2}\int_{\mathbb{R}^{N}}\int_{\mathbb{R}^{N}}\frac{(u_{n}(x)-u_{n}(y))((u_{n}+\epsilon)^{\mathfrak{S}}(x)-(u_{n}+\epsilon)^{\mathfrak{S}}(y))}{|x-y|^{N+2s}}\leavevmode\nobreak\
dx\leavevmode\nobreak\ dy\\\
&\leq\int_{\Omega}\frac{f_{n}(x)}{(u_{n}+\epsilon)^{\gamma-\mathfrak{S}}}\leavevmode\nobreak\
dx\leq
C\int_{\Omega}\frac{f(x)}{\delta^{\gamma-\mathfrak{S}}(x)}\leavevmode\nobreak\
dx\\\ &\leq
C\|f\|_{L^{r}(\Omega)}\left(\int_{\Omega}\delta^{\frac{-r(\gamma-\mathfrak{S})}{r-1}}(x)\leavevmode\nobreak\
dx\right)^{\frac{r-1}{r}}\\\ &\leq
C(\|f\|_{L^{r}(\Omega)},N,\Omega,r)\quad\text{if}\quad\frac{-r(\gamma-\mathfrak{S})}{r-1}>-1\
\text{and}\ \mathfrak{S}>0\end{split}$ (4.12)

where the last condition is equivalent to

$\frac{-r(\gamma-\mathfrak{S})}{r-1}>-1\ \text{and}\
\mathfrak{S}>0\quad\text{if and only
if}\quad\mathfrak{S}\in\left\\{\begin{array}[]{ll}(0,\gamma]&\text{ if }\
\gamma+\frac{1}{r}<1,\vspace{0.1cm}\\\
\left(\gamma-1+\frac{1}{r},\gamma\right]&\text{ if }\ \gamma+\frac{1}{r}\geq
1.\end{array}\right.$

Now, by using the same arguments as in proof of Lemma 4.1(i) and passing
$\epsilon\to 0$, we obtain our claim. ∎

###  4.2 Green’s function estimates

In this section, we prove series of Lemma involving the upper and lower
estimate of the action of Green operator on the logarithmic perturbation of
the distance function. Similar type of Green estimates for large class of
nonlocal operator is proved in [9]. For estimates near the boundary i.e. in
$\Omega_{\eta}$, we partition the set into the following five components: For
this, we partition the set $\Omega_{\eta}$ into the following five components

$\Omega_{1}:=B(x,\delta(x)/2),\quad\Omega_{2}:=\Omega_{\eta}\setminus B(x,1),$
$\Omega_{3}:=\\{y:\delta(y)<\delta(x)/2\\}\cap
B(x,1),\quad\Omega_{4}:=\left\\{y:\frac{3\delta(x)}{2}<\delta(y)<\eta\right\\}\cap
B(x,1),$
$\Omega_{5}:=\left\\{y:\frac{\delta(x)}{2}<\delta(y)<\frac{3\delta(x)}{2}\right\\}\cap\left(B(x,1)\setminus
B(x,\delta(x)/2)\right),$

and set $\ell(t)=\ln\left(\frac{\mathcal{D}_{\Omega}}{t}\right).$ For
$x\in\Omega_{\eta}$, we denote $\phi_{x}:B(x,1)\to B(0,1)$ be a diffeomorphism
such that

$\displaystyle\phi_{x}(\Omega\cap
B(x,1))=B(0,1)\cap\\{y\in\mathbb{R}^{N}:y\cdot e_{N}>0\\},$
$\displaystyle\phi_{x}(y)\cdot e_{N}=\delta(y)\ \text{for}\ y\in B(x,1)\
\text{and}\ \phi_{x}(x)=\delta(x)e_{N}.$

######  Lemma 4.3.

For $\Xi\in[0,1)$, we have

$\mathbb{G}^{\Omega}\left[\frac{\ell^{-\Xi}(\delta(\cdot))}{\delta(\cdot)}\right](x)\gtrsim\delta(x)\ell^{1-\Xi}(\delta(x))\quad\forall\
x\in\Omega.$

###### Proof.

Let $\eta>0$ small. We begin by splitting the integrals over two regions
$\Omega_{\eta}$ and $\Omega\setminus\Omega_{\eta}$ as follows

$\displaystyle\mathbb{G}^{\Omega}\left[\frac{\ell^{-\Xi}(\delta(\cdot))}{\delta(\cdot)}\right](x)$
$\displaystyle=\sum_{i=1}^{5}\int_{\Omega_{i}}\frac{G^{\Omega}(x,y)}{\delta(y)}\ell^{-\Xi}(\delta(y))\leavevmode\nobreak\
dy+\int_{\Omega\setminus\Omega_{\eta}}\frac{G^{\Omega}(x,y)}{\delta(y)}\ell^{-\Xi}(\delta(y))\leavevmode\nobreak\
dy$ (4.13) $\displaystyle:=\sum_{i=1}^{5}I_{1}^{(i)}(x)+I_{2}(x).$

As we know that $I_{1}^{(i)}(x),I_{2}(x)\geq 0$ for every $x\in\Omega$ and
$i=1,2,\dots,5$, so it is enough find the lower estimate of the term
$I_{1}^{(4)}$ when $x\in\Omega_{\frac{\eta}{2}}$, and the lower estimate of
the term $I_{2}$ when $x\in\Omega\setminus\Omega_{\frac{\eta}{2}}.$ For
$y\in\Omega_{4}$ and $x\in\Omega_{\frac{\eta}{2}}$, we have

$\left(\frac{\delta(x)\delta(y)}{|x-y|^{2}}\wedge
1\right)\asymp\frac{\delta(x)\delta(y)}{|x-y|^{2}},\quad\ell^{-\Xi}(\delta(x))\leq\ell^{-\Xi}(\delta(y)).$

Therefore

$\displaystyle
I_{1}(x)\geq\int_{\Omega_{4}}\frac{G^{\Omega}(x,y)}{\delta(y)}\ell^{-\Xi}(\delta(y))\leavevmode\nobreak\
dy\geq\ell^{-\Xi}(\delta(x))\int_{\Omega_{4}}\frac{G^{\Omega}(x,y)}{\delta(y)}\leavevmode\nobreak\
dy:=\ell^{-\Xi}(\delta(x))\mathcal{J}$ (4.14)

Now, by using estimates on [2, Proof of Lemma 3.3, Page 40] and performing
change of variables via diffeomorphism $\phi_{x}$, we obtain

$\begin{split}&\mathcal{J}\gtrsim\delta(x)\int_{\\{\frac{3\delta(x)}{2}<w_{N}<\eta\\}\cap
B(0,1)}\frac{1}{(|\delta(x)-w_{N}|+|w^{\prime}|)^{N}}\leavevmode\nobreak\
dw_{N}dw^{\prime}\\\
&=\delta(x)\int_{3/2}^{\eta/\delta(x)}\int_{0}^{1/\delta(x)}\frac{t^{N-2}}{(|1-h|+t)^{N}}\leavevmode\nobreak\
dt\leavevmode\nobreak\ dh\\\
&=\delta(x)\int_{3/2}^{\eta/\delta(x)}\frac{1}{(h-1)}\int_{0}^{1/(h-1)\delta(x)}\frac{r^{N-2}}{(1+r)^{N}}\leavevmode\nobreak\
dr\leavevmode\nobreak\ dh\\\
&\gtrsim\delta(x)\int_{3/2}^{\eta/\delta(x)}\frac{1}{(h-1)}\int_{1}^{1/(h-1)\delta(x)}\frac{1}{(1+r)^{2}}\leavevmode\nobreak\
dr\leavevmode\nobreak\
dh\gtrsim\delta(x)\int_{3/2}^{\eta/\delta(x)}\frac{1}{(h-1)}\leavevmode\nobreak\
dh\\\
&\gtrsim\delta(x)\int_{3/2}^{\eta/\delta(x)}\frac{1}{h}\leavevmode\nobreak\
dh=\delta(x)\left(\ln\left(\frac{\eta}{\delta(x)}\right)-\ln\left(\frac{3}{2}\right)\right).\end{split}$
(4.15)

By combining (4.13), (4.14) and (4.15), we obtain

$\mathbb{G}^{\Omega}\left[\frac{1}{\delta(\cdot)}\ell^{-\Xi}(\delta(\cdot))\right](x)\geq
I_{1}(x)\geq C\delta(x)\ell^{1-\Xi}(\delta(x))\ \text{for}\
x\in\Omega_{\frac{\eta}{2}}$ (4.16)

where the constant is independent of the parameter $\Xi.$ Let
$x\in\Omega\setminus\Omega_{\frac{\eta}{2}}>0.$ Since the operator
$\mathbb{G}^{\Omega}$ maps
$L^{\infty}_{c}(\Omega)\to\delta^{\gamma}C(\overline{\Omega})$ ([2, Theorem
2.10]), therefore we have the following estimates:

$\begin{split}I_{2}(x)&=\int_{\Omega\setminus\Omega_{\eta}}\frac{G^{\Omega}(x,y)}{\delta(y)}\ell^{-\Xi}(\delta(y))\leavevmode\nobreak\
dy\geq\ell^{-\Xi}(\eta)\mathbb{G}^{\Omega}\left[\frac{\chi_{\Omega\setminus\Omega_{\eta}}}{\delta}\right](x)\\\
&\geq C\ell^{-\Xi}\left(\frac{\eta}{2}\right)\delta(x)\geq
C(\eta)\ell^{1-\Xi}\left(\frac{\eta}{2}\right)\delta(x)\geq
C(\eta)\ell^{1-\Xi}\left(\delta(x)\right)\delta(x).\end{split}$ (4.17)

Finally, by combining (4.16) and (4.17), we obtain our claim. ∎

######  Lemma 4.4.

For $\Xi\in(0,1)$ following estimate holds:

$\mathbb{G}^{\Omega}\left[\frac{1}{\delta(\cdot)}\ell^{-\Xi}(\delta(\cdot))\right](x)\lesssim\delta(x)\ell^{1-\Xi}(\delta(x))\quad\forall\
x\in\Omega.$

###### Proof.

Let $I_{1}(x)$ and $I_{2}(x)$ as in (4.13). To derive the estimate for
$I_{1}$, we divide the proof into two cases depending upon the location of the
point $x.$
Case 1: $x\in\Omega_{\frac{\eta}{2}}.$
By partitioning the domain of integral $I_{1}(x)$ over
$\\{\Omega_{i}\\}_{i=1}^{5}$, we find the upper estimates over each subdomain
$\Omega_{i}$. Observing, for $y_{1}\in\Omega_{1}$ and
$y_{2}\in\cup_{i=2}^{5}\Omega_{i}$, we have

$\left(\frac{\delta(x)\delta(y_{1})}{|x-y_{1}|^{2}}\wedge 1\right)\asymp 1\
\text{and}\ \left(\frac{\delta(x)\delta(y_{2})}{|x-y_{2}|^{2}}\wedge
1\right)\asymp\frac{\delta(x)\delta(y_{2})}{|x-y_{2}|^{2}}.$

Now, by again using the change of variables via diffeomorphism $\phi_{x}$ and
[2, Proof of Lemma 3.3], we get estimates in each domain.
Estimate over $\Omega_{1}$: Choosing $\eta$ small enough such that
$0<\eta<\frac{2\mathcal{D}_{\Omega}}{\exp(1)}$, for some $c\in(0,1)$, we have,
for any $y\in\Omega_{1}$,

$\delta(y)\leq\frac{3}{2}\delta(x),\quad\ell^{-\Xi}(\delta(y))\leq
c^{-\Xi}\ell^{-\Xi}(\delta(x))\leq c^{-1}\ell^{-\Xi}(\delta(x)).$

Therefore

$\begin{split}\int_{\Omega_{1}}\frac{G^{\Omega}(x,y)}{\delta(y)}\ell^{-\Xi}(\delta(y))\leavevmode\nobreak\
dy&\leq
c^{-1}\ell^{-\Xi}(\delta(x))\frac{1}{\delta(x)}\int_{B(x,\delta(x)/2)}\frac{1}{|x-y|^{N-2}}\leavevmode\nobreak\
dy\\\ &\leq c^{-1}\ell^{-\Xi}(\delta(x))\delta(x)\leq
C\ell^{1-\Xi}(\delta(x))\delta(x),\end{split}$

where $C$ is independent of parameter $\Xi.$
Estimate over $\Omega_{2}$:

$\begin{split}\int_{\Omega_{2}}\frac{G^{\Omega}(x,y)}{\delta(y)}\ell^{-\Xi}(\delta(y))\leavevmode\nobreak\
dy&\leq\ell^{-\Xi}\left(\eta\right)\delta(x)\int_{\Omega_{2}}\frac{1}{|x-y|^{n}}\leavevmode\nobreak\
dy\leq C(\eta)\ \ell^{1-\Xi}(\delta(x))\delta(x),\end{split}$

where $C$ is independent of parameter $\Xi.$
Estimate over $\Omega_{3}$: We note that, for any $y\in\Omega_{3}$,
$\ell^{-\Xi}(\delta(y))\leq\ell^{-\Xi}(\delta(x)).$ Therefore

$\begin{split}\int_{\Omega_{3}}&\frac{G^{\Omega}(x,y)}{\delta(y)}\ell^{-\Xi}(\delta(y))\leavevmode\nobreak\
dy\leq\ell^{-\Xi}(\delta(x))\delta(x)\int_{\Omega_{3}}\frac{1}{|x-y|^{N}}\leavevmode\nobreak\
dy\\\
&\lesssim\ell^{-\Xi}(\delta(x))\delta(x)\int_{|w^{\prime}|<1}\int_{0}^{\delta(x)/2}\frac{1}{(|\delta(x)-w_{N}|+|w^{\prime}|)^{N}}\leavevmode\nobreak\
dw_{N}dw^{\prime}\\\
&\lesssim\ell^{-\Xi}(\delta(x))\delta(x)\int_{0}^{1/\delta(x)}\int_{0}^{1/2}\frac{h^{N-2}}{((1-t)+h)^{N}}\leavevmode\nobreak\
dt\leavevmode\nobreak\ dh\\\
&\lesssim\ell^{-\Xi}(\delta(x))\delta(x)\int_{0}^{1/\delta(x)}\frac{h^{N-2}}{(1+h)^{N}}\leavevmode\nobreak\
dh\leq C(\eta)\ell^{1-\Xi}(\delta(x))\delta(x).\end{split}$

Estimate over $\Omega_{4}$: We have

$\begin{split}\int_{\Omega_{4}}&\frac{G^{\Omega}(x,y)}{\delta(y)}\ell^{-\Xi}(\delta(y))\leavevmode\nobreak\
dy\leq\delta(x)\int_{\Omega_{4}}\frac{1}{|x-y|^{N}\ell^{\Xi}(\delta(y))}\leavevmode\nobreak\
dy\\\ &\lesssim\delta(x)\int_{\\{\frac{3\delta(x)}{2}<w_{N}<\eta\\}\cap
B(0,1)}\frac{1}{(|\delta(x)-w_{N}|+|w^{\prime}|)^{N}\ell^{\Xi}(w_{N})}\leavevmode\nobreak\
dw_{N}dw^{\prime}\\\
&=\delta(x)\int_{3/2}^{\eta/\delta(x)}\int_{0}^{1/\delta(x)}\frac{t^{N-2}}{(|1-h|+t)^{N}\ell^{\Xi}(h\delta(x))}\leavevmode\nobreak\
dt\leavevmode\nobreak\ dh\\\
&=\delta(x)\int_{3/2}^{\eta/\delta(x)}\frac{1}{(h-1)\ell^{\Xi}(h\delta(x))}\int_{0}^{1/(h-1)\delta(x)}\frac{r^{N-2}}{(1+r)^{N}}\leavevmode\nobreak\
dr\leavevmode\nobreak\ dh\\\
&\lesssim\delta(x)\int_{3/2}^{\eta/\delta(x)}\frac{1}{(h-1)\ell^{\Xi}(h\delta(x))}\int_{0}^{1/(h-1)\delta(x)}\frac{1}{(1+r)^{2}}\leavevmode\nobreak\
dr\leavevmode\nobreak\ dh\\\
&\lesssim\delta(x)\int_{3/2}^{\eta/\delta(x)}\frac{1}{(h-1)\ell^{\Xi}(h\delta(x))}\leavevmode\nobreak\
dh\lesssim\delta(x)\int_{3/2}^{\eta/\delta(x)}\frac{\ell^{-\Xi}(h\delta(x))}{h}\leavevmode\nobreak\
dh\\\
&\lesssim\delta(x)\int^{\ln(2\mathcal{D}_{\Omega}/3\delta(x))}_{\ln(\mathcal{D}_{\Omega}/\eta)}\frac{1}{t^{\Xi}}dt\leq\frac{C(\eta)}{(1-\Xi)}\ell^{1-\Xi}(\delta(x))\delta(x)\end{split}$

Estimate over $\Omega_{5}$: Again, by choosing $\eta$ small enough such that
$\eta<\frac{2A}{\exp(1)}$, for some $c\in(0,1)$, we have, for
$y\in\Omega_{5}$,

$\delta(y)\leq\frac{3}{2}\delta(x)\ \text{and}\ \ell^{-\Xi}(\delta(y))\leq
c^{-\Xi}\ell^{-\Xi}(\delta(x))\leq c^{-1}\ell^{-\Xi}(\delta(x)).$

Therefore

$\begin{split}\int_{\Omega_{5}}&\frac{G^{\Omega}(x,y)}{\delta(y)}\ell^{-\Xi}(\delta(y))\leavevmode\nobreak\
dy\lesssim\delta^{\gamma}(x)\int_{\Omega_{5}}\frac{1}{|x-y|^{N}\ell^{-\Xi}(\delta(y))}\leavevmode\nobreak\
dy\\\
&\lesssim\delta(x)\ell^{-\Xi}(\delta(x))\int_{\Omega_{5}}\frac{1}{|x-y|^{N}}\leavevmode\nobreak\
dy\lesssim\delta(x)\ell^{-\Xi}(\delta(x))\int_{\delta(x)/2}^{3\delta(x)/2}\int_{\delta(x)/2}^{1}\frac{t^{N-2}}{((\delta(x)-h)+t)^{N}}\leavevmode\nobreak\
dt\leavevmode\nobreak\ dh\\\
&\lesssim\delta(x)\ell^{-\Xi}(\delta(x))\int_{\delta(x)/2}^{1}t^{-1}\int_{-\delta(x)/2t}^{\delta(x)/2t}\frac{1}{(|r|+1)^{N}}\leavevmode\nobreak\
dr\leavevmode\nobreak\ dt\\\
&\lesssim\delta(x)\ell^{-\Xi}(\delta(x))\int_{1/2}^{1/\delta(x)}\rho^{-1}\int_{0}^{1/\rho}\frac{1}{(r+1)^{N}}drd\rho\leq
C(\eta)\ell^{1-\Xi}(\delta(x))\delta(x).\end{split}$

Case 2: $x\in\Omega\setminus\Omega_{\frac{\eta}{2}}.$
Using the estimates from the proof of Lemma 3.2 [2, p.37], for $\eta$ small
enough, we obtain

$\begin{split}I_{1}(x)&\leq\ell^{-\Xi}(\eta)\delta(x)\left(\int_{\Omega_{\eta/4}}\frac{G^{\Omega}(x,y)}{\delta(y)}\leavevmode\nobreak\
dy+\int_{\Omega_{\eta}\setminus\Omega_{\eta/4}}\frac{G^{\Omega}(x,y)}{\delta(y)}\leavevmode\nobreak\
dy\right)\\\ &\leq
C\ell^{-\Xi}(\eta)\delta(x)\left(1+\frac{1}{\eta}\right)\leq
C(\mathcal{D}_{\Omega},\eta)\ell^{1-\Xi}(\delta(x))\delta(x).\end{split}$

For any $x\in\Omega$, we have the following estimate for $I_{2}(x):$

$\begin{split}I_{2}(x)&\leq\frac{\ln^{-\Xi}\left(2\right)}{\eta}\mathbb{G}^{\Omega}[\chi_{\Omega}](x)\leq
C(\mathcal{D}_{\Omega})\ell^{1-\Xi}(\delta(x))\delta(x)\frac{1}{\eta\ell^{1-\Xi}\left(\frac{\eta}{2}\right)}\\\
&\leq
C(\mathcal{D}_{\Omega},\eta)\ell^{1-\Xi}(\delta(x))\delta(x).\end{split}$
(4.18)

Finally, by collecting all the estimates in $\\{\Omega_{i}\\}$ in Case 1, and
Case 2 and (4.18), we get the desired upper estimate. ∎

Denote

$\beta:=\frac{2\gamma+\zeta}{\gamma+1}\quad\kappa:=2-\beta=\frac{2-\zeta}{\gamma+1}$
(4.19)

######  Lemma 4.5.

For $0<\epsilon,\eta<1$, $N\geq 3$ and $\gamma>1$, there exist positive
constant $C_{1},C_{2}>0$ such that the following hold:

$\mathbb{G}^{\Omega}\left[\frac{1}{(\delta+\epsilon^{1/\kappa})^{\beta}}\chi_{\Omega_{\eta}}\right](x)\gtrsim\left(\frac{1}{2}(\delta(x)+\epsilon^{1/\kappa})^{\kappa}-\epsilon\right)\quad\forall\
x\in\Omega_{\frac{\eta}{2}}$ (4.20)

and

$\mathbb{G}^{\Omega}\left[\frac{1}{(\delta+\epsilon^{1/\kappa})^{\beta}}\right](x)\gtrsim(\delta(x)+\epsilon^{1/\kappa})^{\kappa}-\epsilon\quad\forall\
x\in\Omega_{\frac{\eta}{2}}^{c}.$ (4.21)

###### Proof.

Denote $\epsilon_{1}=\epsilon^{1/\kappa}$ and
$\mathcal{B}_{\delta}^{x}:=\\{y\in\Omega:|x-y|<\frac{\delta(x)}{2}\\}\subset\Omega_{\eta}.$
Fix $x\in\Omega_{\frac{\eta}{2}}$. Then, for $y\in\mathcal{B}_{\delta}^{x}$,
we have

$\left(\frac{\delta(x)}{|x-y|}\wedge 1\right)\geq
1,\left(\frac{\delta(y)}{|x-y|}\wedge 1\right)\geq 1\quad\text{and }\
\frac{1}{(\delta(y)+\epsilon_{1})^{\beta}}\geq\left(\frac{2}{3}\right)^{\beta}\frac{1}{(\delta(x)+\epsilon_{1})^{\beta}}.$
(4.22)

Using (4.22), we obtain the following

$\begin{split}\mathbb{G}^{\Omega}&\left[\frac{1}{(\delta+\epsilon_{1})^{\beta}}\chi_{\Omega_{\eta}}\right](x)=\int_{\Omega_{\eta}}\frac{\mathcal{G}^{\Omega}(x,y)}{(\delta(y)+\epsilon_{1})^{\beta}}\leavevmode\nobreak\
dy\\\ &\geq
C\int_{\mathcal{B}_{\delta}^{x}}\frac{1}{(\delta(y)+\epsilon_{1})^{\beta}}\frac{1}{|x-y|^{N-2}}\left(\frac{\delta(x)\delta(y)}{|x-y|^{2}}\wedge
1\right)\,dy,\\\ &\geq
C\int_{\mathcal{B}_{\delta}^{x}}\frac{1}{(\delta(y)+\epsilon_{1})^{\beta}}\frac{1}{|x-y|^{N-2}}\left(\frac{\delta(x)}{|x-y|}\wedge
1\right)\left(\frac{\delta(y)}{|x-y|}\wedge 1\right)\,dy,\\\
&\geq\frac{C^{\prime}}{(\delta(x)+\epsilon_{1})^{\beta}}\int_{\mathcal{B}_{\delta}^{x}}\frac{1}{|x-y|^{N-2}}\leavevmode\nobreak\
dy=\frac{C^{\prime\prime}\delta(x)^{2}}{(\delta(x)+\epsilon_{1})^{\beta}}\geq
C^{\prime\prime}\left(\frac{1}{2}(\delta(x)+\epsilon_{1})^{\kappa}-\frac{\epsilon_{1}^{2}}{(\delta(x)+\epsilon_{1})^{\beta}}\right)\\\
&\geq
C_{1}\left(\frac{1}{2}(\delta(x)+\epsilon^{1/\kappa})^{\kappa}-\epsilon\right).\end{split}$
(4.23)

For the second claim: let $x\in\Omega_{\frac{\eta}{2}}^{c}$. Then by using [2,
Theorem 3.4], we get

$\begin{split}\mathbb{G}^{\Omega}\left[\frac{1}{(\delta+\epsilon_{1})^{\beta}}\right](x)&\geq\frac{1}{(\mathcal{D}_{\Omega}+1)^{\beta}}\mathbb{G}^{\Omega}[\chi_{\Omega}](x)\geq
c\delta(x)\geq
c_{2}(\delta(x)+\epsilon^{1/\kappa})^{\kappa}-\epsilon,\end{split}$

where $c_{2}=c_{2}(\eta,\kappa,\mathcal{D}_{\Omega})$. ∎

###  4.3 Singular weights

####  4.3.1 Boundary behavior

Let $\varphi_{1}$ and $\lambda_{1}$ are first eigenfunction and eigenvalue for
the mixed local-nonlocal operator $(-\Delta)+(-\Delta)^{s}$ such that ([24,
Proposition 3.7] and [17, Proposition 5.3, 5.4])

$\varphi_{1}\asymp\delta^{\gamma}\ \text{and}\
\varphi_{1}=\lambda_{1}\mathbb{G}^{\Omega}[\varphi_{1}]\quad\text{in}\ \Omega$
(4.24)

######  Lemma 4.6.

Let $u_{n}$ be the weak solution of the problem $(P_{n})$, then for any
$\gamma>0$ and $\zeta\geq 0$, there exist a constant $C_{0}>0$ independent of
$n$ and $x\in\Omega$ such that $C_{0}\delta(x)\leq u_{n}(x)$. Moreover, for
$\zeta\in[0,2)$, we have

  1. 1.

if $\zeta+\gamma\leq 1$, there exist two constants $C_{1},C_{2}>0$ independent
of $n$ and $x\in\Omega$ such that

$C_{1}\delta(x)\leq u_{n}(x)\leq
C_{2}\delta(x)\left\\{\begin{array}[]{ll}1&\text{ if }\
\zeta+\gamma<1,\vspace{0.1cm}\\\
\ln\left(\frac{\mathcal{D}_{\Omega}}{\delta(x)}\right)&\text{ if }\
\zeta+\gamma=1,\end{array}\right.\ \text{for}\ x\in\Omega.$ (4.25)

  2. 2.

if $\zeta+\gamma>1$, there exist two constants $C^{1},C^{2}$ and $C^{3}>0$
independent of $n$ and $x\in\Omega$ such that

$C^{1}(\delta(x)+n^{\frac{-(1+\gamma)}{(2-\zeta)}})^{\frac{2-\zeta}{\gamma+1}}-C^{2}n^{-1}\leq
u_{n}(x)\leq C^{3}\delta^{\frac{2-\zeta}{\gamma+1}}(x)\ \text{for}\
x\in\Omega.$ (4.26)

###### Proof.

To prove the boundary behavior, we divide our study into three cases:
Case 1: $\zeta+\gamma\leq 1$
Using the integral representation in (3.10), lower boundary behavior in
(3.11), $f_{n}\leq f$, and Lemma 6.2, we get

$\begin{split}u_{n}(x)&=\int_{\Omega}\frac{G(x,y)f_{n}(y)}{(u_{n}+\frac{1}{n})^{\gamma}}\leavevmode\nobreak\
dy\leq\int_{\Omega}\frac{G(x,y)f(y)}{u_{n}^{\gamma}}\leavevmode\nobreak\ dy\\\
&\lesssim\mathbb{G}^{\Omega}\left[\frac{1}{\delta^{\gamma+\zeta}}\right]\lesssim\delta(x)\left\\{\begin{array}[]{ll}1&\text{
if }\ \zeta+\gamma<1,\vspace{0.1cm}\\\
\ln\left(\frac{\mathcal{D}_{\Omega}}{\delta(x)}\right)&\text{ if }\
\zeta+\gamma=1,\end{array}\right.\end{split}$ (4.27)

Using (4.27) and estimates in Theorem 3.10, we obtain our first claim.
Case 2: $\zeta+\gamma>1$.
For $n\in\mathbb{N}$ and $\kappa>0$, we define
$h_{\epsilon_{n}}:=\frac{1}{(\delta+\epsilon_{n}^{1/\kappa})^{2-\kappa}}\in
L^{\infty}(\Omega)$ with $\epsilon_{n}=\frac{1}{n}.$ Since $\gamma>1$ and
$\zeta<2$, we choose

$\kappa=\frac{2-\zeta}{\gamma+1}\quad\text{such
that}\quad\kappa\gamma+\zeta=2-\kappa.$

Then, from Lemma 3.1, there exists a unique positive weak solution
$v_{\epsilon_{n}}\in H_{0}^{1}(\Omega)\cap L^{\infty}(\Omega)$ of the
following problem

$\left\\{\begin{aligned} -\Delta
v_{\epsilon_{n}}+(-\Delta)^{s}v_{\epsilon_{n}}&=h_{\epsilon_{n}},\quad
v_{\epsilon_{n}}>0\quad&&\text{ in }\Omega,\\\ v_{\epsilon_{n}}&=0&&\text{ in
}\ \mathbb{R}^{N}\setminus\Omega,\end{aligned}\right.$ ($S_{\epsilon_{n}}$)

Set $\epsilon_{1}^{n}={\epsilon_{n}}^{1/\kappa}>0$. Using the uniqueness
property of $v_{\epsilon_{n}}$, integral representation of the solution via
Green function and Lemma 6.2, we obtain

$\begin{split}v_{\epsilon_{n}}=\mathbb{G}^{\Omega}\left[\frac{1}{(\delta+\epsilon_{1}^{n})^{2-\kappa}}\right]\leq\mathbb{G}^{\Omega}\left[\frac{1}{\delta^{2-\kappa}}\right]\leq
C\delta^{\kappa}\quad\text{in }\Omega\quad\text{if}\ \kappa>1,\end{split}$
(4.28)

and on the other hand for $0<\eta<1$, Lemma 4.5 gives

$\begin{split}v_{\epsilon_{n}}&(x)=\mathbb{G}^{\Omega}\left[\frac{1}{(\delta+\epsilon_{1}^{n})^{2-\kappa}}\right](x)\\\
&\begin{aligned}
\gtrsim\left\\{\begin{array}[]{ll}\mathbb{G}^{\Omega}\left[\frac{1}{(\delta+\epsilon_{1}^{n})^{2-\kappa}}\chi_{\Omega_{\eta}}\right](x)\gtrsim\left(\frac{1}{2}(\delta(x)+\epsilon_{1}^{n})^{\kappa}-\epsilon_{n}\right)&\text{
if }x\in\Omega_{\frac{\eta}{2}},\\\
(\delta(x)+\epsilon_{1}^{n})^{\kappa}-\epsilon_{n}&\text{ if
}x\in\Omega\setminus\Omega_{\frac{\eta}{2}},\\\
\end{array}\right.\end{aligned}\\\
&\gtrsim\left(\frac{1}{2}(\delta(x)+\epsilon_{1}^{n})^{\kappa}-\epsilon_{n}\right),\quad
x\in\Omega\quad\text{if}\ \gamma>1\ \text{and}\ \zeta<2.\\\ \end{split}$
(4.29)

Collecting the estimates in (4.28)-(4.29), there exists a constant
$C_{3},C_{4}>0$ independent of $n$ such that

$C_{3}\left(\frac{1}{2}(\delta(x)+\epsilon_{1}^{n})^{\kappa}-\epsilon_{n}\right)\leq
v_{\epsilon_{n}}(x)\leq C_{4}\delta(x)^{\kappa},\ x\in\Omega\quad\text{if}\
\gamma>1\ \text{and}\ \zeta<2.$ (4.30)

Define

$\underline{u}^{\epsilon_{n}}=c_{\eta}v_{\epsilon_{n}}\ \text{with}\
0<c_{\eta}<\frac{C_{1}}{C_{4}}\left(\frac{\eta}{2}\right)^{1-\kappa},$

where $C_{1},C_{4}$ are defined in (4.25) and (4.30) respectively. We note
that $c_{\eta}$ is independent of $\epsilon_{n}$ such that

$\underline{u}^{\epsilon_{n}}\leq u_{n}\ \text{in
}\Omega\setminus\Omega_{\frac{\eta}{2}}$

and

$\left(\underline{u}^{\epsilon_{n}}+\epsilon_{n}\right)\leq
2C_{4}c_{\eta}(\delta+\epsilon_{1}^{n})^{\kappa}+(1-C_{4}c_{\eta})\epsilon_{n}\quad\text{
in }\Omega.$

If
$2C_{4}c_{\eta}(\delta(x)+\epsilon_{1}^{n})^{\kappa}\geq(1-C_{4}c_{\eta})\epsilon_{n}$
then by choosing $\eta$ small enough such that
$c_{\eta}<\mathcal{G}_{1}^{\frac{1}{\gamma+1}}(4C_{4})^{-\frac{\gamma}{\gamma+1}}$,
we have

$\begin{split}(-\Delta)\underline{u}^{\epsilon_{n}}+(-\Delta)^{s}\underline{u}^{\epsilon_{n}}=c_{\eta}(\delta+\epsilon_{1}^{n})^{-(2-\kappa)}<\mathcal{G}_{1}(4C_{4}c_{\eta})^{-\gamma}(\delta+\epsilon_{1}^{n})^{-\kappa\gamma-\zeta}\leq
f_{\epsilon_{n}}(x)\left(\underline{u}^{\epsilon_{n}}+\epsilon_{n}\right)^{-\gamma}\quad\text{in}\
\Omega.\end{split}$

where $\mathcal{G}_{1}$ is defined in (3.4).
If
$2C_{4}c_{\eta}(\delta(x)+\epsilon_{1}^{n})^{\kappa}\leq(1-C_{4}c_{\eta})\epsilon_{n}$
then again by choosing $\eta$ small enough such that $C_{4}c_{\eta}<1$ and
$c_{\eta}<\mathcal{G}_{1}(2(1-C_{4}c_{\eta}))^{-q}$ , we have

$\begin{split}(-\Delta)\underline{u}^{\epsilon_{n}}+(-\Delta)^{s}\underline{u}^{\epsilon_{n}}&=c_{\eta}(\delta+\epsilon_{1}^{n})^{-(2-\kappa)}\leq\mathcal{G}_{1}(2(1-C_{4}c_{\eta}))^{-q}(\delta+\epsilon_{1}^{n})^{-\kappa\gamma-\zeta}\\\
&\leq f_{\epsilon_{n}}(x)(1-C_{4}c_{\eta})^{-q}(2{\epsilon_{n}})^{-\gamma}\leq
f_{\epsilon_{n}}(x)\left(\underline{u}^{\epsilon_{n}}+\epsilon_{n}\right)^{-\gamma}\quad\text{in}\
\Omega.\end{split}$

Now, by considering both the cases and applying the weak comparison principle
in $\Omega_{\frac{\eta}{2}}$ for $u_{n}$ and $\underline{u}^{\epsilon_{n}}$,
we get $\underline{u}^{\epsilon_{n}}\leq u_{n}$ in $\Omega$, namely there
exist constants $0<C^{1},C^{2}<\frac{1}{2}$ (by taking $\eta$ small enough)
such that

$C^{1}(\delta+\epsilon_{1}^{n})^{\kappa}-C^{2}{\epsilon_{n}}\leq
u_{n}\quad\text{in }\Omega.$ (4.31)

Finally, by using the integral representation, we obtain

$\begin{split}u_{n}&=\mathbb{G}^{\Omega}\left[\frac{f_{\epsilon_{n}}(x)}{(u_{n}+\epsilon_{n})^{\gamma}}\right]\leq\mathbb{G}^{\Omega}\left[\frac{f_{\epsilon_{n}}(x)}{(C^{1}(\delta+\epsilon_{1}^{n})^{\kappa}+(1-C^{2}){\epsilon_{n}})^{\gamma}}\right]\\\
&\lesssim\mathbb{G}^{\Omega}\left[\frac{1}{\delta^{\kappa\gamma+\zeta}}\right]=\mathbb{G}^{\Omega}\left[\frac{1}{\delta^{2-\kappa}}\right]\lesssim\delta^{\kappa}\quad\text{if}\
\gamma>1\ \text{and}\ \zeta<2.\end{split}$ (4.32)

∎

####  4.3.2 Sobolev regularity estimates

######  Lemma 4.7.

Let $\gamma>0$, $\zeta\in[0,2)$ and $u_{n}$ be the weak solution of the
problem $(P_{n})$. Then,

$u_{n}^{\frac{\mathfrak{L}+1}{2}}\ \text{is uniformly bounded in}\
H_{0}^{1}(\Omega)\ \text{for any}\
\mathfrak{L}>\left\\{\begin{array}[]{ll}0&\text{ if }\ \zeta+\gamma\leq
1,\vspace{0.1cm}\\\ \mathfrak{L}^{*}&\text{ if }\
\zeta+\gamma>1,\end{array}\right.\ \text{where}\
\mathfrak{L}^{*}:=\frac{\gamma+\zeta-1}{2-\zeta}.$

###### Proof.

Let $n\in\mathbb{N}$ and $u_{n}$ be the weak solution of the problem $(P_{n})$
given by Lemma 3.2. Since, $u_{n}\in L^{\infty}(\Omega)\cap H_{0}^{1}(\Omega)$
and positive, then for any $\epsilon>0$ and $\mathfrak{L}>0$,
$(u_{n}+\epsilon)^{\mathfrak{L}}-\epsilon^{\mathfrak{L}}$ belongs to
$H_{0}^{1}(\Omega)$, therefore, an admissible test function in (3.5). Taking
it so for $\epsilon\in(0,\frac{1}{n})$ and passing $\epsilon\to 0$ as in the
proof of Lemma 4.1(i), we obtain

$\begin{split}\frac{4\mathfrak{L}}{(\mathfrak{L}+1)^{2}}\int_{\Omega}|\nabla
u_{n}^{\frac{\mathfrak{L}+1}{2}}|^{2}\leavevmode\nobreak\
dx&+\frac{2C(N,s)\mathfrak{L}}{(\mathfrak{L}+1)^{2}}\int_{\mathbb{R}^{N}}\int_{\mathbb{R}^{N}}\frac{((u_{n}^{\frac{\mathfrak{L}+1}{2}}(x)-u_{n}^{\frac{\mathfrak{L}+1}{2}}(y))^{2}}{|x-y|^{N+2s}}\leavevmode\nobreak\
dx\leavevmode\nobreak\ dy\\\
&\leq\int_{\Omega}f_{n}(x)u_{n}^{\mathfrak{L}-\gamma}\leavevmode\nobreak\
dx:=\mathfrak{G}(u_{n}).\end{split}$ (4.33)

Now to estimate the right hand side term in (4.33), we divide the proof into
three cases. Using (3.4) and (4.25), we obtain the following:
Case 1: $\zeta+\gamma<1$

$\begin{split}\mathfrak{G}(u_{n})&\leq\int_{\Omega}f_{n}(x)u_{n}^{\mathfrak{L}-\gamma}\leavevmode\nobreak\
dx\leq\mathcal{G}_{2}\int_{\Omega}\frac{u_{n}^{\mathfrak{L}-\gamma}}{\delta^{\zeta}(x)}\leavevmode\nobreak\
dx\lesssim\int_{\Omega}u_{n}^{\mathfrak{L}-\gamma-\zeta}\leavevmode\nobreak\
dx\\\
&\lesssim\int_{\Omega}\delta^{\mathfrak{L}-\gamma-\zeta}(x)\leavevmode\nobreak\
dx\leq C\quad\text{if}\quad\mathfrak{L}>\gamma+\zeta-1.\end{split}$ (4.34)

where $\mathfrak{L}>\gamma+\zeta-1$ holds trivially.
Case 2: $\zeta+\gamma=1$
Since $\mathfrak{L}>0$, we can choose $\chi\in(0,\mathfrak{L})$ small enough
such that

$\max\left\\{\ln^{\zeta}\left(\frac{\mathcal{D}_{\Omega}}{\delta(x)}\right),\ln^{\mathfrak{L}-\gamma}\left(\frac{\mathcal{D}_{\Omega}}{\delta(x)}\right)\right\\}\leq
C(\mathcal{D}_{\Omega},\zeta,\mathfrak{L},\gamma)\delta^{-\chi}(x)\quad\text{for
all}\ x\in\Omega.$

Then, we have

$\begin{split}\mathfrak{G}(u_{n})&\leq\int_{\Omega}f_{n}(x)u_{n}^{\mathfrak{L}-\gamma}\leavevmode\nobreak\
dx\leq\mathcal{G}_{2}\int_{\Omega}\frac{u_{n}^{\mathfrak{L}-\gamma}}{\delta^{\zeta}(x)}\leavevmode\nobreak\
dx\lesssim\int_{\Omega}\ln^{\zeta}\left(\frac{\mathcal{D}_{\Omega}}{\delta(x)}\right)u_{n}^{\mathfrak{L}-1}\leavevmode\nobreak\
dx\\\
&\lesssim\left\\{\begin{array}[]{ll}\int_{\Omega}\ln^{\mathfrak{L}-\gamma}\left(\frac{\mathcal{D}_{\Omega}}{\delta(x)}\right)\delta^{\mathfrak{L}-1}(x)\leavevmode\nobreak\
dx&\text{ if }\ \mathfrak{L}\geq 1,\vspace{0.1cm}\\\
\int_{\Omega}\ln^{\zeta}\left(\frac{\mathcal{D}_{\Omega}}{\delta(x)}\right)\delta^{\mathfrak{L}-1}(x)\leavevmode\nobreak\
dx&\text{ if }\ \mathfrak{L}<1,\end{array}\right.\\\
&\lesssim\int_{\Omega}\delta^{\mathfrak{L}-1-\chi}(x)\leavevmode\nobreak\
dx\leq C\quad\text{if}\ \mathfrak{L}>\chi.\end{split}$ (4.35)

Case 3: $\zeta+\gamma>1$

$\begin{split}\mathfrak{G}(u_{n})&\leq\int_{\Omega}f_{n}(x)u_{n}^{\mathfrak{L}-\gamma}\leavevmode\nobreak\
dx\leq\mathcal{G}_{2}\int_{\Omega}\frac{u_{n}^{\mathfrak{L}-\gamma}}{\delta^{\zeta}(x)}\leavevmode\nobreak\
dx=\mathcal{G}_{2}\int_{\Omega}u_{n}^{\mathfrak{L}-\gamma}\left(\delta^{\frac{2-\zeta}{(\gamma+1)}}(x)\right)^{\frac{-\zeta(\gamma+1)}{2-\zeta}}\leavevmode\nobreak\
dx\\\
&\lesssim\int_{\Omega}u_{n}^{\mathfrak{L}-\gamma-\frac{\zeta(\gamma+1)}{2-\zeta}}\leavevmode\nobreak\
dx\lesssim\int_{\Omega}\delta^{\left(\mathfrak{L}-\gamma\right)\frac{2-\zeta}{(\gamma+1)}-\zeta}(x)\leavevmode\nobreak\
dx\leq
C\quad\text{if}\quad\mathfrak{L}>\frac{\gamma+\zeta-1}{2-\zeta}.\end{split}$
(4.36)

Collecting the estimates in (4.34)-(4.36), we obtain:

$\int_{\Omega}|\nabla
u_{n}^{\frac{\mathfrak{L}+1}{2}}|^{2}\leavevmode\nobreak\
dx\leq\frac{(\mathfrak{L}+1)^{2}}{4\mathfrak{L}}\mathfrak{G}(u_{n})\leq
C(\mathfrak{L},\gamma,\zeta)\quad\text{for all}\
\mathfrak{L}>\left\\{\begin{array}[]{ll}0&\text{ if }\ \zeta+\gamma\leq
1,\vspace{0.1cm}\\\ \frac{\gamma+\zeta-1}{2-\zeta}&\text{ if }\
\zeta+\gamma>1.\end{array}\right.$

∎

##  5 Existence, summability and Sobolev regularity results

###  5.1 Proof of Theorem 2.1

Let
$(r,\gamma)\in\mathcal{P}_{r,\gamma}\cap\\{(r,\gamma):r\in[1,r^{\sharp}),0\leq\gamma<1\\}$
and $u_{n}$ be the weak solution of the problem $(P_{n})$ in the sense that

$\begin{split}\int_{\Omega}\nabla u_{n}\cdot\nabla\psi\leavevmode\nobreak\
dx+\frac{C(N,s)}{2}\int_{\mathbb{R}^{N}}\int_{\mathbb{R}^{N}}&\frac{(u_{n}(x)-u_{n}(y))(\psi(x)-\psi(y))}{|x-y|^{N+2s}}\leavevmode\nobreak\
dx\leavevmode\nobreak\ dy\\\
&=\int_{\Omega}\frac{f_{n}(x)}{(u_{n}+\frac{1}{n})^{\gamma}}\psi\leavevmode\nobreak\
dx\quad\forall\psi\in H_{0}^{1}(\Omega).\end{split}$ (5.1)

From Lemma 4.1(iv), we know that $u_{n}$ is uniformly bounded in
$W_{0}^{1,q}(\Omega)$ with $q:=\frac{Nr(1+\gamma)}{N-r(1-\gamma)}.$ Then,
there exists a $u\in W_{0}^{1,q}(\Omega)$ such that

$\begin{split}&u_{n}\rightharpoonup u\ \text{in}\ W_{0}^{1,q}(\Omega),\ \
u_{n}\to u\ \text{in}\ L^{j}(\Omega)\ \ \text{for}\ 1\leq j<\sigma_{r}\
\text{and a.e. in}\ \mathbb{R}^{N}.\end{split}$ (5.2)

Since, for any $\psi\in W_{0}^{1,q^{\prime}}(\Omega)$,

$W_{0}^{1,q}(\Omega)\ni f\mapsto\int_{\Omega}\nabla
f\cdot\nabla\psi\leavevmode\nobreak\
dx+\frac{C(N,s)}{2}\int_{\mathbb{R}^{N}}\int_{\mathbb{R}^{N}}\frac{(f(x)-f(y))(\psi(x)-\psi(y))}{|x-y|^{N+2s}}\leavevmode\nobreak\
dx\leavevmode\nobreak\ dy$

is a bounded linear functional on $W_{0}^{1,q}(\Omega).$ Then for every
$\psi\in W_{0}^{1,q^{\prime}}(\Omega)$, we obtain

$\begin{split}\int_{\Omega}\nabla u_{n}&\cdot\nabla\psi\leavevmode\nobreak\
dx+\frac{C(N,s)}{2}\int_{\mathbb{R}^{N}}\int_{\mathbb{R}^{N}}\frac{(u_{n}(x)-u_{n}(y))(\psi(x)-\psi(y))}{|x-y|^{N+2s}}\leavevmode\nobreak\
dx\leavevmode\nobreak\ dy\\\ &\to\int_{\Omega}\nabla
u\cdot\nabla\psi\leavevmode\nobreak\
dx+\frac{C(N,s)}{2}\int_{\mathbb{R}^{N}}\int_{\mathbb{R}^{N}}\frac{(u(x)-u(y))(\psi(x)-\psi(y))}{|x-y|^{N+2s}}\leavevmode\nobreak\
dx\leavevmode\nobreak\ dy\quad\text{as}\ n\to\infty\end{split}$ (5.3)

and using the fact that $f_{n}\leq f$, $C_{1}(\omega)\leq u_{1}\leq u_{n}$
a.e. in $\omega\Subset\Omega$ and Lebesgue dominated convergence theorem we
get

$\begin{split}\int_{\Omega}\frac{f_{n}(x)}{(u_{n}+\frac{1}{n})^{\gamma}}\psi\leavevmode\nobreak\
dx\to\int_{\Omega}\frac{f(x)}{u^{\gamma}}\psi\leavevmode\nobreak\
dx\quad\text{as}\ n\to\infty\quad\text{for $\psi\in L^{r^{\prime}}(\Omega)$
with $\operatorname{supp}(\psi)\Subset\Omega$}.\end{split}$ (5.4)

Passing limit $n\to\infty$ in (5.1) and using (5.3)-(5.4), we obtain our
claim. ∎

###  5.2 Proof of Theorem 2.2

Let
$(r,\gamma)\in\mathcal{P}_{r,\gamma}\setminus\\{(r,\gamma):r\in[1,r^{\sharp}),0\leq\gamma<1\\}$
and $u_{n}$ be the weak solution of the problem $(P_{n})$ and satisfies (5.1).
From Lemma 4.1, we know that $u_{n}^{\frac{\mathfrak{S}_{r}+1}{r}}$ is
uniformly bounded in $H_{0}^{1}(\Omega).$ The condition
$(r,\gamma)\in\mathcal{P}_{r,\gamma}\setminus\\{(r,\gamma):r\in[1,r^{\sharp}),0\leq\gamma<1\\}$
implies $\mathfrak{S}_{r}\geq 1$. Together with the fact that for every
compact subset $\omega\Subset\Omega$ there exists $C=C(\omega)$ independent of
$n$ such that $0<C(\omega)\leq u_{n}(x)$ for $x\in\omega$, we get $u_{n}$ is
uniformly bounded in $H_{loc}^{1}(\Omega)$. Precisely,

$\begin{split}\int_{\omega}|\nabla u_{n}|^{2}\leavevmode\nobreak\ dx\leq
C^{-(\mathfrak{S}_{r}-1)}\int_{\omega}u_{n}^{(\mathfrak{S}_{r}-1)}|\nabla
u_{n}|^{2}\leavevmode\nobreak\
dx\leq\frac{4C^{-(\mathfrak{S}_{r}-1)}}{(\mathfrak{S}_{r}+1)}\int_{\omega}|\nabla
u_{n}^{\frac{\mathfrak{S}_{r}+1}{2}}|^{2}\leavevmode\nobreak\ dx\leq
C_{1}\end{split}$

where $C_{1}$ is independent of $n.$ Then, there exists a $u\in
H_{loc}^{1}(\Omega)$ such that

$\begin{split}&u_{n}\rightharpoonup u\ \text{in}\ H_{loc}^{1}(\Omega),\ \
u_{n}\to u\ \text{in}\ L^{j}_{loc}(\Omega)\ \text{for}\ 1\leq j<2^{\ast}\
\text{and a.e. in}\ \mathbb{R}^{N}.\end{split}$ (5.5)

Therefore, by using the weak convergence property and adopting the same
arguments from [21, Theorem 3.6], we are able to pass limits in the left hand
side of (5.1), i.e. for any $\psi\in H^{1}_{loc}(\Omega)$ with
$\operatorname{supp}(\psi)\Subset\Omega$,

$\begin{split}\int_{\Omega}\nabla u_{n}&\cdot\nabla\psi\leavevmode\nobreak\
dx+\frac{C(N,s)}{2}\int_{\mathbb{R}^{N}}\int_{\mathbb{R}^{N}}\frac{(u_{n}(x)-u_{n}(y))(\psi(x)-\psi(y))}{|x-y|^{N+2s}}\leavevmode\nobreak\
dx\leavevmode\nobreak\ dy\\\ &\to\int_{\Omega}\nabla
u\cdot\nabla\psi\leavevmode\nobreak\
dx+\frac{C(N,s)}{2}\int_{\mathbb{R}^{N}}\int_{\mathbb{R}^{N}}\frac{(u(x)-u(y))(\psi(x)-\psi(y))}{|x-y|^{N+2s}}\leavevmode\nobreak\
dx\leavevmode\nobreak\ dy\quad\text{as}\ n\to\infty.\end{split}$ (5.6)

Now, by repeating the same arguments as from the proof of Theorem 2.1 for the
limit passage in the right hand side term and by passing limits $n\to\infty$
using (5.6), we obtain our claim. ∎

###  5.3 Proof of Theorem 2.3

Let $u$ be the weak solution of the problem (1.2) obtained from the
approximable solution $u_{n}$ of the problem $(P_{n})$.

  1. (i)

Let $r\in[1,\frac{N}{2})$ and $N>2$. Then from Lemma 4.1, we know that
$u_{n}^{\frac{\mathfrak{S}_{r}+1}{2}}$ is uniformly bounded in
$H_{0}^{1}(\Omega)$ which further implies that there exists a $v\in
H_{0}^{1}(\Omega)$ such that
$u_{n}^{\frac{\mathfrak{S}_{r}+1}{2}}\rightharpoonup v\ \text{in}\
H_{0}^{1}(\Omega)$ and by Sobolev’s embedding, we get
$u_{n}^{\frac{\mathfrak{S}_{r}+1}{2}}\to v\ \text{in}\ L^{j}(\Omega)$ for
every $1\leq j<2^{\ast}$ and a.e. in $\mathbb{R}^{N}.$ Together with
convergence estimate in (5.2) and (5.5), we obtain
$v=u^{\frac{\mathfrak{S}_{r}+1}{2}}$ a.e. in $\mathbb{R}^{N}.$ The remaining
estimates follows from using the embedding $L^{r}(\Omega)\hookrightarrow
L^{j}(\Omega)$ and repeating the same proof by replacing $\mathfrak{S}_{r}$ by
$\mathfrak{S}_{j}.$

  2. (ii)

Let $r=\frac{N}{2}$ and $N\geq 2.$ Again from Lemma 4.1, we imply that
$\mathfrak{H}\left(\frac{u_{n}}{2}\right)$ in case of $0\leq\gamma\leq 1$ and
$\mathfrak{D}\left(\frac{u_{n}}{2}\right)$ in case of $\gamma>1$ are uniformly
bounded in $H_{0}^{1}(\Omega)$ which further implies that there exist
$v_{1},v_{2}\in H_{0}^{1}(\Omega)$ such that

$\begin{split}&0\leq\gamma\leq 1:\
\mathfrak{H}\left(\frac{u_{n}}{2}\right)\rightharpoonup v_{1}\ \text{in}\
H_{0}^{1}(\Omega)\ \text{and}\ \mathfrak{H}\left(\frac{u_{n}}{2}\right)\to
v_{1}\ \text{in}\ L^{j}(\Omega)\ \text{for every}\ 1\leq j<2^{\ast}\ \text{and
a.e. in}\ \mathbb{R}^{N},\\\ &\gamma>1:\
\mathfrak{D}\left(\frac{u_{n}}{2}\right)\rightharpoonup v_{2}\ \text{in}\
H_{0}^{1}(\Omega)\ \text{and}\ \mathfrak{D}\left(\frac{u_{n}}{2}\right)\to
v_{2}\ \text{in}\ L^{j}(\Omega)\ \text{ for every}\ 1\leq j<2^{\ast}\
\text{and a.e. in}\ \mathbb{R}^{N}.\end{split}$

Using the above estimate with the convergence properties in (5.2) and (5.5),
we identify the limit functions $v_{1}$ and $v_{2}$ as

$v_{1}=\mathfrak{H}\left(\frac{u}{2}\right)\quad\text{and}\quad
v_{2}=\mathfrak{D}\left(\frac{u}{2}\right)\quad\text{a.e. in}\
\mathbb{R}^{N}.$

  3. (iii)

From Lemma 4.1, we know that $u_{n}$ is uniformly bounded in
$L^{\infty}(\Omega)$ when $r>\frac{N}{2}$ and hence $u\in L^{\infty}(\Omega).$

  4. (iv)

The first part in the claim follows from the proof of Theorem 2.1 and the
second part of claim follows using the embedding $L^{r}(\Omega)\hookrightarrow
L^{1}(\Omega)$ for $1\leq r\leq\infty$ and by repeating the proof of Theorem
2.1 and 2.2 with $\mathfrak{S}_{r}=1$.

∎

###  5.4 Proof of Theorem 2.4

The convergence estimate in the proof of Theorem 2.1 and Theorem 2.2, and the
uniform a priori estimates in Lemma 4.2 implies the required claim. The only
if statement in proving optimal Sobolev regularity follows from the Hardy
inequality and the boundary behavior of the weak solution. Precisely, if

$\mathfrak{S}\leq 0\ \text{and}\ \gamma+\frac{1}{r}<1,\quad\text{then}\quad
u^{\frac{\mathfrak{S}+1}{2}}\notin H_{0}^{1}(\Omega).$

Indeed, in this case, we have

$\|u^{\frac{\mathfrak{S}+1}{2}}\|_{H_{0}^{1}(\Omega)}^{2}\geq
C\int_{\Omega}\frac{|u^{\frac{\mathfrak{S}+1}{2}}(x)|^{2}}{|d(x)|^{2}}\leavevmode\nobreak\
dx\geq C\int_{\Omega}d^{\mathfrak{S}-1}(x)\leavevmode\nobreak\ dx=\infty$

and hence, we deduce $u^{\frac{\mathfrak{S}+1}{2}}\notin H_{0}^{1}(\Omega).$ ∎

###  5.5 Proof of Theorem 2.5

Let $u$ and $v$ be two weak solutions of the problem (1.2)-(1.3) for datum $f$
and $g$ respectively, obtained in Theorem 2.1 and Theorem 2.2. Let
$\\{u_{n}\\}_{n\in\mathbb{N}},\\{v_{n}\\}_{n\in\mathbb{N}}\subset
H_{0}^{1}(\Omega)\cap l^{\infty}(\Omega)$ are sequence of the solutions of the
non-singular approximating problem $(P_{n})$ for $n\in\mathbb{N}$ where
$f_{n}$ and $g_{n}$ are increasing sequences such that $f_{n}\to f$ and
$g_{n}\to g$ in $L^{r}(\Omega).$ By taking
$[(u_{n}-v_{n})_{+}+\epsilon]^{\mathfrak{S}_{r}}-\epsilon^{\mathfrak{S}_{r}}$
as a test function in (3.5), we obtain

$\begin{split}&\mathfrak{S}_{r}\int_{\Omega}[\epsilon+(u_{n}-v_{n})_{+}]^{\mathfrak{S}_{r}-1}\nabla(u_{n}-v_{n})\cdot\nabla(u_{n}-v_{n})_{+}\leavevmode\nobreak\
dx+\int_{\Omega}\frac{g_{n}}{\left(v_{n}+\frac{1}{n}\right)^{\gamma}}([(u_{n}-v_{n})_{+}+\epsilon]^{\mathfrak{S}_{r}}-\epsilon^{\mathfrak{S}_{r}})\leavevmode\nobreak\
dx\\\
&+\frac{C(N,s)}{2}\int_{\mathbb{R}^{N}}\int_{\mathbb{R}^{N}}\frac{((u_{n}-v_{n})(x)-(u_{n}-v_{n})(y))([(u_{n}-v_{n})_{+}+\epsilon]^{\mathfrak{S}_{r}}(x)-[(u_{n}-v_{n})_{+}+\epsilon]^{\mathfrak{S}_{r}}(y))}{|x-y|^{N+2s}}\leavevmode\nobreak\
dx\leavevmode\nobreak\ dy\\\
&=\int_{\Omega}\frac{f_{n}}{\left(u_{n}+\frac{1}{n}\right)^{\gamma}}([(u_{n}-v_{n})_{+}+\epsilon]^{\mathfrak{S}_{r}}-\epsilon^{\mathfrak{S}_{r}})\leavevmode\nobreak\
dx\end{split}$

In order to pass limits in the integrals on the left hand side of the above
equality we use Fatou’s Lemma and for the integral on the right hand side, we
use Lebesgue dominated convergence, since for $\epsilon<\frac{1}{n}$, and
using (4.3) and (4.6), we have the dominating function

$\begin{split}\frac{f_{n}}{\left(u_{n}+\frac{1}{n}\right)^{\gamma}}([(u_{n}-v_{n})_{+}+\epsilon]^{\mathfrak{S}_{r}}-\epsilon^{\mathfrak{S}_{r}})\leq
f_{n}(u_{n}+1/n)^{\mathfrak{S}_{r}-\gamma}\leq
f(u+1)^{\mathfrak{S}_{r}-\gamma}\leq\frac{|f|^{r}}{r}+\frac{(r-1)}{r}|u+1|^{\sigma_{r}}.\end{split}$

Therefore, we have

$\begin{split}I_{1}+I_{2}&+I_{3}:=\mathfrak{S}_{r}\int_{\Omega}(u_{n}-v_{n})_{+}^{\mathfrak{S}_{r}-1}\nabla(u_{n}-v_{n})\cdot\nabla(u_{n}-v_{n})_{+}\leavevmode\nobreak\
dx+\int_{\Omega}\frac{g_{n}}{\left(v_{n}+\frac{1}{n}\right)^{\gamma}}(u_{n}-v_{n})_{+}^{\mathfrak{S}_{r}}\leavevmode\nobreak\
dx\\\
&+\frac{C(N,s)}{2}\int_{\mathbb{R}^{N}}\int_{\mathbb{R}^{N}}\frac{((u_{n}-v_{n})(x)-(u_{n}-v_{n})(y))((u_{n}-v_{n})_{+}^{\mathfrak{S}_{r}}(x)-(u_{n}-v_{n})_{+}^{\mathfrak{S}_{r}}(y))}{|x-y|^{N+2s}}\leavevmode\nobreak\
dx\leavevmode\nobreak\ dy\\\
&\leq\int_{\Omega}\frac{f_{n}}{\left(u_{n}+\frac{1}{n}\right)^{\gamma}}(u_{n}-v_{n})_{+}^{\mathfrak{S}_{r}}\leavevmode\nobreak\
dx:=I_{4}\end{split}$ (5.7)

Now, we separately estimate the integrals in the above inequality.
Estimate for $I_{4}$: With the choice of $\mathfrak{S}_{r}$ in (4.3) and
(4.6),

$0\leq\frac{f_{n}}{\left(u_{n}+\frac{1}{n}\right)^{\gamma}}(u_{n}-v_{n})_{+}^{\mathfrak{S}_{r}}\leq
fu^{\mathfrak{S}_{r}-\gamma}\leq\frac{|f|^{r}}{r}+\frac{(r-1)}{r}|u|^{\sigma_{r}},$
(5.8)

and the fact that $u_{n}\to u$ and $v_{n}\to v$ a.e. in $\Omega$, Lebesgue
dominated convergence theorem implies

$I_{4}\to\int_{\Omega}\frac{f}{u^{\gamma}}(u-v)_{+}^{\mathfrak{S}_{r}}\leavevmode\nobreak\
dx$ (5.9)

Estimate for $I_{2}+I_{3}$: It is easy to see that

$((u_{n}-v_{n})(x)-(u_{n}-v_{n})(y))((u_{n}-v_{n})_{+}^{\mathfrak{S}_{r}}(x)-(u_{n}-v_{n})_{+}^{\mathfrak{S}_{r}}(y))\geq
0\quad\text{for}\quad(x,y)\in\mathbb{R}^{N}\times\mathbb{R}^{N}.$ (5.10)

Now by using the Fatou’s lemma and (5.10) we obtain,

$\int_{\Omega}\frac{g}{v^{\gamma}}(u-v)_{+}^{\mathfrak{S}_{r}}\leavevmode\nobreak\
dx\leq\lim\inf_{n\to\infty}I_{2}+I_{3}$ (5.11)

Estimate for $I_{1}$: For this we divide the proof into two cases:
$\mathfrak{S}_{r}\leq 1$ and $\mathfrak{S}_{r}>1.$ In the first case, for
every $\delta>0$, we have

$\begin{split}\frac{4\mathfrak{S}_{r}}{(\mathfrak{S}_{r}+1)^{2}}\int_{\Omega}\left|\nabla\left[(u_{n}-v_{n})_{+}+\delta\right]^{\frac{\mathfrak{S}_{r}+1}{2}}\right|^{2}\leavevmode\nobreak\
dx=\mathfrak{S}_{r}\int_{\Omega}\frac{\left|\nabla(u_{n}-v_{n})_{+}\right|^{2}}{\left((u_{n}-v_{n})_{+}+\delta\right)^{1-\mathfrak{S}_{r}}}\leavevmode\nobreak\
dx\leq I_{1}.\end{split}$

Together with the estimate in (5.8) and (5.10), we obtain
$\left[(u_{n}-v_{n})_{+}+\delta\right]^{\frac{\mathfrak{S}_{r}+1}{2}}$ is
uniformly bounded in $L^{2}$ and which combined with the a.e. convergence
$u_{n}-v_{n}\to u-v$ in $H_{0}^{1}(\Omega)$, we obtain

$\left[(u_{n}-v_{n})_{+}+\delta\right]^{\frac{\mathfrak{S}_{r}+1}{2}}\rightharpoonup\left[(u-v)_{+}+\delta\right]^{\frac{\mathfrak{S}_{r}+1}{2}}\
\text{in}\ H_{0}^{1}(\Omega)$

and the weak lower semicontinuity of the norm gives

$\begin{split}\mathfrak{S}_{r}\int_{\Omega}\frac{\left|\nabla(u-v)_{+}\right|^{2}}{\left((u-v)_{+}+\delta\right)^{1-\mathfrak{S}_{r}}}\leavevmode\nobreak\
dx&=\frac{4\mathfrak{S}_{r}}{(\mathfrak{S}_{r}+1)^{2}}\int_{\Omega}\left|\nabla\left[(u-v)_{+}+\delta\right]^{\frac{\mathfrak{S}_{r}+1}{2}}\right|^{2}\leavevmode\nobreak\
dx\\\
&\leq\lim\inf_{n\to\infty}\frac{4\mathfrak{S}_{r}}{(\mathfrak{S}_{r}+1)^{2}}\int_{\Omega}\left|\nabla\left[(u_{n}-v_{n})_{+}+\delta\right]^{\frac{\mathfrak{S}_{r}+1}{2}}\right|^{2}\leavevmode\nobreak\
dx\leq\lim\inf_{n\to\infty}I_{1}.\end{split}$

Using Monotone convergence theorem by taking $\delta$ decreasing to $0$, we
get

$\begin{split}\mathfrak{S}_{r}\int_{\Omega}\frac{\left|\nabla(u-v)_{+}\right|^{2}}{(u-v)_{+}^{1-\mathfrak{S}_{r}}}\leavevmode\nobreak\
dx\leq\lim\inf_{n\to\infty}I_{1}.\end{split}$ (5.12)

In the second case, $\mathfrak{S}_{r}>1$, taking $\delta=0$ and then reasoning
as above, we get the same estimate (5.12). Now, by passing $n\to\infty$ in
(5.7) and using the above convergence estimates in (5.9), (5.11) and (5.12),
we obtain

$\begin{split}\frac{4\mathfrak{S}_{r}}{(\mathfrak{S}_{r}+1)^{2}}\int_{\Omega}|\nabla(u-v)_{+}^{\frac{\mathfrak{S}_{r}+1}{2}}|^{2}\leavevmode\nobreak\
dx&\leq\int_{\Omega}\left(\frac{f(x)}{u^{\gamma}}-\frac{g(x)}{v^{\gamma}}\right)(u-v)_{+}^{\mathfrak{S}_{r}}\leavevmode\nobreak\
dx\leq\int_{\Omega}(f(x)-g(x))\frac{(u-v)_{+}^{\mathfrak{S}_{r}}}{u^{\gamma}}\leavevmode\nobreak\
dx\\\
&\leq\int_{\Omega}(f(x)-g(x))(u-v)_{+}^{\mathfrak{S}_{r}-\gamma}\leavevmode\nobreak\
dx.\end{split}$ (5.13)

Now, by using Hölder inequality, Sobolev embeddings and using the relation of
$\mathfrak{S}_{r}$ in (4.3) and (4.6), we obtain

$\begin{split}\frac{4\mathfrak{S}_{r}}{(\mathfrak{S}_{r}+1)^{2}}&\int_{\Omega}|\nabla(u-v)_{+}^{\frac{\mathfrak{S}_{r}+1}{2}}|^{2}\leavevmode\nobreak\
dx\leq\|f-g\|_{L^{r}(\Omega_{1})}\left(\int_{\Omega}\left((u-v)_{+}^{\frac{(\mathfrak{S}_{r}+1)}{2}}\right)^{\frac{2N}{N-2}}\leavevmode\nobreak\
dx\right)^{\frac{1}{r^{\prime}}}\\\
&\leq\|f-g\|_{L^{r}(\Omega_{1})}\left(S(N)\int_{\Omega}|\nabla(u-v)_{+}^{\frac{\mathfrak{S}_{r}+1}{2}}|^{2}\leavevmode\nobreak\
dx\right)^{\frac{N}{r^{\prime}(N-2)}}\ \text{where}\ \Omega_{1}:=\\{u\geq
v\\}\end{split}$

which further implies

$\begin{split}\|\nabla(u-v)_{+}^{\frac{\mathfrak{S}_{r}+1}{2}}\|^{2}_{L^{2}(\Omega)}&\leq
C\|f-g\|_{L^{r}(\Omega_{1})}^{\frac{r(N-2)}{N-2r}}\quad\text{where}\quad
C=\left(\frac{(\mathfrak{S}_{r}+1)^{2}}{4\mathfrak{S}_{r}}\right)^{\frac{r(N-2)}{N-2r}}\left(S(N)\right)^{\frac{N(r-1)}{N-2r}}\end{split}$

and analogously

$\begin{split}\|\nabla(v-u)_{+}^{\frac{\mathfrak{S}_{r}+1}{2}}\|^{2}_{L^{2}(\Omega)}&\leq
C\|f-g\|_{L^{r}(\Omega\setminus\Omega_{1})}^{\frac{r(N-2)}{N-2r}}.\end{split}$

∎
Proof of Corollary 2.1 The proof follows from the inequality relation (5.13)
by taking $f\leq g.$ ∎

###  5.6 Proof of Theorem 2.6:

Let $\zeta+\gamma\leq 1$ and $u_{n}$ be the weak solution of the problem
$(P_{n})$ in the sense that

$\begin{split}\int_{\Omega}\nabla u_{n}\cdot\nabla\psi\leavevmode\nobreak\
dx+\frac{C(N,s)}{2}\int_{\mathbb{R}^{N}}\int_{\mathbb{R}^{N}}&\frac{(u_{n}(x)-u_{n}(y))(\psi(x)-\psi(y))}{|x-y|^{N+2s}}\leavevmode\nobreak\
dx\leavevmode\nobreak\ dy\\\
&=\int_{\Omega}\frac{f_{n}(x)}{(u_{n}+\frac{1}{n})^{\gamma}}\psi\leavevmode\nobreak\
dx\quad\forall\psi\in H_{0}^{1}(\Omega).\end{split}$ (5.14)

From Lemma 4.7 and Lemma 4.6, we know that, for any $\mathfrak{L}>0$,
$u_{n}^{\frac{\mathfrak{L}+1}{2}}$ is uniformly bounded in $H_{0}^{1}(\Omega)$
and $u_{n}\geq C\delta(x)$ for $x\in\Omega$ and $C>0$ independent of $n$.
Therefore, by taking $\mathfrak{L}=1$, we obtain

$\begin{split}&u_{n}\rightharpoonup u\ \text{in}\ H_{0}^{1}(\Omega),\ \
u_{n}\to u\ \text{in}\ L^{j}(\Omega)\ \text{for}\ 1\leq j<2^{\ast}\ \text{and
a.e. in}\ \mathbb{R}^{N}.\end{split}$ (5.15)

Using Hardy’s inequality, for any $\psi\in H_{0}^{1}(\Omega)$,
$\frac{f_{n}\psi}{(u_{n}+\frac{1}{n})^{\gamma}}$ is uniformly integrable.
Indeed,

$\begin{split}\int_{\Omega}\frac{f_{n}(x)}{(u_{n}+\frac{1}{n})^{\gamma}}\psi\leavevmode\nobreak\
dx\leq\int_{\Omega}\frac{f(x)}{u_{n}^{\gamma}}\psi\leavevmode\nobreak\
dx\lesssim\int_{\Omega}\delta^{1-\gamma-\zeta}(x)\frac{\psi}{\delta(x)}\leavevmode\nobreak\
dx\lesssim\|\frac{\psi}{\delta}\|_{L^{2}(\Omega)}\leq\|\psi\|_{H_{0}^{1}(\Omega)}.\end{split}$

Finally, by using Vitali convergence theorem and convergence properties in
(5.15), we are able to pass limits in (5.14) and obtain

$\begin{split}\int_{\Omega}\nabla u\cdot\nabla\psi\leavevmode\nobreak\
dx+\frac{C(N,s)}{2}\int_{\mathbb{R}^{N}}\int_{\mathbb{R}^{N}}&\frac{(u(x)-u(y))(\psi(x)-\psi(y))}{|x-y|^{N+2s}}\leavevmode\nobreak\
dx\leavevmode\nobreak\
dy=\int_{\Omega}\frac{f(x)}{u^{\gamma}}\psi\leavevmode\nobreak\
dx\quad\forall\psi\in H_{0}^{1}(\Omega).\end{split}$

Finally, by comparison principle (see [21, Theorem 4.2]) for any
$n\in\mathbb{N}$, $u_{n}\leq v$ a.e. in $\Omega$ where $v$ is another weak
solution of (1.2)-(4.6). Passing to the limit $n\to\infty$ gives that u is a
minimal solution. ∎

###  5.7 Proof of Theorem 2.7

Let $\zeta+\gamma>1$ and $u_{n}$ be the weak solution of the problem
$(P_{n})$. From Lemma 4.7 and Lemma 3.2, we have, for any
$\mathfrak{L}>\mathfrak{L}^{*}$, $u_{n}^{\frac{\mathfrak{L}+1}{2}}$ is
uniformly bounded in $H_{0}^{1}(\Omega)$ and $u_{n}\geq C(\omega)$ for
$x\in\omega\Subset\Omega$ and $C>0$ independent of $n$. Now we divide the
proof into two cases:
Case 1:  $\mathfrak{L}^{*}>1$
In this case, $u_{n}$ is uniformly bounded in $H_{loc}^{1}(\Omega)$.
Precisely,

$\begin{split}\int_{\omega}|\nabla u_{n}|^{2}\leavevmode\nobreak\ dx\leq
C^{-(\mathfrak{L}-1)}(\omega)\int_{\omega}u_{n}^{(\mathfrak{L}-1)}|\nabla
u_{n}|^{2}\leavevmode\nobreak\
dx\leq\frac{4C^{-(\mathfrak{L}-1)}}{(\mathfrak{L}+1)}\int_{\omega}|\nabla
u_{n}^{\frac{\mathfrak{L}+1}{2}}|^{2}\leavevmode\nobreak\ dx\leq
C_{0}\end{split}$ (5.16)

where $C_{0}$ is independent of $n.$ Then, there exists a $u\in
H_{loc}^{1}(\Omega)$ such that

$\begin{split}&u_{n}\rightharpoonup u\ \text{in}\ H_{loc}^{1}(\Omega),\ \
u_{n}\to u\ \text{in}\ L^{j}_{loc}(\Omega)\ \text{for}\ 1\leq j<2^{\ast}\
\text{and a.e. in}\ \mathbb{R}^{N}.\end{split}$ (5.17)

Using Hardy’s inequality we have for any $\psi\in H_{0}^{1}(\Omega)$ with
$\operatorname{supp}(\psi)\Subset\Omega$,
$\frac{f_{n}\psi}{(u_{n}+\frac{1}{n})^{\gamma}}$ is uniformly integrable in
$L^{1}(\omega)$. Then, by using Vitali convergence theorem, convergence
properties in (5.17) and adopting the same arguments from [21, Theorem 3.6],
we are able to pass limits in (5.14), i.e. for any
$\psi\in\bigcup_{\omega\Subset\Omega}H_{0}^{1}(\omega)$,

$\begin{split}\int_{\Omega}\nabla u\cdot\nabla\psi\leavevmode\nobreak\
dx+\frac{C(N,s)}{2}\int_{\mathbb{R}^{N}}\int_{\mathbb{R}^{N}}&\frac{(u(x)-u(y))(\psi(x)-\psi(y))}{|x-y|^{N+2s}}\leavevmode\nobreak\
dx\leavevmode\nobreak\
dy=\int_{\Omega}\frac{f(x)}{u^{\gamma}}\psi\leavevmode\nobreak\
dx.\end{split}$

Case 2:  $\mathfrak{L}^{*}\leq 1$
In this case, $u_{n}$ is uniformly bounded in $H_{0}^{1}(\Omega)$, precisely,
by taking $\mathfrak{L}=1$. Now, by repeating the proof of Theorem 2.6, we
obtain our claim.

###  5.8 Proof of Theorem 2.8

From Theorem 2.6, Theorem 2.7, Lemma 4.6 and using the fact that
$u:=\lim_{n\to\infty}u_{n}$ is a weak solution of the main problem (1.2), we
obtain that $u$ satisfies

$C_{1}\delta(x)\leq u(x)\leq C_{2}\delta(x)\left\\{\begin{array}[]{ll}1&\text{
if }\ \zeta+\gamma<1,\vspace{0.1cm}\\\
\ln\left(\frac{\mathcal{D}_{\Omega}}{\delta(x)}\right)&\text{ if }\
\zeta+\gamma=1,\end{array}\right.$

and

$C^{1}\delta^{\frac{2-\zeta}{\gamma+1}}(x)\leq u(x)\leq
C^{2}\delta^{\frac{2-\zeta}{\gamma+1}}(x).$

Now, it only remains to prove the optimal boundary behavior in case of
$\zeta+\gamma=1.$ Let $w_{\Xi}$ be the solution to the problem

$\left\\{\begin{aligned} -\Delta
w_{\Xi}+(-\Delta)^{s}w_{\Xi}&=\frac{1}{\delta}\ln^{-\Xi}\left(\frac{\mathcal{D}_{\Omega}}{\delta}\right),\quad
w_{\Xi}>0\quad&&\text{ in }\Omega,\\\ w_{\Xi}&=0&&\text{ in }\
\mathbb{R}^{N}\setminus\Omega.\end{aligned}\right.$ (5.18)

where $\Xi\in[0,1).$ Then, by using the integral representation of the
solution, Lemma 4.3 and Lemma 4.4, we obtain

$\displaystyle
w_{\Xi}(x):=\left\\{\begin{array}[]{ll}\mathbb{G}^{\Omega}\left[\frac{1}{\delta}\ln^{-\Xi}(\frac{\mathcal{D}_{\Omega}}{\delta})\right](x)&\text{
if }x\in\Omega,\\\ 0&\text{ if }x\in\mathbb{R}^{N}\setminus\Omega.\\\
\end{array}\right.$ (5.19)

and there exists $c_{1},c_{2}>0$ such that

$c_{1}\delta(x)\ell^{1-\Xi}(\delta(x))\leq w_{\Xi}(x)\leq
c_{2}\delta(x)\ell^{1-\Xi}(\delta(x))\quad\forall\ x\in\Omega.$ (5.20)

Moreover, for any $\Xi_{0}<1$, $c_{1}$ and $c_{2}$ are uniform for any
$0\leq\Xi\leq\Xi_{0}.$ Now we divide the proof into two cases:
Case 1: $\zeta>0$
Since $u$ satisfies (4.25), $\gamma+\zeta=1$ and $f\in\mathcal{A}_{\zeta}$,
there exists $d_{0},d_{1}>0$ such that $d_{0}\delta^{-\zeta}(x)\leq f(x)\leq
d_{1}\delta^{-\zeta}(x)$ for $x\in\Omega$ and $u$ satisfies

$(-\Delta)(d_{0}w_{\gamma})+(-\Delta)^{s}(d_{0}w_{\gamma})=d_{0}\delta^{-1}(x)\ln^{-\gamma}\left(\frac{\mathcal{D}_{\Omega}}{\delta(x)}\right)\leq\frac{f(x)}{u^{\gamma}}=(-\Delta)u+(-\Delta)^{s}u.$

Now, by using the comparison principle and (5.20) with $\Xi_{1}:=\gamma$, we
obtain

$c_{1}d_{0}\delta(x)\ln^{1-\Xi_{1}}\left(\frac{\mathcal{D}_{\Omega}}{\delta(x)}\right)\leq
d_{0}w_{\Xi_{1}}(x)\leq u(x)\ \text{for}\ x\in\Omega$

and

$(-\Delta)u+(-\Delta)^{s}u=\frac{f(x)}{u^{\gamma}}\leq\frac{d_{1}}{(c_{1}d_{0})^{\gamma}}\delta^{-1}(x)\ln^{-\Xi_{2}}\left(\frac{\mathcal{D}_{\Omega}}{\delta(x)}\right)=(-\Delta)\left(\frac{d_{1}}{(c_{1}d_{0})^{\gamma}}w_{\Xi_{2}}\right)+(-\Delta)^{s}\left(\frac{d_{1}}{(c_{1}d_{0})^{\gamma}}w_{\Xi_{2}}\right)$

where $\Xi_{2}:=\gamma(1-\Xi_{1}).$ Again, using the comparison principle, we
obtain

$u(x)\leq\frac{d_{1}}{(c_{1}d_{0})^{\gamma}}w_{\Xi_{2}}(x)\leq\frac{d_{1}c_{2}}{(c_{1}d_{0})^{\gamma}}\delta(x)\ln^{1-\Xi_{2}}\left(\frac{\mathcal{D}_{\Omega}}{\delta(x)}\right)\
\text{for}\ x\in\Omega$

Iterating these estimates, we obtain for any $j\in\mathbb{N}$

$\begin{split}\frac{(c_{1}d_{0})^{1+\gamma^{2}+\dots+\gamma^{2j}}}{(d_{1}c_{2})^{\gamma+\gamma^{3}+\dots+\gamma^{2j-1}}}\delta(x)&\ln^{1-\gamma+\gamma^{2}-\gamma^{3}\dots-\gamma^{2j+1}}\left(\frac{\mathcal{D}_{\Omega}}{\delta(x)}\right)\leq
u(x)\\\
&\leq\frac{(c_{1}d_{0})^{1+\gamma^{2}+\dots+\gamma^{2j}}}{(d_{1}c_{2})^{\gamma+\gamma^{3}+\dots+\gamma^{2j+1}}}\delta(x)\ln^{1-\gamma+\gamma^{2}\dots+\gamma^{2j+2}}\left(\frac{\mathcal{D}_{\Omega}}{\delta(x)}\right)\end{split}$

Passing to the limit $j\to\infty$, we obtain

$d_{2}\delta(x)\ln^{\frac{1}{1+\gamma}}\left(\frac{\mathcal{D}_{\Omega}}{\delta(x)}\right)\leq
u(x)\leq
d_{3}\delta(x)\ln^{\frac{1}{1+\gamma}}\left(\frac{\mathcal{D}_{\Omega}}{\delta(x)}\right)$

where $d_{2},d_{3}>0$ are constants depending upon $\gamma,d_{0},d_{1},c_{1}$
and $c_{2}.$
Case 2: $\zeta=0$
In this case, by taking $\Xi=\frac{1}{2}$ in (5.18) and (5.19), there exists
$c_{3},c_{4}>0$ such that

$c_{3}\delta(x)\ln^{\frac{1}{2}}\left(\frac{\mathcal{D}_{\Omega}}{\delta(x)}\right)\leq
w_{\frac{1}{2}}(x)\leq
c_{4}\delta(x)\ln^{\frac{1}{2}}\left(\frac{\mathcal{D}_{\Omega}}{\delta(x)}\right)\quad\forall\
x\in\Omega.$ (5.21)

and for a positive constant $C$ large enough, we have

$(-\Delta)(Cw_{\frac{1}{2}})+(-\Delta)^{s}(Cw_{\frac{1}{2}})\geq\frac{1}{Cw_{\frac{1}{2}}}\
\text{and}\
(-\Delta)\left(\frac{w_{\frac{1}{2}}}{C}\right)+(-\Delta)^{s}\left(\frac{w_{\frac{1}{2}}}{C}\right)\leq\frac{C}{w_{\frac{1}{2}}}.$

Therefore again by comparison principle, we obtain

$\frac{c_{3}}{C}\delta(x)\ln^{\frac{1}{2}}\left(\frac{\mathcal{D}_{\Omega}}{\delta(x)}\right)\leq
u(x)\leq\frac{c_{4}}{C}\delta(x)\ln^{\frac{1}{2}}\left(\frac{\mathcal{D}_{\Omega}}{\delta(x)}\right)\quad\forall\
x\in\Omega.$

The only if statement in proving optimal Sobolev regularity follows from the
Hardy inequality and the boundary behavior of the weak solution. Precisely, if

$\mathfrak{L}\leq\left\\{\begin{array}[]{ll}0&\text{ if }\ \zeta+\gamma\leq
1,\vspace{0.1cm}\\\ \mathfrak{L}^{*}&\text{ if }\
\zeta+\gamma>1.\end{array}\right.\quad\text{then}\quad
u^{\frac{\mathfrak{L}+1}{2}}\notin H_{0}^{1}(\Omega).$

Indeed, in case of $\zeta+\gamma\leq 1$ and $\mathfrak{L}\leq 0$, we have

$\|u^{\frac{\mathfrak{L}+1}{2}}\|_{H_{0}^{1}(\Omega)}^{2}\geq
C\int_{\Omega}\frac{|u^{\frac{\mathfrak{L}+1}{2}}(x)|^{2}}{|d(x)|^{2}}\leavevmode\nobreak\
dx\geq C\int_{\Omega}d^{\mathfrak{L}-1}(x)\leavevmode\nobreak\ dx=\infty.$

In the same way, if $\zeta+\gamma>1$ and $\mathfrak{L}\leq\mathfrak{L}^{*}$,
then

$\|u^{\frac{\mathfrak{L}+1}{2}}\|_{H_{0}^{1}(\Omega)}^{2}\geq
C\int_{\Omega}\frac{|u^{\frac{\mathfrak{L}+1}{2}}(x)|^{2}}{|d(x)|^{2}}\leavevmode\nobreak\
dx\geq
C\int_{\Omega}d^{\frac{(\mathfrak{L}+1)(2-\zeta)}{(1+\gamma)}-2}(x)\leavevmode\nobreak\
dx=\infty$

and we deduce $u^{\frac{\mathfrak{L}+1}{2}}\notin H_{0}^{1}(\Omega).$ ∎

###  5.9 Proof of Theorem 2.9

Let $\zeta\geq 2$ and $f\in\mathcal{A}_{\zeta}$. We choose $\Gamma\in(0,1)$
and $\zeta_{0}<2$ such that $g(x)\leq f(x)$ for $g\in\mathcal{A}_{\zeta_{0}}$
and the constant $\Gamma$ is independent of $\zeta_{0}$ for
$\zeta_{0}\geq\zeta_{0}^{*}$ with $\zeta_{0}^{*}>0$. To prove our claim, we
proceed by contradiction. Assume there exist a weak solution $w\in
H_{loc}^{1}(\Omega)$ of the problem (1.2) and $\mathfrak{L}_{0}\geq 1$ such
that $w^{\mathfrak{L}_{0}}\in H_{0}^{1}(\Omega)$.
For $n\in\mathbb{N}$, let $v_{n}\in H_{0}^{1}(\Omega)\cap
C^{0,\ell}(\overline{\Omega})$ be the unique weak solution of

$\int_{\Omega}\nabla v_{n}\cdot\nabla\phi\leavevmode\nobreak\
dx+\frac{C(N,s)}{2}\iint_{\mathbb{R}^{2N}}\frac{(v_{n}(x)-v_{n}(y))(\phi(x)-\phi(y))}{|x-y|^{N+2s}}\leavevmode\nobreak\
dx\leavevmode\nobreak\ dy=\int_{\Omega}\frac{\Gamma
g_{n}(x)}{(v_{n}+\frac{1}{n})^{\gamma}}\phi\leavevmode\nobreak\ dx$ (5.22)

for any $\phi\in H_{0}^{1}(\Omega)$. By the continuity of $v_{n}$, for given
$\theta>0$, there exists a $\eta{=\eta(n,\theta)}>0$ such that
$v_{n}\leq\frac{\theta}{2}$ in $\Omega_{\eta}$. Since $w\geq 0$, then
$u:=v_{n}-w-\theta\leq-\frac{\theta}{2}<0$ in $\Omega_{\eta}$ and

$\mbox{supp}({u}^{+})\subset\mbox{supp}((v_{n}-\theta)^{+})\subset\Omega\setminus\Omega_{\eta}.$

We have $u^{+}\in H_{0}^{1}(\tilde{\Omega})\subset H_{0}^{1}(\Omega)$ for some
$\tilde{\Omega}$ such that
$\Omega\setminus\Omega_{\eta}\subset\tilde{\Omega}\Subset\Omega$. Hence,
choosing $u^{+}$ as a test function in (5.22), we get

$\begin{split}\int_{\Omega}\nabla v_{n}\cdot\nabla u^{+}\leavevmode\nobreak\
dx&+\frac{C(N,s)}{2}\iint_{\mathbb{R}^{2N}}\frac{(v_{n}(x)-v_{n}(y))(u^{+}(x)-u^{+}(y))}{|x-y|^{N+2s}}\leavevmode\nobreak\
dx\leavevmode\nobreak\ dy\\\ &=\int_{\Omega}\frac{\Gamma
g_{n}(x)}{(v_{n}+\epsilon)^{\gamma}}u^{+}\leavevmode\nobreak\
dx\leq\int_{\Omega}\frac{\Gamma
g_{n}(x)}{v_{n}^{\gamma}}u^{+}\leavevmode\nobreak\ dx.\end{split}$ (5.23)

Moreover, $w$ is a weak solution of $(P)$ and taking $u^{+}\in
H_{0}^{1}(\tilde{\Omega})$ as test function, we have

$\begin{split}\int_{\Omega}\nabla w\cdot\nabla u^{+}\leavevmode\nobreak\
dx&+\frac{C(N,s)}{2}\iint_{\mathbb{R}^{2N}}\frac{(w(x)-w(y))(u^{+}(x)-u^{+}(y))}{|x-y|^{N+2s}}\leavevmode\nobreak\
dx\leavevmode\nobreak\ dy\\\
&=\int_{\Omega}\frac{f(x)}{w^{\gamma}}w^{+}\leavevmode\nobreak\
dx\geq\int_{\Omega}\frac{\Gamma g_{n}(x)}{w^{\gamma}}w^{+}\leavevmode\nobreak\
dx.\end{split}$ (5.24)

By subtracting (5.24) and (5.23), we get

$\begin{split}0\leq\int_{\Omega}(\nabla v_{n}-\nabla w)\cdot\nabla
u^{+}\leavevmode\nobreak\
dx&+\frac{C(N,s)}{2}\iint_{\mathbb{R}^{2N}}\frac{((v_{n}(x)-v_{n}(y))-(w(x)-w(y)))(u^{+}(x)-u^{+}(y))}{|x-y|^{N+2s}}\leavevmode\nobreak\
dx\leavevmode\nobreak\ dy\\\ &\leq\int_{\Omega}\left(\frac{\Gamma
g_{n}(x)}{v_{n}^{\gamma}}-\frac{\Gamma
g_{n}(x)}{w^{\gamma}}\right)u^{+}\leavevmode\nobreak\ dx\leq 0.\end{split}$

which further implies $w^{+}=(v_{n}-w-\theta)^{+}=0$ a.e. in $\Omega$. Since
$\theta$ is arbitrary, we deduce $v_{n}\leq w$ in $\Omega.$ Using Lemma 4.6,
we obtain

$c_{1}(\delta(x)+n^{\frac{-(1+\gamma)}{(2-\zeta_{0})}})^{\frac{2-\zeta_{0}}{\gamma+1}}-c_{2}n^{-1}\leq
v_{n}\leq w\ \text{in}\ \Omega.$

Now, by using the Hardy inequality and $w^{\mathfrak{L}_{0}}\in
H_{0}^{1}(\Omega)$, we obtain

$\int_{\Omega}\left(c_{1}(\delta(x)+n^{\frac{-(1+\gamma)}{(2-\zeta_{0})}})^{\frac{2-\zeta_{0}}{\gamma+1}}-c_{2}n^{-1}\right)^{2\mathfrak{L}_{0}}d^{-2}(x)\leavevmode\nobreak\
dx\leq\int_{\Omega}\left|\frac{w^{\mathfrak{L}_{0}}}{d(x)}\right|^{2}\leavevmode\nobreak\
dx<\infty.$

Now, by choosing $\zeta_{0}$ close enough to $2$ and by taking $n\to\infty$,
we obtain that the left hand side is not finite, which is a contradiction and
hence claim. ∎

##  6 Appendix

In this part, we recall some preliminary results comprised of upper and lower
estimates of Green kernel, lower Hopf type estimate, and action of the Green
kernel on the distance power functions.

######  Proposition 6.1 (Theorem 1, [26]).

Let $\Omega$ be a $C^{1,1}$ open set in $\mathbb{R}^{N}$ and
$\mathcal{G}^{\Omega}$ denote the Green kernel for the mixed operator
$(-\Delta)+(\Delta)^{s}$ defined on
$\operatorname{dom}(\mathcal{G}^{\Omega}):=\Omega\times\Omega\setminus\\{(x,x):x\in\Omega\\}$.
Then, for any $(x,y)\in\operatorname{dom}(\mathcal{G}^{\Omega}):$

$\mathcal{G}^{\Omega}(x,y)\asymp\frac{1}{|x-y|^{N-2}}\left(\frac{\delta(y)\delta(x)}{|x-y|^{2}}\wedge
1\right)\ \text{ if }\ N\geq 3.$

######  Lemma 6.1 (Theorem 2.6, [2]).

There exists $C>0$ such that for any $f\geq 0$,

$\mathbb{G}^{\Omega}\left[f\right](x):=\int_{\Omega}\mathcal{G}^{\Omega}(x,y)f(y)\leavevmode\nobreak\
dy\geq
C\delta(x)^{\gamma}\int_{\Omega}\delta^{\gamma}(y)f(y)\leavevmode\nobreak\
dy,\ x\in\Omega.$

######  Lemma 6.2 (Theorem 3.4, [2]).