Continuity of solutions to a nonlinear fractional diffusion equation

We study a parabolic equation for the fractional $p-$Laplacian of order $s$, for $p\ge 2$ and $0<s<1$. We provide space-time H\"older estimates for weak solutions, with explicit exponents. The proofs are based on iterated discrete differentiation of the equation in the spirit of J. Moser.

1. Introduction 1.1. The problem. In this paper, we study the regularity of weak solutions to the nonlinear and nonlocal parabolic equation where 2 ≤ p < ∞, 0 < s < 1 and (−∆ p ) s is the fractional p-Laplacian of order s, i.e. the operator formally defined by is still a solution. By combining the last two facts, we also get that u λ,µ = µ u λ x, µ p−2 λ s p t , for λ, µ > 0, still solves (1.1). We will make a repeated use of this simple fact.
In this paper, we are concerned with the Hölder regularity for weak solutions of (1.1). More precisely, we prove that local weak solutions (see Definition  , if s ≥ p − 1 p .
To the best of our knowledge, our result is the first pointwise continuity estimate for solutions of this equation.
1.2. Background and recent developments. In recent years there has been a surge of interest around the operator (1.2), after its introduction in [18]. In particular, equation (1.1) has been studied in [1,25,26,31,33] and [34]. References [26], [25] and [33] dealt with existence and uniqueness of solutions, together with their long time asymptotic behaviour. Similar properties for (1.1) with a general right-hand side in place of 0 are studied in [1]. In [34], some regularity of the semigroup operator generated by (−∆ p ) s was studied. In [31], the local boundedness of weak solutions of (1.1) is proved. Recently, in [15], a weaker pointwise regularity result was obtained for viscosity solutions of the doubly nonlinear equation (1.3) |∂ t u| p−2 ∂ t u + (−∆ p ) s u = 0, by using completely different methods. This equation and its large time behavior is related to the eigenvalue problem for the fractional p−Laplacian. A crucial difference between this equation and (1.1), is that the former is homogeneous, a feature which is not shared by our equation, as already observed in Remark 1.1. Moreover, the nonlinearity in the time derivative in (1.3) makes the notion of weak solutions less useful. It is not clear whether the methods in [15] can be adapted to the present situation or not.
In the linear or non-degenerate case, corresponding to p = 2, the literature on regularity is vast. We mention only a fraction of it, namely [7,22,23,29,30] and [32]. However, we point out that neither of these results apply to our setting.
The stationary version of (1.1), i.e., (−∆ p ) s u = 0, has attracted a lot of attention, as well. The regularity of solutions has been studied for instance in [3,4,6,12,13,17,16,19,20,21,24,27] and [34]. In particular, the regularity result proved in the present paper can be seen as the parabolic version of that obtained by the first two authors and Schikorra in [4] for the stationary equation.
The local counterpart of (1.1) is the parabolic equation for the p−Laplacian This has been intensively studied and only in the last decades has its theory reached a rather complete state. We refer to [10] and [11] for a complete account on the regularity results for this equation and some of its generalizations. At present, the best local regularity known is spatial C 1,α −regularity for some α > 0 (see [10, Chapter IX]) and C 0,1/2 −regularity in time (see [2,Theorem 2.3]). None of these exponents is known to be sharp. However, due to the explicit solution it is clear that solutions cannot be better than C 1,1/(p−1) in space.

Main result.
The main result of our paper is the following Hölder regularity for local weak solutions of (1.1). Here, we use the following notation for parabolic cylinders with B r (x 0 ) denoting the N −dimensional ball of radius r centered at the point x 0 . For the precise definition of local weak solution, as well as of the spaces C δ x,loc (Ω × I) and C γ t,loc (Ω × I), we refer the reader to Sections 3.1 and 2.3, respectively , if s ≥ p − 1 p .

Remark 1.3 (Comment on the time regularity)
. The regularity in time is almost sharp for s p ≤ (p − 1). Indeed, our result in this case gives Hölder continuity for any exponent less than 1. The following example from [23] shows that solutions are not C 1 in time in general. Let v(x, t) = 0, if t < −1/2, C (1/2 + t) + 1 B3\B2 (x), if t ≥ −1/2, where C = 0 is chosen so that v is a local weak subsolution (see Definition 3.1)  u ∈ L ∞ loc (I; L ∞ loc (Ω)), and (1.8) u ∈ L ∞ loc (I; L p−1 s p (R N )), where the tail space L p−1 s p (R N ) is defined by L p−1 s p (R N ) = u ∈ L p−1 loc (R N ) :ˆR N |u| p−1 1 + |x| N +s p dx < +∞ .
We point out that by [31,Lemma 2.6], condition (1.8) is a natural one in order to guarantee the local boundedness (1.7). However, it is not known apriori if the quantity (1.8) is finite whenever u is a weak solution. Indeed, even if u solves the initial boundary value problem    ∂ t u + (−∆ p ) s u = 0, in Ω × I, u = g, on (R N \ Ω) × I, u = u 0 , on Ω × {t = t 0 }, with the boundary data g satisfying g ∈ L ∞ loc (I; L p−1 s p (R N )), it is not evident that this is sufficient to entail (1.8). For this reason, and to not overburden an already technical proof, we have chosen to assume the simpler condition (1.4). For completeness, in Appendix A we give some sufficient conditions assuring that our weak solutions verify (1.4), see Corollary A.5 below.
1.4. Main ideas of the paper. The idea we use to prove Theorem 1.2 is very similar to the method employed in [4] for the elliptic case: we differentiate equation (1.1) in a discrete sense and then test the differentiated equation against functions of the form For suitable choices of ϑ > 0 and β ≥ 1, this gives an integrability gain (see Proposition 4.1) of the form for −1/2 ≤ T ≤ 0 and an arbitrary µ > 0. By first fixing T = 0 and ignoring the second term in the left-hand side of (1.9), this can be iterated finitely many times in order to obtain δ h u |h| s ∈ L q ([−1/2, 0]; L q loc ), for every q < ∞, uniformly in |h| ≪ 1.
We can then use the second term in the left-hand side of (1.9), so to get Thus, by using a Morrey-type embedding result, we can conclude that u ∈ C δ loc spatially for any 0 < δ < s. After this, we prove Proposition 5.1, which comprises a refined version of the scheme (1.9). Namely, an estimate of the form dt.
Also (1.10) can be iterated, where now both the differentiability ϑ and the integrability β change. The result is that u ∈ C δ loc spatially, for every 0 < δ < Θ(s, p), again uniformly in time. The last part of the paper, where we obtain the regularity in time, is quite standard for this kind of diffusion equations (see for example [8, page 118]). It amounts to using the already established spatial regularity and the information given by the equation. However, due to the fractional character of the spatial part of our equation, some care is needed in order to properly handle the time regularity. In particular, we have to treat the cases separately. This is done in Proposition 6.2 and it yields the γ−Hölder continuity in time for any given that the solution is δ−Hölder continuous in the x variable. In particular, by the possible choice of δ, this yields that we may choose any γ < Γ(s, p), where the latter exponent is the one defined in (1.5).
1.5. Plan of the paper. The plan of the paper is as follows. In Section 2, we introduce the expedient spaces and notation used in this paper. In Section 3, we define local weak solutions and justify that we can insert certain test functions in the differentiated equation (see Lemma 3.3 below). This is followed by Section 4, where we prove that weak solutions are almost s−Hölder continuous in the spatial variable. In Section 5, we improve this result up to the exponent Θ(s, p) defined in (1.5). This result is then used in Section 6, where we prove the corresponding Hölder regularity in time. Finally, in Section 7 we prove our main theorem. The paper is complemented by an appendix, where for completeness we prove existence and uniqueness of weak solutions for the initial boundary value problem related to our equation. A comparison principle is also presented.

Preliminaries
2.1. Notation. We denote by B r (x 0 ) the N −dimensional open ball of radius r centered at the point x 0 . The ball of radius r centered at the origin is denoted by B r . Its Lebesgue measure is given by We use the following notation for the parabolic cylinder Again, when x 0 = 0 and t 0 = 0, we simply write Q R,r .
Let 1 < p < ∞, we denote by p ′ = p/(p − 1) the dual exponent of p. For every β > 1, we define the monotone function For a function ψ : R N × R → R and a vector h ∈ R N , we define . It is not difficult to see that the following discrete Leibniz rule holds 2.2. Sobolev spaces. We now recall the main notations and definitions for the relevant fractional Sobolev-type spaces throughout the paper.
Let 1 ≤ q < ∞ and let ψ ∈ L q (R N ), for 0 < β ≤ 1 we set , and for 0 < β < 2 We then introduce the two Besov-type spaces We also need the Sobolev-Slobodeckiȋ space We endow these spaces with the norms A few times we will also work with the space W β,q (Ω) for a subset Ω ⊂ R N , The space W β,q 0 (Ω) is the subspace of W β,q (R N ) consisting of functions that are identically zero in the complement of Ω.
2.3. Parabolic Banach spaces. Let I ⊂ R be an interval and let V be a separable, reflexive Banach space, endowed with a norm · V . We denote by V * its topological dual space. Let us suppose that v is a mapping such that for almost every t ∈ I, v(t) belongs to V . If the function t → v(t) V is measurable on I and 1 ≤ p ≤ ∞, then v is an element of the Banach space L p (I; V ) if and only if By [28,Theorem 1.5], the dual space of L p (I; V ) can be characterized according to We write v ∈ C(I; V ) if the mapping t → v(t) is continuous with respect to the norm on V . We say that u is locally α−Hölder continuous in space (respectively, locally β−Hölder continuous in time) on Ω × I and write That is, if u ∈ C α x (K × J) (respectively, u ∈ C β t (K × J)).

Tail spaces.
We recall the definition of tail space q ≥ 1 and α > 0, which is endowed with the norm For every x 0 ∈ R N , R > 0 and u ∈ L q α (R N ), the following quantity plays an important role in regularity estimates for solutions of fractional problems. We recall the following result, see for example [  • we have the continuous inclusion L m α (R N ) ⊂ L q α (R N ); • for every 0 < r < R and x 0 ∈ R N we have Tail q,α (u; x 0 , R) q .

Weak formulation
3.1. Local weak solutions. In the following, we assume that Ω ⊂ R N is a bounded open set in R N .
We say that u is a local weak solution to the equation if for any closed interval J = [T 0 , T 1 ] ⊂ I, the function u is such that u ∈ L p (J; W s,p loc (Ω)) ∩ L p−1 (J; L p−1 s p (R N )) ∩ C(J; L 2 loc (Ω)), and it satisfies for any φ ∈ L p (J; W s,p (Ω)) ∩ C 1 (J; L 2 (Ω)) which has spatial support compactly contained in Ω. In equation (3.2), the symbol ·, · stands for the duality pairing between W s,p (Ω) and its dual space (W s,p (Ω)) * . We also say that u is a local weak subsolution if instead of the equality above, we have the ≤ sign, for any non-negative φ as above. A local weak supersolution is defined similarly.
Remark 3.2. We observe that L ∞ (R N ) ⊂ L p−1 s p (R N ). This in turn implies that L ∞ (J; L ∞ (R N )) ⊂ L p−1 (J; L p−1 s p (R N )). We will use this fact repeatedly.

3.2.
Regularization of test functions. Let ζ : R → R be a nonnegative, even smooth function with compact support in (−1/2, 1/2), satisfying´R ζ(τ ) dτ = 1. If f ∈ L 1 ((a, b)), we define the convolution The following result justifies that we may take powers of differential quotients of a solution, as test functions. This is needed in the sequel. In the rest of the paper, we will use the abbreviated notation dµ(x, y) = dx dy |x − y| N +s p .
Then, for any locally Lipschitz function F : R → R and any h ∈ R N such that 0 < h < dist (supp η, ∂B 2 )/4, we haveˆT Proof. We take J = [T 0 , T 1 ] ⊂ (−1, 0) and φ ∈ L p (J; W s,p (B 2 )) ∩ C 1 (J; L 2 (B 2 )), whose spatial support is compactly contained in B 2 . We want to use the time-regularization φ ε as test function in (3.1). For this, we take Then, we preliminary observe that from elementary properties of convolutions, Fubini's Theorem and integration by parts, we have For simplicity, we have set Thus from (3.2) it follows that for 0 < ε < ε 0 Before proceeding further, we observe that by using an integration by parts, the term Σ(ε) can be rewritten as where we also used that ζ has compact support in (−1/2, 1/2). By further using a suitable change of variables, we can also write and then changing variables, we get (3.7) The quantity Σ h (ε) is defined as in (3.6), with u h in place of u. We subtract (3.5) from (3.7), so to get ), whose spatial support is compactly contained in B 2 . We take F as in the statement and use (3.8) with the test function and η and τ are as in the statement. By observing that Observe that we used the properties of τ ε . In order to deal with the integral containing the time derivative of δ h u ε , we first observe that . Thus we can use integration by parts, which yieldŝ By inserting this into (3.9), we get (3.10) We recall that this is valid for 0 < |h| < h 0 4 and 0 < ε < ε 0 .
Before taking the limit as ε goes to 0, we first observe that for t ∈ This shows that we have the uniform L ∞ estimate Finally, we pass to the limit in (3.10) as ε goes to 0. We start from the right-hand side: by using the local Lipschitz regularity of F and (3.11), we have where C > 0 does not depend on ε. Thus, by using that η has compact support in B 2 and 0 < |h| < h 0 /4, we get from the last estimate (after a change of variable) The constant C is still independent of 0 < ε < ε 0 . If we now use that u ∈ C((−2, 0]; L 2 loc (B 2 )), we get that the last quantity converges to 0, as ε goes to 0.
For the termˆB we proceed similarly as above. We observe that We can now use again that u ∈ C((−2, 0]; L 2 loc (B 2 )) and obtain that the last quantity converges to 0, as ε goes to 0.
As for the term we can proceed exactly as before, we omit the details. In a similar fashion, we can also show that This is still similar to the previous limits. It is sufficient to use the expression (3.6), the uniform L ∞ estimate (3.11) and the fact u ∈ C((−2, 0]; L 2 loc (B 2 )), in order to apply the Lebesgue Dominated Convergence Theorem. Finally, the convergence of the double integral requires quite lengthy computations and thus we prefer to postpone them to Appendix B below.
Remark 3.4. We observe that the global L ∞ bound on the weak solution is not needed in the previous result. It is sufficient to know that the weak solution is locally bounded. We refer to [32, Theorem 1.1] for local boundedness of weak solutions.

Spatial almost C s -regularity
The following result is an integrability gain for the discrete derivative of order s of a local weak solution. This is the parabolic version of [4, Proposition 4.1], to which we refer for all the missing details.
and that, for some q ≥ p and 0 < h 0 < 1/10, we havê Proof. We divide the proof into seven steps.
Step 1: Discrete differentiation of the equation. We take for the moment T 1 < 0, then we will show at the end of the proof how to include the case T 1 = 0. We already introduced the notation For notational simplicity, we also set r = R − 4 h 0 . Let β ≥ 2 and ϑ ∈ R be such that 0 < 1 + ϑ β < β, and use (3.4) for 0 < |h| < h 0 , where: • • η is a non-negative standard Lipschitz cut-off function supported in B (R+r)/2 , such that • τ is a smooth function such that 0 ≤ τ ≤ 1 and Here µ is as in the statement, i.e. any positive number such that µ < T 1 − T 0 . Note that the assumptions on η imply After dividing by |h| 1+ϑ β , we obtain from Lemma 3.3, The triple integral is now divided into three pieces: where and where we used that η vanishes identically outside B (R+r)/2 . We also suppressed the t−dependence, for notational simplicity. We also have the term in the right-hand side By proceeding exactly as in Step 1 of the proof of [4, Proposition 4.1], we get the following lower bound for I 1 (t) where c = c(p, β) > 0 and C = C(p, β) > 0. We use that |h| 1+ϑβ η p dx = I 4 and the estimate for I 1 (t). This entails that where we set I 11 =´T 1 T0 I 11 τ dt, I 12 =´T 1 T0 I 12 τ dt and and Step 2: Estimates of the local terms I 11 and I 12 . Here we can follow the same computations as in Step 2 of the proof of [4, Proposition 4.1], so to get Step 3: Estimates of the nonlocal terms I 2 and I 3 . Both nonlocal terms I 2 and I 3 can be treated in the same way. We only estimate I 2 for simplicity. We can use that |u| ≤ 1 on R N × [−1, 0] to infer that Hence, we obtain by Young's inequality. Here C = C(h 0 , N, s, p, q, β) > 0 as before.
Step 4: Estimates of I 4 . By using that |u| ≤ 1 in R N × [−1, 0] and the properties of τ , we get In the last inequality we further used Young's inequality. By inserting the estimates (4.6) and (4.7) in (4.5), using that τ is non-negative and such that τ = 1 on This is the parabolic counterpart of [4, equation (4.10)]. Observe that the constant C now depends on 1/µ, as well.
Step 6: Conclusion for T 1 < 0. As in the final step of the step of [4, Proposition 4.1], we now fix where q ≥ p is as in the statement. These choices assure that where C = C(N, h 0 , p, q, s) > 0. Up to a suitable modification of the constant C, we obtain in particular as desired. Observe that we used that r = R − 4 h 0 .
Step 7: Conclusion for T 1 = 0. In this case, the previous proof does not directly work because it relies on Lemma 3.3, which needed T 1 < 0. However, the constant C in (4.1) does not depend on T 1 , we can thus use a limit argument. By assumption, we have that for some q ≥ p and 0 We then observe that by the Dominated Convergence Theorem. As for the second term on the left-hand side, we know by definition of local weak solution that is a continuous function on (−2, 0], with values in L 2 (B R−4 h0 ), for every fixed 0 < |h| < h 0 . Thus This in turn implies that 1 (4.14) lim inf , for every 0 < |h| < h 0 . By using (4.13) and (4.14) in (4.12), we get the desired conclusion for T 1 = 0, as well.
As in [4,Theorem 4.2], by iterating the previous result, we can obtain the following regularity estimate. 1 We use the following standard fact: if {fn} n∈N converges to f in L α (E), then for any β = α.
such that u ∈ L ∞ loc (I; L ∞ (R N )). Then u ∈ C δ x,loc (Ω × I) for every 0 < δ < s. More precisely, for every 0 < δ < s, R > 0 and every (x 0 , T 0 ) such that there exists a constant C = C(N, s, p, δ) > 0 such that Proof. We assume for simplicity that x 0 = 0 and T 0 = 0, then we set By taking into account the scaling properties of our equation (see Remark 1.1), the function u R,α is a local weak solution of u t + (−∆ p ) s u = 0, in B 2 × (−2, 0], and satisfies We will prove that u R,α satifies the estimate for C = C(N, s, p, δ) > 0 independent of α. By scaling back, this would give which is the desired result. In what follows, we suppress the subscript R, α and simply write u in place of u R,α , in order not to overburden the presentation.
We fix 0 < δ < s and choose i ∞ ∈ N \ {0} such that Then we define the sequence of exponents We define also We note that By applying Proposition 4.1 (ignoring the second term in the left-hand side of (4.1)) with 2 and R = R i and q = q i = p + i, for i = 0, . . . , i ∞ − 1, and observing that R i − 4 h 0 = R i+1 + 4 h 0 , we obtain the iterative scheme of inequalities: Here C = C(N, δ, p, s) > 0 as always. We note that by using the relation and then appealing to [3, Proposition 2.6], we havê where we also have used the assumptions (4.16) on u. Hence, the iterative scheme of inequalities leads us tô (N, δ, p, s). 2 We observe that by construction we have Thus these choices are admissible in Proposition 4.1.
It is now time to exploit the full power of Proposition 4.1: we apply it once more, with We obtain (ignoring the first term in the left-hand side of (4.1), this time) Since this is valid for every −1/2 ≤ T 1 ≤ 0, this in turn implies that (4.18) sup (N, δ, p, s).
In particular, we have for all |h| > 0 |δ h χ| We also recall that Hence, for 0 < |h| < h 0 and any t ∈ [−5/8, 0] by (4.18). Finally, by noting that thanks to the choice of i ∞ we have we may invoke the Morrey-type embedding of [4, Theorem 2.8] with Thus we obtain for any t ∈ [−1/2, 0], where we used (4.19). This concludes the proof.
with C now depending on σ as well (and blowing-up as σ ր 7/8). Indeed, if σ ≤ 1/2 then this is immediate. If 1/2 < σ < 7/8, then we can cover Q σR,σR s p (x 0 , T 0 ) with a finite number of cylinders By using (4.15) on each of these cylinders, we get By taking the supremum over 1 ≤ i ≤ k and 1 ≤ j ≤ m, we get the desired conclusion.

Improved spatial Hölder regularity
Once we know that solutions are locally spatially δ−Hölder continuous for any 0 < δ < s, we can obtain the following improvement of Proposition 4.1. The latter provided a recursive gain of integrability. In contrast, the next result provides a gain which is interlinked between differentiability and integrability. Assume further that for some 0 < h 0 < 1/10 and ϑ < 1, β ≥ 2 such that (1 + ϑ β)/β < 1, we havê for a radius 4 h 0 < R ≤ 1 − 5 h 0 and two time instants Proof. This is analogous to the proof of [4, Proposition 5.1]. As above, we will refer to [4] for the main computations and only list the major changes. We first notice that it sufficient to prove (5.1) for T 1 < 0, with a constant independent of T 1 . Then the same argument of Step 7 in Proposition 4.1 will be enough to handle the case T 1 = 0, as well.
We go back to the estimates in the proof of Proposition 4.1. The acquired knowledge on the spatial regularity of u permits to improve the estimate on the term I 11 (t) defined in ( [u(·, t)] C s−ε (B R+h 0 ) ≤ C(N, h 0 , p, s).
Using this together with the assumed regularity of η, we have for (x, y) ∈ B R+h0 and t ∈ [T 0 , As usual, we are suppressing the time dependence. Thanks to the choice of ε, the last exponent is strictly larger than −N and we may concludê for any x ∈ B R . Therefore, by suppressing as before the t−dependence for simplicity, we have the estimate for some C = C(N, h 0 , p, s) > 0.
As for I 12 , by going back to its definition (4.4) and using the properties of the cut-off function η, we get where we used the local L ∞ bound on u, as above. In addition, from the first inequality in (4.6) together with the properties of the cut-off function τ , we have Combining these new estimates with (4.7) and (4.2), we can reproduce the last part of [4, Proposition 5.1] and arrive atˆT for some C = C(N, h 0 , p, s, β) > 0. By appealing again to [4, Lemma 2.6] and using that 1 + ϑ β β < 1, we may replace the first order differential quotients in the right-hand side by second order ones. This leads tô for some C = C(N, h 0 , p, s, β) > 0. By recalling again that r = R − 4 h 0 , we eventually conclude the proof.
We are now ready to prove the claimed Hölder regularity in space.
By induction, we see that {ϑ i } i∈N is explicitely given by the increasing sequence and thus lim The proof is now split into two different cases.

Define also
We note that By applying 3 Proposition 5.1 (ignoring the second term of the left-hand side of (5.1)) with we obtain the iterative scheme of inequalities: Here C = C(N, p, s, δ) > 0 as always. As in (4.17) we havê (N, δ, s, p). (N, δ, p, s). Now we apply Proposition 5.1 once more, this time with

Hence, the previous iterative scheme of inequalities implieŝ
We obtain (now ignoring the first term in the left-hand side of (5.1)) Since this is valid for every −1/2 ≤ T 1 ≤ 0, we obtain (N, δ, p, s).
From here, we may repeat the arguments at the end of the proof of Theorem 4.2 (see (4.19)) and use the Morrey-type embedding of [4, Theorem 2.8], with to obtain sup which concludes the proof in this case.
Once we land here, as before we can repeat the arguments at the end of the proof of Theorem 4.2 and use the Morrey-type embedding, this time with

Regularity in time
In this section, we prove Hölder regularity in time using the previously obtained regularity in space. This approach uses energy estimates to control the growth of local integrals which yields a Campanato-type estimate. We will use the notation When the center x 0 is clear from the context, we often simply write u R . For u ∈ L 1 (Q R,r (x 0 , t 0 )), we set Again, when the center (x 0 , t 0 ) is clear from the context, we simply write u R,r .
The following simple Poincaré-type inequality will be useful.
Lemma 6.1. Let 1 ≤ p < ∞ and let B r = B r (x 0 ). Suppose that u ∈ W s,p (B r ), then for any nonnegative η ∈ C ∞ 0 (B r ) such that η r = 1, there holds Proof. By using the fact that´B r η dx = |B r | and Jensen's inequality, we obtain This concludes the proof.

Proposition 6.2. Suppose that u is a local weak solution of
where Θ(s, p) is the exponent defined in (1.5). Then there is a constant C = C(N, s, p, K δ , δ) > 0 such that In particular, u ∈ C γ t (Q 1 4 , 1 4 ) for any γ < Γ(s, p), where Γ(s, p) is the exponent defined in (1.5).
Proof. We take (x 0 , t 0 ) ∈ Q 1/4,1/4 and choose Consider the parabolic cylinder Observe that by construction we have Let η ∈ C ∞ 0 (B r/2 (x 0 )) be a non-negative cut-off function, such that for some constant C = C( η L ∞ (B r/2 (x0)) , N ) > 0. Observe that, thanks to the condition on its average, we have Thus the constant appearing in (6.1) will only depend on N, s and p.
We now write where we have set u(y, t) η(y) dy.
We first note that Thus it suffices to estimate A 1 and A 3 . In view of Lemma 6.1, we have , for some C = C(N, s, p) > 0. Recalling that δ > s and using the spatial Hölder continuity of u, we find that Indeed, by observing that for every where we used spherical coordinates to compute the last integral. Observe that the width θ of the time does not come into play here. We now turn to A 3 and first note that If T 0 , T 1 ∈ (t 0 − θ, t 0 ] with T 0 < T 1 , we use the weak formulation (3.2) with φ(x, t) = η(x), to obtain T0¨Br (x0)×Br (x0) J p (u(x, τ ) − u(y, τ )) (η(x) − η(y)) dµ(x, y) dτ In order to control J 2 , we claim that for t ∈ [−1/2, 0], x ∈ B r (x 0 ) and y ∈ R N , Indeed, if y ∈ B 1/2 this follows directly from the assumption. On the other hand, if y ∈ R N \ B 1/2 , then by construction Additionally, if y ∈ R N \ B r (x 0 ) and x ∈ B r/2 (x 0 ), we have Thus, by using this and (6.7), we get for some C = C(δ, N, s, p, K δ ) > 0. Observe that we used that δ (p − 1) − s p < 0, in order to assure that the integral on R N \ B r (x 0 ) converges. As for J 1 , we have for δ > s for some C = C(N, s, p, δ) > 0. By recalling (6.5) and using the estimates on J 1 and J 2 in (6.6), we have thus shown that A 3 ≤ C K p−1 δ θ r δ (p−1)−s p , for some C = C(δ, N, s, p) > 0.
Hence, by also using (6.4) and (6.3), we get (6.8) We now have to distinguish two cases: • Case s p ≥ (p − 1). We now choose θ as follows Observe that since s p ≥ (p − 1), then Θ(s, p) = 1 and we always have 5 We thus obtain from (6.8) |u − u r,θ | dx dt ≤ C r δ , for some C = C(δ, K δ , N, s, p) > 0. 5 Indeed, observe that thanks to the fact that 0 < δ < 1. This in turn implies as claimed.
By the characterization of Campanato spaces on R N +1 with respect to a general metric (see [9,Teorema 3.I] and also [14,Theorem 3.2]), this implies that u is δ−Hölder continuous in Q 1/4,1/4 with respect to the metric By keeping (6.9) into account, we can infer that d is a true metric. Thus, in particular, we have the estimate where C = C(K δ , N, s, p) > 0. Observe that the continuous function is increasing and that .
• Case s p < (p − 1). In this case, we revert the hierarchy between time and space and choose r as follows Observe that the exponent on θ is positive: indeed, for p = 2 this is straightforward, while for p > 2 we use that We further notice that now (6.10) s p − (p − 2) δ ≤ 1, up to choose δ sufficiently close 6 to s p/(p − 1). This time, we obtain from (6.8) |u − u r,θ | dx dt ≤ C θ δ s p−(p−2) δ , for some C = C(δ, K δ , N, s, p) > 0.
Again by the Campanato-type theorem of [9, Teorema 3.I], this shows that u is (δ/(s p − (p − 2) δ))−Hölder continuous in Q 1/4,1/4 with respect to the metric Observe that this is indeed a metric, thanks to (6.10). In particular, we have the estimate 6 More precisely, it is sufficient to take with 0 < ε < s/(p − 1) such that Such a choice is feasible, since now s p < (p − 1).
where C = C(δ, K δ , N, s, p) > 0. We now use that the continuous function is increasing and that lim δր s p Thus, for every γ < 1, there exists s < δ < s p/(p − 1) such that This concludes the proof in this case, as well.

Proof of the main theorem
Before proving our main result, we will need the following lemma, which allows us to control the parabolic Sobolev-Slobodeckiȋ seminorm of a local weak solution u in terms of its L ∞ norm.
Lemma 7.1. Let p ≥ 2 and 0 < s < 1. Let u be a local weak solution of for some C = C(N, s, p) > 0.
Proof. Without loss of generality, we may suppose that x 0 = 0. Let us set Then u is a local weak solution in B 2 × (−2 R s p , 0] and u ≥ 1 in R N × [−R s p , 0]. For all φ(x, t) = η(x) ψ(t) with ψ ∈ C ∞ such that ψ(t) = 0 for t ≤ −R s p and ψ(0) = 1, and η ∈ C ∞ 0 (B R ), we get from a slight modification of We choose η such that and ψ such that It is then a routine matter to show that where C = C(N, s, p) > 0 and we used that p ≥ 2 and k ≥ 1. This proves the lemma.
We are now in the position to prove Theorem 1.2.
Proof of Theorem 1.2. The continuity in space is contained in Theorem 5.2, thus we only need to prove the continuity in time. We take for simplicity T 0 = 0. If u is a local weak solution in B 2 × (−2 R s p , 0], we obtain from (5.2) An application of Lemma 7.1 gives We set Appendix A. Existence for an initial boundary value problem In order to give the definition of weak solution for an initial boundary value problem, we need to define a suitable functional space. We assume that Ω ⋐ Ω ′ ⊂ R N , where Ω ′ is a bounded open set in R N . Given a function ψ ∈ W s,p (Ω ′ ) ∩ L p−1 s p (R N ), we define as in [19] (see also [4,Proposition 2.12]) the space When ψ ≡ 0, the boundedness of Ω ′ entails that We endow the space X s,p 0 (Ω, Ω ′ ) with the norm W s,p (Ω ′ ), then this is a reflexive Banach space. Thanks to the previous inclusion, we also have that (W s,p (Ω ′ )) * ⊂ (X s,p 0 (Ω, Ω ′ )) * . Definition A.1. Let I = [t 0 , t 1 ] and p ≥ 2. With the notation above, assume that the functions u 0 , f and g satisfy u 0 ∈ L 2 (Ω), f ∈ L p ′ (I; (W s,p (Ω ′ )) * ), g ∈ L p (I; W s,p (Ω ′ )) ∩ L p−1 (I; L p−1 s p (R N )) and ∂ t g ∈ L p ′ (I; (W s,p (Ω ′ )) * ). We say that u is a weak solution of the initial boundary value problem if the following properties are verified: • u ∈ L p (I; W s,p (Ω ′ )) ∩ L p−1 (I; L p−1 s p (R N )) ∩ C(I; L 2 (Ω)); • u ∈ X g(t) (Ω, Ω ′ ) for almost every t ∈ I, where (g(t))(x) = g(x, t); • lim t→t0 u(·, t) − u 0 L 2 (Ω) = 0; The starting point for proving the existence of weak solutions is an abstract theorem for parabolic equations in Banach spaces. Before stating the theorem, we will briefly explain its framework. Let V be a separable reflexive Banach space and let H be a Hilbert space that we identify with its dual, i.e. H * = H. Suppose that V is dense and continuously embedded in H. If v ∈ V and h ∈ H, we identify h as an element of V * through the relation 7 Here ·, · denotes the duality pairing between V and V * and (·, ·) H denotes the scalar product in H. Let I be an interval and 1 < p < ∞. By [28, Proposition 1.2, Chapter], we have , v(t) . More generally, by [28, Corollary 1.1, Chapter III], for every u, v ∈ W p (I) the scalar product t → (u(t), v(t)) H is an absolutely continuous function and there holds d dt for a. e. t ∈ I.
We recall that an operator A : V → V * is said to be 7 With these identifications, we have V ⊂ H ⊂ V * . This is sometimes called in the literature Gelfand triple.
• monotone if for every u, v ∈ V , • hemicontinuous if the real function λ → A(u + λ v), v is continuous, for every u, v ∈ V . (ii) for almost every t ∈ I, the operator A(t, ·) : V → V * is monotone, hemicontinuous and bounded by for v ∈ V and k ∈ L p ′ (I), (iii) there exist a real number β > 0 and a function ℓ ∈ L 1 (I) such that Then for each f ∈ V * = L p ′ (I; V * ) and u 0 ∈ H, there exists a unique u ∈ W p (I) satisfying This means that u ∈ V, u ′ ∈ V * and Proof. The existence of a unique solution u ∈ V is contained in [28,  In order to prove existence for our problem (A.1), we will use Theorem A.2 with the choice V = X s,p 0 (Ω, Ω ′ ). This is the content of the next result, which generalizes [25,Theorem 2.5]. The latter only deals with the case f ≡ g ≡ 0.
Observe that this is well-defined, since v + g(t) ∈ X s,p g(t) (Ω, Ω ′ ), for every v ∈ X s,p 0 (Ω, Ω ′ ). We next show that the operator A, together with the spaces V = X s,p 0 (Ω, Ω ′ ), V = L p (I; X s,p 0 (Ω, Ω ′ )) and H = L 2 (Ω), fits into the framework of Theorem A.2. Since p ≥ 2 and Ω ′ is bounded, X s,p 0 (Ω, Ω ′ ) is dense and continuously embedded in L 2 (Ω). This follows from Hölder's inequality and the fact that smooth functions are dense in both spaces. Note that A inherits the property of monotonicity from A t since We next claim that We have The first term on the right-hand side of (A.5) can be bounded by using Hölder's inequality. For the second term we observe that, when x ∈ Ω and y ∈ R N \ Ω ′ , where C > 1 depends only on the distance between Ω and Ω ′ . Since 1/(1 + |y| N +s p ) ∈ L 1 (R N ), the second term in the right-hand side of (A.5) can be estimated by where we used the continuous inclusion L p (Ω) ⊂ W s,p (Ω ′ ). This finally shows (A.4). Observe that thanks to the assumptions on g. Thus in order to verify (ii) of Theorem A.2, we are left with proving hemicontinuity. For this, fixed t ∈ I and λ, λ 0 ∈ R, we consider for u, v ∈ X s,p 0 (Ω, Ω ′ ). In order to show that this differences goes to 0 as λ goes to λ 0 , it is sufficient to write and then use [19,Lemma 3]. This proves that A is hemicontinuous for almost every t ∈ I.

Thus we obtained
for every J = [T 0 , T 1 ] ⊂ I and every φ ∈ L p (J; X s,p 0 (Ω, Ω ′ )) ∩ C 1 (J; L 2 (Ω)). By recalling the definition of A t , this shows u is a weak solution of (A.1).
Given an initial datum u 0 ∈ L 2 (Ω), we consider the unique weak solution u to the initial boundary value problem   If there exists M ∈ R such that Proof. We take J = [T 0 , T 1 ] ⋐ (t 0 , t 1 ), by proceeding as in the first part of Lemma 3.3, we obtain for every φ ∈ L p (J; X s,p 0 (Ω, Ω ′ )) ∩ C 1 (J; L 2 (Ω)). We still use the notation φ ε and u ε for the convolution in the time variable, as defined in (3.3). Moreover, we still indicate by Σ(ε) the error term (3.6). We now take the test function 8 φ(x, t) = (u ε (x, t) − M ) + . Observe that this function is only Lipschitz in time, but it is not difficult to see that Lipschitz functions are still feasible test functions (by a simple density argument). This giveŝ By taking the limit as ε goes to 0, we thus get T1 T0¨R N ×R N J p (u(x, t) − u(y, t)) (u(x, t) − M ) + − (u(y, t) − M ) + dµ(x, y) dt By using that (see [5,Lemma A.2]) we thus getˆT This is valid for every t 0 < T 0 < T 1 < t 1 . By using that u ∈ L p (I; W s,p (Ω ′ )) ∩ L p−1 (I; L p−1 s p (R N )) ⊂ L p (I; W s,p (R N )), and that u ∈ C(I; L 2 (Ω)), we can pass to the limit as T 0 goes to t 0 and obtain We used that u 0 ≤ M , by construction. This implies that u(x, T 1 ) ≤ M, for a. e. x ∈ Ω.
Since T 1 is arbitrary, we finally get that u(x, t) ≤ M, for a. e. x ∈ Ω, for t ∈ I.
This concludes the proof.
As a straightforward consequence of the previous result, we get the following We also include the following comparison principle with bounded subsolutions.
Then u(x, t) ≥ v(x, t), in R N × I.

By [4, Lemma A.3], we have
for some C = C(p) > 0. Then for every t 0 < T 0 < T 1 < t 1 . We can now let T 0 converge to t 0 and obtain This implies u(x, T 1 ) ≥ v(x, T 1 ), for a. e. x ∈ Ω.
Since T 1 is arbitrary, this entails the desired result.
We start by splitting the integral as followŝ We now observe that where we used the properties of convolutions, the fact that F is locally Lipschitz and the uniform L ∞ bound (3.11). Thus, up to extracting a subsequence, we can infer weak convergence in L p ([T 0 , T 1 ]; W s,p (B 2−2 h )), of F (δ h u ε (x, t)) τ ε (t) ε η p , to the function F (δ h u(x, t)) τ (t) η p .
By definition, this is the same as saying that the function , belongs to L p ′ ([T 0 , T 1 ]; L p ′ (B 2−2 h × B 2−2 h )).
This in turn permits to infer that Θ 2 (ε) goes to 0, as well. This concludes the proof of Lemma 3.3.