Global higher integrability of weak solutions of porous medium systems

We establish higher integrability up to the boundary for the gradient of solutions to porous medium type systems, whose model case is given by \begin{equation*} \partial_t u-\Delta(|u|^{m-1}u)=\mathrm{div}\,F\,, \end{equation*} where $m>1$. More precisely, we prove that under suitable assumptions the spatial gradient $D(|u|^{m-1}u)$ of any weak solution is integrable to a larger power than the natural power $2$. Our analysis includes both the case of the lateral boundary and the initial boundary.


INTRODUCTION
We are concerned with the boundary regularity of solutions to Cauchy-Dirichlet problems of the form ∂ t u − div A(x, t, u, D(|u| m−1 u)) = div F in Ω T := Ω × (0, T ), for vector-valued solutions u : Ω T → R N , where m > 1, the domain Ω ⊂ R n is bounded with dimension n ≥ 2, the dimension of the target space is N ∈ N, and ∂ par Ω T denotes the parabolic boundary of the space-time cylinder Ω T , where T > 0. We cover a large class of vector fields A that we only require to satisfy growth and ellipticity conditions corresponding to the model case A(x, t, u, ζ) = ζ of the porous medium system. The assumptions on the data are made precise in Section 1.1 below. Our starting point are weak solutions, by which we mean in particular that the spatial gradient satisfies D(|u| m−1 u) ∈ L 2 (Ω T ). Our goal is to establish the self-improving property of integrability up to the boundary in the sense that D(|u| m−1 u) ∈ L 2+ε (Ω T ) holds true for some ε > 0.
The question for higher integrability of solutions has a long history that starts with the classical work by Elcrat & Meyers [31] on elliptic systems of p-Laplace type, which in turn is based on the work of Gehring [14]. Since then, similar results have been established for a variety of other elliptic problems, and the higher integrability of solutions has proved to be a very useful tool for the derivation of further regularity results. We refer to [18,19,17,21] and the references therein. The question of higher integrability up to the boundary for equations of p-Laplace type has been answered positively by Kilpeläinen & Koskela [25]. They observed that the natural condition to impose on the regularity of the domain Ω ⊂ R n is the property of uniform p-thickness of the complement R n \ Ω, see [25,Rem. 3.3].
The first higher integrability result for a parabolic problem is due to Giaquinta & Struwe [20], who treated the quasilinear case. However, it turned out that the techniques of Elcrat & Meyers could not directly be extended to the case of the parabolic p-Laplace system due to the anisotropic scaling behaviour of this system. This problem was solved by Kinnunen & Lewis in [26] for weak solutions to p-Laplace type systems. The much more intricate case of very weak solutions was settled by the same authors in [27]. Their approach relies on the idea of intrinsic cylinders by DiBenedetto, see [9,8,10]. The heuristic idea is to compensate for the inhomogeneity of the parabolic p-Laplace operator ∂ t u − div (|Du| p−2 Du) by working with cylinders that depend on the size of |Du|. More precisely, for a parameter λ > 0 that is in some sense comparable to |Du|, the idea by DiBenedetto is to consider cylinders of the form Q (λ) The boundary version of the higher integrability result for the parabolic p-Laplacian has been established by Parviainen [32,33], see also Bögelein & Parviainen [2,7] for the higher order case. The required regularity of the boundary is the same as in the case of the elliptic p-Laplacian, i.e. the complement of the domain is assumed to be uniformly pthick. Finally, we note that Adimurthi & Byun [1] proved global higher integrability even for very weak solutions of parabolic p-Laplace equations.
Even after the case of the parabolic p-Laplace equation had been quite well understood, the corresponding question for porous medium type equations stayed open for a long time. This case turned out to pose additional challenges, which stem from the fact that the differential operator ∂ t u − ∆u m = ∂ t u − mdiv(u m−1 Du) can degenerate depending on the size of u, and not on the size of the gradient as for the parabolic p-Laplace. This type of degeneracy makes it much more involved to derive gradient estimates, because both the size of the solution and of the gradient have to be taken into account. In particular, it is natural to work with intrinsic cylinders of the type where θ m corresponds to 1 u m . The construction of a family of such intrinsic cylinders that is suitable for the derivation of gradient estimates has first been established by Gianazza and the third author in [15], using an idea from [34]. The article [15] contains the first result on higher integrability of the gradient for porous medium type equations and opened the path to further results in this direction. The higher integrability result was already extended to systems in [4], to singular porous medium equations and systems, i.e. the case m < 1, in [16,6], and to a doubly nonlinear system in [3]. All of the mentioned results are restricted to the interior case. The present article is devoted to the question whether the higher integrability of the gradient can be extended up to the boundary. As to be expected from the p-Laplacian case, we have to assume that the complement of the domain is uniformly 2-thick. However, it turns out that we need a further assumption on the domain in the case of the porous medium equation. The additional problem stems from the fact that the degeneracy of the porous medium equation depends on the values of the solution itself rather than on the gradient. This means that close to the boundary, the degeneracy also depends on the value of the boundary values. In order to rebalance this nonlinearity with the help of intrinsic cylinders of the type (1.1), we need to estimate the difference of the boundary values and the constant θ by means of a suitable Poincaré inequality on cylinders centred on the boundary. In order to obtain such an inequality for arbitrary boundary data, we have to restrict ourselves to Sobolev extension domains. The exact assumptions will be given in the following section.
with u : Ω T → R N , where A : Ω T × R N × R N n → R N n is a Carathéodory function satisfying A(x, t, u, ζ) · ζ ≥ ν|ζ| 2 |A(x, t, u, ζ)| ≤ L|ζ| (1.3) for a.e. (x, t) ∈ Ω T and any (u, ζ) ∈ R n × R N n . Note that for u ∈ R N we used the short hand notation u α = |u| α−1 u for α > 0, where we interpret u α as zero if u is zero. For the inhomogeneity F : Ω T → R N n we assume that and for the boundary datum g : Ω T → R N we suppose that g m ∈ L 2+ε 0, T ; W 1,2+ε (Ω, R N ) , g ∈ C 0 [0, T ], L m+1 (Ω, R N ) , for some ε > 0. We consider weak solutions in the following sense.
is called a global weak solution to the Cauchy-Dirichlet problem holds true for every test-function ϕ ∈ C ∞ 0 (Ω T , R N ) and, moreover (u m − g m )(·, t) ∈ W 1,2 0 (Ω, R N ) for almost every t ∈ (0, T ) and for a given function g satisfying (1.5).
In order to state our assumptions on the boundary of the domain, we recall the following two definitions. The first one is already familiar from corresponding results for p-Laplace equations.
for all x o ∈ E and for all 0 < < o .
For the treatment of the porous medium equation, we rely on a suitable Poincaré inequality for the boundary values, see Lemma 4.3. In order to achieve our main result for arbitrary boundary values, we need to assume that Ω is a Sobolev extension domain in the following sense. Definition 1.3. A domain Ω ⊂ R n is called a W 1,p -extension domain if there exists a linear operator E : W 1,p (Ω) → W 1,p (R n ) such that Eu(x) = u(x) for a.e. x ∈ Ω and for any u ∈ W 1,p (Ω) and a constant c E ∈ R ≥0 .
In [22] it was shown that every W 1,p -extension domain satisfies the measure density condition, i.e. there exists α > 0 such that for all x o ∈ Ω and 0 < ≤ 1 This allows us to formulate the main result of our paper. In order to state the local estimate, we consider parabolic cylinders of the form

Moreover, for any parabolic cylinder
where we abbreviated The constant ε o depends at most on m, n, N, ν, L, µ, o , and α, and c depends on the same data and additionally on c E . Here, the parameters µ, o are introduced in Definition 2.7 with p = 2, c E is the constant from Definition 1.3 with p = 2 + ε and α is given by (1.9). Remark 1.5. A close inspection of the proof shows that the constants in the preceding theorem actually depend continuously on m > 1 and remain bounded when m ↓ 1.

1.2.
Technical novelties and plan of the paper. It has been observed by Gianazza and the third author in [15] that higher integrability in the interior of the domain can be derived by working with cylinders Q (θ) that are intrinsic in the sense More precisely, a key step in [15] is the construction of a suitable system of cylinders that are sub-intrinsic in the sense that the above integral is bounded from above by θ 2m . Close to the lateral boundary, however, it turns out that the suitable choice of the scaling parameter θ has to depend on the boundary values as well. In the boundary situation, we work with cylinders that satisfy a coupling of the type Both of these coupling conditions have to be taken into account for the construction of a system of sub-intrinsic cylinders as in [15]. In fact, when considering a point z o close to the lateral boundary, it is not clear a priori if the mentioned construction yields a cylinder for which the doubled cylinder Q (θ) 2 (z o ) touches the boundary or not. This is the reason why both the interior scaling (1.11) and the boundary scaling (1.12) enter in the construction of the cylinders, cf. Section 6.2. As a matter of course, the derivation of the desired reverse Hölder inequalities on these cylinders requires a much more extensive case-bycase analysis than in the interior case.
At the initial boundary, we use an extension argument in order to avoid the occurrence of a third type of coupling condition. More precisely, we extend the solution by the reflected boundary values, cf. (3.2) below. Then we use a scaling as in (1.11) with u replaced by its extension. This enables us to treat the initial boundary case with a coupling condition analogous to the interior.
This article is organized as follows. In the preliminary Section 2, we collect some technical tools that will be crucial for the proof. In Section 3 , we derive suitable Caccioppoli type estimates and Section 4 is devoted to Sobolev-Poincaré type inequalities for the solutions. Both estimates are combined in Section 5 to establish reverse Hölder type inequalities on sub-intrinsic cylinders. Each of the three last-mentioned sections is subdivided into one subsection that is concerned with the case of the lateral boundary and another one that deals with the initial boundary. Moreover, for the derivation of the reverse Hölder inequality, we have to consider two different types of coupling conditions for the sub-intrinsic cylinders that can be understood as the non-degenerate case (see (5.1) for the lateral boundary and (5.4) for the initial boundary) and the degenerate case (cf. (5.3) and (5.6), respectively). The final Section 6 contains the construction of a suitable system of cylinders, which can be shown to satisfy one of the mentioned coupling conditions that lead to a reverse Hölder inequality. By a Vitali type covering argument, the reverse Hölder estimates on the cylinders can be extended to estimates on the super-level sets. Then, a standard Fubini type argument yields the result.
In the case θ = 1 we use the shorter notation Q (z o ) := Q (1) (z o ). From the definition of the cylinders it becomes clear that the parabolic dimension associated to our problem is d := n + 1 + 1 m . Moreover, we will use the notations For the mean value of a function f ∈ L 1 (A) over a set A ⊂ R k of finite positive measure we write (f ) A := − A f dx, and for a function v ∈ L 1 (Ω T ), we abbreviate moreover where t ∈ [0, T ]. Finally, we define the boundary term as b[u m , a m ] := m m+1 |a| m+1 − |u| m+1 − u · a m − u m .

Auxiliary tools.
In order to prove energy estimates we have to use a mollification in time. For this purpose we define for v ∈ L 1 (Ω T , For the basic properties of the mollification · h we refer to [ Lemma 2.1. Let α > 1. There exists a constant c = c(α) such that for any u, a ∈ R N the following holds true: There exists a constant c = c(m) such that for every u, a ∈ R N we have There exists a constant c = c(m) such that for any bounded A ∈ R n , any u ∈ L m+1 (A, R N ), and any a ∈ R N there holds The proof of the following lemma can be found in [3,Lemma 3.5], see also [11,Lemma 6.2] for an earlier version in a special case.
Lemma 2.4. Let p ≥ 1 and α ≥ 1 p . Then there exists a constant c = c(α, p) such that for any bounded sets of positive measure satisfying A ⊂ B ⊂ R k , k ∈ N and any u ∈ L αp (B, R N ) and constant a ∈ R N there holds Finally, we state a well-known absorption Lemma, that can be found in [21, Lemma 6.1] for instance. Lemma 2.5. Let 0 < ϑ < 1, A, C ≥ 0 and α, β > 0. Then, there exists a constant c = c(β, ϑ) such that there holds: For any 0 < r < and any nonnegative bounded

2.3.
Variational p-capacity. Let 1 < p < ∞ and D ⊂ R n be an open set. The variational p-capacity of a compact set C ⊂ D is defined by where the infimum is taken over all functions f ∈ C ∞ 0 (D) such that f ≡ 1 in C. In order to define the variational p-capacity of an open set U ⊂ E, we are taking the supremum over the capacities of compact sets contained in U . The variational p-capacity for an arbitrary set E is defined by taking the infimum over the capacities of the open sets containing E. The capacity of a ball is For more details we refer to [12,Ch. 4] or [24,Ch. 2].
At this point we introduce the uniform capacity density condition, which is essential for proving a boundary version of a Sobolev-Poincaré type inequality, where we note that this condition is essentially sharp in the context of higher integrability. For the elliptic setting we see [25], whereas the equations of parabolic p-Laplacian type were treated in [29].
We recall the definition of uniform p-thickness introduced in Definition 1.2. The following consequences of this property are well-known, see e.g. [32,Lemma 3.8].
Lemma 2.6. Let Ω ⊂ R n be a bounded open set and assume that R n \ Ω is uniformly p-thick. Choose y ∈ Ω such that B 4 /3 (y) \ Ω = ∅. Then there exists a constantμ = µ(n, µ, o , p) > 0 such that Lemma 2.7. If a compact set E is uniformly p-thick, then E is uniformly ϑ-thick for any ϑ ≥ p.
The next theorem shows that a uniformly p-thick set has a self-improving property, see [30].
Before we proceed, let us recall that u ∈ W 1,p (Ω) is called p-quasicontinuous if for each ε > 0 there exists an open set U ⊂ Ω ⊂ B R such that cap p (U, B 2R ) ≤ ε and the restriction of u to the set Ω \ U is finite valued and continuous. Note that every function u ∈ W 1,p (Ω) has a p-quasicontinuous representative. A proof of the next lemma can be found in [23]. Lemma 2.9. Let B (x o ) be a ball in R n and fix a q-quasicontinuous representative of u ∈ W 1,q (B (x o )). Denote Then there exists a constant c = c(n, q) > 0 such that The following Lemma can be found for instance in [32,Lemma 3.13].
Then, forq ∈ [q, q * ] with q * = nq n−q there exists a constant c = c(n, q) > 0 such that

ENERGY ESTIMATES
In this section, we will prove energy estimates that are required to prove a reverse Hölder inequality.

3.1.
Estimates near the lateral boundary. We begin with a Caccioppoli type estimate at the lateral boundary.
Proof. The mollified version of the system (1.6) reads as as testing function in the mollified weak formulation (3.1). We start with the parabolic part of the equation and estimatë where we also used that We are now able to pass to the limit h ↓ 0 in the right-hand side of the previous estimate and obtain Now, we pass to the limit ε ↓ 0 and obtain for the first term where we note that the integral at the time t = 0 vanishes by assumption (1.7) in connection with Lemma 2.2. The second term can be estimated as follows whereas the third term is estimated with the help of Young's inequality and Lemma 2.1 (i) Next we will treat the diffusion term. After passing to the limit h ↓ 0 we use the ellipticity and growth condition (1.3) and Young's inequality and hence we arrive aẗ for a constant depending on m, ν and L. Let us now consider the right hand side in (3.1). Note that the second term vanishes in the limit h ↓ 0, since which follows from (1.7), Lemma 2.1(ii) and Hölder's inequality. In the term containing F we also pass to the limit h ↓ 0 and use Young's inequality afterwards to obtain We combine all these estimates and pass to the limit ε ↓ 0. This showŝ in the first term on the left-hand side and then pass to the limit in the second term. This proves the lemma.

3.2.
Estimates near the initial boundary and in the interior. Up next we prove the corresponding Caccioppoli estimate near the initial boundary Ω × {0}. For the initial datum we use the abbreviation We do not impose an additional regularity assumption on the initial datum except from g 0 ∈ L m+1 (Ω, R N ). However, we exploit the fact that there is an extension g : Ω T → R N with g(·, 0) = g 0 and g m ∈ L 2+ε (0, T ; W 1,2+ε (Ω, R N )) as well as ∂ t g m ∈ L m(2+ε) At the initial boundary, we begin with a Caccioppoli type estimate for the extended function We note that the following result also contains the interior case 0 < ≤ 1 and θ > 0, the following holds. For every r ∈ [ /2, ) and every a ∈ R N , the energy estimate holds true, whereû is defined according to (3.2).
Proof. We start with arguments similar to the proof of Lemma 3.1. We consider the mollified version (3.1) of the equation and use now the test-function with η, ζ, and ψ ε defined as in Lemma 3.1 and g m replaced by a m . Observe that ∂ t a m = 0 and Da m = 0. For the parabolic part we obtain By first passing to the limit h ↓ 0, then ε ↓ 0 and using the same estimates as in Lemma 3.1 we arrive at Here we also used the fact that which follows from Lemma 2.1 (i) and assumption (1.7). The diffusion term and the term containing F are treated exactly in the same way as in Lemma 3.1 with a m instead of g m (with obvious simplifications as Da m = 0). The second integral on the right-hand side of the mollified equation (3.1) vanishes in the limit h ↓ 0 because of assumption (1.7). By combining these estimates we obtain the bound It remains to estimate the last integral. We start with the observation that two applications of Lemma 2.
and moreover, we have the identity This enables us to estimate where we have abbreviated Q (θ) ,− := Q (θ) ∩ {t < 0}. Next, we use Young's inequality, the facts ≤ 1 andû(t) = g(−t) for t < 0, as well as Lemmas 2.1 and 2.2, with the result Plugging this estimate into (3.3), we arrive at It remains to estimate the terms on the left-hand side for negative times t ∈ Λ (θ) For the estimate of the first term, we observe that |Λ To the remaining term, we apply Young's inequality and Fubini's theorem, which leads to the estimate In the last step, we used Lemma 2.1 (i), the fact ≤ 1 and the definition ofû. Moreover, from the definition ofû, we immediately obtain the estimate Combining the estimates (3.5) and (3.6) with (3.4), we deduce the claim.
Next we prove a lemma that allows us to compare slice-wise values of the solution between the initial time and any given point of time. This type of lemma is termed gluing lemma, and we will use it later in the proof of a Sobolev-type inequality near the initial boundary.
We start by recalling the gluing lemma from the interior case, see [4,Lemma 3.2]. By applying this result to the cylinder Q (θ) ,+ (z o ) and using initial condition (1.7) in the case t = 0, we infer the following lemma.
with a constant c = c(L).
We extend this result to a version adapted to the initial boundary.
Proof. Throughout the proof, we omit the reference to the center z o in the notation. We choose the radiusˆ ∈ [ 2 , ] that is provided by Lemma 3.3. We follow different strategies depending on whether the considered times are positive or negative. In the case t, τ ≥ 0, we combine Lemma 3.3 with Lemma 2.1 (ii) to obtain Next, we consider the case t, τ ≤ 0, in which we can estimate We apply Lemma 2.4, the definition ofû and Poincaré's inequality in order to estimate Obviously, the same estimate holds true for τ in place of t. From the two preceding estimates, we deduce for any τ, t ∈ Λ (θ) with τ, t ≤ 0. It remains to consider the case t < 0 < τ . In this case we combine the estimates (3.7) with t = 0 and (3.9) with τ = 0 and deduce Using (3.8) with τ = 0, the term Θ m−1 τ,0 can be bounded as follows.
Plugging this into the preceding estimate and applying Young's inequality with exponents m m−1 and m, we arrive at In view of Estimates (3.7) and (3.9), this estimates holds in any case, i.e. for arbitrary times t, τ ∈ Λ (θ) . We multiply the preceding estimate with where in the last step we applied Young's inequality, once with exponents m m−1 and m and a second time with m 2 (m−1) 2 and m 2 2m−1 . We re-absorb the first term of the right-hand side into the left and take the mth root of both sides. This yields the asserted estimate.

SOBOLEV POINCARÉ TYPE INEQUALITIES
4.1. Estimates near the lateral boundary. The next lemma is an adoption of Lemma 4.2 of [7]. However, for the sake of completeness we will state a proof.
We can extend u m − g m outside of Ω T by zero (still denoted in the same way) and define for fixed t ∈ (t o − s, t o + s) ∩ (0, T ) the set Using Lemma 2.9 showŝ , with a constant c depending only on n, N, ϑ. Since R n \ Ω is uniformly ϑ-thick, Lemma 2.6 and (2.1) imply whereμ =μ(n, µ, o , ϑ). Combining the previous estimates leads tô Finally, integrating this inequality with respect to t over (t o − s, t o + s) ∩ (0, T ) finishes the proof of the Lemma.
Next, we are going to prove a different version of a Sobolev-type inequality. To this end, we assume that the boundary values are extended to a function g ∈ L 2+ε (0, T ; W 1,2+ε (Ω, R N )), which is possible since Ω is an extension domain. Moreover, we extend the solution and the boundary values across the initial boundary by letting Note thatû m −ĝ m = 0 outside of Ω T . For the proof of the Sobolev-type inequality, we assume that the cylinders Q (θ) satisfy the sub-intrinsic scaling We observe that (4.2) implies that holds true. . Then there exists q = q(n, µ) ∈ (1, 2) such that for every ε ∈ (0, 1) where c = c(n, m, N, µ, o ).
Proof. To shorten the notation, we will omit z o as the reference point for the cylinder. Note that the condition dist(B (x o ), ∂Ω) = 0 implies that B 4 3 (x o ) \ Ω = ∅. With a similar argument as in the proof of Lemma 4.1, where we use Lemma 2.10 instead of Lemma 2.9, we obtain an exponent ϑ = ϑ(n, µ) ∈ (1, 2) so that for every q ∈ [ϑ, 2) we have 1 for a constant c = c(n, N, µ, o , q). For α = α(m, q) ∈ (0, 2) to be chosen later we estimate with the help of Lemma 2.2 (iii), Hölder's inequality and the sub-intrinsic scaling (4.3) where we used the short-hand notation for exponents and and p , r are the Hölder conjugates of p and r. Let us note that p, r > 1 holds true when we choose α suitably, as we do below. Next, we apply Hölder's inequality and then estimate (4.4). In this way, we deduce At this point we are choosing α such that Next, we observe that we obtain in the limit q ↑ 2 that α → 0, p → 1 and r → n 2 . Therefore, we can choose q close to 2 such that p > 1, r > 1, and α 2 p r < 1 hold true and we are able to use Hölder's and Young's inequalities to deduce Noting that 2 (2−α)p = 2 q finishes the proof.
Next, we are going to prove a Poincaré inequality for the boundary function g that will be very useful in the course of the paper. This is the point in the proof at which the Sobolev extension property of the domain is crucial. We recall that the extension property in particular implies the measure density condition (1.9), which in turn implies the lower bound for a constant c depending only on m, n and α.
Proof. For simplicity we omit the center of the cylinder in the notation. Adding and subtracting the slice-wise mean value integral leads to Using Lemma 2.4, the measure density condition (1.9) and the Sobolev-inequality shows for the first term for a constant c = c(α, m, n). In order to treat the second term we may assume that t > 0 because otherwise we use the identityĝ(t) =ĝ(−t). This allows us to estimate This proves the following estimate Using the sub-intrinsic scaling of the cylinders and m > 1, we obtain This in connection with the last estimate shows We multiply this estimate with , and take both sides to the power m 2m−1 , which leads to the bound Absorbing the first term on the right-hand side, using the measure density condition (4.6) and applying Hölder's inequality finishes the proof of the lemma.

Estimates near the initial boundary.
Here we prove a Sobolev-type inequality near the initial boundary. First we prove one auxiliary lemma, since we cannot use the Sobolev inequality in both space and time directions simultaneously. That is due to the fact that less regularity is assumed in the time direction. That is why we first use the gluing lemma to treat the time direction and then use the Sobolev inequality slice-wise in space. In this section we assume that the considered cylinders ,+ (zo) Proof. Letˆ ∈ [ 2 , ] be the radius in Lemma 3.4. For simplicity we omit the reference point z o in the notation. We start by decomposing The first integral we can estimate by using Lemma 2.4 slice-wise and the factˆ ∈ [ 2 , ] to obtain in which c = c(n, m). For the second term II we use Gluing Lemma 3.4, Hölder's inequality and the sub-intrinsic scaling (4.8) such that we have For the third term, Hölder's inequality and the estimate for I imply which completes the proof. Now we are able to prove a suitable Sobolev-type inequality near the initial boundary.

REVERSE HÖLDER INEQUALITIES
In this section we will prove reverse Hölder inequalities. Since the construction of our cylinders does not ensure that we always have intrinsic coupling, we have to distinguish between two cases here. Additionally, we have to treat the lateral boundary in a different way than the initial boundary.
Proof. Let 0 < ≤ r < s ≤ 2 . To shorten the notation, we will again omit the reference point z o . Utilizing Lemma 3.1 shows For the second term we use the intrinsic coupling (5.1), noting that

and Lemma 2.2 to obtain
In order to estimate the second term on the right-hand side, we first use the Poincaré inequality (4.7) to obtain where we abbreviated Using Lemma 2.2 and Young's inequality, we further estimate Inserting the estimates for the terms I and II, and using Lemma 4.2 shows for every ε ∈ (0, 1) that We are now in position to apply Lemma 2.5 and obtain This finishes the proof of the Lemma.
Proof. We consider again radii r, s > 0 with ≤ r < s ≤ 2 and take the Caccioppoli inequality from Lemma 3.1 as starting point. We use the same short-hand notation as in the proof of Lemma 5.1. The first term in (5.2) can be estimated in the same way as before, whereas the second term will be treated in a different way. By using Young's inequality and Lemma 2.2 (ii) we obtain Using the intrinsic coupling (5.3) allows us to absorb the term involving Du m and moreover, exploiting Lemma 4.2 leads to Proceeding as in the proof of Lemma 5.1 completes the proof.

5.2.
The initial boundary. Our next goal is the proof of reverse Hölder inequalities at the initial boundary. Again, we have to distinguish between two cases.
for some 0 < ≤ 1 and θ ≥ 1 we have the following reverse Hölder inequality for a constant c = c(n, m, ν, L) and for q := 2n d < 2. Proof. We omit the reference point z o in notation, and consider radii ≤ r < s ≤ 2 . From the Caccioppoli estimate in Lemma 3.2 we obtain Choosing ε = 1 2cR 4m m+1 r,s and using the Iteration Lemma 2.5 in order to reabsorb the supterm into the left-hand side, we deduce the asserted estimate.
Next we prove the reverse Hölder inequality in the degenerate case.
for some 0 < ≤ 1 and θ ≥ 1 we have the following reverse Hölder inequality for a constant c = c(n, m, ν, L, K) and for q := 2n d < 2. Proof. Similarly as in the proof of the preceding lemma, we start with estimate (5.5). The term I is estimated in the same way as before, but now we estimate II by Using assumption (5.6) 2 to bound the first term and Lemma 4.5 for the estimate of the second, we deduce Choosing first δ and then ε small in the form δ = 1 4K and ε = ∈ Ω T ∪ ∂ par Ω T . Since the center will be fixed throughout this section, we will simply write Q := Q (y o , τ o ) for > 0. We fix a specific extension of the boundary values in order to derive an estimate on the cylinder Q R . To this end, we choose a standard cut-off function η ∈ C ∞ 0 (B 8R , [0, 1]) with η ≡ 1 in B 4R and |Dη| ≤ 1 R in B 8R . We assume that the extension of the boundary values is given by g m = E(ηg m ) on Q 8R,+ \ Ω T , where the extension operator E from Definition 1.3 is applied separately on each time slice. Then, for each fixed time t ∈ Λ 4R (τ o ) ∩ (0, T ) we have the estimateŝ In the case n > 2, we use Hölder's inequality and Sobolev's embedding to infer where the last estimate follows from (6.1). In dimension n = 2, we use the Sobolev where we used the measure density property (1.9) and (6.1) in the last step. From the three preceding estimates, we deduce the bound Using the extension of the boundary values specified above, we now define whereû andĝ are defined in (4.1) and We use (6.2) and Lemma 3.1 with θ = 1 in order to estimate .
For the estimates we also used the measure density condition (1.9), which implies |Q 8R ∩ Ω T | ≥ c|Q 8R |.

6.2.
Construction of a non-uniform system of cylinders. The following construction is inspired by the one in [4,15,34]. However, the boundary case becomes much more involved due to the fact that the notion of intrinsic cylinder is different at the lateral boundary compared to the interior of the domain. The transition between both cases requires additional carefulness. For whenever ∈ (0, R] and θ ≥ 1.
For ∈ (0, R], we define the parameterθ ≡θ zo; bỹ Observe thatθ is well defined, since the integral condition is satisfied for some θ ≥ λ o . This follows from the fact that in the limit θ → ∞, the integral on the left-hand side converges to zero, while the right-hand side blows up. Note that we can rewrite the condition for the integral in the definition ofθ as By the very definition ofθ we either havẽ In any case we haveθ R ≥ λ o ≥ 1. If λ o <θ R and R ≥ d o then we obtaiñ where we used the fact that Q (θ R ) R ⊂ Q 4R . If R < d o we argue similarly. In any case, we obtain Next, we prove that the functionθ is piecewise continuous: Proof. Without loss of generality, we may assume that d o ∈ (0, R). If ∈ (0, d o ) the proof works as in [4,Section 6.1]. If ∈ [d o , R] the idea still remains the same, but we will present the proof for convenience. Therefore, we consider ∈ [d o , R] and ε > 0 and define θ + :=θ + ε. Then there exists δ = δ(ε, ) > 0 such that with |r − | < δ. This can be verified by observing that the strict inequality holds true if r = and that both sides are continuous with respect to the radius. This showsθ r ≤ θ + =θ + ε if |r − | < δ. To prove the reverse inequality we set θ − :=θ − ε. If θ − ≤ λ o the desired estimate follows directly from the construction. In the other case we obtain Then, by Lemma 6.1 and the construction, the map → θ is continuous and monotonically decreasing. This construction can be considered as a rising sun construction (see Figure 1).
(ii) For any s ∈ ( , R] we have θ ≤ s d+2 m+1 θ s . (iii) For any 0 < ≤ R there holds Illustration of the rising sun construction.
Proof. (i) By definition we haveθ s ≤ θ so that Q (iii) Combining (ii) for s = R with estimate (6.5) yields the claim, since θ R =θ R .
We note that due to the monotonicity of the map → θ ,zo the above constructed cylinders are nested in the sense that However, these cylinders in general only fulfill the sub-intrinsic coupling condition from Lemma 6.2 (i).
6.3. Covering property. Next, we will present a Vitali type covering property for the above constructed cylinders. Using the just established bounds from Lemma 6. Proof. For j ∈ N define a sub-collection of F as Next we construct a countable collection of disjoint cylinders G ⊂ F inductively as follows. Let G 1 be a maximal disjoint collection of cylinders in F 1 . Observe that the measure of every cylinder in G 1 is bounded from below by Lemma 6.2 (iii), which implies that G 1 is finite. For k ≥ 2, define G k to be any maximal sub-collection of Collection G k is again finite, and we can define Now G is a countable disjoint sub-collection of F. To conclude the result we show that for any Q ∈ F there exists a cylinder Q * ∈ G such that Q ⊂Q * .
Fix Q = Q in whichĉ is the constant from the Vitali covering lemma 6.3. For radii s with we obtain by using Lemma 6.2 (iii) . By absolute continuity of the integral and the continuity of s → θ s , there exists a maximal radius zo ∈ (0, R2−R1 The maximality of zo implies for any s ∈ ( zo , R]. 6.5. Reverse Hölder inequalities. As before we consider z o ∈ E(R 1 , λ) and abbreviate θ zo ≡ θ zo, zo . From now on we denote the exponent q < 2 as the maximum of the Sobolev exponents q in Lemma 4.2 and 2n d in Lemma 4.5. We distinguish between the non-degenerate and the degenerate case, which correspond to the cases˜ zo ≤ 2 zo and zo > 2 zo . 6.5.1. The non-degenerate case˜ zo ≤ 2 zo . In this case, we note that the cylinder Q (θ zo ) zo (z o ) is intrinsic since˜ zo < R. We first consider the boundary case˜ zo ≥ d o . Lemma 6.2 (i) and the fact that Q 2˜ zo (z o ) intersects or touches the lateral boundary. Hence, we are in position to use Lemma 5.1 on this cylinder to obtain which is the desired reverse Hölder inequality. In the remaining case˜ zo < d o , we are either in the interior case (if Q (θ zo ) 2˜ zo Ω T ) or we might intersect with the initial boundary.
Therefore, Lemma 6.2 (i) and the fact that Q On the other hand, inequality (6.13), the monotonicity of the mapping → θ and Lemma 6.2 (ii) imply that 6.7. Proof of the gradient estimate. With estimate (6.18) on super-level sets and using Fubini-type arguments we are finally able to prove the higher integrability for the gradient of the solution. In order to ensure that quantities we end up re-absorbing are finite we consider truncations. For k > λ 1 we define Du m k := min Du m , k m , and the corresponding super-level set as E k (r, λ) := z ∈ Q r ∩ Ω T : Du m k (z) > λ m . With this notation and estimate (6.18) we havë for k > λ 1 . Here we exploited the facts that Du m k ≤ Du m a.e., E k (r, λ) = E(r, λ) if k > λ and E k (r, λ) = ∅ if k ≤ λ.
Let ε ∈ (0, 1]. We multiply the inequality above by λ εm−1 and integrate over the interval (λ 1 , ∞). By using Fubini's theorem, on the left-hand side we havê For the first term on the right-hand side we obtain Combining the estimates and multiplying by εm we obtain Du m q χ Ω T dxdt.
Adding the two previous estimates we deducë Du m 2 χ Ω T dxdt + c¨Q and assume that ε ≤ ε o . Now λ ε 1 = (ηBλ o ) ε ≤ Bλ ε o since η ≤ 1, B > 1 and ε < 1. We obtain Du m 2 χ Ω T dxdt + c¨Q 2R,+ G 2+ε dxdt, for any R 1 , R 2 satisfying R ≤ R 1 < R 2 ≤ 2R. By using Iteration Lemma 2.5 we can re-absorb the first term into the left-hand side. Then by passing to the limit k → ∞ and using Fatou's Lemma we can concludë Estimating λ o by means of (6.3) and the last integral by (6.2) proves the estimatë with c = c(m, n, N, ν, L, α, µ, o , c E ). Finally, we note that we can replace the integrals over Q 8R by integrals over Q 2R by a standard covering argument. More precisely, we cover the cylinder Q R by smaller cylinders Q R/8 (z i ) with centers z i ∈ Q R , apply the preceding estimate on each of the smaller cylinders and sum up the resulting inequalities. This procedure leads to the asserted estimate (1.10). The local estimate implies |Du m | ∈ L 2+ε (Ω τ ) for every τ < T . However, we can assume that the solution is given on the larger cylinder Ω 2T by reflecting the boundary values across the time slice Ω × {T } and solving a Cauchy-Dirichlet problem on Ω × [T, 2T ]. Applying the preceding result on Ω 2T , we deduce the remaining assertion |Du m | ∈ L 2+ε (Ω T ). This completes the proof of Theorem 1.4.