Proof. Step 1. We collect some useful estimates.

We denote \[p^{(2)}(x,t)=\frac{1}{(4\pi t)^{\frac{n}{2}}}\exp(-\frac{|x|^2}{4t})\] and \(p^{(2s)}(x,t)\) are the heat kernels associated with \(-\Delta\) and \((-\Delta)^s\), respectively. We also define \[\begin{align} \tilde{p}^{(2)}(x,t)=\frac{1}{(4\pi t)^{\frac{n}{2}}}\exp(-\frac{|x|^2}{16t}),\quad\hat{p}^{(2)}(x,t)=\frac{1}{(4\pi t)^{\frac{n}{2}}}\exp(-\frac{|x|^2}{t}), \end{align}\] and \[\begin{align} q_1(x,t)&=\begin{cases} \hat{p}^{(2)}(x,t),\quad\text{for }|x|^2<t<|x|^{2s}\leq 1,\\ \max\{\hat{p}^{(2)}(x,t),{p}^{(2s)}(x,t)\},\quad\text{for }t<|x|^{2}\leq 1,\\ {p}^{(2)}(x,t),\quad\text{for }|x|^{2s}\leq t\leq 1,\\ {p}^{(2s)}(x,t),\quad\text{for } t\geq 1\text{ or }|x|\geq 1, \end{cases}\\ q_2(x,t)&=\begin{cases} p^{(2)}(x,t),\quad\text{for }|x|^2<t<|x|^{2s}\leq 1,\\ \max\{\tilde{p}^{(2)}(x,t),{p}^{(2s)}(x,t)\},\quad\text{for }t<|x|^{2}\leq 1,\\ {p}^{(2)}(x,t),\quad\text{for }|x|^{2s}\leq t\leq 1,\\ {p}^{(2s)}(x,t),\quad\text{for } t\geq 1\text{ or }|x|\geq 1. \end{cases} \end{align}\] Theorem 2.13 in [47] shows that there exists a positive constant \(C_1\) such that for any \((x,t)\in\mathbb{R}^n\times (0,+\infty)\), \[\label{mix-kernel-estimate} \frac{1}{C_1}q_1(x,t)\leq H(x,t)\le C_1 q_2(x,t).\tag{5}\] Note that by Blumenthal and Getoor [49], there exists a positive constant \(C_2\) such that for any \((x,t)\in\mathbb{R}^n\times(0,+\infty)\), \[\label{s-heat32kernel-estimate} \frac{1}{C_2}\min\{t^{-\frac{n}{2s}},t|x|^{-n-2s}\}\leq p^{(2s)}(x,t)\leq C_2\min\{t^{-\frac{n}{2s}},t|x|^{-n-2s}\}.\tag{6}\] We observe that \[\label{min-cases} \min\{t^{-\frac{n}{2s}},t|x|^{-n-2s}\}= \begin{cases} t|x|^{-n-2s},\quad\text{if }|x|^{2s}>t,\\ t^{-\frac{n}{2s}},\quad\text{if }|x|^{2s}\leq t. \end{cases}\tag{7}\]

Step 2. We prove ?? .

Assume that \(0<|x|\leq 1\), then \[\begin{align} \Gamma(x)&\leq C_1\int^\infty_0q_2(x,t)dt\\ &= C _1\left(\int^{|x|^2}_0 \max\{\tilde{p}^{(2)}(x,t),{p}^{(2s)}(x,t)\}dt+\int^{1}_{|x|^{2}}{p}^{(2)}(x,t)dt+\int^{\infty}_1 p^{(2s)}(x,t)dt\right) \end{align}\] First, note that \(t<|x|^2\) implies that \(t<|x|^{2s}\), hence by 6 7 , we have \[\int^{|x|^2}_0 {p}^{(2s)}(x,t)dt\leq C \int^{|x|^2}_0 t|x|^{-n-2s}dt=C|x|^{4-n-2s}.\] Next, by the variable substitution \(z=|x|^2/t\), we obtain \[\begin{align} \int^{|x|^2}_0 \tilde{p}^{(2)}(x,t)dt= C\int^{|x|^2}_0 t^{-\frac{n}{2}}e^{-\frac{|x|^2}{16t}} dt=C\int^{\infty}_1|x|^{2-n}z^{\frac{n}{2}-2}e^{-\frac{z}{16}}dz=C|x|^{2-n}, \end{align}\] and \[\begin{align} \int^{1}_{|x|^{2}}{p}^{(2)}(x,t)dt= C \int^{1}_{|x|^{2}}t^{-\frac{n}{2}}e^{-\frac{|x|^2}{4t}} dt = C\int^{1}_{|x|^2}|x|^{2-n}z^{\frac{n}{2}-2}e^{-\frac{z}{4}}dz\leq C |x|^{2-n}. \end{align}\] In addition, 6 and 7 yield that \[\int^{\infty}_1 p^{(2s)}(x,t)dt \leq C\int^{\infty}_1 t^{-\frac{n}{2s}}dt\leq C.\] Combining the above facts, we get \(\Gamma(x)\leq C|x|^{2-n}\) for all \(0<|x|\leq 1\). The argument for proving \(\Gamma(x)\geq C'|x|^{2-n}\) for all \(0<|x|\leq 1\) is analogous.

Step 3. We prove ?? .

Assume that \(|x|\geq 1\), then by 7 , we have \[\begin{align} \Gamma(x)&\leq C_1\int^\infty_0q_2(x,t)dt= C_1\int^\infty_0p^{(2s)}(x,t)dt\leq C\left(\int_0^{|x|^{2s}}t|x|^{-n-2s}dt+\int^{\infty}_{|x|^{2s}}t^{-\frac{n}{2s}}dt\right)\\ &\leq C|x|^{2s-n}. \end{align}\] Similar discussion yields \(\Gamma(x)\geq C'|x|^{2s-n}\) for all \(|x|\geq 1\).

Now, combining ?? and ?? , \(\Gamma\) is well defined in \(\mathbb{R}^n\backslash\{0\}\). Then, by Theorem 3.1 in [48], we know that \(H\) is nonnegative, radially symmetric, and nonincreasing with respect to \(r=|x|\), hence so is \(\Gamma\).

Step 4. We verify that \(\Gamma\) is a fundamental solution of \(-\Delta +(-\Delta)^s.\) Taking the Fourier transform, equation ?? becomes equivalent to \[\mathcal{F}(\Gamma)(\xi)=\frac{1}{|\xi|^2+|\xi|^{2s}},\] where \(\mathcal{F}(\cdot)\) denotes the Fourier transform. Hence, it suffices to show that for any \(\phi\in C^\infty_c(\mathbb{R}^n)\), it holds \[\label{equ-dis-test} \langle \mathcal{F}^{-1}\left(\frac{1}{|\xi|^2+|\xi|^{2s}}\right),\phi\rangle= \langle\Gamma,\phi\rangle.\tag{8}\]

First, by ?? , we have \(\Gamma\phi\in L^1(\mathbb{R}^n)\). Next, for a fixed \(t>0\), we observe that \(A(x,\xi):=e^{-t(|\xi|^2+|\xi|^{2s})+ix\cdot\xi}\phi(x)\in L^1(\mathbb{R}^n\times\mathbb{R}^n)\). Moreover, since \(\mathcal{F}^{-1}(\phi)\) is a Schwartz function, then \(|\mathcal{F}^{-1}(\phi)(\xi)|(1+|\xi|)^{N}\) is bounded for any \(N\in\mathbb{N}\) and thus \(\mathcal{F}^{-1}(\phi)(\xi)\frac{1}{|\xi|^2+|\xi|^{2s}}\in L^1(\mathbb{R}^n)\). With these facts, employing 4 and applying Fubini’s theorem repeatedly, we derive \[\begin{align} \langle\Gamma,\phi\rangle&=\int_{\mathbb{R}^n}\int^{\infty}_{0}H(x,t)\phi(x) dtdx\\ &=\int^{\infty}_{0}\int_{\mathbb{R}^n}H(x,t)\phi(x) dxdt\\ &=\int^{\infty}_{0}\int_{\mathbb{R}^n}\int_{\mathbb{R}^n}e^{-t(|\xi|^2+|\xi|^{2s})+ix\cdot\xi}\phi(x)d\xi dx dt\\ &=\int^{\infty}_{0}\int_{\mathbb{R}^n}e^{-t(|\xi|^2+|\xi|^{2s})}\left(\int_{\mathbb{R}^n}e^{ix\cdot\xi}\phi(x)dx\right) d\xi dt\\ &=\int^{\infty}_{0}\int_{\mathbb{R}^n}e^{-t(|\xi|^2+|\xi|^{2s})}\mathcal{F}^{-1}(\phi)(\xi) d\xi dt\\ &=\int_{\mathbb{R}^n}\mathcal{F}^{-1}(\phi)(\xi)\left(\int^{\infty}_{0}e^{-t(|\xi|^2+|\xi|^{2s})} dt\right)d\xi\\ &=\int_{\mathbb{R}^n}\mathcal{F}^{-1}(\phi)(\xi)\frac{1}{|\xi|^2+|\xi|^{2s}}d\xi=\langle \mathcal{F}^{-1}\left(\frac{1}{|\xi|^2+|\xi|^{2s}}\right),\phi\rangle. \end{align}\] Therefore, we obtain 8 and complete the proof. ◻

Proof. We use the smooth approximation. Fix \(\eta\in C^\infty_c(B_1)\) is radial about 0 satisfying \(\eta\geq 0\) and \(\int_{\mathbb{R}^n}\eta dx=1\). For \(\varepsilon\in (0,1)\) and a function \(v\), we denote \(v_{\varepsilon}(x)=\int_{\mathbb{R}^n}v(x-y)\varepsilon^{-n}\eta(\varepsilon^{-1}y)dy\).

Now, We claim that \(u_\varepsilon\in C^{\infty}(\mathbb{R}^n)\) satisfies the following properties:

\(\mathrm{(i)}\) \(u_\varepsilon\in \mathcal{L}^{2s}(\mathbb{R}^n)\). Note that for any \(y\in B_{\varepsilon}\) and \(x\in\mathbb{R}^n\), we have \[\frac{1+|x-y|^{n+2s}}{1+|x|^{n+2s}}\leq C,\] hence \[\begin{align} \int_{\mathbb{R}^n}\frac{|u_\varepsilon(x)|}{1+|x|^{n+2s}}dx&=\int_{\mathbb{R}^n}\int_{\mathbb{R}^n}\frac{|u(x-y)|}{1+|x|^{n+2s}}dx\varepsilon^{-n}\eta(\varepsilon^{-1}y)dy\\ &\leq C\int_{\mathbb{R}^n}\int_{\mathbb{R}^n}\frac{|u(x-y)|}{1+|x-y|^{n+2s}}dx\varepsilon^{-n}\eta(\varepsilon^{-1}y)dy\\&=C\int_{\mathbb{R}^n}\frac{|u(x)|}{1+|x|^{n+2s}}dx<\infty. \end{align}\]

\(\mathrm{(ii)}\) \(-\Delta u_{\varepsilon}+(-\Delta)^s u_{\varepsilon}\geq 0\) in \(\mathbb{R}^n\) in the classical sense. Take any \(\phi\in C^{\infty}_c(\mathbb{R}^n)\) with \(\phi\geq 0\). Note that \[(-\Delta)^s\phi(x)=\frac{C_{n,s}}{2}\int_{\mathbb{R}^n}\frac{2\phi(x)-\phi(x+y)-\phi(x-y)}{|y|^{n+2s}}dy,\] we easily verify that \(((-\Delta)^s\phi)_\varepsilon=(-\Delta)^s\phi_\varepsilon\). We also have \((-\Delta\phi)_\varepsilon=-\Delta\phi_\varepsilon\). Hence, \[\begin{align} \int_{\mathbb{R}^n}[-\Delta+(-\Delta)^s]u_\varepsilon\phi dx&=\int_{\mathbb{R}^n}u_\varepsilon[-\Delta+(-\Delta)^s]\phi dx=\int_{\mathbb{R}^n}u([-\Delta+(-\Delta)^s]\phi)_\varepsilon dx\\&=\int_{\mathbb{R}^n}u[-\Delta+(-\Delta)^s]\phi_\varepsilon dx\geq 0, \end{align}\] where the last inequality follows from 9 . By the arbitrariness of \(\phi\), we obtain \(-\Delta u_{\varepsilon}+(-\Delta)^s u_{\varepsilon}\geq 0\) pointwise in \(\mathbb{R}^n\).

\(\mathrm{(iii)}\) \(\liminf_{|x|\rightarrow\infty}u_\varepsilon(x)\geq 0.\) By 10 , for any \(\delta>0\), there exists some \(R>0\) such that \(u(x)\geq -\delta\) for a.e. \(x\in B_R^c\). Note that for any \(x\in B_{R+\varepsilon}^c\) and \(y\in B_\varepsilon\), we have \(|x-y|\geq R\), hence \[u_{\varepsilon}(x)=\int_{B_\varepsilon}u(x-y)\varepsilon^{-n}\eta(\varepsilon^{-1}y)dy\geq -\delta\int_{B_\varepsilon}\varepsilon^{-n}\eta(\varepsilon^{-1}y)dy=-\delta.\] Therefore, we obtain that \(\liminf_{|x|\rightarrow\infty}u_\varepsilon(x)\geq 0.\)

From the above properties of \(u_{\varepsilon}\), we have \(u_\varepsilon\geq 0\) in \(\mathbb{R}^n\). Otherwise, \(\inf_{\mathbb{R}^n}u_\varepsilon<0\). By \(\mathrm{(iii)}\), there exists an \(x_0\in\mathbb{R}^n\) such that \(u_\varepsilon(x_0)=\min_{\mathbb{R}^n}u_{\varepsilon}<0\). However, we have \[[-\Delta +(-\Delta)^s]u_\varepsilon(x_0)\leq (-\Delta)^s u_\varepsilon(x_0)=\int_{\mathbb{R}^n}\frac{ u_\varepsilon(x_0)- u_\varepsilon(y)}{|x_0-y|^{n+2s}}dy<0,\] which contradicts \(\mathrm{(ii)}\). Therefore, we obtain \(u_\varepsilon\geq 0\) in \(\mathbb{R}^n\). Let \(\varepsilon\rightarrow 0\), we get \(u\geq 0\) a.e. in \(\mathbb{R}^n\). ◻

Proof. For \(R>0\), denote \[u_R= \begin{cases} u\quad\text{in } B_R,\\ 0\quad\text{in } \mathbb{R}^n\backslash B_R. \end{cases}\] Let \[w_R(x)=\int_{\mathbb{R}^n} \Gamma(x-y)u^q_R(y)dy.\] Then by Proposition 4, one can easily check that \(w_R\in\mathcal{L}^{2s}(\mathbb{R}^n)\) is a distributional solution to \[[-\Delta +(-\Delta)^s]w_R=u_R^q\quad\text{in }\mathbb{R}^n.\] with the asymptotic behavior\[\label{wR-decay} w_R(x)\rightarrow 0,\;as\;|x|\rightarrow\infty.\tag{11}\] Set \(v_R=u-w_R\), then \(v_R\) is a distributional solution of \[\label{vR-equation} [-\Delta +(-\Delta)^s]v_R=u^q-u_R^q\geq 0\quad\text{in }\mathbb{R}^n,\tag{12}\] with \(\liminf_{|x|\rightarrow\infty}v_R(x)\geq 0\). Applying Lemma 1 to \(v_R\), we have \(v_R\geq 0\) a.e. in \(\mathbb{R}^n\). Letting \(R\rightarrow+\infty\) and applying the monotone convergence theorem, we obtain ?? . ◻

Proof of Theorem 1. Case 1: \(1<q< \frac{n}{n-2s}\).

We present two different proofs.

Proof by test functions: Let \(u\) be a positive supersolution of 1 and \(\eta\in C^\infty_c(B_2)\) satisfy \(0\leq \eta\leq 1\), \(\eta|_{B_1}=1\). For any \(R>0\), we set \(\eta_{R}=\eta(\frac{\cdot}{R})\) For some \(m\in\mathbb{Z}_{+}\) to be determined, we use \(\eta^m_{R}\) as a test function for 1 . Note that \[|\nabla \eta_R|\leq \frac{C}{R},\quad |\nabla^2\eta_R|\leq \frac{C}{R^2}.\] Hence, we have \[\int_{\mathbb{R}^n}u(-\Delta\eta^m_R)dx=-\int_{\mathbb{R}^n}u(m\eta^{m-1}_R\Delta\eta_R+m(m-1)\eta^{m-2}_R|\nabla\eta_R|^2)dx\leq \frac{C}{R^2}\int_{\mathbb{R}^n}u\eta^{m-1}_Rdx.\] Now we set \(m=\lfloor\frac{q}{q-1}\rfloor+1\), then \(m-1-\frac{m}{q}>0\). By Hölder inequality, we obtain that \[\label{Delta-eta} \begin{align} \int_{\mathbb{R}^n}u(-\Delta\eta^m_R)dx&\leq \frac{C}{R^2}\int_{\mathbb{R}^n}u\eta^{\frac{m}{q}}_R\eta^{m-1-\frac{m}{q}}_Rdx\\ &\leq \frac{C}{R^2}\left(\int_{\mathbb{R}^n}u^q\eta^m_Rdx\right)^{\frac{1}{q}}\left(\int_{\mathbb{R}^n}\eta^{(m-1-\frac{m}{q})(\frac{q}{q-1})}_Rdx\right)^{1-\frac{1}{q}}\\ &\leq CR^{n(1-\frac{1}{q})-2}\left(\int_{\mathbb{R}^n}u^q\eta^m_Rdx\right)^{\frac{1}{q}}. \end{align}\tag{17}\] Next, by the convexity of \(t\mapsto t^m\), we have \[(a^m-b^m)\leq ma^{m-1}(a-b),\;for all\;a,b\geq 0,\] hence \[\begin{align} \int_{\mathbb{R}^n}u(-\Delta)^s(\eta^m_R)dx&=\int_{\mathbb{R}^n}\int_{\mathbb{R}^n}u(x)\frac{\eta^m_R(x)-\eta^m_R(y)}{|x-y|^{n+2s}}dydx\\ &\leq \int_{\mathbb{R}^n}\int_{\mathbb{R}^n}u(x)\frac{m\eta^{m-1}_R(x)(\eta_R(x)-\eta_R(y))}{|x-y|^{n+2s}}dydx\\ &= \int_{\mathbb{R}^n}u(x)m\eta^{m-1}_R(x)\frac{1}{R^{2s}}\int_{\mathbb{R}^n}\frac{(\eta(\frac{x}{R})-\eta(\frac{y}{R}))}{|\frac{x}{R}-\frac{y}{R}|^{n+2s}}\frac{dy}{R^n}dx\\ &= \frac{1}{R^{2s}}\int_{\mathbb{R}^n}u(x)m\eta^{m-1}_R(x)[(-\Delta)^s\eta]\left(\frac{x}{R}\right)dx. \end{align}\] Note that \(supp(\eta_R)\subset B_{2R}\) and \(|(-\Delta)^s\eta|\leq C\) in \(B_2\) yield that \[\int_{\mathbb{R}^n}u(-\Delta)^s(\eta^m_R)dx\leq \frac{C}{R^{2s}}\int_{\mathbb{R}^n}u\eta^{m-1}_Rdx.\] Following a discussion similar to 17 , we obtain \[\label{s-Delta-eta} \int_{\mathbb{R}^n}u(-\Delta)^s(\eta^m_R)dx\leq CR^{n(1-\frac{1}{q})-2s}\left(\int_{\mathbb{R}^n}u^q\eta^m_Rdx\right)^{\frac{1}{q}}.\tag{18}\] On the other hand, by \(u\) is a supersolution of 1 , we get \[\label{mix-super} \int_{\mathbb{R}^n}u[(-\Delta)+(-\Delta)^s](\eta^m_R)dx \geq \int_{\mathbb{R}^n}u^p\eta^m_R dx.\tag{19}\] Combining 17 19 , we have \[\label{decay-estimate-inte} \left(\int_{\mathbb{R}^n}u^q\eta^m_Rdx\right)^{1-\frac{1}{q}}\leq C(R^{n(1-\frac{1}{q})-2s}+R^{n(1-\frac{1}{q})-2}).\tag{20}\] By \(q<\frac{n}{n-2s}\), we have \[n(1-\frac{1}{q})-2<n(1-\frac{1}{q})-2s<0.\] Hence, let \(R\rightarrow+\infty\) in 20 , we get \(u\equiv0\) in \(\mathbb{R}^n\), which contradicts \(u>0\) in \(\mathbb{R}^n\).

Proof by equivalent integral characterization: Let \(u\) be a positive supersolution of 1 . By ?? and ?? , we have, for \(|x|\geq 2\), \[u(x)\geq \int_{B_1}\Gamma(x-y)u^q(y)dy\geq\frac{1}{C}\int_{B_1}\frac{1}{|x-y|^{n-2s}}u^q(y)dy\geq \frac{C_1}{|x|^{n-2s}}.\] for some constant \(C_1>0\). Hence, for \(|x|\geq 4\), \[\begin{align} u(x)&\geq \int_{B_{\frac{|x|}{2}}(x)\backslash B_{\frac{|x|}{4}}(x)}\Gamma(x-y)u^q(y)dy\geq \frac{1}{C}\int_{B_{\frac{|x|}{2}}(x)\backslash B_{\frac{|x|}{4}}(x)}\frac{1}{|x-y|^{n-2s}}\frac{C^q_1}{|y|^{q(n-2s)}}dy\\ &\geq\frac{C_2}{|x|^{q(n-2s)-2s}},\dots \end{align}\] By induction, for \(j=1,2,\dots\), when \(|x|\geq 2^j\), we obtain \[\label{decay-alphaj} u(x)\geq \frac{C_j}{|x|^{\alpha_j}},\tag{21}\] where \(C_j>0\) and \[\alpha_1=n-2s,\quad\alpha_j=q\alpha_{j-1}-2s~(j\geq 2).\] Note that \(1<q<\frac{n}{n-2s}\) implies \(n-2s+\frac{2s}{1-q}<0\), hence, for \(j\) large, \[\alpha_j=q^{j-1}\alpha_1-2s(1+q+q^2+\cdots +q^{j-1})=\left(n-2s+\frac{2s}{1-q}\right)q^{j-1}-\frac{2s}{1-q}<0.\] We choose a \(j_0\) such that \(\alpha_{j_0}<0\). Then, by ?? and 21 , we get, for \(|x|\leq 2^{j_0-1}\), \[u(x)\geq \int_{\mathbb{R}^n\backslash B_{2^{j_0}}}\Gamma (x-y)u^q(y)dy\geq \frac{1}{C}\int_{\mathbb{R}^n\backslash B_{2^{j_0}}}\frac{1}{|x-y|^{n-2s}}\frac{C^q_{j_0}}{|y|^{q\alpha_{j_0}}}dy=+\infty.\] Therefore we obtain a contradiction, which shows that 1 has no positive supersolution.

Case 2: \(q=\frac{n}{n-2s}\).

Let \(R\geq 1\). By ?? , ?? , and ?? , we have \[\label{u-inf-estimate} \begin{align} u(x)&\geq \int_{B_R}\Gamma (x-y)u^q(y)dy\\ &=\left(\int_{B_R\cap\{|x-y|\leq 1\}}+\int_{B_R\cap\{|x-y|\geq 1\}}\right)\Gamma (x-y)u^q(y)dy\\ &\geq\frac{1}{C}\left[\int_{B_R\cap\{|x-y|\leq 1\}}\frac{u^q(y)}{|x-y|^{n-2}}dy+\int_{B_R\cap\{|x-y|\geq 1\}}\frac{u^q(y)}{|x-y|^{n-2s}}dy\right] \\ &\geq \frac{1}{C}\left[\int_{B_R\cap\{|x-y|\leq 1\}}u^q(y)dy+\frac{1}{(|x|+R)^{n-2s}}\int_{B_R\cap\{|x-y|\geq 1\}}{u^q(y)}dy\right]\\ &\geq \frac{1}{C}\frac{1}{(|x|+R)^{n-2s}}\int_{B_R}u^q(y)dy. \end{align}\tag{22}\] Hence, \[\int_{B_R}u^q(x)dx\geq \frac{1}{C^q}\int_{B_R}\frac{1}{(|x|+R)^{n}}dx\left(\int_{B_R}u^q(y)dy\right)^q\geq c\left(\int_{B_R}u^q(y)dy\right)^q,\] where \(c\) is independent of \(R\). Let \(R\rightarrow +\infty\), we get \(u\in L^q(\mathbb{R}^n)\). Furthermore, by 22 , we have \[\int_{B_{2R}\backslash B_R} u^q(x)dx\geq \frac{1}{C^q}\int_{B_{2R}\backslash B_R}\frac{1}{(|x|+R)^{n}}dx\left(\int_{B_R}u^q(y)dy\right)^q\geq c'\left(\int_{B_R}u^q(y)dy\right)^q,\] where \(c'\) is independent of \(R\). Let \(R\rightarrow +\infty\), we obtain \(\int_{\mathbb{R}^n}u^q(y)dy=0\) and hence \(u\equiv0\). This contradicts the fact that \(u>0\).

Case 3: \(q>\frac{n}{n-2s}\).

We will construct a positive supersolution of 1 . Let \(f\) be defined by 13 with 14 15 . We first compute \((-\Delta)^s f\). Fix \(x,y\in\mathbb{R}^n\) with \(x\neq0\). We set \(r=|x|\) and write \(x=rA_x e_1\) for some rotation \(A_x\) in \(\mathbb{R}^n\). We also write \(y=r A_x z\) for some \(z\in\mathbb{R}^n\). Let \(g_r\) be defined by 16 , we have \[f(x)-f(y)=r^{2s-n+\sigma}(g_r(e_1)-g_r(z)).\] Hence, we obtain \[\label{s-delta-f} \begin{align} (-\Delta)^s f(x)&=\mathrm{p.v.}\int_{\mathbb{R}^n}\frac{f(x)-f(y)}{|x-y|^{n+2s}}dy= \mathrm{p.v.}\int_{\mathbb{R}^n}\frac{r^{2s-n+\sigma}(g_r(e_1)-g_r(z))}{r^{n+2s}|e_1-z|^{n+2s}}r^ndz\\ &=r^{\sigma-n}(-\Delta)^s g_r(e_1)=f(x)^q\frac{r^{\sigma-n}}{(1+r^2)^{\frac{(2s-n+\sigma)q}{2}}}(-\Delta)^s g_r(e_1). \end{align}\tag{23}\] On the other hand, a direct computation yields \[\label{delta-f} \begin{align} (-\Delta)f(x)&=-(1+r^2)^{\frac{2s-n+\sigma}{2}-1}(2s-n+\sigma)\left((2s-n+\sigma-2)\frac{|x|^2}{1+|x|^2}+n\right)\\ &=f(x)^q\frac{r^{\sigma-n}}{(1+r^2)^{\frac{(2s-n+\sigma)q}{2}}}h(r), \end{align}\tag{24}\] where \[h(r):=-\frac{(1+r^2)^{\frac{2s-n+\sigma}{2}-1}}{r^{2s-n+\sigma-2}}\cdot r^{2s-2}(2s-n+\sigma)\left((2s-n+\sigma-2)\frac{r^2}{1+r^2}+n\right).\] Note that \[\label{h-infty-0} \lim_{r\rightarrow+\infty}h(r)=0.\tag{25}\] By 14 , we have \[\label{com-infty} \lim_{r\rightarrow\infty}\frac{r^{\sigma-n}}{(1+r^2)^{\frac{(2s-n+\sigma)q}{2}}}=+\infty.\tag{26}\] Combining Lemma 2 and 23 26 , we can find a \(R>0\) such that for any \(|x|\geq R\), \[\label{exter-equation} [-\Delta +(-\Delta)^s] f(x)\geq f(x)^q.\tag{27}\] Furthermore, for \(\delta>0\) and \(u\) a positive supersolution of 1 , we consider \[u_\delta(x)=\delta^\frac{2}{q-1}u(\delta x),\quad x\in\mathbb{R}^n.\] Then \(u_\delta\) satisfies \[\label{u95delta-equation} -\Delta u_\delta+\delta^{2(1-s)}(-\Delta)^su_\delta\geq u_\delta^q\quad\text{in }\mathbb{R}^n.\tag{28}\] Hence, to construct a positive solution of 1 , we only need to construct a positive supersolution of 28 for some \(\delta>0\). By Lemma 2, we know that \((-\Delta)^s f\) is a continuous positive function in \(\mathbb{R}^n\). Now, we take \[\delta=\max\left\{1,\max_{|x|\leq R}\left(\frac{(f^q+\Delta f)_+}{(-\Delta)^sf}\right)^{\frac{1}{2(1-s)}}\right\},\] then \(-\Delta f+\delta^{2(1-s)}(-\Delta)^sf\geq f^q\) in \(\{|x|\leq R\}\). Combining this with 27 shows that \(u_\delta=f\) is a supersolution of 28 . Consequently, \(\delta^{-\frac{2}{q-1}} f( \cdot/\delta)\) is a supersolution of 1 , which completes the proof. ◻