Combined effects in mixed local–nonlocal stationary problems

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 nonexistence properties, power- and exponential-type Sobolev regularity results, and the boundary behaviour of the weak solution, in the light of the interplay between the summability of the datum and the power exponent in singular nonlinearities.

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

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

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 as the sum 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}}\,{\rm d}y. \]

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}}\,{\rm 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 phenomena and motivation behind studying problems involving such nonlocal operators, we refer the reader to review the famous Hitchhiker's guide [46] and references 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 nonlocal 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 thorough discussion, see [39] for mixed dispersal movement strategy, [42] for nonlocal diffusion strategy, [25] for conditional dispersal strategy). These type of operators also arises in the models obtained from the superposition of two different scaled stochastic processes. For a detailed presentation on this, we refer the reader to [43].

Very recently, a great amount of attention has been paid to studying elliptic problems involving a 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, 36, 37], symmetry results [13], Faber–Krahn-type inequality [12], Neumann problems [43], Green functions estimates [26, 27] have been answered.

The study of elliptic or integral equations involving singular terms started in the early 1960s with 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, 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 ${f(x)}/{u^\gamma }$, and the temperature distribution in the conducting medium satisfies the local counterpart of the equation (see (1.2)). For interested readers, we refer to [30, 44, 48] 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 a weight function $f$ and singular nonlinearities:

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

subject to the homogeneous Dirichlet boundary conditions:

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

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 problems (1.2)–(1.3). An elementary way to define the notion of a solution of 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:

   (1.4)\begin{align} & \int_{\Omega} \nabla u \cdot \nabla \psi\,{\rm d}x + \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}}\,{\rm d}x\,{\rm d}y\nonumber\\ & \quad = \int_{\Omega} \frac{f(x)}{u^\gamma}\psi\,{\rm d}x, \end{align}
   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 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).

A ‘$natural$’ space to look for the solution $u$ of problems (1.2)–(1.3), more accurately, to define the integrals on the left-hand side of (1.4), is as follows:

\[ \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 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 the above identification and employing the classical embedding theorem [46, Proposition 2.2], both the integrals on the 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 our problems (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 [40]) 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 [8, 21]):

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.$

To provide a meaningful interpretation for the integral in equation (2.5), we carefully select a suitable class of test functions denoted $\mathcal {T}(\Omega )$. The choice of this class depends on crucial factors, such as the exponent value $\gamma \geq 0$ and the regularity of both the datum $f$ and the solution $u$. This selection is pivotal for the rigorous analysis and understanding of the problem at hand. Motivating from the above discussion, the exact notion of a weak solution to our main problems (1.2)–(1.3) with different choices of a class of test functions and sufficient regularity of the solution $u$, is detailed in § 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 et al. [29], which majorly set-up this direction of research. Afterwards, 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, e.g. [30, 31, 48] 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 problems (1.2)–(1.3) with Laplace operator and $f=1$. By using the classical method of sub-supersolutions on the nonsingular approximating problem, they proved the existence and uniqueness results of the classical solution of our original problem. In addition, by exploiting the second-order ordinary differential equation techniques and localization near the boundary, the boundary behaviour of the solution is deduced. By Stuart [49], 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, 49] provide results for more general differential operators with smooth coefficients, not necessarily in divergence form, and for nonmonotone 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 [40] in which they simplified the proof of boundary behaviour of classical by constructing appropriate sub and super solutions. In addition to that, they also obtained the optimal power related to the existence of finite-energy solutions. In fact, a solution $u \in H_0^1(\Omega )$ exists if and only if $\gamma <3.$ In [50], Yijing and Zhang analysed the threshold value $3$ when the datum $f \in L^1(\Omega )$ and a positive function, and provided 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:

\[ \left\{\begin{array}{@{}ll} u \in W_0^{1,({Nr(1+\gamma)})/({N-r(1-\gamma)})}(\Omega) & \text{ if } \ 0<\gamma <1 \ \text{and} \ f \in L^r(\Omega) \ \text{with} \ r \\ & \quad\in [1,(2^\ast{/}(1-\gamma))'),\\ u \in H_0^{1}(\Omega) & \text{ if } \ 0<\gamma <1 \ \text{and} \ f \in L^r(\Omega) \ \text{with} \ r\\ & \quad =(2^\ast{/}(1-\gamma))',\\ u \in H_0^{1}(\Omega) & \text{ if } \ \gamma =1 \ \text{and} \ f \in L^1(\Omega),\\ u^{{(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 Moreno-Mérida in [5] studied some particular cases of strongly singular elliptic equations, i.e. $1< \gamma < ({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 Casado-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 [18, 30, 33, 48, 52] and for singular nonlinearities with source terms or absorption terms, we refer to [20, 28, 45, 49] 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 the presence of singular nonlinearities and weight function $f$. In particular, they proved the existence and regularity of solution in a very weak sense depending upon the regularity of the datum $f$ and the singular exponent $\gamma$:

\[ \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))',\\ u^{{(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 [51], Youssfi and Mahmoud extended the result of [5] to the fractional Laplacian and established the following:

\begin{align*} & u \in W_0^{1,{(Nr(1+\gamma))}/{(N-rs(1-\gamma))}}(\Omega) \ \text{if} \ 0<\gamma <1 \ \text{and} \ f\\ & \quad\in L^r(\Omega) \ \text{with} \ r \in [1,(2_s^\ast{/}(1-\gamma))'), \end{align*}

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 the 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 nonlinearities 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 the Green function and maximum principle, they proved the sharp boundary behaviour of the weak solution. For further issues on nonlocal and nonlinear singular problems, the interested reader can consult to the bibliographic references in [6, 8, 21, 23, 38].

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'$ the conjugate exponent of $q$, namely $q'={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)' \ \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 classes of weight functions $f$, and in both the presence or absence of singular nonlinearities, i.e. $\gamma >0$ or $\gamma =0$, respectively. In this section, we give a short description of our main results and for a detailed presentation, we refer the reader to § 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 the classical approach of regularizing the singular nonlinearities $u^{-\gamma }$ by $(u+{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 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\},\\ 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:= ({Nr(1+\gamma )})/({N-r(1-\gamma )}).$ Here, the notion of solution in two disjoint 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 complement the results for classical Laplacian in [16] and fractional Laplacian in [10, 51].

• • 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 an appropriate choice of test functions when $\alpha$ is large and the second by using the lower-boundary behaviour of the approximating solution when $\alpha$ is small.

  Type 1: When $1 \leq r < {N}/{2}$, $\gamma \geq 0$, we show that

  \[ u^{\alpha} \in H_0^1(\Omega) \ \text{for any}\ \alpha \in \left[\frac{\gamma+1}{2}, \frac{\mathfrak{S}_r+1}{2}\right] \ \text{where}\ \mathfrak{S}_r:=\frac{N(r-1) + \gamma r(N-2)}{N-2r}. \]
  We notice that as $r \to ({N}/{2})^-$, $\mathfrak {S}_r \to \infty$ and so it is natural to expect the exponential-type Sobolev regularity when $r={N}/{2}.$ In this regard, for $r= {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 the 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 is even new for the classical Laplacian and fractional Laplacian singular and nonsingular problems.

  Type 2: In this case, to handle the singular nonlinearities, we exploit the boundary behaviour 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} \displaystyle\left(\dfrac{1}{2},\dfrac{\gamma+1}{2}\right] & \text{if}\ \gamma+\dfrac{1}{r} < 1,\\ \left(\dfrac{r\gamma+ 1}{2r},\dfrac{\gamma+1}{2}\right] & \text{if}\ \gamma+\dfrac{1}{r} \geq 1. \end{array}\right. \]
  In addition to that, we also have shown that if $\alpha \leq {1}/{2}$, then $u^\alpha \not \in H_0^1(\Omega )$ if $\gamma +{1}/{r} < 1$, which highlight the optimality of this result.

• • Continuity with respect to datum: As an application of the type 1 power-type Sobolev regularity results, we show that for any $1 \leq r < {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|^{{(\mathfrak{S}+1)}/{2}}\|^2_{L^2(\Omega)} \leq C \|f-g\|_{L^r(\Omega)}^{{(r(N-2))}/{(N-2r)}}. \]
  The above inequality also implies the continuity of the solution with respect to the 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 Green's kernel associated with 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,\\ H^1_{loc}(\Omega) & \text{ if } \ \gamma+ \zeta > 1. \end{array}\right. \]

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

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

• • Optimal power-type Sobolev regularity and nonexistence results: As an application of the above optimal-boundary behaviour, we prove:

  \[ u^{{(\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, \\ \mathfrak{L}^* & \text{ if }\ \zeta+ \gamma > 1. \end{array}\right. \]
  and for $\zeta \geq 2$ and $\gamma >0$, nonexistence results is established. These results for mixed-type operator in this class complement the results for fractional Laplacian in [4, 7].

The article is organized as follows: In § 2, we provide a detailed statement of the main results of this work. In § 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 behaviour of the solution of the approximated problem. In § 5, we provide the proof of our main results. At the end, we provide a short appendix.

2. Statement of main results

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

2.1 Lebesgue weights

We start by 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 problems (1.2)-(1.3) in the following sense$:$

1. (i) $u \in W_0^{1,q}(\Omega )$ with $q:= ({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'}(\Omega ) \cap L^{r'}_c(\Omega )$

   (2.1)\begin{align} & \int_{\Omega} \nabla u \cdot \nabla \psi\,{\rm d}x + \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}}\,{\rm d}x\,{\rm d}y\nonumber\\ & \quad = \int_{\Omega} \frac{f(x)}{u^\gamma}\psi\,{\rm d}x. \end{align}

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 problems (1.2)-(1.3) in the following sense$:$

1. (i) $u \in H^1_{loc}(\Omega )$ and $u=0$ on $\Omega$ in the sense of definition 1.1,

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'}_c(\Omega )$

   (2.2)\begin{align} & \int_{\Omega} \nabla u \cdot \nabla \psi\,{\rm d}x + \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}}\,{\rm d}x\,{\rm d}y\nonumber\\ & \quad = \int_{\Omega} \frac{f(x)}{u^\gamma}\psi\,{\rm d}x. \end{align}

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 problems (1.2)-(1.3) obtained in theorems 2.1 and 2.2. Then,

1. (i) (Weak case) For $r \in [1, {N}/{2})$ and $N > 2$:

   \[ u^{{(\mathfrak{S}+1)}/{2}} \in H_0^1(\Omega) \ \text{for any} \ \mathfrak{S} \in [\gamma,\mathfrak{S}_r]\ \text{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:= ({Nr(1+\gamma )})/({N-2r}).$

2. (ii) (Limit case) For $r= {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
      (2.3)\begin{equation} \frac{2}{S(N) \|f\|_{{N}/{2}}} \geq \left\{\begin{array}{@{}ll} \beta \max\{1, \left(\beta\gamma^{{-}1}\right)^\gamma\} & \text{ if }\ 0 < \gamma \leq 1,\\ \beta & \text{ if }\ \gamma=0. \end{array}\right. \end{equation}
      Moreover, there exist constants $C_1, C_2$ depending upon $\beta, S(N),f, |\Omega |$ such that
      \[ \int_{\Omega} \exp\left(\frac{ \beta N u}{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 $({\alpha }/{(2^{3-2\tau })}) \max \{1, \alpha ^\tau \}< {1}/({\tau S(N) \|f\|_{{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
      \begin{align*} & \int_{\Omega} \exp\left(\frac{ 2N \alpha \tau u}{N-2}\right) \leq C_3 \quad \text{when}\ N >2 \quad \text{and} \\ & \|u_n\|_{L^\infty(\Omega)} \leq C_4 \quad \text{when} \ N =2. \end{align*}

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

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

4. (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:= ({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 problems (1.2)-(1.3) obtained in theorems 2.1 and 2.2. Then,

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

and

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

Theorem 2.5 Let $u$ and $v$ be two solutions of problems (1.2)–(1.3) obtained in theorems 2.1 and 2.2 with respect to datum $f$ and $g$ in $L^r(\Omega )$ for $1 \leq r< {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

\[ \|\nabla |u-v|^{{(\mathfrak{S}+1)}/{2}}\|^2_{L^2(\Omega)} \leq C\ \|f-g\|_{L^r(\Omega)}^{{(r(N-2))}/{(N-2r)}} . \]

In addition to the above, we have

\[ \int_{\Omega} |\nabla (u-v)_+^{{(\mathfrak{S}_r+1)}/{2}}|^2\,{\rm d}x \leq \int_{\Omega} (f(x)-g(x)) (u-v)_+^{\mathfrak{S}_r-\gamma}\,{\rm d}x. \]

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

Remark 2.7 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.8 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 ).$

Remark 2.9 We note that the regularity result in Theorem 2.13 of [35], where the exponent $p$ is set to 2, can be seen as a specific subcase of the results established in theorem 2.3(iv). This alignment becomes evident when considering the partial case $r = r^\ast$ and $0 \leq \gamma < 1$ of the exact Sobolev regularity results for the weak solution $u$ in problem (1.2) proven in theorem 2.3(iv).

2.2 Singular weights

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

Theorem 2.10 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 problems (1.2)-(1.3) in the following sense:

1. (i) $u \in H_0^{1}(\Omega ),$

2. (ii) for every $\psi \in H_0^{1}(\Omega )$

   (2.4)\begin{align} & \int_{\Omega} \nabla u \cdot \nabla \psi\,{\rm d}x + \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}}\,{\rm d}x\,{\rm d}y\nonumber\\ & \quad = \int_{\Omega} \frac{f(x)}{u^\gamma}\psi\,{\rm d}x. \end{align}

Theorem 2.11 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 problems (1.2)-(1.3) in the following sense:

1. (i) $u \in H^1_{loc}(\Omega )$ and $u=0$ on $\Omega$ in the sense of definition 1.1,

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$:

   (2.5)\begin{align} & \int_{\Omega} \nabla u \cdot \nabla \psi\,{\rm d}x + \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}}\,{\rm d}x\,{\rm d}y\nonumber\\ & \quad = \int_{\Omega} \frac{f(x)}{u^\gamma}\psi\,{\rm d}x. \end{align}

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

Theorem 2.12 Let $\gamma >0$, $\zeta \in [0,2)$ and $u$ be a weak solution of problems (1.2)-(1.3) obtained in theorems 2.10 and 2.11. Then,

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

and

\[ u^{{(\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, \\ \mathfrak{L}^* & \text{ if}\ \zeta+ \gamma > 1. \end{array}\right. \]

As an application of the above theorem, we prove the following nonexistence result:

Theorem 2.13 Let $\zeta \geq 2.$ Then, there does not exist a weak solution of problems (1.2)-(1.3) in the sense of theorems 2.10 and 2.11.

Remark 2.14 Repeating the same proof [8, Theorem 1.1], we can obtain the uniqueness of weak solution proved in theorems 2.10 and 2.11 when $\gamma >0$ and $\beta \in [0, {3}/{2})$.

Remark 2.15 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

(S)\begin{equation} \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. \end{equation}

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\,{\rm d}x + \frac{C(N,s)}{4} \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|u(x)-u(y)|^2}{|x-y|^{N+2s}}\,{\rm d}x\,{\rm d}y - \int_{\Omega} h u\,{\rm d}x. \]

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

\begin{align*} & \mathcal{I}_h(u) \geq \frac{1}{2} \|\nabla u\|_2^2 - |\Omega|^{{1}/{2}} \|h\|_{\infty} \|u\|_{2} \geq \|\nabla u\|_2 \left(\frac{1}{2} \|\nabla u\|_2 - S(N) |\Omega|^{{1}/{2}} \|h\|_{\infty}\right)\nonumber\\ & \quad \to \infty\ \text{as}\ \|\nabla u\|_2 \to \infty \end{align*}

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.:

(3.1)\begin{align} & \int_{\Omega} \nabla u \cdot \nabla \psi\,{\rm d}x + \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}}\,{\rm d}x\,{\rm d}y\nonumber\\ & \quad = \int_{\Omega} h \psi\,{\rm d}x \quad \forall \psi \in H_0^1(\Omega). \end{align}

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\,{\rm d}x \leq 0. \]

Hence, $u \geq 0$ a.e. in $\Omega.$ Now, we show that problem $(S)$ has a unique solution. Let $u_1, u_2 \in H_0^1(\Omega )$ be two solutions of problem $(S)$. Therefore, for all $\psi \in H_0^1(\Omega )$, we have

(3.2)\begin{equation} \int_{\Omega} \nabla u_1 \cdot \nabla \psi\,{\rm d}x + \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}}\,{\rm d}x\,{\rm d}y \!=\! \int_{\Omega} h \psi\,{\rm d}x, \end{equation}

(3.3)\begin{equation} \int_{\Omega} \nabla u_2 \cdot \nabla \psi\,{\rm d}x + \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}}\,{\rm d}x\,{\rm d}y \!=\! \int_{\Omega} h \psi\,{\rm d}x. \end{equation}

By subtracting the above two equations and inserting $\psi = u_1-u_2$, we obtain:

\[ \int_{\Omega} |\nabla (u_1- u_2)|^2\,{\rm d}x + \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}}\,{\rm d}x\,{\rm d}y =0, \]

which further gives $u_1=u_2$ a.e. in $\Omega.$ The boundedness and regularity of the solution follow by employing the classical method of Stampacchia, see 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 the second case, when $f \in \mathcal {A}_\zeta$, we consider:

\[ f_n(x):=\left\{\begin{array}{@{}ll} \left(f^{{-1}/{\zeta}}(x) + \left(\dfrac{1}{n}\right)^{{(\gamma+1)}/{(2-\zeta)}}\right)^{-\zeta} & \text{ if }\ x \in \Omega,\\ 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$:

(3.4)\begin{equation} \frac{\mathcal{G}_1}{\left(d(x) + \left(\frac{1}{n}\right)^{{(\gamma+1)}/{(2-\zeta)}}\right)^{\zeta}} \leq f_n(x) \leq \frac{\mathcal{G}_2}{\left(d(x) + \left(\frac{1}{n}\right)^{{(\gamma+1)}/{(2-\zeta)}}\right)^{\zeta}}. \end{equation}

By considering the above choice of $f_n$, we study the following approximated singular problem:

(P _n)\begin{equation} \left\{\!\begin{aligned} \Delta u + (-\Delta)^s u =\frac{f_n}{(u+{1}/{n})^\gamma},\quad u>0 & \text{ in } \Omega, \\\quad u = 0 & \text{ in } \ \mathbb{R}^N {\setminus} \Omega, \end{aligned}\right. \end{equation}

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 problem $(P_n)$ in the sense that:

(3.5)\begin{align} & \int_{\Omega} \nabla u_n \cdot \nabla \psi\,{\rm d}x + \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}}\,{\rm d}x\,{\rm d}y\nonumber\\ & \quad = \int_{\Omega} \frac{f_n(x)}{(u_n + {1}/{n})^\gamma} \psi\,{\rm d}x \quad \forall \psi \in H_0^1(\Omega). \end{align}

Moreover,

1. (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).$

2. (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}.$

3. (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^+{+} {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 problem $(S)$ (see statement of lemma 3.1) with $h$ defined above. Therefore, we define an operator $T: L^{2}(\Omega ) \to L^2(\Omega )$ such that $T(g)=w$ where $w$ satisfies:

(3.6)\begin{align} & \int_{\Omega} \nabla w \cdot \nabla \psi\,{\rm d}x + \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}}\,{\rm d}x\,{\rm d}y\nonumber\\ & \quad = \int_{\Omega} \frac{f_n(x)}{(g^+{+} {1}/{n})^\gamma} \psi\,{\rm d}x \quad \forall \psi \in H_0^1(\Omega). \end{align}

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

(3.7)\begin{equation} \int_{\Omega} |\nabla w|^2\,{\rm d}x \leq \int_{\Omega} \frac{f_n(x) w}{(g^+{+} {1}/{n})^\gamma}\,{\rm d}x \leq R \|\nabla w\|_{2} \Longrightarrow \|\nabla w\|_{2} \leq R \end{equation}

where $R:= S(N) \|f_n\|_{L^\infty (\Omega )} n^\gamma |\Omega |^{{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 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 sequence $\{w\}_{k \in \mathbb {N}}$ and choosing $w_k-w$ as a test function, we get:

(3.8)\begin{align} & \int_{\Omega} |\nabla w_k-w|^2\,{\rm d}x = \int_{\Omega} \left(\frac{f_n(x)}{(h_k^+{+} {1}/{n})^\gamma} - \frac{f_n(x)}{(h^+{+} {1}/{n})^\gamma}\right) (w_k-w)\,{\rm d}x\nonumber\\ & \quad \leq\! S(N) \left(\int_{\Omega} \left(\frac{f_n(x)}{(h_k^+{+} {1}/{n})^\gamma} - \frac{f_n(x)}{(h^+{+} {1}/{n})^\gamma}\right)^{{2N}/{(N\!+2)}}\,{\rm d}x\right)^{{(N\!+2)}/{2N}}\|\nabla w_k \!-\!\nabla w\|_{2}. \end{align}

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 )$:

(3.9)\begin{align} & \int_{\Omega} \nabla w_k \cdot \nabla \psi\,{\rm d}x + \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}}\,{\rm d}x\,{\rm d}y\nonumber\\ & \quad = \int_{\Omega} \frac{f_n(x)}{(h_k^+{+} {1}/{n})^\gamma} \psi\,{\rm d}x. \end{align}

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

\[ \frac{(w_k(x)- w_k(y))(\psi(x)- \psi(y))}{|x-y|^{{(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|^{{(N+2s)}/{2}}} \to \frac{(w(x)- w(y))(\psi(x)- \psi(y))}{|x-y|^{{(N+2s)}/{2}}}\ \text{a.e. in} \ \Omega. \]

Then, since

\[ \frac{(\psi(x)- \psi(y))(\psi(x)- \psi(y))}{|x-y|^{{(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

\begin{align*} & \int_{\Omega} \nabla w \cdot \nabla \psi\,{\rm d}x + \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}}\,{\rm d}x\,{\rm d}y\nonumber\\ & \quad = \int_{\Omega} \frac{f_n(x)}{(h^+{+} {1}/{n})^\gamma} \psi\,{\rm d}x. \end{align*}

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 point $u_n$ such that $T(u_n)= u_n$ for all $n \in \mathbb {N}$, i.e.:

\begin{align*} & \int_{\Omega} \nabla u_n \cdot \nabla \psi\,{\rm d}x + \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}}\,{\rm d}x\,{\rm d}y\nonumber\\ & \quad = \int_{\Omega} \frac{f_n(x)}{(u_n^+{+} {1}/{n})^\gamma} \psi\,{\rm d}x. \end{align*}

The boundedness and regularity of solution follow by employing the classical method of Stampacchia, see 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 ${f_n}/{((u^+_n + {1}/{n} )^\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 problem $(P_n)$ and $(P_{n+1})$ respectively, i.e. for any $\psi \in H_0^1(\Omega )$, we have

\begin{align*} & \int_{\Omega} \nabla u_n \cdot \nabla \psi\,{\rm d}x + \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}}\,{\rm d}x\,{\rm d}y\nonumber\\ & \quad = \int_{\Omega} \frac{f_n(x)}{(u_n + {1}/{n})^\gamma} \psi\,{\rm d}x \end{align*}

and

\begin{align*} & \int_{\Omega} \nabla u_{n+1} \cdot \nabla \psi\,{\rm d}x + \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}}\,{\rm d}x\,{\rm d}y\nonumber\\ & \quad = \int_{\Omega} \frac{f_{n+1}(x)}{(u_{n+1} + ({1}/{(n+1)}))^\gamma} \psi\,{\rm d}x. \end{align*}

Subtracting the above equalities by taking the test function $\psi = (u_n - u_{n+1})^+$ and using the following inequality [41, 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{align*} & \int_{u_n \geq u_{n+1}} |\nabla (u_n - u_{n+1})|^2\,{\rm d}x \\ & \quad \leq \int_{\Omega} \left( \frac{f_{n}(x)}{(u_{n} + {1}/{n})^\gamma} -\frac{f_{n+1}(x)}{\left(u_{n+1} + ({1}/{(n+1)})\right)^\gamma}\right) (u_n - u_{n+1})^+ \,{\rm d}x\\ & \quad \leq \int_{\Omega} f_{n+1}(x) \frac{(u_{n} + {1}/{n})^\gamma- (u_{n+1} + ({1}/{(n+1)}))^\gamma}{\left((u_{n} + {1}/{n})(u_{n+1} + ({1}/{(n+1)}))\right)^\gamma} (u_n - u_{n+1})^+\,{\rm d}x \leq 0 \end{align*}

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 behaviour of the approximated sequence $u_n$ with the help of the integral representation of the solution via Green's operator.

Theorem 3.3 Let $u_n$ be the weak solution of 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 a sequence of weak solutions of problem $(P_n)$ and $G(x,y)$ the Green function associated with the mixed operator with homogeneous Dirichlet-boundary conditions in $\Omega.$ Then, using the integral representation of the solution, we have

(3.10)\begin{equation} u_n(x) := \mathbb{G}^\Omega\left[\frac{f_n(y)}{(u_n + {1}/{n})^\gamma}\right](x)= \int_{\Omega} \frac{G(x,y) f_n(y)}{(u_n + {1}/{n})^\gamma} ~{\rm d}y \quad \text{for} \ n \in \mathbb{N}. \end{equation}

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

(3.11)\begin{align} \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)\nonumber\\ & \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)\nonumber\\ & \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{align}

4. Uniform a priori 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 problem $(P_n)$:

1. (i) For $r \in [1, {N}/{2})$ and $N > 2$:

   \[ u_n^{{(\mathfrak{S}_r+1)}/{2}}\ \text{is uniformly bounded in}\ H_0^1(\Omega)\ \text{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:= ({Nr(1\!+\!\gamma )})/({N\!-\!2r}).$

2. (ii) For $r= {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
      \begin{align*} & \int_{\Omega} \exp\left(\frac{ \beta N u_n}{N-2}\right) \leq C_1 \quad \text{when} \ N >2 \quad \text{and} \quad \|u_n\|_{L^\infty(\Omega)}\\ & \quad \leq C_2 \quad \text{when} \ N =2. \end{align*}

   2. 2. For $\tau \geq \gamma > 1$ and $\alpha >0$ such that $({\alpha }/{(2^{3-2\tau })}) \max \{1, \alpha ^\tau \}< {1}/(\tau S(N) \|f\|_{{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
      \begin{align*} & \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)}\\ & \quad \leq C_4 \quad \text{when} \ N =2. \end{align*}

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

3. (iii) For $r > {N}/{2}$ and $N \geq 2$, $u_n$ is uniformly bounded in $L^\infty (\Omega ).$

4. (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:= ({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 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, {1}/{n})$ and $\mathfrak {S} \in [\gamma, \infty )$, it yields

(4.1)\begin{align} & \int_{\Omega} \nabla u_n \cdot \nabla (u_n + \epsilon)^{\mathfrak{S}}\,{\rm d}x + \frac{C(N,s)}{2} \int_{\mathbb{R}^N} \nonumber\\ & \quad\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}}\,{\rm d}x\,{\rm d}y \nonumber\\ & \quad \leq \int_{\Omega} \frac{f_n(x)}{(u_n + {1}/{n})^\gamma} (u_n + \epsilon)^{\mathfrak{S}}\,{\rm d}x \leq \int_{\Omega} f_n(x) (u_n + \epsilon)^{\mathfrak{S} - \gamma}\,{\rm d}x. \end{align}

Passing $\epsilon \to 0$ in the above estimate via Fatou's theorem and using lemma A.4(i), we obtain:

(4.2)\begin{align} & \frac{4 \mathfrak{S}}{(\mathfrak{S} +1)^2}\int_{\Omega} |\nabla u_n^{{(\mathfrak{S} +1)}/{2}}|^2\,{\rm d}x + \frac{2 C(N,s) \mathfrak{S}}{(\mathfrak{S} +1)^2} \int_{\mathbb{R}^N}\nonumber\\ & \quad \int_{\mathbb{R}^N} \frac{((u_n^{{(\mathfrak{S} +1)}/{2}}(x)- u_n^{{(\mathfrak{S} +1)}/{2}}(y))^2}{|x-y|^{N+2s}}\,{\rm d}x\,{\rm d}y \nonumber\\ & \quad \leq \int_{\Omega} f_n(x) u_n^{\mathfrak{S} - \gamma}\,{\rm d}x. \end{align}

In order to estimate the right-hand side term of (4.2), we choose $\mathfrak {S} = \mathfrak {S}_r$ such that

(4.3)\begin{equation} \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. \end{equation}

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

(4.4)\begin{align} \int_{\Omega} f_n(x) u_n^{\mathfrak{S}_r - \gamma}\,{\rm d}x & \leq \|f_n\|_{L^r(\Omega)} \left(\int_{\Omega} u_n^{(\mathfrak{S}_r - \gamma) r'}\,{\rm d}x\right)^{{1}/{r'}}\nonumber\\ & \quad = \|f_n\|_{L^r(\Omega)} \left(\int_{\Omega} \left(u_n^{{(\mathfrak{S}_r +1)}/{2}}\right)^{{2N}/{(N-2)}}\,{\rm d}x\right)^{{1}/{r'}}\nonumber\\ & \leq \|f_n\|_{L^r(\Omega)} \left( S(N) \int_{\Omega} |\nabla u_n^{{(\mathfrak{S}_r +1)}/{2}}|^2\,{\rm d}x \right)^{{N}/{(r'(N-2))}} \end{align}

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:

(4.5)\begin{align} & \int_{\Omega} |\nabla u_n^{{(\mathfrak{S}_r +1)}/{2}}|^2\,{\rm d}x \leq \left( S(N)\right)^{{(N(r-1))}/{(N-2r)}} \nonumber\\ & \quad \left(\|f_n\|_{L^r(\Omega)}\frac{(\mathfrak{S}_r +1)^2}{4 \mathfrak{S}_r}\right)^{{((N-2)r)}/{(N-2r)}} \leq C(\|f\|_{L^r(\Omega)}, \mathfrak{S}_r, N). \end{align}

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:

(4.6)\begin{equation} \frac{N(\mathfrak{S}_r +1)}{N-2} = \frac{Nr(1+\gamma)}{N-2r} = \sigma_r \end{equation}

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= {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 computation leads to:

(4.7)\begin{equation} \int_{\Omega} \nabla u_n \cdot \nabla \mathfrak{H}(u_n)\,{\rm d}x= \int_{\Omega} \mathfrak{H}'(u_n) |\nabla u_n|^2\,{\rm d}x= \frac{4}{\beta} \int_{\Omega} \left| \nabla \mathfrak{H} \left(\frac{u_n}{2}\right)\right|^2\,{\rm d}x \end{equation}

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

(4.8)\begin{align} \frac{4}{\beta} & \int_{\Omega} \left| \nabla \mathfrak{H} \left(\frac{u_n}{2}\right)\right|^2\,{\rm d}x \nonumber\\ & \leq \int_{\Omega} \nabla u_n \cdot \nabla \mathfrak{H}(u_n)\,{\rm d}x + \frac{C(N,s)}{2} \int_{\mathbb{R}^N}\nonumber\\ & \quad \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}}\,{\rm d}x\,{\rm d}y \nonumber\\ & \leq \int_{\Omega} \frac{f_n(x)}{(u_n + {1}/{n})^\gamma} \mathfrak{H}(u_n)\,{\rm d}x \leq \int_{\Omega} \frac{f(x) \mathfrak{H}(u_n)}{u_n^\gamma}\,{\rm d}x\nonumber\\ & \quad = \int_{\Omega} \frac{f(x) \left(\exp(\beta u_n)-1\right)}{u_n^\gamma}\,{\rm d}x \end{align}

Case 1: $0 < \gamma \leq 1.$

To estimate the last term in (4.8), we use lemma A.4(iii) and $\exp (t)-1 \leq t \exp (t)$ for every $t \geq 0$,

\begin{align*} \int_{\Omega} \left| \nabla \mathfrak{H} \left(\frac{u_n}{2}\right)\right|^2\,{\rm d}x & \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\,{\rm d}x \\ & \quad + \frac{\beta}{4} \int_{\Omega \cap \{u_n \geq 1\}} f(x) \exp\left(\beta u_n\right)\,{\rm d}x\\ & \leq C^\sharp \int_{\Omega} f(x) \left(\exp\left(\frac{\beta u_n}{2}\right)\right)^2 ~{\rm d}x\nonumber\\ & \quad = C^\sharp \int_{\Omega} f(x) \left(\mathfrak{H} \left(\frac{u_n}{2}\right) +1\right)^2\,{\rm d}x \end{align*}

where $C^\sharp := ({\beta }/{4}) \max \{1, (\beta \gamma ^{-1})^\gamma \}.$

Case 2: $\gamma =0.$

In this case, the estimate in (4.8) takes the following form:

\begin{align*} & \int_{\Omega} \left| \nabla \mathfrak{H} \left(\frac{u_n}{2}\right)\right|^2\,{\rm d}x \leq C^\sharp \int_{\Omega} f(x) \left(\exp\left(\frac{\beta u_n}{2}\right)\right)^2\,{\rm d}x\nonumber\\ & \quad = C^\sharp \int_{\Omega} f(x) \left(\mathfrak{H} \left(\frac{u_n}{2}\right) +1\right)^2\,{\rm d}x \end{align*}

with $C^\sharp ={\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 ${N}/{2}$ and ${N}/{(N-2)}$, we further estimate:

\begin{align*} \int_{\Omega} \left| \nabla \mathfrak{H} \left(\frac{u_n}{2}\right)\right|^2\,{\rm d}x & \leq 2 C^\sharp \left(\int_{\Omega} f(x)\,{\rm d}x + \int_{\Omega} f(x) \left(\mathfrak{H} \left(\frac{u_n}{2}\right)\right)^2\,{\rm d}x \right)\\ & \leq 2 C^\sharp \left( \|f\|_1 + \|f\|_{{N}/{2}} \left\| \mathfrak{H} \left(\frac{u_n}{2}\right)\right\|_{{2N}/{(N-2)}}^2 \right). \end{align*}

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

\[ \frac{2}{S(N) \|f\|_{{N}/{2}}} \geq \left\{\begin{array}{ll} \beta \max\{1, \left(\beta\gamma^{{-}1}\right)^\gamma\} & \text{ if } 0 < \gamma \leq 1,\\ \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{2 C^\sharp \|f\|_1}{\left(1 - 2 C^\sharp S(N) \|f\|_{{N}/{2}}\right)} \]

and

\[ \int_{\Omega} \exp \left(\frac{ \beta N u_n}{N-2}\right)\,{\rm d}x \leq C \]

where $C$ depends upon $\beta, S(N), \|f\|_{{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 A.4(iii) yields $\mathfrak {D}(u_n) \in H_0^1(\Omega )$ and

(4.9)\begin{align} \int_{\Omega} \nabla u_n \cdot \nabla \mathfrak{D}(u_n)\,{\rm d}x & = \int_{\Omega} \mathfrak{D}'(u_n) |\nabla u_n|^2\,{\rm d}x\nonumber\\ & \quad = \alpha \tau \int_{\Omega} \left| \exp\left(\frac{\alpha u_n}{2}\right) \left(\exp(\alpha u_n) -1 \right)^{{(\tau-1)}/{2}} \nabla u_n \right|^2\,{\rm d}x \nonumber\\ & \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\,{\rm d}x\nonumber\\ & = \frac{4}{\alpha \tau} \int_{\Omega} \left|\nabla \left(\exp\left(\frac{\alpha u_n}{2}\right) -1 \right)^{\tau} \right|^2\,{\rm d}x \!=\! \frac{4}{\alpha \tau} \int_{\Omega} \left|\nabla \mathfrak{D}\left(\frac{u_n}{2}\right) \right|^2\,{\rm d}x \end{align}

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

\begin{align*} \frac{4}{\alpha \tau} & \int_{\Omega} \left|\nabla \mathfrak{D}\left(\frac{u_n}{2}\right) \right|^2\,{\rm d}x \\ & \leq \int_{\Omega} \nabla u_n \cdot \nabla \mathfrak{D}(u_n)\,{\rm d}x + \frac{C(N,s)}{2} \int_{\mathbb{R}^N}\nonumber\\ & \quad \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}}\,{\rm d}x\,{\rm d}y \\ & \leq \int_{\Omega} \frac{f_n(x)}{(u_n + {1}/{n})^\gamma} \mathfrak{D}(u_n)\,{\rm d}x \leq \int_{\Omega} \frac{f(x) \mathfrak{D}(u_n)}{u_n^\gamma}\,{\rm d}x = \int_{\Omega} \frac{f(x) \left(\exp(\alpha u_n)-1\right)^\tau}{u_n^\gamma}\,{\rm d}x \\ & \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\,{\rm d}x + \int_{\Omega \cap \{u_n \geq 1\}} f(x) \exp\left(\alpha \tau u_n\right)\,{\rm d}x\\ & \leq \max\{1, \alpha^\tau\} \int_{\Omega } f(x) \exp\left(\alpha \tau u_n\right)\,{\rm d}x \end{align*}

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) = (\mathfrak {D}^{{1}/{\tau }}(t) +1)^{2\tau } \leq 2^{2\tau -1} (\mathfrak {D}(t))^2+1)$ and Hölder inequality with exponents ${N}/{2}$ and ${N}/{(N-2)}$, we obtain:

\begin{align*} \int_{\Omega} \left| \nabla \mathfrak{D} \left(\frac{u_n}{2}\right)\right|^2\,{\rm d}x & \leq C^* \int_{\Omega} f(x) \left(\exp\left(\frac{\alpha u_n}{2}\right)\right)^{2\tau} \,{\rm d}x\nonumber\\ & \quad \leq 2^{2\tau-1} C^* \int_{\Omega} f(x) \left(\mathfrak{D}^{{1}/{\tau}} \left(\frac{u_n}{2}\right) +1\right)^{2\tau}\,{\rm d}x\\ & \leq 2^{2\tau-1} C^* \left(\int_{\Omega} f(x)\,{\rm d}x + \int_{\Omega} f(x) \left(\mathfrak{D} \left(\frac{u_n}{2}\right)\right)^2\,{\rm d}x \right)\\ & \leq 2^{2\tau-1} C^* \left( \|f\|_1 + \|f\|_{{N}/{2}} \left\| \mathfrak{D} \left(\frac{u_n}{2}\right)\right\|_{{2N}/{(N-2)}}^2 \right) \\ & \qquad C^*:= \frac{\alpha \tau}{4} \max\{1, \alpha^\tau\}. \end{align*}

Now, by choosing $\alpha$ small enough such that $\alpha \max \{1, \alpha ^\tau \}< {(2^{3-2\tau })}/{(\tau S(N) \|f\|_{{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\|_{{N}/{2}}\right)}\quad \text{and} \quad \int_{\Omega} \exp \left(\frac{2N \alpha \tau u_n}{N-2}\right)\,{\rm d}x \leq C \]

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

Step 3: To prove the boundedness result when $r > {N}/{2}$, we use the classical arguments from the seminal paper of Stampacchia [47]. 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:

(4.10)\begin{align} & \left(\int_{A_k} |G_k(u_n)|^{{2N}/{(N-2)}}\,{\rm d}x \right)^{{(N-2)}/{N}} \nonumber\\ & \quad \leq \int_{A_k} |\nabla G_k(u_n)|^2 \,{\rm d}x = \int_{\Omega} \nabla u_n \cdot \nabla G_k(u_n)\,{\rm d}x\nonumber\\ & \leq \int_{\Omega} \frac{f_n}{\left(u_n + {1}/{n}\right)^\gamma} G_k(u_n)\,{\rm d}x \leq \int_{A_k} f_n(x) G_k(u_n)\,{\rm d}x\nonumber\\ & \leq C \|f\|_{L^r(\Omega)} \left( \int_{A_k} (G_k(u_n))^{{2N}/{(N-2)}}\right)^{{(N-2)}/{2N}} |A_k|^{1-({(N-2)}/{2N})- {1}/{r}} \end{align}

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{align*} |h-k| |A_h|^{{(N-2)}/{2N}} & \leq \left(\int_{A_h} |G_k(u_n)|^{{2N}/{(N-2)}}\,{\rm d}x \right)^{{(N-2)}/{2N}}\nonumber\\ & \leq \left(\int_{A_k} |G_k(u_n)|^{{2N}/{(N-2)}}\,{\rm d}x \right)^{{(N-2)}/{2N}} \\ & \leq C \|f\|_{L^r(\Omega)} |A_k|^{1-({(N-2)}/{2N})- {1}/{r}} \end{align*}

which further implies:

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

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

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

Hence, we apply lemma A.5 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

(4.11)\begin{equation} {\frac{(1-\mathfrak{S}_r)q}{2-q}} = \frac{Nr (1+\gamma)}{N-2r}=\sigma_r. \end{equation}

By applying Hölder inequality with exponent ${2}/{q}$ and ${2}/{(2-q)}$ and using claim in step 1, we get:

\begin{align*} \int_{\Omega} |\nabla u_n|^q\,{\rm d}x & = \int_{\Omega} \frac{|\nabla u_n|^q}{u_n^{{((1-\mathfrak{S}_r)q)}/{2}}} u_n^{{((1-\mathfrak{S}_r)q)}/{2}}\,{\rm d}x \\ & \leq \left(\int_{\Omega} u_n^{\mathfrak{S}_r-1} |\nabla u_n|^2\,{\rm d}x \right)^{{q}/{2}} \left(\int_{\Omega} u_n^{{((1-\mathfrak{S}_r)q)}/{(2-q)}}\,{\rm d}x \right)^{{(2-q)}/{2}}\\ & \leq \left(\frac{4}{(\mathfrak{S}_r+1)^2}\int_{\Omega} |\nabla u_n^{{(\mathfrak{S}_r+1)}/{2}}|^2\,{\rm d}x \right)^{{q}/{2}} \left(\int_{\Omega} u_n^{\sigma_r}\right)^{{(2-q)}/{2}}\nonumber\\ & \quad \leq C \ \text{(independent of $n$)}. \end{align*}

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.2

• • 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.3 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 problem $(P_n)$. Then,

\[ u_n^{{(\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+\dfrac{1}{r} < 1,\\ \left(\gamma-1\!+\!\dfrac{1}{r}, \gamma \right] & \text{ if } \ \gamma\!+\!\dfrac{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, {1}/{n})$ and $\mathfrak {S} \in (0,\gamma ]$ and using lower-boundary estimates in theorem 3.3, we obtain:

(4.12)\begin{align} & \int_{\Omega} \nabla u_n \cdot \nabla (u_n + \epsilon)^{\mathfrak{S}}\,{\rm d}x + \frac{C(N,s)}{2} \int_{\mathbb{R}^N} \nonumber\\ & \quad\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}}\,{\rm d}x\,{\rm d}y \nonumber\\ & \quad \leq \int_{\Omega} \frac{f_n(x)}{(u_n + \epsilon)^{\gamma-\mathfrak{S}}}\,{\rm d}x \leq C \int_{\Omega} \frac{f(x)}{\delta^{\gamma-\mathfrak{S}}(x)}\,{\rm d}x \nonumber\\ & \quad \leq C \|f\|_{L^r(\Omega)} \left(\int_{\Omega} \delta^{{({-}r(\gamma-\mathfrak{S}))}/{(r-1)}}(x)\,{\rm d}x\right)^{{(r-1)}/{r}} \nonumber\\ & \quad \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{align}

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+\dfrac{1}{r} \!<\! 1,\\ \left(\gamma-1\!+\!\dfrac{1}{r}, \gamma \right] & \text{ if} \ \gamma\!+\!\dfrac{1}{r} \!\geq\! 1. \end{array}\right. \]

Now, by using the same arguments as in the 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 a series of lemmas involving the upper and lower estimate of the action of the Green operator on the logarithmic perturbation of the distance function. A similar type of Green estimate for a large class of nonlocal operators 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 ({\mathcal {D}_\Omega }/{t}).$ For $x \in \Omega _\eta$, we denote $\phi _x: B(x,1) \to B(0,1)$ be a diffeomorphism such that:

\begin{align*} & \phi_x(\Omega \cap B(x,1)) = B(0,1) \cap \{y \in \mathbb{R}^N : y \cdot e_N >0\}, \\ & \phi_x(y) \cdot e_N = \delta(y) \ \text{for} \ y \in B(x,1)\ \text{and} \ \phi_x(x)= \delta(x) e_N. \end{align*}

Lemma 4.4 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:

(4.13)\begin{align} \mathbb{G}^\Omega\left[ \frac{\ell^{-\Xi}(\delta({\cdot}))}{\delta({\cdot})} \right](x) & = \sum_{i=1}^5 \int_{\Omega_i} \frac{G^{\Omega}(x,y)}{\delta(y)} \ell^{-\Xi}(\delta(y))\,{\rm d}y + \int_{\Omega {\setminus} \Omega_{\eta}} \frac{G^{\Omega}(x,y)}{\delta(y)} \ell^{-\Xi}(\delta(y))\,{\rm d}y \nonumber\\ & := \sum_{i=1}^5 I_1^{(i)}(x) + I_2(x). \end{align}

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 to find the lower estimate of the term $I_1^{(4)}$ when $x \in \Omega _{{\eta }/{2}}$, and the lower estimate of the term $I_2$ when $x \in \Omega {\setminus} \Omega _{{\eta }/{2}}.$ For $y \in \Omega _4$ and $x \in \Omega _{{\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,

(4.14)\begin{equation} I_1(x) \geq \int_{\Omega_4} \frac{G^{\Omega}(x,y)}{\delta(y)} \ell^{-\Xi}(\delta(y))\,{\rm d}y \geq \ell^{-\Xi}(\delta(x)) \int_{\Omega_4} \frac{G^\Omega(x,y)}{\delta(y)}\,{\rm d}y := \ell^{-\Xi}(\delta(x)) \mathcal{J} \end{equation}

Now, by using estimates on [2, proof of Lemma 3.3, p. 40] and performing change of variables via diffeomorphism $\phi _x$, we obtain:

(4.15)\begin{align} & \mathcal{J} \gtrsim \delta(x) \int_{\{{3\delta(x)}/{2} < w_N < \eta\} \cap B(0,1)} \frac{1}{(|\delta(x)-w_N| + |w'|)^{N}}\,{\rm d}w_N\,{\rm d}w' \nonumber\\ & \quad = \delta(x) \int_{3/2}^{\eta/\delta(x)} \int_0^{1/\delta(x)} \frac{t^{N-2} }{(|1-h|+t)^{N}}\,{\rm d}t\,{\rm d}h \nonumber\\ & \quad = \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}}\,{\rm d}r\,{\rm d}h \nonumber\\ & \quad \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} }\,{\rm d}r\,{\rm d}h \gtrsim \delta(x) \int_{3/2}^{\eta/\delta(x)} \frac{1}{(h-1)}\,{\rm d}h \nonumber\\ & \quad \gtrsim \delta(x) \int_{3/2}^{\eta/\delta(x)} \frac{1}{h} \,{\rm d}h = \delta(x) \left(\ln\left(\frac{\eta}{\delta(x)}\right) - \ln\left(\frac{3}{2} \right) \right). \end{align}

By combining (4.13)–(4.15), we obtain:

(4.16)\begin{equation} \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_{{\eta}/{2}} \end{equation}

where the constant is independent of the parameter $\Xi.$ Let $x \in \Omega {\setminus} \Omega _{{\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:

(4.17)\begin{align} I_2(x) & = \int_{\Omega {\setminus} \Omega_{\eta}} \frac{G^{\Omega}(x,y)}{\delta(y)} \ell^{-\Xi}(\delta(y))\,{\rm d}y \geq \ell^{-\Xi}(\eta) \mathbb{G}^\Omega\left[ \frac{\chi_{\Omega {\setminus} \Omega_{\eta}}}{\delta}\right](x)\nonumber\\ & \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{align}

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

Lemma 4.5 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 _{{\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 < {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{align*} \int_{\Omega_1} \frac{G^{\Omega}(x,y)}{\delta(y)} \ell^{-\Xi}(\delta(y))\,{\rm d}y & \leq c^{{-}1} \ell^{-\Xi}(\delta(x)) \frac{1}{\delta(x)} \int_{B(x, \delta(x)/2)} \frac{1}{|x-y|^{N-2}}\,{\rm d}y \\ & \leq c^{{-}1} \ell^{-\Xi}(\delta(x)) \delta(x) \leq C \ell^{1-\Xi}(\delta(x)) \delta(x), \end{align*}

where $C$ is independent of parameter $\Xi$.

Estimate over $\Omega _2$:

\[ \int_{\Omega_2} \frac{G^{\Omega}(x,y)}{\delta(y)} \ell^{-\Xi}(\delta(y))\,{\rm d}y \leq \ell^{-\Xi}\left(\eta\right) \delta(x) \int_{\Omega_2} \frac{1}{|x-y|^{n}}\,{\rm d}y \leq C(\eta) \ \ell^{1-\Xi}(\delta(x)) \delta(x), \]

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{align*} \int_{\Omega_3} & \frac{G^{\Omega}(x,y)}{\delta(y)} \ell^{-\Xi}(\delta(y))\,{\rm d}y \leq \ell^{-\Xi}(\delta(x)) \delta(x) \int_{\Omega_3} \frac{1}{|x-y|^{N}}\,{\rm d}y \\ & \lesssim \ell^{-\Xi}(\delta(x)) \delta(x) \int_{|w'| <1} \int_0^{\delta(x)/2} \frac{1}{(|\delta(x)-w_N| + |w'|)^{N}}\,{\rm d}w_N\,{\rm d}w'\\ & \lesssim \ell^{-\Xi}(\delta(x)) \delta(x) \int_0^{1/\delta(x)} \int_0^{1/2} \frac{h^{N-2} }{((1-t)+h)^{N}}\,{\rm d}t\,{\rm d}h\\ & \lesssim \ell^{-\Xi}(\delta(x)) \delta(x) \int_0^{1/\delta(x)} \frac{h^{N-2}}{(1+h)^{N}}\,{\rm d}h \leq C(\eta) \ell^{1-\Xi}(\delta(x)) \delta(x). \end{align*}

Estimate over $\Omega _4$: We have

\begin{align*} \int_{\Omega_4} & \frac{G^{\Omega}(x,y)}{\delta(y)} \ell^{-\Xi}(\delta(y))\,{\rm d}y \leq \delta(x) \int_{\Omega_4} \frac{1}{|x-y|^{N} \ell^{\Xi}(\delta(y))}\,{\rm d}y \\ & \lesssim \delta(x) \int_{\{{3\delta(x)}/{2} < w_N < \eta\} \cap B(0,1)} \frac{1}{(|\delta(x)-w_N| + |w'|)^{N} \ell^{\Xi}(w_N)}\,{\rm d}w_N\,{\rm d}w'\\ & = \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)) } \,{\rm d}t\,{\rm d}h \\ & = \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} }\,{\rm d}r\,{\rm d}h \\ & \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} }\,{\rm d}r\,{\rm d}h \\ & \lesssim \delta(x) \int_{3/2}^{\eta/\delta(x)} \frac{1}{(h-1) \ell^{\Xi}(h\delta(x)) }\,{\rm d}h \lesssim \delta(x) \int_{3/2}^{\eta/\delta(x)} \frac{ \ell^{-\Xi}(h\delta(x)) }{h}\,{\rm d}h \\ & \lesssim \delta(x) \int^{\ln(2\mathcal{D}_{\Omega}/3\delta(x))}_{\ln(\mathcal{D}_{\Omega}/\eta)} \frac{1}{t^\Xi}\,{\rm d}t \leq \frac{C(\eta)}{(1-\Xi)} \ell^{1-\Xi}(\delta(x)) \delta(x) \end{align*}

Estimate over $\Omega _5$: Again, by choosing $\eta$ small enough such that $\eta < {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{align*} \int_{\Omega_5} & \frac{G^{\Omega}(x,y)}{\delta(y)} \ell^{-\Xi}(\delta(y)) ~{\rm d}y \lesssim \delta^{\gamma}(x) \int_{\Omega_5} \frac{1}{|x-y|^{N} \ell^{-\Xi}(\delta(y)) } \,{\rm d}y \\ & \lesssim \delta(x) \ell^{-\Xi}(\delta(x)) \int_{\Omega_5} \frac{1}{|x-y|^{N}}\,{\rm d}y \lesssim \delta(x) \ell^{-\Xi}(\delta(x)) \int_{\delta(x)/2}^{3 \delta(x)/2}\nonumber\\ & \quad \int_{\delta(x)/2}^1 \frac{t^{N-2}}{((\delta(x)-h)+t)^{N}} \,{\rm d}t\,{\rm d}h\\ & \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}}\,{\rm d}r\,{\rm d}t\\ & \lesssim \delta(x) \ell^{-\Xi}(\delta(x)) \int_{1/2}^{1/\delta(x)} \rho^{{-}1} \int_{0}^{1/\rho} \frac{1}{(r+1)^{N}}\,{\rm d}r\,{\rm d}\rho \leq C(\eta) \ell^{1-\Xi}(\delta(x)) \delta(x). \end{align*}

Case 2: $x \in \Omega {\setminus} \Omega _{{\eta }/{2}}$.

Using the estimates from the proof of Lemma 3.2 [2, p. 37], for $\eta$ small enough, we obtain:

\begin{align*} I_1(x) & \leq \ell^{-\Xi}(\eta) \delta(x) \left( \int_{\Omega_{\eta/4}}\frac{G^{\Omega}(x,y)}{\delta(y)}\,{\rm d}y + \int_{\Omega_{\eta} {\setminus} \Omega_{\eta/4}} \frac{G^{\Omega}(x,y)}{\delta(y)}\,{\rm d}y \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{align*}

For any $x \in \Omega$, we have the following estimate for $I_2(x)$:

(4.18)\begin{align} 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({\eta}/{2}\right)}\nonumber\\ & \leq C(\mathcal{D}_{\Omega},\eta) \ell^{1-\Xi}(\delta(x)) \delta(x). \end{align}

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

(4.19)\begin{equation} \beta:= \frac{2\gamma+ \zeta}{\gamma+1} \quad \kappa:=2-\beta = \frac{2-\zeta}{\gamma+1}. \end{equation}

Lemma 4.6 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:

(4.20)\begin{equation} \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_{{\eta}/{2}} \end{equation}

and

(4.21)\begin{equation} \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_{{\eta}/{2}}^c. \end{equation}

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

(4.22)\begin{equation} \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}}. \end{equation}

Using (4.22), we obtain the following:

(4.23)\begin{align} \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}\,{\rm d}y \nonumber\\ & \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) \,{\rm d}y,\nonumber\\ & \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) \,{\rm d}y,\nonumber\\ & \geq \frac{C'}{(\delta(x)+ \epsilon_1)^\beta} \int_{\mathcal{B}_\delta^x} \frac{1}{|x-y|^{N-2}}\,{\rm d}y = \frac{C'' \delta(x)^{2}}{(\delta(x)+ \epsilon_1)^\beta}\nonumber\\ & \geq C'' \left(\frac{1}{2} (\delta(x)+ \epsilon_1)^\kappa - \frac{\epsilon_1^2}{(\delta(x)+ \epsilon_1)^\beta}\right)\nonumber\\ & \geq C_1 \left( \frac{1}{2} (\delta(x)+\epsilon^{1/\kappa})^{\kappa} -\epsilon\right) . \end{align}

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

\[ \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, \]

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

4.3 Singular weights

4.3.1 Boundary behaviour

Let $\varphi _1$ and $\lambda _1$ are first eigenfunction and eigenvalue for the mixed local–nonlocal operator $(-\Delta ) + (-\Delta )^s$ such that (see [24, Proposition 3.7] and [17, Propositions 5.3, 5.4]):

(4.24)\begin{equation} \varphi_1 \asymp \delta^\gamma\ \text{and}\ \varphi_1= \lambda_1 \mathbb{G}^\Omega[\varphi_1] \quad \text{in} \ \Omega . \end{equation}

Lemma 4.7 Let $u_n$ be the weak solution of 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

   (4.25)\begin{equation} C_1 \delta(x) \leq u_n(x) \leq C_2 \delta(x) \left\{ \begin{array}{@{}ll} 1 & \text{ if } \ \zeta + \gamma < 1, \\ \ln\left(\dfrac{\mathcal{D}_\Omega}{\delta(x)}\right) & \text{ if} \ \zeta + \gamma = 1, \end{array} \right.\ \text{for} \ x \in \Omega. \end{equation}

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

   (4.26)\begin{align} & C^1 (\delta(x) + n^{({-(1+\gamma)})/({(2-\zeta)})})^{{(2-\zeta)}/{(\gamma+1)}} - C^2 n^{{-}1} \leq u_n(x)\nonumber\\ & \quad \leq C^3 \delta^{{(2-\zeta)}/{(\gamma+1)}}(x) \ \text{for} \ x \in \Omega. \end{align}

Proof. To prove the boundary behaviour, we divide our study into three cases:

Case 1: $\zeta + \gamma \leq 1$.

Using the integral representation in (3.10), lower-boundary behaviour in (3.11), $f_n \leq f$, and lemma A.3, we get:

(4.27)\begin{align} u_n(x) & =\int_{\Omega} \frac{G(x,y) f_n(y)}{(u_n + {1}/{n})^\gamma}\,{\rm d}y \leq \int_{\Omega} \frac{G(x,y) f(y)}{u_n^\gamma}\,{\rm d}y \nonumber\\ & \lesssim \mathbb{G}^{\Omega}\left[\frac{1}{\delta^{\gamma+\zeta}}\right] \lesssim \delta(x) \left\{ \begin{array}{@{}ll} 1 & \text{ if } \ \zeta + \gamma < 1, \\ \ln\left(\dfrac{\mathcal{D}_\Omega}{\delta(x)}\right) & \text{ if } \ \zeta + \gamma = 1. \end{array}\right. \end{align}

Using (4.27) and estimates in theorem 3.3, we obtain our first claim.

Case 2: $\zeta + \gamma > 1$.

For $n \in \mathbb {N}$ and $\kappa >0$, we define $h_{\epsilon _n}:= {1}/{((\delta +\epsilon _n^{1/\kappa })^{2-\kappa })} \in L^\infty (\Omega )$ with $\epsilon _n={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:

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

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 A.3, we obtain:

(4.28)\begin{equation} 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{equation}

and on the other hand for $0< \eta <1$, lemma 4.6 gives

(4.29)\begin{align} v_{\epsilon_n}(x)& = \mathbb{G}^\Omega\left[\frac{1}{(\delta+\epsilon_1^n)^{2-\kappa}}\right](x) \nonumber\\ & \gtrsim \left\{\begin{array}{@{}ll} \displaystyle \mathbb{G}^{\Omega}\left[\dfrac{1}{(\delta+\epsilon_1^n)^{2-\kappa}} \chi_{\Omega_{\eta}}\right](x) \gtrsim \left(\dfrac{1}{2}(\delta(x)+ \epsilon_1^n)^\kappa - \epsilon_n \right) & \text{ if } x \in \Omega_{{\eta}/{2}},\\ (\delta(x)+ \epsilon_1^n)^\kappa - \epsilon_n & \text{ if } x \in \Omega {\setminus} \Omega_{{\eta}/{2}}, \\ \end{array}\right.\nonumber\\ & \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{align}

Collecting the estimates in (4.28)–(4.29), there exists a constant $C_3, C_4>0$ independent of $n$ such that

(4.30)\begin{equation} 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. \end{equation}

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_{{\eta}/{2}} \]

and

\[ \left(\underline u^{\epsilon_n } + \epsilon_n\right) \leq 2 C_4 c_\eta (\delta + \epsilon_1^n)^\kappa + (1-C_4 c_\eta) \epsilon_n \quad \text{ in }\Omega. \]

If $2 C_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^{{1}/{(\gamma +1)}} (4C_4)^{-({\gamma }/{(\gamma +1)})}$, we have

\begin{align*} & (-\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}\nonumber\\ & \quad \leq f_{\epsilon_n}(x) \left(\underline u^{\epsilon_n } + \epsilon_n\right)^{-\gamma} \quad \text{in}\ \Omega, \end{align*}

where $\mathcal {G}_1$ is defined in (3.4).

If $2 C_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{align*} (-\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{align*}

Now, by considering both the cases and applying the weak comparison principle in $\Omega _{{\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 < {1}/{2}$ (by taking $\eta$ small enough) such that

(4.31)\begin{equation} C^1 (\delta + \epsilon_1^n)^\kappa - C^2 {\epsilon_n} \leq u_n \quad \text{in } \Omega. \end{equation}

Finally, by using the integral representation, we obtain:

(4.32)\begin{align} 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] \nonumber\\ & \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{align}

4.3.2 Sobolev regularity estimates

Lemma 4.8 Let $\gamma >0$, $\zeta \in [0,2)$ and $u_n$ be the weak solution of problem $(P_n)$. Then,

\begin{align*} & 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, \\ \mathfrak{L}^* & \text{ if } \ \zeta+ \gamma > 1, \end{array}\right. \\ & \quad \text{where} \ \mathfrak{L}^*:= \frac{\gamma + \zeta-1}{2-\zeta}. \end{align*}

Proof. Let $n \in \mathbb {N}$ and $u_n$ be the weak solution of 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, {1}/{n})$ and passing $\epsilon \to 0$ as in the proof of lemma 4.1(i), we obtain:

(4.33)\begin{align} \frac{4 \mathfrak{L}}{(\mathfrak{L} +1)^2}\int_{\Omega} |\nabla u_n^{{(\mathfrak{L} +1)}/{2}}|^2\,{\rm d}x & + \frac{2 C(N,s) \mathfrak{L}}{(\mathfrak{L} +1)^2} \int_{\mathbb{R}^N} \nonumber\\ & \quad\int_{\mathbb{R}^N} \frac{((u_n^{{(\mathfrak{L} +1)}/{2}}(x)- u_n^{{(\mathfrak{L} +1)}/{2}}(y))^2}{|x-y|^{N+2s}}\,{\rm d}x\,{\rm d}y\nonumber\\ & \leq \int_{\Omega} f_n(x) u_n^{\mathfrak{L} - \gamma}\,{\rm d}x:= \mathfrak{G}(u_n). \end{align}

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$.

(4.34)\begin{align} \mathfrak{G}(u_n) & \leq \int_{\Omega} f_n(x) u_n^{\mathfrak{L} - \gamma}\,{\rm d}x \leq \mathcal{G}_2 \int_{\Omega} \frac{u_n^{\mathfrak{L} - \gamma}}{\delta^\zeta(x)}\,{\rm d}x \lesssim \int_{\Omega} u_n^{\mathfrak{L} - \gamma -\zeta}\,{\rm d}x \nonumber\\ & \lesssim \int_{\Omega} \delta^{\mathfrak{L} - \gamma -\zeta}(x)\,{\rm d}x \leq C \quad \text{if} \quad \mathfrak{L} > \gamma + \zeta-1. \end{align}

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

(4.35)\begin{align} \mathfrak{G}(u_n) & \leq \int_{\Omega} f_n(x) u_n^{\mathfrak{L} - \gamma}\,{\rm d}x \leq \mathcal{G}_2 \int_{\Omega} \frac{u_n^{\mathfrak{L} - \gamma}}{\delta^\zeta(x)}\,{\rm d}x \lesssim \int_{\Omega} \ln^\zeta\left(\frac{\mathcal{D}_\Omega}{\delta(x)}\right) u_n^{\mathfrak{L} - 1}\,{\rm d}x \nonumber\\ & \lesssim \left\{ \begin{array}{@{}ll} \int_{\Omega} \ln^{\mathfrak{L}-\gamma}\left(\dfrac{\mathcal{D}_\Omega}{\delta(x)}\right) \delta^{\mathfrak{L} - 1}(x)\,{\rm d}x & \text{ if } \ \mathfrak{L} \geq 1, \\ \int_{\Omega} \ln^\zeta\left(\dfrac{\mathcal{D}_\Omega}{\delta(x)}\right) \delta^{\mathfrak{L} - 1}(x)\,{\rm d}x & \text{ if }\ \mathfrak{L} < 1, \end{array}\right.\nonumber\\ & \lesssim \int_{\Omega} \delta^{\mathfrak{L} - 1-\chi}(x)\,{\rm d}x \leq C \quad \text{if}\ \mathfrak{L} > \chi. \end{align}

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

(4.36)\begin{align} \mathfrak{G}(u_n) & \leq \int_{\Omega} f_n(x) u_n^{\mathfrak{L} - \gamma}\,{\rm d}x \leq \mathcal{G}_2 \int_{\Omega} \frac{u_n^{\mathfrak{L} - \gamma}}{\delta^\zeta(x)}\,{\rm d}x = \mathcal{G}_2 \nonumber\\ & \quad \int_{\Omega} u_n^{\mathfrak{L} - \gamma} \left(\delta^{({2-\zeta})/({(\gamma+1)})}(x)\right)^{({-\zeta(\gamma+1)})/{(2-\zeta)}}\,{\rm d}x\nonumber\\ & \lesssim \int_{\Omega} u_n^{\mathfrak{L} - \gamma -(({\zeta(\gamma+1)})/{(2-\zeta)}})\,{\rm d}x \lesssim \int_{\Omega} \delta^{\left(\mathfrak{L} - \gamma\right) ({(2-\zeta)}/({(\gamma+1)} -\zeta}))(x)\,{\rm d}x \nonumber\\ & \quad\leq C \quad \text{if}\quad \mathfrak{L} > \frac{\gamma + \zeta-1}{2-\zeta}. \end{align}

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

\begin{align*} & \int_{\Omega} |\nabla u_n^{{(\mathfrak{L} +1)}/{2}}|^2\,{\rm d}x \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}\nonumber\\ & \quad > \left\{ \begin{array}{ll} 0 & \text{ if } \ \zeta+ \gamma \leq 1, \\ \dfrac{\gamma + \zeta-1}{2-\zeta} & \text{ if} \ \zeta+ \gamma > 1. \end{array}\right. \end{align*}

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 problem $(P_n)$ in the sense that

(5.1)\begin{align} & \int_{\Omega} \nabla u_n \cdot \nabla \psi\,{\rm d}x + \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}}\,{\rm d}x\,{\rm d}y \nonumber\\ & \quad = \int_{\Omega} \frac{f_n(x)}{(u_n + {1}/{n})^\gamma} \psi\,{\rm d}x \quad \forall \psi \in H_0^1(\Omega). \end{align}

From lemma 4.1(iv), we know that $u_n$ is uniformly bounded in $W_0^{1,q}(\Omega )$ with $q:= ({Nr(1+\gamma )})/({N-r(1-\gamma )}).$ Then, there exists a $u \in W_0^{1, q}(\Omega )$ such that

(5.2)\begin{equation} 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{equation}

Since, for any $\psi \in W_0^{1,q'}(\Omega )$:

\begin{align*} & W_0^{1,q}(\Omega) \ni f \mapsto \int_{\Omega} \nabla f \cdot \nabla \psi\,{\rm d}x + \frac{C(N,s)}{2} \int_{\mathbb{R}^N}\\ & \quad \int_{\mathbb{R}^N} \frac{(f(x)- f(y))(\psi(x)- \psi(y))}{|x-y|^{N+2s}}\,{\rm d}x\,{\rm d}y \end{align*}

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

(5.3)\begin{align} & \int_{\Omega} \nabla u_n \cdot \nabla \psi\,{\rm d}x + \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}}\,{\rm d}x\,{\rm d}y \nonumber\\ & \quad \to \int_{\Omega} \nabla u \cdot \nabla \psi\,{\rm d}x + \frac{C(N,s)}{2} \nonumber\\ & \quad \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{(u(x)- u(y))(\psi(x)- \psi(y))}{|x-y|^{N+2s}}\,{\rm d}x\,{\rm d}y \quad \text{as} \ n \to \infty \end{align}

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:

(5.4)\begin{align} & \int_{\Omega} \frac{f_n(x)}{(u_n + {1}/{n})^\gamma} \psi\,{\rm d}x \to \int_{\Omega} \frac{f(x)}{u^\gamma} \psi\,{\rm d}x \quad \text{as} \ n\nonumber\\ & \quad \to \infty \quad \text{for }\psi \in L^{r'}(\Omega)\ \text{with}\ \operatorname{supp}(\psi) \Subset \Omega. \end{align}

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 problem $(P_n)$ and satisfies (5.1). From lemma 4.1, we know that $u_n^{{(\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{align*} & \int_{\omega} |\nabla u_n|^2\,{\rm d}x \leq C^{-(\mathfrak{S}_r-1)} \int_{\omega} u_n^{(\mathfrak{S}_r-1)} |\nabla u_n|^2\,{\rm d}x \leq \frac{4 C^{-(\mathfrak{S}_r-1)}}{(\mathfrak{S}_r +1)} \nonumber\\ & \quad \int_{\omega} |\nabla u_n^{{(\mathfrak{S}_r+1)}/{2}}|^2\,{\rm d}x \leq C_1 \end{align*}

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

(5.5)\begin{equation} 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{equation}

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$:

(5.6)\begin{align} & \int_{\Omega} \nabla u_n \cdot \nabla \psi\,{\rm d}x + \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}}\,{\rm d}x\,{\rm d}y\nonumber\\ & \quad \to \int_{\Omega} \nabla u \cdot \nabla \psi\,{\rm d}x + \frac{C(N,s)}{2} \int_{\mathbb{R}^N}\nonumber\\ & \quad \int_{\mathbb{R}^N} \frac{(u(x)- u(y))(\psi(x)- \psi(y))}{|x-y|^{N+2s}}\,{\rm d}x\,{\rm d}y \quad \text{as}\ n \to \infty. \end{align}

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 the limit $n \to \infty$ using (5.6), we obtain our claim.

5.3 Proof of theorem 2.3

Let $u$ be the weak solution of problem (1.2) obtained from the approximable solution $u_n$ of problem $(P_n)$.

1. (i) Let $r \in [1, {N}/{2})$ and $N > 2$. Then, from lemma 4.1, we know that $u_n^{{(\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^{{(\mathfrak {S}_r+1)}/{2}} \rightharpoonup v \ \text {in} \ H_0^1(\Omega )$ and by Sobolev's embedding, we get $u_n^{{(\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^{{(\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={N}/{2}$ and $N \geq 2$. Again from lemma 4.1, we imply that $\mathfrak {H} ({u_n}/{2})$ in case of $0 \leq \gamma \leq 1$ and $\mathfrak {D} ({u_n}/{2})$ 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{align*} & 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) \nonumber\\ & \quad \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)\nonumber\\ & \quad \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{align*}
   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> {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 the claim follows using the embedding $L^r(\Omega ) \hookrightarrow L^1(\Omega )$ for $1 \leq r \leq \infty$ and by repeating the proof of theorems 2.1 and 2.2 with $\mathfrak {S}_r=1$.

5.4 Proof of theorem 2.4

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

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

Indeed, in this case, we have

\[ \|u^{{(\mathfrak{S}+1)}/{2}}\|_{H_0^1(\Omega)}^2 \geq C\int_{\Omega} \frac{|u^{{(\mathfrak{S}+1)}/{2}}(x)|^2}{|d(x)|^2}\,{\rm d}x \geq C\int_{\Omega} d^{\mathfrak{S}-1}(x)\,{\rm d}x =\infty \]

and hence, we deduce $u^{{(\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 problems (1.2)–(1.3) for datum $f$ and $g$, respectively, obtained in theorems 2.1 and 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 nonsingular 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{align*} & \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)_+\,{\rm d}x\nonumber\\ & \quad + \int_{\Omega} \frac{g_n}{\left(v_n+ {1}/{n}\right)^\gamma} ([(u_n-v_n)_+{+} \epsilon]^{\mathfrak{S}_r} -\epsilon^{\mathfrak{S}_r})\,{\rm d}x \\ & \qquad+ \frac{C(N,s)}{2} \int_{\mathbb{R}^N}\nonumber\\ & \quad \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}}\,{\rm d}x \,{\rm d}y \\ & \quad= \int_{\Omega} \frac{f_n}{\left(u_n+ {1}/{n}\right)^\gamma} ([(u_n-v_n)_+{+} \epsilon]^{\mathfrak{S}_r} -\epsilon^{\mathfrak{S}_r})\,{\rm d}x \end{align*}

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 < {1}/{n}$, and using (4.3) and (4.6), we have the dominating function:

\begin{align*} & \frac{f_n}{\left(u_n+ {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}\nonumber\\ & \quad \leq f (u + 1)^{\mathfrak{S}_r-\gamma} \leq \frac{|f|^r}{r} + \frac{(r-1)}{r} |u+1|^{\sigma_r}. \end{align*}

Therefore, we have

(5.7)\begin{align} 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)_+ ~{\rm d}x \nonumber\\ & \quad + \int_{\Omega} \frac{g_n}{\left(v_n+ {1}/{n}\right)^\gamma} (u_n-v_n)_+^{\mathfrak{S}_r}\,{\rm d}x\nonumber\\ & \quad+ \frac{C(N,s)}{2} \int_{\mathbb{R}^N}\nonumber\\ & \quad \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}}\,{\rm d}x\,{\rm d}y \nonumber\\ & \leq \int_{\Omega} \frac{f_n}{\left(u_n+ {1}/{n}\right)^\gamma} (u_n-v_n)_+^{\mathfrak{S}_r}\,{\rm d}x: =I_4 \end{align}

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):

(5.8)\begin{equation} 0 \leq \frac{f_n}{\left(u_n+ {1}/{n}\right)^\gamma} (u_n-v_n)_+^{\mathfrak{S}_r} \leq f u^{\mathfrak{S}_r-\gamma} \leq \frac{|f|^r}{r} + \frac{(r-1)}{r} |u|^{\sigma_r}, \end{equation}

and the fact that $u_n \to u$ and $v_n \to v$ a.e. in $\Omega$, Lebesgue-dominated convergence theorem implies:

(5.9)\begin{equation} I_4 \to \int_{\Omega} \frac{f}{u^\gamma} (u-v)_+^{\mathfrak{S}_r} ~{\rm d}x \end{equation}

Estimate for $I_2 + I_3$: It is easy to see that

(5.10)\begin{align} & ((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 \nonumber\\ & \quad\quad \text{for} \quad (x,y) \in \mathbb{R}^N \times \mathbb{R}^N. \end{align}

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

(5.11)\begin{equation} \int_{\Omega} \frac{g}{v^\gamma} (u-v)_+^{\mathfrak{S}_r}\,{\rm d}x \leq \lim\inf_{n \to \infty} I_2 + I_3 \end{equation}

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

\[ \frac{4 \mathfrak{S}_r}{(\mathfrak{S}_r\!+\!1)^2} \int_{\Omega} \left|\nabla \left[(u_n\!-\!v_n)_+{+} \delta\right]^{{(\mathfrak{S}_r+1)}/{2}}\right|^2\,{\rm d}x\!=\! \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}}\,{\rm d}x \!\leq\! I_1. \]

Together with the estimate in (5.8) and (5.10), we obtain $[(u_n-v_n)_+ + \delta ]^{{(\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]^{{(\mathfrak{S}_r+1)}/{2}} \rightharpoonup \left[(u-v)_+{+} \delta\right]^{{(\mathfrak{S}_r+1)}/{2}} \ \text{in} \ H_0^1(\Omega) \]

and the weak lower semicontinuity of the norm gives

\begin{align*} \mathfrak{S}_r \int_{\Omega} \frac{\left|\nabla (u-v)_+\right|^2}{\left((u-v)_+{+} \delta\right)^{1-\mathfrak{S}_r}}\,{\rm d}x & = \frac{4 \mathfrak{S}_r}{(\mathfrak{S}_r+1)^2} \int_{\Omega} \left|\nabla \left[(u-v)_+{+} \delta\right]^{{(\mathfrak{S}_r+1)}/{2}}\right|^2\,{\rm d}x\\ & \leq \lim\inf_{n \to \infty} \frac{4 \mathfrak{S}_r}{(\mathfrak{S}_r+1)^2} \nonumber\\ & \quad\int_{\Omega} \left|\nabla \left[(u_n-v_n)_+{+} \delta\right]^{{(\mathfrak{S}_r+1)}/{2}}\right|^2\,{\rm d}x \leq \lim\inf_{n \to \infty} I_1. \end{align*}

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

(5.12)\begin{equation} \mathfrak{S}_r \int_{\Omega} \frac{\left|\nabla (u-v)_+\right|^2}{(u-v)_+^{1-\mathfrak{S}_r}}\,{\rm d}x \leq \lim\inf_{n \to \infty} I_1. \end{equation}

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:

(5.13)\begin{align} \frac{4 \mathfrak{S}_r}{(\mathfrak{S}_r+1)^2} \int_{\Omega} |\nabla (u-v)_+^{{(\mathfrak{S}_r+1)}/{2}}|^2\,{\rm d}x & \leq \int_{\Omega} \left(\frac{f(x)}{u^\gamma} - \frac{g(x)}{v^\gamma} \right) (u-v)_+^{\mathfrak{S}_r}\,{\rm d}x\nonumber\\ & \leq \int_{\Omega} (f(x)-g(x)) \frac{(u-v)_+^{\mathfrak{S}_r}}{u^\gamma}\,{\rm d}x\nonumber\\ & \leq \int_{\Omega} (f(x)-g(x)) (u-v)_+^{\mathfrak{S}_r-\gamma}\,{\rm d}x. \end{align}

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{align*} & \frac{4 \mathfrak{S}_r}{(\mathfrak{S}_r+1)^2} \int_{\Omega} |\nabla (u-v)_+^{{(\mathfrak{S}_r+1)}/{2}}|^2\,{\rm d}x \nonumber\\ & \quad\leq \|f-g\|_{L^r(\Omega_1)} \left(\int_{\Omega} \left((u-v)_+^{{(\mathfrak{S}_r +1)}/{2}}\right)^{{2N}/{(N-2)}}\,{\rm d}x\right)^{{1}/{r'}}\\ & \quad \leq \|f-g\|_{L^r(\Omega_1)} \left( S(N) \int_{\Omega} |\nabla (u-v)_+^{{(\mathfrak{S}_r +1)}/{2}}|^2\,{\rm d}x \right)^{{N}/{(r'(N-2))}} \nonumber\\ & \qquad \text{where} \ \Omega_1:= \{u \geq v\} \end{align*}

which further implies

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

and analogously

\[ \|\nabla (v-u)_+^{{(\mathfrak{S}_r+1)}/{2}}\|^2_{L^2(\Omega)} \leq C \|f-g\|_{L^r(\Omega {\setminus} \Omega_1)}^{{(r(N-2))}/{(N-2r)}}. \]

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

5.6 Proof of theorem 2.10

Let $\zeta +\gamma \leq 1$ and $u_n$ be the weak solution of problem $(P_n)$ in the sense that

(5.14)\begin{align} & \int_{\Omega} \nabla u_n \cdot \nabla \psi\,{\rm d}x + \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}}\,{\rm d}x\,{\rm d}y \nonumber\\ & \quad = \int_{\Omega} \frac{f_n(x)}{(u_n + {1}/{n})^\gamma} \psi\,{\rm d}x \quad \forall \psi \in H_0^1(\Omega). \end{align}

From lemmas 4.8 and 4.7, we know that, for any $\mathfrak {L}>0$, $u_n^{{(\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:

(5.15)\begin{equation} 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{equation}

Using Hardy's inequality, for any $\psi \in H_0^1(\Omega )$, ${f_n \psi }/({(u_n+{1}/{n})^\gamma })$ is uniformly integrable. Indeed,

\[ \int_{\Omega} \frac{f_n(x)}{(u_n \!+\! {1}/{n})^\gamma} \psi\,{\rm d}x \!\leq\! \int_{\Omega} \frac{f(x)}{u_n^\gamma} \psi\,{\rm d}x \lesssim\! \int_{\Omega} \delta^{1-\gamma-\zeta}(x) \frac{\psi}{\delta(x)}\,{\rm d}x \!\lesssim\! \|\frac{\psi}{\delta}\|_{L^2(\Omega)} \!\leq\! \|\psi\|_{H_0^1(\Omega)}. \]

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

\begin{align*} & \int_{\Omega} \nabla u\cdot \nabla \psi\,{\rm d}x + \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}}\,{\rm d}x\,{\rm d}y\nonumber\\ & \quad = \int_{\Omega} \frac{f(x)}{u^\gamma} \psi\,{\rm d}x \quad \forall \psi \in H_0^1(\Omega). \end{align*}

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.7). Passing to the limit $n \to \infty$ gives that $u$ is a minimal solution.

5.7 Proof of theorem 2.11

Let $\zeta +\gamma > 1$ and $u_n$ be the weak solution of problem $(P_n)$. From lemmas 4.8 and 3.2, we have, for any $\mathfrak {L}>\mathfrak {L}^*$, $u_n^{{(\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,

(5.16)\begin{align} & \int_{\omega} |\nabla u_n|^2\,{\rm d}x \leq C^{-(\mathfrak{L}-1)}(\omega) \int_{\omega} u_n^{(\mathfrak{L}-1)} |\nabla u_n|^2\,{\rm d}x \leq \frac{4 C^{-(\mathfrak{L}-1)}}{(\mathfrak{L} +1)}\nonumber\\ & \quad \int_{\omega} |\nabla u_n^{{(\mathfrak{L}+1)}/{2}}|^2\,{\rm d}x \leq C_0 \end{align}

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

(5.17)\begin{equation} 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{equation}

Using Hardy's inequality, we have for any $\psi \in H_0^1(\Omega )$ with $\operatorname {supp}(\psi ) \Subset \Omega$, ${f_n \psi }/({(u_n+{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{align*} & \int_{\Omega} \nabla u\cdot \nabla \psi\,{\rm d}x + \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}}\,{\rm d}x\,{\rm d}y \nonumber\\ & \quad= \int_{\Omega} \frac{f(x)}{u^\gamma} \psi\,{\rm d}x. \end{align*}

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.10, we obtain our claim.

5.8 Proof of theorem 2.12

From theorems 2.10, 2.11, lemma 4.7 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, \\ \ln\left(\dfrac{\mathcal{D}_\Omega}{\delta(x)}\right) & \text{ if } \ \zeta + \gamma = 1, \end{array}\right. \]

and

\[ C^1\delta^{{(2-\zeta)}/{(\gamma+1)}}(x)\leq u(x) \leq C^2 \delta^{{(2-\zeta)}/{(\gamma+1)}}(x). \]

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

(5.18)\begin{equation} \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 & \text{ in } \Omega, \\ w_\Xi = 0 & \text{ in } \ \mathbb{R}^N {\setminus} \Omega, \end{aligned}\right. \end{equation}

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

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

and there exists $c_1, c_2>0$ such that

(5.20)\begin{equation} 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. \end{equation}

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

\begin{align*} & (-\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) \nonumber\\ & \quad = (-\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) \end{align*}

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{align*} & \frac{(c_1 d_0)^{1+\gamma^2 + \dots + \gamma^{2j}}}{(d_1c_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) \\ & \quad \leq \frac{(c_1 d_0)^{1+\gamma^2 + \dots + \gamma^{2j}}}{(d_1c_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{align*}

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

\[ d_2 \delta(x) \ln^{{1}/{(1+\gamma)}}\left(\frac{\mathcal{D}_\Omega}{\delta(x)}\right) \leq u(x) \leq d_3 \delta(x) \ln^{{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 ={1}/{2}$ in (5.18) and (5.19), there exists $c_3, c_4>0$ such that

(5.21)\begin{equation} c_3 \delta(x) \ln^{{1}/{2}}\left(\frac{\mathcal{D}_{\Omega}}{\delta(x)}\right) \leq w_{{1}/{2}}(x) \leq c_4 \delta(x) \ln^{{1}/{2}}\left(\frac{\mathcal{D}_{\Omega}}{\delta(x)}\right) \quad \forall\ x \in \Omega, \end{equation}

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

\begin{align*} & (-\Delta) (C w_{{1}/{2}}) + (-\Delta)^s (C w_{{1}/{2}}) \geq \frac{1}{C w_{{1}/{2}}} \ \text{and} \ (-\Delta) \left(\frac{w_{{1}/{2}}}{C}\right) \nonumber\\ & \quad + (-\Delta)^s \left(\frac{w_{{1}/{2}}}{C}\right) \leq \frac{C}{ w_{{1}/{2}}}. \end{align*}

Therefore, again by comparison principle, we obtain:

\[ \frac{c_3}{C} \delta(x) \ln^{{1}/{2}}\left(\frac{\mathcal{D}_{\Omega}}{\delta(x)}\right) \leq u(x) \leq \frac{c_4}{C} \delta(x) \ln^{{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 behaviour of the weak solution. Precisely, if

\[ \mathfrak{L} \leq \left\{\begin{array}{@{}ll} 0 & \text{ if} \ \zeta+ \gamma \leq 1,\\ \mathfrak{L}^* & \text{ if} \ \zeta+ \gamma > 1, \end{array} \right. \quad \text{then} \quad u^{{(\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^{{(\mathfrak{L}+1)}/{2}}\|_{H_0^1(\Omega)}^2 \geq C\int_{\Omega} \frac{|u^{{(\mathfrak{L}+1)}/{2}}(x)|^2}{|d(x)|^2}\,{\rm d}x \geq C\int_{\Omega} d^{\mathfrak{L}-1}(x)\,{\rm d}x =\infty. \]

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

\[ \|u^{{(\mathfrak{L}+1)}/{2}}\|_{H_0^1(\Omega)}^2 \geq C\int_{\Omega} \frac{|u^{{(\mathfrak{L}+1)}/{2}}(x)|^2}{|d(x)|^2}\,{\rm d}x \geq C\int_{\Omega} d^{({((\mathfrak{L}+1)(2-\zeta))}/{(1+\gamma)})-2}(x)\,{\rm d}x =\infty \]

and we deduce $u^{{(\mathfrak {L}+1)}/{2}} \notin H_0^{1}(\Omega ).$

5.9 Proof of theorem 2.13

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 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

(5.22)\begin{align} & \int_{\Omega} \nabla v_n \cdot \nabla \phi\,{\rm d}x + \frac{C(N,s)}{2} \iint_{\mathbb{R}^{2N}} \frac{(v_n(x)- v_n(y)) (\phi(x)- \phi(y))}{|x-y|^{N+2s}}\,{\rm d}x\,{\rm d}y \nonumber\\ & \quad= \int_{\Omega} \frac{\Gamma g_n(x)}{(v_n+ {1}/{n})^\gamma} \phi\,{\rm d}x \end{align}

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 {\theta }/{2}$ in $\Omega _\eta$. Since $w \geq 0$, then $u:= v_n- w- \theta \leq - ({\theta }/{2}) <0$ in $\Omega _\eta$ and

\[ \text{supp}({u}^+) \subset \text{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

(5.23)\begin{align} & \int_{\Omega} \nabla v_n \cdot \nabla u^+\,{\rm d}x + \frac{C(N,s)}{2} \iint_{\mathbb{R}^{2N}} \frac{(v_n(x)- v_n(y)) (u^+(x)- u^+(y))}{|x-y|^{N+2s}}\,{\rm d}x\,{\rm d}y\nonumber\\ & \quad = \int_{\Omega} \frac{\Gamma g_n(x)}{(v_n+ \epsilon)^\gamma} u^+\,{\rm d}x \leq \int_{\Omega} \frac{\Gamma g_n(x)}{v_n^\gamma} u^+\,{\rm d}x. \end{align}

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

(5.24)\begin{align} & \int_{\Omega} \nabla w \cdot \nabla u^+\,{\rm d}x + \frac{C(N,s)}{2} \iint_{\mathbb{R}^{2N}} \frac{(w(x)- w(y)) (u^+(x)- u^+(y))}{|x-y|^{N+2s}}\,{\rm d}x\,{\rm d}y \nonumber\\ & \quad = \int_{\Omega} \frac{f(x)}{w^\gamma} w^+\,{\rm d}x \geq \int_{\Omega} \frac{ \Gamma g_n(x)}{w^\gamma} w^+\,{\rm d}x. \end{align}

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

\begin{align*} 0 & \leq \int_{\Omega} (\nabla v_n - \nabla w) \cdot \nabla u^+\,{\rm d}x + \frac{C(N,s)}{2} \nonumber\\ & \quad \iint_{\mathbb{R}^{2N}} \frac{((v_n(x)- v_n(y))-(w(x)- w(y))) (u^+(x)- u^+(y))}{|x-y|^{N+2s}}\,{\rm d}x\,{\rm d}y\\ & \leq \int_{\Omega} \left(\frac{\Gamma g_n(x)}{v_n^\gamma} - \frac{\Gamma g_n(x)}{w^\gamma} \right) u^+\,{\rm d}x \leq 0, \end{align*}

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.7, we obtain:

\[ c_1 (\delta(x) + n^{{(-(1+\gamma))}/{((2-\zeta_0))}})^{{(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:

\begin{align*} & \int_{\Omega} \left(c_1 (\delta(x) + n^{{(-(1+\gamma))}/{((2-\zeta_0))}})^{{(2-\zeta_0)}/{(\gamma+1)}} - c_2 n^{{-}1}\right)^{2\mathfrak{L}_0} d^{{-}2}(x)\,{\rm d}x \nonumber\\ & \quad\leq \int_{\Omega} \left|\frac{w^{\mathfrak{L}_0}}{{\rm d}(x)}\right|^2\,{\rm d}x < \infty. \end{align*}

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.