Nonlocal Harnack inequalities in the Heisenberg group

We deal with a wide class of nonlinear integro-differential problems in the Heisenberg-Weyl group $\mathbb{H}^n$, whose prototype is the Dirichlet problem for the $p$-fractional subLaplace equation. These problems arise in many different contexts in quantum mechanics, in ferromagnetic analysis, in phase transition problems, in image segmentations models, and so on, when non-Euclidean geometry frameworks and nonlocal long-range interactions do naturally occur. We prove general Harnack inequalities for the related weak solutions. Also, in the case when the growth exponent is $p=2$, we investigate the asymptotic behavior of the fractional subLaplacian operator, and the robustness of the aforementioned Harnack estimates as the differentiability exponent $s$ goes to $1$.


Introduction
We deal with a very general class of nonlinear nonlocal operators, which include, as a particular case, the fractional subLaplacian. Precisely, let Ω be a bounded domain in the Heisenberg-Weyl group H n , and let g be in the fractional Sobolev space W s,p (H n ), for any s ∈ (0, 1) and any p > 1. We shall prove general Harnack inequalities for the weak solutions to the following class of nonlinear integro-differential problems, where f = f (·, u) belongs to L ∞ loc (H n ) uniformly in Ω, and L is the operator defined by with d o being a homogeneous norm on H n , and Q = 2n + 2 the usual homogeneous dimension of H n . The symbol P. V. in the display above stands for "in the principal value sense". We immediately refer the reader to Section 2 below for the precise definitions of the involved quantities and related properties, as well as for further observations in order to relax some of the assumptions listed in the present section.
Integral-differential operators in the form as in (1.2) do arise as a generalization of the fractional subLaplacian on the Heisenberg group, naturally defined in the fractional Sobolev space H s (H n ) for any s ∈ (0, 1) as follows (−∆ H n ) s u(ξ) := C(n, s) P. V.
where | · | H n is the standard homogeneous norm of H n , and C(n, s) is a positive constant which depends only on n and s. In this fashion, the prototype of the wide class of problems in (1.1) reads as follows, In the last decades, a great attention has been focused on the study of problems involving fractional equations, both from a pure mathematical point of view and for concrete applications since they naturally arise in many different contexts. Despite its relatively short history, the literature is really too wide to attempt any comprehensive treatment in a single paper; we refer for instance to the paper [19] for an elementary introduction to fractional Sobolev spaces and for a quite extensive (but still far from exhaustive) list of related references. For what concerns specifically the family of equations in (1.1) and the corresponding energy functionals, both in the nonlocal and in the local framework, the link with several concrete models arises from many different contexts in Probability (e.g., in non-Markovian coupling for Brownian motions [5]), in Physics (e. g., in group theory in quantum mechanics [46], in ferromagnetic trajectories [34], in image segmentation models [13], in phase transition problems described by Ising models [38], and many others), where the analysis in sub-Riemannian geometry revealed to be decisive. In this respect, as proven in the literature, Harnack-type inequalities constitute a fundamental tool of investigation.
Let us focus now merely on regularity and related results in the fractional panorama in the Heisenberg group. It is firstly worth stressing that one can find various definitions of the involved operator and related extremely different approaches. In the linear case when p = 2, an explicit integral definition can be found in the relevant paper [41], where several Hardy inequalities for the conformally invariant fractional powers of the sublaplacian are proven, and [11] for related Hardy and uncertainty inequalities on general stratified Lie groups involving fractional powers of the Laplacian; we also refer to [2], where, amongst other important results, Morrey and Sobolev-type embeddings are derived for fractional order Sobolev spaces. Still in the linear case when p = 2, very relevant results have been obtained based on the construction of fractional operators via a Dirichlet-to-Neumann map associated to degenerate elliptic equations, as firstly seen for the Euclidean framework in the celebrated Caffarelli-Silvestre s-harmonic extension. For this, we would like to mention the very general approach in [23]; the Liouville-type theorem in [12]; the Harnack and Hölder results in Carnot groups in [20]; the connection with the fractional perimeters of sets in Carnot group in [21].
For what concerns the more general situation as in (1.2) when a p-growth exponent is considered, in our knowledge, a regularity theory is very far from be complete; nonetheless, very interesting estimates have been recently proven, as, e. g., in [26,45], and in our recent paper [33] where local boundedness and Hölder estimates have been proven for the weak solutions to (1.1).
In order to state our main results, we need to introduce a special quantity which plays a central role when dealing with nonlocal operators. Namely, we define the nonlocal tail Tail(u; ξ 0 , R) of a function u centered in ξ 0 ∈ H n of radius R > 0, (1.5) We immediately notice that the quantity above is finite whenever u ∈ L q (H n ), with q ≥ p − 1. The nonlocal tail in (1.5) can be seen as the natural generalization in the Heisenberg setting of that originally introduced in [18,17], and subsequently revealed to be decisive in the analysis of many different nonlocal problems when a fine quantitative control of the long-range interactions is needed; see for instance the subsequent results proven in [30,28,31,29] and the references therein.
Our main result reads as follows, Theorem 1 (Nonlocal Harnack inequality) For any s ∈ (0, 1) and where Tail(·) is defined in (1.5), u − := max{−u, 0} is the negative part of the function u, and c depends only on n, s, p, and the structural constant Λ defined in (2.2).
Notice that in the case when u is nonnegative in the whole H n , the formulation in (1.6) does reduce to that of the classical Harnack inequality.
In the particular situation when u is merely a weak supersolution to Problem (1.1), still in analogy with the classical case when s = 1, a weak Harnack inequality can be proven, as stated in the following Theorem 2 (Nonlocal weak Harnack inequality) For any s ∈ (0, 1) and any p ∈ (1, ∞), let u ∈ W s,p (H n ) be a weak supersolution to (1.1) such that u ≥ 0 in B R ⊆ Ω. Then, for for any s − Q/p < ε < s and t < (p−1)s Tail(·) is defined in (1.5), and u − := max{−u, 0} is the negative part of the function u.
The constant c depends only on n, s, p, and the structural constant Λ defined in (2.2).
As expected, a nonlocal tail contribution should be still taken into account, and, again, such a contribution in (1.7) will disappear in the case when the function u is nonnegative in the whole H n .
It is now worth noticing that the main difficulty into the treatment of the equation in (1.1) lies in the very definition of the leading operator L defined in (1.2), which combines the typical issues given by its nonlocal feature together with the ones given by its nonlinear growth behaviour and with those naturally arising from the non-Euclidean geometrical structure.
For this, some very important tools recently introduced in the nonlocal theory and successfully applied in the fractional sublaplacian on the Heisenberg group, as the aforementioned Caffarelli-Silvestre s-harmonic extension, and the approach via Fourier representation, as well as other successful tools, like for instance the commutator estimates in [43], the pseudo-differential commutator compactness in [36], and many others, seem not to be adaptable to the framework we are dealing with. However, even in such a nonlinear non-Euclidean framework, we will be able to extend part of the strategy developed in [17] where nonlocal Harnack inequalities have been proven for the homogeneous version of the analogue of problem (1.1) in the Euclidean framework. Further efforts are also needed due to the presence of the non homogeneous datum f , as well as in order to deal with the limit case when sp = Q, both of them are novelty even with respect to the results proven in the Euclidean framework in [17].
Let now focus on the linear case when p = 2 when the datum f in the right-hand side of the equation in (1.1) is zero. It is worth mentioning that the necessity of the presence of the tail term in the Harnack inequalities stated above is a very recent achievement. Indeed, during the last decades, the validity of the classical Harnack inequality without extra positivity assumptions on the solutions has been an open problem in the nonlocal setting, and more in general for integro-differential operators of the form in (1.1) even in the Euclidean framework. An answer has been eventually given by Kassmann in his breakthrough papers [24,25], where a simple counter-example is provided in order to show that positivity cannot be dropped nor relaxed even in the most simple case when L does coincide with the fractional Laplacian operator (−∆) s ; see Theorem 1.2 in [24]. The same author proposed a new formulation of the Harnack inequality without requiring the additional positivity on solutions by adding an extra term, basically a natural tailtype contribution on the right-hand side, in accordance with the result presented here; see Theorem 3.1 in [25], where the robustness of the estimates as s goes to 1 is also presented. In the same spirit, we also investigate the special linear case in which L does reduce to the pure fractional subLaplacian, namely problem (1.4). Firstly, as expected, we prove that the fractional subLaplacian (−∆ H n ) s effectively converges to the standard subLaplacian −∆ H n as s goes to 1, as stated in Proposition 1 below. For this, we shall carefully estimate the weighted second order integral form of the fractional subLaplacian with the aim of a suitable Mac-Laurin-type expansion in the Heisenberg group.
Proposition 1 For any u ∈ C ∞ 0 (H n ) the following statement holds true Notice that a proof of the result above for fractional subLaplacian on Carnot groups can be found in the relevant paper [21], via heat kernel characterization.
Secondly, we revisit the proofs of Theorems 1-2 by taking care of the dependance of the differentiability exponent s in all the estimates, so that we are eventually able to obtain the results below in clear accordance with the analogous ones in the Euclidean framework [25], by proving that the nonlocal tail term will vanish when s goes to 1, in turn recovering the classical Harnack formulation.
Theorem 3 For any s ∈ (0, 1) let u ∈ H s (H n ) be a weak solution to (1.4) such that u ≥ 0 in B R (ξ 0 ) ⊂ Ω. Then, the following estimate holds true for any B r such that B 6r ⊂ B R , where Tail(·) is defined in (1.5) by taking p = 2 there, u − := max{−u, 0} is the negative part of the function u, and c = c (n, s).
Theorem 4 For any s ∈ (0, 1) let u ∈ H s (H n ) be a weak supersolution to (1.4), such that u ≥ 0 in B R (ξ 0 ) ⊆ Ω. Then, the following estimate holds true for any B r such that B 6r ⊂ B R , and any t < Q/(Q − 2s), where Tail(·) is defined in (1.5) by taking p = 2 there, u − := max{−u, 0} is the negative part of the function u, and c = c (n, s).
Further developments. Starting from the results proven in the present paper, several questions naturally arise.
• Firstly, it is worth remarking that we treat general weak solutions, namely by truncation and dealing with the resulting error term as a right hand-side, in the same flavour of the papers [17,18,29], in the spirit of De Giorgi-Nash-Moser. However, one could approach the same family of problems by focusing solely to bounded viscosity solutions in the spirit of Krylov-Safonov.
• Still in clear accordance with the Euclidean counterpart, one would expect selfimproving properties of the solutions to (1.1). For this, one should extend the recent nonlocal Gehring-type theorems proven in [30,43,44].
• One could expect nonlocal Härnack inequalities and other regularity results for the solutions to a strictly related class of problems; that is, by adding in (1.1) a second integral-differential operators, of differentiability exponent t > s and summability growth q > 1, controlled by the zero set of a modulating coefficient: the so-called nonlocal double phase problem, in the same spirit of the Euclidean case treated in [16,7], starting from the pioneering results in the local case, when s = 1, by Colombo and Mingione; see for instance [14,15] and the references therein.
In the same spirit, it could be interesting to understand if our methods do apply in non-Euclidean setting for even more general nonlocal nonstandard growth equations, as the one recently considered in [10,40].
• Recently, mean value properties for solutions to fractional equations have been of great interest. It could be interesting to generalize such an investigation in a fractional non-Euclidean framework as the one considered in the present paper; we refer to the relevant results in [8,9] and the references therein.
• Moreover, to our knowledge, nothing is known about Harnack inequalities and more in general about the regularity for solutions to parabolic nonlocal integro-differential equations involving the nonlinear operators in (1.2).
• Finally, by starting from the estimates proven in the present paper, in [39] regularity results up to the boundary have been proven for very general boundary data, and for the related obstacle problem. As expected, a tail contribution naturally appears in those estimates in order to control the nonlocal contributions coming from far. Many subsequent related problems are still open, and not for free because of the possible degeneracy and singularity of L, as for instance boundary Harnack inequalities or Carleson estimates for the homogeneous case. The boundary Hölder estimates and the comparison results in the aforementioned paper [39] together with the Harnack estimates presented here could be a starting point for such a delicate investigation. Still for what concerns boundary Harnack inequalities, we also refer the reader to [42] for general strategy for equations with possibly unbounded right hand-side data.
To summarize. The results in the present paper seem to be the first ones concerning Harnack estimates for nonlinear nonlocal equations in the Heisenberg group. We prove that one can extend to the Heisenberg setting the strategy successfully applied in the fractional Euclidean case ( [17,18,30]), by attacking even a more general equation which applies to non-zero data and also to the case when sp ≥ Q. From another point of view, our results can be seen as the (nonlinear) nonlocal extension of the Heisenberg counterpart of the celebrated classical Harnack inequality ( [1,6,32]). Moreover, since we derive all our results for a general class of nonlinear integro-differential operators, via our approach by taking into account all the nonlocal tail contributions in a precise way, we obtain alternative proofs that are new even in the by-now classical case of the pure fractional sublaplacian operator (−∆ H n ) s ; also, in such a case, we are able to prove the robustness of the Harnack estimates with respect to s in the limit as s goes to 1.
We would guess that our estimates will be important in a forthcoming nonlinear nonlocal theory in the Heisenberg group.
The paper is organized as follows. In Section 2 below we set up notation, and we briefly recall our underlying geometrical structure, by also recalling the involved functional spaces, and providing a few remarks on the assumptions on the data. A few classical technical tools are also stated. In Section 3 we present very recent results for fractional equations in the Heisenberg group. In Section 4, we firstly carry out a suitable positivity expansion and some tail estimate. Then we complete the proof of the Harnack inequality with tail, and the weak Harnack inequality with tail, respectively. Section 5 is devoted to the asymptotic of the fractional subLaplacian operator, and the robustness of the Harnack inequalities in the linear case.

Preliminaries
In this section we state the general assumptions on the quantity we are dealing with. We keep these assumptions throughout the paper. Firstly, notice that we will follow the usual convention of denoting by c a general positive constant which will not necessarily be the same at different occurrences and which can also change from line to line. For the sake of readability, dependencies of the constants will be often omitted within the chains of estimates, therefore stated after the estimate. Relevant dependencies on parameters will be emphasized by using parentheses.

The Heisenberg-Weyl group
We start by very briefly recalling a few well-known facts about the Heisenberg group; see for instance [6] for a more exhaustive treatment.
We denote points in R 2n+1 by ξ := (z, t) = (x 1 , . . . , x n , y 1 , . . . , y n , t). For any ξ, ξ ′ ∈ R 2n+1 , consider the group multiplication • defined by For any λ > 0, the automorphism group (Φ λ ) λ>0 on R 2n+1 is defined by ξ → Φ λ (ξ) := (λx, λy, λ 2 t), and, as customary, Q ≡ 2n + 2 is the homogeneous dimension of R 2n+1 with respect to (Φ λ ) λ>0 , so that the Heisenberg-Weyl group H n : The Jacobian base of the Heisenberg Lie algebra h n of H n is given by Since [X j , X n+j ] = −4∂ t for every 1 ≤ j ≤ n, it plainly follows that rank Lie{X 1 , . . . , X 2n , T }(0, 0) = 2n + 1, so that H n is a Carnot group with the following stratification of the algebra We have now the following Definition 1 A homogeneous norm on H n is a continuous function (with respect to the Euclidean topology ) d o : H n → [0, +∞) such that: , for every λ > 0 and every ξ ∈ H n ; If d o is a homogeneous norm on H n , then the function Ψ defined by is a pseudometric on H n . We recall that the standard homogeneous norm | · | H n on H n is given by For any fixed ξ 0 ∈ H n and R > 0 we denote with B R (ξ 0 ) the ball with center ξ 0 and radius R defined by It is now worth noticing that for any homogeneous norm d o on H n one can prove the existence of a positive constant Λ such that As a consequence, in most of the estimates in the forthcoming proofs, one can simply take into account the pure homogeneous norm defined in (2.1) with no modifications at all.
In the analysis of the special case when the integro-differential operator (1.2) does reduce to the standard fractional subLaplacian, we will need to obtain fine estimates by taking into account the differentiability exponent s near 1, and thus several modifications with respect to the proof of similar estimates in the Euclidean framework are needed. In particular, in order to obtain the desired characterization of the asymptotic behaviour as s goes to 1 of the fractional subLaplacian, and consequently proving the consistency of our Harnack estimates with tail in the limit (see Section 5), we are able to overcome some difficulties coming from the non-Euclidean structure considered here by making use of a suitable MacLaurin-type expansion.
In the case of the Heisenberg group one can explicitly write the MacLaurin polynomial of Φ λ -degree 2; we have where the subgradient ∇ H n u is given by ∇ H n u(ξ) := X 1 u(ξ), . . . , X 2n u(ξ) , and D 2, * H n is the symmetrized horizontal Hessian matrix; that is, is the Taylor polynomial of H n -degree m centered at ξ 0 associated to u.

The fractional framework
We now introduce the natural fractional framework; we point out that our setup is in align with Section 2 in [33].
Let p > 1 and s ∈ (0, 1), the Gagliardo seminorm of a measurable function u : H n → R is given by 4) and the fractional Sobolev space W s,p (H n ) given by is equipped with the natural norm Similarly, given a domain Ω ⊂ H n , the fractional Sobolev space W s,p (Ω) is given by : endowed with the norm We denote by W s,p 0 (Ω) the closure of C ∞ 0 (Ω) in W s,p (H n ). We conclude this section by recalling the natural definition of weak solutions to the class of problem we are deal with; that is, where the datum f ≡ f (·, u) ∈ L ∞ loc (H n ) uniformly in Ω, the nonlocal boundary datum g belongs to W s,p (H n ), and the leading operator L is an integro-differential operator of differentiability exponent s ∈ (0, 1) and summability exponent p > 1 given by with d o being an homogeneous norm according to Definition 1.
For any g ∈ W s,p (H n ), consider the classes of functions We have the following Definition 4 A function u ∈ K − g (Ω) (respectively, K + g (Ω)) is a weak subsolution (respectively, supersolution) to (2.9) if for any nonnegative ψ ∈ W s,p 0 (Ω). A function u is a weak solution to problem (2.9) if it is both a weak sub-and supersolution.

Remark 1
The requirement on the boundary datum g to be in the whole W s,p (H n ) can be weakened by assuming a local fractional differentiability, namely g ∈ W s,p loc (Ω), in addition to the boundedness of its nonlocal tail; i. e., Tail(g; ξ 0 , R) < ∞, for some ξ 0 ∈ H n and some R > 0. This is not restrictive, and it does not bring modifications in the rest of the paper. For further details on the related "Tail space", we refer the interested reader to papers [28,29].

Classical technical tools
As in the classical variational approach to local Harnack estimates, the following wellknown iteration lemmata are needed.
Lemma 1 Let β > 0 and let {A j } j∈N be a sequence of real positive numbers such that where c 1 , c 2 , θ and ζ < 1 are nonnegative constant. Then, there exists a constant c depending only on θ and ζ, such that for every ρ, R, T 0 ≤ ρ < R ≤ T 1 , we have

Some recent results on fractional operators in the Heisenberg group
Similarly to what happens in the Euclidean case, a fractional Sobolev embedding can be proved in the non-Euclidean setting of the Heisenberg group. Indeed, the following result holds true, Theorem 5 Let p > 1 and s ∈ (0, 1) such that sp < Q. For any measurable compactly supported function u : H n → R there exists a positive constant c = c(n, p, s) such that where p * := Qp Q−sp is the critical Sobolev exponent. For the proof of the previous statement we refer to Theorem 2.5 in [26], where the authors extend the same strategy used in [19,37] in order to prove the fractional Sobolev embedding in the Euclidean setting.
Before proving the weak Harnack inequality, namely Theorem 2, we also need to recall a boundedness estimate and a Caccioppoli-type one for the weak sub-and super-solution to (1.1).

Theorem 7 (Caccioppoli estimates with tail) [Theorem 1.3 in
where w := (u + d) p−q p , and the constant c depends only on n, p, q, φ L ∞ (supp φ) and the structural constant Λ defined in (2.2).

Proofs of the Harnack estimates with tail
For the sake of readability, from now on we adopt the following notation, (4.1)

Expansion of positivity
In this section we prove a careful estimate of the weak supersolutions to (1.1), by generalizing the original strategy applied in the local framework as well as that in the fractional Euclidean one. Clearly, we have to take into account the needed modifications to handle the difficulties given in the fractional Heisenberg framework; also, the presence of the nonhomogeneous datum f would require further care when some iterative argument will be called.
Proof. We firstly notice that with no loss of generality one can suppose Tail(u − ; ξ 0 , R) > 0, otherwise (4.3) plainly follows from (4.2) by choosing the constantc large enough. Take a smooth function φ ∈ C ∞ 0 (B 7r ) such that Now, choose ψ :=ũ 1−p φ p as a test function in Definition 4 by making use of the fact that u is a supersolution. It follows where we have also used the definition ofũ and the translation invariant property of the fractional seminorms. We start by estimating the integral contribution in the left-hand side. Thanks to the hypothesis on f and in view of the definition of ψ, we have get Now, we focus on the terms on the right-hand side. We can treat the first integral I 1 as in the proof of the Logarithmic Lemma 1.4 in [33], so that We now proceed to estimating the second integral I 2 in (4.5), in turn obtaining an estimate for I 3 as well. We split I 2 as follows, The contribution in I 2,1 can be estimated as follows where we have used that, for any ξ ∈ B 7r and any η ∈ H n B 8r , For I 2,2 , notice thatũ ≥ 0 in B 7r . Then,ũ(ξ) −ũ(η) ≤ũ(ξ), for every η ∈ H n B 8r ∩ ũ(η) ≥ 0 and every ξ ∈ B 7r ; we get Thus, by combining (4.6) with (4.7) and (4.8), we arrive at Since v is a truncation of log(k + d) − log(ũ), the estimate in (4.9) yields Therefore, the inequality in (4.2) yields Also, notice that where we also used (4.10) and the Hölder inequality. From all the previous estimates we finally arrive at so that the desired inequality plainly follows by inserting the definition of d given in (4.4).
We are now in the position to refine the estimates in order to prove the main result of this section; i. e., Lemma 4 Let s ∈ (0, 1) and p ∈ (1, ∞) and let u ∈ W s,p (H n ) be a weak supersolution to Let k ≥ 0 and suppose that there exists σ ∈ (0, 1] such that for some r satisfying 0 < 6r < R. Then, there exists a constant δ ∈ (0, 1/4) depending on n, s, p, σ and Λ, such that Proof. We immediately notice that in the case when k = 0, the inequality (4.11) does trivially hold, since u ≥ 0 in B R . Also, with no loss of generality, we can assume that for Now, for any r ≤ ρ ≤ 6r, take a smooth function φ ∈ C ∞ 0 (B ρ ) such that 0 ≤ φ ≤ 1, and consider the test function ψ := w − φ p , where w − := (ℓ − u) + , for any ℓ ∈ (δk, 2δk). Testing Definition 4 with such a smooth function ψ yields (4.13) We begin to estimate the term on the left-hand side. As done in our proof of Lemma 3, we obtain that Now we focus on the right-hand side of (4.13). It is convenient to split J 2 as follows Now, notice that which yields For J 2,2 , since u ≥ 0 in B ρ , we can write By reasoning as above for the integral J 3 , we finally arrive at It remains to estimate the contribution J 1 , and for this one can proceed as seen in the Euclidean setting (see Theorem 1.4 in [18]); it follows Combining all the above estimates, we get At this point an iteration argument is needed. To this aim, define ℓ ≡ ℓ j := δk + 2 −j−1 δk, ρ ≡ ρ j := 4r + 2 1−j r, andρ j := ρ j+1 + ρ j 2 , ∀j = 0, 1, ...
Note that both 4r < ρ j ,ρ j < 6r and Moreover, in view of (4.12), we have We are now in the position to apply the estimate in Lemma 3; we arrive at Denote by B j := B ρj (ξ 0 ) and let φ j ∈ C ∞ 0 (Bρ j ) be such that We have where in the last inequality we have used the Sobolev embedding with p * = Qp/(Q − sp). Let us estimate (4.17) with the aid of (4.14). Firstly, by the particular choice of φ, we have where, proceeding as in the proof of the Logarithmic Lemma 1.4 in [33], we have that Now, notice that, for any η ∈ H n B j and any ξ ∈ supp φ j ⊂ Bρ j , it holds Thus, where we have also used the fact that u ≥ 0 in B R , the estimate in (4.12), and the fact that δk < ℓ j .
From (4.14), (4.17), and (4.18), we arrive at We are finally in the position to apply the classic iteration Lemma 1. We set the previous estimate can be rewritten as follows and using the estimate in (4.16), we arrive at

from (4.16) we get that
Then which implies inf B4r u ≥ δk and hence (4.11) when sp = Q.
The case when sp > Q can be deduced as in the latter, without relevant modifications, choosing the parameter ε > (s−Q/p). With such a choice in hand, we can use the Sobolev embedding for W sε,p in Theorem 5 and the desired result plainly follows as in the previous case when sp = Q.

Proof of Theorem 1
In order to derive the Harnack inequalities with tail, we firstly need the estimate (4.20) below, which is a straightforward consequence of the refined positivity expansion proven in the previous section, together with the classical Krylov-Safonov covering lemma (whose proof can be found for instance in [27,Lemma 7.2]), which can be adjusted to our framework thanks to the role of the nonlocal tail, as shown in the Euclidean framework in the proof of [17,Lemma 4.1].
In the next lemma, we prove that the tail of the positive part of the weak solutions to (1.1) can be controlled in a precise way.
Proof. Set k := sup Br u and choose a cut-off function φ ∈ C ∞ 0 (B r ) such that 0 ≤ φ ≤ 1, φ ≡ 1 on B r/2 and |∇ H n φ| ≤ 8/r. Take now the test function ψ := (u − 2k)φ p . We have (4.22) The last two integral in the identity above can be estimated as in the proof of Lemma 4.2 in [17]; we have For what concerns the contribution H 1 in (4.22), we have where we argued as in the proof of Lemma 1.4 in [33]. The contribution given by the datum f can be easily estimated as follows, Armed with the tail estimate in Lemma 6, and the interpolative inequality given by Theorem 6, we are ready to complete the proof of the Harnack inequality with tail in (1.6). The strategy does generalize that successfully applied in [17] in the analysis of the homogeneous case in the Euclidean framework.
We firstly deal with the case when sp < Q. In such a case one can apply Theorem 5 to the function wφ, with w :=ũ where we have also used the estimate below Collecting the estimates (4.27) and (4.28) with the Caccioppoli inequality of Theorem 7, we obtain where we have used the fact that Now, we choose d as in the proof of Lemma 3; see (4.4) there. It follows where c depends only on n, p, s and the structural constant Λ defined in (2.2). Let t = (p − q)Q/(Q − sp) for any q ∈ (1, p). Thanks to a standard finite Moser iteration, the inequality in (4.30) becomes for any 0 < t ′ < t < (p − 1)Q/(Q − sp). Since 6r < R, we can apply Lemma 5 with α = t ′ there; it follows which provides the desired inequality, up to relabelling r.
We investigate now the case when sp = Q. Fix 0 < ε < s, and set s ε := s − ε. Since s ε p < sp = Q, we can make use of the Sobolev inequality for the function wφ, to get where we have obtained the last inequality by following by the same argument used at the end of the proof of Lemma 4. Now, we choose d as in (4.4), so that the inequality in (4.32) becomes Thus, recalling the definition of w and that of the cut-off function φ, we have that Set t = (p − q)s/ε, for any q ∈ (1, p); a standard application of the finite Moser iteration yields for any 0 < t ′ < t < (p − 1)s/ε.
Finally, we can conclude as in the proof in the case when sp < Q; that is, it suffices to apply Lemma 5, with α = t ′ there, in order to get which provides the desired inequality, up to relabelling r.
The case when sp > Q can be deduced as in the latter, without relevant modifications, choosing the parameter ε > (s − Q/p). With such a choice in hand, it plainly follows that s ε p < Q, and thus one can use the Sobolev embedding in Theorem 5 and proceed as done in the limit case when sp = Q. where P 2 (u, ξ) is the Taylor polynomial of H n -degree 2 associated to u and centered at ξ presented in Section 2.1. Also, by the very definition of Taylor polynomial, it follows 1 r 2s−1 dr = c 2 (n, s) 2(1 − s) X 2 i u(ξ).

Robustness of the nonlocal Harnack estimates
The proofs of Theorem 3 and Theorem 4 can be plainly deduced from the ones in Section 4 for the subcritical case, by taking there p = 2, f ≡ 0, and the Korani-Folland norm in place of the generic homogeneous norm d o . Below we stated the related needed lemmata, by indicating only the modifications in the estimates where a special care on the involved quantities is needed in order to successfully obtaining the desired robustness in the limit as s goes to 1.
Firstly, we need the related positivity expansion, which can be condensed in the following two lemmata.
The proof will then follow with no further modifications at all.
As well as in the proof in the general nonlinear framework presented in Section 4, we can obtain for the pure subLaplacian case the analogue of the estimate in Lemma 5 and that of the Tail control estimate stated in Lemma 6. For the sake of the reader, we prefer to restate these results by stressing the novelty of the dependance on s here. No modifications in the related proofs are essentially needed, thanks to the results obtained in Lemma 7 and Lemma 8.