{'fraction_non_alphanumeric': 0.06559611429137858, 'fraction_numerical': 0.021411047505761654, 'mean_word_length': 3.555141663844678, 'pattern_counts': {'":': 0, '<': 20, '<?xml version=': 0, '>': 29, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On regularity of weak solutions to linear parabolic systems with measurable coefficients

We establish a new regularity property for weak solutions of parabolic systems with coefficients depending measurably on time as well as on all spatial variables. Namely, weak solutions are locally H{\"o}lder continuous Lp valued functions for some p>2.


Introduction
This work is concerned with local regularity of weak solutions to linear parabolic equations or systems in divergence form, (1.1) in absence of any regularity of the coefficients besides measurability. The system is considered in an open parabolic cylinder (t, x) ∈ Ω ⊆ R × R n , n ≥ 1, ellipticity is imposed in the sense of a weak Gårding inequality and weak solutions belong to the usual Lions class, that is to say, u and ∇ x u are locally square-integrable. We note at this stage that we do not impose solutions to be locally bounded in any sense.
The topic has a long history, probably starting with the famous results of Nash [31] and Moser [30] that weak solutions to parabolic equations with real coefficients are Hölder continuous with respect to the parabolic distance. This is not true for equations with complex coefficients let alone systems, even in dimension n = 2. First counterexamples were constructed by Frehse-Meinel [14] in dimension n ≥ 3 and very recently by Mooney [29] in dimension 2. Results obtaining Hölder estimates beyond the Nash-Moser theorem have mostly focused on either quasilinear equations and systems with more regular coefficients, notably C 1 -regular in time and space [7] or absolutely continuous in time [15], and systems with a very specially structured real coefficient matrix such as the diagonal systems in [8]. As for most general linear parabolic systems, which are considered in this work, Lions [27] showed continuity in time valued in spatial L 2 loc . This was improved to L p loc for some p > 2 in the case of real coefficients by Nečas-Sverák [33], and Giaquinta-Struwe obtained the local higher integrability of ∇ x u in [17]. In this paper, we reveal a new regularity property of weak solutions. In its simplest form it can be stated as follows. Theorem 1.1. If u is a local weak solution to the homogeneous system ∂ t u − div x A(t, x)∇ x u = 0, then u is locally bounded and Hölder continuous in time with values in spatial L p loc for some p > 2. Most of the regularity properties of solutions to parabolic systems have been established through the local variational methods emerging from the Lions theory [27]. It does not seem that those methods give access to our result. Instead we rely on a global variational approach based on this simple observation: We can extend the local solution u via multiplication with a smooth cut-off to a global function v and study the corresponding inhomogeneous problem (1.2) now on all of R n+1 . Indeed, any local information of v carries over to u, but the global setup with the real line describing the time enables to bring powerful tools such as singular integral operators and the Fourier transform into play. Most notably, splitting ∂ t = D 1/2 according to the Fourier decomposition −iτ = |τ | 1/2 (−i sgn(τ ))|τ | 1/2 , there is a sesquilinear form (1.3) corresponding to (1.2) which, in contrast to the Lions theory on finite time intervals, admits a hidden coercivity on a natural energy space. This uses the algebraic properties of the Hilbert transform on the real line in a crucial way. See Lemma 3.2 below for details. These observations are not new but it does not seem they have been fully exploited for obtaining local regularity of solutions up to now. To the best of our knowledge, the idea first appeared in the work of Kaplan [24]. This idea was rediscovered by Hofmann-Lewis [21] in the context of parabolic boundary value problems (see also [5] and the references therein) and has recently played a role in the realm of non-autonomous maximal regularity [4,13]. Half-order derivatives on finite time intervals were also studied in [26,Ch. VII].
Compared to the local variational approach, where the t-derivative is understood in a weak sense only through the equation (see Section 2), we can now study the exact amount of differentiability that v should admit on the global level through a locally integrable function -the fractional derivative D 1/2 t v(t, x). In fact, we have a priori a parabolic differential |Dv| := |∇ x v| + |H t D 1/2 Since time is a one-dimensional variable, square-integrability of D 1/2 t v is already the borderline case from the view-point of Sobolev embeddings, not enough for continuity in time though, which probably explains why D 1/2 t v has not been exploited before. On the other hand, higher integrability of D 1/2 t v would yield Hölder continuity in time. Indeed most of this paper is devoted to establishing the selfimprovement of integrability for the spatial gradient and the half-order time derivative simultaneously, that is to say, |Dv| ∈ L p (R; L p (R n )), for some p > 2. (1.4) We present two proofs relying on rather different methods, both using the global variational formulation explained above. We think they each have their own interest, with potential applicability to non-linear systems for the first one and to other types of parabolic equations as well as fractional elliptic equations for the second one.
1.1. Strategy of the proofs. In Section 6 we present a real analysis proof of (1.4). The idea is to use, as in the analogous result for elliptic equations [28], self-improvement properties of reverse Hölder inequalities known as Gehring's lemma. First, we prove a new and delicate reverse Hölder inequality for Dv, by extending ideas from [5]. The non-locality of the fractional derivative reflects in local L 2 averages of Dv being controlled only by a weighted infinite sum of L 1 averages. Hence, we need a substantial extension of the classical Gehring lemma, which we shall prove in Section 5 and could be of independent interest. An unrelated Gehring type lemma "with tail" recently obtained in the context of fractional elliptic equations [25] has been inspiring to us. We mention that we shall explore further such extensions in a forthcoming work [2]. For other modifications of the local parabolic Gehring lemma we refer to [6] and references therein.
In Section 7 we present an operator theoretic proof of (1.4). We consider the operator L associated with the sesquilinear form (1.3) in virtue of the Lax-Milgram lemma and use an analytic perturbation argument. More precisely, exploiting the hidden coercivity in a crucial way, L plus a large constant turns out to be invertible from the natural L 2 energy space to its dual and extends boundedly to the corresponding L p -based spaces. The higher integrability of Dv then follows from the fact that invertibility of a bounded operator between complex interpolation scales extrapolates. The latter is known asŠneȋberg's theorem [35].

Main results.
All this is for homogeneous equations so far, that is f = 0 and F = 0. However, as the extension to R n+1 forces us to work with inhomogeneous equations anyway, there is no real obstacle to start with an inhomogeneous equation right away. Here we give an informal description of our main results. Precise theorems and their proofs are found in Sections 4 and 6 -8.
We consider weak solutions in Ω = I 0 × Q 0 to (1.1) under the assumptions f ∈ L q ′ loc (Ω) and F ∈ L 2 loc (Ω), where q = 2 + 4/n and q ′ is its Hölder conjugate. It is worth mentioning that q ′ < 2 and hence we go beyond the usual Lions variational approach that uses ∂ t u ∈ L 2 loc (I 0 ; W −1,2 loc (Q 0 )) as a starting point to obtain in-time continuity of u valued in spatial L 2 loc . We could have made u ∈ L ∞ loc (I 0 ; L 2 loc (Q 0 )) an assumption as for instance in [32] but in fact -and this is an observation we have not found in the literature -this uniform local L 2 boundedness in space and the local L q integrability both follow from the hypotheses, again thanks to the global variational approach and the use of half-order derivatives. More precisely, we are still able to obtain the "classical" cornerstones of the Lions theory: Next, if p > 2 is sufficiently close to 2, depending only on ellipticity and dimensions, then we have under the assumptions f ∈ L p * loc (Ω) and F ∈ L p loc (Ω) the following improvements. Here, p * = (n + 2)p/(n + 2 + p).
As explained before, (iv) and (v) are completely new in this generality and the main contribution of this work. Moreover, (vi) was first obtained in [17] when f = 0, F = 0 by means of the classical Gehring lemma and was generalized to non zero right-hand side in [12], but with stronger requirements on f and F = 0. Such results have impact on partial regularity of nonlinear systems [12,16,17].
In the following Section 2 we introduce relevant notation. The global variational setup using (1.3) is discussed afterwards in Section 3.

Notation and basic definitions
Most of our notation is standard. One exception is that for X a Banach space we let X * be the (anti-)dual space of conjugate linear functionals on X. For exponents p ∈ (1, ∞) we define the upper and lower Sobolev conjugate p * and p * with respect to parabolic scaling and the Hölder conjugate p ′ through the relations whenever they belong again to the interval (1, ∞). Hence, for the exponent q used above we have q = 2 * and q ′ = 2 * . With regard to parabolic systems and their weak solutions, we use the following notions.

Ellipticity.
In what follows we assume bounded, measurable, complex valued coefficients x)) α,β=1,...,m i,j=1,...,n ∈ L ∞ (R n+1 ; L(C mn )) and that there exist λ > 0 and κ ≥ 0 such that the (weak) Gårding inequality holds for all u ∈ W 1,2 (R n ; C m ), uniformly in t ∈ R. Our notation is where for the sake of readability we shall usually stick to the summation convention on repeated indices and do not write out sums explicitly. We shall refer to λ, κ and an upper bound for the L ∞ -norm of A as ellipticity and to n and the number m ≥ 1 of equations as dimensions.
Let us remark that for the local results we are after, it is no restriction to define A on all of R n+1 . Indeed, if, for some open interval I 0 ⊂ R and ball Q 0 ⊂ R n , coefficients A ∈ L ∞ (I 0 × Q 0 ; L(C mn )) satisfy (2.2) only for u ∈ W 1,2 0 (Q 0 ; C m ) uniformly in t ∈ I 0 , then given ε ∈ (0, 1), there exists A with We denote by ℓ(I 0 ) the length of I 0 and by r(Q 0 ) the radius of Q 0 . Given f ∈ L 1 loc (Ω; C m ) and F ∈ L 1 loc (Ω; C mn ), we say that u is a weak solution to in Ω if u ∈ L 2 loc (I 0 , W 1,2 loc (Q 0 ; C m )) and for all smooth functions with compact support φ ∈ C ∞ 0 (Ω; C m ), Having posed the setup, whenever the context is clear we are going to ignore the target spaces C m or C mn in our notation and do not write the Lebesgue measures dx and dt. We abbreviate ∇ := ∇ x and div := div x for the gradient and divergence in the spatial variables x, respectively.

Fractional time derivatives and related spaces.
In the following D 1/2 t and H t denote the half-order time derivative and Hilbert transform in time defined on S ′ (R)/C, the tempered distributions modulo constants, through the Fourier symbols |τ | 1/2 and −i sgn(τ ), respectively. For summarizing properties see for example Section 3 in [5]. In particular, the time derivative factorizes as ∂ t = D 1/2 We shall use the space H 1/2 (R; L 2 (R n )) of functions in L 2 (R n+1 ) such that D Here, we identify L 2 (R n+1 ) with L 2 (R; L 2 (R n )), and having said this, we extend D 1/2 t and H t to R n+1 by acting only on the time variable.
More generally, for 1 < p < ∞ we introduce the spaces For the sake of completeness only, we remark that up to equivalent norms these are the (vector-valued) Bessel potential spaces usually denoted by the same symbol [9]. We also note that C ∞ 0 (R n+1 ), the space of smooth and compactly supported functions, is dense in these spaces using smooth convolution and truncation. Lebesgue space norms are denoted with the usual symbol · p .

The global variational setup
Our starting point is the parabolic problem on the whole space R n+1 . We define the Hilbert space V := L 2 (R; W 1,2 (R n )) with norm u V := ( u 2 2 + ∇u 2 2 ) 1/2 . This is the natural space for global weak solutions to parabolic problems. We recall that 2 * = 2 + 4 n is the upper Sobolev conjugate of 2 and that (2 * ) ′ = 2 * is its dual exponent.
The following proposition summarizes the properties of global weak solutions.
The implicit constant depends only on dimensions and ellipticity.
We need a few short lemmas to prepare its proof. The key tool is the following definition of the parabolic operator through a sesquilinear form and its accretivity on the parabolic energy space The result is basically that of [21,24] but we repeat the short argument for the reader's convenience.

This operator is invertible and its norm as well as the norm of the inverse depend only on ellipticity and dimensions.
Proof. The E → E * boundedness of L is clear by definition. Next, for the invertibility, the form for u, v ∈ E, is bounded and satisfies an accretivity bound for δ > 0 sufficiently small. Indeed, from the ellipticity condition and the fact that the Hilbert transform is L 2 -isometric, skew-adjoint and commutes with D 1/2 t and ∇, The following lemma is well-known on the Hilbert space V , see [34, Prop. III.1.2], but we need a slight variant involving the smaller spaces Then v ∈ C(R; L 2 (R n )) and t → v(t, · ) 2 2 is absolutely continuous on R, vanishes at ±∞ and satisfies Proof. We can identify V p and V * p canonically with subspaces of the (tempered) distributions on R 1+n . Indeed, C ∞ 0 (R n+1 ) is dense V p through approximation via smooth convolution and truncation. We approximate v through convolution with smooth compactly supported kernels in the t-variable only, say v ε = v * t ϕ ε , where ε > 0. By construction v ε is in the class C ∞ (R; L 2 (R n ) ∩ L p (R n )) and vanishes at ±∞. Hence we can differentiate in the classical sense Integrating this identity over an interval I, we obtain . As for ∂ t v ε , we have boundedness in V * p uniformly in ε and weak * -convergence towards ∂ t v. Hence, we can pass to the limit as ε → 0 in the equation above.
In fact, we can derive such an estimate also for differences . This means that t → v(t, ·) 2 2 belongs to W 1,1 (I) and therefore can be re-defined as an absolutely continuous function.
If we take I = (−∞, T ], with T such that at least for a subsequence v ε (T, · ) 2 2 → v(T, · ) 2 2 , then the fundamental theorem of calculus reveals that the left-hand side of (3.1) tends to v(T, · ) 2 2 whereas the right-hand side tends to 2 Re We may obtain the final statement Re ∂ t v, v = 0 by taking I = R and passing to the limit as ε → 0.
Finally, we need a slight variant of the usual parabolic Sobolev embeddings. For background, we refer to [20].
We can now give the Proof of Proposition 3.1. In virtue of the canonical identifications we have the continuous inclusion E ⊂ E * and a bounded mapping div : L 2 (R n+1 ) → E * . It follows from Lemma 3.4 and density of C ∞ 0 (R n+1 ) in E (standard mollification and truncation), that E embeds into L 2 * (R n+1 ). Hence, It follows from Lemma 3.2 that there exists v ∈ E such that L v = f + (κ + 1)v + div F . By definition of the respective embeddings, this means that for all φ ∈ E, We may now apply the Caccioppoli inequality, for any parabolic cylinder I × Q with ℓ(I) = r(Q) 2 , see Remark 4.7 below for convenience. Since we have u ∈ L 2 (R n+1 ), we obtain ∇u = 0 on passing to the limit r(Q) → ∞. Hence, u depends only on t. Again, as u ∈ L 2 (R n+1 ), u must be 0. It follows that v = v ∈ E, hence (i) is proved and follows by Lemma 3.2. Applying Lemma 3.4 again yields (ii) by density. As for (iii) we have seen v ∈ V 2 * , which in turn implies ∂ t v ∈ V * 2 * by the equation for v. Hence, we can apply Lemma 3.3 to v and obtain the statements on continuity. We also obtain sup t∈R v(t, · ) 2 2 ≤ 2 v Vp ∂ t v V * p and so the required bound follows on controlling the right-hand side by means of (i) and (ii).

Local estimates
As a first application of the global results obtained in the previous section we present the "classical" local estimates for weak solution within our setting. We recall that by definition a weak solution u to a parabolic problem in a parabolic cylinder Ω = I 0 × Q 0 satisfies u ∈ L 2 loc (Ω) and ∇u ∈ L 2 loc (Ω). The following lemma is nothing but a simple calculation. Nevertheless, it is of fundamental importance for all subsequent considerations. Here, we suggestively use the notation for the scalar case m = 1 (even when m > 1), as we are only interested in norm estimates later on.
With this at hand we can prove the local higher integrability and absolute continuity of weak solutions.

Theorem 4.2. Assume that u is a weak solution to
Using the assumption on f, F , the local square-integrability of u and ∇u, and 2 * < 2, we see that f ∈ L 2 * (R n+1 ) and F ∈ L 2 (R n+1 ). The conclusion follows from Proposition 3.1.
We continue with the Caccioppoli inequality. It will be convenient to formulate it with an additional zero-order term on the right-hand side.

Proposition 4.3 (Caccioppoli inequality). Let u be a weak solution in
The implicit constant depends on ellipticity, dimensions, γ, constants controlling the ratio r 2 /ℓ and B ∞ .
Before we give the proof of Caccioppoli's inequality, let us conclude the reverse Hölder estimate of Nečas-Šverák [33] for |u| in our setting.
where the implicit constants depend only on ellipticity, dimensions, γ and the constants controlling the ratio r 2 /ℓ.
Proof. The equation (4.4) with − − γ 2 I×γQ |u| 2 1/2 on the right-hand side follows from (4.2) and Proposition 4.3 -at least when r = 1, which suffices since we can rescale as before. The improvement to − − γ 2 I×γQ |u| follows from a classical self-improvement feature of reverse Hölder inequalities, see Theorem 2 in [22]. A simple proof that applies in our situation with parabolic scaling can be found in Theorem B.1 of [10].

Remark 4.5.
Under suitable assumptions on f and F the classical Gehring lemma (with parabolic scaling) can be used to improve the exponent of integrability on the left-hand side to 2 * + ε, where ε > 0 depends on ellipticity, dimensions, γ and the constants controlling the ratio r 2 /ℓ. For f = 0 and F = 0, the argument is found in the textbook [11].

Concatenating (4.3) and (4.4) yields
We complete the section with the proof of Caccioppoli's inequality. The argument follows the traditional one and can be omitted on a first reading.
Proof of Proposition 4.3. For the argument set q := 2 * with Hölder conjugate q ′ = 2 * . By scaling we may assume that r = 1 as our hypotheses are invariant under dilations. We pick χ ∈ C ∞ 0 (R n+1 ), real-valued, with χ = 1 on I × Q and support contained in γ 2 I × γQ. As in Lemma 4.1 we write the equation satisfied by v := uχ as We have v ∈ V q thanks to Proposition 3.1 and from the assumptions on f , F , B and Hölder's inequality we can infer ∂ t v ∈ V * q . Thus, Lemma 3.3 yields First, we note that the final integral on the right has positive real part by assumption. We then isolate Re(A∇v · ∇v) and use Young's inequality to give Now, we apply Gårding's inequality (2.2) on the left and cancel λ|∇v| 2 on both sides to obtain We pick yet another real-valued function φ ∈ C ∞ 0 (R n+1 ), supported in γ 2 I × γQ with φ = 1 on the support of χ. On recalling v = χu, we can hence write We insert this expression along with the definitions of v and F into (4.8). Then we can estimate all terms appearing on the right but f φuχ 2 simply by Cauchy-Schwarz' and Young's inequality and use the uniform bound for A whenever convenient. This results in λ |∇(uχ)| 2 I + |uχA∇u · ∇χ| + |f φ||uχ 2 |, (4.9) where I := |u| 2 (|∇χ| 2 + |∂ t χ| + |χ| 2 ) + |F | 2 (|χ| 2 + |∇χ| 2 ).
As for the second term on the right-hand side of (4.9), we write This allows us to apply Young's inequality with parameters chosen such that the contribution of ∇(uχ) appearing on the right-hand side of (4.9) can be absorbed into the left-hand side. Invoking again the uniform bound for A, all other terms created in this step will only increase the implicit constant in front of I. Thus, we can note |∇(uχ)| 2 I + |f φ||uχ 2 |. (4.10) Had we assumed f ∈ L 2 loc (Ω), then a simple application of the Cauchy-Schwarz inequality would complete the argument. But we only assumed f ∈ L q ′ loc (Ω). In order to master the situation, we introduce w = uχ 2 , which is of the same nature as v except that the cut-off function changed from χ to χ 2 . We also define f ′ , F ′ just as f , F but with χ replaced by χ 2 . Hölder's inequality followed by the global bound provided by Proposition 3.1 (iii) for w yields Crudely using Hölder's inequality on the definition of f ′ , we find We see that up to changing the cut-off function from χ to χ 2 , the terms F ′ 2 and w 2 on the right of (4.11) have already been estimated before when passing from (4.8) to (4.9). Repeating these arguments,

Remark 4.7.
Under the stronger assumption f ∈ L 2 loc (Ω) the proof given above yields without making use of Proposition 3.1. This observation is important in order to make clear that there is no circular reasoning going on in the proof of the latter. Indeed, we can replace q by 2 so that u ∈ V q = V 2 is by definition of weak solutions and then we follow the proof verbatim until we reach (4.10), where now we only have to apply the Cauchy-Schwarz inequality to conclude.

A Gehring type lemma with tail
We provide here the main real analysis lemma to obtain our estimates. For a ball Q ⊂ R n and an interval I ⊂ R with ℓ(I) = r(Q) := r we write B := I × Q. If (t, x) is the center of B, we also use the notation B = B ((t, x), r) and r = r(B) for such a parabolic cylinder (that is, a cylinder which is a ball in the parabolic (quasi-)metric d((t, x), (s, y)) := max( 2|t − s|, |x − y|)). For u ≥ 0 locally integrable and γ > 1 we define where for this section µ denotes the Lebesgue measure on R n+1 and we use the single integral notation for simplicity. The functional a u is an approximate identity indexed over radii of parabolic cylinders when u ∈ L p (R n+1 ) for some p ∈ (1, ∞) in the sense that a u (B((t, x), r)) → u(t, x) as r → 0 for almost every (t, x). Indeed, introduce the maximal operators M x and M t on space and time variables separately. For each j ≥ 0 we have and as M t M x is bounded on L p (R n+1 ), this average converges to u(t, x) for almost every point. So, the claim follows from the dominated convergence theorem for series and ∞ j=0 2 −j−1 = 1. In addition, we have a u (B((t, x), r)) → 0 when r → ∞ by Hölder's inequality. This last point also holds when u ∈ L 1 (R n+1 ).

Lemma 5.1.
Let g, f, h be be non-negative functions with g 2 , f 2 , h s ∈ L 1 (R n+1 ) for some 1 < s < n + 2, and suppose that for some A ≥ 1, a h s (B)) 1/s holds for all parabolic cylinders B. Let p > 2 and suppose there are α, β ≥ 0 and q α , q β > 1 (depending on p) such that If |p − 2| is sufficiently small (depending on A and dimension), then The implicit constant depends on A, α, β, q α , q β and dimension. Proof. Let m > 0. Denote g m := min(g, m). By the Cavalieri principle we have where g 2 (A) := A g 2 dµ. We define three functions G(t, x, r) := a g (B((t, x), r)), (B((t, x), r)) 1/2 , (B((t, x), r)) 1/s and for λ > 0, we denote E λ := {g > λ}. We have for almost every (t, x) by the discussion before the statement of the lemma and we define E λ as the subset of E λ where this holds. We also note lim r→∞ G(t, x, r) + F (t, x, r) + H(t, x, r) = 0 for all (t, x), using the global assumptions on g, f, h and s < n + 2.
By definition, if (t, x) ∈ E λ , then lim r→0 G(t, x, r) > λ, and thus for (t, x) ∈ E λ we can define the stopping time radius We readily see that sup (t,x)∈ E λ r t,x < ∞. Indeed, since G, F, H are continuous functions of r > 0 for fixed (t, x), we have at r = r t,x equality G(t, x, r) + F (t, x, r) + H(t, x, r) = λ and thus either G(t, x, r) ≥ λ/3 or F (t, x, r) ≥ λ/3 or H(t, x, r) ≥ λ/3. In the last case for example, we obtain Now, using the hypothesis for each B i and pairwise disjointness of the balls 1 K B i , we find Let M β x be the fractional maximal function with respect to the x-variable: Similarly, define M α t with respect to the t-variable. Since 2α + β = s, the parabolic scaling r(B) = r(Q) = ℓ(I) yields r(B) s = r(Q) β × ℓ(I) α . Thus, it follows from the definition of r i that We thus have established Going back to the start of the proof, so far we have found where the integrals I, II, III correspond to the decomposition of S λ above. By the Cavalieri principle, we obtain for p > 2, by iterating the two maximal function L p/2 bounds, so that the implicit constant depends only on the dimension n. Note that p > 2 and that p is determined by the other parameters in (5.1). Similarly, By hypothesis, we have exponents q α , q β > 1 such that With a slight abuse in our notation, ignoring the other variable, these are precisely the conditions guaranteeing that M α t : L qα (R) → L p/s (R) and M β x : L q β (R n ) → L p/s (R n ) are bounded, see Theorem 3.1.4 in [1]. Now, using this and Minkowski's inequality along with sq α /p = 1 − αq α ≤ 1 in the second step, we see that with implicit constant depending on α, β, q α , q β and dimension. The remaining term is To handle I, we first notice that From this inclusion and the weak type ( 3 2 , 3 2 )-bound for the iterated maximal function (which follows from the strong type ( 3 2 , 3 2 )), we obtain for a dimensional constant C. Using this bound in I yields Choosing p − 2 > 0 small enough, depending on A and dimension, we see from (5.2) that R n+1 for some constant C depending on α, β, q α , q β and dimension. This finishes the proof after simplifying the first term and taking the limit m → ∞.
Remark 5.2. The same estimate also holds with the mixed norm in different order as we are free to interchange the fractional maximal functions. If we want q α = q β , then (5.1) reveals that α, β are uniquely determined by s, p and n. Hence, for each p there is at most one such pair.
In the application to our parabolic PDE, we consider special values for the auxiliary parameters in Lemma 5.1.
6. Higher integrability of the parabolic differential: Real analysis proof Let Ω = I 0 × Q 0 be the ambient parabolic cylinder and u a weak solution to (1.1) in Ω. Given χ ∈ C ∞ 0 (Ω) we know from Section 3 that under suitable assumptions on f and F the (localized) parabolic differential belongs to L 2 (R n+1 ). In this section, we give a first proof of the following higher integrability result, which lies at the heart of our considerations. Since the Hilbert transform is isometric on L 2 (R) it does not matter whether or not we include H t D 1/2 t (uχ) here but for L p -results it seems appropriate to treat both half-order derivatives. Theorem 6.1. Suppose p > 2 is sufficiently close to 2, depending only on ellipticity and dimensions. If u is a weak solution in Ω to ∂ t u − div A∇u = f + div F , where f ∈ L p * loc (Ω; C m ) and F ∈ L p loc (Ω; C mn ), then . In this section we want to give a proof using the Gehring lemma with tail. The required non-local reverse Hölder estimate is provided by the following key lemma. Naturally its proof is somewhat technical and will be postponed until the end of this section. It can be skipped on a first reading.
Let γ > 1 and let I × Q be a parabolic cylinder with ℓ(I) ∼ r(Q) 2 . Then Here, I k := kℓ(I) + I are the disjoint translates of I covering the real line up to a countable set. The implicit constant depends only on ellipticity, dimensions, γ and the constants controlling the ratio r(Q) 2 /ℓ(I).
For the moment, let us admit the lemma and record its consequences.
Let g be as in Lemma 6.2. If p > 2 is sufficiently close to 2, then The implicit constant as well as p depends only on ellipticity and dimensions.
Proof. Rearranging unions of translates of an interval I into unions of its dilates, and vice versa, reveals that for any positive function h on the real line we have with absolute implicit constants. Lemma 6.2 together with this observation and Hölder's inequality yields Thus, we have the setup of Lemma 5.1 with s = 2 * and we conclude by Corollary 5.3.
Theorem 6.1 is obtained through a by now well-known localization procedure.
Proof of Theorem 6.1. The function v := uχ ∈ V is a weak solution to ∂ t v−div A(t, x)∇v = f +div F on R n+1 with the relations (4.1). We can a priori assume 2 < p ≤ 2 * and hence have v ∈ L p (R n+1 ) from Theorem 3.1. Then we have L 2 loc (R n+1 ) ⊆ L p * loc (R n+1 ) since p * ≤ 2. This being said, f ∈ L p * (R n+1 ), F ∈ L p (R n+1 ) follows from (4.1) and the hypotheses on f , F . Hence, Corollary 6.3 applies and the claim follows.
We turn to the proof of Lemma 6.2. We follow the argument presented in Section 8 of [5] for f = 0 and F = 0. We omit duplicated arguments but give all other details so that the reader does not have to work through any other section of [5]. In this reference, the order of variables was (x, t) and an additional spatial dimension was carried through the argument, both for the purpose of treating boundary value problems. The latter plays no role here and can be ignored. Next, u in [5] has become v here and the extra property D 1/2 t v ∈ L 2 (R n+1 ) provided by Proposition 3.1 means that v is a reinforced weak solution in the terminology there.
Proof of Lemma 6.2. We remark that g ∈ L 2 (R n+1 ) due to Proposition 3.1 and the fact that H t is isometric on L 2 (R n+1 ). It suffices to prove the claim for γ = 8 since a posteriori a covering argument, which we leave to the reader, gives us the inequality with any γ > 1.
For simplicity, we are also going to assume r(Q) ∼ 1 and that I × Q is centered at (0, 0) as scaling and translating give us back the general estimate. Having normalized to scale 1, averages are integrals (up to numerical constants). For the time being it will be enough to work with γ = 4, so that the parabolic enlargement is 16I × 4Q. We fix a smooth cut-off η : R n+1 → [0, 1] with support in 4I × 2Q that is 1 on an enlargement 9 4 I × 3 2 Q. For a reason which will become clear later on, we choose η to have the product form where η I is symmetric about 0 (the midpoint of I). We also give a name to the translation sums |h|.
Step 1: The spatial average. The estimate (4.5) with v − c and c : and apply the L 1 -Poincaré's inequality to the first term and treat the second term (which does not depend on x) via the fractional Poincaré inequality in the following lemma with p = q = 1. Details are written out in the proof of Lemma 8.4 in [5].

The analogous inequality with
t h on the right-hand side also holds. The resulting estimate is Thus, we obtain a bound of the required type 4I×2Q |∇v| 2 Step 2: Local terms. For the local term w := η(v − c) we have using that H t is isometric on L 2 (R n+1 ). Since w solves an equation ∂ t w − div A∇w = f ′ + div F ′ , Proposition 3.1 implies as we can read off from (4.1). Using the formulae for f ′ , F ′ and Hölder's inequality to bound L 2 * averages by L 2 averages whenever convenient, we arrive at For the first term on the right we use (4.4) and then (6.2). For the second one we use (6.3). This leads to and decomposing 16I into translates of I as before gives an estimate of the required type.
Step 3: First error term. We come to the delicate steps in [5]. The non-locality of the operators D 1/2 t and H t D 1/2 t cannot be circumvented anymore. As η Q = 1 on Q, we have For the terms involving w 1 : we can use a kernel representation for D 1/2 t and then the fractional Poincaré inequality of Lemma 6.4. This is (the proof of) Lemma 8.6 in [5]. As a result, Inserting (6.3) for each 4I j × 2Q instead of 4I × 2Q and using the convolution inequality j∈Z 1 1 + |j| 3/2 1 1 + |k − j| 3/2 we obtain the desired bound by Step 4: Second error term. The remaining average of |H t D 1/2 is treated independently of knowing that v is a solution. Indeed, since v ∈ H 1/2 (R; L 2 (R n )) we have h, (1 − η I )h ∈ H 1/2 (R) and we had chosen η I in such a way that the following lemma applies. We note that in its proof the symmetry of η I is used to control the Hilbert transform, which has an odd kernel.
Lemma 6.5 ([5,Lem. 8.7]). Let I be a bounded interval and η I : R → [0, 1] be a smooth cut-off function with support in 4I that is identically 1 on 9 4 I. Suppose furthermore that η I is symmetric about the midpoint of I. If h, (1 − η I )h ∈ H 1/2 (R), then almost everywhere on I, We take the L 2 (I × Q) average on both side of (6.4) and use that H t D 1/2 t commutes with averages in the spatial variable to give This completes the required bound for the L 2 (I × Q) average of |H t D 1/2 t v| 2 with γ = 4 on the right-hand side.
It only remains to consider the bound for the L 2 (I × Q) average of D 1/2 . This can be done by another lemma on real functions from [5], with a weaker conclusion since D 1/2 t has an even kernel.
Indeed, for h = − 2Q v as before, we obtain The upshot is that we have already completed the reverse Hölder estimate for |H t D 1/2 t v| 2 on any parabolic cylinder, and in particular on 4I × 2Q, with spatial enlargement by a factor 4 on the righthand side. Hence, if we finally use γ = 8 on the right-hand side of (6.1), then the above estimate completes the reverse Hölder bound for D 7. Higher integrability of the parabolic differential: Operator theoretic proof We provide a second proof for the higher integrability of the parabolic differential using a completely different method. For 1 < p < ∞, we set E p := L p (R; W 1,p (R n )) ∩ H 1/2,p (R; L p (R n )) with norm u Ep := ( u p p + ∇u p p + D 1/2 t u p p ) 1/p , so that in particular E = E 2 is as in Section 3. These are Banach spaces with C ∞ 0 (R n+1 ) as a common dense subspace as is seen again by approximation via smooth convolution and truncation.
We shall use some results on complex interpolation of Banach spaces. The reader will find necessary background information in the textbook [9]. For the understanding of this paper it will be enough to know the complex interpolation identity in the sense of Banach spaces with equivalent norms. We shall say that (E p ) 1<p<∞ is a complex interpolation scale and the same holds true for the dual scale (E * p ′ ) 1<p<∞ . These assertions were proved in [5, Lem. 6.1].
We then have the following extension of Lemma 3.2.
Lemma 7.1. The operator L = ∂ t − div A(t, x)∇ + κ + 1 : E 2 → E * 2 extends by density to a bounded operator from E p to (E p ′ ) * for 1 < p < ∞. This extension is invertible for |p − 2| small enough and its inverse agrees with the one calculated when p = 2 on (E 2 ) * ∩ (E p ′ ) * . The norm of the inverse and the smallness of |p − 2| depend only on ellipticity and dimensions.
Proof. By definition, L : E 2 → E * 2 acts via Thus, E p → (E p ′ ) * boundedness of L follows from Hölder's inequality and the norm depends only on ellipticity and dimension. As L is invertible when p = 2 by Lemma 3.2, the invertibility for |p − 2| small enough follows fromŠneȋberg's result on bounded operators acting on interpolation scales, see [35] or [3, Thm. A.1] for a qualitative version revealing that the smallness of |p − 2| and the bound for the inverse depend only on ellipticity and dimensions. Finally, the compatibility of the inverses is an abstract feature of complex interpolation, see Theorem 8.1 in [23].
A simple but important consequence is sufficiently small (depending only on ellipticity and dimensions) and in this case with an implicit constant depending only on ellipticity and dimensions.
With this at hand, we are ready to give the second Proof of Theorem 6.1. Let χ ∈ C ∞ 0 (Ω). As before, v := uχ ∈ V is a weak solution to ∂ t v − div A(t, x)∇v = f + div F on R n+1 with f , F given by (4.1). By Proposition 3.1 we know that v ∈ E 2 and Lv = (κ + 1)v + f + div F . Let now p > 2 be such that we have Lemma 7.2 at our disposal. We may also suppose p ≤ 2 * , which is equivalent to p * ≤ 2. If we assume f ∈ L p * loc (Ω) and F ∈ L p loc (Ω), then f + (κ + 1)v ∈ L p * (R n+1 ) and F ∈ L p (R n+1 ), using also Theorem 6.1(i) to control uχ 2 * in the formula for F . Hence, by compatibility of the inverses (Lemma 7.1) and Lemma 7.2, we obtain The left-hand side controls uχ p + ∇(uχ) p + D 1/2 t (uχ) p and we are done.

Local higher regularity estimates
Eventually, we shall use the previously obtained qualitative information of higher integrability for the parabolic differential of the localized solution to obtain scale-invariant local higher regularity estimates. This is summarized in the following theorem. As usual, Ω = I 0 × Q 0 denotes the ambient parabolic cylinder.
Proof. We assume again r(Q) = 1 as rescaling will give us the right powers of r(Q). We follow the usual strategy and pick χ ∈ C ∞ 0 (Ω), χ = 1 on I ×Q, with support in γ 2 I ×γQ. Then x)∇v = f + div F on R n+1 with the relations (4.1) and f ∈ L 2 * (R n+1 ), F ∈ L 2 (R n+1 ). By Corollary 6.3 we have if p − 2 > 0 is small enough, Alternatively, we could have used (7.1) here at the expense of a term v p * v p on the right, which turns out to be harmless. Indeed, we have from (4.4) if p ≤ 2 * , as we may assume, We have used 2 * ≤ p * and 2 < p in the second step.
We have shown that v, D 1/2 t v are controlled in L p (R n+1 ) and in L 2 (R n+1 ). Since p > 2, a Hölder norm estimate on v will follow from classical embeddings. Indeed, for almost every x ∈ R n we can apply the fractional Poincaré inequality from Lemma 6.4 to v( · , x). Hence, for any interval J ⊂ R, and α = 1/2 − 1/p, we have where the implicit constant depends also on ℓ(J), see Theorem 2.9 in [18]. Taking J = I and averaging the p-th power of this estimate over x ∈ Q reveals that the left hand side of (8.1) is bounded by v p + ∇v p + D 1/2 t v p . (The reader should recall the normalization r = 1 ∼ ℓ and the construction of χ.) In view of (8.2) and (8.3), we see that it remains to control F p + f p * from above by the right-hand side of (8.1).
We begin with F . Let 1 < δ < γ be such that the support of χ is contained in δ 2 I × δQ. By (4.1) we have, Similarly, we use the definition of f in (4.1) to infer f p * δ 2 I×δQ |f | p * + |u| p * + |∇u| p * + |F | p * 1/p * and we are done as p * ≤ 2 and since the term integral over |∇u| p * can be treated using Proposition 4.4.
Finally, we obtain a true reverse Hölder estimate for ∇u, that is to say, an analogous estimate without u on the right-hand side. The implicit constant depends on ellipticity, dimensions, γ and the constants controlling the ratio r(Q) 2 /ℓ(I).
Proof. As in the proof of Proposition 4.4 the self-improving properties [10,22] of reverse Hölder estimates yield the conclusion once we have managed to prove (8.4) with an L 2 average of |∇u| on the right-hand side. By scaling it is also enough to assume r(Q) = 1. We use the "weighted means trick" introduced by Struwe in [36]. We choose χ real-valued, equal to 1 on I × Q and supported in γ 2 I × γQ of the form χ(t, x) = η(t)ϕ(x). Then define the weighted mean We set w(t, x) := (u(t, x) − u(t))η(t). We remark that ∇u = ∇w. It is thus enough to estimate I×Q |∇w| p 1/p . We proceed as follows. It follows from the equation for u that u is absolutely continuous on I with ∂ t u(t) = a −(A∇u + F ) · ∇ϕ + f ϕ dx (8.5) almost everywhere. Since η and u depend only on t, we have ∂ t w − div A(t, x)∇w = f ′ + div(ηF ) in Ω, where, omitting the variables except for the integration, Of course, we may assume p ≤ 2 * , which is equivalent to p * ≤ 2 and hence allows us bound L p * averages of |∇w| = |∇u| and |w| by the corresponding L 2 averages. As mentioned previously, it suffices to prove (8.4) with an L 2 average of |∇u| on the right hand side. So, we are left with controlling the L 2 average of w. Since taking weighted averages u → u defines a projection from W 1,2 (γQ) onto C, a variant of Poincaré's inequality on L 2 (γQ) discussed for example in [1,Lem. 8 The proof is complete. Remark 8.3. Once Theorem 6.1 is established, it is also possible to prove directly (8.4) under our assumptions by adapting the original argument in [17] and invoking the usual Gehring lemma.