# Existence and global behaviour of weak solutions to a doubly nonlinear evolution equation involving fractional p-Laplacian equation

Jacques Giacomoni, Abdelhamid Gouasmia, Abdelhafid Mokrane

► To cite this version:

Jacques Giacomoni, Abdelhamid Gouasmia, Abdelhafid Mokrane. Existence and global behaviour of weak solutions to a doubly nonlinear evolution equation involving fractional p-Laplacian equation. Electronic Journal of Differential Equations, 2021, 2021 (01-104), pp.09. 10.58997/ejde.2021.09 . hal-03150699

HAL Id: hal-03150699

https://hal.science/hal-03150699v1

Submitted on 12 Jul 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# EXISTENCE AND GLOBAL BEHAVIOR OF WEAK SOLUTIONS TO A DOUBLY NONLINEAR EVOLUTION FRACTIONAL p-LAPLACIAN EQUATION

JACQUES GIACOMONI, ABDELHAMID GOUASMIA, ABDELHAFID MOKRANE

**ABSTRACT.** In this article, we study a class of doubly nonlinear parabolic problems involving the fractional *p*-Laplace operator. For this problem, we discuss existence, uniqueness and regularity of the weak solutions by using the time-discretization method and monotone arguments. For global weak solutions, we also prove stabilization results by using the accretivity of a suitable associated operator. This property is strongly linked to the Picone identity that provides further a weak comparison principle, barrier estimates and uniqueness of the stationary positive weak solution.

## 1. INTRODUCTION AND STATEMENT OF MAIN RESULTS

Let $1 < q \le p < \infty$, $0 < s < 1$, $Q_T := (0, T) \times \Omega$, where $\Omega \subset \mathbb{R}^N$, with $N > sp$, is an open bounded domain with $C^{1,1}$ boundary. $\Gamma_T := (0, T) \times \partial\Omega$ denotes the lateral boundary of the cylinder $Q_T$. In this work, we deal with the existence, uniqueness and other qualitative properties of the weak solution to the following doubly nonlinear parabolic equation:

$$ 
\begin{array}{l}
\displaystyle \frac{q}{2q-1} \partial_t (u^{2q-1}) + (-\Delta)_p^s u = f(x, u) + h(t, x) u^{q-1} \quad \text{in } Q_T, \\
\displaystyle \qquad u > 0 \quad \text{in } Q_T, \\
\displaystyle \qquad u = 0 \quad \text{on } \Gamma_T, \\
\displaystyle \qquad u(0, \cdot) = u_0 \quad \text{in } \Omega.
\end{array} 
\quad (1.1) $$

Here $(-\Delta)_p^s u$ is the fractional $p$-Laplace operator, defined for $1 < p < \infty$, as

$$ (-\Delta)_p^s u(x) := 2 \operatorname{P.V.} \int_{\mathbb{R}^N} \frac{|u(x) - u(y)|^{p-2} (u(x) - u(y))}{|x - y|^{N+sp}} dy, $$

where P.V. denotes the Cauchy principal value. We refer to [21, 29, 38] for the main properties of this nonlinear fractional elliptic operator.

Throughout this article we assume the following hypothesis:

(H1) $f: \bar{\Omega} \times \mathbb{R}^+ \to \mathbb{R}^+$ is a continuous function, such that $f(x, 0) \equiv 0$ and $f$ is positive on $\Omega \times \mathbb{R}^+ \setminus \{0\}$.

(H2) For a.e. $x \in \Omega$, $s \mapsto \frac{f(x,s)}{s^{q-1}}$ is non-increasing in $\mathbb{R}^+ \setminus \{0\}$.

<footnote n="2010">2010 Mathematics Subject Classification. 35B40, 35K59, 35K55, 35K10, 35R11.</footnote>

<footnote>Key words and phrases. Fractional p-Laplace equation; doubly nonlinear evolution equation; Picone identity; stabilization; nonlinear semi-group theory.</footnote>

<footnote>©2021 Texas State University.</footnote>

<footnote>Submitted June 6, 2020. Published February 23, 2021.</footnote>
(H3) If $q = p$, $s \mapsto \frac{f(x,s)}{s^{p-1}}$ is decreasing in $\mathbb{R}^+ \setminus \{0\}$ for a.e. $x \in \Omega$ and
$$ \lim_{r \to +\infty} \frac{f(x,r)}{r^{p-1}} = 0 $$ uniformly in $x \in \Omega$.

(H4) There exists $\underline{h} \in L^\infty(\Omega) \setminus \{0\}$, $\underline{h} \ge 0$ such that $h(t, x) \ge \underline{h}(x)$ a.e. in $Q_T$.

(H5) If $q = p$,

$$ \|h\|_{L^\infty(Q_T)} < \lambda_{1,s,p} := \inf_{\phi \in W_0^{s,p}(\Omega) \setminus \{0\}} \frac{\|\phi\|_{W_0^{s,p}(\Omega)}^p}{\|\phi\|_{L^p(\Omega)}^p}. $$

(H6) If $q = p$, $\underline{h}$, $f$ fulfills the condition

$$ \inf_{x \in \Omega} \left( \underline{h}(x) + \lim_{s \to 0^+} \frac{f(x, s)}{s^{p-1}} \right) > \lambda_{1,s,p}. $$

1.1. **State of the art.** The study of nonlocal elliptic operators arouse more and more interest in mathematical modeling, see e.g. [8, 11, 12, 14, 27, 34, 42] and the references cited therein. Concerning the investigation on parabolic equations involving nonlocal operators, we refer to [1, 5, 15, 16, 18, 19, 24, 25, 30, 31, 32, 33, 35, 37, 38, 39, 41] without giving an exhaustive list. These types of operators arise in several contexts: in finance, physics, fluid dynamics, image processing and in various fields like continuum mechanics, stochastic processes of Lévy type, phase transitions, population dynamics, optimal control and game theory, see for further discussion [15, 17, 21, 29, 39] and the references therein. In particular [15] shows some non-local diffusion models coming from game theory. In connection to our doubly nonlinear problem (1.1), [37] shows different methods (entropy method and contraction semi-group theory) two evolution models of flows in porous media involving fractional operators:

* The first model is based on Darcy's law and is given by

$$ \begin{align*} \partial_t u &= \nabla \cdot (u \nabla P) \quad \text{in } (0, \infty) \times \mathbb{R}^N, \\ P &= (-\Delta)^{-s} u \quad \text{in } (0, \infty) \times \mathbb{R}^N, \\ u(0, x) &= u_0(x) \quad \text{in } \mathbb{R}^N, \end{align*} $$

where $u$ is the particle density of the fluid, $P$ is the pressure and $(-\Delta)^{-s}$ is the inverse of the fractional Laplace operator (i.e. $p = 2$). The initial data $u_0$ is a nonnegative, bounded and integrable function in $\mathbb{R}^N$ (see also [13] for further explanations).

* The second model in analogy to classical models of transport through porous media (see [22]) is described in the non local case by

$$ \partial_t u + (-\Delta)^s (u^m) = 0. \qquad (1.2) $$

For $s \to 1^-$ and $m = 1$, the limiting model (1.2) is the well known heat equation. Furthermore if $m > 1$, (1.2) is known as the porous media equation (PME for short) whereas in case $m < 1$ it is referred as the fast diffusion equation (FDE for short). Existence and global behaviour of solutions are described in [37] for the two types of models. We refer again to [39] for further explanations about the physical background and the adequacy of nonlocal diffusion operators (see also [19] for related issues). The paper [18] deals with the problem (1.2) in the special case $s = \frac{1}{2}$, and $p = 2$ and investigates the local existence, uniqueness and regularity of the weak solution. We highlight here that few results are available about the parabolic equation involving fractional p-Laplacian operator in contrast with the stationary elliptic equation.
In [25], considering the more general case $1 < p < \infty$, authors obtain the existence, uniqueness, and regularity of the weak solution to the fractional reaction diffusion equation

$$ 
\begin{array}{l}
\partial_t u + (-\Delta)_p^s u + g(x, u) = f(x, u) \quad \text{in } Q_T; \\
u = 0 \quad \text{in } \mathbb{R}^N \setminus Q_T; \\
u(0, \cdot) = u_0 \quad \text{in } \mathbb{R}^N,
\end{array}
\quad (1.3) $$

with $f$ and $g$, satisfying suitable growth and homogeneity conditions. In addition, the authors prove that global solutions converge to the unique positive stationary solution as $t \to \infty$. Previously, [1] has dealt with the case where the nonlinearity $f$ depends only on $x$ and $t$ and have established the existence and some properties of nonnegative entropy solutions. In [24], the authors have studied (1.3), under similar conditions about $f$ and $g(x, u) := -|u(t, x)|^{q-2}u(t, x)$, with $q \ge 2$. They prove the existence of locally-defined strong solutions to the problem with any initial data $u_0 \in L^r(\Omega)$ and $r \ge 2$. They also investigate the occurrence of finite time blow up behavior. In [30, 38] the results about existence, uniqueness and $T$-accretivity in $L^1$ of strong solutions to the fractional $p$-Laplacian heat equation with Dirichlet or Neumann boundary conditions, are obtained through the theory of nonlinear accretive operators. The asymptotic decay of solutions and the study of asymptotic models as $p \to 1^+$ are also investigated. In [26], authors extend the results obtained in [4] in case of singular nonlinearities and fractional diffusion. We refer the reader to [28, 33, 36, 40, 41] for further investigations of above issues.

The aim of this article is to discuss similar issues about local existence, uniqueness, regularity and global behavior of solutions to the doubly nonlinear and non local equation (1.1). Up to our knowledge, (1.1) which covers several PME and FDE models in the fractional setting has not been investigated in the literature. By using the semi-discretization in time method applied to an auxiliary evolution problem, we prove the local existence of weak energy solutions. The uniqueness of weak solutions are obtained via the fractional version of the Picone identity (see below) which leads to a new comparison principle and $T$-accretivity of an associated operator in $L^2$. Using the comparison principle, we also prove the existence of barrier functions from which we derive that weak solutions are global. We then show that weak solutions converge to the unique non trivial stationary solution as $t \to \infty$. To achieve this goal, our approach borrows techniques from the contraction semi-group theory.

## 1.2. Preliminaries and functional setting.
First, we recall some notation which will be used throughout the paper. Considering a measurable function $u : \mathbb{R}^N \to \mathbb{R}$, we adopt

• Let $p \in [1; +\infty[$, the norm in the space $L^p(\Omega)$ is denoted by

$$ \|u\|_{L^p(\Omega)} := \left( \int_{\Omega} |u|^p dx \right)^{1/p}. $$

• Set $0 < s < 1$ and $p > 1$, we recall that the fractional Sobolev space $W^{s,p}(\mathbb{R}^N)$ is defined as

$$ W^{s,p}(\mathbb{R}^N) := \{ u \in L^p(\mathbb{R}^N) : \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|u(x) - u(y)|^p}{|x - y|^{N+sp}} dx dy < \infty \}, $$
endowed with the norm

$$ \|u\|_{W^{s,p}(\mathbb{R}^N)} := \left( \|u\|_{L^p(\mathbb{R}^N)}^p + \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|u(x) - u(y)|^p}{|x - y|^{N+sp}} dx dy \right)^{1/p} $$

• The space $W_0^{s,p}(\Omega)$ is the set of functions

$$ W_0^{s,p}(\Omega) := \{u \in W^{s,p}(\mathbb{R}^N) : u = 0 \text{ a.e. in } \mathbb{R}^N \setminus \Omega\}, $$

and the norm is given by the Gagliardo semi-norm

$$ \|u\|_{W_0^{s,p}(\Omega)} := \left( \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|u(x) - u(y)|^p}{|x - y|^{N+sp}} dx dy \right)^{1/p} $$

We recall that by the fractional Poincaré inequality (e.g., in [29, Theorem 6.5]; see also Theorem 1.3 below), $\| \cdot \|_{W^{s,p}(\mathbb{R}^N)}$ and $\| \cdot \|_{W_0^{s,p}(\Omega)}$ are equivalent norms on $W_0^{s,p}(\Omega)$. From the results in [21], [29], we have that $W_0^{s,p}(\Omega)$ is continuously embedded in $L^r(\Omega)$ when $1 \le r \le \frac{Np}{N-sp}$ and compactly for $1 \le r < \frac{Np}{N-sp}$.

• Let $\alpha \in (0, 1]$, we consider the space of Hölder continuous functions:

$$ C^{\alpha}(\bar{\Omega}) = \{u \in C(\bar{\Omega}), \sup_{x,y \in \bar{\Omega}, x \neq y} \frac{|u(x) - u(y)|}{|x - y|^{\alpha}} < \infty\}, $$

endowed with the norm

$$ \|u\|_{C^{\alpha}(\bar{\Omega})} = \|u\|_{L^{\infty}(\bar{\Omega})} + \sup_{x,y \in \bar{\Omega}, x \neq y} \frac{|u(x) - u(y)|}{|x - y|^{\alpha}}. $$

• Let $T > 0$, and consider a measurable function

$$ u:]0, T[ \to W_0^{s,p}(\Omega), $$

and we denote $u(t)(x) := u(t, x)$. Let $C([0, T], W_0^{s,p}(\Omega))$ the space of continuous functions in $[0, T]$ with vector values in $W_0^{s,p}(\Omega)$, endowed with the norm

$$ \|u\|_{C([0,T], W_0^{s,p}(\Omega))} := \sup_{t \in [0,T]} \|u(t)\|_{W_0^{s,p}(\Omega)}. $$

• We denote by $d(\cdot)$ the distance function up to the boundary $\partial\Omega$. That means

$$ d(x) := \operatorname{dist}(x, \partial\Omega) = \inf_{y \in \partial\Omega} |x - y|. $$

• We define for $r > 0$, the sets

$$ \begin{gathered} \mathcal{M}_{ds}^{r}(\Omega) := \{u : \Omega \to \mathbb{R}^{+} : u \in L^{\infty}(\Omega) \text{ and } \exists c > 0 \text{ s.t.} \\ c^{-1}ds^{(x)} \le u^{r}(x) \le cd^{s}(x)\}, \quad (1.4) \\ \dot{V}_{+}^{r} := \{u : \Omega \to (0, \infty) : u^{1/r} \in W_{0}^{s,p}(\Omega)\}. \end{gathered} $$

• We define the weighted space

$$ L_{ds}^{\infty}(\Omega) := \{u : \Omega \to \mathbb{R} : u \in L^{\infty}(\Omega) \text{ s.t. } \frac{u}{ds(\cdot)} \in L^{\infty}(\Omega)\}. $$

Let $\phi_{1,s,p}$ be the positive normalized eigenfunction ($\|\phi_{1,s,p}\|_{L^\infty(\Omega)} = 1$) of $(-\Delta)_p^s$ in $W_0^{s,p}(\Omega)$ associated to the first eigenvalue $\lambda_{1,s,p}$. We recall that $\phi_{1,s,p} \in C^{0,\alpha}(\bar{\Omega})$ for some $\alpha \in (0, s]$ (see Theorem 1.1 in [27]) and $\phi_{1,s,p} \in \mathcal{M}_{ds}^1(\Omega)$ (see [27, Theorem 4.4] and [20, Theorem 1.5]). Next, we recall some results that will be used in the sequel.
**Proposition 1.1** (Discrete hidden convexity [9, Proposition 4.1]). Let $1 < p < \infty$ and $1 < q \le p$. For every $u_0, u_1 \ge 0$, we define

$$ \sigma_t(x) = [(1-t)u_0^q(x) + tu_1^q(x)]^{1/q}, \quad t \in [0,1], x \in \mathbb{R}^N. $$

Then

$$ |\sigma_t(x) - \sigma_t(y)|^p \le (1-t)|u_0(x) - u_0(y)|^p + t|u_1(x) - u_1(y)|^p, \quad t \in [0,1], x, y \in \mathbb{R}^N. $$

**Proposition 1.2** (Discrete Picone inequality [9, Proposition 4.2]). Let $1 < p < \infty$ and $1 < r \le p$. Let $u, v$ be two Lebesgue-measurable functions with $v \ge 0$ and $u > 0$. Then

$$ \begin{aligned} & |u(x) - u(y)|^{p-2}(u(x) - u(y)) \left[ \frac{v(x)^r}{u(x)^{r-1}} - \frac{v(y)^r}{u(y)^{r-1}} \right] \\ & \le |v(x) - v(y)|^r |u(x) - u(y)|^{p-r}. \end{aligned} $$

As we will see, Proposition 1.2 provides a comparison principle, barrier estimates and uniqueness of weak solutions.

**Theorem 1.3** ([21, Theorem 6.5]). Let $s \in (0,1)$, $p \ge 1$ with $N > sp$. Then, there exists a positive constant $C = C(N, p, s)$ such that, for any measurable and compactly supported $u : \mathbb{R}^N \to \mathbb{R}$ function, we have

$$ \|u\|_{L^{p_s^*}(\mathbb{R}^N)}^p \le C \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|u(x) - u(y)|^p}{|x - y|^{N+sp}} dx dy, $$

where $p_s^* = \frac{Np}{N-sp}$. Consequently, the space $W^{s,p}(\mathbb{R}^N)$ is continuously embedded in $L^q(\mathbb{R}^N)$ for $q \in [p, p_s^*]$.

**Theorem 1.4** (Aubin-Lions-Simon, [7, Theorem II.5.16]). Let $B_0 \subset B_1 \subset B_2$ be three Banach spaces. We assume that the embedding of $B_1$ in $B_2$ is continuous and that the embedding of $B_0$ in $B_1$ is compact. Let $p, r$ such that $1 \le p, r \le \infty$. For $T > 0$, we define

$$ E_{p,r} = \{v \in L^p([0, T]; B_0) : \frac{dv}{dt} \in L^r([0, T]; B_2)\}. $$

Then the following holds:

(a) If $p < \infty$, then the embedding of $E_{p,r}$ in $L^p([0, T]; B_1)$ is compact.

(b) If $p = \infty$ and $r > 1$, then the embedding of $E_{p,r}$ in $C([0, T]; B_1)$ is compact.

We now recall the definition of the strict ray-convexity.

**Definition 1.5.** Let $X$ be a real vector space. Let $C$ be a non empty convex cone in $X$. A functional $W: C \to \mathbb{R}$ will be called *ray-strictly convex* (strictly convex, respectively) if it satisfies

$$ W((1-t)v_1 + tv_2) \le (1-t)W(v_1) + tW(v_2), $$

for all $v_1, v_2 \in C$ and for all $t \in (0,1)$, where the inequality is always strict unless $\frac{v_1}{v_2} \equiv c > 0$ (always strict unless $v_1 \equiv v_2$, respectively).

**Remark 1.6.** We observe that by Proposition 1.1, the set $\dot{V}_+^r$ defined in (1.4) is a convex cone, i.e. for $\lambda \in (0, \infty)$, $f, g \in \dot{V}_+^r$ implies $\lambda f + g \in \dot{V}_+^r$.
**Proposition 1.7** (Convexity). Let $1 < p < \infty$ and $1 < r \le p$. The functional $W: \dot{V}_+^r \to \mathbb{R}_+$ defined by

$$W(w) := \frac{1}{p} \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|w(x)^{1/r} - w(y)^{1/r}|^p}{|x - y|^{N+sp}} dx dy,$$

is ray-strictly convex on $\dot{V}_+^r$. Furthermore, if $p \ne r$, then $W$ is even strictly convex on $\dot{V}_+^r$.

*Proof.* According to Definition 1.5, let us consider any $w_1, w_2 \in \dot{V}_+^r$ and $t \in [0, 1]$. Let us denote $w = tw_1 + (1-t)w_2$, we obtain by Proposition 1.1

$$W(w) \le t W(w_1) + (1-t) W(w_2). \quad (1.5)$$

If the equality holds, then

$$|w(x)^{1/r} - w(y)^{1/r}|^p = t |w_1(x)^{1/r} - w_1(y)^{1/r}|^p + (1-t) |w_2(x)^{1/r} - w_2(y)^{1/r}|^p$$

a.e. $x, y \in \mathbb{R}^N$. If $p = r$, we obtain

$$\|a\|_{\ell^r} - \|b\|_{\ell^r}^r = \|a - b\|_{\ell^r}^r \quad \text{a.e. } x, y \in \mathbb{R}^N,$$

where $\| \cdot \|_{\ell^r}$ denotes the $\ell^r$-norm in $\mathbb{R}^2$, and

$$a = ((tw_1(x))^1/r, ((1-t)w_2(x))^1/r), \quad b = ((tw_1(y))^1/r, ((1-t)w_2(y))^1/r).$$

Since $r > 1$, there exists a constant $c > 0$ such that $w_1 = cw_2$ a.e. $x \in \mathbb{R}^N$. Then, $W$ is ray-strictly convex on $\dot{V}_+^r$. On the other hand, if $p \ne r$ thanks to the strict convexity of $\tau \mapsto \tau^p$ on $\mathbb{R}^+$, we obtain $w_1 = w_2$ a.e. $x \in \mathbb{R}^N$ and $W$ is strictly convex on $\dot{V}_+^r$. $\square$

**Lemma 1.8.** Let $1 < p < \infty$. Then, for $1 < r \le p$ and for any $u, v$ two measurable and positive functions in $\Omega$:

$$|u(x) - u(y)|^{p-2} (u(x) - u(y)) \left[ \frac{u(x)^r - v(x)^r}{u(x)^{r-1}} - \frac{u(y)^r - v(y)^r}{u(y)^{r-1}} \right] \\ + |v(x) - v(y)|^{p-2} (v(x) - v(y)) \left[ \frac{v(x)^r - u(x)^r}{v(x)^{r-1}} - \frac{v(y)^r - u(y)^r}{v(y)^{r-1}} \right] \ge 0 \quad (1.6)$$

for a.e. $x, y \in \Omega$. Moreover, if $u, v \in W_0^{s,p}(\Omega)$ and if the equality occurs in (1.6) for a.e. $x, y \in \Omega$, then we have the following two statements:

(1) $u/v \equiv \text{const} > 0$ a.e. in $\Omega$.

(2) If also $p \ne r$, then $u \equiv v$ a.e. in $\Omega$.

*Proof.* Let $u, v$ be two measurable functions such that $u, v > 0$ in $\Omega$ and $1 < r \le p$. Then by using Proposition 1.2, we obtain for $x, y \in \Omega$,

$$|u(x) - u(y)|^{p-2} (u(x) - u(y)) \left[ \frac{v(x)^r}{u(x)^{r-1}} - \frac{v(y)^r}{u(y)^{r-1}} \right] \\ \le |v(x) - v(y)|^r |u(x) - u(y)|^{p-r}. \quad (1.7)$$

Let us start with the case $r = p$. By using the above inequality, in this case, we obtain

$$|u(x) - u(y)|^{p-2} (u(x) - u(y)) \left[ \frac{u(x)^p - v(x)^p}{u(x)^{p-1}} - \frac{u(y)^p - v(y)^p}{u(y)^{p-1}} \right] \\ \ge |u(x) - u(y)|^p - |v(x) - v(y)|^p. \quad (1.8)$$
By exchanging the roles of $u$ and $v$, we obtain

$$ \begin{aligned} & |v(x) - v(y)|^{p-2}(v(x) - v(y)) \left[ \frac{v(x)^p - u(x)^p}{v(x)^{p-1}} - \frac{v(y)^p - u(y)^p}{v(y)^{p-1}} \right] \\ & \geq |v(x) - v(y)|^p - |u(x) - u(y)|^p. \end{aligned} \quad (1.9) $$

Combining (1.8) and (1.9), we obtain

$$ \begin{aligned} & |u(x) - u(y)|^{p-2}(u(x) - u(y)) \left[ \frac{u(x)^p - v(x)^p}{u(x)^{p-1}} - \frac{u(y)^p - v(y)^p}{u(y)^{p-1}} \right] \\ & + |v(x) - v(y)|^{p-2}(v(x) - v(y)) \left[ \frac{v(x)^p - u(x)^p}{v(x)^{p-1}} - \frac{v(y)^p - u(y)^p}{v(y)^{p-1}} \right] \geq 0 \end{aligned} $$

which concludes the proof of (1.6) for $r = p$.

We deal finally with the case $1 < r < p$. By using Young's inequality, (1.7) implies

$$ \begin{aligned} & |u(x) - u(y)|^{p-2}(u(x) - u(y)) \left[ \frac{u(x)^r - v(x)^r}{u(x)^{r-1}} - \frac{u(y)^r - v(y)^r}{u(y)^{r-1}} \right] \\ & \geq \frac{r}{p} \left[ |u(x) - u(y)|^p - |v(x) - v(y)|^p \right]. \end{aligned} \quad (1.10) $$

Reversing the role of $u$ and $v$:

$$ \begin{aligned} & |v(x) - v(y)|^{p-2}(v(x) - v(y)) \left[ \frac{v(x)^r - u(x)^r}{v(x)^{r-1}} - \frac{v(y)^r - u(y)^r}{v(y)^{r-1}} \right] \\ & \geq \frac{r}{p} \left[ |v(x) - v(y)|^p - |u(x) - u(y)|^p \right]. \end{aligned} \quad (1.11) $$

Adding the above inequalities, we obtain (1.6).

Now, let us consider $u, v \in W_0^{s,p}(\Omega)$, such that $u > 0, v > 0$ a.e. in $\Omega$ and $\theta \in (0,1)$. Setting $w := (1-\theta)u^r + \theta v^r$, one can easily check that $w \in \dot{V}_+^r$. Thus, by Proposition 1.7, it is easy to prove that the function, defined in $[0,1]$,

$$ \theta \mapsto \Phi(\theta) := W(w) = W((1-\theta)u^r + \theta v^r) $$

is convex, differentiable and for $\theta \in (0,1)$:

$$ \Phi'(\theta) = \int_{\mathbb{R}^{2N} \setminus (\Omega^c \times \Omega^c)} \frac{|w(x)^{1/r} - w(y)^{1/r}|^{p-2} (w(x)^{1/r} - w(y)^{1/r})}{|x-y|^{N+sp}} \times \left( \frac{v(x)^r - u(x)^r}{w(x)^{1-1/r}} - \frac{v(y)^r - u(y)^r}{w(y)^{1-1/r}} \right) dx dy. $$

Finally, let us assume that the equality in (1.6) holds. By the monotonicity of $\Phi' : (0,1) \to \mathbb{R}$, we deduce that $\Phi'(\theta) = \text{const}$ in $(0,1)$. It follows that $\Phi : [0,1] \to \mathbb{R}$ must be linear, i.e.

$$ \Phi(\theta) = W(w) = (1 - \theta)\Phi(0) + \theta\Phi(1) = (1 - \theta)W(u^r) + \theta W(v^r), $$

for all $\theta \in [0,1]$. We conclude that $u \equiv \text{const.}v$ with $\text{const} > 0$ and if $p \neq r$, then $u \equiv v$, thanks to Proposition 1.7. $\square$
1.3. Main results. We consider the associated problem of (1.1),

$$
\begin{align}
v^{q-1} \partial_t (v^q) + (-\Delta)_p^s v &= h(t, x)v^{q-1} + f(x, v) \quad \text{in } Q_T, \nonumber \\
v &> 0 \quad \text{in } Q_T, \nonumber \\
v &= 0 \quad \text{on } \Gamma_T, \tag{1.12} \\
v(0, \cdot) &= v_0 \quad \text{in } \Omega. \nonumber
\end{align}
$$

**Claim 1.9.** *Any bounded weak solution of the above problem is also a weak solution to problem (1.1).*

To this aim, we introduce the notion of the weak solution to problem (1.12) as
follows.

**Definition 1.10.** Let $T > 0$. A weak solution to problem (1.12) is any nonnegative function $v \in L^\infty(0, T; W_0^{s,p}(\Omega)) \cap L^\infty(Q_T)$ such that $v > 0$ in $\Omega$, $\partial_t(v^q) \in L^2(Q_T)$ and satisfying for any $t \in (0, T]$:

$$
\begin{align*}
& \int_0^t \int_{\Omega} \partial_t (v^q) v^{q-1} \varphi \, dx \, ds \\
& + \int_0^t \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|v(s, x) - v(s, y)|^{p-2} (v(s, x) - v(s, y)) (\varphi(s, x) - \varphi(s, y))}{|x - y|^{N+sp}} \, dx \, dy \, ds \\
& = \int_0^t \int_{\Omega} (h(s, x) v^{q-1} + f(x, v)) \, \varphi \, dx \, ds,
\end{align*}
$$

for any $\varphi \in L^2(Q_T) \cap L^1(0, T; W_0^{s,p}(\Omega))$, with $v(0, \cdot) = v_0$ a.e. in $\Omega$.

**Remark 1.11.** According to Definition 1.10, a weak solution of (1.12) belongs to $L^\infty(Q_T)$. Then, we obtain

$$
\frac{q}{2q-1} \partial_t (v^{2q-1}) = v^{q-1} \partial_t (v^q)
$$

weakly, and we deduce that a weak solution to (1.12) is a weak solution to (1.1).

Our main result about existence and properties of solutions to (1.12) is as follows.

**Theorem 1.12.** Let $T > 0$ and $q \in (1, p]$. Assume that $f$ satisfies (H1)–(H3), (H6) and

(H7) The map $x \mapsto \phi_{1,s,p}^{1-q}(x)f(x, \phi_{1,s,p}(x))$ belongs to $L^2(\Omega)$.

Assume in addition that $h \in L^\infty(Q_T)$ satisfies (H4), (H5) and that $v_0 \in \mathcal{M}_{ds}^1(\Omega) \cap W_0^{s,p}(\Omega)$. Then there exists a unique weak solution $v$ to (1.12). Furthermore,

(i) $v \in C([0, T]; W_0^{s,p}(\Omega))$ and satisfies for any $t \in [0, T]$ the energy estimate

$$
\begin{align*}
& \int_0^t \int_{\Omega} \left(\frac{\partial v^q}{\partial t}\right)^2 dx \, ds + \frac{q}{p} \|v(t)\|_{W_0^{s,p}(\Omega)}^p \\
& = \int_0^t \int_{\Omega} h\left(\frac{\partial v^q}{\partial t}\right) dx \, ds + \int_0^t \int_{\Omega} \frac{f(x, v)}{v^{q-1}} \frac{\partial v^q}{\partial t} dx \, ds + \frac{q}{p} \|v_0\|_{W_0^{s,p}(\Omega)}^p.
\end{align*}
$$

(ii) If $w$ is a weak solution to (1.12) associated to the initial data $w_0 \in \mathcal{M}_{ds}^1(\Omega) \cap W_0^{s,p}(\Omega)$ and the right hand side $g \in L^\infty(Q_T)$ satisfying (H4) and (H5), then the following estimate (T-accretivity in $L^2(\Omega)$) holds:

$$
\| (v^q(t) - w^q(t))^+ \|_{L^2(\Omega)} \le \|(v_0^q - w_0^q)^+\|_{L^2(\Omega)} + \int_0^t \|(h(s) - g(s))^+\|_{L^2(\Omega)} ds \quad (1.13)
$$

for any $t \in [0, T]$.
The T-acretivity in L² stated in (1.13) was proved for p-Laplace operators in [22] with a different approach (by the study of properties of the associated subdifferential via the potential theory) and for quasilinear elliptic operators with variable exponents in [2] (see also [6] and [3] for related issues). The uniqueness of the solution in Theorem 1.12 can be also obtained by the following theorem under less restrictive assumptions about $v_0$ and $h$.

**Theorem 1.13.** Let $v, w$ be two solutions of the problem (1.12) in sense of Definition 1.10, with respect to the initial data $v_0, w_0 \in L^{2^q}(\Omega)$, $v_0, w_0 \ge 0$ and $h, \tilde{h} \in L^2(Q_T)$. Then, for any $t \in [0, T]$,

$$ \|v^q(t) - w^q(t)\|_{L^2(\Omega)} \le \|v_0^q - w_0^q\|_{L^2(\Omega)} + \int_0^t \|h(s) - \tilde{h}(s)\|_{L^2(\Omega)} ds. \quad (1.14) $$

Using the theory of maximal accretive operators, we introduce the nonlinear operator $\mathcal{T}_q: L^2(\Omega) \supset D(\mathcal{T}_q) \to L^2(\Omega)$ defined by

$$ \mathcal{T}_q u = u^{\frac{1-q}{q}} \left( 2\mathrm{P.V.} \int_{\mathbb{R}^N} \frac{|u^{1/q}(x) - u^{1/q}(y)|^{p-2} (u^{1/q}(x) - u^{1/q}(y))}{|x-y|^{N+sp}} dy - f(x, u^{1/q}) \right) \quad (1.15) $$

with

$$ D(\mathcal{T}_q) = \{ w : \Omega \to \mathbb{R}^+, \ w^{1/q} \in W_0^{s,p}(\Omega), \ w \in L^2(\Omega), \mathcal{T}_q w \in L^2(\Omega) \}. $$

Using the T-accretive property of $\mathcal{T}_q$ in $L^2(\Omega)$ proved below and under additional assumptions on regularity of initial data, we obtain the following stabilization result for the weak solutions to the problem (1.12).

**Theorem 1.14.** Assume that the hypothesis in Theorem 1.12 hold for any $T > 0$. Let $v$ be the weak solution of the problem (1.12) with the initial data $v_0 \in \mathcal{M}_d^1(\Omega) \cap W_0^{s,p}(\Omega)$. Assume in addition that there exists $h_\infty \in L^\infty(\Omega)$ such that

$$ l(t) \|h(t, \cdot) - h_\infty\|_{L^2(\Omega)} = O(1) \quad \text{as } t \to \infty \quad (1.16) $$

with $l$ continuous and positive on $]s_0; +\infty[$ and $\int_s^{+\infty} \frac{dt}{l(t)} < +\infty$, for some $s > s_0 \ge 0$. Then, for any $r \ge 1$,

$$ \|v^q(t, \cdot) - v_\infty^q\|_{L^r(\Omega)} \to 0 \quad \text{as } t \to \infty, $$

where $v_\infty$ is the unique stationary solution to (1.12) associated to the potential $h_\infty$.

This article is organized as follows: In Section 2, we study the stationary nonlinear problem

$$ v^{2q-1} + \lambda(-\Delta)_p^s v = h_0(x)v^{q-1} + \lambda f(x, v) \quad \text{in } \Omega, $$

$$ v > 0 \quad \text{in } \Omega, $$

$$ v = 0 \quad \text{in } \mathbb{R}^N \setminus \Omega, $$

related to the parabolic problem (1.12) and establish the existence and the uniqueness results in case $h_0 \in L^\infty(\Omega)$ [Theorem 2.2, Corollary 2.4] and in case $h_0 \in L^2(\Omega)$ [Theorem 2.5, Corollary 2.6]. Section 3 is devoted to prove Theorem 1.12. The proof is divided into three main steps. First, by using a semi-discretization in time with
implicit Euler method, we prove the existence of a weak solution in sense of Definition 1.10 (see Theorem 3.1). Next, we prove the contraction property given in Theorem 1.13 which implies the uniqueness of the weak solution stated in Corollary 3.2. The regularity of weak solutions is established in Theorem 3.4 that brings the completion of the proof of Theorem 1.12. In Section 4, we show the stabilization result (see Theorem 1.14) for problem (1.12) via classical arguments of the semi-group theory. Finally in the appendix 5.1, we establish some new regularity results ($L^\infty$ bound) for a class of quasilinear elliptic equations involving fractional $p$-Laplace operator. Via the Picone identity, we also obtain a new weak comparison principle that provides existence of barrier functions for stationary problems of (1.12).

## 2. $p$-FRACTIONAL ELLIPTIC EQUATION ASSOCIATED WITH PROBLEM (1.1)

The aim of this section is to study the elliptic problem corresponding to (1.12). For this, we have several cases.

### 2.1. Potential $h_0 \in L^\infty(\Omega)$.

$$ 
\begin{array}{l}
v^{2q-1} + \lambda(-\Delta)_p^s v = h_0(x)v^{q-1} + \lambda f(x, v) \quad \text{in } \Omega, \\
v > 0 \quad \text{in } \Omega \\
v = 0 \quad \text{in } \mathbb{R}^N \setminus \Omega,
\end{array} 
\quad (2.1) $$

where $\lambda$ is a positive parameter and $h_0 \in (L^\infty(\Omega))^+$ satisfying the hypothesis (H8) $h_0(x) \ge \lambda \underline{h}(x)$ for a.e. in $\Omega$, where $\underline{h}$ is defined in (H4).

We have the following notion of weak solutions.

**Definition 2.1.** A weak solution of the problem (2.1) is any nonnegative and nontrivial function $v \in \mathbf{W} := W_0^{s,p}(\Omega) \cap L^{2q}(\Omega)$ such that for any $\varphi \in \mathbf{W}$,

$$ 
\begin{aligned}
& \int_{\Omega} v^{2q-1} \varphi dx + \lambda \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|v(x) - v(y)|^{p-2} (v(x) - v(y)) (\varphi(x) - \varphi(y))}{|x-y|^{N+sp}} dx dy \\
& = \int_{\Omega} h_0 v^{q-1} \varphi dx + \lambda \int_{\Omega} f(x, v) \varphi dx.
\end{aligned} 
\quad (2.2) $$

We first investigate the existence and uniqueness of the weak solution to (2.1).

**Theorem 2.2.** Assume that $f$ satisfies (H1), (H2), (H6). In addition suppose that $h_0 \in L^\infty(\Omega)$ and satisfies (H8). Then, for any $1 < q \le p$ and $\lambda > 0$, there exists a positive weak solution $v \in C(\bar{\Omega}) \cap M_{ds}^1(\Omega)$ to (2.1).

Moreover, let $v_1, v_2$ be two weak solutions to (2.1) with $h_1, h_2 \in L^\infty(\Omega)$ satisfy (H8), respectively, we have (with the notation $t^+ = \max\{0, t\}$),

$$ \|(v_1^q - v_2^q)^+\|_{L^2} \le \|(h_1 - h_2)^+\|_{L^2}. \quad (2.3) $$

*Proof.* We divided the proof into 3 steps.

**Step 1: Existence of a weak solution.** Consider the energy functional $\mathcal{J}$ corresponding to the problem (2.1), defined on $\mathbf{W}$ equipped with the Cartesian norm $\| \cdot \|_{\mathbf{W}} = \| \cdot \|_{W_0^{s,p}(\Omega)} + \| \cdot \|_{L^{2q}(\Omega)}$ by

$$ 
\begin{aligned}
\mathcal{J}(v) = & \frac{1}{2q} \int_{\Omega} v^{2q} dx + \frac{\lambda}{p} \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|v(x) - v(y)|^p}{|x - y|^{N+ps}} dx dy \\
& - \frac{1}{q} \int_{\Omega} h_0(v^+)^q dx - \lambda \int_{\Omega} F(x, v) dx
\end{aligned} 
\quad (2.4) $$
where

$$F(x,t) = \begin{cases} \int_{0}^{t} f(x,s)ds & \text{if } 0 \le t < +\infty, \\ 0 & \text{if } -\infty < t < 0. \end{cases}$$

We extend accordingly the domain of f to all of $\Omega \times \mathbb{R}$ by setting

$$f(x,t) = \frac{\partial F}{\partial t}(x,t) = 0 \quad \text{for } (x,t) \in \Omega \times (-\infty, 0).$$

From (H1) and (H2) there exists $C > 0$ large enough such that for any $(x, s) \in \Omega \times \mathbb{R}^+$,

$$0 \le f(x, s) \le C(1 + s^{q-1}). \quad (2.5)$$

Thus, we infer that:

* $\mathcal{J}$ is well defined and weakly lower semi-continuous on $\mathbf{W}$.

* From (2.5), the Hölder inequality and Theorem 1.3, we obtain

$$\begin{aligned} \mathcal{J}(v) &\ge \frac{1}{2q} \|v\|_{L^{2q}(\Omega)}^{2q} + \frac{\lambda}{p} \|v\|_{W_0^{s,p}(\Omega)}^{p} - \frac{1}{q} \|h_0\|_{L^2(\Omega)} \|v\|_{L^{2q}(\Omega)}^{q} - C\lambda \int_{\Omega} |v|dx \\ &\quad - \lambda \frac{C}{q} \int_{\Omega} |v|^q dx \\ &\ge \|v\|_{L^{2q}(\Omega)}^{q} (c_1 \|v\|_{L^{2q}(\Omega)}^{q} - c_2) + \|v\|_{W_0^{s,p}(\Omega)} (c_3 \|v\|_{W_0^{s,p}(\Omega)}^{p-1} - c_4), \end{aligned}$$

where the constants $c_1, c_2, c_3$ and $c_4$ do not depend on $v$. Therefore, we obtain that $\mathcal{J}(v)$ is coercive on $\mathbf{W}$. Therefore, $\mathcal{J}$ admits a global minimizer on $\mathbf{W}$, denoted by $v_0$. Thus, adopting the notation $t = t^+ - t^-$, we have

$$\begin{aligned} \mathcal{J}(v_0) = \mathcal{J}(v_0^+) &+ \frac{1}{2q} \int_{\Omega} (v^-)^{2q} dx + \frac{\lambda}{p} \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|(v^-)(x) - (v^-)(y)|^p}{|x - y|^{N+ps}} dx dy \\ &+ \frac{2\lambda}{p} \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|(v^-)(x) - (v^+)(y)|^p}{|x - y|^{N+ps}} dx dy \ge \mathcal{J}(v_0^+). \end{aligned}$$

Therefore, $v_0 \ge 0$. In order to show that $v_0 \ne 0$ in $\Omega$, we find a suitable function $v$ in $\mathbf{W}$ such that $\mathcal{J}(v) < 0 = \mathcal{J}(0)$. For that, we start by dealing with the case $q < p$. Let $\phi \in C_c^1(\Omega)$ be nonnegative and non trivial with $\text{supp}(\phi) \subset \text{supp}(\underline{h})$. Then, for any $t > 0$,

$$\mathcal{J}(t\phi) \le c_1 t^{2q} + c_2 t^p - c_3 t^q,$$

where the constants $c_1, c_2$ and $c_3$ are independent of $t$ and $c_3 > 0$ thanks to $h_0 \ge \lambda \underline{h} \ne 0$. Hence for $t > 0$ small enough, $\mathcal{J}(t\phi) < 0$. We now consider the remaining case $q = p$. Assumption (H6) implies that for $c > 0$ small enough there exists $s_0 = s_0(c) > 0$ such that

$$\lambda \underline{h}(x) s^{p-1} + \lambda f(x, s) > \lambda (\lambda_{1,p,s} + c) s^{p-1}$$

for all $s \le s_0$ and uniformly in $x \in \Omega$. Hence, for $\epsilon$ small enough, we deduce that

$$\begin{aligned} \mathcal{J}(\epsilon\phi_{1,p,s}) &< \frac{1}{2p} \|\phi_{1,p,s}\|_{L^{2p}(\Omega)}^{2p} \epsilon^{2p} + \frac{\lambda}{p} \|\phi_{1,p,s}\|_{W_0^{s,p}(\Omega)}^{p} \epsilon^p \\ &\quad - \frac{\lambda}{p} (\lambda_{1,p,s} + c) \|\phi_{1,p,s}\|_{L^p(\Omega)}^{p} \epsilon^p \\ &\quad = \epsilon^p \left( \frac{1}{2p} \|\phi_{1,p,s}\|_{L^{2p}(\Omega)}^{2p} \epsilon^p - \frac{c\lambda}{p} \|\phi_{1,p,s}\|_{L^p(\Omega)}^{p} \epsilon^p \right) < 0. \end{aligned}$$

Since $\mathcal{J}(0) = 0$, we deduce $v_0 \ne 0$. From the Gâteaux differentiability of $\mathcal{J}$, we obtain that $v_0$ satisfies (2.2).
**Step 2: Regularity and positivity of weak solutions.** We first claim that all weak solutions to (2.1) belongs to $L^\infty(\Omega)$. To this aim, we adapt arguments from [[23], Theorem 3.2]. Precisely, let $v_0$ be a weak solution. Then, it is enough to prove that

$$ \|v_0\|_{L^\infty(\Omega)} \le 1 \quad \text{if} \quad \|v_0\|_{L^p(\Omega)} \le \delta \quad \text{for some } \delta > 0 \text{ small enough.} \quad (2.6) $$

For this purpose, we consider the function $w_k$ defined as follows

$$ w_k(x) := (v_0(x) - (1 - 2^{-k}))^+ \quad \text{for } k \ge 1. $$

We first state the following straightforward observations about $w_k(x)$,

$$ w_k \in W_0^{s,p}(\Omega) \quad \text{and} \quad w_k = 0 \quad \text{a.e. in} \quad \partial\Omega, $$

and

$$ \begin{gathered} w_{k+1}(x) \le w_k(x) \quad \text{a.e. in } \mathbb{R}^N, \\ v_0(x) < (2^{k+1} + 1)w_k(x) \quad \text{for } x \in \{w_{k+1} > 0\}. \end{gathered} \quad (2.7) $$

Also the inclusion

$$ \{w_{k+1} > 0\} \subseteq \{w_k > 2^{-\frac{p}{k+1}}\} \quad (2.8) $$

holds for all $k \in \mathbb{N}$.

Setting $V_k := \|w_k\|_{L^p(\Omega)}^p$, using (2.5), (2.7) and the inequality

$$ x^+ - y^+|^p \le |x - y|^{p-2}(x^+ - y^+)(x - y) $$

for any $x, y \in \mathbb{R}$, we obtain

$$ \begin{aligned} \lambda \|w_{k+1}\|_{W_0^{s,p}(\Omega)}^p &= \lambda \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|w_{k+1}(x) - w_{k+1}(y)|^p}{|x - y|^{N+sp}} dx dy \\ &\le \lambda \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|v_0(x) - v_0(y)|^{p-2}(w_{k+1}(x) - w_{k+1}(y))(v_0(x) - v_0(y))}{|x - y|^{N+sp}} dx dy \\ &\le \int_{\Omega} (h_0(x)v_0^{q-1} + \lambda f(x, v_0))w_{k+1} dx \\ &\le C_1 [ \int_{\{w_{k+1}>0\}} w_{k+1} dx + \int_{\{w_{k+1}>0\}} v_0^{q-1} w_{k+1} dx ] \\ &\le C_1 [|w_{k+1}>0|^{1-\frac{1}{p}} V_k^{1/p} + (2^{k+1} + 1)^{q-1} |\{w_{k+1}>0\}|^{1-\frac{q}{p}} V_k^{\frac{q}{p}}] \end{aligned} $$

where $C_1 > 0$ is a constant. Now, from (2.8) we have

$$ V_k = \int_{\Omega} w_k^p dx \ge \int_{\{w_{k+1}>0\}} w_k^p dx \ge 2^{-(k+1)p} |\{w_{k+1}>0\}|. \quad (2.9) $$

Therefore,

$$ \|w_{k+1}\|_{W_0^{s,p}(\Omega)}^p \le C_2(2^{k+1} + 1)^{p-1}V_k $$

where $C_2 > 0$ is a constant. On the other hand, by the Hölder's inequality, fractional Sobolev imbeddings (Theorem 1.3) and (2.9), we obtain

$$ V_{k+1} = \int_{\{w_{k+1}>0\}} w_{k+1}^p dx \le C_3 \|w_{k+1}\|_{W_0^{s,p}(\Omega)}^p (2^{(k+1)p} V_k)^{\frac{sp}{N}}, $$

where $C_3 > 0$ is a constant. Hence, the above inequality

$$ V_{k+1} \le C^k V_k^{1+\alpha}, \quad \text{for all } k \in \mathbb{N} $$
holds for a suitable constant $C > 1$ and $\alpha = \frac{sp}{N}$. This implies that

$$ \lim_{k \to \infty} V_k = 0 \qquad (2.10) $$

provided that

$$ \|v_0\|_{L^p(\Omega)}^p = V_0 \le C^{-\frac{1}{\alpha^2}} =: \delta^p $$

as it can be easily checked. Since $w_k$ converges to $(v_0 - 1)^+$ a.e. in $\mathbb{R}^N$, from (2.10) we infer that (2.6) holds as desired. Then, we deduce that $v_0 \in L^\infty(\Omega)$ and [27, Theorem 1.1] provides the $C^{0,\alpha}(\bar{\Omega})$-regularity of $v_0$, for some $\alpha \in (0, s]$. Now, we show that $v_0 > 0$ in $\Omega$. We argue by contradiction: Suppose that there exists $x_0 \in \Omega$, where $v_0(x_0) = 0$, then it follows that

$$ \begin{aligned} & 0 > 2\lambda \int_{\mathbb{R}^N} \frac{|v_0(x_0) - v_0(y)|^{p-2}(v_0(x_0) - v_0(y))}{|x_0 - y|^{N+sp}} dy \\ & = h_0(x)v_0(x_0)^{q-1} + \lambda f(x_0, v_0(x_0)) - v_0(x_0)^{2q-1} = 0 \end{aligned} $$

from which we obtain a contradiction. Thus $v_0 > 0$ in $\Omega$. Finally, starting with the case $q = p$, the Hopf lemma (see [20, Theorem 1.5]) implies that $v_0 \ge k d^s(x)$ for some $k > 0$. Next, supposing $q < p$, we have that for $\epsilon > 0$ small enough, $\epsilon\phi_{1,s,p}$ is a subsolution to problem (2.1). Indeed, for a constant $\epsilon > 0$ small enough, we have

$$ (\epsilon\phi_{1,s,p})^{2q-1} + \lambda(-\Delta)_p^s (\epsilon\phi_{1,s,p}) \le h_0(x)(\epsilon\phi_{1,s,p})^{q-1} + \lambda f(x, \epsilon\phi_{1,s,p}) \quad \text{in } \Omega. $$

From the comparison principle (Theorem 5.4), we obtain $\epsilon\phi_{1,s,p} \le v_0$. Then, we deduce that $v_0 \ge kd^s(x)$ for some $k > 0$. Again by using [[27], Theorem 4.4], we obtain that $v_0 \in \mathcal{M}_{d^s}^1(\Omega)$.

**Step 3: Contraction property** (2.3) Let $v_1, v_2 \in \mathcal{M}_{d^s}^1(\Omega)$ be two weak solutions of (2.1) associated to $h_1$ and $h_2$ respectively. Namely, for any $\Phi, \Psi \in \mathbf{W}$ we have

$$ \begin{aligned} & \int_{\Omega} v_1^{2q-1} \Phi dx + \lambda \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|v_1(x) - v_1(y)|^{p-2} (v_1(x) - v_1(y)) (\Phi(x) - \Phi(y))}{|x-y|^{N+sp}} dx dy \\ & = \int_{\Omega} h_1 v_1^{q-1} \Phi dx + \lambda \int_{\Omega} f(x, v_1) \Phi dx \end{aligned} $$

and

$$ \begin{aligned} & \int_{\Omega} v_2^{2q-1} \Psi dx + \lambda \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|v_2(x) - v_2(y)|^{p-2} (v_2(x) - v_2(y)) (\Psi(x) - \Psi(y))}{|x-y|^{N+sp}} dx dy \\ & = \int_{\Omega} h_2 v_2^{q-1} \Psi dx + \lambda \int_{\Omega} f(x, v_2) \Psi dx. \end{aligned} $$

Since $v_1, v_2 \in \mathcal{M}_{d^s}^1(\Omega) \cap W_0^{1,s}(\Omega)$, we obtain that

$$ \Phi = \frac{(v_1^q - v_2^q)^+}{v_1^{q-1}}, \quad \Psi = \frac{(v_2^q - v_1^q)^-}{v_2^{q-1}} $$

are well-defined and belong to $\mathbf{W}$.

Subtracting the two expressions above and using (H2) and Lemma 1.8, we obtain

$$ \int_{\Omega} ((v_1^q - v_2^q)^+)^2 dx \le \int_{\Omega} (h_1 - h_2)(v_1^q - v_2^q)^+ dx. $$

Finally, applying the Hölder inequality we obtain (2.3). □

**Remark 2.3.** Inequality (2.3) implies the uniqueness of the weak solution to the problem (2.1) in the sense of Definition 2.2 in $\mathcal{M}_{d^s}^1(\Omega)$.
From Theorem 2.2, we deduce the T-accretivity of $\mathcal{T}_q$ (see (1.15)) as follows.

**Corollary 2.4.** Let $\lambda > 0$, $q \in (1, p]$, $f : \Omega \times \mathbb{R}^+ \to \mathbb{R}^+$ satisfies (H1), (H2), (H6). Assume in addition that $h_0 \in L^\infty(\Omega)$ satisfies (H8). Then, there exists a unique solution $u \in C(\bar{\Omega})$ of the problem

$$ 
\begin{aligned}
& u + \lambda \mathcal{T}_q u = h_0 \quad \text{in } \Omega, \\
& u > 0 \quad \text{in } \Omega, \\
& u \equiv 0 \quad \text{in } \mathbb{R}^N \setminus \Omega.
\end{aligned} 
\quad (2.11) $$

Namely, $u$ belongs to $\dot{V}_+^q \cap \mathcal{M}_{ds}^{1/q}(\Omega)$, and satisfies

$$ 
\begin{aligned}
& \int_{\Omega} u \Psi \, dx + \lambda \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} |u^{1/q}(x) - u^{1/q}(y)|^{p-2} (u^{1/q}(x) - u^{1/q}(y)) \\
& \times \left( (u^{\frac{1-q}{q}} \Psi)(x) - (u^{\frac{1-q}{q}} \Psi)(y) \right) / |x-y|^{N+sp} \, dx dy \\
& = \int_{\Omega} h_0 \Psi \, dx + \lambda \int_{\Omega} f(x, u^{1/q}) u^{\frac{1-q}{q}} \Psi \, dx
\end{aligned} 
\quad (2.12) $$

for any $\Psi$ such that

$$ |\Psi|^{1/q} \in L_{ds}^{\infty}(\Omega) \cap W_0^{s,p}(\Omega). \quad (2.13) $$

Moreover, if $u_1$ and $u_2$ are two solutions of (2.11), corresponding to $h_1$ and $h_2$ respectively, then

$$ \|(u_1 - u_2)^+\|_{L^2} \le \|(u_1 - u_2 + \lambda(\mathcal{T}_q(u_1) - \mathcal{T}_q(u_2)))^+\|_{L^2}. \quad (2.14) $$

**Proof.** We define the energy functional $\xi$ on $\dot{V}_+^q \cap L^2(\Omega)$ as $\xi(u) = \mathcal{J}(u^{1/q})$, where $\mathcal{J}$ is defined in (2.4). Let $v_0$ be the weak solution of (2.1) and the global minimizer of (2.4). We set $u_0 = v_0^q$. Then

$$ u_0 \in \dot{V}_+^q \cap \mathcal{M}_{ds}^{1/q}(\Omega). $$

Let $\Psi \ge 0$ satisfy (2.13), then there exists $t_0 = t_0(\Psi) > 0$ such that for $t \in (0, t_0)$, $u_0 + t\Psi > 0$. Hence, we have

$$ 
\begin{aligned}
& 0 \le \xi(u_0 + t\Psi) - \xi(u_0) \\
& = \frac{1}{2q} \left( \int_{\Omega} (t\Psi)^2 dx + 2t \int_{\Omega} u_0 \Psi dx \right) - \frac{1}{q} \int_{\Omega} t h_0 \Psi dx \\
& \quad + \frac{\lambda}{p} \left( \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|(u_0 + t\Psi)^{1/q}(x) - (u_0 + t\Psi)^{1/q}(y)|^p}{|x-y|^{N+ps}} dx dy \right. \\
& \quad \left. - \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|(u_0)^{1/q}(x) - (u_0)^{1/q}(y)|^p}{|x-y|^{N+ps}} dx dy \right) \\
& \quad - \lambda \left( \int_{\Omega} F(x, (u_0 + t\Psi)^{1/q}) dx - \int_{\Omega} F(x, (u_0)^{1/q}) dx \right).
\end{aligned} $$

Then dividing by $t$ and passing to the limit $t \to 0$, we obtain that $u_0$ satisfies (2.12). On the other hand, consider $u_1 \in \dot{V}_+^q \cap \mathcal{M}_{ds}^{1/q}(\Omega)$ a solution satisfying (2.12). Thus $v_1 = u_1^{1/q}$ satisfies (2.2), by Remark 2.3, we deduce $v_1 = v_2$. Finally, (2.14) follows from (2.3). $\square$
2.2. Potential $h_0 \in L^2(\Omega)$. In this subsection, we extend the existence results above.

**Theorem 2.5.** Assume that $f$ satisfies (H1), (H2), (H6). Then, for any $1 < q \le p$, $\lambda > 0$ and $h_0 \in L^2(\Omega)$ satisfies (H8), there exists a positive weak solution $v \in W$ to (2.1). Moreover assuming that $h_0$ belongs to $L^r(\Omega)$ for some $r > \frac{N}{sp}$, $v \in L^\infty(\Omega)$. Moreover, let $v_1, v_2$ be two weak solutions to (2.1) associated with $h_1, h_2 \in L^2(\Omega)$, respectively, satisfy (H8). Then, we have

$$ \|(v_1^q - v_2^q)^+\|_{L^2} \le \|(h_1 - h_2)^+\|_{L^2}. \quad (2.15) $$

*Proof.* Let $\tilde{h}_n \in C_c^1(\Omega)$, $\tilde{h}_n \ge 0$ with $\tilde{h}_n \to h_0$ in $L^2(\Omega)$, we take $h_n = \max(\tilde{h}_n, \lambda \underline{h})$. By Theorem 2.2, for any $n \ge n_0$, define $v_n \in C^{0,\alpha}(\bar{\Omega}) \cap \mathcal{M}_{ds}^1(\Omega)$ as the unique positive weak solution of (2.1). Then, for any $\varphi \in W$,

$$ \begin{aligned} & \int_{\Omega} v_n^{2q-1} \varphi dx \\ & + \lambda \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|v_n(x) - v_n(y)|^{p-2} (v_n(x) - v_n(y)) (\varphi(x) - \varphi(y))}{|x-y|^{N+sp}} dx dy \quad (2.16) \\ & = \int_{\Omega} h_n v_n^{q-1} \varphi dx + \lambda \int_{\Omega} f(x, v_n) \varphi dx. \end{aligned} $$

One has

$$ (a-b)^{2r} \le (a^r - b^r)^2 \quad \text{for any } r \ge 1, a, b \ge 0 \quad (2.17) $$

from which together with (2.3) it follows for any $n, m \in \mathbb{N}^*$,

$$ \|(v_n - v_m)^+\|_{L^{2q}} \le \|(v_n^q - v_m^q)^+\|_{L^2}^{1/q} \le \|(h_n - h_m)^+\|_{L^2}^{1/q}. $$

Thus we deduce that $(v_n)$ converges to some $v \in L^{2q}(\Omega)$. We infer that the limit $v$ does not depend on the choice of the sequence $(h_n)$. Indeed, consider $\tilde{h}_n \ne h_n$ such that $\tilde{h}_n \to h_0$ in $L^2(\Omega)$ and $\tilde{v}_n$ the positive solution to (2.1) corresponding to $\tilde{h}_n$ which converges to $\tilde{v}$. Then, for any $n \in \mathbb{N}$, (2.3) implies

$$ \|(v_n^q - \tilde{v}_n^q)^+\|_{L^2} \le \|(h_n - \tilde{h}_n)^+\|_{L^2} $$

and passing to the limit we obtain $\tilde{v} \ge v$ and then by reversing the role of $v$ and $\tilde{v}$, we obtain $\tilde{v} = v$.

For $n \in \mathbb{N}^*$, let $h_n = \min\{h_0, n\lambda \underline{h}\}$. So, it is easy to check by (2.3), $(v_n)_{n \in \mathbb{N}}$ is nondecreasing and for any $n \in \mathbb{N}^*$, $v_n \le v$ a.e. in $\Omega$ which implies

$$ v(x) \ge v_1(x) \ge c d^s(x) > 0 \quad \text{in } \Omega \quad (2.18) $$

for some $c$ independent of $n$. We choose $\varphi = v_n$ in (2.16), by the Hölder inequality and (2.5), we obtain

$$ \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|v_n(x) - v_n(y)|^p}{|x - y|^{N+sp}} dx dy \le C[\|v_n\|_{L^{2q}(\Omega)}^q (\|h_n\|_{L^2(\Omega)} + 1) + \|v_n\|_{L^{2q}(\Omega)}] \quad (2.19) $$

where $C$ does not depend on $n$. Then, we deduce that $(v_n)_{n \in \mathbb{N}}$ is uniformly bounded in $W_0^{s,p}(\Omega)$. Hence,

$$ \left\{ \frac{|v_n(x) - v_n(y)|^{p-2} (v_n(x) - v_n(y))}{|x - y|^{\frac{N+sp}{p'}}} \right\} \text{ is bounded in } L^{p'}(\mathbb{R}^N \times \mathbb{R}^N) $$
where $p' = \frac{p}{p-1}$ and by the pointwise convergence of $v_n$ to $v$, we obtain

$$ \frac{|v_n(x) - v_n(y)|^{p-2}(v_n(x) - v_n(y))}{|x-y|^{\frac{N+sp}{p'}}} \rightarrow \frac{|v(x) - v(y)|^{p-2}(v(x) - v(y))}{|x-y|^{\frac{N+sp}{p'}}} $$

a.e. in $\mathbb{R}^N \times \mathbb{R}^N$. It follows that

$$ \frac{|v_n(x) - v_n(y)|^{p-2}(v_n(x) - v_n(y))}{|x-y|^{\frac{N+sp}{p'}}} \rightarrow \frac{|v(x) - v(y)|^{p-2}(v(x) - v(y))}{|x-y|^{\frac{N+sp}{p'}}} $$

weakly in $L^{p'}(\mathbb{R}^N \times \mathbb{R}^N)$. Then, since $\varphi \in \mathbf{W} = W_0^{s,p}(\Omega) \cap L^{2q}(\Omega)$, we obtain

$$ \lim_{n \to \infty} \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|v_n(x) - v_n(y)|^{p-2}(v_n(x) - v_n(y))(\varphi(x) - \varphi(y))}{|x-y|^{N+sp}} dx dy \\ = \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|v(x) - v(y)|^{p-2}(v(x) - v(y))(\varphi(x) - \varphi(y))}{|x-y|^{N+sp}} dx dy. $$

With similar arguments, by the Hölder inequality, $(v_n^{2q-1})_{n \in \mathbb{N}}$ and $(h_n v_n^{q-1})_{n \in \mathbb{N}}$ are uniformly bounded in $L^{\frac{2q}{2q-1}}(\Omega)$. By (2.5), we infer that $f(x, v_n)$ are uniformly bounded in $L^{\frac{2q}{q-1}}(\Omega)$ and $f(x, v_n) \rightarrow f(x, v)$ a.e. in $\Omega$. Since $\varphi \in \mathbf{W} = W_0^{s,p}(\Omega) \cap L^{2q}(\Omega)$, we obtain

$$ \lim_{n \to \infty} \int_{\Omega} v_n^{2q-1} \varphi dx = \int_{\Omega} v^{2q-1} \varphi dx, \quad \lim_{n \to \infty} \int_{\Omega} h_n v_n^{q-1} \varphi dx = \int_{\Omega} h v^{q-1} \varphi dx, \\ \lim_{n \to \infty} \int_{\Omega} f(x, v_n) \varphi dx = \int_{\Omega} f(x, v) \varphi dx. $$

By passing to the limit in (2.16), $v$ is a weak solution to (2.1). Finally, the fact that $v \in L^{\infty}(\Omega)$ follows from Corollary 5.3. $\square$

From Theorem 5.4, we obtain the following result.

**Corollary 2.6.** Let $\lambda > 0$, $q \in (1, p]$, $f: \Omega \times \mathbb{R}^+ \rightarrow \mathbb{R}^+$ satisfy (H1), (H2), (H6). In addition suppose that $h_0 \in L^2(\Omega) \cap L^r(\Omega)$, for some $r > \frac{N}{sp}$ and satisfies (H8). Then, there exists a unique solution $u$ of problem (2.11). Namely, $u$ belongs to $\dot{V}_+^q \cap L^\infty(\Omega)$, satisfies (2.12) for any $\Psi$ satisfying (2.13) and there exists $c > 0$ such that $u(x) \ge cd^{sq}(x)$ a.e. in $\Omega$.

Moreover, if $u_1$ and $u_2$ are two solutions to the problem (2.11) associated with $h_1, h_2 \in L^2(\Omega)$ satisfy (H8), then

$$ \|(u_1 - u_2)^+\|_{L^2} \le \|(u_1 - u_2 + \lambda(\mathcal{T}_q(u_1) - \mathcal{T}_q(u_2)))^+\|_{L^2}. \quad (2.20) $$

*Proof.* The existence of a solution $v$ in Theorem 2.5 can be obtained by a global minimization argument as in Step 1 of the proof of Theorem 2.2. Therefore, we deduce from Theorem 5.4 that $v$ is a global minimizer of $\mathcal{J}$ defined in (2.4).

As in the proof of Corollary 2.4, we can define the energy functional $\xi$ on $\dot{V}_+^q \cap L^2(\Omega)$ as $\xi(u) = \mathcal{J}(u^{1/q})$. We set $u_0 = v_0^q$. Then, $u_0$ belongs to $\dot{V}_+^q \cap L^\infty(\Omega)$. By (2.18) we obtain $u_0(x) \ge cd^{sq}(x)$ a.e. in $\Omega$. Let $\Psi$ satisfy (2.13), then for $t$ small enough, $\xi(u_0 + t\Psi) - \xi(u_0) \ge 0$. By using the Taylor expansion, we deduce that $u_0$ satisfies (2.12). Finally, (2.15) gives (2.20). $\square$
### 3. EXISTENCE OF A WEAK SOLUTION TO PARABOLIC PROBLEM (1.1)

In light of Remark 1.11, we consider problem (1.12) and establish the existence of weak solution when $v_0 \in \mathcal{M}_{d^s}^1(\Omega) \cap W_0^{s,p}(\Omega)$. In this section, we prove Theorem 1.12. We begin the next subsection with some auxiliary results.

**3.1. Existence and regularity of a weak solution.** We divided the subsection into three main parts concerning: existence, uniqueness, and regularity of solutions.

*Existence of a weak solution.*

**Theorem 3.1.** Under the assumptions of Theorem 1.12, there exists a weak solution $v$ to the problem (1.12) (in sense of Definition 1.10). Furthermore, $v$ belongs to $C([0, T]; L^r(\Omega))$ for any $1 \le r < \infty$ and there exists $C > 0$ such that, for any $t \in [0, T]$:

$$C^{-1}d^s(x) \le v(t, x) \le Cd^s(x) \quad \text{a.e. in } \Omega. \quad (3.1)$$

*Proof.* We use the time semi-discretization method: Let $n^* \in \mathbb{N}^*$ and $T > 0$. We set $\Delta_t = \frac{T}{n^*}$ and for $n \in \{1, \dots, n^*\}$, we define $t_n = n\Delta_t$. We perform the proof along four main steps.

**Step 1: Approximation of $h$.** For $n \in \{1, \dots, n^*\}$, we define for $(t, x) \in [t_{n-1}, t_n) \times \Omega$,

$$h_{\Delta_t}(t, x) = h^n(x) := \frac{1}{\Delta_t} \int_{t_{n-1}}^{t_n} h(s, x)ds.$$

The Jensen's inequality implies that

$$\|h_{\Delta_t}\|_{L^2(Q_T)} \le \|h\|_{L^2(Q_T)}.$$

Hence $h_{\Delta_t} \in L^2(Q_T)$, $h^n \in L^2(\Omega)$. It is easy to prove by density arguments that

$$h_{\Delta_t} \to h \quad \text{in } L^2(Q_T).$$

On the other hand, we obtain

$$\|h_{\Delta_t}\|_{L^\infty(Q_T)} \le \|h\|_{L^\infty(Q_T)}.$$

**Step 2: Time discretization of problem (1.12).** We define the following implicit Euler scheme: $v^0 = v_0$ and for $n \ge 1$, $v_n$ is the weak solution of

$$\begin{gathered} \left(\frac{v_n^q - v_{n-1}^q}{\Delta_t}\right) v_n^{q-1} + (-\Delta)^s v_n = h^n v_n^{q-1} + f(x, v_n) \quad \text{in } \Omega, \tag{3.2} \\ v_n > 0 \quad \text{in } \Omega, \\ v_n = 0 \quad \text{in } \mathbb{R}^N \setminus \Omega. \end{gathered}$$

The sequence $(v_n)_{n=1,2,\dots,n^*}$ is well-defined. Indeed, existence and uniqueness of $v_1 \in C(\bar{\Omega}) \cap \mathcal{M}_{d^s}^1(\Omega)$ follow from Theorem 2.2 with $h_0 = \Delta_t h^1 + v_0^q \in L^\infty(\Omega)$ and $\Delta_t h^1 + v_0^q \ge \Delta_t \underline{h}$. Hence by induction we obtain in the same way the existence and the uniqueness of the solution $v_n$ for any $n = 2, 3, \dots, n^*$ where $v_n \in C(\bar{\Omega}) \cap \mathcal{M}_{d^s}^1(\Omega)$.

**Step 3: Existence of subsolutions and supersolutions.** In this step, we establish the existence of a sub-solution $w$ and a super-solution $\bar{w}$ such that $v_n \in [w, \bar{w}]$ for all $n \in \{0, 1, 2, \dots, n^*\}$. First, we rewrite (3.2) as

$$v_n^{2q-1} + \Delta_t (-\Delta)^s v_n = (\Delta_t h^n + v_{n-1}^q) v_n^{q-1} + \Delta_t f(x, v_n). \quad (3.3)$$
As in Theorem 2.2, we prove that for any $\mu \in (0, 1]$, the problem below admits a unique weak solution $\underline{w}_{\mu} \in C(\bar{\Omega}) \cap \mathcal{M}_{d^s}^1(\Omega)$,

$$ \begin{aligned} (-\Delta)_p^s w &= \mu(\underline{h}w^{q-1} + f(x, w)) \quad \text{in } \Omega, \\ w &\ge 0 \quad \text{in } \Omega, \\ w &= 0 \quad \text{in } \mathbb{R}^N \setminus \Omega, \end{aligned} \quad (3.4) $$

where $\underline{h}$ is defined in (H4).

Let $\mu_1 < \mu_2 \le 1$ and $\underline{w}_{\mu_1}, \underline{w}_{\mu_2} \in C(\bar{\Omega}) \cap \mathcal{M}_{d^s}^1(\Omega)$ be two weak solutions of (3.4). Then

$$ \begin{aligned} & \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|\underline{w}_{\mu_1}(x) - \underline{w}_{\mu_1}(y)|^{p-2}(\underline{w}_{\mu_1}(x) - \underline{w}_{\mu_1}(y))(\Phi(x) - \Phi(y))}{|x-y|^{N+sp}} dx dy \\ &= \mu_1 \int_{\Omega} (\underline{h} \underline{w}_{\mu_1}^{q-1} + f(x, \underline{w}_{\mu_1})) \Phi dx \end{aligned} $$

and

$$ \begin{aligned} & \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|\underline{w}_{\mu_2}(x) - \underline{w}_{\mu_2}(y)|^{p-2}(\underline{w}_{\mu_2}(x) - \underline{w}_{\mu_2}(y))(\Psi(x) - \Psi(y))}{|x-y|^{N+sp}} dx dy \\ &= \mu_2 \int_{\Omega} (\underline{h} \underline{w}_{\mu_2}^{q-1} + f(x, \underline{w}_{\mu_2})) \Psi dx. \end{aligned} $$

Subtracting the above expressions and taking

$$ \Phi = \frac{(\underline{w}_{\mu_1}^q - \underline{w}_{\mu_2}^q)^+}{\underline{w}_{\mu_1}^{q-1}}, \quad \Psi = \frac{(\underline{w}_{\mu_2}^q - \underline{w}_{\mu_1}^q)^-}{\underline{w}_{\mu_2}^{q-1}}, $$

we deduce that $(\underline{w}_{\mu})_{\mu}$ is nondecreasing. From [27, Corollary 4.2 and Theorem 1.1], we obtain for some $\mu_0 > 0$ and $0 < \alpha \le s$ that

$$ \|\underline{w}_{\mu}\|_{C^{0,\alpha}(\bar{\Omega})} \le C(\mu_0) \text{ for any } \mu \le \mu_0 \text{ and } \|\underline{w}_{\mu}\|_{L^{\infty}(\Omega)} \to 0 \text{ as } \mu \to 0. $$

Furthermore, by using [27, Theorem 4.4], we can choose $\mu < 1$ small enough such that there exists $\underline{w} \in C(\bar{\Omega}) \cap \mathcal{M}_{d^s}^1(\Omega)$ satisfies $0 < \underline{w} := \underline{w}_{\mu} \le v_0$. We infer that $\underline{w}$ is the subsolution of the problem (3.3) for $n = 1$, i.e.

$$ \begin{aligned} & \int_{\Omega} \underline{w}^{2q-1} \varphi dx + \Delta_t \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|\underline{w}(x) - \underline{w}(y)|^{p-2}(\underline{w}(x) - \underline{w}(y))(\varphi(x) - \varphi(y))}{|x-y|^{N+sp}} dx dy \\ & \le \Delta_t \int_{\Omega} (h^1 \underline{w}^{q-1} + f(x, \underline{w})) \varphi dx + \int_{\Omega} v_0^q \underline{w}^{q-1} \varphi dx \end{aligned} $$

for all $\varphi \in \mathbf{W}$ and $\varphi \ge 0$. We also recall that $v_1$ satisfies

$$ \begin{aligned} & \int_{\Omega} v_1^{2q-1} \psi dx + \Delta_t \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|v_1(x) - v_1(y)|^{p-2}(v_1(x) - v_1(y))(\psi(x) - \psi(y))}{|x-y|^{N+sp}} dx dy \\ &= \Delta_t \int_{\Omega} (h^1 v_1^{q-1} + f(x, v_1)) \psi dx + \int_{\Omega} v_0^q v_1^{q-1} \psi dx \end{aligned} $$

for all $\psi \in \mathbf{W}$. By Theorem 5.4, we obtain $\underline{w} \le v_1$ and then by induction $0 < \underline{w} \le v_n$ in $\Omega$ for $n = 0, 1, 2, \dots, n^*.$

Next, we construct a uniform supersolution. We start with the case $q < p$ for which we consider the problem

$$ \begin{aligned} (-\Delta)_p^s w &= 1 \quad \text{in } \Omega, \\ w &= 0 \quad \text{in } \mathbb{R}^N \setminus \Omega. \end{aligned} \quad (3.5) $$
As above, we can prove that there exists a unique weak solution $w \in C(\bar{\Omega}) \cap \mathcal{M}_{ds}^1(\Omega)$ to (3.5). We easily check that for some $K > 0$ fixed, $w_K = K^{\frac{1}{p-1}}w$ is the unique weak solution of the problem

$$(-\Delta)_p^s w_K = K \quad \text{in } \Omega, \\ w_K = 0 \quad \text{in } \mathbb{R}^N \setminus \Omega$$

and

$$c^{-1} d(x)^s K^{\frac{1}{p-1}} \le w_K(x) \le c d(x)^s K^{\frac{1}{p-1}}, \quad (3.6)$$

where $c > 0$ is a constant. Again by using [27, Theorem 4.4], we obtain $\bar{w} = w_K \ge v_0$ for $K$ large enough. By (2.5) and (3.6), it is easy to prove that $\bar{w}$ is the supersolution of the problem

$$(-\Delta)_p^s w = \|h\|_{L^\infty(\Omega)} w^{q-1} + f(x, w) \quad \text{in } \Omega, \\ w > 0 \quad \text{in } \Omega, \\ w = 0 \quad \text{in } \mathbb{R}^N \setminus \Omega. \quad (3.7)$$

We now study the case $q = p$. Using (H3), we can choose for any $\epsilon > 0$, $r_0 = r_0(\epsilon) > 0$ large enough, such that for $r \ge r_0$,

$$f(x, r) \le \epsilon r^{p-1}. \quad (3.8)$$

Let $w$ be the solution of the problem

$$(-\Delta)_p^s w = C + \beta w^{p-1} \quad \text{in } \Omega, \\ w > 0 \quad \text{in } \Omega, \\ w = 0 \quad \text{in } \mathbb{R}^N \setminus \Omega$$

with $C > 0$ and $\beta < \lambda_{1,p,s}$. Then, by a similar proof as in Theorem 2.2 step 2, we obtain $w \in L^\infty(\Omega)$. On the other hand, by [20, Theorems 1.4 and 1.5, p. 768], we obtain that $w > 0$ in $\Omega$ and satisfies $w \ge kd^s(x)$, for some $k = k(C, \beta) > 0$. Finally, using [27, Theorem 4.4], we obtain that $w \in \mathcal{M}_{ds}^1(\Omega)$. By (3.8), (H5) and for $C > 0$ large enough and $\beta$ close enough to $\lambda_{1,p,s}$, we obtain

$$(-\Delta)_p^s(w) = C + \beta w^{p-1} \ge \|h\|_{L^\infty(\Omega)} w^{p-1} + f(x, w).$$

Hence, $\bar{w} = w$ is supersolution of (3.7). Again using [27, Theorem 4.4]] and taking $C > 0$ large enough, we obtain $v_0 \le \bar{w}$.

Then, since $\bar{w} \ge v_0$, $\bar{w}$ is the supersolution to (3.3) for $n = 1$, i.e.

$$\int_{\Omega} \bar{w}^{2q-1} \varphi dx + \Delta_t \int_{\mathbb{R}^n} \int_{\mathbb{R}^n} \frac{|\bar{w}(x) - \bar{w}(y)|^{p-2} (\bar{w}(x) - \bar{w}(y))(\varphi(x) - \varphi(y))}{|x-y|^{N+sp}} dx dy \\ \ge \Delta_t \int_{\Omega} (h^1 \bar{w}^{q-1} + f(x, \bar{w})) \varphi dx + \int_{\Omega} v_0^q \bar{w}^{q-1} \varphi dx$$

for all $\varphi \in \mathbf{W}$ and $\varphi \ge 0$. From Theorem 5.4, we obtain $\bar{w} \ge v_1$ and then by induction we have $\bar{w} \ge v_n$ for all $n = 1, 2, 3, \dots, n^*$. Finally, we conclude that $\bar{w} \le v_n \le \bar{w}$ for $n = 0, 1, 2, 3, \dots, n^*$, i.e. $c_1 d^s(x) \le v_n(x) \le c_2 d^s(x)$ in $\Omega$, where $c_1, c_2$ are positive constants independent of $n$.

**Step 3: A priori estimates.** For $n \in \{1, 2, 3, \dots, n^*\}$ and $t \in [t_{n-1}, t_n)$ let the functions $v_{\Delta_t}(t)$ and $\tilde{v}_{\Delta_t}(t)$ be as follows:

$$v_{\Delta_t}(t) = v_n,$$
$$ \tilde{v}_{\Delta_t}(t) = \frac{(t - t_{n-1})}{\Delta_t} (v_n^q - v_{n-1}^q) + v_{n-1}^q. $$

One can easily check that

$$ v_{\Delta_t}^{q-1} \frac{\partial \tilde{v}_{\Delta_t}}{\partial t} + (-\Delta)_p^s v_{\Delta_t} = h^n v_{\Delta_t}^{q-1} + f(x, v_{\Delta_t}). \quad (3.9) $$

We observe now that as $\Delta_t \to 0$, the discrete equation (3.9) converges to (1.12). We further point out that there exists $c > 0$ independent of $\Delta_t$ such that for any $(t, x) \in Q_T$,

$$ c^{-1} d^s(x) \le v_{\Delta_t}, \quad \tilde{v}_{\Delta_t}^{1/q} \le cd^s(x). \quad (3.10) $$

Now, multiplying (3.2) by $\frac{v_n^q - v_{n-1}^q}{v_n^{q-1}} \in W_0^{s,p}(\Omega) \cap L^\infty(\Omega)$ and summing from $n=1$ to $n' \le n^*$, we obtain

$$ \sum_{n=1}^{n'} \int_{\Omega} \Delta_t \left( \frac{v_n^q - v_{n-1}^q}{\Delta_t} \right)^2 dx + \sum_{n=1}^{n'} \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|v_n(x) - v_n(y)|^{p-2}(v_n(x) - v_n(y))}{|x-y|^{N+sp}} \\ \times \left[ \left( \frac{v_n^q - v_{n-1}^q}{v_n^{q-1}} \right)(x) - \left( \frac{v_n^q - v_{n-1}^q}{v_n^{q-1}} \right)(y) \right] dx dy \\ = \sum_{n=1}^{n'} \int_{\Omega} h^n (v_n^q - v_{n-1}^q) dx + \sum_{n=1}^{n'} \int_{\Omega} \frac{f(x, v_n)}{v_n^{q-1}} (v_n^q - v_{n-1}^q) dx. $$

Since $v_n \in [\underline{w}, \overline{w}] \subset \mathcal{M}_{ds}^1(\Omega)$, we have that $(\frac{f(x, v_n)}{v_n^{q-1}}(v_n^q - v_{n-1}^q))$ is uniformly bounded. By Young's inequality, we have

$$ \sum_{n=1}^{n'} \int_{\Omega} \Delta_t \left( \frac{v_n^q - v_{n-1}^q}{\Delta_t} \right)^2 dx + \sum_{n=1}^{n'} \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|v_n(x) - v_n(y)|^{p-2}(v_n(x) - v_n(y))}{|x-y|^{N+sp}} \\ \times \left[ \left( \frac{v_n^q - v_{n-1}^q}{v_n^{q-1}} \right)(x) - \left( \frac{v_n^q - v_{n-1}^q}{v_n^{q-1}} \right)(y) \right] dx dy \\ \le \frac{1}{2} \sum_{n=1}^{n'} \Delta_t \|h^n\|_{L^2(\Omega)}^2 + \frac{1}{2} \sum_{n=1}^{n'} \int_{\Omega} \Delta_t \left( \frac{v_n^q - v_{n-1}^q}{\Delta_t} \right)^2 dx + C, $$

i.e.

$$ \frac{1}{2} \sum_{n=1}^{n'} \int_{\Omega} \Delta_t \left( \frac{v_n^q - v_{n-1}^q}{\Delta_t} \right)^2 dx + \sum_{n=1}^{n'} \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|v_n(x) - v_n(y)|^{p-2}(v_n(x) - v_n(y))}{|x-y|^{N+sp}} \\ \times \left[ \left( \frac{v_n^q - v_{n-1}^q}{v_n^{q-1}} \right)(x) - \left( \frac{v_n^q - v_{n-1}^q}{v_n^{q-1}} \right)(y) \right] dx dy \\ \le \frac{1}{2} \sum_{n=1}^{n'} \Delta_t \|h^n\|_{L^2(\Omega)}^2 + C, $$

where $C$ is independent of $n'$. Then by step 1, we obtain

$$ \left( \frac{\partial \tilde{v}_{\Delta_t}}{\partial t} \right) \text{ is bounded in } L^2(Q_T) \text{ uniformly in } \Delta_t. \quad (3.11) $$
Now, from Proposition 1.2 and by Young's inequality in the case $q < p$, we have

$$ 
\begin{aligned}
& |v_n(x) - v_n(y)|^{p-2}(v_n(x) - v_n(y)) \left[ \frac{v_{n-1}(x)^q}{v_n(x)^{q-1}} - \frac{v_{n-1}(y)^q}{v_n(y)^{q-1}} \right] \\
& \le |v_{n-1}(x) - v_{n-1}(y)|^q |v_n(x) - v_n(y)|^{p-q} \quad (3.12) \\
& \le \frac{q}{p} |v_{n-1}(x) - v_{n-1}(y)|^p + \frac{p-q}{p} |v_n(x) - v_n(y)|^p.
\end{aligned}
$$

Next, for $p = q$ we obtain

$$ 
\begin{aligned}
& |v_n(x) - v_n(y)|^{p-2}(v_n(x) - v_n(y)) \left[ \frac{v_{n-1}(x)^p}{v_n(x)^{p-1}} - \frac{v_{n-1}(y)^p}{v_n(y)^{p-1}} \right] \\
& \le |v_{n-1}(x) - v_{n-1}(y)|^p. \quad (3.13)
\end{aligned}
$$

Then, for any $n' \ge 1$ and $p \ne q$ we obtain

$$ 
\begin{aligned}
& \frac{1}{2} \sum_{n=1}^{n'} \Delta_t \|h^n\|_{L^2(\Omega)}^2 + C \\
& \ge \sum_{n=1}^{n'} \left[ \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|v_n(x) - v_n(y)|^p}{|x-y|^{N+sp}} dx dy - \frac{q}{p} \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|v_{n-1}(x) - v_{n-1}(y)|^p}{|x-y|^{N+sp}} dx dy \right. \\
& \qquad \left. - \frac{p-q}{p} \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|v_n(x) - v_n(y)|^p}{|x-y|^{N+sp}} dx dy \right].
\end{aligned}
$$

For $p = q$, we have

$$ 
\begin{aligned}
& \frac{1}{2} \sum_{n=1}^{n'} \Delta_t \|h^n\|_{L^2(\Omega)}^2 + C \ge \sum_{n=1}^{n'} \left[ \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|v_n(x) - v_n(y)|^p}{|x-y|^{N+sp}} dx dy \right. \\
& \qquad \left. - \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|v_{n-1}(x) - v_{n-1}(y)|^p}{|x-y|^{N+sp}} dx dy \right].
\end{aligned}
$$

Finally, we obtain

$$ 
\begin{aligned}
& \frac{1}{2} \sum_{n=1}^{n'} \Delta_t \|h^n\|_{L^2(\Omega)}^2 + C \ge \frac{q}{p} \left[ \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|v_{n'}(x) - v_{n'}(y)|^p}{|x-y|^{N+sp}} dx dy \right. \\
& \qquad \left. - \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|v_0(x) - v_0(y)|^p}{|x-y|^{N+sp}} dx dy \right]
\end{aligned}
$$

which implies that

$$ (v_{\Delta t}) \text{ is bounded in } L^{\infty}(0, T; W_0^{s,p}(\Omega)) \text{ uniformly in } \Delta t. \quad (3.14) $$

Since $\tilde{v}_{\Delta t} = \xi v_n^q + (1-\xi)v_{n-1}^q$, where $\xi = \frac{t-t_{n-1}}{\Delta t}$, by Proposition 1.1, we obtain

$$ 
\begin{aligned}
& \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|\tilde{v}_{\Delta t}^{1/q}(x) - \tilde{v}_{\Delta t}^{1/q}(y)|^p}{|x-y|^{N+sp}} dx dy \\
& \le \xi \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|v_n(x) - v_n(y)|^p}{|x-y|^{N+sp}} + (1-\xi) \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|v_{n-1}(x) - v_{n-1}(y)|^p}{|x-y|^{N+sp}}.
\end{aligned}
$$

Then, we conclude that

$$ (\tilde{v}_{\Delta t}^{1/q}) \text{ is bounded in } L^{\infty}(0, T; W_0^{s,p}(\Omega)) \text{ uniformly in } \Delta t. \quad (3.15) $$
Thus, $v_{\Delta t} \stackrel{*}{\sim} v$ and $\tilde{v}_{\Delta t}^{1/q} \stackrel{*}{\sim} \tilde{v}$ in $L^{\infty}(0, T; W_0^{s,p}(\Omega))$. Furthermore using (2.17), (3.11),

$$ \sup_{t \in [0, T]} \| \tilde{v}_{\Delta t}^{1/q} - v_{\Delta t} \|_{L^{2q}(\Omega)}^{2q} \le \sup_{t \in [0, T]} \| \tilde{v}_{\Delta t} - v_{\Delta t}^q \|_{L^2(\Omega)}^2 \le C \Delta t \to 0 \quad \text{as} \quad \Delta t \to 0. \quad (3.16) $$

It follows that $v = \tilde{v}$.

Now, from (3.11), (3.15) and since $W_0^{s,p}(\Omega) \hookrightarrow L^r(\Omega)$ compactly for all $1 \le r < p_s^*$, using Theorem 1.4 we obtain that $(\tilde{v}_{\Delta t})$ is compact in $C([0, T]; L^r(\Omega))$. Then from (3.16),

$$ \tilde{v}_{\Delta t} \to v^q \quad \text{in} \quad C([0, T]; L^r(\Omega)), \quad \text{for} \quad 1 \le r < p_s^*. $$

Using $\tilde{v}_{\Delta t} \in L^{\infty}(\Omega)$ with the interpolation inequality with $p_s^* \le r < \infty$,

$$ \| \cdot \|_r \le \| \cdot \|_{\infty}^{\alpha} \| \cdot \|_{p_s^*}^{1-\alpha}, \quad \text{with} \quad \alpha \in [0, 1], $$

we obtain that

$$ \tilde{v}_{\Delta t} \to v^q \quad \text{in} \quad C([0, T]; L^r(\Omega)), \quad \text{for all} \ r \ge 1. \quad (3.17) $$

Hence, from the estimate

$$ \sup_{t \in [0, T]} \|v_{\Delta t}^q - \tilde{v}_{\Delta t}\|_{L^2(\Omega)} \le C(\Delta t)^{1/2}, \quad (3.18) $$

we have

$$ v_{\Delta t} \to v \quad \text{in} \quad L^{\infty}([0, T]; L^r(\Omega)), \quad \text{for all} \ r \ge 1. \quad (3.19) $$

Hence, (3.10) implies (3.1). From (3.11) and (3.17), we obtain

$$ \frac{\partial \tilde{v}_{\Delta t}}{\partial t} \to \frac{\partial v^q}{\partial t} \quad \text{in} \quad L^2(Q_T). \quad (3.20) $$

**Step 4: $v$ satisfies (1.12).** • First, from (3.14), we have

$$ \left\{ \frac{|v_{\Delta t}(t, x) - v_{\Delta t}(t, y)|^{p-2} (v_{\Delta t}(t, x) - v_{\Delta t}(t, y))}{|x-y|^{\frac{N+s_p}{p'}}} \right\} $$

is bounded in $L^{\infty}(0, T; L^p(\mathbb{R}^N \times \mathbb{R}^N))$, where $p' = \frac{p}{p-1}$, and by the pointwise convergence of $v_{\Delta t}$ to $v$, we obtain as $\Delta t \to 0^+$ and for a.e. $t \in [0, T]$,

$$ \frac{|v_{\Delta t}(t, x) - v_{\Delta t}(t, y)|^{p-2} (v_{\Delta t}(t, x) - v_{\Delta t}(t, y))}{|x-y|^{\frac{N+s_p}{p'}}} \to \frac{|v(t, x) - v(t, y)|^{p-2} (v(t, x) - v(t, y))}{|x-y|^{\frac{N+s_p}{p'}}} $$

a.e. in $\mathbb{R}^N \times \mathbb{R}^N$, it follows that as $\Delta t \to 0^+$,

$$ \frac{|v_{\Delta t}(t, x) - v_{\Delta t}(t, y)|^{p-2} (v_{\Delta t}(t, x) - v_{\Delta t}(t, y))}{|x-y|^{\frac{N+s_p}{p'}}} \to \frac{|v(t, x) - v(t, y)|^{p-2} (v(t, x) - v(t, y))}{|x-y|^{\frac{N+s_p}{p'}}} $$
weakly in $L^p((0,T) \times \mathbb{R}^{2N})$. Then, we conclude that for any $\phi \in C_c^\infty(Q_T)$,

$$ \begin{aligned} & \int_0^T \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \left( |v_{\Delta_t}(t,x) - v_{\Delta_t}(t,y)|^{p-2} (v_{\Delta_t}(t,x) - v_{\Delta_t}(t,y)) \right. \\ & \qquad \left. \times (\varphi(t,x) - \varphi(t,y)) \right) / |x-y|^{N+sp} dx dy dt \\ & \rightarrow \int_0^T \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \left( |v(t,x) - v(t,y)|^{p-2} (v(t,x) - v(t,y)) \right. \\ & \qquad \left. \times (\varphi(t,x) - \varphi(t,y)) \right) / |x-y|^{N+sp} dx dy dt. \end{aligned} \quad (3.21) $$

• Next, from (2.17), (3.16) and (3.18) we have

$$ \begin{aligned} & \|v_{\Delta_t}^{q-1} - v^{q-1}\|_{L^2(Q_T)} \\ & \le T^{\frac{1}{2}} \|v_{\Delta_t}^{q-1} - v^{q-1}\|_{L^\infty(0,T;L^2(\Omega))} \\ & \le T^{\frac{1}{2}} |\Omega|^{\frac{1}{2q}} \|v_{\Delta_t}^{q-1} - v^{q-1}\|_{L^\infty(0,T;L^{\frac{2q}{q-1}}(\Omega))} \quad (3.22) \\ & \le T^{\frac{1}{2}} |\Omega|^{\frac{1}{2q}} \|v_{\Delta_t}^{q} - v^{q}\|_{L^\infty(0,T;L^2(\Omega))} \\ & \le T^{\frac{1}{2}} |\Omega|^{\frac{1}{2q}} \left[ \|v_{\Delta_t}^{q} - \tilde{v}_{\Delta_t}\|_{L^\infty(0,T;L^2(\Omega))} + \|\tilde{v}_{\Delta_t} - v^{q}\|_{L^\infty(0,T;L^2(\Omega))} \right]^{\frac{q-1}{q}} \rightarrow 0 \end{aligned} $$

as $\Delta_t \rightarrow 0$. By the Hölder inequality, for all $\varphi \in C_c^\infty(Q_T)$ we have

$$ \begin{aligned} & \left| \int_0^T \int_{\Omega} (v_{\Delta_t}^{q-1} \frac{\partial \tilde{v}_{\Delta_t}}{\partial t} - v^{q-1} \frac{\partial v^q}{\partial t}) \varphi(t,x) dx dt \right| \\ & \le \left| \int_0^T \int_{\Omega} v^{q-1} \left( \frac{\partial \tilde{v}_{\Delta_t}}{\partial t} - \frac{\partial v^q}{\partial t} \right) \varphi(t,x) dx dt \right| + \|v_{\Delta_t}^{q-1} - v^{q-1}\|_{L^2(Q_T)} \left\| \frac{\partial \tilde{v}_{\Delta_t}}{\partial t} \varphi \right\|_{L^2(Q_T)} \end{aligned} $$

and

$$ \begin{aligned} & \int_0^T \int_{\Omega} (h^n v_{\Delta_t}^{q-1} - h v^{q-1}) \varphi \, dx \, dt \\ & = \int_0^T \int_{\Omega} h^n (v_{\Delta_t}^{q-1} - v^{q-1}) \varphi \, dx \, dt + \int_0^T \int_{\Omega} (h^n - h) v^{q-1} \varphi \, dx \, dt \\ & \le \|h^n \varphi\|_{L^2(Q_T)} \|v_{\Delta_t}^{q-1} - v^{q-1}\|_{L^2(Q_T)} + \|v^{q-1} \varphi\|_{L^2(Q_T)} \|h^n - h\|_{L^2(Q_T)}. \end{aligned} $$

Then from (3.11), (3.19), (3.20), (3.22) and Step 1, we obtain

$$ \int_0^T \int_{\Omega} (v_{\Delta_t}^{q-1} \frac{\partial \tilde{v}_{\Delta_t}}{\partial t} - v^{q-1} \frac{\partial v^q}{\partial t}) \varphi(t,x) dx dt \rightarrow 0, \quad (3.23) $$

$$ \int_0^T \int_{\Omega} (h^n v_{\Delta_t}^{q-1} - h v^{q-1}) \varphi(t,x) dx dt \rightarrow 0 \quad (3.24) $$

as $\Delta_t \rightarrow 0$. From (3.19), we have $f(x, v_{\Delta_t})\varphi \rightarrow f(x, v)\varphi$ a.e. in $Q_T$, (up to a subsequence). Furthermore from (2.5) and (3.10), $(f(x, v_{\Delta_t}))$ is bounded in $L^2(Q_T)$ uniformly in $\Delta_t$. Then, by the dominated convergence Theorem, we obtain

$$ \int_0^T \int_{\Omega} f(x, v_{\Delta_t}) \varphi \, dx \, dt \rightarrow \int_0^T \int_{\Omega} f(x, v) \varphi \, dx \, dt, \quad \text{as } \Delta_t \rightarrow 0. \quad (3.25) $$
Finally, gathering (3.21), (3.23), (3.24), (3.25) and passing to the limit in (3.9) as $\Delta_t \to 0^+$, we conclude that $v$ satisfies (1.12), i.e.

$$ 
\begin{aligned}
& \int_0^T \int_{\Omega} \partial_t (v^q) v^{q-1} \varphi \, dx \, ds \\
& + \int_0^T \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|v(s, x) - v(s, y)|^{p-2} (v(s, x) - v(s, y)) (\varphi(s, x) - \varphi(s, y))}{|x-y|^{N+sp}} \, dx \, dy \, ds \\
& = \int_0^T \int_{\Omega} (h(s, x) v^{q-1} + f(x, v)) \varphi \, dx \, ds
\end{aligned} 
\quad (3.26) $$

for any $\varphi \in C_c^\infty(Q_T)$. Since $C_c^\infty(Q_T)$ is dense in $L^2(Q_T) \cap L^1(0, T; W_0^{s,p}(\Omega))$. Hence, we conclude that (3.26) is satisfied for any $\varphi \in L^2(Q_T) \cap L^1(0, T; W_0^{s,p}(\Omega))$. $\square$

Uniqueness.

*Proof of Theorem 1.13.* We again use the Picone identity. Let $v$ and $w$ be two weak solutions to (1.12) with $h$ and $\tilde{h}$ respectively. For $\epsilon \in (0; 1)$, we set

$$ \Phi := \frac{(v + \epsilon)^q - (w + \epsilon)^q}{(v + \epsilon)^{q-1}}, \quad \Psi := \frac{(w + \epsilon)^q - (v + \epsilon)^q}{(w + \epsilon)^{q-1}}. \quad (3.27) $$

$\Phi$ and $\Psi$ belong to $L^2(Q_T) \cap L^1(0, T; W_0^{s,p}(\Omega))$ and for any $t \in (0, T]$,

$$ 
\begin{aligned}
& \int_0^t \int_{\Omega} \partial_t (v^q) v^{q-1} \Phi \, dx \, ds \\
& + \int_0^t \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|v(s, x) - v(s, y)|^{p-2} (v(s, x) - v(s, y)) (\Phi(s, x) - \Phi(s, y))}{|x - y|^{N+sp}} \, dx \, dy \, ds \\
& = \int_0^t \int_{\Omega} (h(s, x) v^{q-1} + f(x, v)) \Phi \, dx \, ds
\end{aligned} 
$$

and

$$ 
\begin{aligned}
& \int_0^t \int_{\Omega} \partial_t (w^q) w^{q-1} \Psi \, dx \, ds \\
& + \int_0^t \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|w(s, x) - w(s, y)|^{p-2} (w(s, x) - w(s, y)) (\Psi(s, x) - \Psi(s, y))}{|x - y|^{N+sp}} \, dx \, dy \, ds \\
& = \int_0^t \int_{\Omega} (\tilde{h}(s, x) w^{q-1} + f(x, w)) \Psi \, dx \, ds.
\end{aligned} 
$$

Summing the above equalities, we obtain $\mathbf{I}_\epsilon = \mathbf{J}_\epsilon$ where

$$ 
\begin{aligned}
\mathbf{I}_{\epsilon} &= \int_{0}^{t} \int_{\Omega} \left( \frac{\partial_t (v^q) v^{q-1}}{(v+\epsilon)^{q-1}} - \frac{\partial_t (w^q) w^{q-1}}{(w+\epsilon)^{q-1}} \right) ((v+\epsilon)^q - (w+\epsilon)^q) \, dx \, ds \\
&\quad + \int_{0}^{t} \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|v(s, x) - v(s, y)|^{p-2} (v(s, x) - v(s, y))}{|x-y|^{N+sp}} \\
&\quad \times \left[ \frac{(v+\epsilon)^q(s, x) - (w+\epsilon)^q(s, x)}{(v+\epsilon)^{q-1}(s, x)} - \frac{(v+\epsilon)^q(s, y) - (w+\epsilon)^q(s, y)}{(v+\epsilon)^{q-1}(s, y)} \right] \, dx \, dy \, ds \\
&\quad + \int_{0}^{t} \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|w(s, x) - w(s, y)|^{p-2} (w(s, x) - w(s, y))}{|x-y|^{N+sp}} \\
&\quad \times \left[ \frac{(w+\epsilon)^q(s, x) - (v+\epsilon)^q(s, x)}{(w+\epsilon)^{q-1}(s, x)} - \frac{(w+\epsilon)^q(s, y) - (v+\epsilon)^q(s, y)}{(w+\epsilon)^{q-1}(s, y)} \right] \, dx \, dy \, ds
\end{aligned} 
$$
and

$$ \mathbf{J}_{\epsilon} = \int_{0}^{t} \int_{\Omega} \left( \frac{hv^{q-1}}{(v+\epsilon)^{q-1}} - \frac{\tilde{h}w^{q-1}}{(w+\epsilon)^{q-1}} \right) ((v+\epsilon)^q - (w+\epsilon)^q) dx ds \\ + \int_{0}^{t} \int_{\Omega} \left( \frac{f(x,v)}{(v+\epsilon)^{q-1}} - \frac{f(x,w)}{(w+\epsilon)^{q-1}} \right) ((v+\epsilon)^q - (w+\epsilon)^q) dx ds. $$

First, we deal with $\mathbf{I}_{\epsilon}$. Since $\frac{v}{v+\epsilon}$, $\frac{w}{w+\epsilon} < 1$ and $v, w \in L^{\infty}(Q_T)$,

$$ \left| \frac{\partial_t(v^q)v^{q-1}}{(v+\epsilon)^{q-1}} - \frac{\partial_t(w^q)w^{q-1}}{(w+\epsilon)^{q-1}} \right| \left| (v+\epsilon)^q - (w+\epsilon)^q \right| \le C \left( |\partial_t(v^q)| + |\partial_t(w^q)| \right) $$

where C does not depend on $\epsilon$. Moreover as $\epsilon \to 0^+$,

$$ \left( \frac{\partial_t(v^q)v^{q-1}}{(v+\epsilon)^{q-1}} - \frac{\partial_t(w^q)w^{q-1}}{(w+\epsilon)^{q-1}} \right) \left| (v+\epsilon)^q - (w+\epsilon)^q \right| \to \frac{1}{2} \partial_t(v^q - w^q)^2. $$

Therefore, by the dominated convergence Theorem and Lemma 1.8, we obtain

$$ \liminf_{\epsilon \to 0} \mathbf{I}_{\epsilon} \ge \frac{1}{2} \int_{0}^{t} \int_{\Omega} \partial_t (v^q - w^q)^2 dx ds. $$

Next, dealing with $\mathbf{J}_{\epsilon}$, dominated convergence Theorem implies

$$ \lim_{\epsilon \to 0} \int_{0}^{t} \int_{\Omega} \left( \frac{hv^{q-1}}{(v+\epsilon)^{q-1}} - \frac{\tilde{h}w^{q-1}}{(w+\epsilon)^{q-1}} \right) ((v+\epsilon)^q - (w+\epsilon)^q) dx ds \\ = \int_{0}^{t} \int_{\Omega} (h - \tilde{h})(v^q - w^q) dx ds. $$

Moreover, by using Fatou's Lemma, we have

$$ \liminf_{\epsilon \to 0} \int_{0}^{t} \int_{\Omega} \frac{f(x,v)}{(v+\epsilon)^{q-1}} (w+\epsilon)^q dx ds \ge \int_{0}^{t} \int_{\Omega} \frac{f(x,v)}{v^{q-1}} w^q dx ds, $$

$$ \liminf_{\epsilon \to 0} \int_{0}^{t} \int_{\Omega} \frac{f(x, w)}{(w+\epsilon)^{q-1}} (v+\epsilon)^q dx ds \ge \int_{0}^{t} \int_{\Omega} \frac{f(x, w)}{w^{q-1}} v^q dx ds. $$

Hence gathering the three limits above and from (H2), we obtain

$$ \liminf_{\epsilon \to 0} \mathbf{J}_{\epsilon} \le \int_{0}^{t} \int_{\Omega} (h - \tilde{h})(v^q - w^q) dx ds. $$

Since $\mathbf{I}_{\epsilon} = \mathbf{J}_{\epsilon}$, using Hölder inequality we conclude that for any $t \in [0, T]$,

$$ \frac{1}{2} \int_{0}^{t} \int_{\Omega} \partial_t (v^q - w^q)^2 dx ds \le \int_{0}^{t} \|h - \tilde{h}\|_{L^2(\Omega)} \|v^q - w^q\|_{L^2(\Omega)} ds $$

and by Grönwall Lemma [10, Lemma A.5], we deduce (1.14). □

The uniqueness of the weak solution in sense of Definition 1.10 in Theorem 1.12 is a consequence of Theorem 1.13. Precisely, we have the following Corollary.

**Corollary 3.2.** Let $v, w$ be weak solutions of (1.12) in sense of Definition 1.10 with the initial data $v_0 \in L^{2q}(\Omega)$, $v_0 \ge 0$ and $h \in L^2(Q_T)$. Then, $v \equiv w$.

We use Theorem 3.1 and Corollary 3.2 to infer the existence result concerning the parabolic problem involving the operator $\mathcal{T}_q$.
**Theorem 3.3.** Under the assumptions of Theorem 1.12, for any the initial data $u_0$ such that $u_0^{1/q} \in \mathcal{M}_d^1(\Omega) \cap W_0^{s,p}(\Omega)$, there exists a unique weak solution $u \in L^\infty(Q_T)$ of the problem

$$ 
\begin{array}{l}
\partial_t u + \mathcal{T}_q u = h \quad \text{in } Q_T, \\
u > 0 \quad \text{in } Q_T, \\
u = 0 \quad \text{on } \Gamma_T, \\
u(0, \cdot) = u_0 \quad \text{in } \Omega.
\end{array} 
\qquad (3.28) $$

In particular,

* $u^{1/q} \in L^\infty(0, T; W_0^{s,p}(\Omega))$, $\partial_t u \in L^2(Q_T)$;

* there exists $c > 0$ such that for any $t \in [0, T]$;

$$ c^{-1}d^s(x) \le u^{1/q}(t, x) \le cd^s(x) \quad \text{a.e. in } \Omega; $$

* for any $t \in [0, T]$, $u$ satisfies

$$ 
\begin{aligned}
& \int_0^t \int_{\Omega} \partial_t u \Psi \, dx \, ds + \int_0^t \int_{\mathbb{R}^{2N}} |u^{1/q}(s, x) - u^{1/q}(s, y)|^{p-2} (u^{1/q}(s, x) - u^{1/q}(s, y)) \\
& \times \left( \left( u^{\frac{1-q}{q}} \Psi \right)(s, x) - \left( u^{\frac{1-q}{q}} \Psi \right)(s, y) \right) / |x - y|^{N+sp} \, dx \, dy \, ds \\
& = \int_0^t \int_{\Omega} h(s, x) \Psi \, dx \, ds + \int_0^t \int_{\Omega} f(x, u^{1/q}) u^{\frac{1-q}{q}} \Psi \, dx \, ds
\end{aligned}
$$

for any $\Psi \in L^2(Q_T)$ such that

$$ |\Psi|^{1/q} \in L^1(0, T; W_0^{s,p}(\Omega)) \cap L^\infty(0, T; L_{ds}^\infty(\Omega)). $$

Moreover, for any $1 \le r < \infty$, $u$ belongs to $C([0, T]; L^r(\Omega))$.

The proof of the above theorem follows straightforward from Theorem 3.1 and corollary 3.2.

Regularity of weak solutions.

**Theorem 3.4.** Under the assumptions of Theorem 1.12, the weak solution $v$, of (1.12) obtained by Theorem 3.1, belongs to $C(0, T; W_0^{s,p}(\Omega))$ and for any $t \in [0, T]$ satisfies

$$ 
\begin{aligned}
& \int_0^t \int_{\Omega} \left(\frac{\partial v^q}{\partial t}\right)^2 dx \, ds + \frac{q}{p} \|v(t)\|_{W_0^{s,p}(\Omega)}^p \\
& = \int_0^t \int_{\Omega} h\left(\frac{\partial v^q}{\partial t}\right) dx \, ds + \int_0^t \int_{\Omega} \frac{f(x, v)}{v^{q-1}} \frac{\partial v^q}{\partial t} dx \, ds + \frac{q}{p} \|v_0\|_{W_0^{s,p}(\Omega)}^p.
\end{aligned}
$$

*Proof.* Since $v \in L^\infty(0, T; W_0^{s,p}(\Omega)) \cap L^\infty(Q_T)$ and $\partial_t v^q \in L^2(Q_T)$, by Theorem 1.4, we obtain that $v$ belongs to $C([0, T]; L^r(\Omega))$ for any $r \ge 1$. From the Sobolev embedding (Theorem 1.3), we have that $W_0^{s,p}(\Omega)$ is compactly embedded in $L^p(\Omega)$. So we deduce that $v: [0, T] \to W_0^{s,p}(\Omega)$ is weakly continuous. Therefore, for any $t_0 \in [0, T]$,

$$ \|v(t_0)\|_{W_0^{s,p}(\Omega)} \le \liminf_{t \to t_0} \|v(t)\|_{W_0^{s,p}(\Omega)}. \qquad (3.29) $$
Multiplying (3.2) by $\frac{v_n^q - v_{n-1}^q}{v_n^{q-1}} \in W$, integrating over $\mathbb{R}^N$ and summing from $1 \le n = N'$ to $N'' \le n^*$, we obtain

$$ \begin{aligned} & \sum_{n=N'}^{N''} \int_{\Omega} \Delta_t \left( \frac{v_n^q - v_{n-1}^q}{\Delta_t} Big \right)^2 dx + \sum_{n=N'}^{N''} \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|v_n(x) - v_n(y)|^{p-2}(v_n(x) - v_n(y))}{|x-y|^{N+sp}} \\ & \times \left[ \left( \frac{v_n^q - v_{n-1}^q}{v_n^{q-1}} \right)(x) - \left( \frac{v_n^q - v_{n-1}^q}{v_n^{q-1}} \right)(y) \right] dx dy \\ &= \sum_{n=N'}^{N''} \Delta_t \int_{\Omega} h^n \left( \frac{v_n^q - v_{n-1}^q}{\Delta_t} \right) dx + \Delta_t \sum_{n=N'}^{N''} \int_{\Omega} \frac{f(x, v_n)}{v_n^{q-1}} \left( \frac{v_n^q - v_{n-1}^q}{\Delta_t} \right) dx. \end{aligned} $$

Now, from (3.12) and (3.13), we obtain

$$ \begin{aligned} & \sum_{n=N'}^{N''} \int_{\Omega} \Delta_t \left( \frac{v_n^q - v_{n-1}^q}{\Delta_t} \right)^2 dx + \frac{q}{p} \left( \|v_{N'}\|_{W_0^{s,p}(\Omega)}^p - \|v_{N''-1}\|_{W_0^{s,p}(\Omega)}^p \right) \\ & \le \sum_{n=N''}^{N'} \Delta_t \int_{\Omega} h^n \left( \frac{v_n^q - v_{n-1}^q}{\Delta_t} \right) dx + \Delta_t \sum_{n=N'}^{N''} \int_{\Omega} \frac{f(x, v_n)}{v_n^{q-1}} \left( \frac{v_n^q - v_{n-1}^q}{\Delta_t} \right) dx. \end{aligned} \quad (3.30) $$

For any $t \in [t_0, T]$, we choose $N'$ and $N''$ such that $N'\Delta_t \to t$ and $N''\Delta_t \to t_0$. By (H7), then (3.30) gives

$$ \begin{aligned} & \int_{t_0}^{t} \int_{\Omega} \left(\frac{\partial v^q}{\partial t}\right)^2 dx ds + \frac{q}{p} \|v(t)\|_{W_0^{s,p}(\Omega)}^p \\ & \le \int_{t_0}^{t} \int_{\Omega} h\left(\frac{\partial v^q}{\partial t}\right) dx ds + \int_{t_0}^{t} \int_{\Omega} \frac{f(x, v)}{v^{q-1}} \frac{\partial v^q}{\partial t} dx ds + \frac{q}{p} \|v(t_0)\|_{W_0^{s,p}(\Omega)}^p. \end{aligned} \quad (3.31) $$

Taking lim sup in (3.31) as $t \to t_0^+$ and by (3.29), we obtain

$$ \|v(t_0)\|_{W_0^{s,p}(\Omega)} = \lim_{t \to t_0^+} \|v(t)\|_{W_0^{s,p}(\Omega)} $$

and hence the right-continuity of $v : [0, T] \to W_0^{s,p}(\Omega)$ follows.

Now, for proving the left continuity, consider $0 < \eta \le t - t_0$, multiply (1.12) by

$$ \tau_{\eta} v = \frac{v^q(\cdot + \eta, \cdot) - v^q(\cdot, \cdot)}{\eta v^{q-1}} \in L^2(Q_T) \cap L^1(0, T; W_0^{s,p}(\Omega)) $$

and integrate over $(t_0, t) \times \Omega$. Using Proposition 1.2 and Young's inequality again, we obtain

$$ \begin{aligned} & \int_{t_0}^{t} \int_{\Omega} v^{q-1} \partial_t (v^q) \tau_{\eta} v \, dx \, ds + \frac{q}{p\eta} \int_{t_0}^{t} \left( \|v(s+\eta)\|_{W_0^{s,p}(\Omega)}^p - \|v(s)\|_{W_0^{s,p}(\Omega)}^p \right) ds \\ & \ge \int_{t_0}^{t} \int_{\Omega} h v^{q-1} \tau_{\eta} v \, dx \, ds + \int_{t_0}^{t} \int_{\Omega} f(x, v) \tau_{\eta} v \, dx \, ds. \end{aligned} $$
It follows that

$$ 
\begin{aligned}
& \int_{t_0}^{t} \int_{\Omega} v^{q-1} \partial_t (v^q) \tau_{\eta} v \, dx \, ds \\
& + \frac{q}{p\eta} \left( \int_{t}^{t+\eta} \|v(s)\|_{W_0^{s,p}(\Omega)}^p ds - \int_{t_0}^{t_0+\eta} \|v(s)\|_{W_0^{s,p}(\Omega)}^p ds \right) && (3.32) \\
\geq & \int_{t_0}^{t} \int_{\Omega} h v^{q-1} \tau_{\eta} v \, dx \, ds + \int_{t_0}^{t} \int_{\Omega} f(x, v) \tau_{\eta} v \, dx \, ds.
\end{aligned}
$$

By the right continuity of $v : [0, T] \to W_0^{s,p}(\Omega)$ and by dominated convergence Theorem, as $\eta \to 0^+$ we have

$$ 
\begin{aligned}
\frac{q}{p\eta} \int_{t}^{t+\eta} \|v(s)\|_{W_0^{s,p}(\Omega)}^p ds &\to \frac{q}{p} \|v(t)\|_{W_0^{s,p}(\Omega)}^p, \\
\frac{q}{p\eta} \int_{t_0}^{t_0+\eta} \|v(s)\|_{W_0^{s,p}(\Omega)}^p ds &\to \frac{q}{p} \|v(t_0)\|_{W_0^{s,p}(\Omega)}^p.
\end{aligned}
$$

Hence as $\eta \to 0^+$, (3.32) yields

$$ 
\begin{aligned}
& \int_{t_0}^{t} \int_{\Omega} \left(\frac{\partial v^q}{\partial t}\right)^2 dx \, ds + \frac{q}{p} \|v(t)\|_{W_0^{s,p}(\Omega)}^p \\
\geq & \int_{t_0}^{t} \int_{\Omega} h\left(\frac{\partial v^q}{\partial t}\right) dx \, ds + \int_{t_0}^{t} \int_{\Omega} \frac{f(x, v)}{v^{q-1}} \frac{\partial v^q}{\partial t} dx \, ds + \frac{q}{p} \|v(t_0)\|_{W_0^{s,p}(\Omega)}^p.
\end{aligned}
$$

From the above inequality, we deduce that the equality in (3.31) holds and the left-continuity of $v : [0, T] \to W_0^{s,p}(\Omega)$ follows. $\square$

## 4. STABILIZATION

### 4.1. Existence and uniqueness of the solution of the stationary problem related to (1.12).
In this subsection we deal with the stationary problem in order to determine the asymptotic behavior of trajectories to (1.1). Precisely, we consider the problem

$$ 
\begin{aligned}
(-\Delta)_p^s v &= b(x)v^{q-1} + f(x, v) \quad \text{in } \Omega, \\
v &> 0 \quad \text{in } \Omega, \\
v &= 0 \quad \text{in } \mathbb{R}^N \setminus \Omega
\end{aligned}
\quad (4.1) $$

where $b \in L^\infty(\Omega)$ and nonnegative. We define the notion of a weak solution as follows.

**Definition 4.1.** A positive function $v \in W_0^{s,p}(\Omega) \cap L^\infty(\Omega)$ is called a weak solution to problem (4.1) if

$$ 
\begin{aligned}
& \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|v(x) - v(y)|^{p-2} (v(x) - v(y)) (\varphi(x) - \varphi(y))}{|x - y|^{N+sp}} dx dy \\
& = \int_{\Omega} (b(x)v^{q-1} + f(x, v))\varphi dx
\end{aligned}
\quad (4.2) $$

for any $\varphi \in W_0^{s,p}(\Omega)$.

**Theorem 4.2.** Assume that $f$ satisfies (H1)–(H3). Let $q \in (1, p]$. In addition if $q = p$ suppose that $\|b\|_\infty < \lambda_{1,p,s}$. Then, there exists a unique weak solution $v \in C(\bar{\Omega}) \cap M_{ds}^1(\Omega)$ to problem (4.1).
**Proof.** By following the same arguments as in Theorem 2.2, we deduce the existence of a nonnegative global minimizer to the following energy functional

$$ \mathcal{L}(v) = \frac{1}{p} \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|v(x) - v(y)|^p}{|x - y|^{N+ps}} dx dy - \frac{1}{q} \int_{\Omega} b(v^+)^q dx - \int_{\Omega} F(x, v)dx $$

where $F$ is defined in (2.4). Then, as in the proof of Theorem 2.2 step 2, we infer that $v \in L^\infty(\Omega)$. Furthermore, by using [27, Theorem 1.1], there is $\alpha \in (0, s]$ such that $v \in C^{0,\alpha}(\bar{\Omega})$. Next, by [20, Theorems 1.4 and 1.5, p. 768], we obtain that $v > 0$ in $\Omega$ and satisfies $v \ge k d^s(x)$ for some $k > 0$. Finally, [[27, Theorem 4.4]] implies that $v \in \mathcal{M}_{d^s}^1(\Omega)$.

Let $v_1, v_2 \in C(\bar{\Omega}) \cap \mathcal{M}_{d^s}^1(\Omega)$ be two solutions of (4.1), we choose $\frac{v_1^q - v_2^q}{v_1^{q-1}}$ and $\frac{v_2^q - v_1^q}{v_2^{q-1}}$ as test functions in (4.1) satisfied by $v_1, v_2$ respectively. Then adding the equations, we deduce from Lemma 1.8 and (H2),

$$ \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \left( \frac{|v_1(x) - v_1(y)|^{p-2}(v_1(x) - v_1(y))}{|x - y|^{N+sp}} \left[ \left( \frac{v_1^q - v_2^q}{v_1^{q-1}} \right)(x) - \left( \frac{v_1^q - v_2^q}{v_1^q} \right)(y) \right] + \frac{|v_2(x) - v_2(y)|^{p-2}(v_2(x) - v_2(y))}{|x - y|^{N+sp}} \left[ \left( \frac{v_2^q - v_1^q}{v_2^{q-1}} \right)(x) - \left( \frac{v_2^q - v_1^q}{v_2^q} \right)(y) \right] \right) dx dy = 0. $$

Again by Lemma 1.8, for $1 < q < p$, we obtain $v_1 \equiv v_2$ in $\mathbb{R}^N$. While for $q = p$, we have $v_1(x) = k v_2(x)$ a.e. in $\mathbb{R}^N$, for some $k > 0$. Without loss of generality, we can assume that $k \le 1$. Then from (H2) we obtain

$$ \begin{aligned} (-\Delta)_p^s (kv_2) &= k^{p-1} (-\Delta)_p^s (v_2) = k^{p-1} (b(x)v_2^{p-1} + f(x, v_2)) \\ &< b(x)(kv_2)^{p-1} + f(x, kv_2) \\ &= (-\Delta)_p^s (v_1) \end{aligned} $$

which yields a contradiction. Hence $k = 1$ and $v_1 \equiv v_2$. $\square$

<figure><img src="image_2.png" /></figure>

Next, as in the proof of Corollary 2.4, we obtain the following result.

**Corollary 4.3.** Under the conditions of Theorem 4.2, there exists one and only one weak solution $u \in \dot{V}_+^q \cap \mathcal{M}_{d^s}^{1/q}(\Omega)$ to the problem

$$ \begin{aligned} & \mathcal{T}_q u = b \quad \text{in } \Omega, \\ & u > 0 \quad \text{in } \Omega, \\ & u = 0 \quad \text{in } \mathbb{R}^N \setminus \Omega. \end{aligned} \tag{4.3} $$

Furthermore,

$$ \int_{\mathbb{R}^n} \int_{\mathbb{R}^n} \frac{|u^{1/q}(x) - u^{1/q}(y)|^{p-2} (u^{1/q}(x) - u^{1/q}(y)) \left( (u^{\frac{1-q}{q}} \Psi)(x) - (u^{\frac{1-q}{q}} \Psi)(y) \right)}{|x-y|^{N+sp}} dx dy \\ = \int_{\Omega} b \Psi dx + \int_{\Omega} f(x, u^{1/q}) u^{\frac{1-q}{q}} \Psi dx $$

for all $\Psi$ satisfies (2.13).

**Proof of Theorem 1.14.** We are ready now to prove our stabilization result by using the same approach as in the proof of [25, Theorem 3.10].

**Proof of Theorem 1.14.** We consider two cases.
**Case 1:** $h = h_\infty$. We introduce the family of operators $\{S(t) : t \ge 0\}$ defined on $\dot{V}_+^q \cap \mathcal{M}_{ds}^{1/q}(\Omega)$ as $w(t) = S(t)w_0$ where $w$ is the unique solution (obtained in Theorem 3.3) to

$$ 
\begin{array}{l}
\partial_t w + \mathcal{T}_q w = h_\infty \quad \text{in } Q_T, \\
w > 0 \quad \text{in } Q_T, \\
w = 0 \quad \text{on } \Gamma, \\
w(0, \cdot) = w_0 \quad \text{in } \Omega.
\end{array} 
\quad (4.4) $$

Thus, we claim that $\{S(t) : t \ge 0\}$ defines a semi-group of contractions in $L^2(\Omega)$. Indeed, from the uniqueness and above properties of solutions to problem (4.4) we infer that for any $w_0 \in \dot{V}_+^q \cap \mathcal{M}_{ds}^{1/q}(\Omega)$,

$$ 
\begin{array}{l}
S(t+s)w_0 = S(t)S(s)w_0, \\
S(0)w_0 = w_0.
\end{array} 
\quad (4.5) $$

From (3.17) and (3.19) the map $[0, \infty) \ni t \mapsto S(t)w_0$ is continuous and T-accretive $L^2(\Omega)$. Note that $\tilde{v} = (S(t)w_0)^{1/q}$ is the solution of (1.12) with $h = h_\infty$ and the initial data $w_0^{1/q}$.

Let us denote $v$ the solution of (1.12) with $h = h_\infty$ and the initial data $v_0$. Hence we obtain $u(t) = v(t)^q = S(t)u_0$ with $u_0 = v_0^q$. Let $\underline{w} = w_\mu$ be the solution of (3.4) and $\overline{w} = w_K$ or the solution to (3.7) if $q = p$. Then, $\underline{w}, \overline{w} \in \mathcal{M}_{ds}^{1/q}(\Omega)$ and for $\mu$ small enough and $K$ large enough, $w$ is a subsolution and $\overline{w}$ a supersolution to (4.1) with $b = h_\infty$ such that $\underline{w} \le v_0 \le \overline{w}$. We then define $\underline{u}(t) = S(t)\underline{w}^q$ and $\overline{u}(t) = S(t)\overline{w}^q$ the solutions to (4.4). Therefore, $\underline{u} := (\underline{v})^q$ and $\overline{u} := (\overline{v})^q$ are obtained by the iterative scheme (3.2) with $v_0 = \underline{w}$ and $v_0 = \overline{w}$. Hence, by comparison principle the maps $t \mapsto \underline{u}(t)$ and $t \mapsto \overline{u}(t)$ are respectively nondecreasing and non-increasing. In the other hand, (1.13) ensures that for any $t \ge 0$,

$$ \underline{w} \le \underline{u}(t) \le u(t) \le \overline{u}(t) \le \overline{w}. \quad (4.6) $$

We set $\underline{u}_\infty = \lim_{t \to \infty} \underline{u}(t)$ and $\overline{u}_\infty = \lim_{t \to \infty} \overline{u}(t)$. Then from (4.5), we obtain

$$ 
\begin{aligned}
\underline{u}_\infty &= \lim_{s \to \infty} S(t+s)\underline{w}^q = S(t) \lim_{s \to \infty} (S(s)(\underline{w}^q)) = S(t)\underline{u}_\infty, \\
\overline{u}_\infty &= \lim_{s \to \infty} S(t+s)\overline{w}^q = S(t) \lim_{s \to \infty} (S(s)(\overline{w}^q)) = S(t)\overline{u}_\infty.
\end{aligned} $$

This implies that $\underline{u}_\infty$ and $\overline{u}_\infty$ are the stationary solutions to (4.3) with $b = h_\infty$. By uniqueness, we have $u_{\text{stat}} := \underline{u}_\infty = \overline{u}_\infty$ where $u_{\text{stat}}$ is the stationary solution to (4.4). Therefore from (4.6) and by dominated convergence Theorem, we obtain

$$ \|u(t) - u_{\text{stat}}\|_{L^2(\Omega)} \to 0 \quad \text{as} \quad t \to \infty. $$

Thus using (4.6) and the interpolation inequality with $2 < r < \infty$,

$$ \| \cdot \|_r \le \| \cdot \|_{\infty}^{\alpha} \| \cdot \|_2^{1-\alpha}, $$

we obtain, the above convergence for any $r \ge 1$.

**Case 2:** $h \ne h_\infty$. From (1.16), for any $\epsilon > 0$ there exists $t_0 > 0$ large enough such that $\int_{t_0}^{+\infty} \frac{1}{l(t)}dt < \epsilon$ and for any $t \ge t_0$,

$$ l(t)\|h(t, \cdot) - h_\infty\|_{L^2(\Omega)} \le M \text{ for some } M > 0. $$

Let $T > 0$ and $v$ be the solution of the problem (1.12) obtained by Theorem 3.1 with $h$ and the initial data $v_0 = u_0^{1/q}$ and set $u = v^q$. Since $v$ satisfies (3.1), we can
define $\tilde{u}(t) = S(t + t_0)u_0 = S(t)u(t_0)$. Then, by (1.13) and uniqueness argument, we have for any $t > 0$,

$$ \begin{aligned} \|u(t + t_0, \cdot) - \tilde{u}(t, \cdot)\|_{L^2(\Omega)} &\le \int_0^t \|h(s + t_0, \cdot) - h_\infty\|_{L^2(\Omega)} ds \\ &\le M \int_{t_0}^{+\infty} \frac{1}{l(s)} ds \le M\epsilon. \end{aligned} $$

By Case 1, we have $\tilde{u}(t) \to u_{\text{stat}}$ in $L^2(\Omega)$ as $t \to \infty$. Therefore, we obtain

$$ \|u(t) - u_{\text{stat}}\|_{L^2(\Omega)} \to 0 \quad \text{as } t \to \infty. $$

Using again the interpolation inequality above, we conclude the proof of Theorem 1.14. □

## 5. Appendix

### 5.1. Regularity results.
The first one is obtained by a similar proof as in [23] (see also [25]).

**Proposition 5.1.** Let $u \in W_0^{s,p}(\Omega)$ satisfying

$$ \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|u(x) - u(y)|^{p-2}(u(x) - u(y))(\Psi(x) - \Psi(y))}{|x - y|^{N+sp}} dx dy = \int_{\Omega} f(x, u) \Psi dx \quad (5.1) $$

for all $\Psi \in W_0^{s,p}(\Omega)$, where $f$ satisfies for all $(x, t) \in \Omega \times \mathbb{R}$,

$$ |f(x,t)| \le C (1 + |t|^{r-1}), \quad \forall x \in \bar{\Omega}, \quad 1 < r \le p. $$

Then $u \in L^\infty(\Omega)$.

**Proposition 5.2.** Let $1 < q \le p$. Assume that $u \in \mathbf{W}$ and nonnegative satisfying for any $\Psi \in \mathbf{W}$,

$$ \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|u(x) - u(y)|^{p-2}(u(x) - u(y))(\Psi(x) - \Psi(y))}{|x - y|^{N+sp}} dx dy = \int_{\Omega} h u^{q-1} \Psi dx \quad (5.2) $$

where $h \in L^2(\Omega) \cap L^r(\Omega)$ with $r > \max\{1, \frac{N}{sp}\}$ and $h \ge 0$ a.e. in $\Omega$. Then $u \in L^\infty(\Omega)$.

**Proof.** We follow the main steps in the proof of [8, Theorem 3.1]. For every $\delta > 0$, we define $u_\delta = u + \delta$. Given $\beta \ge 1$, we insert the test function $\psi = u_\delta^\beta - \delta^\beta$ in (5.2), then we obtain

$$ \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|u(x) - u(x)|^{p-2}(u(x) - u(x))(u_\delta(x)^\beta - u_\delta(y)^\beta)}{|x - y|^{N+sp}} dx dy \le \int_{\Omega} h u^{q-1} u_\delta^\beta dx. $$

By using the inequality in [8, Lemma A.2], we obtain

$$ \frac{\beta p^p}{(\beta + p - 1)^p} \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|u_\delta(x)^{\frac{\beta+p-1}{p}} - u_\delta(y)^{\frac{\beta+p-1}{p}}|^p}{|x-y|^{N+sp}} dx dy \le \int_{\Omega} h u^{q-1} u_\delta^\beta dx. $$

By Theorem 1.3, we obtain

$$ \begin{aligned} & \left( \int_{\Omega} \left( u_{\delta}(x)^{\frac{\beta+p-1}{p}} - \delta^{\frac{\beta+p-1}{p}} \right) p_s^* dx \right)^{p/p_s^*} \\ & \le C_{N,s,p} \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|u_{\delta}(x)^{\frac{\beta+p-1}{p}} - u_{\delta}(y)^{\frac{\beta+p-1}{p}}|^p}{|x-y|^{N+sp}} dx dy. \end{aligned} $$
By the triangle inequality, the left-hand side of the bovd inequality, can be estimated as

$$ \left( \int_{\Omega} \left( u_{\delta}^{\frac{\beta+p-1}{p}} \right)^{p_s^*} dx \right)^{p/p_s^*} \le \left( \int_{\Omega} \left( u_{\delta}^{\frac{\beta+p-1}{p}} - \delta^{\frac{\beta+p-1}{p}} \right)^{p_s^*} dx \right)^{p/p_s^*} + \delta^{\beta+p-1} |\Omega|^{p/p_s^*}. $$

On the other hand, we use the inequality $u_{\delta}^{\beta+p-1} \ge \delta^{p-q}u_{\delta}^{\beta+q-1}$, Hölder and interpolation inequalities, for $r > \frac{N}{sp}$ and with the observation that $p < pr' < p_s^*$, where $r' = \frac{r}{r-1}$ to obtain

$$ \begin{aligned} \int_{\Omega} h u^{q-1} u_{\delta}^{\beta} dx &\le \delta^{q-p} \int_{\Omega} h u_{\delta}^{p+\beta-1} dx \le \delta^{q-p} \|h\|_{L^r} \left( \int_{\Omega} \left( u_{\delta}^{\frac{p+\beta-1}{p}} \right)^{pr'} dx \right)^{1/r'} \\ &\le \delta^{q-p} \|h\|_{L^r} \left( \int_{\Omega} \left( u_{\delta}^{\frac{p+\beta-1}{p}} \right)^{p_s^*} dx \right)^{\frac{p}{p_s^*}} \left( \int_{\Omega} u_{\delta}^{p+\beta-1} dx \right)^{1-\alpha} \end{aligned} \quad (5.3) $$

where $\frac{1}{pr'} = \frac{\alpha}{p} + \frac{1-\alpha}{p_s^*}$ and $0 \le \alpha \le 1$. Using Young's inequality,

$$ \int_{\Omega} h u^{q-1} u_{\delta}^{\beta} dx \le \delta^{q-p} \|h\|_{L^r} \left[ \epsilon \left( \int_{\Omega} \left( u_{\delta}^{\frac{p+\beta-1}{p}} \right)^{p_s^*} dx \right)^{p/p_s^*} + C_{\epsilon} \int_{\Omega} u_{\delta}^{p+\beta-1} dx \right] $$

with $C_{\epsilon} = \epsilon^{-\frac{1}{\alpha-1}}$, it is easy to see that

$$ \delta^{p+\beta-1} |\Omega|^{p/p_s^*} \le \frac{1}{\beta} \left( \frac{p+\beta-1}{p} \right)^p |\Omega|^{\frac{p}{p_s^*}-1} \int_{\Omega} u_{\delta}^{p+\beta-1} dx. $$

Taking

$$ \epsilon = \frac{\beta \delta^{p-q}}{2C_{N,s,p} \|h\|_{L^r}} \left( \frac{p}{p+\beta-1} \right)^p > 0, $$

we obtain

$$ \begin{aligned} & \left( \int_{\Omega} \left( u_{\delta}^{\frac{\beta+p-1}{p}} \right)^{p_s^*} dx \right)^{p/p_s^*} \\ &\le \frac{C_{N,s,p}}{\beta} \left( \frac{p+\beta-1}{p} \right)^p \left[ \delta^{q-p} \|h\|_{L^r} C_{\epsilon} + |\Omega|^{\frac{p}{p_s^*}-1} \right] \int_{\Omega} u_{\delta}^{p+\beta-1} dx. \end{aligned} $$

We then choose

$$ \delta = (C_{\epsilon} \|h\|_{L^r})^{\frac{1}{q-p}} |\Omega|^{\frac{1}{q-p} (\frac{p}{p_s^*} - 1)} > 0 $$

and set $v = \beta + p - 1$. Then the previous inequality can be written as

$$ \left( \int_{\Omega} u_{\delta}^{\left(\frac{p_s^*}{p}\right)v} dx \right)^{\frac{1}{\left(\frac{p_s^*}{p}\right)v}} \le \left[ C \left|\Omega\right|^{\frac{p}{p_s^*}-1} \right]^{1/v} \left( \frac{v}{p} \right)^{p/v} \left( \int_{\Omega} u_{\delta}^{v} dx \right)^{1/v} $$

with $C = C(N, s, p) > 0$. We now iterate the previous inequality, by taking the sequence of exponents

$$ v_0 = 1 \quad \text{and} \quad v_{n+1} = \left(\frac{p_s^*}{p}\right) v_n = \left(\frac{p_s^*}{p}\right)^{n+1}. $$

We have

$$ \begin{aligned} & \sum_{n=0}^{\infty} \frac{1}{v_n} = \sum_{n=0}^{\infty} \left(\frac{p}{p_s^*}\right)^n = \frac{p_s^*}{p_s^* - p}, \\ & \prod_{n=0}^{\infty} \left(\frac{v_n}{p}\right)^{\frac{p}{v_n}} < \infty. \end{aligned} $$
By starting from 0 at the step $n$ we have

$$ \|u_{\delta}\|_{L^{\nu_{n+1}}(\Omega)} \leq [C |\Omega|^{\frac{p}{p_s^*}-1}]^{\sum_{i=0}^{n} \frac{1}{v_i}} \prod_{i=0}^{n} \left(\frac{v_i}{p}\right)^{\frac{p}{v_i}} \|u_{\delta}\|_{L^1(\Omega)}. $$

By taking the limit as $n$ approaches $\infty$, we finally obtain

$$ \|u_{\delta}\|_{L^{\infty}(\Omega)} \leq \frac{C'}{|\Omega|} \|u_{\delta}\|_{L^1(\Omega)} \leq \frac{C'}{|\Omega|} (\|u\|_{L^1(\Omega)} + \delta|\Omega|) $$

for some constant $C' = C'(N, p, s) > 0$. $\square$

Combining Proposition 5.1 with Proposition 5.2, we have the following corollary:

**Corollary 5.3.** Let $1 < q \leq p$. Assume $u \in W$, nonnegative and satisfying for any nonnegative $\Psi \in W$

$$ \int_{\Omega} u^{2q-1} \Psi dx + \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|u(x) - u(y)|^{p-2}(u(x) - u(y))(\Psi(x) - \Psi(y))}{|x-y|^{N+sp}} dx dy \\ \leq \int_{\Omega} (f(x, u) + h u^{q-1}) \Psi dx $$

where $f$ satisfies for all $(x, t) \in \Omega \times \mathbb{R}$, $|f(x, t)| \leq C(1+t^{q-1})$ and $h \in L^2(\Omega) \cap L^r(\Omega)$ with $r > \max\{1, \frac{N}{sp}\}$. Then $u \in L^\infty(\Omega)$.

5.2. **Comparison principle.** Following the proof of [2, Theorem 4.3] and using Lemma 1.8, we have the following new comparison principle.

**Theorem 5.4.** Assume $f$ satisfies (H1), (H2). Let $\underline{v}, \bar{v} \in W \cap L^\infty(\Omega)$ be nonnegative functions respectively subsolution and supersolution to (2.1) for some $h_0 \in (L^r(\Omega))^+$ with $r \geq 2$. Then $\underline{v} \leq \bar{v}$.

*Proof.* For any nonnegative pair $\Phi, \Psi \in W$ we have

$$ \int_{\Omega} \underline{v}^{2q-1} \Phi dx + \lambda \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|\underline{v}(x) - \underline{v}(y)|^{p-2} (\underline{v}(x) - \underline{v}(y))(\Phi(x) - \Phi(y))}{|x-y|^{N+sp}} dx dy \\ \leq \int_{\Omega} h_0 \underline{v}^{q-1} \Phi dx + \lambda \int_{\Omega} f(x, \underline{v}) \Phi dx $$

and

$$ \int_{\Omega} \bar{v}^{2q-1} \Psi dx + \lambda \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} \frac{|\bar{v}(x) - \bar{v}(y)|^{p-2} (\bar{v}(x) - \bar{v}(y))(\Psi(x) - \Psi(y))}{|x-y|^{N+sp}} dx dy \\ \geq \int_{\Omega} h_0 \bar{v}^{q-1} \Psi dx + \lambda \int_{\Omega} f(x, \bar{v}) \Psi dx. $$

Subtracting the above inequalities with test functions

$$ \Phi = \left( \frac{(\underline{v} + \epsilon)^q - (\bar{v} + \epsilon)^q}{(\underline{v} + \epsilon)^{q-1}} \right)^{+}, \quad \Psi = \left( \frac{(\bar{v} + \epsilon)^q - (\underline{v} + \epsilon)^q}{(\bar{v} + \epsilon)^{q-1}} \right)^{-} \in W, $$
with $\epsilon \in (0, 1)$, we obtain

$$ 
\begin{aligned}
& \int_{\{\underline{v} > \bar{v}\}} \left( \frac{\underline{v}^{2q-1}}{(\underline{v} + \epsilon)^{q-1}} - \frac{\bar{v}^{2q-1}}{(\bar{v} + \epsilon)^{q-1}} \right) ((\underline{v} + \epsilon)^q - (\bar{v} + \epsilon)^q) dx \\
& + \lambda \int_{\{\underline{v} > \bar{v}\}} \int_{\{\bar{v} > \underline{v}\}} \frac{|\underline{v}(x) - \underline{v}(y)|^{p-2} (\underline{v}(x) - \underline{v}(y))}{|x - y|^{N+sp}} \\
& \times \left[ \frac{(\underline{v}(x) + \epsilon)^q - (\bar{v}(x) + \epsilon)^q}{(\underline{v}(x) + \epsilon)^{q-1}} - \frac{(\underline{v}(y) + \epsilon)^q - (\bar{v}(y) + \epsilon)^q}{(\underline{v}(y) + \epsilon)^{q-1}} \right] dx dy \\
& + \lambda \int_{\{\underline{v} > \bar{v}\}} \int_{\{\bar{v} > \underline{v}\}} \frac{|\bar{v}(x) - \bar{v}(y)|^{p-2} (\bar{v}(x) - \bar{v}(y))}{|x - y|^{N+sp}} \\
& \times \left[ \frac{(\bar{v}(x) + \epsilon)^q - (\underline{v}(x) + \epsilon)^q}{(\bar{v}(x) + \epsilon)^{q-1}} - \frac{(\bar{v}(y) + \epsilon)^q - (\underline{v}(y) + \epsilon)^q}{(\bar{v}(y) + \epsilon)^{q-1}} \right] dx dy \\
& \leq \int_{\{\underline{v} > \bar{v}\}} h_0 \left( \frac{\underline{v}^{q-1}}{(\underline{v} + \epsilon)^{q-1}} - \frac{\bar{v}^{q-1}}{(\bar{v} + \epsilon)^{q-1}} \right) ((\underline{v} + \epsilon)^q - (\bar{v} + \epsilon)^q) dx \\
& + \lambda \int_{\{\underline{v} > \bar{v}\}} \left( \frac{f(x, \underline{v})}{(\underline{v} + \epsilon)^{q-1}} - \frac{f(x, \bar{v})}{(\bar{v} + \epsilon)^{q-1}} \right) ((\underline{v} + \epsilon)^q - (\bar{v} + \epsilon)^q) dx.
\end{aligned} 
\quad (5.4) $$

Since $\frac{\underline{v}}{\underline{v}+\epsilon} \leq \frac{\bar{v}}{\bar{v}+\epsilon} < 1$ in $\{\underline{v} > \bar{v}\}$, we obtain

$$ 
\begin{aligned}
& \left( \frac{\underline{v}^{2q-1}}{(\underline{v} + \epsilon)^{q-1}} - \frac{\bar{v}^{2q-1}}{(\bar{v} + \epsilon)^{q-1}} \right) ((\underline{v} + \epsilon)^q - (\bar{v} + \epsilon)^q) \\
& = \left( \underline{v}^q \left( \frac{\underline{v}}{(\underline{v} + \epsilon)} \right)^{q-1} - \bar{v}^q \left( \frac{\bar{v}}{(\bar{v} + \epsilon)} \right)^{q-1} \right) ((\underline{v} + \epsilon)^q - (\bar{v} + \epsilon)^q) \\
& \leq \underline{v}^q ((\underline{v} + \epsilon)^q - (\bar{v} + \epsilon)^q) \leq \underline{v}^q (\bar{v} + 1)^q.
\end{aligned} $$

In the same spirit, we infer that

$$ 0 \leq h_0 \left( \frac{\underline{v}^{q-1}}{(\underline{v} + \epsilon)^{q-1}} - \frac{\bar{v}^{q-1}}{(\bar{v} + \epsilon)^{q-1}} \right) ((\underline{v} + \epsilon)^q - (\bar{v} + \epsilon)^q) \leq h_0 (\bar{v} + 1)^q. $$

Moreover, as $\epsilon \to 0$, we have

$$ 
\begin{aligned}
& \left( \frac{\underline{v}^{2q-1}}{(\underline{v} + \epsilon)^{q-1}} - \frac{\bar{v}^{2q-1}}{(\bar{v} + \epsilon)^{q-1}} \right) ((\underline{v} + \epsilon)^q - (\bar{v} + \epsilon)^q) \to (\underline{v}^q - \bar{v}^q)^2, \\
& h_0 \left( \frac{\underline{v}^{q-1}}{(\underline{v} + \epsilon)^{q-1}} - \frac{\bar{v}^{q-1}}{(\bar{v} + \epsilon)^{q-1}} \right) ((\underline{v} + \epsilon)^q - (\bar{v} + \epsilon)^q) \to 0.
\end{aligned} $$

Then, by the dominated convergence Theorem, we have

$$ 
\begin{aligned}
& \lim_{\epsilon \to 0} \int_{\{\underline{v} > \bar{v}\}} \left( \frac{\underline{v}^{2q-1}}{(\underline{v} + \epsilon)^{q-1}} - \frac{\bar{v}^{2q-1}}{(\bar{v} + \epsilon)^{q-1}} \right) ((\underline{v} + \epsilon)^q - (\bar{v} + \epsilon)^q) dx \\
& = \int_{\{\underline{v} > \bar{v}\}} (\underline{v}^q - \bar{v}^q)^2 dx
\end{aligned} 
\quad (5.5) $$

and

$$ \lim_{\epsilon \to 0} \int_{\{\underline{v} > \bar{v}\}} h_0 \left( \frac{\underline{v}^{q-1}}{(\underline{v} + \epsilon)^{q-1}} - \frac{\bar{v}^{q-1}}{(\bar{v} + \epsilon)^{q-1}} \right) ((\underline{v} + \epsilon)^q - (\bar{v} + \epsilon)^q) dx = 0. \quad (5.6) $$

Then by Fatou's Lemma and (H1), we have

$$ - \liminf_{\epsilon \to 0} \int_{\{\underline{v} > \bar{v}\}} \frac{f(x, \underline{v})}{(\underline{v} + \epsilon)^{q-1}} (\bar{v} + \epsilon)^q dx \leq - \int_{\{\underline{v} > \bar{v}\}} \frac{f(x, \underline{v})}{\underline{v}^{q-1}} \bar{v}^q dx, \quad (5.7) $$
$$ -\liminf_{\epsilon \to 0} \int_{\{\bar{v} > v\}} \frac{f(x, \bar{v})}{(\bar{v} + \epsilon)^{q-1}} (\underline{v} + \epsilon)^q dx \le -\int_{\{\bar{v} > v\}} \frac{f(x, \bar{v})}{\bar{v}^{q-1}} \underline{v}^q dx \quad (5.8) $$

$$ \lim_{\epsilon \to 0} \int_{\{\bar{v} > v\}} f(x, \underline{v})(\underline{v} + \epsilon) dx = \int_{\{\bar{v} > v\}} f(x, \underline{v}) \underline{v} dx, \quad (5.9) $$

$$ \lim_{\epsilon \to 0} \int_{\{\bar{v} > v\}} f(x, \bar{v})(\bar{v} + \epsilon) dx = \int_{\{\bar{v} > v\}} f(x, \bar{v}) \bar{v} dx. \quad (5.10) $$

From Lemma 1.8, we have

$$ \begin{aligned} & \int_{\{\underline{v} > \bar{v}\}} \int_{\{\underline{v} > \bar{v}\}} \frac{|\underline{v}(x) - \underline{v}(y)|^{p-2} (\underline{v}(x) - \underline{v}(y))}{|x-y|^{N+sp}} \\ & \times \left[ \frac{(\underline{v}(x) + \epsilon)^q - (\bar{v}(x) + \epsilon)^q}{(\underline{v}(x) + \epsilon)^{q-1}} - \frac{(\underline{v}(y) + \epsilon)^q - (\bar{v}(y) + \epsilon)^q}{(\underline{v}(y) + \epsilon)^{q-1}} \right] dx dy \\ & + \int_{\{\underline{v} > \bar{v}\}} \int_{\{\underline{v} > \bar{v}\}} \frac{|\bar{v}(x) - \bar{v}(y)|^{p-2} (\bar{v}(x) - \bar{v}(y))}{|x-y|^{N+sp}} \\ & \times \left[ \frac{(\bar{v}(x) + \epsilon)^q - (\underline{v}(x) + \epsilon)^q}{(\bar{v}(x) + \epsilon)^{q-1}} - \frac{(\bar{v}(y) + \epsilon)^q - (\underline{v}(y) + \epsilon)^q}{(\bar{v}(y) + \epsilon)^{q-1}} \right] dx dy \ge 0. \end{aligned} \quad (5.11) $$

Therefore, plugging (5.5)-(5.11) and taking lim sup in (5.4), we obtain from (H2),
$\epsilon \to 0$

$$ 0 \le \int_{\{\underline{v} > \bar{v}\}} (\underline{v}^q - \bar{v}^q)^2 dx \le \lambda \int_{\{\underline{v} > \bar{v}\}} \left( \frac{f(x, \underline{v})}{\underline{v}^{q-1}} - \frac{f(x, \bar{v})}{\bar{v}^{q-1}} \right) (\underline{v}^q - \bar{v}^q) dx \le 0 $$

from which $\underline{v} \le \bar{v}$ follows. □

## REFERENCES

[1] B. Abdellaoui, A. Attar, R. Bentifour, I. Peral; On fractional p-Laplacian parabolic problem with general data, *Ann. Mat. Pura Appl.*, (4) **197** (2018), 329-356.

[2] R. Arora, J. Giacomoni, G. Warnault; A Picone identity for variable exponent operators and applications, *Adv. Nonlinear Anal.*, **9** (2020), no. 1, 327-360.

[3] R. Arora, J. Giacomoni, G. Warnault; Doubly nonlinear equation involving p(x)-homogeneous operators: local existence, uniqueness and global behaviour, *J. Math. Anal. Appl.*, **487** (2020), no. 2, 27 p.

[4] M. Badra, K. Bal, J. Giacomoni; A singular parabolic equation: Existence, stabilization, *J. Differential Equations*, **252** (2012), 5042-5075.

[5] B. Barrios, A. Figalli, X. Ros-Oton; Free boundary regularity in the parabolic fractional obstacle problem, *Comm. Pure Appl. Math.*, **71** (2018), no. 10, 2129-2159.

[6] Benilan P., Picard C; *Quelques aspects non lineaires du principe du maximum*. In: Hirsch F., Mokobodzki G. (eds) Séminaire de Théorie du Potentiel Paris, No. 4. Lecture Notes in Mathematics, vol 713. Springer, Berlin, Heidelberg, 1979, 1-37.

[7] F. Boyer, P. Fabrie; *Mathematical tools for the study of the incompressible Navier-Stokes equations and related models*, volume 183 of Applied Mathematical Sciences. Springer, New York, 2013.

[8] L. Brasco, E. Parini; The second eigenvalue of the fractional p-Laplacian, *Adv. Calc. Var.*, **9** (2016), no. 4, 323-355.

[9] L. Brasco, G. Franzina; Convexity properties of Dirichlet integrals and Picone-type inequalities, *Kodai Math. J.*, **37** (2014), 769-799.

[10] H. Brezis; *Opérateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert*, Elsevier, 1973.

[11] X. Cabré, J. Tan; Positive solutions of nonlinear problems involving the square root of the Laplacian, *Adv. Math.*, **224** (2010), 2052-2093.

[12] X. Cabré, Y. Sire; Nonlinear equations for fractional Laplacians, I: Regularity, maximum principles, and Hamiltonian estimates, *Ann. Inst. H. Poincaré Anal. Non Linéaire*, **31** (2014), 23-53.
[13] L. Caffarelli, J. L. Vázquez; Nonlinear porous medium flow with fractional potential pressure, *Arch. Ration. Mech. Anal.*, **202** (2011), 537–565.

[14] L. Caffarelli, L. Silvestre; An extension problem related to the fractional Laplacian, *Comm. Partial Differential Equations*, **32** (2007), no. 8, 1245-1260.

[15] L. Caffarelli; Non-local diffusions, drifts and games, *Nonlinear Partial Differential Equations*, volume 7 of Abel Symposia, (2012), 37–52.

[16] L. Caffarelli; A. Figalli; Regularity of solutions to the parabolic fractional obstacle problem, *J. Reine Angew. Math.*, **680** (2013), 191–233.

[17] R. Cont, P. Tankov; *Financial Modelling with Jump Processes*, Chapman & Hall/CRC Financial. Math. Ser. 2004.

[18] A. De Pablo, F. Quirós, A. Rodriguez, J. L. Vázquez; A fractional porous medium equation. *Adv. Math.*, **226** (2011), no.2, 1378–1409.

[19] A. De Pablo, F. Quirós, A. Rodriguez, J. L. Vázquez; A general fractional porous medium equation, *Comm. Pure Applied Math.*, **65** (2012), 1242–1284.

[20] L. M. Del Pezzo, A. Quaas; A hopf's lemma and a strong maximum principle for the fractional $p$-Laplacian, *J. Differential Equations*, **263** (2017), 765–778.

[21] E. Di Nezza, G. Palatucci, E. Valdinoci; Hitchhiker's guide to the fractional Sobolev spaces, *Bull. Sci. Math.*, **136** (2012), 521-573.

[22] J. I. Díaz; New applications of monotonicity methods to a class of non-monotone parabolic quasilinear sub-homogeneous problems, to appear in *Pure and Applied Functional Analysis*.

[23] G. Franzina, G. Palatucci; Fractional $p$-eigenvalues, *Riv. Mat. Univ. Parma*, **5** (2014), 315–328.

[24] C. G. Gal, M. Warma; On some degenerate non-local parabolic equation associated with the fractional $p$-Laplacian, *Dyn. Partial Differ. Equ.*, **14** (2017), 47–77.

[25] J. Giacomoni, S. Tiwari; Existence and global behavior of solutions to fractional $p$-Laplacian parabolic problems, *Electron. J. Differential Equations*, 2018, no. 44, 20 p.

[26] J. Giacomoni, T. Mukherjee, K. Sreenadh; Existence and stabilization results for a singular parabolic equation involving the fractional Laplacian, *Discrete Contin. Dyn. Syst. Ser. S*, **12**, 2 (2019), 311–337.

[27] A. Iannizzotto, S. Mosconi, M. Squassina; Global Hölder regularity for the fractional $p$-Laplacian, *Rev. Mat. Iberoam.*, **32** (2016), 1353–1392.

[28] M. Liao, Q. Liu, H. Ye; Global existence and blow-up of weak solutions for a class of fractional $p$-Laplacian evolution equations, *Adv. Nonlinear Anal.*, **9** (2020), no. 1, 1569–1591.

[29] G. Molica Bisci, V. Radulescu, R. Servadei; *Variational Methods for Nonlocal Fractional Problems*, Cambridge University Press, Encyclopedia Math. Appl. (2016).

[30] J. M. Mazón, J. D. Rossi, J. Toledo; Fractional $p$-Laplacian Evolution Equations, *J. Math. Pures Appl.* (9), **105** (2016), 810–844.

[31] N. Pan, B. Zhang, J. Cao; Degenerate Kirchhoff-type diffusion problems involving the fractional $p$-Laplacian, *Nonlinear Anal. Real World Appl.*, **37** (2017), 56–70.

[32] P. Pucci, M. Xiang, B. Zhang; A diffusion problem of Kirchhoff type involving the nonlocal fractional $p$-Laplacian, *Discrete Contin. Dyn. Syst.*, **37** (2017), 4035–4051.

[33] D. Puhst; On the evolutionary fractional $p$-Laplacian, *Appl. Math. Res. Express. AMRX*, **2** (2015), 253–273.

[34] X. Ros-Oton; Nonlocal elliptic equations in bounded domains: a survey, *Publ. Mat.*, **60**, no. 1, 3-26.

[35] M. Strömqvist; Local boundedness of solutions to non-local parabolic equations modeled on the fractional $p$-Laplacian, *J. Differential Equations*, **266** (2019), no. 12, 7948–7979.

[36] L. Tang; Random homogenization of $p$-Laplacian with obstacles in perforated domain and related topics. Ph.D Dissertation, The University of Texas at Austin, 2011.

[37] J. L. Vázquez; Nonlinear diffusion with fractional Laplacian operators, *Nonlinear Partial Differential Equations*, Abel Symposia **7** (2012) 271–298.

[38] J. L. Vázquez; The Dirichlet problem for the fractional $p$-Laplacian evolution equation, *J. Differential Equations*, **260** (2016), 6038–6056.

[39] J. L. Vázquez; The evolution fractional $p$-Laplacian equation in $\mathbb{R}^N$. Fundamental solution and asymptotic behaviour. arXiv preprint arXiv:2004.05799, (2020).

[40] M. Warma; The fractional Neumann and Robin type boundary conditions for the regional fractional $p$-Laplacian, *NoDEA Nonlinear Differential Equations Appl.*, **23** (2016), no. 1, 46 p.
[41] M. Warma; Local Lipschitz continuity of the inverse of the fractional $p$-Laplacian, Hölder type continuity and continuous dependence of solutions to associated parabolic equations on bounded domains, *Nonlinear Anal.*, **135** (2016), 129-157.

[42] M. Xiang, B. Zhang, V. D. Rădulescu; Existence of solutions for perturbed fractional $p$-Laplacian equations, *J. Differential Equations*, **260** (2016), no. 2, 1392-1413.

JACQUES GIACOMONI
UNIVERSITÉ DE PAU ET DES PAYS DE L'ADOUR, CNRS, E2S LMAP (UMR 5142), IPRA, AVENUE DE L'UNIVERSITÉ-F-64013 PAU, FRANCE

Email address: jacques.giacomoni@univ-pau.fr

ABDELHAMID GOUASMIA
LABORATOIRE D'ÉQUATIONS AUX DÉRIVÉES PARTIELLES NON LINÉAIRES ET HISTOIRE DES MATHÉMATIQUES, ÉCOLE NATIONALE SUPÉRIEURE, B.P. 92, VIEUX KOUBA, 16050 ALGIERS, ALGERIA

Email address: gouasmia.abdelhamid@gmail.com

ABDELHAFID MOKRANE
LABORATOIRE D'ÉQUATIONS AUX DÉRIVÉES PARTIELLES NON LINÉAIRES ET HISTOIRE DES MATHÉMATIQUES, ÉCOLE NATIONALE SUPÉRIEURE, B.P. 92, VIEUX KOUBA, 16050 ALGIERS, ALGERIA

Email address: abdelhafid.mokrane@ens-kouba.dz