\betagin{document}

\betagin{frontmatter}

\title{Boundary higher integrability for very weak solutions of quasilinear parabolic equations.}

\author[myaddress]{Karthik Adimurthi\corref{mycorrespondingauthor}\tnoteref{thanksfirstauthor}}
\cortext[mycorrespondingauthor]{Corresponding author}
\ead{karthikaditi@gmail.com and kadimurthi@snu.ac.kr}
\tnotetext[thanksfirstauthor]{Supported by the National Research Foundation of Korea grant NRF-2015R1A2A1A15053024.}

\author[myaddress,myaddresstwo]{Sun-Sig Byun\tnoteref{thankssecondauthor}}
\ead{byun@snu.ac.kr}
\tnotetext[thankssecondauthor]{Supported by the National Research Foundation of Korea grant  NRF-2015R1A4A1041675. }

\address[myaddress]{Department of Mathematical Sciences, Seoul National University, Seoul 08826, Korea.}
\address[myaddresstwo]{Research Institute of Mathematics, Seoul National University, Seoul 08826, Korea.}

\betagin{abstract}
We prove boundary higher integrability for the (spatial) gradient of \emph{very weak} solutions of quasilinear parabolic equations of the form $$u_t - \text{div}\,\mathcal{A}(x,t, \nabla u)=0 \quad \text{on} \ \Omegaega \times \mathbb{R},$$ where the non-linear structure $\text{div}\,\mathcal{A}(x, t,\nabla u)$ is modelled after the  $p$-Laplace operator.    To this end, we prove that the gradients satisfy a reverse H\"older inequality near the boundary. In order to do this, we construct a suitable test function which is Lipschitz continuous and  preserves the boundary values. \emph{These results are new even for linear parabolic equations on domains with smooth boundary and make no assumptions on the smoothness of $\mathcal{A}(x,t,\nabla u)$}. These results are also applicable for systems as well as higher order parabolic  equations.

\vspace*{0.5cm}
\noindent {\bf R\'esum\'e}\\

Nous montrons l'int\'egrabilit\'e des limites pour le gradient (spatial) de solutions tr\`es faibles d'\'equations paraboliques quasi-lin\'eaires de la forme
$$
u_t-\text{div}\mathcal{A}(x,t,\nabla u)=0 \quad \text{ sur } \Omegaega\times \mathbb{R},
$$
o\`u la structure non lin\'eaire $\text{div}\mathcal{A}(x,t,\nabla u)$ est model\'ee d'apr\`es l'op\'erateur
$p$-Laplace. Nous prouvons que les gradients satisfont l'in\'egalit\'e de H\"older inverse pr\`es du bord. Pour ce faire, nous construisons une fonction test appropri\'ee qui est lipschitzienne et pr\'eserve les valeurs limites. Ces r\'esultats sont nouveaux, m\^eme pour les \'equations  paraboliques lin\'eaires \`a fronti\`ere lisse et ne font aucune hypoth\`ese sur la r\'egularit\'e de $\text{div}\mathcal{A}(x,t,\nabla u)$ et sont \'egalement applicables aux syst\`emes ainsi qu'aux \'equations paraboliques d'ordre sup\'erieur. 
\end{abstract}

\betagin{keyword}
quasilinear parabolic equations, boundary higher integrability, very weak solutions
\MSC[2010] 35K10\sep  99-00\sep 35K59\sep 35K65\sep 35K67.
\end{keyword}

\end{frontmatter}

\section{Introduction}
In this paper, we consider the boundary regularity  of quasilinear parabolic partial differential equations of the form 
\betagin{equation}
\lambdabel{main}
 \left\{
 \betagin{array}{ll}
u_t - \dv  \mathcal{A}(x,t,\nabla u) = 0 & \text{on} \ \Omega \times (-T,T) \\
u = 0 & \text{on} \ \partialrtial \Omega \times (-T,T),
 \end{array}
\right.
\end{equation}
modelled after the well studied $p$-Laplacian operator in a bounded domain $\Omega\subset \mathbb{R}^n$,  potentially with non-smooth boundary $\partial \Omega$.

Weak solutions to \eqref{main} are in the space $u \in L^2(-T,T; L^2(\Omega)) \cap L^p(-T,T; W_0^{1,p}(\Omega))$ which allows one to use $u$ as a test function. But from the definition of weak solution, we see that the expression (see Definition \ref{very_weak_solution}) makes sense if we only assume $u \in L^2(-T,T; L^2(\Omega)) \cap L^{s}(-T,T; W_0^{1,s}(\Omega))$ for some $s > p-1$. But under this milder notion of solution called \emph{very weak solution}, we lose the ability to use $u$ as a test function. This difficulty was overcome in \cite{KL} in the interior case by constructing a suitable test function after modifying $u$ on a \emph{bad} set. The result was further extended to higher order systems in \cite{Bog} {\color{black}(see also \cite[Chapter 3]{bogelein2013regularity} where a more general form of Lipschitz truncation method was developed)}. In this paper, we show the boundary higher integrability for the gradient of very weak solutions to \eqref{main} in the sense of \cite{KL}. \emph{The results presented here are new even for linear parabolic equations on smooth domains.}

Regarding the assumption on the structure of the boundary of the domain, we assume that the complement of the domain $\Omega^c$ satisfies a uniform capacity density condition (see Definition \ref{p_thick_domain}) . This condition is essentially optimal for our main results, as this condition is needed to ensure that Poincar\'e's inequality is applicable (see \cite[Chapter 11]{AH} for example).  The proofs are based on the Caccioppoli and Sobolev-Poincar\'e type inequalities as well as on the careful analysis of the associated strong Maximal function. As with regularity results concerning quasilinear parabolic equations, we make use of intrinsic scaling and covering techniques. 

The first nonlinear parabolic higher integrability results apparently date back to a 1982 paper \cite{Gstr} in which the author studied the local higher integrability for systems of parabolic equations with quadratic growth conditions. However for more general quasilinear systems, the problem remained open until the papers \cite{KL, KL1} in which interior higher integrability was shown for both \emph{weak} and \emph{very weak} solutions to equations of the form considered in \eqref{main}.  For higher order quasilinear parabolic systems, interior higher integrability for very weak solutions was proved in \cite{Bog}. The boundary higher integrability for \emph{weak} solutions under the condition that $\Omega^c$ satisfies a uniform capacity density condition was done in \cite{Par, Mik1}.  

We shall make precise all the notations in the following section, but let us first state the main theorem:
\betagin{theorem}
 \lambdabel{main_theorem}
 Let $n\geq 2$ and $\frac{2n}{n+2} < p<\infty$ and $\Omega$ be a bounded domain whose complement is uniformly $p$-thick with constants $r_0,b_0$, then there exists $\beta_0 = \beta_0(n,p,\Lambda_0,\Lambda_1,b_0) \in (0,1)$ such that the following holds: For any $\beta \in (0,\beta_0]$ and any very weak solution $u \in L^2(-T,T; L^2(\Omega)) \cap L^{p-\beta}(-T,T; W_0^{1,p-\beta}(\Omega))$ solving \eqref{main} in $\Omega \times (-T,T)$,  we have the following improved integrability $u \in L^2(-T,T; L^2(\Omega)) \cap L^{p+\beta}(-T,T; W_0^{1,p+\beta}(\Omega))$ along with the following qualitative estimates: for any $Q_{z_0}(\rhoo,s)$ with $z_0 \in (\overline{\Omega} \times (-T,T))$, there holds

\betagin{equation}
\lambdabel{main_est}
  \Yint\tiltlongdash_{Q_{z_0}(\rhoo,s)} |\nabla u|^{p} \ dz \apprle \lbr \Yint\tiltlongdash_{Q_{z_0}(2\rhoo, 2^2s)} \lbr |\nabla u| + \theta\rbr^{p-\beta} \ dz\rbr^{1+\frac{\beta}{d}}  + \Yint\tiltlongdash_{Q_{z_0}(2\rhoo, 2^2 s)} (1 + \theta^p) \ dz. 
 \end{equation}
where $\theta$ is from \eqref{bound_b}, $C = C(n,p,b_0,\Lambda_0,\Lambda_1)$  and 
\betagin{equation}
 \lambdabel{de_d}
 d:= \left\{ \betagin{array}{ll}
              2-\beta & \text{if} \ p \geq 2 \\
              p-\beta - \frac{(2-p)n}{2} & \text{if} \ p<2
             \end{array}\right.
\end{equation}

\end{theorem}

\betagin{remark}
     In Theorem \ref{main_theorem}, we claim $u \in L^2(-T,T; L^2(\Omega)) \cap L^{p+\beta}(-T,T; W_0^{1,p+\beta}(\Omega))$ whereas the estimate in \eqref{main_est} only gives $u \in L^2(-T,T; L^2(\Omega)) \cap L^{p}(-T,T; W_0^{1,p}(\Omega))$. The current paper only proves higher integrability till the natural exponent as in \eqref{main_est}, but to reach $p+\beta$, we can now make use of the results from \cite{Par,Mik1}. 
\end{remark}

The plan of the paper is as follows: In Section \ref{Preliminaries}, we shall collect all the assumptions and some well known, but useful lemmas. In Section \ref{Construction_of_test_function}, we shall give the construction of the test function using Lipschitz truncation and collect its properties. In Section \ref{Caccioppoli_inequality}, we will use results from Section \ref{Construction_of_test_function} to obtain Caccioppoli inequality and Reverse H\"older inequality and finally in Section \ref{main_estimates}, we will provide the proof of Theorem \ref{main_theorem}. 

\section{Preliminaries}\lambdabel{Preliminaries}

In this section, we shall collect some well known theories and estimates that will be used in the latter parts of the paper. 

\subsection{Variational \texorpdfstring{$p$}{pc}-Capacity}
Let $1<p<\infty$, then the \emph{variational $p$-capacity} of a compact set $K \Subset \Omega$ is defined to be 
\betagin{equation*}
 \lambdabel{var_cap_def}
 \cpt(K, \mathbb{R}^n) = \inf \left\{ \int_{\mathbb{R}^n} |\nabla \phii|^p \ dx \ : \ \phii \in C_c^{\infty}(\mathbb{R}^n),\ \lsb{\chii}{K}(x) \leq \phii(x) \leq 1 \right\},
\end{equation*}
where $\lsb{\chii}{K} (x) = 1 $ for $x \in K$ and $\lsb{\chii}{K}(x) = 0$ for $x \notin K$. To define the variational $p$-capacity of an open set $O \subset \mathbb{R}^n$, we take the supremum over the capacities of the compact sets belonging to $O$. The variational $p$-capacity of an arbitrary set $E \subset \mathbb{R}^n$ is defined by taking the infimum over the capacities of the open sets containing $E$. For further details, see \cite{AH,HKM}.

Let us now introduce the capacity density conditions which we later impose on the complement of the domain. For higher integrability results, this condition is essentially sharp (see  \cite[Remark 3.3]{KK94} for the linear elliptic case and \cite{AP1} for the quasilinear elliptic case). 
\betagin{definition}[Uniform $p$-thickness] \lambdabel{p_thick}
Let $\tilde{\Omega}\subset\mathbb{R}^n$ be a bounded domain. We say that the complement $\tilde{\Omega}^c:=\mathbb{R}^n\setminus \tilde{\Omega}$ is 
uniformly $p$-thick for some $1< p \leq n$ with constants  $ \tilde{b}_0,\tilde{r}_0>0$, if  the inequality  
$$ {\rm cap}_{p} ( \overline{B_r(x)} \cap \tilde{\Omega}^c, B_{2r}(x)) \geq  \tilde{b}_0\, {\rm cap}_{p} (\overline{B_{r}(x)}, B_{2r}(x))$$ 
holds for any $x \in \partialrtial \tilde{\Omega}$ and $r\in(0, \tilde{r}_0]$. 
\end{definition}

It is well-known that the class of  domains with uniform $p$-thick complements
is very large. They include all domains with Lipschitz boundaries or even those that satisfy a uniform exterior corkscrew condition, where the latter means that there exist constants $c_0, r_0>0$ such that for all $0<t\leq r_0$ and all $x\in \mathbb{R}^n\setminus \Omega$, there is $y\in B_t(x)$
such that $B_{t/c_0}(y)\subset \mathbb{R}^n\setminus \Omega$.

If we replace the capacity with the Lebesgue measure in Definition \ref{p_thick}, then we obtain a measure density condition. A set $E$ satisfying the measure density condition is uniformly $p$-thick for all $p>1$. If $p>n$, then every non-empty set is uniformly $p$-thick. The following lemma from \cite[Lemma 3.8]{Par} extends the capacity estimate in Definition \ref{p_thick} to make precise the notion of being \emph{uniformly $p$-thick}:
\betagin{lemma}[\cite{Par}]
 Let $\tilde{\Omega}$ be a bounded open set, and suppose that $\mathbb{R}^n \setminus \tilde{\Omega}$ is uniformly $p$-thick with constant $\tilde{b}_0,\tilde{r}_0$. Choose $y \in \tilde{\Omega}$ such that $B_{\frac34 r}(y)\setminus \tilde{\Omega} \neq \emptyset$, then there exists a constant $\tilde{b}_1 = \tilde{b}_1(b_0,r_0,n,p)>0$ such that 
 \[
  \cpt(\overline{B_{2r}(y)}\setminus \tilde{\Omega}, B_{4r}(y) ) \geq \tilde{b}_1 \cpt( \overline{B_{2r}(y)}, B_{4r}(y) ). 
 \]
\end{lemma}

Following the definition of $p$-thickness, a simple consequence of H\"older's and Young's inequality gives the following result (for example, see \cite[Lemma 3.13]{Par} for the proof):
\betagin{lemma}
 Let $1 < p \leq n$ be given and suppose  a set $E \subset \mathbb{R}^n$ is uniformly $p$-thick with constants $\tilde{b}_0,\tilde{r}_0$. Then $E$ is uniformly $q$-thick for all $q \geq p$ with constants $\tilde{b}_1, \tilde{r}_1$.
\end{lemma}

We remark that the 

notion of $p$-thickness is self improving (see \cite{Le88} or  \cite{Anc,Mik} for the details):
\betagin{theorem}[\cite{Le88}]
\lambdabel{self_improv_cap}
 Let $1 < p \leq n$ be given and suppose a set $E \subset \mathbb{R}^n$ is uniformly $p$-thick with constants $\tilde{b}_0,\tilde{r}_0$. Then there exists an exponent $q = q(n,p,b_0)$ with $1<q<p$ for which $E$ is uniformly $q$-thick with constants $\tilde{b}_1, \tilde{r}_1$. 
\end{theorem}
We next state a generalized Sobolev-Poincar\'e's inequality which was originally obtained by V. Maz'ya \cite[Sec. 10.1.2]{Maz} (see also    \cite[Sec. 3.1]{KK94} and \cite[Corollary 8.2.7]{AH}). 
\betagin{theorem} \lambdabel{sobolev-poincare} Let $B$ be a ball and $\phii \in W^{1,p}(B)$ be $p$-quasicontinuous, with $p>1$. Let $\kappappa=n/(n-p)$ if $1<p<n$ and $\kappappa=2$ if 
$p=n$. Then there exists a constant $c = c(n, p) > 0$ such that 
\betagin{equation}\lambdabel{sob_poin_est}
 \left( \Xint\tilthlongdash_B |\phii|^{\kappappa p} \, dx\right)^{\frac{1}{\kappappa p}} \leq c \left( 
\frac{1}{{\rm cap}_{1,p} (\overline{N(\phii)}, 2B ) } \int_B |\nabla \phii|^p \, dx\right)^{\frac{1}{p}},
\end{equation}
where $N(\phii)=\{x \in B: \phii (x) = 0 \}$.
\end{theorem}

\betagin{lemma}
\lambdabel{interpolation_capacity_version_main}
Let $B_{\rhoo}$ be a ball and let   $1 < \sigma < \infty$ ,  $1< \vartheta \leq n$  be any two exponents. Let $g $ be any function such that the following assumption holds: 
$$\cpt[\vartheta] (\overline{N(g)},2B) \geq c_0 \cpt[\vartheta] (\overline{B},2B) = c_0 \rhoo^{n-p}, \quad \text{where} \ N(g) = \{ x \in B : g(x) = 0\}. $$
Let $\tilde{\alphapha} \in (0,1)$ and $1<\tilde{q},\tilde{r}<\infty$  be given such that the exponents are all related by the relations:
\betagin{gather}
 \tilde{\alphapha} \sigma < \tilde{q} \leq \vartheta^* = \frac{n \vartheta}{n-\vartheta}   \qquad \text{and} \qquad \tilde{r}:={\frac{(1-\tilde{\alphapha}) \sigma \tilde{q}}{ \tilde{q} - \tilde{\alphapha} \sigma}} \Longleftrightarrow \tilde{q}:= \frac{\tilde{\alphapha} \sigma \tilde{r}}{\tilde{r} - (1-\tilde{\alphapha})\sigma}. \lambdabel{lemma1.7.1}

\end{gather}
which is equivalent to the relation
\betagin{gather*}
     - \frac{n}{\sigma} \leq \tilde{\alphapha}\lbr 1 - \frac{n}{\vartheta} \rbr  - (1-\tilde{\alphapha}) \frac{n}{\tilde{r}}. \nonumber 
\end{gather*}
Then the following interpolation type estimate holds:
\[
 \Xint\tilthlongdash_B \lbr \frac{|g|}{\rhoo}\rbr^{\sigma} \ dx  \leq C_{(n,\sigma,\vartheta,\tilde{r},c_0)} \lbr \Xint\tilthlongdash_B|\nabla g|^{ \vartheta}  \ dx \rbr^{\frac{\tilde{\alphapha} \sigma}{\vartheta}} \lbr \Xint\tilthlongdash_B \lbr \frac{|g|}{\rhoo}\rbr^{\tilde{r}} \ dx \rbr^{\frac{(1-\tilde{\alphapha})\sigma}{\tilde{r}}}.
\]
\end{lemma}
\betagin{proof}
Let us first consider the case $\vartheta<n$. The proof follows by applying Theorem \ref{sobolev-poincare} as follows: let $\tilde{\alphapha} \sigma < \tilde{q} \leq \vartheta^* = \frac{n \vartheta}{n-\vartheta}$, then making use of \eqref{lemma1.7.1},  we obtain
\betagin{equation*}
 \betagin{array}{ll}
  \Xint\tilthlongdash_B \lbr \frac{|g|}{\rhoo}\rbr^{\sigma} \ dx 

  & \apprle \lbr \Xint\tilthlongdash_B \lbr \frac{|g|}{\rhoo} \rbr^{ \vartheta^*} \ dx \rbr^{\frac{\tilde{\alphapha} \sigma}{\vartheta^*}} \lbr \Xint\tilthlongdash_B \lbr \frac{|g|}{\rhoo}\rbr^{\tilde{r}} \ dx \rbr^{\frac{(1-\tilde{\alphapha})\sigma}{\tilde{r}}} \\
  & \apprle \lbr \Xint\tilthlongdash_B |\nabla g|^{ \vartheta} \ dx \rbr^{\frac{\tilde{\alphapha} \sigma}{\vartheta}} \lbr \Xint\tilthlongdash_B \lbr \frac{|g|}{\rhoo}\rbr^{\tilde{r}} \ dx \rbr^{\frac{(1-\tilde{\alphapha})\sigma}{\tilde{r}}}.

 \end{array}
\end{equation*}

 In the case $\vartheta =n$ (and also $\vartheta >n$), we observe that $W^{1,\vartheta}$ embeds into $L^t$ for all $t< \infty$ and the proof becomes simpler.
\end{proof}

\betagin{definition}
 \lambdabel{p_thick_domain}
 In this paper, we shall assume that the domain $\Omega$ is bounded and that it's complement $\Omega^c$ is uniformly $p$-thick with constants $b_0,r_0$ in the sense of Definition of \ref{p_thick}.
\end{definition}

\subsection{Assumptions on the Nonlinear structure}

We shall now collect the assumptions on the nonlinear structure in \eqref{main}. 
We assume that $\mathcal{A}(x,t,\nabla u)$ is an Carath\'eodory function, i.e., we have $(x,t) \mapsto \mathcal{A}(x,t,\zetata)$  is measurable for every $\zetata \in \mathbb{R}^n$ and 
$\zetata \mapsto \mathcal{A}(x,t,\zetata)$ is continuous for almost every  $(x,t) \in \Omega \times (-T,T)$.

We further assume that for a.e. $(x,t) \in \Omega \times (-T,T)$ and for any $\zetata \in \mathbb{R}^n$, there exists two given positive constants $\Lambda_0,\Lambda_1$ such that  the following bounds are satisfied   by the nonlinear structures :
\betagin{gather}
 \iprod{\mathcal{A}(x,t,\zetata)}{\zetata} \geq \Lambda_0 |\zetata |^p  - h_1  \txt{and} |\mathcal{A}(x,t,\zetata)| \leq \Lambda_1  |\zetata|^{p-1} + h_2, \lambdabel{abounded}

\end{gather}
where, the functions $h_1, h_2: \Omega \times (-T,T) \mapsto \mathbb{R}$ are assumed to be  measurable  with bounded norm 
\betagin{equation}
 \lambdabel{bound_b}
 \theta^p:=|h_1| + |h_2|^{\frac{p}{p-1}} \txt{and} \| \theta\|_{L^{\hat{q}}(\Omega \times (-T,T))} < \infty \quad \text{ for some  } \hat{q} \geq p+\beta_0,
\end{equation}
where $\beta_0$ is from Theorem \ref{main_theorem} (for example, $\beta_0 = 1$ is admissible in  \eqref{bound_b}).  \emph{An important aspect of the estimates obtained in this paper is that we do not make any  assumptions regarding the smoothness of $\mathcal{A}(x,t,\zetata)$ with respect to $x,t,\zetata$}.

As the basic sets for our estimates, we will use parabolic cylinders where the radii in space and time are coupled. This is due to the fact that in the case that $p \neq 2$, the size of the cylinders implicitly depend on the solution itself and cannot be taken independently. This difficulty extends to the problems dealing with \emph{very weak} solutions also. 

In what follows, we will always assume 
\betagin{equation}
 \lambdabel{restriction_p}
 \frac{2n}{n+2} < p < \infty.
\end{equation}

\betagin{remark}The restriction in \eqref{restriction_p} is necessary when dealing with parabolic problems because of the embedding $W^{1,\frac{2n}{n+2}} \hookrightarrow L^2$ due to which,  we will invariably have to deal with the $L^2$-norm of the solution which comes from the time-derivative. On the other hand, if we assume $u \in L^{r}(\Omega \times (-T,T))$ for some $r\geq 1$ such that $n(p-2) + rp >0$ (see \cite[Chapter 5]{DiB1} for more on this), then we can obtain analogous result as to Theorem \ref{main_theorem}. This extension of Theorem \ref{main_theorem} to the case $1< p \leq \frac{2n}{n+2}$ is beyond the scope of this paper and will be pursued elsewhere. \end{remark}

\subsection{Function Spaces}\lambdabel{function_spaces}

Let $1\leq \vartheta < \infty$, then $W_0^{1,\vartheta}(\Omega)$ denotes the standard Sobolev space which is the completion of $C_c^{\infty}(\Omega)$ under the $\|\cdot\|_{W^{1,\vartheta}}$ norm. 

The parabolic space $L^{\vartheta}(-T,T; W^{1,\vartheta}(\Omega))$ for any $\vartheta \in (1,\infty)$ is the collection of measurable functions $f(x,t)$ such that for almost every $t \in (-T,T)$, the function $x \mapsto f(x,t)$ belongs to $W^{1,\vartheta}(\Omega)$ with the following norm being finite:
\[
 \| f\|_{L^{\vartheta}(-T,T; W^{1,\vartheta}(\Omega)} := \lbr \int_{-T}^T \| u(\cdot, t) \|_{W^{1,\vartheta}(\Omega)}^{\vartheta} \ dt \rbr^{\frac{1}{\vartheta}} < \infty.
\]

Analogously, the parabolic space $L^{\vartheta}(-T,T; W_0^{1,\vartheta}(\Omega))$ is the collection of measurable functions $f(x,t)$ such that for almost every $t \in (-T,T)$, the function $x \mapsto f(x,t)$ belongs to $W_0^{1,\vartheta}(\Omega)$.

\subsection{Notations}

We shall clarify the notation that will be used throughout the paper: 
\betagin{enumerate}[(i)]
    
 \item\lambdabel{not1} We shall use $\nabla$ to denote derivatives with respect the space variable $x$.
\item\lambdabel{not2} We shall sometimes alternate between using $\ddt{f}$, $\partial_t f$ and $f'$ to denote the time derivative of a function $f$.
 \item\lambdabel{not3} We shall use $D$ to denote the derivative with respect to both the space variable $x$ and time variable $t$ in $\mathbb{R}^{n+1}$. 
 \item\lambdabel{not4}  Let $z_0 = (x_0,t_0) \in \mathbb{R}^{n+1}$ be a point and $\rhoo, s >0$ be two given parameters and let $\alpha \in (0,\infty)$. We shall use the following symbols to denote the following regions:
 \betagin{equation}\lambdabel{notation_space_time}
\deltaf2.2{1.5}
 \betagin{array}{ll}
 I_s(t_0) := (t_0 - s, t_0+s) \subset \mathbb{R},
& \qquad Q_{\rhoo,s}(z_0) := B_{\rhoo}(x_0) \times I_{s}(t_0) \subset \mathbb{R}^{n+1},\\ 
 \alpha Q_{\rhoo,s}(z_0) := B_{\alpha \rhoo}(x_0) \times I_{\alpha^2s}(t_0)  \subset \mathbb{R}^{n+1},
 &\qquad \mathcal{H}_s(t_0) := \mathbb{R}^n \times I_s(t_0) \subset \mathbb{R}^{n+1},\\ 
 \alphapha \mathcal{H}_s(t_0) := \mathbb{R}^n \times I_{\alphapha^2 s}(t_0) \subset \mathbb{R}^{n+1},
 & \qquad\mathcal{C}_{\rhoo}(x_0) := \Omega \cap B_{\rhoo}(x_0) \times \mathbb{R} \subset \mathbb{R}^{n+1},\\

\Omega_{\rhoo,s}(z_0) := \Omega \cap B_{\rhoo}(x_0) \times I_s(t_0) \subset \mathbb{R}^{n+1}, 
& \qquad\Omega_{\rhoo}(x_0) := \Omega \cap B_{\rhoo}(x_0) \subset \mathbb{R}^{n}.
 \end{array}
\end{equation}

\item\lambdabel{not5} We shall use $\int$ to denote the integral with respect to either space variable or time variable and use $\iint$ to denote the integral with respect to both space and time variables simultaneously. 

Analogously, we will use $\fint$ and $\Yint\tiltlongdash$ to denote the average integrals as defined below: for any set $A \times B \subset \mathbb{R}^n \times \mathbb{R}$, we define
\betagin{gather*}

\avgs{A}{f}:= \fint_A f(x) \ dx = \frac{1}{|A|} \int_A f(x) \ dx,\\
\avgs{A\times B}{f}:=\Yint\tiltlongdash_{A\times B} f(x,t) \ dx \ dt = \frac{1}{|A\times B|} \iint_{A\times B} f(x,t) \ dx \ dt.

\end{gather*}

\item\lambdabel{not6} Given any positive function $\mu$, we shall denote $\avgs{\mu}{f} := \int f\frac{\mu}{\|\mu\|_{L^1}}dm$ where the domain of integration is the domain of definition of $\mu$ and $dm$ denotes the associated measure. 

\item\lambdabel{not7} Given any $\lambda> 0$, we can convert $\mathbb{R}^{n+1}$ into a metric space where the parabolic cylinders correspond to \emph{balls} under the parabolic metric given by:
\betagin{equation}
 \lambdabel{parabolic_metric}
d_{\lambda} (z_1,z_2) := \max \left\{ |x_2-x_1|, \sqrt{\lambda^{p-2}|t_2-t_1|} \right\}.
\end{equation}

\item\lambdabel{not8} In what follows, $r_0$ and $b_0$ will denote the constants arising from the assumption that  $\Omega^c$ is uniformly $p$-thick (see Definition \ref{p_thick_domain}).
\end{enumerate}

\subsection{Very weak solution}

There is a well known difficulty in defining the notion of solution for \eqref{main} due to a lack of time derivative of $u$. To overcome this, one can either use Steklov average or convolution in time. In this paper, we shall use the former approach (see also \cite[Page 20, Eqn (2.5)]{DiB1} for further details).

Let us first define Steklov average as follows: let $h \in (0, 2T)$ be any positive number, then we define
\betagin{equation}\lambdabel{stek1}

  u_{h}(\cdot,t) := \left\{ \betagin{array}{ll}
                              \Xint\tilthlongdash_t^{t+h} u(\cdot, \tauu) \ d\tauu \quad & t\in (-T,T-h), \\
                              0 & \text{else}.
                             \end{array}\right.
 \end{equation}

\betagin{definition}[Very weak solution] 
\lambdabel{very_weak_solution}

 

 Let $ \beta \in (0,1)$ and $h \in (0,2T)$ be given and suppose $p-\beta > 1$. We  then say $u \in L^2(-T,T; L^2(\Omega)) \cap L^{p-\beta}(-T,T; W_0^{1,p-\beta}(\Omega))$ is a very weak solution of \eqref{main} if for any $\phii \in W_0^{1,\frac{p-\beta}{1-\beta}}(\Omega) \cap L^{\infty}(\Omega)$, the following holds:
 \betagin{equation}
 \lambdabel{def_weak_solution}
  \int_{\Omega \times \{t\}} \frac{d [u]_{h}}{dt} \phii + \iprod{[\mathcal{A}(x,t,\nabla u)]_{h}}{\nabla \phii} \ dx = 0 \txt{for any}-T < t < T-h.
 \end{equation}

 
\end{definition}

 \subsection{Maximal Function}
 
 
For any $f \in L^1(\mathbb{R}^{n+1})$, let us now define the strong maximal function in $\mathbb{R}^{n+1}$ as follows:
\betagin{equation*}
 \lambdabel{par_max}
 \mathcal{M}(|f|)(x,t) := \sup_{\tilde{Q} \ni(x,t)} \Yint\tiltlongdash_{\tilde{Q}} |f(y,s)| \ dy \ ds
\end{equation*}
where the supremum is  taken over all parabolic cylinders $\tilde{Q}_{a,b}$ with $a,b \in \mathbb{R}^+$ such that $(x,t)\in \tilde{Q}_{a,b}$. An application of the Hardy-Littlewood maximal theorem in $x-$ and $t-$ directions shows that the Hardy-Littlewood maximal theorem still holds for this type of maximal function (see \cite[Lemma 7.9]{Gary} for details):
\betagin{lemma}
\lambdabel{max_bound}
 If $f \in L^1(\mathbb{R}^{n+1})$, then for any $\alpha >0 $, there holds
 \[
  |\{ z \in \mathbb{R}^{n+1} : \mathcal{M}(|f|)(z) > \alpha\}| \leq \frac{5^{n+2}}{\alpha} \|f\|_{L^1(\mathbb{R}^{n+1})}.
 \]

 and if $f \in L^{\vartheta}(\mathbb{R}^{n+1})$ for some $1 < \vartheta \leq \infty$, then there holds
 \[
  \| \mathcal{M}(|f|) \|_{L^{\vartheta}(\mathbb{R}^{n+1})} \leq C_{(n,\vartheta)} \| f \|_{L^{\vartheta}(\mathbb{R}^{n+1})}.
 \]

\end{lemma}

 \subsection{A few well known lemmas}\lambdabel{few_well_known_lemmas}

 We shall recall the following well known lemmas that will be used throughout the paper. 
 The first one is a standard lemma regarding integral averages (for a proof in this setting, see for example \cite[Chapter 8.2]{bogelein2007thesis} for the details).
\betagin{lemma}
\lambdabel{time_average}
Let $\gamma >0$ be any fixed number and suppose  $[f]_h(x,t) : = \Xint\tilthlongdash_{t-\gamma h^2}^{t+\gamma h^2} f(x,\tauu) \ d\tauu$. Then we have the following properties:
\betagin{enumerate}[(i)]
 \item $[f]_h \rightarrow f$ a.e $(x,t) \in \mathbb{R}^{n+1}$ as $h \searrow 0$, $[f]_h$ is continuous and  bounded  in time for a.e. $x \in \mathbb{R}^n$.
 \item For any cylinder $Q_{r, \gamma r^2} \subset \mathbb{R}^{n+1}$ with $r >0$, there holds
\[
  \Yint\tiltlongdash_{Q_{r,\gamma r^2}} [f]_h(x,t) \ dx \ dt \apprle_n \Yint\tiltlongdash_{Q_{r,\gamma(r+h)^2}}f(x,t) \ dx \ dt. 
 \]
 \item The function $[f]_h(x,t)$ is differentiable with respect to $t \in \mathbb{R}$, moreover $[f]_h(\cdot, t) \in C^{1}(\mathbb{R})$ for a.e. $x \in \mathbb{R}^n$. 
 \end{enumerate}
\end{lemma}

The next Lemma we recall is a standard iteration lemma (see for example  \cite[Lemma 5.1]{Gia}):
\betagin{lemma}
\lambdabel{iteration_lemma}	
 Let $0< \thetaeta < 1$, $B \geq 0$, $A \geq 0$, $\tilde{\alphapha} >0$ and $0 < r < \rhoo < \infty$ and let $f \geq 0$ be a bounded measurable function satisfying 
 \[
  f(r_1) \leq \thetaeta f(r_2) + A(r_2-r_1)^{-\tilde{\alphapha}} + B \qquad \forall \ r < r_1 < r_2 < \rhoo,
 \]
then there exists a constant $C=C(\tilde{\alphapha},A,B)$ such that the following holds:
\[
 f(r) \leq C (A (\rhoo - r)^{-\tilde{\alphapha}} + B), \qquad \text{ for every }\ r < \rhoo.
\]

\end{lemma}

 We could not find the proof of the following version of the  Parabolic Poincar\'e Lemma in literature, hence for the sake of completeness, we also present the proof.
\betagin{lemma}
 \lambdabel{lemma_crucial_1}
 Let $f \in L^{\vartheta} (-T,T; W^{1,\vartheta}(\Omega))$ with $\vartheta \in (1,\infty)$  and suppose that $\mathcal{B}_{r} \Subset \Omega$ be compactly contained ball of radius $r>0$. Let $I \subset (-T,T)$ be a time interval  and $\rhoo(x,t) \in L^1(\mathcal{B}_r \times I)$ be any positive function such that $\|\rhoo\|_{L^{\infty}(\mathcal{B}_r\times I)} \apprle_n \frac{\|\rhoo\|_{L^1(\mathcal{B}_r\times I)}}{|\mathcal{B}_r\times I|} $ and $\mu(x) \in C_c^{\infty}(\mathcal{B}_r)$ such that $\int_{\mathcal{B}_r} \mu(x) \ dx = 1$ with $|\mu| \leq \frac{C{(n)}}{r^n}$ and $|\nabla \mu| \leq  \frac{C(n)}{r^{n+1}}$. Then there holds
 \betagin{equation*}
 \betagin{array}{ll}
  \Yint\tiltlongdash_{\mathcal{B}_r \times I} \left|\frac{f - \avgs{\rhoo}{f}}{r}\right|^{\vartheta} \ dz & \apprle_{n,\vartheta} \Yint\tiltlongdash_{\mathcal{B}_r \times I} |\nabla f|^{\vartheta} \ dz + \sup_{t_1,t_2 \in I} \left| \frac{\avgs{\mu}{f}(t_2) - \avgs{\mu}{f}(t_1)}{r} \right|^{\vartheta},
  \end{array}
 \end{equation*}
where $\avgs{\rhoo}{f}:= \int_{\mathcal{B}_r\times I} f(z) \frac{\rhoo(z)}{\|\rhoo\|_{L^1(\mathcal{B}_r\times I)}} \ dz $ and $\avgs{\mu}{f}(t_i) := \int_{\mathcal{B}_r} f(x,t_i) \mu(x) \ dx$ for $i = 1,2$. 
\end{lemma}
\betagin{proof}
 Let us first consider the case of $\rhoo(x,t) = \mu(x)\lsb{\chii}{I}(t)$. In this case, we get
 \betagin{equation}
 \lambdabel{lem5.1.2}
  \betagin{array}{rcl}
  \Yint\tiltlongdash_{\mathcal{B}_r\times I} \left| \frac{f(z)  - \avgs{\mu \times \lsb{\chii}{I}}{f}}{r} \right|^{\vartheta} \ dz 

 & \overset{\redlabel{lemma2.14.a}{a}}{\apprle} &  \Yint\tiltlongdash_{\mathcal{B}_r\times I} |\nabla f|^{\vartheta} \ dz + \sup_{t_1,t_2 \in {I}} \left|\frac{\avgs{\mu}{f}(t_2) - \avgs{\mu}{f}(t_1)}{r}\right|^{\vartheta}.
  \end{array}
 \end{equation}
Here to obtain $\redref{lemma2.14.a}{a}$ above, we made of the standard Poincar\'e's inequality.

For the general case, we observe that 
\betagin{equation}
  \lambdabel{lem5.1.3}
  \Yint\tiltlongdash_{\mathcal{B}_r\times I} \left|\frac{f - \avgs{\rhoo}{f}}{r}\right|^{\vartheta} \ dz \apprle \Yint\tiltlongdash_{\mathcal{B}_r\times I} \left|\frac{f - \avgs{\mu \times \lsb{\chii}{I}}{f}}{r}\right|^{\vartheta} \ dz  + \Yint\tiltlongdash_{\mathcal{B}_r\times I} \left|\frac{\avgs{\rhoo}{f} - \avgs{\mu \times \lsb{\chii}{I}}{f}}{r}\right|^{\vartheta} \ dz.
 \end{equation}
The first term of \eqref{lem5.1.3} can be controlled as in \eqref{lem5.1.2} and  to control the second term, we observe that
\betagin{equation*}
 \lambdabel{lem5.1.4}
 \betagin{array}{ll}
 |\avgs{\rhoo}{f} - \avgs{\mu \times \lsb{\chii}{I}}{f}| & \leq \frac{\|\rhoo\|_{L^\infty(\mathcal{B}_r\times I)}}{\|\rhoo\|_{L^1(\mathcal{B}_r\times I)}} \iint_{\mathcal{B}_r\times I} |f - \avgs{\mu \times \lsb{\chii}{I}}{f}|\ dz  \apprle \Yint\tiltlongdash_{\mathcal{B}_r\times I} |f - \avgs{\mu \times \lsb{\chii}{I}}{f}|\ dz.\\

 \end{array}
\end{equation*}
This completes the proof of the Lemma. 
\end{proof}

We will use the following result which can be found in \cite[Theorem 3.1]{PG} (see also \cite{Prato}):
\betagin{lemma}
\lambdabel{metric_lipschitz}
 Let $\lambda >0$ and  $\mathcal{Z} \subset \mathbb{R}^{n+1}$ be given. For any $z \in \mathcal{Z}$ and $r>0$, let $Q_{r,\gamma r^2}(z)$ be the parabolic cylinder centred  at $z$ with radius $r$.  Suppose  there exists a constant $C>0$ independent of $z$ and $r$ such that the following bound holds:
 \[
  \frac{1}{|Q_{r,\gamma r^2}(z)\cap \mathcal{Z}|} \iint_{Q_{r,\gamma r^2}(z)\cap \mathcal{Z}} \left|\frac{f(x,t) - \avgs{Q_{r,\gamma r^2}(z)\cap \mathcal{Z}}{f}}{r}\right| \ dx \ dt \leq C \qquad \forall\  z \in \mathcal{Z} \text{ and } r>0.
 \]
 Then $f \in C^{0,1} (\mathcal{Z})$. 
\end{lemma}

Let us also recall a parabolic version of Gehring's lemma (for example, see \cite[Lemma 6.4]{Bog}).
\textcolor{black}{
\betagin{lemma}
\lambdabel{gehring_lemma}
 Let $\alpha_0 \geq 1$, $\kappa \geq 1$, $\varepsilon_0>0$, $p>1$ and $\beta_0>0$ be given. Let $q$ be given such that  $1<p-\varepsilon_0 \leq q < p-2\beta<p-\beta$ for some $\beta \in (0,\beta_0)$.  Furthermore, for a cylinder $Q_2 = Q_{2,2^2}$, let  $f \in L^{p-\beta}(Q_2)$ and $g \in L^{\tilde{p}}(Q_2)$ for some $\tilde{p} \geq p$ be given. Suppose for each $\lambda \geq \alpha_0$ and almost every $z \in Q_2$ with $f(z) >\lambda$, there exists a parabolic cylinder $Q = Q_{\rhoo,s}(z)\subset Q_2$ such that 
 \[
  \frac{\lambda^{p-\beta}}{\kappa} \leq \Yint\tiltlongdash_Q f^{p-\beta}(x,t) \ dx \ dt \leq \kappa \lbr \Yint\tiltlongdash_Q f^q(x,t) \ dx\ dt \rbr^{\frac{p-\beta}{q}} + \kappa \Yint\tiltlongdash_Q g^{p-\beta}(x,t) \ dx \ dt \leq \kappa^2 \lambda^{p-\beta},
 \]
then there exists $\delta_0 = \delta_0(\kappa,p,\beta,q,\varepsilon_0)$ and $C = C(\kappa,p,\beta,q,\varepsilon_0)$, such that $f \in L^{p-\beta + \delta_1}(Q_2)$ with $\delta_1 = \min \{ \delta_0, \tilde{p}-p+\beta\}$. This improved higher integrability comes with the following bound:
\[
 \iint_{Q_2} f^{p-\beta + \delta}(x,t) \ dx \ dt\apprle \alpha_0^{\delta} \iint_{Q_2} f^{p-\beta}(x,t) \ dx\ dt + \iint_{Q_2} g^{p-\beta+\delta}(x,t) \ dx\ dt \quad \text{for all} \ \delta \in (0,\delta_1]. 
\]
\end{lemma}}

\section{Construction of test function.}\lambdabel{Construction_of_test_function}

Since we are interested in estimates near the boundary, for the subsequent calculations, we will fix a parabolic cylinder of the form $Q:=Q_{\rhoo, s} (0,0) = B \times I$, where $(0,0) \in \partial \Omega \times (-T,T)$. We will further assume 
\betagin{equation}
 \lambdabel{def_s}
 s = \alpha_0^{2-p} \rhoo^2 \txt{for some} \alpha_0 >0 \ \text{to be chosen}.
\end{equation}
\emph{Since $(0,0)$ and $(\rhoo,s)$ are fixed, we shall henceforth simplify the notation in  \eqref{notation_space_time} by dropping obvious references to  $\rhoo$ and $s$.}

\betagin{definition}
 \lambdabel{exponent_choice} 
 Since we assume that domain $\Omega^c$ is uniformly $p$-thick with constants $r_0,b_0$ (see Definition \ref{p_thick_domain}), by the self improving property (see Theorem \ref{self_improv_cap}), there exists an $\varepsilon_0 > 0$ such that $\Omega^c$ is uniformly  $p-\varepsilon_0$ thick with constants $\tilde{b}_0,\tilde{r}_0$.  Let us now take $\beta$ and $q$ such that the following holds:
 \betagin{equation}
  \lambdabel{def_q}
  p-\varepsilon_0 < q \leq p-2\beta < p-\beta < p. 
 \end{equation}
Note that we will eventually obtain $\beta_0 = \beta_0(n,p,b_0,\Lambda_0,\Lambda_1)\in (0,\varepsilon_0/2)$  such that all the calculations will hold for any $\beta \in (0,\beta_0]$. 
\end{definition}

Consider the following cut off functions:
\betagin{equation}
\lambdabel{cut_off_function}\betagin{array}{c}
 \eta(x) \in C_c^{\infty} (8B) , \ \eta(x) \equiv 1 \ \text{for} \ x \in 4B,  \qquad  \zetata(t) \in C_c^{\infty} (8I) , \ \zetata(t) \equiv 1 \ \text{for} \ t \in 4I , \\
 |\eta(x)| \leq 1, \qquad |\eta(x)| \leq \frac{C(n)}{\rhoo}, \txt{and} |\zetata(t)| \leq 1, \qquad | \zetata'(t)| \leq \frac{C}{s}.
 \end{array}
\end{equation}

\betagin{definition}
 Since we assume that there exists a very weak solution $u \in L^2(-T,T; L^2(\Omega)) \cap L^{p-\beta} (-T,T; W_0^{1,p-\beta}(\Omega))$, we can now multiply the solution by the cut off functions to get:
 \betagin{equation}
  \lambdabel{def_u_tilde}
  \tilde{u}h(x,t) := [u]_h(x,t) \eta(x) \zetata(t).
 \end{equation}
Here $[u]_h$ denotes the usual Steklov average defined in \eqref{stek1}. 
\end{definition}

Let us define the following function with $\theta$ as in \eqref{bound_b}:
\betagin{equation}
 \lambdabel{def_g}
 g(x,t) := \max \left\{ \mathcal{M}(|\nabla \tilde{u}|^q)^{\frac{1}{q}}(x,t), \mathcal{M}((|\nabla u|+|\theta|)^q\lsb{\chii}{8Q})^{\frac{1}{q}}(x,t)\right\} \quad \text{for any} \ (x,t) \in \mathbb{R}^{n+1}.
\end{equation}
We have the following bound for $g(x,t)$:
\betagin{lemma}
 \lambdabel{bound_g_x_t}
 Under the assumption that $\Omega$ satisfies Definition \ref{p_thick_domain}, the following bound holds:
 \betagin{equation*}
  \lambdabel{bound_g_x_t_1}
  \| g\|_{L^{p-\beta} (\mathbb{R}^{n+1})} \leq C_{(n,p,b_0)} \| (|\nabla u| + |\theta|) \lsb{\chii}{8Q} \|_{L^{p-\beta}(\mathbb{R}^{n+1})}.
 \end{equation*}
\end{lemma}
\betagin{proof}
 Making use of  Lemma \ref{max_bound}, we get
 \betagin{equation*}
  \betagin{array}{ll}
   \iint_{\mathbb{R}^{n+1}} |g(x,t)|^{p-\beta} \ dz & \apprle \iint_{\mathbb{R}^{n+1} }\lbr \mathcal{M}(|\nabla \tilde{u}|^q)^{\frac{p-\beta}{q}}(z)+ \mathcal{M}((|\nabla u|+|\theta|)^q\lsb{\chii}{8Q})^{\frac{p-\beta}{q}}(z)\rbr \ dz \\
   & \apprle \iint_{\mathbb{R}^{n+1} }\lbr |\nabla \tilde{u}|^{p-\beta}+ \lbr |\nabla u|+|\theta|\rbr^{p-\beta} \lsb{\chii}{8Q} \rbr \ dz.
  \end{array}
 \end{equation*}
 Since $\tilde{u}(x,t) = u(x,t) \eta(x) \zetata(t)$ and $\spt(\tilde{u}) \subset 8Q$, we use Theorem \ref{sobolev-poincare} to get
 \[

   \iint_{\mathbb{R}^{n+1} } |\nabla \tilde{u}|^{p-\beta}\ dz  \apprle \iint_{8Q}\lbr |\nabla u|^{p-\beta}  + \left|\frac{u}{\rhoo}\right|^{p-\beta} \rbr \ dz 
    \apprle \iint_{8Q} |\nabla u|^{p-\beta}\ dz.

 \]
This completes the proof of the Lemma. 
\end{proof}

We now define the good set
\betagin{equation*}
 \lambdabel{def_elam}
\elam:=  \left\{ (x,t) \in \mathbb{R}^{n+1} : g(x,t) \leq \lambda \right\} \quad \text{for} \ \lambda > c_e\alpha_0. 
 \end{equation*}
 \emph{Note that $\alpha_0>0$  and $c_e$ are constants that will be quantified explicitly in Section \ref{Caccioppoli_inequality}. }

We now have the following parabolic Whitney type decomposition of $\elam^c$ which is taken from  \cite[Section 2.3]{Seb} {\color{black}(see also \cite[Chapter 3]{bogelein2013regularity} where this was proved in the stated form or \cite[Lemma 3.1]{diening2010existence} for details):}

\betagin{lemma}
\lambdabel{whitney_decomposition}

Let $\gamma := \lambda^{2-p}$, then there exists a $\gamma$-parabolic Whitney covering  $\{Q_i(z_i)\}$ of $\elam^c$ in the following sense:
 \betagin{description}
 \setcounter{whitney}{1}
  \deltascitem{W\thetaewhitney}{W7} $Q_j(z_j) = B_j(x_j) \times I_j(t_j)$ where $B_j(x_j) = B_{r_j}(x_j)$ and $I_j(t_j) = (t_j - \gamma r_j^2, t_j + \gamma r_j^2)$. 
  \stepcounter{whitney}\deltascitem{W\thetaewhitney}{W2} Since $\gamma = \lambda^{2-p}$, recalling the notation from \eqref{parabolic_metric}, we have $d_{\lambda}(z_j,\elam) = 16r_j$.
  \stepcounter{whitney}\deltascitem{W\thetaewhitney}{W1} $\bigcup_j \frac12 Q_j(z_j) = \elam^c$.
  \stepcounter{whitney}\deltascitem{W\thetaewhitney}{W3} for all $j \in \mathbb{N}$, we have $8Q_j \subset \elam^c$ and $16Q_j \cap \elam \neq \emptyset$.
  \stepcounter{whitney}\deltascitem{W\thetaewhitney}{W4} if $Q_j \cap Q_k \neq \emptyset$, then $\frac12 r_k \leq r_j \leq 2r_k$.
  \stepcounter{whitney}\deltascitem{W\thetaewhitney}{W5} $\frac14 Q_j \cap \frac14Q_k = \emptyset$ for all $j \neq k$.
  \stepcounter{whitney}\deltascitem{W\thetaewhitney}{W6} $\sum_j \lsb{\chii}{4Q_j}(z) \leq c(n)$ for all $z \in \elam^c$.
  \end{description}
  Subject to this Whitney covering, we have an associated partition of unity denoted by $\{ \omega_j\}_{j\in\mathbb{N}} \in C_c^{\infty}(\mathbb{R}^{n+1})$ such that the following holds:
  \betagin{description}
  \stepcounter{whitney}\deltascitem{W\thetaewhitney}{W8} $\lsb{\chii}{\frac12Q_j} \leq \omega_j \leq \lsb{\chii}{\frac34Q_j}$.
  \stepcounter{whitney}\deltascitem{W\thetaewhitney}{W9} $\|\omega_j\|_{\infty} + r_j \| \nabla \omega_j\|_{\infty} + r_j^2 \| \nabla^2 \omega_j\|_{\infty} + \lambda r_j^2 \| \partial_t \omega_j\|_{\infty} \leq C(n)$.
  \end{description}
  For a fixed $k \in \mathbb{N}$, let us define 
  \betagin{equation*}\lambdabel{Ak}A_k := \left\{ j \in \mathbb{N}: \frac34Q_k \cap \frac34Q_j \neq \emptyset\right\},\end{equation*} then we have
  \betagin{description}
  \stepcounter{whitney}\deltascitem{W\thetaewhitney}{W10} Let $i \in \mathbb{N}$ be given, then $\sum_{j \in A_i} \omega_j(z) = 1$  for all $z \in \frac34Q_i$.
  \stepcounter{whitney}\deltascitem{W\thetaewhitney}{W11} Let $i \in \mathbb{N}$ be given and let  $j \in A_i$, then $\max \{ |Q_j|, |Q_i|\} \leq C(n) |Q_j \cap Q_i|.$
  \stepcounter{whitney}\deltascitem{W\thetaewhitney}{W12}  Let $i \in \mathbb{N}$ be given and let  $j \in A_i$, then $ \max \{ |Q_j|, |Q_i|\} \leq \left|\frac34Q_j \cap \frac34Q_i\right|.$
  \stepcounter{whitney}\deltascitem{W\thetaewhitney}{W13} For any $i \in \mathbb{N}$, we have $\# A_i \leq c(n)$.
  \stepcounter{whitney}\deltascitem{W\thetaewhitney}{W14} Let $i \in \mathbb{N}$ be given, then for any $j \in A_i$, we have $\frac34Q_j \subset 4Q_i$.
 \end{description}
\end{lemma}

Once we have the Whitney-type parabolic decomposition as in Lemma \ref{whitney_decomposition}, we can now define the following extension:
\betagin{equation}
 \lambdabel{lipschitz_extension}
\vlh(z) := \tilde{u}h(z) - \sum_i \omega_i(z) \lbr \tilde{u}h(z) - \tilde{u}h^i \rbr. 
\end{equation}
where 
\betagin{equation}
\lambdabel{def_tuh_i}
 \tilde{u}h^i := \left\{ \betagin{array}{ll}
                   \frac1{\|\omega_i\|_{L^1(\frac34Q_i)}}\iint_{\frac34Q_i} \tilde{u}h(z) \omega_i (z) \ dz & \text{if} \ \frac34Q_i \subset \Omega_{8\rhoo}\times \mathbb{R} = \mathcal{C}_{8\rhoo}, \\
                   0 & \text{else}.
                  \end{array}\right. 
\end{equation}

For each $\lambda > c_e \alpha_0$, we will show that $\vlh$ constructed in \eqref{lipschitz_extension} is Lipschitz continuous with respect to the parabolic metric $d_{\lambda}$ on $2\mathcal{H} = \mathbb{R}^n \times (-4s,4s)$. In order to do this, we shall make use of Lemma \ref{metric_lipschitz}. In order to apply Lemma \ref{metric_lipschitz} with $Y = 2\mathcal{H}$, we will consider the following two cases, one where the cylinders are contained inside $4\mathcal{H} = \mathbb{R}^n \times (-16s,16s)$ and the other where the cylinders are not contained inside $4\mathcal{H}$. 
In this regard, let us define the following:
\betagin{gather}
 \Thetaeta := \left\{ i \in \mathbb{N} : \frac34Q_i \cap 2\mathcal{H} \neq \emptyset\right\}, \lambdabel{theta} \\
 \Thetaeta_1 := \left\{ i \in \Thetaeta : 8Q_i \subset 4\mathcal{H} \right\}, \lambdabel{theta_1} \\
 \Thetaeta_2 := \left\{ i \in \Thetaeta : 8Q_i \cap (4\mathcal{H})^c \neq \emptyset\right\}  = \Thetaeta \setminus \Thetaeta_1.\lambdabel{theta_2} 
\end{gather}

We will use the following important lemma throughout  the paper:
\betagin{lemma}
\lambdabel{lemma_crucial_2}
 Let $\gamma >0$ be given and suppose that $u \in L^2(-T,T;L^2(\Omega)) \cap L^{p-\beta}(-T,T; W_0^{1,p-\beta}(\Omega))$ with $0 \leq \beta \leq \min\{1,p-1\}$ be a very weak solution of \eqref{main}. Let $\mathcal{B} \subset \Omega$ be a compactly contained region and $(t_1,t_2) \subset (-T, T-\gamma h^2)$ for some $h \in (0,T)$ be a time interval. Let $\phii(x) \in C_c^{\infty}(\mathcal{B})$, $\varphi(t) \in C_c^{\infty}(t_1,t_2)$ be two non-negative functions and $[u]_h$ be the Steklov average as defined in \eqref{stek1}. Then   the following estimate holds:
 \betagin{equation*}
  \lambdabel{lemma_crucial_2_est}
  \betagin{array}{ll}
  |\avgs{\phii}{{[u]_h\varphi}} (t_2) - \avgs{\phii}{{[u]_h\varphi}}(t_1)| & \leq \Lambda_1 \|\nabla \phii\|_{L^{\infty}{(\mathcal{B})}} \|\varphi\|_{L^{\infty}(t_1,t_2)} \iint_{\mathcal{B} \times (t_1,t_2)} {[|\nabla u|^{p-1} + |\theta|^{p-1}]_h} \ dz + \\

   & \qquad +  \|\phii\|_{L^{\infty}{(\mathcal{B})}} \|\varphi'\|_{L^{\infty}(t_1,t_2)} \iint_{\mathcal{B} \times (t_1,t_2)} |[u]_h| \ dz.
  \end{array}
 \end{equation*}
\end{lemma}
\betagin{proof}
 Let us use $\phii(x)\varphi(t)$ as a test function in \eqref{def_weak_solution} to get 
 \[
  \int_{\Omega\times\{t\}} \ddt{[u]_h} (x,t) \phii(x)\varphi(t) \ dx +  \iprod{[\mathcal{A}(x,t,\nabla u)]_h}{\nabla \phii}(x,t)\varphi(t) \ dx =0
 \]
 Using the Fundamental theorem of calculus, we get
 \betagin{equation*}
  \betagin{array}{ll}
   |\avgs{\phii}{{[u]_h\varphi}} (t_2) - \avgs{\phii}{{[u]_h\varphi}}(t_1)|

   &= \left|\iint_{\mathcal{B} \times (t_1,t_2)} \frac{d}{dt} \lbr[(] [u]_h(x,t) \phii(x)\varphi(t) \rbr[)] \ dz \right|\\
   & \leq \left|\iint_{\mathcal{B} \times (t_1,t_2)} \iprod{[\mathcal{A}(x,t,\nabla u)]_h}{\nabla \phii}(x,t) \varphi(t)\ dz \right| + \\

   & \qquad + \left| \iint_{\mathcal{B} \times (t_1,t_2)} [u]_h(x,t) \phii(x) \ddt{\varphi(t)} \ dz \right|  \\
   & \overset{\redlabel{lemma3.4.1}{a}}{\leq} C(\Lambda_1,p)\|\nabla \phii\|_{L^{\infty}{(\mathcal{B})}} \|\varphi\|_{L^{\infty}(t_1,t_2)} \iint_{\mathcal{B} \times (t_1,t_2)} {[|\nabla u|^{p-1} + |\theta|^{p-1}]_h} \ dz \\

   & \qquad + \|\phii\|_{L^{\infty}{(\mathcal{B})}} \|\varphi'\|_{L^{\infty}(t_1,t_2)} \iint_{\mathcal{B} \times (t_1,t_2)} |[u]_h| \ dz.
  \end{array}
 \end{equation*}
 To obtain \redref{lemma3.4.1}{a} above, we made use of \eqref{abounded} and  \eqref{bound_b} which completes the proof.
\end{proof}

\subsection{Properties of the test function}
In this subsection, we shall prove all the properties that the constructed function in \eqref{lipschitz_extension} satisfies. Recall that $\Omega^c$ is uniformly $p$-thick with constants $b_0,r_0$. 
The first lemma gives a point-wise bound for $\vlh$ on $\elam^c$. 

\betagin{lemma}
\lambdabel{lemma3.7}
 For any $z \in \elam^c$, we have
 \betagin{equation*}
  \lambdabel{bound_v_l_h}
  |\vlh(z)| \leq C_{(n,p,q,\Lambda_0,\Lambda_1,b_0)} \rhoo \lambda.
 \end{equation*}

\end{lemma}

\betagin{proof}
From \eqref{lipschitz_extension},  for any  $z \in \elam^c$, we see that $\vlh(z) = \sum_j \omega_j(z) \tilde{u}h^j$ with $\tilde{u}h^j\neq 0$ if and only if $\frac34Q_j \subset \mathcal{C}_{8\rhoo}$, which automatically implies $\vlh(z) = 0$ for all $z \in \mathcal{C}_{8\rhoo}^c$. 

In order to prove the Lemma, making use of \deltascref{W8}, we see that \eqref{bound_v_l_h} follows if the following holds:
\betagin{equation}
\lambdabel{claim_bound}
|\tilde{u}h^j| \leq C_{(n,p,q,\Lambda_0,\Lambda_1,b_0)}\  \rhoo \lambda,  \txt{for all}\ j \in \mathbb{N}.
\end{equation}

We shall now proceed with proving \eqref{claim_bound}.    Since we only have to consider the case  $\frac34Q_j \subset \mathcal{C}_{8\rhoo}$, we have the following two sub-cases:
 \betagin{description}[leftmargin=*]
  \item[Case $r_j \geq \rhoo$:] In this case, we observe that $8B \subset 32B_j$ which gives the following sequence of estimates:
  \betagin{equation}
  \lambdabel{est_1}
   \betagin{array}{ll}
    |\tilde{u}h^j| 

    & \overset{\redlabel{3.13.a}{a}}{\apprle} \rhoo \Xint\tilthlongdash_{32I_j} \lbr \Xint\tilthlongdash_{32B_j} \left| \frac{\tilde{u}h(x,t)}{r_j} \right|^q \ dx\rbr^{\frac1{q}} \ dt 

     \overset{\redlabel{3.13.b}{b}}{\apprle} \rhoo \lbr \Yint\tiltlongdash_{32Q_j}  \left| \nabla \tilde{u}h(x,t) \right|^q \ dx\ dt\rbr^{\frac1{q}}  
     \overset{\redlabel{3.13.c}{c}}{\apprle} \rhoo \lambda. 
   \end{array}
  \end{equation}

  In order to obtain \redref{3.13.a}{a}, we have used the fact that $r_j \leq \frac{32}{3} \rhoo$ since $\frac34Q_j \subset \mathcal{C}_{8\rhoo}$.  To get \redref{3.13.b}{b}, we  used Theorem \ref{sobolev-poincare} since $32B_j \cap \Omega^c \neq \emptyset$. Finally to obtain \redref{3.13.c}{c}, we made use of \deltascref{W4}.
  
  \item[Case $r_{j} \leq \rhoo$:] Let us define the following constant $k_0 := \min\{ \tilde{k}_1,\tilde{k}_2\}$, where $\tilde{k}_1$ and $\tilde{k}_2$ satisfy
   \betagin{equation}\lambdabel{def_k_0}
2^{\tilde{k}_1 - 1} r_j < \rhoo \leq 2^{\tilde{k}_1} r_j \txt{and} 2^{\tilde{k}_2 -1} Q_j \cap (\Omega \times \mathbb{R})^c = \emptyset \ \text{  but  } \ 2^{\tilde{k}_2}Q_j \cap (\Omega \times \mathbb{R})^c \neq \emptyset.    
   \end{equation}

Using \eqref{cut_off_function} and  triangle inequality, we get
  \betagin{equation}
   \lambdabel{est_2} 
   \betagin{array}{ll}
    |\tilde{u}h^j| 

    & \overset{\text{}}{\apprle} \sum_{m=0}^{k_0 -2} \lbr \avgs{2^m Q_j}{{[u]_h}} - \avgs{2^{m+1}Q_j}{{[u]_h}}\rbr + \avgs{2^{k_0-2}Q_j}{{[u]_h}} := \sum_{m=0}^{k_0-2}S_1^m + S_2.
   \end{array}
  \end{equation}

  \betagin{description}[leftmargin=*]
   \item[Estimate for $S_1^m$ for $m \leq {k_0-2}$:]  In this case,  we see that  $2^{m+1}Q_j \subset \mathcal{C}_{8\rhoo}$. Thus applying Lemma \ref{lemma_crucial_1} for any $\mu \in C_c^{\infty}(B(2^{m+1}r_j, x_j))$ satisfying $| \mu(x)| \leq \frac{C(n)}{(2^{m+1}r_j)^{n}}$ and $|\nabla \mu(x)| \leq \frac{C(n)}{(2^{m+1}r_j)^{n+1}}$, we get
   \betagin{equation}
   \lambdabel{S_1_1}
    \betagin{array}{ll}

     S_1^m

     & \apprle  (2^{m+1}r_j) \lbr \Yint\tiltlongdash_{2^{m+1}Q_j} |[\nabla u]_h|^q \ dz + \sup_{t_1,t_2 \in {2^{m+1}I_j}} \left|\frac{\avgs{\mu}{{[u]_h}}(t_2)-\avgs{\mu}{{[u]_h}}(t_1)}{2^{m+1}r_j} \right|^q\rbr^{\frac1{q}}.
    \end{array}
   \end{equation}

Since $B(x_j,2^{m+1}r_j) \subset \Omega$, we can apply Lemma \ref{lemma_crucial_2} with the test function $\phii(x)=\mu(x)$ and $\varphi(t) \equiv 1$, which  for any $t_1,t_2 \in \frac34 I_j$, gives
\betagin{equation}
\lambdabel{S_1_2}
 \betagin{array}{ll}
  |\avgs{\mu}{{[u]_h}}(t_2)-\avgs{\mu}{{[u]_h}}(t_1)| & \apprle  \|\nabla \mu \|_{L^{\infty}} \iint_{2^{m+1}Q_j} [|\nabla u|^{p-1}+|\theta|^{p-1}]_h \ dz\\

  & \apprle  2^{m+1} r_j \gamma \Yint\tiltlongdash_{2^{m+1}Q_j} [|\nabla u|^{p-1}+|\theta|^{p-1}]_h \ dz.

 \end{array}
\end{equation}

If $m \leq 3$, then we rescale to reach $16Q_j$ and use \deltascref{W4}, and  if $m \geq 4$, then trivially we can use \deltascref{W4}. Thus combining \eqref{S_1_1} and \eqref{S_1_2} along with \deltascref{W4}, we get
\betagin{equation}
 \lambdabel{S_1_3}
 S_1^m \apprle 2^{m+1} r_j \lbr \lambda^q + (\lambda^{p-1}\gamma)^q +  \rbr^{\frac{1}{q}} \apprle 2^{m+1}r_j \lambda ,
\end{equation}
where we have used $\gamma = \lambda^{2-p}$.

   \item[Estimate for $S_2$:] This term can be easily estimated as in the case $r_j \geq \rhoo$ to obtain
   \betagin{equation}
   \lambdabel{S_2_1}
    \betagin{array}{ll}
     \avgs{2^{k_0-2}Q_j}{{[u]_h}} 

     & \apprle \rhoo \lambda. 
    \end{array}
   \end{equation}
   
  \end{description}
Substituting \eqref{S_1_3} and \eqref{S_2_1} into \eqref{est_2} and making use of \eqref{def_k_0}, we get
\betagin{equation}
 \lambdabel{east_2}
 \betagin{array}{ll}
  |\tilde{u}h^j|  \apprle\lambda  \sum_{m=0}^{k_0-2} 2^{m+1}r_j + \lambda \rhoo 
 \apprle \lambda  \sum_{m=0}^{k_0-2} 2^{m+1}2^{1-k_0}\rhoo + \lambda \rhoo 

  \leq C_{(n,p,q,\Lambda_0,\Lambda_1,b_0)} \rhoo \lambda. 
 \end{array}
\end{equation}
\end{description}
Thus \eqref{est_1} and \eqref{east_2} imply  \eqref{claim_bound} holds, which proves the lemma.
\end{proof}

Now we prove a sharper estimate over $\Theta_1$ from \eqref{theta_1}.
\betagin{lemma}
\lambdabel{lemma3.8}
 For any $i \in \Theta_1$ and any $j \in A_i$, there holds
 \betagin{equation*}
  \lambdabel{3.26}
  |\tilde{u}h^i - \tilde{u}h^j| \leq C_{(n,p,q,\Lambda_0,\Lambda_1,b_0)} \min\{\rhoo, r_i\} \lambda.
 \end{equation*}

\end{lemma}

\betagin{proof}
We only have to consider the case $r_i \leq \rhoo$, because if $\rhoo \leq r_i$, we can directly use Lemma \ref{lemma3.7} to get the desired conclusion.

If either $\tilde{u}h^i = 0$ or $\tilde{u}h^j = 0$, then using \deltascref{W4}, \deltascref{W14} and Theorem \ref{sobolev-poincare}, we get

\betagin{equation}
\lambdabel{lemma3.8.1}
 \betagin{array}{ll}
  |\tilde{u}h^j| & \apprle r_j \lbr \Yint\tiltlongdash_{16Q_j} \left| \frac{\tilde{u}h}{r_j} \right|^q \ dz \rbr^{\frac{1}{q}} \apprle r_j \lbr \Yint\tiltlongdash_{16Q_j} \left| {\nabla \tilde{u}h} \right|^q \ dz \rbr^{\frac{1}{q}} \apprle r_j \lambda. 
 \end{array}
\end{equation}
Combining \eqref{lemma3.8.1} and \eqref{claim_bound}, the lemma follows in this case.

Now let us consider the case of $\tilde{u}h^i \neq 0$ and $\tilde{u}h^j \neq 0$. This implies $\frac34Q_i \subset \mathcal{C}_{8\rhoo}$ and $\frac34Q_j \subset \mathcal{C}_{8\rhoo}$ along with $j \in \Theta_1$. Using triangle inequality and \deltascref{W12}, we get
\betagin{equation}
 \lambdabel{3.29}
 \betagin{array}{ll}
 |\tilde{u}h^i - \tilde{u}h^j| 

 & \apprle \frac{|\frac34Q_i|}{|\frac34Q_i \cap \frac34Q_j|}\Yint\tiltlongdash_{\frac34Q_i} |\tilde{u}h(z) - \tilde{u}h^i| \ dz + \frac{|\frac34Q_j|}{|\frac34Q_i \cap \frac34Q_j|}\Yint\tiltlongdash_{ \frac34Q_j} |\tilde{u}h(z) - \tilde{u}h^j| \ dz \\
 & \apprle \Yint\tiltlongdash_{\frac34Q_i} |\tilde{u}h(z) - \tilde{u}h^i| \ dz + \Yint\tiltlongdash_{ \frac34Q_j} |\tilde{u}h(z) - \tilde{u}h^j| \ dz.
 \end{array}
\end{equation}

Let us now estimate each of the terms in \eqref{3.29} as follows: we apply H\"older's inequality followed by  Lemma \ref{lemma_crucial_1} with $\mu \in C_c^{\infty}(\frac34B_i)$ satisfying $|\mu(x)| \apprle \frac{1}{r_i^n}$ and $|\nabla \mu(x)| \apprle \frac{1}{r_i^{n+1}}$,  to get
\betagin{equation}
 \lambdabel{3.30}
 \betagin{array}{ll}
 \Yint\tiltlongdash_{\frac34Q_i} | \tilde{u}h(z) - \tilde{u}h^i| \ dz 
 & \apprle r_i \lbr \Yint\tiltlongdash_{\frac34Q_i} |\nabla \tilde{u}h|^q \ dz + \sup_{t_1,t_2 \in \frac34I_i} \left|\frac{\avgs{\mu}{\tilde{u}h}(t_2) - \avgs{\mu}{\tilde{u}h}(t_1)}{r_i} \right|^q \rbr^{\frac{1}{q}}.
 \end{array}
\end{equation}

The first term on the right of \eqref{3.30} can be controlled using \deltascref{W4}.  To control the second term of \eqref{3.30}, recalling  \eqref{def_u_tilde}, we apply Lemma \ref{lemma_crucial_2} with $\phii(x) = \eta(x) \mu(x)$ and $\varphi(t) = \zetata(t) \equiv 1$ on $\frac34I_i$ as a test function. Further making  use of \eqref{bound_b}, we get
\betagin{equation}
 \lambdabel{3.31}
 \betagin{array}{ll}
  |\avgs{\mu}{\tilde{u}h}(t_2) - \avgs{\mu}{\tilde{u}h}(t_1)| & \apprle \| \nabla (\eta \mu) \|_{\infty} \iint_{\frac34Q_i} [|\nabla u|^{p-1}+ |\theta|^{p-1}]_h \ dz \\

  & \overset{\redlabel{3.24.a}{a}}{\apprle} \lbr \frac{1}{\rhoo r_i^n} + \frac{1}{r_i^{n+1}}  \rbr |16Q_i| \lambda^{p-1} 

   \overset{\redlabel{3.24.b}{b}}{\apprle}   \lambda r_i.
 \end{array}
\end{equation}
To obtain \redref{3.24.a}{a}, we again made use of \eqref{bound_b}, \deltascref{W4},  \eqref{cut_off_function} and the structure of $\mu$. To obtain \redref{3.24.b}{b}, we used $r_i \leq \rhoo$ and $\gamma = \lambda^{2-p}$. 

Combining \eqref{3.30} and \eqref{3.31}, we get in the case $i \in \Theta_1$, the bound
\betagin{equation}
 \lambdabel{bound_when_i_1}
 \Yint\tiltlongdash_{\frac34Q_i} | \tilde{u}h(z) - \tilde{u}h^i| \ dz  \apprle r_i \lambda.
\end{equation}

This completes the proof of the lemma. 
\end{proof}

\betagin{lemma}
\lambdabel{lemma3.9}
 Given any  $z \in \elam^c$, we have $z \in \frac34Q_i$ for some $i \in \mathbb{N}$. If either  $i \in \Theta_1$ or $i \in \Theta_2$ with $\rhoo \leq r_i$, then there holds
 \betagin{equation}
  \lambdabel{3.34}
  |\nabla \vlh(z)| \leq C_{(n,p,q,\Lambda_0,\Lambda_1,b_0)} \lambda.
 \end{equation}

\end{lemma}

\betagin{proof}
 Let us prove each of the cases separately:
 \betagin{description}[leftmargin=*]
  \item[Case $i\in \Theta_1$:] In this case, we observe that $\sum_{j} \omega_j(z) = \sum_{j\in A_i} \omega_j(z)=1$ for any $z \in \elam^c$, which implies
  \betagin{equation}
   \lambdabel{3.35}
  \sum_j \nabla \omega_j(z) =  \nabla \lbr \sum_j \omega_j(z) \rbr = 0 \qquad \forall \ z \in \elam^c.
  \end{equation}

  Thus using \eqref{lipschitz_extension} and  \eqref{3.35}, we get
  \betagin{equation*}
   \lambdabel{3.36}
   \betagin{array}{ll}
   \nabla \vlh(z) = \sum_j \nabla \omega_j(z) \tilde{u}h^j 
   & = \sum_{ j\in A_i} \nabla \omega_j(z) \lbr \tilde{u}h^j - \tilde{u}h^i\rbr. 
   \end{array} 
  \end{equation*}

  Now making use of \deltascref{W9} and  \deltascref{W13} along with \eqref{theta} and Lemma \ref{lemma3.8}, we get
  \betagin{equation*}
  \lambdabel{3.37}
   |\nabla \vlh(z)| \apprle \sum_{ j\in A_i} \frac{1}{r_j} \min\{r_i, \rhoo\} \lambda \apprle \lambda.
  \end{equation*}

  \item[Case $i \in \Theta_2$ and $\rhoo \leq r_i$:] This is a simpler case, as we can directly use \eqref{lipschitz_extension} and  \eqref{claim_bound} along with \deltascref{W13} to get
  \betagin{equation*}
   \lambdabel{3.38}

    |\nabla \vlh(z)|  \apprle \sum_{ j\in A_i} |\nabla \omega_j(z)| |\tilde{u}h^j| \apprle \sum_{ j\in A_i} \frac{1}{r_j} \rhoo \lambda  \apprle \lambda.

  \end{equation*}

 \end{description}
 This completes the proof of the Lemma. 
\end{proof}

\subsection{Estimates on the spatial derivative of \texorpdfstring{$\vlh$}.}
In this subsection, we shall prove a few useful estimates on the spatial derivative of $\vlh(z)$ for $z \in \elam^c$. 

\betagin{lemma}
 \lambdabel{lemma3.10.1}
 Let $z \in \elam^c$ and $\varepsilon >0$ be any number, then there exists a constant $C(n,p,q,\Lambda_0,\Lambda_1,b_0)$ such that the following holds:
 \betagin{align}
  |\vlh(z)| & \leq C \Yint\tiltlongdash_{4Q_i}|\tilde{u}h(y,s)| \ dy\ ds \leq  \frac{Cr_i\lambda}{\varepsilon} + \frac{C\varepsilon}{\lambda r_i} \Yint\tiltlongdash_{4Q_i}|\tilde{u}h(y,s)|^2 \ dy\ ds, \lambdabel{lemma3.10_bound1}\\
  |\nabla \vlh(z)| &\leq C \frac{1}{r_i} \Yint\tiltlongdash_{4Q_i}|\tilde{u}h(y,s)| \ dy\ ds \leq  \frac{C \lambda}{\varepsilon} + \frac{C\varepsilon}{\lambda r_i^2} \Yint\tiltlongdash_{4Q_i}|\tilde{u}h(y,s)|^2 \ dy\ ds\lambdabel{lemma3.10_bound2}.

 \end{align}
\end{lemma}
\betagin{proof} Let us prove each of the estimates as follows:
 \betagin{description}[leftmargin=*]
  \item[Proof of \eqref{lemma3.10_bound1}:] Making use of \deltascref{W9} along with \eqref{lipschitz_extension}, \deltascref{W13} and \deltascref{W14}, we get
  \betagin{equation*}

    |\vlh(z)| 
     \apprle \sum_{j \in A_i} |\omega_j(z)| \Yint\tiltlongdash_{\frac34Q_j} |\tilde{u}h(y,s)|  \ dy \ ds 
     \apprle \Yint\tiltlongdash_{4Q_i} |\tilde{u}h(y,s)|  \ dy \ ds.

  \end{equation*}
  To get the second inequality in \eqref{lemma3.10_bound1}, we use Young's inequality to get
  \betagin{equation*}
   \Yint\tiltlongdash_{4Q_i} |\tilde{u}h(y,s)|  \sqrt{\frac{\lambda r_i}{\varepsilon}} \sqrt{\frac{\varepsilon}{\lambda r_i}}\ dy \ ds \apprle \frac{r_i\lambda}{\varepsilon} + \frac{\varepsilon}{\lambda r_i} \Yint\tiltlongdash_{4Q_i}|\tilde{u}h(y,s)|^2 \ dy\ ds.
  \end{equation*}
\item[Proof of \eqref{lemma3.10_bound2}:] This again follows exactly as the bound for \eqref{lemma3.10_bound1}, but in this case we additionally make use of the fact that $|\nabla \omega_j| \apprle \frac{1}{r_j}$ along with \deltascref{W5} to get the desired conclusion. 
 \end{description}
 This proves the lemma.
\end{proof}

\betagin{lemma}
 \lambdabel{lemma3.10.2}
 Let $z \in \elam^c$ and $\varepsilon \in (0,1)$ be given.  From Lemma \ref{whitney_decomposition}, we have $z \in \frac34Q_i$ for some $i \in \mathbb{N}$. Suppose $i \in \Theta_1$, then there holds:
 \betagin{align}
  |\vlh(z)|& \leq C_{(n,p,q,\Lambda_0,\Lambda_1,b_0)} \lbr \min\{ \rhoo, r_i\} \lambda + |\tilde{u}h^i| \rbr \leq  C_{(n,p,q,\Lambda_0,\Lambda_1,b_0)} \lbr \frac{ r_i\lambda}{\varepsilon} + \frac{\varepsilon}{r_i \lambda} |\tilde{u}h^i|^2 \rbr  \lambdabel{lemma3.10_bound3}\\
  |\nabla \vlh(z)| &\leq C_{(n,p,q,\Lambda_0,\Lambda_1,b_0)} \frac{\lambda}{\varepsilon}.\lambdabel{lemma3.10_bound4}
 \end{align}
\end{lemma}
\betagin{proof}
 Let us prove each of the estimates as follows:
 \betagin{description}[leftmargin=*]
  \item[Proof of \eqref{lemma3.10_bound3}:] Using triangle inequality, \deltascref{W9}, \deltascref{W13} and Lemma \ref{lemma3.8}, we get
  \betagin{equation*}
    |\vlh(z)| 
     \apprle  \sum_{j : j \in A_i} \omega_j(z)|\tilde{u}h^j -\tilde{u}h^i| + |\tilde{u}h^i| 
     \apprle \min\{ \rhoo, r_i\} \lambda + |\tilde{u}h^i| \overset{\redlabel{3.33.a}{a}}{\apprle}\frac{r_i\lambda}{\varepsilon} + \frac{\varepsilon}{r_i \lambda} |\tilde{u}h^i|^2.
  \end{equation*}
  To obtain \redref{3.33.a}{a}, we made use of Young's inequality and this proves the estimate. 
  \item[Proof of \eqref{lemma3.10_bound4}:] This follows exactly as \eqref{lemma3.10_bound3}, but in this case we additionally make use of the fact that $|\nabla \omega_j| \apprle \frac{1}{r_j}$ along with \deltascref{W5} to get the desired conclusion. 
 \end{description}
 This proves the lemma.
\end{proof}

\betagin{lemma}
 \lambdabel{lemma3.10.3}
 Let $z \in \elam^c$, then from Lemma \ref{whitney_decomposition}, we have $z \in \frac34Q_i$ for some $i \in \mathbb{N}$. Suppose $i \in \Theta_2$, then there holds:
  \betagin{align}  
  | \vlh(z)| &\leq C_{(n,p,q,\Lambda_0,\Lambda_1,b_0)}\lbr  r_i \lambda + \frac{\lambda^{1-p}r_i}{s} \Yint\tiltlongdash_{4Q_i}|\tilde{u}h(y,s)|^2 \ dy\ ds \rbr,\lambdabel{lemma3.10_bound5}\\
  |\nabla \vlh(z)| &\leq C_{(n,p,q,\Lambda_0,\Lambda_1,b_0)}\lbr  \lambda + \frac{\lambda^{1-p}}{s} \Yint\tiltlongdash_{4Q_i}|\tilde{u}h(y,s)|^2 \ dy\ ds\rbr.\lambdabel{lemma3.10_bound6} 
 \end{align}
\end{lemma}
\betagin{proof}
 Since $i \in \Theta_2$, we must have from the definition of $\Theta_2$ in \eqref{theta_2} and \deltascref{W5} that $\gamma r_j^2 \apprge s$ for all $j \in A_i$,  where $\gamma = \lambda^{2-p}$. Setting $\varepsilon = 1$ in Lemma \ref{lemma3.10.1} and making use of the bound $\lambda^{2-p} r_i^2 \apprge s$ (which uses \deltascref{W5}) along with \deltascref{W13}, the proof of the lemma follows.

\end{proof}

\subsection{Estimates on the time derivative of \texorpdfstring{$\vlh$}.}

\betagin{lemma}
\lambdabel{lemma3.11}
Let $z \in \elam^c$, then from Lemma \ref{whitney_decomposition}, there exists an $i \in \mathbb{N}$ such that $z \in \frac34Q_i$. Then the following estimates for the time derivative of $\vlh$ holds:
 \betagin{align}
 |\partial_t \vlh(z)| & \leq C_{(n,p,q,\Lambda_0,\Lambda_1,b_0)} \frac{1}{\gamma r_i^2} \Yint\tiltlongdash_{4Q_i} |\tilde{u}h(y,s)| \ dy\ ds. \lambdabel{lemma3.11.bound1} 
 \end{align}
 If $i \in \Theta_1$, then we have
 \betagin{align}

   |\partial_t \vlh(z)| & \leq C_{(n,p,q,\Lambda_0,\Lambda_1,b_0)} \frac{1}{\gamma r_i^2} \min\{r_i,\rhoo\} \lambda. \lambdabel{lemma3.11.bound2}
  \end{align}
  If $i \in \Theta_2$, then there holds
  \betagin{align}
   |\partial_t \vlh(z)| & \leq C_{(n,p,q,\Lambda_0,\Lambda_1,b_0)} \frac{\rhoo \lambda}{s}   \lambdabel{lemma3.11.bound3}.
 \end{align}

\end{lemma}
\betagin{proof}
 We shall prove each of the estimates as follows:
 \betagin{description}[leftmargin=*]
  \item[Estimate \eqref{lemma3.11.bound1}:] Using \deltascref{W9} and \deltascref{W13}, we can proceed analogous to how  \eqref{lemma3.10_bound1} was obtained to get 
  \betagin{equation*}
  \betagin{array}{ll}
   |\partial_t \vlh(z)| &\apprle \sum_{j:j \in A_i} |\partial_t \omega_j(z)| |\tilde{u}h^j| 
    \apprle \frac{1}{\gamma r_i^2} \Yint\tiltlongdash_{4Q_i} |\tilde{u}h(y,s)| \ dy\ ds .
   \end{array}
  \end{equation*}

  \item[Estimate \eqref{lemma3.11.bound2}:] Since $\sum_{j} \omega_j(\tilde{z})\equiv 1$ for all $\tilde{z} \in \elam^c$, we must have $\sum_j \partial_t \omega_j(\tilde{z}) = 0$, which gives

\betagin{equation*}
\lambdabel{3.54}
 \betagin{array}{ll}
  |\partial_t \vlh(z)| &= \left|\sum_{ j\in A_i} \partial_t\omega_j(z) \tilde{u}h^j\right| \apprle \sum_{ j\in A_i} |\partial_t\omega_j(z)| |\tilde{u}h^j-\tilde{u}h^i|

   \overset{\redlabel{3.54.a}{a}}{\apprle} \frac{1}{\gamma r_i^2} \min\{ r_i,\rhoo\} \lambda.
 \end{array}
\end{equation*}
To get \redref{3.54.a}{a}, we have made use of Lemma \ref{lemma3.8} (which is applicable since $i \in \Theta_1$) along with \deltascref{W5} and \deltascref{W13}. 
  \item[Estimate \eqref{lemma3.11.bound3}] Since $i \in \Theta_2$, we automatically have $\gamma r_i^2 \apprge s$. Now making use of  \eqref{claim_bound}, \deltascref{W5} and \deltascref{W13}, we get
  \betagin{equation*}
   \betagin{array}{ll}
    |\partial_t \vlh(z)| & \apprle \sum_{j\in A_i} |\partial_t \omega_j(z)| |\tilde{u}h^j|  \apprle \sum_{j\in A_i} \frac{1}{\gamma r_j^2} \rhoo \lambda 

     \apprle \frac{\rhoo \lambda}{s}.
   \end{array}
  \end{equation*}
 \end{description}
 This proves the lemma.
\end{proof}

\subsection{Main estimates and Lipschitz continuity of \texorpdfstring{$\vlh$}. } 
\betagin{lemma}
\lambdabel{lemma3.12}
For any $\vartheta \geq 1$,  we have the following bound:
 \betagin{equation*}
  \lambdabel{3.56}
  \iint_{8Q\setminus\elam } |\vlh(z)|^{\vartheta} \ dz \leq C_{(n,p,q,\Lambda_0,\Lambda_1,b_0)} \iint_{8Q\setminus\elam } |\tilde{u}h(z)|^{\vartheta} \ dz.
 \end{equation*}

\end{lemma}

\betagin{proof}
 Since $\elam^c$ is covered by Whitney cylinders (see Lemma \ref{whitney_decomposition}), let us pick some $i \in \mathbb{N}$ and consider the corresponding parabolic Whitney cylinder. Using the construction from \eqref{lipschitz_extension} along with \deltascref{W5}, \deltascref{W9} and \deltascref{W13}, we get
 \betagin{equation}
 \lambdabel{3.57}

   \iint_{\frac34Q_i} |\vlh(z)|^{\vartheta} \ dz  \apprle \sum_{j:j \in A_i} \iint_{\frac34Q_i} \omega_j(z)^{\vartheta} |\tilde{u}h^j|^{\vartheta} 

    \apprle \iint_{4Q_i} |\tilde{u}h(z)|^{\vartheta} \ dz. 

 \end{equation}
 Summing \eqref{3.57} over all $i \in \mathbb{N}$ and using $\spt(\tilde{u}h)\subset 8Q$ along with \deltascref{W4} and \deltascref{W7}, we get
\betagin{equation}
 \iint_{8Q\setminus \elam} |\vlh(z)|^{\vartheta} \ dz \leq \sum_i\iint_{\frac34Q_i} |\vlh(z)|^{\vartheta} \ dz \apprle \sum_i \iint_{4Q_i} |\tilde{u}h(z)|^{\vartheta} \ dz \apprle \iint_{8Q\setminus \elam} |\tilde{u}h(z)|^{\vartheta} \ dz.
\end{equation}
This proves the Lemma.
\end{proof}

\betagin{lemma}
\lambdabel{lemma3.13}
 For any  $i \in \Theta_1$,  we have the following estimate: 
 \betagin{equation}
 \lambdabel{3.59}
  \Yint\tiltlongdash_{\frac34Q_i} \left|\frac{\tilde{u}h(z) - \vlh(z)}{r_i}\right|^q \ dz \leq C_{(n,p,q,\Lambda_0,\Lambda_1,b_0)} \lambda^q.
 \end{equation}

\end{lemma}

\betagin{proof}
 For any  $z \in \frac34Q_i$, using \eqref{lipschitz_extension} along with triangle inequality and  \deltascref{W9}, we get
 \betagin{equation}
 \lambdabel{3.60}

    \Yint\tiltlongdash_{\frac34Q_i}|\tilde{u}h(z) - \vlh(z)|^q \ dz 

   \leq \Yint\tiltlongdash_{\frac34Q_i}|\tilde{u}h(z) - \tilde{u}h^i|^q \ dz +  \sum_{j:j \in A_i}\Yint\tiltlongdash_{\frac34Q_i} \left|\tilde{u}h^j - \tilde{u}h^i \right|^q\  dz := J_1 + J_2.

 \end{equation}

 We shall estimate each of the terms of \eqref{3.60} as follows (note that $i \in \Theta_1$): 
 \betagin{description}[leftmargin=*]
  \item[Estimate for $J_1$:] If $\tilde{u}h^i \neq 0$, then  $J_1$ is exactly as in \eqref{3.30}, which implies 

  \betagin{equation}
  \lambdabel{3.61}
   J_1 \apprle (r_i \lambda)^q. 
  \end{equation}
   If $\tilde{u}h^i = 0$, then we can proceed as \eqref{lemma3.8.1} to again bound $J_1$ by \eqref{3.61}.

  \item[Estimate for $J_2$:] To bound this term, we can directly use Lemma \ref{lemma3.8} to get 
  \betagin{equation}
  \lambdabel{3.62}
   J_2 \apprle (r_i \lambda)^q. 
  \end{equation}
 \end{description}
 Substituting \eqref{3.61} and \eqref{3.62} into \eqref{3.60} and making use of  \deltascref{W13}, the lemma follows.

\end{proof}

\betagin{lemma}
\lambdabel{lemma2.5}
Given any  $i \in \Theta_2$ and $\alpha_0$ as in \eqref{def_s}, there holds:
 \[
   \Yint\tiltlongdash_{\frac34Q_i} \left| \frac{\tilde{u}h - \tilde{u}h^i}{r_i} \right|^q \ dz \leq C_{(n,p,q,\Lambda_0,\Lambda_1,b_0)} \lambda^q \lbr 1  + \lbr \frac{\lambda^{2-p}}{\alpha_0^{2-p}}\rbr^{\frac{n+1}{2}}\rbr^q.
 \]

\end{lemma}

\betagin{proof}
 Since $i \in \Theta_2$, we have $\gamma r_i^2 \apprge s$. Let us split into cases $\tilde{u}h^i \neq 0$ and $\tilde{u}h^i = 0$:
 \betagin{description}[leftmargin=*]
  \item[Case $\tilde{u}h^i \neq 0$:] From the construction in \eqref{lipschitz_extension}, we have  $\frac34Q_i \subset \mathcal{C}_{8\rhoo}$, which also implies the bound $r_i \leq 8\rhoo$. Now applying Lemma  \ref{lemma_crucial_1} with $\mu(x) \in C_c^{\infty} (\frac34B_i)$ satisfying $|\mu| \leq \frac{C(n)}{r_i^n}$ and $|\nabla \mu| \leq \frac{C(n)}{r_i^{n+1}}$,  we get
 \betagin{equation}\lambdabel{lemma2.5.8.1}
  \betagin{array}{rcl}
   \Yint\tiltlongdash_{\frac34Q_i} \left| \frac{\tilde{u}h(z) - \tilde{u}h^i}{r_i} \right|^q \ dz 
   & \apprle& \Yint\tiltlongdash_{\frac34Q_i} |\nabla \tilde{u}h|^q \ dz + \sup_{t_1,t_2 \in \frac34I_i} \left| \frac{\avgs{\mu}{\tilde{u}h}(t_1)-\avgs{\mu}{\tilde{u}h}(t_2)}{r_i} \right|^q .

  \end{array}
 \end{equation}
 The first term of \eqref{lemma2.5.8.1} can be estimated using \deltascref{W4}. To estimate the second term of \eqref{lemma2.5.8.1}, let us apply Lemma \ref{lemma_crucial_2} with $\phii(x) = \mu(x) \eta(x)$ and $\varphi(t) = \zetata(t)$, which combined with \eqref{cut_off_function}, gives for any $t_1,t_2 \in \frac34I_i$, the estimate
  \betagin{equation}
  \lambdabel{lemma2.5.8}
  \betagin{array}{ll}
  |\avgs{\mu}{\tilde{u}h}(t_2)- \avgs{\mu}{\tilde{u}h}(t_1)| 

   & \apprle \|\nabla (\eta \mu) \|_{L^{\infty}(\frac34B_i)}   \iint_{\frac34Q_i} [|\nabla u|^{p-1}+|\theta|^{p-1}]_h\ dz  \\
   & \qquad +\  \|\mu\|_{L^{\infty}(\frac34B_i)}   \iint_{\frac34Q_i\cap 8Q} \frac{[u]_h}{s} \ dz \\
    &:= J_1 + J_2.
  \end{array}
 \end{equation}
  
  
  We can easily control $J_1$ of \eqref{lemma2.5.8} using \deltascref{W4} and  $r_i \leq 8\rhoo$ to get
  \betagin{equation}
   \lambdabel{lemma2.5.9}

  J_1  \apprle \lbr \frac{1}{\rhoo r_i^n} + \frac{1}{ r_i^{n+1}}\rbr|Q_i| \lambda^{p-1}  \apprle r_i \lambda.

  \end{equation}

To control $J_2$ of  \eqref{lemma2.5.8}, we proceed as follows:
\betagin{equation}
\lambdabel{lemma2.5.9.1}

  J_2 \overset{\text{\redlabel{3.50.a}{a}}}{\apprle}  \frac{\rhoo}{r_i^n s} |8Q|  \lbr \Yint\tiltlongdash_{8Q} \lbr \frac{[u]_h}{\rhoo}\rbr^q \ dz\rbr^{\frac1{q}} \overset{\text{\redlabel{3.50.b}{b}}}{\apprle} \frac{\rhoo^{n+1}}{r_i^n}   \lbr \Yint\tiltlongdash_{8Q} |\nabla [u]_h|^q \ dz\rbr^{\frac1{q}} \overset{\text{\redlabel{3.50.c}{c}}}{\apprle} 

   \lbr \frac{\gamma}{\alpha_0^{2-p}}\rbr^{\frac{n+1}{2}} r_i \lambda.

\end{equation}
To obtain \redref{3.50.a}{a}, we made use of H\"older's inequality. To get \redref{3.50.b}{b}, we made use of Theorem \ref{sobolev-poincare}. Finally  to get \redref{3.50.c}{c}, we used \deltascref{W4}, $\gamma r_i^2 \apprge s$ and \eqref{def_s}.

  
  
  Thus substituting \eqref{lemma2.5.8},  \eqref{lemma2.5.9} and \eqref{lemma2.5.9.1} into \eqref{lemma2.5.8.1}, the Lemma follows for case $\tilde{u}h^i \neq 0$.
  \item[Case $\tilde{u}h^i = 0$:]   In this case, we can directly apply Theorem \ref{sobolev-poincare} and make use of \deltascref{W4} to get 
 \betagin{equation*}\lambdabel{lemma2.5.6}
  \betagin{array}{ll}
   \Yint\tiltlongdash_{\frac34Q_i} \left| \frac{\tilde{u}h}{r_i} \right|^q \ dz \apprle \Yint\tiltlongdash_{\frac34Q_i} |\nabla \tilde{u}h|^q \apprle \lambda^q.
  \end{array}
 \end{equation*}
 \end{description}
 This completes the proof of the Lemma. 
 \end{proof}

\betagin{lemma}
\lambdabel{lemma3.14}
 For any $0< \vartheta \leq q$, there exists a positive constant $C(n,p,q,\Lambda_0,\Lambda_1,b_0, \vartheta)>0$ such that the following holds:
 \betagin{equation*}
  \iint_{8B \times 2I \setminus \elam} |\partial_t \vlh(z)  (\vlh(z) - \tilde{u}h(z))|^{\vartheta} \ dz \leq C \lambda^{\vartheta p} |\mathbb{R}^{n+1} \setminus \elam| + \frac{C}{s^{\vartheta}} \iint_{8Q} |\tilde{u}h(z)|^{2\vartheta} \ dz. 
 \end{equation*}

\end{lemma}

\betagin{proof}
From \deltascref{W3} and \eqref{theta}, we see that 

 $8B \times 2I \setminus \elam\subset \sum_{i \in \Theta} 4Q_i$. Thus to prove the lemma, it will be sufficient to  consider the two sub-cases, one where $i \in \Theta_1$ and the other where $i \in \Theta_2$. Thus, for a given $i \in \Theta$, let us define the following: 
 \betagin{equation*}
  J_i:= \iint_{\frac34Q_i} |\partial_t \vlh(z)  (\vlh(z) - \tilde{u}h(z))|^{\vartheta} \ dz .
 \end{equation*}

 \betagin{description}[leftmargin=*]
  \item[Case $i\in \Theta_1$:] Making use of  \eqref{lemma3.11.bound2} and H\"older's inequality applied to \eqref{3.59} (recall $\gamma = \lambda^{2-p}$), we get
  \betagin{equation}
  \lambdabel{3.66}

J_i  \apprle \lbr \frac{1}{\gamma r_i^2} r_i \lambda \rbr^{\vartheta} \iint_{\frac34Q_i} |\vlh(z) - \tilde{u}h(z)|^{\vartheta} \ dz 

 \apprle \lambda^{\vartheta p} \ |\frac34Q_i|.

  \end{equation}

  \item[Case $i\in \Theta_2$:] In this case, we obtain the following sequence of estimates: 
  \betagin{equation}
  \lambdabel{3.67}
   \betagin{array}{ll}
    J_i & \overset{\redlabel{3.52.a}{a}}{\apprle} \lbr \frac{1}{\gamma r_i^2} \rbr^{\vartheta} \lbr \Yint\tiltlongdash_{4Q_i} |\tilde{u}h(z)| \ dz \rbr^{\vartheta} \lbr \iint_{\frac34Q_i} |\vlh(z)| + |\tilde{u}h(z)| \ dz \rbr^{\vartheta} \\
    & \overset{{\redlabel{3.52.b}{b}}}{\apprle} \lbr \frac{1}{\gamma r_i^2} \rbr^{\vartheta} \lbr \Yint\tiltlongdash_{4Q_i} |\tilde{u}h(z)| \ dz \rbr^{\vartheta} \lbr |\frac34Q_i|\Yint\tiltlongdash_{4Q_i} |\tilde{u}h(z)| \ dz + \iint_{4Q_i}|\tilde{u}h(z)| \ dz \rbr^{\vartheta} \\

    & \overset{{\redlabel{3.52.c}{c}}}{\apprle} \frac{1}{s^{\vartheta}}  \lbr \iint_{4Q_i} |\tilde{u}h(z)|^{2\vartheta} \ dz \rbr.
   \end{array}
  \end{equation}
  To get \redref{3.52.a}{a}, we made use of \eqref{lemma3.11.bound1}. To obtain \redref{3.52.b}{b} above, we  made use of \eqref{lemma3.10_bound1}. Finally to get \redref{3.52.c}{c}, we have made use of the fact that $i\in\Theta_2$, which implies $\gamma r_i^2 \apprge s$. 

 \end{description}

 Combining \eqref{3.66} and \eqref{3.67} followed by summing over all $i \in \Theta$ and  using  \deltascref{W7}, we get 
 \betagin{equation*}
  \betagin{array}{ll}
   \iint_{8B \times 2I \setminus \elam} |\partial_t \vlh(z)  (\vlh(z) - \tilde{u}h(z))|^{\vartheta} \ dz 
   & \apprle \sum_{i \in \Theta_1} J_i + \sum_{i \in \Theta_2} J_i \\

   & \overset{}{\apprle} \lambda^{\vartheta p} |\mathbb{R}^{n+1} \setminus \elam| + \frac{1}{s^{\vartheta}} \iint_{8Q} |\tilde{u}h(z)|^{2\vartheta} \ dz. 
  \end{array}
 \end{equation*}

This completes the proof of the lemma. 
\end{proof}

We shall now prove the Lipschitz continuity of $\vlh$ on $2\mathcal{H} = \mathbb{R}^n \times (-4s, 4s)$. {\color{black}The method of making use of Lemma \ref{metric_lipschitz} to prove Lemma \ref{lemma3.15} is due to \cite[Chapter 3]{bogelein2013regularity}.}	
\betagin{lemma}
\lambdabel{lemma3.15}
The extension $\vlh$ from  \eqref{lipschitz_extension} is $C^{0,1}(2\mathcal{H})$ with respect to the parabolic metric $d_{\lambda}$. 
\end{lemma}

\betagin{proof} 
Let us consider a parabolic cylinder $\tilde{Q}_{r}(z) := \tilde{Q}_{r, \gamma r^2} (z) := \tilde{Q}$ for some $z \in 2\mathcal{H}$ and $r>0$. To prove the Lemma, we make use of Lemma \ref{metric_lipschitz} to bound
 \betagin{equation}
 \lambdabel{bound_I_r}
I_r(z) :=   \Yint\tiltlongdash_{\tilde{Q} \cap 2\mathcal{H}} \left|\frac{\vlh(\tilde{z}) - \avgsnoleft{\tilde{Q}\cap2\mathcal{H}}{\vlh}}{r}\right| \ d\tilde{z} \leq {o}(1),
 \end{equation}
where ${o}(1)$ denotes a bound independent of $z \in 2\mathcal{H}$ and $r>0$ only. 

In order to do this, we split the proof into three cases given by

\betagin{gather*}
 (i)\ 2\tilde{Q} \subset 4\mathcal{H} \setminus \elam, \quad (ii) \ 2\tilde{Q} \subset 4\mathcal{H} \quad  \text{with} \ \gamma r^2 \leq 4s, \\
(iii)\ 2\tilde{Q} \subset 4\mathcal{H} \quad \text{with}\  \gamma r^2 > 4s  \quad \text{\bf OR}  \quad 2\tilde{Q} \cap (4\mathcal{H})^c \neq \emptyset. 
\end{gather*}

\betagin{description}[leftmargin=*]
 \item[Case $2\tilde{Q} \subset 4\mathcal{H} \setminus \elam$:]  In this case, we observe the following lower bound must hold:

 \betagin{equation}
  \lambdabel{3.70}
  |\tilde{Q} \cap 2\mathcal{H}| \apprge |\tilde{Q}| \approx r^{n+2} \gamma. 
 \end{equation}
 From the construction of $\vlh$ in \eqref{lipschitz_extension}, we see that  $\vlh$ is smooth on $\elam^c$. 
 
 
 If $2\tilde{Q} \subset  \mathcal{C}_{8\rhoo}$, then setting $\tilde{Q} = \tilde{B}_r \times \tilde{I}_{\gamma r^2}$ and  using the definition of derivatives, we get
  \betagin{equation}
   \lambdabel{3.71}
   \deltaf2.2{2.2}

I_r(z)  \apprle   \frac{1}{r} \Yint\tiltlongdash_{\tilde{Q} \cap 2\mathcal{H}} \Yint\tiltlongdash_{\tilde{Q} \cap 2\mathcal{H}}  |\vlh(\tilde{z}) - \vlh(\tilde{z}_1)| \ d\tilde{z} \ d\tilde{z}_1
 \apprle  \sup_{\tilde{z} \in \tilde{Q} \cap 2\mathcal{H}}\lbr  |\nabla\vlh(\tilde{z})| + r |\partial_t \vlh(\tilde{z})|\rbr .

  \end{equation}
  
  On the other hand, if $2\tilde{Q} \cap \mathcal{C}_{8\rhoo}^c \neq \emptyset$, then applying Theorem \ref{sobolev-poincare} on each times slice, we obtain
  \betagin{equation}
   \lambdabel{3.72}

    I_r(z)  \apprle \Yint\tiltlongdash_{\tilde{Q} \cap 2\mathcal{H}} \left|\frac{\vlh(\tilde{z})}{r}\right| \ d\tilde{z} \apprle \Yint\tiltlongdash_{\tilde{Q} \cap 2\mathcal{H}} \left| \nabla \vlh(\tilde{z})\right| \ d\tilde{z}
    \apprle \sup_{\tilde{z} \in \tilde{Q}\cap 2\mathcal{H}} |\nabla \vlh(\tilde{z})|.

  \end{equation}

From \eqref{3.71} and \eqref{3.72}, we see that in order to  bound $I_r(z)$ independent of $z$ and $r$, we have to bound the following expression independent of $z$ and $r$:

 \betagin{equation}\lambdabel{3.72.1}\sup_{\tilde{z} \in \tilde{Q} \cap 2\mathcal{H}} \lbr |\nabla\vlh(\tilde{z})| + r |\partial_t \vlh(\tilde{z})|\rbr.\end{equation}
 

To this end, let  $\tilde{z} \in \tilde{Q} \cap 2\mathcal{H}$, which implies  $\tilde{z} \in \frac34Q_i$ for some $i \in \mathbb{N}$. Since $\tilde{Q} \subset \elam^c$, we also get
\betagin{equation}
 \lambdabel{3.74}
 r \leq \partialrd (\tilde{z}, \elam) \leq \partialrd (\tilde{z}, z_i) + \partialrd (z_i , \elam) \leq r_i + 16r_i = 17r_i. 
\end{equation}

In order to bound \eqref{3.72.1}, we consider the following three sub-cases:
\betagin{description}
\item[Subcase $i \in \Theta_1$:] Using \eqref{3.34} and  \eqref{lemma3.11.bound2} along with  \eqref{3.74}, we get 
\betagin{equation*}
\lambdabel{3.76}
|\nabla \vlh(\tilde{z})| + r |\partial_t \vlh(\tilde{z})| \apprle \lambda + r \frac{1}{\gamma r_i^2} \min\{ r_i,\rhoo\} \lambda \apprle \lambda + \lambda^{p-1} \leq o(1). 
\end{equation*}

\item[Subcase $i \in \Theta_2$ and $\rhoo \leq r_i$:] In this case, we make use of \eqref{3.34},  \eqref{lemma3.11.bound1} and \eqref{3.74} to get
\betagin{equation*}
\lambdabel{3.77}
\betagin{array}{ll}
|\nabla \vlh(\tilde{z})| + r |\partial_t \vlh(\tilde{z})| & \apprle \lambda + r \frac{1}{\gamma r_i^2} \frac{1}{|Q_i|} \iint_{4Q_i} |\tilde{u}h(y,s)| \ dy\ ds \\

& \overset{}{\apprle } \lambda +  \frac{1}{\gamma^2 \rhoo^{n+3}} \| \tilde{u}h\|_{L^1(\mathbb{R}^{n+1})} \leq o(1).
\end{array}
\end{equation*}

\item[Subcase $i \in \Theta_2$ and $\rhoo \geq r_i$:]  In this case, we use the fact $\gamma r_i^2 \apprge s$ along with \eqref{lemma3.10_bound2} and \eqref{lemma3.11.bound3}  to obtain
\betagin{equation*}
 \lambdabel{3.79}
 \betagin{array}{ll}
  |\nabla \vlh(\tilde{z})| + r |\partial_t \vlh(\tilde{z})| & \apprle \frac{1}{r_i|Q_i|} \iint_{4Q_i} |\tilde{u}h(y,s)| \ dy\ ds + r \frac{\rhoo \lambda}{s} \\

  & \overset{}{\apprle } \frac{1}{\lbr \frac{s}{\gamma}\rbr^{\frac{n+1}{2}} s} \| \tilde{u}h\|_{L^1(\mathbb{R}^{n+1})} + \frac{\rhoo^2 \lambda}{s} \leq o(1).
 \end{array}
\end{equation*}

\end{description}
 

 We have shown that $I_r(z) \leq o(1)$ when $2\tilde{Q} \subset 4\mathcal{H} \setminus \elam$. 

 \item[Case $2\tilde{Q} \subset 4\mathcal{H}$ and $\gamma r^2 \leq 4s$:] If $2\tilde{Q} \subset \elam^c$, then this is just the previous case.  Hence without loss of generality, we can assume 
 \betagin{equation}\lambdabel{second_case}
 2\tilde{Q} \cap \elam \neq \emptyset. 
  \end{equation}
  Noting that \eqref{3.70} must also hold in this case, we apply 

  triangle inequality and estimate $I_r(z)$ by 
 \betagin{equation}
  \lambdabel{3.81}
  \betagin{array}{ll}
   I_r(z) 
   & \leq \Yint\tiltlongdash_{\tilde{Q} \cap 2\mathcal{H}} \left| \frac{\vlh(\tilde{z}) - \tilde{u}h(\tilde{z})}{r}\right|+ \left| \frac{\tilde{u}h(\tilde{z}) - \avgs{\tilde{Q}\cap 2\mathcal{H}}{\tilde{u}h}}{r}\right| + \left| \frac{\avgs{\tilde{Q}\cap 2\mathcal{H}}{\tilde{u}h} - \avgsnoleft{\tilde{Q} \cap 2\mathcal{H}}{\vlh}}{r}\right| \ d\tilde{z} \\

   & \leq 2J_1 + J_2,
  \end{array}
 \end{equation}
 where we have set
 \betagin{gather}
  J_1:= \Yint\tiltlongdash_{\tilde{Q} \cap 2\mathcal{H}} \left| \frac{\vlh(\tilde{z}) - \tilde{u}h(\tilde{z})}{r}\right| \ dz  \txt{and} J_2 := \Yint\tiltlongdash_{\tilde{Q} \cap 2\mathcal{H}} \left| \frac{\tilde{u}h(\tilde{z}) - \avgs{\tilde{Q}\cap 2\mathcal{H}}{\tilde{u}h}}{r}\right|\ dz.\lambdabel{def_J_1_2}
 \end{gather}
We now estimate each of the terms of \eqref{def_J_1_2} as follows:

\betagin{description}[leftmargin=*]
\item[Estimate for $J_1$:] If $\tilde{Q} \subset \elam$, then $\vlh = \tilde{u}h$ which implies $J_1=0$. Hence without loss of generality, along with \eqref{second_case}, we can further assume $\tilde{Q} \cap \elam^c \neq \emptyset$. 
Using the construction of \eqref{lipschitz_extension}, we get
\betagin{equation}
 \lambdabel{3.81.1}
 \betagin{array}{ll}
  J_1 

  & \apprle \sum_{i\in \Theta} \frac{1}{|\tilde{Q}\cap2\mathcal{H}|} \iint_{\tilde{Q} \cap2\mathcal{H}\cap \frac34Q_i} \left| \frac{\tilde{u}h(\tilde{z}) - \tilde{u}h^i}{r}\right| \ d\tilde{z}.
 \end{array}
\end{equation}

Let us fix an $i \in \Theta$ and take two points $\tilde{z}_1 \in \tilde{Q} \cap \frac34Q_i$ and $\tilde{z}_2 \in \elam \cap 2\tilde{Q}$. Making use of  \deltascref{W2} along with the trivial bound $\partialrd (\tilde{z}_1, \tilde{z}_2) \leq  4r$ and $\partialrd (z_i, \tilde{z}_1) \leq 2r_i$,  we get
\betagin{equation}

\lambdabel{3.83.1}
16r_i =\partialrd(z_i,\elam) \leq \partialrd (z_i, \tilde{z}_1) + \partialrd (\tilde{z}_1, \tilde{z}_2) \leq 2r_i + 4r  \ \Longrightarrow \ 2r_i \leq r,
\end{equation}
where $z_i$ denotes the centre of $\frac34Q_i$ as in \deltascref{W1}. 

Note that \eqref{3.70} holds and thus summing over all $i \in \Theta$ such that  $\tilde{Q} \cap2\mathcal{H}\cap \frac34Q_i \neq \emptyset$ in \eqref{3.81.1} and making use of  \eqref{3.83.1}, we get
\betagin{equation*}
 \lambdabel{3.84.1}
 J_1 \apprle \sum_{\substack{i\in\Theta \\ \tilde{Q} \cap2\mathcal{H}\cap \frac34Q_i \neq \emptyset }} \frac{|\frac34Q_i|}{|\tilde{Q}\cap2\mathcal{H}|} \Yint\tiltlongdash_{\frac34Q_i} \left| \frac{\tilde{u}h(\tilde{z}) - \tilde{u}h^i}{r}\right| \ d\tilde{z} \apprle \sum_{i\in \Theta}  \Yint\tiltlongdash_{\frac34Q_i} \left| \frac{\tilde{u}h(\tilde{z}) - \tilde{u}h^i}{r_i}\right| \ d\tilde{z}.
\end{equation*}
To bound \eqref{3.84.1},  in the case $i \in \Theta_2$, we apply H\"older's inequality followed by Lemma \ref{lemma2.5} and in the case $i \in \Theta_1$, we again apply H\"older's inequality and use \eqref{3.61}. Thus we get
\betagin{equation}
\lambdabel{3.84.2}
 J_1 \apprle \lambda \lbr 1  + \lbr \frac{\lambda^{2-p}}{\alpha_0^{2-p}}\rbr^{\frac{n+1}{2}}\rbr.
\end{equation}

 \item[Estimate for $J_2$:] Note that \eqref{3.70} is valid, which we use along with triangle inequality to get

 \betagin{equation}
 \lambdabel{3.82}
  \betagin{array}{ll}
   J_2 
   & \apprle \Yint\tiltlongdash_{\tilde{Q} \cap 2\mathcal{H}}  \left| \frac{\tilde{u}h(\tilde{z}) - \avgs{\tilde{Q}}{\tilde{u}h}}{r}\right| +  \left| \frac{\avgs{\tilde{Q}}{\tilde{u}h} - \avgs{\tilde{Q}\cap 2\mathcal{H}}{\tilde{u}h}}{r}\right| \ d\tilde{z} 

    \apprle 2 \Yint\tiltlongdash_{\tilde{Q}}  \left| \frac{\tilde{u}h(\tilde{z}) - \avgs{\tilde{Q}}{\tilde{u}h}}{r}\right| \ d\tilde{z}.
  \end{array}
 \end{equation}
  If  $\tilde{Q} \subset \mathcal{C}_{8\rhoo}$, then we estimate \eqref{3.82} by first applying Lemma \ref{lemma_crucial_1} for some  $\mu \in C_c^{\infty} (\tilde{B})$ with $|\mu| \leq \frac{C(n)}{\tilde{r}^{n}}$ and  $|\nabla \mu| \leq \frac{C(n)}{\tilde{r}^{n+1}}$ to get
 \betagin{equation*}
  \lambdabel{3.83}
  \betagin{array}{ll}
   J_2 
   & \apprle \lbr \Yint\tiltlongdash_{\tilde{Q}} |\nabla \tilde{u}h(\tilde{z})|^q \ d\tilde{z} + \sup_{t_1,t_2 \in \tilde{I}}\left|\frac{\avgs{\mu}{\tilde{u}h}(t_2)-\avgs{\mu}{\tilde{u}h}(t_1)}{r} \right|^q\rbr^{\frac{1}{q}}.
  \end{array}
 \end{equation*}
 
 Since we have $\tilde{Q} \subset \mathcal{C}_{8\rhoo}$, we have  $\tilde{r} \leq 8\rhoo$. Furthermore, since  $\tilde{Q} \subset 4\mathcal{H}$, we see that  $\zetata(t) \equiv 1$ on $\tilde{Q}$. Thus,  we can apply Lemma \ref{lemma_crucial_2} and make use of the hypothesis \eqref{second_case} and mimic the procedure in the proof of Lemma \ref{lemma3.8} to obtain the bound 

 \betagin{equation}
 \lambdabel{3.88}
  J_2 \apprle \lambda.
 \end{equation}

 On the other hand,  if $\tilde{Q} \nsubseteq \mathcal{C}_{8\rhoo}$, then we can directly apply Theorem \ref{sobolev-poincare} and use  \eqref{3.82},  to get
 \betagin{equation}
 \lambdabel{3.89}
  \betagin{array}{ll}
   J_2 & \apprle \Yint\tiltlongdash_{\tilde{Q}} \left| \frac{\tilde{u}h(\tilde{z})}{r}\right| \ d\tilde{z} \apprle \lbr \Yint\tiltlongdash_{\tilde{Q}} |\nabla \tilde{u}h(\tilde{z})|^q \ d\tilde{z} \rbr^{\frac{1}{q}} \apprle \lambda. 
  \end{array}
 \end{equation}

 
\end{description}

Thus substituting either \eqref{3.88} or \eqref{3.89} along with  \eqref{3.84.2} into \eqref{3.81}, we get  $I_r(z) \leq o(1)$ in the case that $2\tilde{Q} \subset 4\mathcal{H}$ and $\gamma r^2 \leq 4s$. 

 
 \item[Final case $2\tilde{Q} \subset 4\mathcal{H}$ and $\gamma r^2 > 4s$ OR $2\tilde{Q} \cap (4\mathcal{H})^c \neq \emptyset$:] In either case, we have $\gamma r^2 \apprge s$, which gives the estimate
 \betagin{equation*}
  \lambdabel{3.91}
  \betagin{array}{ll}
  \Yint\tiltlongdash_{\tilde{Q} \cap 2\mathcal{H}}  \left| \frac{\vlh(\tilde{z}) - \avgsnoleft{\tilde{Q} \cap 2\mathcal{H}}{\vlh}}{r} \right|\ d\tilde{z} & \apprle \frac{1}{r^{n+3}} \iint_{\tilde{Q} \cap 2\mathcal{H}} |\vlh(\tilde{z})| \ d\tilde{z} \\

  & \apprle \frac{1}{r^{n+3}} \iint_{8Q \cap \elam} |\vlh(\tilde{z})| \ d\tilde{z} + \frac{1}{r^{n+3}} \iint_{8Q\setminus \elam} |\vlh(\tilde{z})| \ d\tilde{z}.
  \end{array}
 \end{equation*}

 By construction of $\vlh$ in \eqref{lipschitz_extension}, we have $\vlh = \tilde{u}h$ on $\elam$.  On $8Q \setminus \elam$, we can apply Lemma \ref{lemma3.12} (with $\vartheta=1$) to obtain the following bound:
 \betagin{equation*}
  \betagin{array}{ll}
   \Yint\tiltlongdash_{\tilde{Q} \cap 2\mathcal{H}}  \left| \frac{\vlh(\tilde{z}) - \avgsnoleft{\tilde{Q} \cap 2\mathcal{H}}{\vlh}}{r} \right|\ d\tilde{z} & \apprle \lbr \frac{\gamma}{s} \rbr^{\frac{n+3}{2}} \iint_{8Q \cap \elam} |\tilde{u}h(\tilde{z})| \ d\tilde{z} +  \lbr \frac{\gamma}{s} \rbr^{\frac{n+3}{2}} \iint_{8Q \setminus \elam} |\tilde{u}h(\tilde{z})| \ d\tilde{z} \\
   & \apprle o(1).
  \end{array}
 \end{equation*}
\end{description}
This completes the proof of the lemma. 

\end{proof}

\betagin{remark}
In the proof of Lemma \ref{lemma3.15}, we can actually show the Lipschitz constant of $\vlh$ on $2\mathcal{H}$  with respect to the parabolic metric \eqref{parabolic_metric} is $C_{(n,p,q,\Lambda_0,\Lambda_1,b_0)}\lambda$. This improved bound will not be needed for the purposes of this paper, but it is worthwhile to mention how to obtain such an improvement. 

 In order to show this improvement, we can either use sharper estimates in the proof of Lemma \ref{lemma3.15} or mimic the proof from \cite{KL}. 

 
 Since we only take Steklov average in the time variable,  a priori we don't have continuity of $\vlh$ in the space variable to be able to implement the argument of \cite{KL} directly. On the other hand, once we have Lemma \ref{lemma3.15}, we know that $\vlh$ is continuous on $2\mathcal{H}$, which enables us to implement the arguments of \cite{KL} and hence obtain the improved bound. 
 
 
\end{remark}

\betagin{lemma}
 \lambdabel{lemma3.18}
 For any $i \in \Theta$ and $k \in \{0,1\}$, there exists a positive  constant $C{(n,p,q,\Lambda_0,\Lambda_1,b_0)}$ such that there holds:
 \betagin{equation}
 \lambdabel{lemma3.18.bound_1}
  \iint_{\frac34Q_i} [|\nabla u|^{p-1}+|\theta|^{p-1}]_h |\nabla^k \vlh| \ dz \leq C \rhoo^{1-k} \lbr \lambda^p |4Q_i| + \frac{\lsb{\chii}{i \in \Theta_2}}{s} \iint_{4Q_i} |\tilde{u}h|^2 \ dz \rbr. 
 \end{equation}
Here we have used the notation  $\lsb{\chii}{i\in \Theta_2} = 1$ if $i \in \Theta_2$ and $\lsb{\chii}{i\in \Theta_2} = 0$ if $i \in \Theta_1$ and $\nabla^0 \vlh := \vlh$. 
\end{lemma}

\betagin{proof}
 We shall split the proof into two cases: 
 \betagin{description}[leftmargin=*]
  \item[Case $i\in \Theta_1$ or $i\in \Theta_2$ with $\rhoo \leq r_i$:] In this case, we make use of Lemma \ref{lemma3.7} when $k=0$ or Lemma \ref{lemma3.9} when $k=1$,  along with \deltascref{W4}, to get
  \betagin{equation*}
   \lambdabel{lemma3.18.1}
   \betagin{array}{ll}
    \iint_{\frac34Q_i}  [|\nabla u|^{p-1}+|\theta|^{p-1}]_h |\nabla^k \vlh| \ dz & \apprle \rhoo^{1-k} \lambda \iint_{\frac34Q_i}  [|\nabla u|^{p-1}+|\theta|^{p-1}]_h  \ dz\\

    &{\apprle} \rhoo^{1-k} \lambda |4Q_i| \lambda^{p-1} = \rhoo^{1-k} \lambda^p |4Q_i|.
   \end{array}
  \end{equation*}

  \item[Case $i \in \Theta_2$ with $r_i \leq \rhoo$:] If $k=0$, we make use of \eqref{lemma3.10_bound6} and in the case $k=1$, we make use of \eqref{lemma3.10_bound5}, along with \deltascref{W4} to  get
  \betagin{equation*}
   \lambdabel{lemma3.18.2}
   \betagin{array}{ll}
    \iint_{\frac34Q_i}  [|\nabla u|^{p-1}+|\theta|^{p-1}]_h |\nabla \vlh| \ dz 
    & \apprle \rhoo^{1-k}\lbr \lambda  + \frac{\lambda^{1-p}}{s} \Yint\tiltlongdash_{4Q_i} |\tilde{u}h|^2 \ dz \rbr |4Q_i| \lambda^{p-1} \\
    & \apprle \rhoo^{1-k} \lbr \lambda^p |4Q_i| + \frac{1}{s} \iint_{4Q_i} |\tilde{u}h|^2 \ dz \rbr.
   \end{array}
  \end{equation*}
 \end{description}
This completes the proof of the lemma. 
\end{proof}

\betagin{corollary}
 \lambdabel{corollary3.20}
There exists a positive constant $C{(n,p,q,\Lambda_0,\Lambda_1,b_0)}$ such that  the following estimate holds for any $k \in \{0,1\}$:
 \betagin{equation*}
  \lambdabel{corollary3.20.1}
  \iint_{8B \times 2I \setminus \elam} [|\nabla u|^{p-1}+|\theta|^{p-1}]_h |\nabla^k \vlh| \ dz \leq C \rhoo^{1-k} \lbr \lambda^p |\mathbb{R}^{n+1} \setminus \elam| + \frac{1}{s} \iint_{8Q} |\tilde{u}h|^2 \ dz \rbr. 
 \end{equation*}

\end{corollary}
\betagin{proof}
 Since  $8B \times 2I \setminus \elam \subset \sum_{i\in \Theta} \frac34Q_i$, all we need to do is to sum over $i\in \Theta$ in \eqref{lemma3.18.bound_1} and make use of the fact that $\spt(\tilde{u}h) \subset 8Q$, to prove the Corollary. 

\end{proof}

We shall now prove a crucial estimate on each time slice. 
\betagin{lemma}
 \lambdabel{pre_crucial_lemma}
 For any $i \in \Theta_1$ and any $0 < \varepsilon \leq 1$, there exists a positive constant $C{(n,p,q, \Lambda_0,\Lambda_1,b_0)}$ such that for  almost every $t \in 2I$, there holds
 \betagin{equation}
 \lambdabel{3.120}
  \left| \int_{8B} (\tilde{u}(x,t) - \tilde{u}^i) \vl(x,t) \omega_i(x,t) \ dx \right| \leq C \lbr  \frac{\lambda^p}{\varepsilon} |4Q_i| + \varepsilon |4B_i| |\tilde{u}^i|^2\rbr. 
 \end{equation}
 
 In the case $i \in \Theta_2$, there exists a positive constant $C{(n,p,q, \Lambda_0,\Lambda_1,b_0)}$ such that for  almost every $t \in 2I$, there holds
 \betagin{equation}
 \lambdabel{3.121}
  \left| \int_{8B} \tilde{u}(x,t) \vl(x,t) \omega_i(x,t) \ dx \right| \leq C\lbr  {\lambda^p} |4Q_i| + \frac{1}{s}\iint_{\frac34Q_i} |u|^2 \lsb{\chii}{8Q}\ dz\rbr. 
 \end{equation}

\end{lemma}
\betagin{proof}
Let us fix any $t \in 2I$,  $i \in \Theta$  and take $\eta(y) \zetata(\tauu) \omega_i(y,\tauu) \vlh(y,\tauu)$ as a test function in \eqref{very_weak_solution}. Further integrating the resulting expression over $ \left(t_i - \gamma \left(\frac34 r_i\right)^2 , t\right)$ along with making use of  the fact that $\omega_i(y,t_i - \gamma (3r_i/4)^2) = 0$, we get for  any $a\in \mathbb{R}$, the equality

\betagin{equation}
 \lambdabel{3.123}
 \betagin{array}{ll}
  \int_{\Omega} \lbr[(]  (\tilde{u}h - a)  \omega_i \vlh \rbr (y,t) \ dy & = \int_{t_i - \gamma \left(\frac34 r_i\right)^2}^t \int_{\Omega} \partial_t \left(  (\tilde{u}h - a) \omega_i \vlh \right) (y,\tauu) \ dy \ d\tauu \\

  & = \int_{t_i - \gamma \left(\frac34 r_i\right)^2}^t \int_{\Omega} \iprod{[\mathcal{A}(y,\tauu,\nabla u)]_h}{\nabla (\eta \zetata \omega_i \vlh)} \ dy \ d\tauu  \\
  & \qquad +\  \int_{t_i - \gamma \left(\frac34 r_i\right)^2}^t \int_{\Omega}   [u]_h  \ \partial_t\lbr \eta \zetata \omega_i \vlh \right) (y,\tauu) \ dy \ d\tauu  \\
  & \qquad -\  \int_{t_i - \gamma \left(\frac34 r_i\right)^2}^t \int_{\Omega} a \partial_t \lbr \omega_i \vlh\rbr \ dy \ d\tauu.
 \end{array}
\end{equation}

We can estimate $|\nabla (\eta \zetata \omega_i \vl)|$ using the chain rule, \eqref{cut_off_function} and \deltascref{W9}, to get
\betagin{equation}
 \lambdabel{3.126}
 \betagin{array}{ll}
 |\nabla (\eta \zetata \omega_i \vl)| 
 & \apprle \frac{1}{\rhoo} |\vl| + \frac{1}{r_i} |\vl| + |\nabla \vl|.
 \end{array}
\end{equation}
Similarly, we can estimate $\left|\partial_t\lbr \eta \zetata \omega_i \vl \right)\right|$ using the chain rule, \eqref{cut_off_function} and \deltascref{W9}, to get
\betagin{align}
 \left|\partial_t\lbr \eta \zetata \omega_i \vl \right)\right|  & \apprle \frac{\lsb{\chii}{i \in \Theta_2}}{s} |\vl| + \frac{1}{\gamma r_i^2} |\vl| + |\partial_t \vl|, \lambdabel{3.127} \\
 \left| \partial_t \lbr \omega_i \vl\rbr\right| & \apprle  \frac{1}{\gamma r_i^2} |\vl| + |\partial_t \vl|,\lambdabel{3.128}
\end{align}
where we have set $\lsb{\chii}{i \in \Theta_2} = 1$ if  $i \in \Theta_2$ and $\lsb{\chii}{i \in \Theta_2}=0$ otherwise. 

Let us now prove each of the assertions of the Lemma.
\betagin{description}[leftmargin=*]
 \item[Proof of \eqref{3.120}:] Note that $i \in \Theta_1$, which implies $\zetata(t) \equiv 1$ on $\frac34I_i$, thus taking $a=\tilde{u}h^i$ in the \eqref{3.123} followed by letting $h \searrow 0$ and making use of \eqref{3.126}, \eqref{abounded} and \eqref{bound_b},  we get

 \betagin{equation}
  \lambdabel{first_1}
  \betagin{array}{ll}
   \left| \int_{8B} \lbr[(] (\tilde{u} - \tilde{u}^i) \omega_i \vl \rbr (y,t) \ dy \right| & \apprle J_1 + J_2 + J_3,

  \end{array}
 \end{equation}
 where we have set 
 \betagin{align}
  J_1& := \frac{1}{\min\{\rhoo, r_i\}} \iint_{8Q} \lbr |\nabla u|^{p-1}+|\theta|^{p-1} \rbr  |\vl|   \lsb{\chii}{\frac34Q_i\cap 8Q} \ dy \ d\tauu, \nonumber \\
  J_2& :=  \iint_{8Q} \lbr |\nabla u|^{p-1}+|\theta|^{p-1} \rbr  |\nabla \vl|   \lsb{\chii}{\frac34Q_i\cap 8Q} \ dy \ d\tauu,\nonumber \\
  J_3&:= \iint_{8Q} |\tilde{u}-\tilde{u}^i| | \partial_t (\omega_i \vl)|  \lsb{\chii}{\frac34Q_i\cap 8Q} \ dy \ d\tauu.\lambdabel{bound_J_3_1}
 \end{align}

 Let us now estimate each of the terms as follows: 
 \betagin{description}
  \item[Bound for $J_1$:] If $\rhoo \leq r_i$,  we can directly use H\"older's inequality, Lemma \ref{lemma3.7} and \deltascref{W4}, to find that for any $\varepsilon \in (0,1]$, there holds 
  \betagin{equation}
   \lambdabel{bound_I_1_rho_leq_r_i}
   \betagin{array}{ll}
    J_1 & \apprle \lambda |Q_i| \lbr \Yint\tiltlongdash_{16Q_i} |\nabla u|^{q}+|\theta|^q   \ dy \ d\tauu\rbr^{\frac{p-1}{q}} \apprle \lambda^p |4Q_i|\leq \frac{\lambda^p}{\varepsilon} |4Q_i|.
   \end{array}
  \end{equation}

  In the case $r_i \leq \rhoo$, we make use of \eqref{lemma3.10_bound3}, \deltascref{W4} along with the fact $|Q_i| = |B_i| \times 2\lambda^{2-p} r_i^2$, to get
  \betagin{equation}
  \lambdabel{bound_I_1_rho_geq_r_i}
   \betagin{array}{ll}
    J_1
    & {\apprle} \frac{1}{r_i} \lbr \frac{r_i \lambda}{\varepsilon} + \frac{\varepsilon}{\lambda r_i} |\tilde{u}^i|^2  \rbr \lbr \iint_{4Q_i} |\nabla u|^{p-1}  +|\theta|^{p-1} \ dy \ d\tauu \rbr \\
    & {\apprle} \frac{1}{r_i} \lbr \frac{r_i \lambda}{\varepsilon} + \frac{\varepsilon}{\lambda r_i}  |\tilde{u}^i|^2  \rbr |4Q_i|\lambda^{p-1} 

     {\apprle} \frac{\lambda^p}{\varepsilon} |4Q_i| +  \varepsilon |4B_i| |\tilde{u}^i|^2.	
   \end{array}
  \end{equation}

Thus combining \eqref{bound_I_1_rho_geq_r_i} and \eqref{bound_I_1_rho_leq_r_i}, we get
\betagin{equation}
 \lambdabel{bound_I_1}
 J_1 \apprle \frac{\lambda^p}{\varepsilon} |4Q_i| +  \lsb{\chii}{r_i \leq \rhoo}\varepsilon |4B_i| |\tilde{u}h^i|^2,
\end{equation}
where we have set $\lsb{\chii}{r_i \leq \rhoo} = 1$ if  $r_i \leq \rhoo$ and $\lsb{\chii}{r_i \leq \rhoo}=0$ else. 

  \item[Bound for $J_2$:] In this case, we can directly use Lemma \ref{lemma3.9} and \deltascref{W4} to get for any $\varepsilon \in (0,1]$, the bound
  \betagin{equation}
   \lambdabel{bound_I_22}

    J_2 
     \apprle \frac{\lambda}{\varepsilon} \iint_{4Q_i}|\nabla u|^{p-1}+ |\theta|^{p-1} \ dy \ d\tauu 
   \apprle \frac{\lambda^p}{\varepsilon} |4Q_i|.

  \end{equation}

  \item[Bound for $J_3$:] Substituting  \eqref{lemma3.10_bound4}, \eqref{lemma3.11.bound2} and \deltascref{W9} into \eqref{3.128},  for any $\varepsilon \in (0,1]$, there holds
  \betagin{equation}
   \lambdabel{bound_3_1}
   \betagin{array}{ll}
    |\partial_t(\omega_i \vl)(z)| 

    & \apprle \frac{1}{\gamma r_i^2} \lbr \frac{r_i \lambda}{\varepsilon}  + \frac{\varepsilon}{r_i \lambda} |\tilde{u}h^i|^2 \rbr. 
   \end{array}
  \end{equation}

Making use of  \eqref{bound_3_1} in the expression for $J_3$ in \eqref{bound_J_3_1}, we get
  \betagin{equation*}
   \lambdabel{bound_I_3}
   \betagin{array}{ll}
    J_3 & \apprle \frac{1}{\gamma r_i^2} \lbr \frac{r_i \lambda}{\varepsilon}  + \frac{\varepsilon}{r_i \lambda} |\tilde{u}h^i|^2 \rbr\iint_{\frac34Q_i} |\tilde{u}-\tilde{u}^i|    \ dy \ d\tauu.

   \end{array}
  \end{equation*}
  We can use \eqref{bound_when_i_1} to get 

\betagin{equation}
 \lambdabel{bound_I33}
 J_3 \apprle \frac{\lambda^2}{r_i^2 \varepsilon} |Q_i| + \frac{\varepsilon}{\gamma r_i^2} |\tilde{u}^i|^2 |Q_i| \apprle \frac{\lambda^p }{\varepsilon} |4Q_i| + \varepsilon |4B_i| |\tilde{u}^i|^2.
\end{equation}

 \end{description}

 Substituting the estimates \eqref{bound_I_1}, \eqref{bound_I_22} and \eqref{bound_I33} into \eqref{first_1}, we obtain the proof of  \eqref{3.120}. 
 
 \item[Proof of \eqref{3.121}:] Note that in this case, we have $\gamma r_i^2 \apprge s$.  Setting $a =0$ in \eqref{3.123} along with making use of the the bounds  \eqref{3.126} and \eqref{3.127}, we get the following estimate:
 \betagin{equation}
  \lambdabel{second_1}
  \betagin{array}{ll}
   \left| \int_{8B} \lbr[(] \tilde{u}  \omega_i \vl \rbr (y,t) \ dy \right| & \apprle II_1 + II_2+II_3+II_4,

  \end{array}
 \end{equation}
where we have set 
\betagin{equation*}
 \betagin{array}{ll}
  II_1 & :=  \frac{1}{\min\{\rhoo, r_i\}} \iint_{8Q}\lbr  |\nabla u|^{p-1}+  |\theta|^{p-1}\rbr |\vl|   \lsb{\chii}{\frac34Q_i\cap 8Q} \ dy \ d\tauu, \\
  II_2 & :=  \iint_{8Q} \lbr |\nabla u|^{p-1} +|\theta|^{p-1}\rbr  |\nabla \vl|   \lsb{\chii}{\frac34Q_i\cap 8Q} \ dy \ d\tauu, \\
  II_3 & :=  \frac1{\min\{s,\gamma r_i^2\}} \iint_{8Q} |u| |\vl|  \lsb{\chii}{\frac34Q_i\cap 8Q} \ dy \ d\tauu,\\
  II_4 & :=   \iint_{8Q} |u| |\partial_t\vl|  \lsb{\chii}{\frac34Q_i\cap 8Q} \ dy \ d\tauu.
 \end{array}
\end{equation*}
We shall now estimate each of the terms as follows:
\betagin{description}[leftmargin=*]
 \item[Bound for $II_1$:] 
 In the case $\rhoo \leq r_i$, we can directly use Lemma \ref{lemma3.7} along with \deltascref{W4}, to get the estimate
 \betagin{equation}
  \lambdabel{3.83_1}

    II_1 \apprle \lambda \iint_{4Q_i} |\nabla u|^{p-1} +|\theta|^{p-1} \ dy \ d\tauu \apprle  |4Q_i| \lambda^p.

  \end{equation}
  
  In the case $r_i \leq \rhoo$, we make use of \eqref{lemma3.10_bound1}, \deltascref{W4} along with the bounds $|\tilde{u}| \leq |u|\lsb{\chii}{8Q}$ and $\gamma r_i^2 \apprge s$, to get 
  \betagin{equation}
  \lambdabel{3.84_1}
   \betagin{array}{ll}
    II_1 & \apprle \frac{1}{r_i} \lbr r_i \lambda + \frac{1}{\lambda r_i} \Yint\tiltlongdash_{4Q_i} |\tilde{u}|^2 \ dy \ d\tauu \rbr \lbr \iint_{4Q_i} |\nabla u|^{p-1}+|\theta|^{p-1} \ dy \ d\tauu \rbr \\

    & \apprle \lambda^p |4Q_i| + \frac{1}{s} \iint_{4Q_i} |u|^2 \lsb{\chii}{8Q} \ dy \ d\tauu.
   \end{array}
  \end{equation}

Thus combining \eqref{3.83_1} and \eqref{3.84_1}, we get 
\betagin{equation}
 \lambdabel{II_1_1}
 II_1 \apprle \lambda^p |4Q_i| + \frac{1}{s} \iint_{4Q_i} |u|^2\lsb{\chii}{8Q} \ dy \ d\tauu.
\end{equation}

 \item[Bound for $II_2$:] To bound this term, we make use of \eqref{lemma3.10_bound6}, \deltascref{W4} along with the observation $|\tilde{u}| \leq |u|\lsb{\chii}{8Q}$ and $\gamma r_i^2 \apprge s$, to get
 \betagin{equation}
  \lambdabel{II_2_1}
  \betagin{array}{ll}
   II_2 & \apprle \lbr \iint_{4Q_i} |\nabla u|^{p-1} +|\theta|^{p-1}\ dy \ d\tauu \rbr \lbr \lambda + \frac{\lambda^{1-p}}{s} \Yint\tiltlongdash_{4Q_i} |\tilde{u}|^2 \ dy \ d\tauu \rbr \\

           & \apprle \lambda^p |4Q_i| + \frac{1}{s} \iint_{4Q_i} |u|^2 \lsb{\chii}{8Q}\ dy \ d\tauu.
  \end{array}
 \end{equation}

 \item[Bound for $II_3$:] For bounding this term, we  use  \eqref{lemma3.10_bound1} along with $\gamma r_i^2 \apprge s$ and $|\tilde{u}|\leq|u| \lsb{\chii}{8Q}$,  to  get
 \betagin{equation}
  \lambdabel{II_3_1}
  \betagin{array}{ll}
   II_3 & \apprle \frac{1}{s} |4Q_i|\lbr \Yint\tiltlongdash_{4Q_i} |u| \lsb{\chii}{8Q} \ dy \ d\tauu\rbr\lbr  \Yint\tiltlongdash_{4Q_i} |\tilde{u}| \ dy \ d\tauu \rbr \\

   & \apprle \frac{1}{s} \iint_{4Q_i} |u|^2 \lsb{\chii}{8Q} \ dy \ d\tauu.
  \end{array}
 \end{equation}

 \item[Bound for $II_4$:] In this case, we make use of \eqref{lemma3.11.bound1} along with  $\gamma r_i^2 \apprge s$ and $|\tilde{u}|\leq|u| \lsb{\chii}{8Q}$, to  get
 \betagin{equation}
  \lambdabel{II_4_1}
  \betagin{array}{ll}
   II_4 & \apprle  |4Q_i|\lbr \Yint\tiltlongdash_{4Q_i} |u| \lsb{\chii}{8Q} \ dy \ d\tauu\rbr\lbr \frac{1}{\gamma r_i^2} \Yint\tiltlongdash_{4Q_i} |\tilde{u}| \ dy \ d\tauu \rbr \\

   & \apprle \frac{1}{s} \iint_{4Q_i} |u|^2 \lsb{\chii}{8Q} \ dy \ d\tauu.
  \end{array}
 \end{equation}

\end{description}
We finally combine the estimates \eqref{II_1_1}, \eqref{II_2_1}, \eqref{II_3_1} and \eqref{II_4_1} and recalling \eqref{second_1},  we obtain \eqref{3.121}.
\end{description}
This completes the proof of the lemma. 
\end{proof}

We now come to  essentially the most important estimate which will be used to obtain the Caccioppoli inequality in Section \ref{Caccioppoli_inequality}.
\betagin{lemma}
 \lambdabel{crucial_lemma}
 There exists a positive constant $C{(n,p,q, \Lambda_0,\Lambda_1,b_0)}$ such that the following estimate holds for every $t \in 2I$:
 \betagin{equation}
 \lambdabel{3122}
  \int_{8B \setminus \elam^t} (|\tilde{u}|^2 - |\tilde{u} - \vl|^2)(x,t) \ dx \geq C \lbr - \lambda^p |\mathbb{R}^{n+1} \setminus \elam|  - \frac{1}{s} \iint_{8Q} |\tilde{u}|^2 \ dz \rbr.
 \end{equation}

\end{lemma}
\betagin{proof}
 Let us fix any $t\in 2I$ and any point $x \in 8B \setminus \elam^t$.  Now define
 \betagin{equation}
  {\Upsilonsilon} := \left\{  i \in \Theta: \spt(\omega_i) \cap 8B \times \{t\} \neq \emptyset \quad \text{and} \quad |\tilde{u}| + |\vl| \neq 0 \quad  \text{on}\ \spt(\omega_i)\cap (8B \times \{t\}) \right\}.
 \end{equation}

 From \eqref{theta}, we see that if $\spt(\omega_i) \cap 8B \times \{t\} \neq \emptyset$, then  $i \in \Theta$.  If $i \neq {\Upsilonsilon}$, then $\tilde{u} = \vl = 0$ on $\spt(\omega_i) \cap 8B \times \{t\}$, which  implies
 \[
  \int_{\spt(\omega_i) \cap 8B \times\{t\}} |\tilde{u}|^2 - |\tilde{u} - \vl|^2 \ dx = 0.
 \]
 Hence we only need to consider $i \in {\Upsilonsilon}$.  We now decompose  ${\Upsilonsilon} = {\Upsilonsilon}_1 \cup {\Upsilonsilon}_2$, where ${\Upsilonsilon}_1 := {\Upsilonsilon} \cap \Theta_1$ and ${\Upsilonsilon}_2 := {\Upsilonsilon} \cap \Theta_2$. Noting that $\sum_{i \in {\Upsilonsilon}} \omega_i(\cdot, t) \equiv 1$ on $\mathbb{R}^n \cap \elam^t$, we can rewrite the left-hand side of \eqref{3122} as 
 \betagin{equation}
  \lambdabel{3124}
  \int_{8B \setminus \elam^t} (|\tilde{u}|^2 - |\tilde{u} - \vl|^2)(x,t) \ dx = J_1 + J_2,
 \end{equation}
where we have set
\betagin{gather*}
 J_1:= \sum_{i \in {\Upsilonsilon}_1} \int_{8B} \omega_i (|\tilde{u}|^2 - |\tilde{u} - \vl|^2) \ dx, \qquad 

 J_2:= \sum_{i \in {\Upsilonsilon}_2} \int_{8B} \omega_i (|\tilde{u}|^2 - |\tilde{u} - \vl|^2) \ dx 

\end{gather*}
We shall now estimate each of the terms as follows:
\betagin{description}[leftmargin=*]
 \item[Estimate of $J_1$:] We can further rewrite this term as follows:
 \betagin{equation}
  \lambdabel{I_1}
  \betagin{array}{ll}
  J_1 & = \sum_{i \in {\Upsilonsilon}_1} \int_{8B} \omega_i(z) \lbr |\tilde{u}^i|^2  + 2 \vl (\tilde{u} - \tilde{u}^i) \rbr \ dx - \sum_{i \in {\Upsilonsilon}_1} \int_{8B} \omega_i(z) |\vl - \tilde{u}^i|^2 \ dx\\

  & := J_1^1 + J_1^2.
  \end{array}
 \end{equation}

 \betagin{description}
  \item[Estimate of $J_1^1$:] Using \eqref{3.120}, we get
  \betagin{equation}
   \lambdabel{I_1_1}
   J_1^1 \apprge \sum_{i \in {\Upsilonsilon}_1} \int_{8B} \omega_i(z)  |\tilde{u}^i|^2 \ dz - \varepsilon \sum_{i \in {\Upsilonsilon}_1} |4B_i| |\tilde{u}^i|^2 - \sum_{i \in {\Upsilonsilon}_1} \frac{\lambda^p}{\varepsilon} |4Q_i|.
  \end{equation}
  From \eqref{def_tuh_i}, we have $\tilde{u}^i = 0$ whenever $\spt(\omega_i) \cap \Omega_{8\rhoo}^c \neq \emptyset$. Hence we only have to sum over all those $i \in {\Upsilonsilon}_1$ for which $\spt(\omega_i) \subset \mathcal{C}_{8\rhoo}$.  In this case, we make use of a suitable choice for $\varepsilon \in (0,1]$, and use \deltascref{W7} along with \deltascref{W8}, to estimate \eqref{I_1_1} from below. We have
  \betagin{equation}
 \lambdabel{bound_I1}
 J_1^1 \apprge -\lambda^p |\mathbb{R}^{n+1} \setminus \elam|. 
\end{equation}

  
  
  
  

  \item[Estimate of $J_1^2$:] For any $x \in 8B \setminus \elam^t$, we have from \deltascref{W10} that $\sum_{j} \omega_j(x,t) = 1$, which gives
  \betagin{equation}
   \lambdabel{I_2_1}
   \betagin{array}{ll}
    \omega_i(z) |\vl(z) - \tilde{u}^i|^2  
    & \apprle \omega_i(z)  \sum_{j \in A_i} |\omega_j(z)|^2  \lbr \tilde{u}^j - \tilde{u}^i\rbr^2 \\
    & \overset{\redlabel{3.104.a}{a}}{\apprle} \min\{ \rhoo, r_i\}^2 \lambda^2.
   \end{array}
  \end{equation}
To obtain \redref{3.104.a}{a} above, we made use of Lemma \ref{lemma3.8} (recall $i \in {\Upsilonsilon}_1 \subset \Theta_1$) along with \deltascref{W13}.  Substituting \eqref{I_2_1} into the expression for $J_1^2$ and using $|Q_i| = |B_i| \times 2\gamma r_i^2$, we get
\betagin{equation}
 \lambdabel{bound_I_2}

J_1^2 
\apprle 

 \sum_{i \in {\Upsilonsilon}_1} |8B \cap \frac34B_i| \frac{\gamma r_i^2}{\gamma} \lambda^2 

 \apprle \lambda^p |\mathbb{R}^{n+1} \setminus \elam|.

\end{equation}
 \end{description}
Substituting \eqref{bound_I1} and \eqref{bound_I_2} into \eqref{I_1}, we get
\betagin{equation}
 \lambdabel{bound_I}
 J_1 \apprge - \lambda^p |\mathbb{R}^{n+1} \setminus \elam|.
\end{equation}

 \item[Estimate of $J_2$:] From the identity $|\tilde{u}|^2 - |\tilde{u}-\vl|^2 = 2\tilde{u} \vl - |\vl|^2$, we get
 \betagin{equation}
  \lambdabel{II1}
  \betagin{array}{ll}
  \sum_{i \in {\Upsilonsilon}_2} \left| \int_{8B} \omega_i (|\tilde{u}|^2 - |\tilde{u} - \vl|^2) \ dx\right| & \apprle \sum_{i \in {\Upsilonsilon}_2} \left| \int_{8B} \omega_i \tilde{u} \vl \ dx\right| + \sum_{i \in {\Upsilonsilon}_2} \left| \int_{8B} \omega_i |\vl|^2 \ dx\right| \\
  & := J_2^1 + J_2^2.
  \end{array}
 \end{equation}
\betagin{description}
 \item[Estimate for $J_2^1$:]  This term can be easily estimated using \eqref{3.121} followed by summing over $i \in {\Upsilonsilon}_2$ and using \deltascref{W13}, we get
 \betagin{equation}
  \lambdabel{II1_1}
  \betagin{array}{ll}
   \sum_{i \in {\Upsilonsilon}_2} \left| \int_{8B} \omega_i \tilde{u} \vl \ dx\right| & \apprle \sum_{i \in {\Upsilonsilon}_2}\lbr  \lambda^p |Q_i| + \frac{1}{s} \iint_{4Q_i} |u|^2 \lsb{\chii}{8Q}\ dz\rbr \\
   & \apprle \lambda^p |\mathbb{R}^{n+1} \setminus \elam| + \frac{1}{s} \iint_{8Q} |u|^2 \ dz.

  \end{array}
 \end{equation}

 \item[Estimate for $J_2^2$:] To estimate this term, we make use of \eqref{lemma3.10_bound1} along with $|Q_i| = |B_i| \times 2\gamma r_i^2$, $\gamma r_i^2 \apprge s$, the bound $|\tilde{u}| \leq |u| \lsb{\chii}{8Q}$ and \deltascref{W13}, to get
 \betagin{equation}
  \lambdabel{II2_1}
  \betagin{array}{ll}
   \sum_{i \in {\Upsilonsilon}_2} \left| \int_{8B} \omega_i |\vl|^2 \ dx\right| & \apprle \sum_{i \in {\Upsilonsilon}_2}  \lbr \Yint\tiltlongdash_{4Q_i} |\tilde{u}| \ dz \rbr^2 |8B \cap \frac34B_i| 

    \apprle \frac{1}{s} \iint_{8Q} |u|^2 \ dz.
  \end{array}
 \end{equation}

\end{description}
 We combine \eqref{II1_1}, \eqref{II2_1} and \eqref{II1} to obtain
 \betagin{equation}
  \lambdabel{bound_II}

  J_2\apprle \lambda^p |\mathbb{R}^{n+1} \setminus \elam| + \frac{1}{s} \iint_{8Q} |u|^2 \ dz.

 \end{equation}
\end{description}
Thus, from \eqref{bound_I}, \eqref{bound_II} and\eqref{3124}, the proof of the Lemma follows.

\end{proof}

\section{Caccioppoli inequality and Reverse H\"older inequality}\lambdabel{Caccioppoli_inequality}
We now have all the estimates need to prove a Caccioppoli type inequality near the boundary. \emph{In subsequent calculations, we shall alternate between $\mathcal{C}_{8\rhoo} =\Omega_{8\rhoo} \times 8I$ and $8Q$ depending on convenience, as $\vlh(z) = \tilde{u}(z) =u(z)= 0$ for all $z \in \mathcal{C}_{8\rhoo}^c$.}

Before we state Lemma \ref{Caccioppoli}, let us make the following remark:
\betagin{remark}
 \lambdabel{def_al_0_1}
 Under the hypothesis of Lemma \ref{Caccioppoli} (more specifically suppose \eqref{hypothesis_1} holds), then we see that there exists a universal constant $c_e = c_e(n)$ such that for any $\lambda \geq c_e \alpha_0$, there holds $\elam \neq \emptyset$. This follows from Lemma \ref{bound_g_x_t} and \eqref{hypothesis_1} below. 
\end{remark}

\subsection{Caccioppoli Inequality}
We shall prove the Caccioppoli inequality in this subsection.
\betagin{lemma}
 \lambdabel{Caccioppoli}
 Let $\kappa \geq 1$ and $\Omega$ satisfy Definition \ref{p_thick_domain} with constants $r_0,b_0$. Then there exists  constants $\beta_0 = \beta_0(n,p,q,\Lambda_0,\Lambda_1,b_0,\kappa) \in (0,1)$ small and $C = C(n,p,q,\Lambda_0,\Lambda_1,\kappa,b_0)$ such that the following holds.

Let $Q := B \times I = Q_{\rhoo,s}(z)$ be a parabolic cylinder  with $s = \alpha_0^{2-p} \rhoo^2$   for some $\alpha_0>0$ as given in \eqref{def_s} and let $u \in L^2(-T,T;L^2(\Omega)) \cap L^{p-\beta}(-T,T; W_0^{1,p-\beta}(\Omega))$ for some $\beta \in (0,\beta_0)$ be any very weak solution of \eqref{main} in the sense of \eqref{very_weak_solution}.  Suppose the following bound holds:
 \betagin{gather}
  \frac{\alpha_0^{p-\beta}}{\kappa} \leq \Yint\tiltlongdash_Q (|\nabla u|+|\theta|)^{p-\beta} \ dz \quad \text{and} \quad \Yint\tiltlongdash_{8Q} (|\nabla u|+|\theta|)^{p-\beta} \ dz \leq \kappa \alpha^{p-\beta}_0 \lambdabel{hypothesis_1},
 \end{gather}
 where $h$ is from \eqref{bound_b}. Then we have
 \betagin{equation*}
  \lambdabel{conclusion_1}
  \alpha_0^{p-\beta} + \alpha_0^{p-2} \sup_{t \inI} \Xint\tilthlongdash_B \left| \frac{u(x,t)}{\rhoo} \right|^2 \mathcal{M}(x,t)^{-\beta} \ dx \leq {C}\Yint\tiltlongdash_{8Q}\lbr  \alpha_0^{p-2-\beta} \left| \frac{u}{\rhoo} \right|^2 + \left| \frac{u}{\rhoo} \right|^{p-\beta} + |\theta|^{p-\beta} \rbr dz
 \end{equation*}
 where we have set $\mathcal{M}(x,t):=\max\{g(x,t), \alpha_0\}$.

\end{lemma}

\betagin{proof}[Proof of Lemma \ref{Caccioppoli}]
Let $(t_1 ,t) \subset 2I$ be a time interval and recall the definition of $\vlh$ from \eqref{lipschitz_extension}.  We  shall make use of  $\eta(x)\vlh(x,\tauu)$ as a test function in \eqref{very_weak_solution} (this is possible since $\spt(\vlh) \subset \Omega_{8\rhoo} \times 8I$ and $\vlh C^{0,1}(2\mathcal{H})$ from Lemma \ref{lemma3.15}). Thus, after integrating the resulting expression of \eqref{very_weak_solution} over $(t_1,t)$, we get

 \betagin{equation}
 \lambdabel{cac1}
  L_1 + L_2:=\int_{t_1}^t \left[ \int_{\Omega_{8\rhoo}} \frac{d{[u]_h}}{d\tauu} \eta(y) \vlh(y,\tauu) + \iprod{[\mathcal{A}(y,\tauu,\nabla u)]_h}{\nabla (\eta \vlh)} \ dy\right]\ d\tauu = 0.
 \end{equation}

 
 
 
 \betagin{description}[leftmargin=*]
  \item[Estimate of $L_1$:] Note that $\zetata(\tauu) = 1$ for all $\tauu \in (t_1,t)$. Thus making use of the standard hole filling technique and  \eqref{def_u_tilde}, we get
  \betagin{equation}
   \lambdabel{6.38}
   \betagin{array}{ll}
    L_1 
    & = \int_{t_1}^t \int_{\Omega_{8\rhoo}} \dds{\tilde{u}_h(y,\tauu)}  \vlh(y,s) \ dy\ d\tauu \\

    & = \int_{t_1}^t \int_{\Omega_{8\rhoo}\setminus \elam^{\tauu}}  \dds{\vlh}  (\vlh-\tilde{u}_h) \ dy \ d\tauu +    \int_{t_1}^t \int_{\Omega_{8\rhoo}}  \frac{d{\lbr (\tilde{u}_h)^2 - (\vlh - \tilde{u}_h)^2 \rbr }}{d\tauu}  \ dy \ d\tauu \\
    & := J_2 + J_1(t) - J_1(t_1),
   \end{array}
  \end{equation}
  where we have set 
  \betagin{equation*}
  \lambdabel{def_i_1}J_1(\tauu) := \frac12 \int_{\Omega_{8\rhoo}} ( (\tilde{u}_h)^2 - (\vlh - \tilde{u}_h)^2 ) (y,\tauu) \ dy.
  \end{equation*}

  \betagin{description}
   \item[Estimate for $J_2$:] Taking absolute values and making use of Lemma \ref{lemma3.14} with $\vartheta =1$, we get
   \betagin{equation}
    \lambdabel{6.39}
    \betagin{array}{ll}
     |J_2| 

     & \apprle \iint_{8Q\setminus \elam}   \left| \dds{\vlh}  (\vlh-\tilde{u}_h)\right| \ dy \ d\tauu 
       \apprle\lambda^p |\mathbb{R}^{n+1} \setminus \elam| + \frac{1}{s} \iint_{8Q} |[u]_h|^2 \ dy \ d\tauu. 
    \end{array}
   \end{equation}

   \item[Estimate for $J_1(t_1)$:] We first claim that we can choose $t_1\in (-2^2s,-s)$ such that
   \betagin{equation}
    \lambdabel{6.40}
    |J_1(t_1)|  

    \leq \frac{1}{s} \int_{-2^2s}^{-s} \left| \int_{\Omega_{8\rhoo}} \lbr |\tilde{u}_h|^2 - |\vlh - \tilde{u}_h|^2 (x,\tauu)\rbr \ dy\right| d\tauu.
   \end{equation}
Suppose not, then for any $t_1 \in (-2^2s,-s)$, we would have 
\betagin{equation}
\lambdabel{J_1_t_11}
 \left| \int_{\Omega_{8\rhoo}} \lbr |\tilde{u}_h|^2 - |\vlh - \tilde{u}_h|^2\rbr (y,t_1) \ dy \right| > \frac{1}{s} \int_{-2^2s}^{-s} \left| \int_{\Omega_{8\rhoo}} \lbr |\tilde{u}_h|^2 - |\vlh - \tilde{u}_h|^2 \rbr(y,\tauu) \ dy\right| d\tauu.
\end{equation}
But we have
\betagin{equation}
\lambdabel{J_1_t_22}
 3 \Xint\tilthlongdash_{-2^2s}^{-s} \left| \int_{\Omega_{8\rhoo}} \lbr |\tilde{u}_h|^2 - |\vlh - \tilde{u}_h|^2 (y,\tauu)\rbr \ dy\right| d\tauu \geq 3 \min_{\tauu \in (-4s,-s)} |J_1(\tauu)|.
\end{equation}
Combining \eqref{J_1_t_11} and \eqref{J_1_t_22}, for any $t_1 \in (-4s,-s)$, we get $|J_1(t_1)| > 3 \min_{\tauu \in (-4s,-s)} |J_1(\tauu)|$,
which is absurd. Hence \eqref{6.40}  must be true for some $t_1 \in (-4s,-s)$. 	 

From the construction \eqref{lipschitz_extension},  we have $\vlh = \tilde{u}_h$ on $\elam$. Furthermore,  $\spt(\tilde{u}_h) \subset 8Q$ and $|\tilde{u}_h| \leq |[u]_h|\lsb{\chii}{8Q}$ holds. For $t_1 \in (-2^2s,-s)$ satisfying \eqref{6.40}, we then get
\betagin{equation}
 \lambdabel{6.44}
 \betagin{array}{ll}
  |J_1(t_1)| 
& \leq   \frac{1}{s} \iint_{8Q\cap \elam } \left| |\tilde{u}_h|^2 - |\vlh - \tilde{u}_h|^2 \right| dz + \frac{1}{s} \iint_{8Q\setminus \elam} \left| |\tilde{u}_h|^2 - |\vlh - \tilde{u}_h|^2 \right| dz \\

& \apprle   \frac{1}{s} \iint_{8Q\cap \elam }  |\tilde{u}_h|^2  dz + \frac{1}{s} \iint_{8Q\setminus \elam}  |\tilde{u}_h|^2 + |\vlh |^2  dz \\

& \overset{\redlabel{4.8.a}{a}}{\apprle} \frac{1}{s} \iint_{8Q} |[u]_h|^2 \ dz.
 \end{array}
\end{equation}
To obtain \redref{4.8.a}{a}, we used Lemma \ref{lemma3.12} (applied with $\vartheta = 2$).

  \end{description}

  \item[Estimate for $L_2$:] 
   We decompose the expression as 
   \betagin{equation}\lambdabel{4.9}
    \betagin{array}{ll}
     L_2 & = \int_{t_1}^t \int_{\elam^{\tauu}} \iprod{[\mathcal{A}(y,\tauu,\nabla u)]_h}{\nabla (\eta \vlh)} \ dy\ d\tauu  + \\
     & \qquad \quad + \int_{t_1}^t \int_{\Omega_{8\rhoo} \setminus \elam^{\tauu}} \iprod{[\mathcal{A}(y,\tauu,\nabla u)]_h}{\nabla (\eta \vlh)} \ dy\ d\tauu \\
     & := L_2^1 + L_2^2.
    \end{array}
   \end{equation}

   \betagin{description}
    \item[Estimate for $L_2^2$:] Using the chain rule, \eqref{abounded}, \eqref{cut_off_function} along with Corollary \ref{corollary3.20}, we get
    \betagin{equation}
    \lambdabel{4.10}
     \betagin{array}{ll}
      L^2_2 
      & \leq \int_{t_1}^t \int_{\Omega_{8\rhoo} \setminus E_{\tauu}(\lambda)} {[|\nabla u|^{p-1}+|\theta|^{p-1}]_h}{|\nabla (\eta \vlh)|} \ dy\ d\tauu  \\
      & \apprle \sum_{k=0}^1  \rhoo^{-k+1}\iint_{(8B\times 2I) \setminus \elam} {[|\nabla u|^{p-1}+|\theta|^{p-1}]_h}{|\nabla^k  \vlh|} \ dy\ d\tauu  \\%+ \frac{1}{\rhoo}\int_{(8B\times 2I) \setminus \elam} {[|\nabla u|^{p-1}]_h}{|\vlh|} \ dy\ d\tauu  \\
&\apprle \lambda^p |\mathbb{R}^{n+1} \setminus \elam| + \frac{1}{s} \iint_{8Q} |\tilde{u}h|^2 \ dy\ d\tauu.
     \end{array}
    \end{equation}

\end{description} 
\end{description}

 
Substituting \eqref{4.10} into \eqref{4.9} and \eqref{6.39}, \eqref{6.44} into \eqref{6.38},  and finally making use of \eqref{cac1} along with the bound $|\tilde{u}h| \leq |[u]_h|\lsb{\chii}{8Q}$, we get

 \betagin{equation}
 \lambdabel{6.45}
  \betagin{array}{l}
    \frac12 \int_{\Omega_{8\rhoo}} | (\tilde{u}_h)^2 - (\vlh - \tilde{u}_h)^2 | (y,t) \ dy +  \int_{t_1}^t \int_{E_{\tauu}(\lambda)} \iprod{[\mathcal{A}(y,\tauu,\nabla u)]_h}{\nabla (\eta \vlh)} \ dy \ d\tauu   \\

   \hspace*{6cm} \apprle \lambda^p |\mathbb{R}^{n+1} \setminus \elam| + \frac{1}{s} \iint_{8Q} |[u]_h|^2 \ dy \ d\tauu.
  \end{array}
 \end{equation}

 Since the estimate in \eqref{6.45} is independent of $h$, we can let $h \searrow 0$ to get
 \betagin{equation}
 \lambdabel{6.46.1}
  \betagin{array}{l}
    \frac12 \int_{\Omega_{8\rhoo}} | (\tilde{u})^2 - (\vl - \tilde{u})^2 | (y,t) \ dy +  \int_{t_1}^t \int_{E_{\tauu}(\lambda)} \iprod{\mathcal{A}(y,\tauu,\nabla u)}{\nabla (\eta \vl)} \ dy \ d\tauu   \\

   \hspace*{6cm} \apprle \lambda^p |\mathbb{R}^{n+1} \setminus \elam| + \frac{1}{s} \iint_{8Q} |u|^2 \ dy \ d\tauu.
  \end{array}
 \end{equation}

  Using the fact that $\tilde{u} = \vl$ on $\elam$ and Lemma \ref{crucial_lemma},  we get
 \betagin{equation}
 \lambdabel{4.13}
  \betagin{array}{ll}
    \int_{\Omega_{8\rhoo}} | (\tilde{u})^2 - (\vl - \tilde{u})^2 | (y,t) \ dy   

   & \apprge \int_{\elam^t} | \tilde{u} (x,t)|^2 \ dx   - \lambda^p |\mathbb{R}^{n+1} \setminus \elam| - \frac{1}{s} \iint_{8Q} |u|^2 \ dz. 
  \end{array}
 \end{equation}

Making use of the bound \eqref{4.13} into \eqref{6.46.1}, followed by multiplying the resulting expression with $\lambda^{-1-\beta}$ and integrating from $(c_e\alpha_0,\infty)$ (recall that  $c_e$ is as in Remark \ref{def_al_0_1}), for almost every $t \in (-s,s)$, we get 
 \betagin{equation}
 \lambdabel{K_expression}
K_1 + K_2 \apprle K_3 + K_4,
 \end{equation}
 where we have set
 \betagin{equation*}
  \betagin{array}{@{}r@{}c@{}l@{}}
   K_1 \ &:=&\  \frac12\int_{c_e \alpha_0}^{\infty} \lambda^{-1-\beta} \int_{\elam^t} | \tilde{u}(y,t)|^2 \ dy \ d\lambda, \\%\lambdabel{def_K_1}\\
  K_2 \ &:=& \ \int_{c_e \alpha_0}^{\infty} \lambda^{-1-\beta}\int_{t_1}^t \int_{\elam^{\tauu}} \iprod{\mathcal{A}(y,\tauu,\nabla u)}{\nabla (\eta u)} \ dy\ d\tauu \ d\lambda, \\%\lambdabel{def_K_2}\\
  K_3\  &:=& \ \int_{c_e \alpha_0}^{\infty} \lambda^{-1-\beta}  \lambda^p |\mathbb{R}^{n+1} \setminus \elam| \  d\lambda, \\%\lambdabel{def_K_3}\\
  K_4 \ &:=& \ \frac{1}{s} \int_{c_e \alpha_0}^{\infty} \lambda^{-1-\beta}   \iint_{8Q} |u|^2 \ dy  \ d\tauu \ d\lambda .
  \end{array}
 \end{equation*}

 We now define the truncated Maximal function $\mathcal{M}(z) := \max \{ g(z), \alpha_0\}$ and then estimate each of the $K_i$ for $i \in \{1,2,3,4\}$ as follows:
 \betagin{description}
  \item[Estimate for $K_1$:]  By applying Fubini, we get
  \betagin{equation}
  \lambdabel{4.19}
   \betagin{array}{ll}
    K_1 & \geq \int_{8B} |\tilde{u}(y,t)|^2 \int_{c_e g(y,t)}^{\infty} \lambda^{-1-\beta} \ d\lambda \ dy 
     \apprge \frac{1}{\beta c_e^{\beta}}  \int_{8B} \mathcal{M}(y,t)^{-\beta} | \tilde{u}(y,t)|^2 \ dy .
   \end{array}
  \end{equation}

  \item[Estimate for $K_2$:] Again applying Fubini, we get
  \betagin{equation}
  \lambdabel{4.20}
   \betagin{array}{ll}
    K_2 & = \frac{1}{\beta c_e^{\beta}} \int_{t_1}^s \int_{\Omega_{8\rhoo}} \mathcal{M}(y,\tauu)^{-\beta} \iprod{\mathcal{A}(y,\tauu,\nabla u)}{\nabla (\eta u)} \ dy \ d\tauu. 
   \end{array}
  \end{equation}
  Applying the chain rule and using \eqref{def_u_tilde} along with \eqref{abounded} and \eqref{bound_b}, we  get
  \betagin{equation}
\lambdabel{4.25}
 \betagin{array}{ll}

  K_2 &= \int_{t_1}^s \int_{\Omega_{8\rhoo}} \mathcal{M}(y,\tauu)^{-\beta} \iprod{\mathcal{A}(y,\tauu,\nabla u)}{\nabla u} \eta^2  \ dy \ d\tauu   \\
  & \qquad + \int_{t_1}^s \int_{\Omega_{8\rhoo}} \mathcal{M}(y,\tauu)^{-\beta} \iprod{\mathcal{A}(y,\tauu,\nabla u)}{\nabla \eta^2 } u \ dy \ d\tauu  \\
  &\apprge \int_{t_1}^s \int_{\Omega_{8\rhoo}} \mathcal{M}(y,\tauu)^{-\beta} |\nabla u|^p \eta^2  \ dy \ d\tauu  - \iint_{8Q} \mathcal{M}(y,\tauu)^{-\beta} |\theta|^p   \ dy \ d\tauu  \\
  & \qquad - \int_{t_1}^s \int_{\Omega_{8\rhoo}} \mathcal{M}(y,\tauu)^{-\beta} \lbr |\nabla u|^{p-1} + |\theta|^{p-1} \rbr \frac{|u|}{\rhoo} \ dy \ d\tauu  \\
  & \apprge  \iint_{Q} \mathcal{M}(y,\tauu)^{-\beta} |\nabla u|^p \eta^2  \ dy \ d\tauu  - \iint_{8Q} \mathcal{M}(y,\tauu)^{-\beta} |\theta|^p   \ dy \ d\tauu \\
  & \qquad - \iint_{8Q} \mathcal{M}(y,\tauu)^{-\beta} \lbr |\nabla u|^{p-1} + |\theta|^{p-1} \rbr \frac{|u|}{\rhoo} \ dy \ d\tauu  \\
  & := \mathcal{A}_1 + \mathcal{A}_2 + \mathcal{A}_3.
 \end{array}
\end{equation}

\betagin{description}
 \item[Estimate for $\mathcal{A}_1$:] Note that $\eta \equiv 1$ on $B$.  Let $S :=\{ z \in Q: |\nabla u(z)| \geq \beta g(z)\},$ then we get
 \betagin{equation}
 \lambdabel{4.26}
  \betagin{array}{@{}l@{}c@{}l@{}}
   \iint_Q |\nabla u|^{p-\beta} \ dz  & =& \iint_S |\nabla u|^{p-\beta} \ dz + \iint_{Q\setminus S} |\nabla u|^{p-\beta} \ dz \\

   & \leq &\beta^{-\beta} \iint_Q \mathcal{M}(z)^{-\beta} |\nabla u|^p \ dz + \beta^{p-\beta} \iint_{Q\setminus S} \mathcal{M}(z)^{p-\beta} \ dz \\

            & \overset{\text{Lemma \ref{bound_g_x_t}}}{\apprle} & \iint_Q \mathcal{M}(z)^{-\beta} |\nabla u|^p \ dz + \beta^{p-\beta} \iint_{8Q} |\nabla u|^{p-\beta} \ dz  + \beta^{p-\beta} |Q| \alpha_0^{p-\beta}.
      \end{array}
 \end{equation}

  \item[Estimate for $\mathcal{A}_2$:]  From \eqref{def_g}, we obtain the bound  $\lsb{\chii}{8Q}( |\nabla u(z)| + |\theta(z)|)  \leq \mathcal{M}(z)$ for a.e $z \in \mathbb{R}^n$. This gives
  \betagin{equation}
   \lambdabel{4.26.1}
 \mathcal{A}_2 =  \iint_{8Q} \mathcal{M}(z)^{-\beta} |\theta|^p   \ dz \apprle  \iint_{8Q}  |\theta|^{p-\beta}   \ dz.
  \end{equation}

 
 \item[Estimate for $\mathcal{A}_3$:]  We use the bound  $\lsb{\chii}{8Q}( |\nabla u(z)| + |\theta(z)|)  \leq \mathcal{M}(z)$ for a.e $z \in \mathbb{R}^n$, along with Young's inequality, Theorem \ref{sobolev-poincare} and  \eqref{hypothesis_1}, to get
 \betagin{equation}
 \lambdabel{4.27}
  \betagin{array}{rcl}
   \mathcal{A}_3 & \apprle& \iint_{8Q}  (|\nabla u|+ |\theta|) ^{p-1-\beta} \frac{|u|}{\rhoo} \ dz

   \apprle  \varepsilon |Q| \kappa \alpha_0^{p-\beta} + C(\varepsilon) \iint_{8Q} \left| \frac{u}{\rhoo}\right|^{p-\beta}\ dz.
  \end{array}
 \end{equation}

\end{description}

  \item[Estimate for $K_3$:] Applying the layer-cake representation (see for example \cite[Chapter 1]{Grafakos}), we get
  \betagin{equation}
  \lambdabel{4.21}
   \betagin{array}{rcl}
    K_3  & =& \frac{1}{p-\beta} \iint_{\mathbb{R}^{n+1}} \mathcal{M}(z)^{p-\beta} \ dz

     \overset{\text{Lemma \ref{bound_g_x_t}}}{\apprle} \frac{1}{p-\beta} \iint_{8Q}  \lbr |\nabla u|^{p-\beta}+ |\theta|^{p-\beta}+\alpha_0^{p-\beta}\rbr  \ dz \\
    & \overset{\text{\eqref{hypothesis_1}}}{\apprle} &\alpha_0^{p-\beta} |Q|.  
   \end{array}
  \end{equation}

  \item[Estimate for $K_4$] Again applying Fubini, we get 
  \betagin{equation}
  \lambdabel{4.22}
   \betagin{array}{ll}
    K_4 & =\frac{1}{s} \int_{{c_e \alpha_0}}^{\infty} \lambda^{-1-\beta}   \iint_{8Q} |u|^2 \ dz \ d\lambda  = \frac{1}{\beta} \iint_{8Q}  (\alpha_0)^{-\beta}  \frac{|u|^2}{s} \ dz.
   \end{array}
  \end{equation}

 \end{description}

Substituting \eqref{4.26}, \eqref{4.26.1} and \eqref{4.27} into \eqref{4.25} followed by combining  \eqref{4.19}, \eqref{4.21} and \eqref{4.22} into \eqref{K_expression}, we get
\betagin{equation*}
\lambdabel{4.28}
 \betagin{array}{ll}
 \frac{1}{2\beta}  \int_{B} \mathcal{M}(y,t)^{-\beta} | \tilde{u}(y,t)|^2 \ dy + \frac{1}{\beta} \iint_Q |\nabla u|^{p-\beta} \ dz \apprle \frac{1}{\beta} \beta^{p-\beta} \iint_{8Q} |\nabla u|^{p-\beta} \ dz +\iint_{8Q} |\theta|^{p-\beta} \ dz\\
    \hspace*{2cm}+ \frac{1}{\beta} \varepsilon |Q| \alpha_0^{p-\beta} + \frac{1}{\beta} C(\varepsilon) \iint_{8Q} \left| \frac{u}{\rhoo}\right|^{p-\beta}\ dz+  \alpha_0^{p-\beta} |Q| + \frac{1}{\beta} \iint_{8Q} \alpha_0^{-\beta}  \frac{|u|^2}{s} \ dz.

 \end{array}
\end{equation*}

Multiplying the above expression  by $\beta$ and then using the intrinsic scaling $s = \rhoo^2 \alpha_0^{2-p}$ from \eqref{def_s} along with \eqref{hypothesis_1}, we get 
\betagin{equation*}
 \betagin{array}{ll}
        \int_{B} \mathcal{M}(y,t)^{-\beta} | \tilde{u}(y,t)|^2 \ dy + |Q| \alpha_0^{p-\beta}  \apprle \beta^{p-\beta} |Q| \alpha_0^{p-\beta} + \varepsilon |Q| \alpha_0^{p-\beta} +\iint_{8Q} |\theta|^{p-\beta} \ dz \\
    \hspace*{2cm}  +  C(\varepsilon) \iint_{8Q} \lbr \frac{|u|}{\rhoo} \rbr^{p-\beta} \ dz + \beta \alpha_0^{p-\beta} |Q| +  \iint_{8Q} \alpha_0^{p-2-\beta}   \lbr \frac{|u|}{\rhoo}\rbr^2 \ dz.
 \end{array}
\end{equation*}

Choosing $\beta \in (0,\beta_0)$  and $\varepsilon \in (0,1)$ small, we then make use of  \eqref{hypothesis_1} and the observation $|Q| \approx |B| s^2 = \rhoo^2 \alpha_0^{2-p}|B|$, to get
\betagin{equation*}
\lambdabel{4.31}

        \sup_{t \in I} \alpha_0^{p-2} \Xint\tilthlongdash_{B} \mathcal{M}(y,t)^{-\beta} \left| \frac{ \tilde{u}(y,t)}{\rhoo}\right|^2 \ dy +  \alpha_0^{p-\beta}   \ \apprle \   \Yint\tiltlongdash_{8Q}\left[ \alpha_0^{p-2-\beta}   \lbr \frac{|u|}{\rhoo}\rbr^2 +  \lbr \frac{|u|}{\rhoo} \rbr^{p-\beta} + |\theta|^{p-\beta} \right]\ dz.

\end{equation*} 
This completes the proof of the Lemma. 
\end{proof}

\subsection{Estimates needed to prove Reverse H\"older inequality}

\betagin{lemma}
\lambdabel{lemma7.4}
Let $\kappa \geq 1$, then there exists $\beta_0(n,p,q,\Lambda_0,\Lambda_1,b_0,\kappa)$. Let  $\beta \in (0,\beta_0)$ be given and consider any very weak solution $u \in L^2(-T,T;L^2(\Omega)) \cap L^{p-\beta}(-T,T; W_0^{1,p-\beta}(\Omega))$. Let $Q_{\rhoo,s}(0,0) = B_{\rhoo}(0) \times I_s(0)$ be the parabolic cylinder with $(0,0) \in \partial \Omega \times (-T,T)$ and $s = \rhoo^2 \alpha_0^{2-p}$ for some $\alpha_0 >0$ as in \eqref{def_s}. Let $\alpha Q$ be a rescaled parabolic cylinder for some $\alpha \in (1,8]$ and let $p-\varepsilon_0$  be as given in \eqref{def_q} and also suppose that 
 \betagin{equation}\lambdabel{hypothesis_2}
  \Yint\tiltlongdash_{ \alpha Q} \lbr |\nabla u|+|\theta| \rbr^{p-\beta} \ dz \leq \kappa \alpha_0^{p-\beta}.
 \end{equation}
 Let us define
\betagin{equation}
 \lambdabel{def_J}
 J:= \sup_{t\in I} \Xint\tilthlongdash_{ B} \lbr \frac{|u|}{\rhoo}\rbr^2 \tilde{g}(x,t)^{-\beta} \ dx,
\end{equation}
where $\tilde{g}(z) := \max\{ \mathcal{M}( |\nabla \tilde{u}|^q \lsb{\chii}{\alpha Q})^{\frac{1}{q}}(z), \mathcal{M}(( |\nabla u|^q+|\theta|^q )\lsb{\chii}{\alpha Q})^{\frac{1}{q}}(z), \alpha_0\}$.
 Then there holds:

\betagin{equation}
\lambdabel{lemma4.3.1}
 \Yint\tiltlongdash_{ Q} \lbr \frac{|u|}{\rhoo} \rbr^{p-\beta} \ dz \leq C_{(n,p,q,\Lambda_0,\Lambda_1,b_0,\kappa)} (\alpha_0^{\beta} J)^{\frac{\varepsilon_0-\beta}{2}}  \lbr \Yint\tiltlongdash_{ Q} |\nabla u|^q \ dz\rbr^{\frac{p-\varepsilon_0}{q}}  .
\end{equation}
 If $\frac{p-\varepsilon_0}{2} \geq 1$, then there exists an $\tilde{\gamma}a = \tilde{\gamma}a(n,p,\varepsilon_0) \in (0,1)$ such that the following holds:
 \betagin{equation}
  \lambdabel{lemma4.3.2}
 \betagin{array}{ll}
\Yint\tiltlongdash_{ Q} \lbr \frac{|u|}{\rhoo} \rbr^{2} \ dz & \leq C_{(n,p,q,\Lambda_0,\Lambda_1,b_0,\kappa)} (\alpha_0^{\beta} J)^{1-\tilde{\gamma}a}\lbr  \Yint\tiltlongdash_{ Q} |\nabla u|^{q} \ dz\rbr^{\frac{2\tilde{\gamma}a}{q}} .
\end{array}
\end{equation}
 If $\frac{p-\varepsilon_0}{2} \in (0,1)$, then there holds
\betagin{equation}

\lambdabel{lemma4.3.3}
 \betagin{array}{ll}
\Yint\tiltlongdash_{ Q} \lbr \frac{|u|}{\rhoo} \rbr^{2} \ dz & \leq C_{(n,p,q,\Lambda_0,\Lambda_1,b_0,\kappa)} (\lambda_0^{\beta} J)^{\frac{2-p+\varepsilon_0}{2}}\lbr  \Yint\tiltlongdash_{ Q} |\nabla u|^{q} \ dz\rbr^{\frac{p-\varepsilon_0}{q}} .
\end{array}
\end{equation}

\end{lemma}

\betagin{proof}
 In order to prove the lemma, we want to make use of Lemma \ref{interpolation_capacity_version_main}. Let us prove each of the assertions as follows:
 

 

\betagin{description}[leftmargin=*]
 \item[Proof of \eqref{lemma4.3.1}:] In this case, in  order to apply Lemma \ref{interpolation_capacity_version_main}, we make the following choice of exponents:
 \betagin{equation}
 \lambdabel{exp_ch}
\sigma =  p-\beta, \quad  \tilde{r} = \frac{2(p-\beta)}{p}, \quad \vartheta = p-\varepsilon_0 \qquad \text{and} \qquad \tilde{\alphapha} = \frac{p-\varepsilon_0}{p-\beta} < 1. 
 \end{equation}
 These choice of exponents are admissible under the additional restriction that $\varepsilon_0 \leq \max \left\{\frac{4}{n+2}, \frac14\right\}$ and  $\beta_0$ is small such that $\frac{(p-\varepsilon_0)^2}{p} \geq \frac{n}{2} (\varepsilon_0-\beta_0)$ holds. Thus applying Lemma \ref{interpolation_capacity_version_main}, we get
 \betagin{equation}
 \lambdabel{7.14}
\Yint\tiltlongdash_Q \lbr \frac{|u|}{\rhoo} \rbr^{p-\beta} \ dz \apprle \Xint\tilthlongdash_{I} \lbr \Xint\tilthlongdash_B |\nabla u|^{p-\varepsilon_0} \ dx \rbr \lbr \Xint\tilthlongdash_B  \lbr \frac{|u|}{\rhoo} \rbr^{\tilde{r}} \ dx \rbr^{\frac{\varepsilon_0 - \beta}{2}\frac{p}{p-\beta}} dt.
\end{equation}
 Making use Lemma \ref{bound_g_x_t} along with \eqref{hypothesis_2}, we get
 \betagin{equation}
\lambdabel{7.15}

 \Yint\tiltlongdash_Q \tilde{g}(z)^{p-\beta} \ dz \apprle  \Yint\tiltlongdash_{\alpha Q} \lbr |\nabla u|+|\theta| +\alpha_0 \rbr ^{{{p-\beta}}} \ dz \apprle \alpha_0^{p-\beta}.

\end{equation}
For almost every $t\in \alpha I$, using H\"older's inequality and  \eqref{def_J}, we also get 
\betagin{equation}\lambdabel{7.16}
\betagin{array}{ll}
 \Xint\tilthlongdash_B  \lbr \frac{|u(x,t)|}{\rhoo} \rbr^{\tilde{r}} \ dx  & \apprle J^{\frac{\tilde{r}}{2}} \lbr \Xint\tilthlongdash_B \tilde{g}(x,t)^{p-\beta} \ dx \rbr^{\frac{2-\tilde{r}}{2}}.
 \end{array}
\end{equation}
Substituting \eqref{7.16} into \eqref{7.14}, followed by an application of H\"older's inequality with exponents $\frac{q}{p-\varepsilon_0}$ and $\tilde{q_0}:=\frac{q}{q-p+\varepsilon_0}$, we get
\betagin{equation*}
\lambdabel{4.44}
 \Yint\tiltlongdash_Q \lbr \frac{|u|}{\rhoo} \rbr^{p-\beta} \ dz \apprle J^{\frac{\varepsilon_0-\beta}{2}}  \lbr \Yint\tiltlongdash_Q |\nabla u|^q \ dz \rbr^{\frac{p-\varepsilon_0}{q}}  \lbr \Xint\tilthlongdash_{I} \lbr \Xint\tilthlongdash_B  \tilde{g}(x,t)^{p-\beta} \ dx \rbr^{\tilde{q_0} \frac{(2-\tilde{r})(\varepsilon_0 - \beta)p}{2(p-\beta)}} dt\rbr^{\frac{1}{\tilde{q_0}}}.
\end{equation*}
We further restrict $\beta_0 \leq \frac{p-\varepsilon_0}{p}$, which ensures $\tilde{q_0} \frac{(2-\tilde{r})(\varepsilon_0 - \beta)p}{2(p-\beta)} = q_0\frac{\beta(\varepsilon_0-\beta)}{2(p-\beta)}<1$ holds.

Thus Jensen's inequality for concave functions becomes applicable and hence we get
\betagin{equation}
\lambdabel{4.45}
 \Yint\tiltlongdash_Q \lbr \frac{|u|}{\rhoo} \rbr^{p-\beta} \ dz \apprle J^{\frac{\varepsilon_0-\beta}{2}}  \lbr \Yint\tiltlongdash_Q |\nabla u|^q \ dz \rbr^{\frac{p-\varepsilon_0}{q}}  \lbr  \Yint\tiltlongdash_Q  \tilde{g}(z)^{p-\beta} \ dz \rbr^{\frac{\beta(\varepsilon_0-\beta)}{2(p-\beta)}}.
\end{equation}
Using \eqref{7.15} in \eqref{4.45}, we get 
\betagin{equation*}
\lambdabel{4.46}
 \Yint\tiltlongdash_Q \lbr \frac{|u|}{\rhoo} \rbr^{p-\beta} \ dz \apprle (\alpha_0^{\beta} J)^{\frac{\varepsilon_0-\beta}{2}}  \lbr \Yint\tiltlongdash_Q |\nabla u|^q \ dz \rbr^{\frac{p-\varepsilon_0}{q}}  .
\end{equation*}

 \item[Proof of \eqref{lemma4.3.2}:] Note that in this case, we have $p-\varepsilon_0 \geq 2$, hence we make the following choice of exponents to apply Lemma \ref{interpolation_capacity_version_main}:
 \betagin{equation}
 \lambdabel{exp_ch_3}
\sigma = 2, \quad  \tilde{r} = \frac{2(p-\beta)}{p}, \quad \vartheta = p-\varepsilon_0.
 \end{equation}
 Let us now choose $\tilde{\gamma}a \in(0,1)$ so as to satisfy the following restrictions:
\betagin{equation}
- \frac{n}{2} \leq \tilde{\gamma}a\lbr 1 - \frac{n}{p-\varepsilon_0} \rbr  -  \frac{(1-\tilde{\gamma}a)np}{2(p-\beta)} \quad \text{and} \quad 
\frac{q (1-\tilde{\gamma}a) \beta}{(q-2\tilde{\gamma}a)(p-\beta)} < 1. \lambdabel{7.24}
\end{equation}
This choice of $\tilde{\gamma}a$ in \eqref{7.24} is possible, due to the following reasons.
\betagin{itemize}
 \item \eqref{7.24} holds  for $\tilde{\gamma}a=1$, since $\frac{n}{p-\varepsilon_0} < 1 + \frac{n}{2}$, as we are in the case $p-\varepsilon_0 \geq 2$.

 \item Since we assumed $p-\varepsilon_0 \geq 2$, this means that $q-2\tilde{\gamma}a >0$ for $\tilde{\gamma}a < 1$. This follows from the fact that  $ q \geq p-\varepsilon_0 \geq 2 \geq 2\tilde{\gamma}a \in (0,2).$ 

\end{itemize}
Thus we can apply Lemma \ref{interpolation_capacity_version_main} with \eqref{exp_ch_3} and $\tilde{\alphapha}=\tilde{\gamma}a$ as the choice of exponents, to get
\betagin{equation}
 \lambdabel{7.20}
\Yint\tiltlongdash_Q \lbr \frac{|u|}{\rhoo} \rbr^{2} \ dz \apprle \Xint\tilthlongdash_{I} \lbr \Xint\tilthlongdash_B |\nabla u|^{p-\varepsilon_0} \ dx \rbr^{\frac{2\tilde{\gamma}a}{p-\varepsilon_0}} \lbr \Xint\tilthlongdash_B  \lbr \frac{|u|}{\rhoo} \rbr^{\tilde{r}} \ dx \rbr^{\frac{2(1-\tilde{\gamma}a)}{r}} dt.
\end{equation}
Since $\tilde{r}<2$, we can substitute \eqref{7.16} into \eqref{7.20} to obtain
\betagin{equation}
 \lambdabel{7.22}
\Yint\tiltlongdash_Q \lbr \frac{|u|}{\rhoo} \rbr^{2} \ dz \apprle J^{1-\tilde{\gamma}a} \Xint\tilthlongdash_{I} \lbr \Xint\tilthlongdash_B |\nabla u|^{p-\varepsilon_0} \ dx \rbr^{\frac{2\tilde{\gamma}a}{p-\varepsilon_0}} \lbr \Xint\tilthlongdash_B \tilde{g}(x,t)^{p-\beta} \ dx \rbr^{\frac{(2-\tilde{r})(1-\tilde{\gamma}a)}{\tilde{r}}} dt.
\end{equation}
We will now apply H\"older inequality with exponents $\frac{q}{2\tilde{\gamma}a}$ and $\frac{q}{q-2\tilde{\gamma}a}$, and then make use of \eqref{7.24}, which will enable us to apply Jensen's inequality. Thus we get
\betagin{equation}
 \lambdabel{7.23}
 \betagin{array}{ll}
\Yint\tiltlongdash_Q \lbr \frac{|u|}{\rhoo} \rbr^{2} \ dz 

& \apprle J^{1-\tilde{\gamma}a}\lbr  \Yint\tiltlongdash_Q |\nabla u|^{q} \ dz\rbr^{\frac{2\tilde{\gamma}a}{q}} \lbr  \Yint\tiltlongdash_Q \tilde{g}(x,t)^{p-\beta} \ dz \rbr^{\frac{(2-\tilde{r})(1-\tilde{\gamma}a)}{\tilde{r}}}.
\end{array}
\end{equation}
Further making use of \eqref{7.15} and \eqref{7.23}, we get
\betagin{equation*}
 \lambdabel{7.25}
 \betagin{array}{ll}
\Yint\tiltlongdash_Q \lbr \frac{|u|}{\rhoo} \rbr^{2} \ dz & \apprle (\alpha_0^{\beta} J)^{1-\tilde{\gamma}a}\lbr  \Yint\tiltlongdash_Q |\nabla u|^{q} \ dz \rbr^{\frac{2\tilde{\gamma}a}{q}} .
\end{array}
\end{equation*}

 \item[Proof of \eqref{lemma4.3.3}:]
 In this case, we  take
\betagin{equation*}
 \lambdabel{exp_ch_4}
\sigma = 2, \quad  \tilde{r} = \frac{2(p-\beta)}{p}, \quad \vartheta = p-\varepsilon_0 \qquad \text{and} \qquad \tilde{\alphapha} = \frac{p-\varepsilon_0}{2} < 1. 
 \end{equation*}
Proceeding as in the proof of \eqref{lemma4.3.1}, we get
\betagin{equation*}
 \lambdabel{7.26}
 \betagin{array}{ll}
\Yint\tiltlongdash_Q \lbr \frac{|u|}{\rhoo} \rbr^{2} \ dz & \apprle (\alpha_0^{\beta} J)^{\frac{2-p+\varepsilon_0}{2}}\lbr  \Yint\tiltlongdash_Q |\nabla u|^{q} \ dz \rbr^{\frac{p-\varepsilon_0}{q}} .
\end{array}
\end{equation*}
\end{description}

\end{proof}

\betagin{lemma}
\lambdabel{lemma6.2}
 Let $\frac{2n}{n+2} < p < 2+\beta$.  Then under the assumptions of Lemma \ref{Caccioppoli}, there holds
 \[
  \Yint\tiltlongdash_Q \lbr \frac{|u|}{\rhoo} \rbr^2 \ dz \leq  C_{(n,p,q,\Lambda_0,\Lambda_1,b_0,\kappa)}\alpha_0^2. 
 \]
\end{lemma}
\betagin{proof}
  Let $1 \leq \tilde{\alphapha}_1 < \tilde{\alphapha}_2 \leq 8$ be any two fixed numbers. We can then apply Lemma \ref{lemma7.4} to get $\tilde{\gamma}a \in (0,1)$ such that the following estimate holds:
  \betagin{equation}
  \lambdabel{eq7.31}

\Yint\tiltlongdash_{ \tilde{\alphapha}_1 Q} \lbr \frac{|u|}{\rhoo} \rbr^{2} \ dz  \apprle (\alpha_0^{\beta} J)^{1-\tilde{\gamma}a}\lbr  \Yint\tiltlongdash_{ \tilde{\alphapha}_2Q} |\nabla u|^{q} \ dz\rbr^{\frac{2\tilde{\gamma}a}{q}}, 

\end{equation}
where $J:= \sup_{t \in \tilde{\alphapha}_1I} \Xint\tilthlongdash_{\tilde{\alphapha}_1 B} \lbr \frac{|u|}{\rhoo} \rbr^2 \tilde{g}(x,t)^{-\beta} \ dx $ with 
$$\tilde{g}(z) := \max\left\{ \mathcal{M}( |\nabla \tilde{u}|^q )^{\frac{1}{q}}(z), \mathcal{M}( \lbr |\nabla u|^q+|\theta|^q\rbr \lsb{\chii}{\tilde{\alphapha}_2 Q})^{\frac{1}{q}}(z), \alpha_0\right\}$$ 
where we have set $\tilde{u}(x,t) = u(x,t) \eta(x)\zetata(t)$ for  $\eta(x) \zetata(t) \in C_c^{\infty}(\tilde{\alphapha}_2 Q)$ with $\eta(x) \zetata(t) \equiv 1$ for  $(x,t) \in \tilde{\alphapha}_1Q$ instead of \eqref{cut_off_function}.

We can now make  use of Lemma \ref{Caccioppoli} and  obtain the following estimate for some $ \beta \in (0,\beta_0]$:
 \betagin{equation}
  \lambdabel{conclusion_1_rescaled}
  \betagin{array}{ll}
  \alpha_0^{p-\beta} + \alpha_0^{p-2} \sup_{t \in\tilde{\alphapha}_1I} \Xint\tilthlongdash_{\tilde{\alphapha}_1B} \left| \frac{u(x,t)}{\rhoo} \right|^2 \tilde{g}(x,t)^{-\beta} \ dx &\leq C \Yint\tiltlongdash_{\tilde{\alphapha}_2Q}  \alpha_0^{p-2-\beta} \left| \frac{u}{(\tilde{\alphapha}_2-\tilde{\alphapha}_1)\rhoo} \right|^2\ dz \\
  &\quad +\Yint\tiltlongdash_{\tilde{\alphapha}_2Q}  \left[\left| \frac{u}{(\tilde{\alphapha}_2-\tilde{\alphapha}_1)\rhoo} \right|^{p-\beta} + |\theta|^{p-\beta}\right] \ dz.
  \end{array}
 \end{equation}

 
 
 
In particular, \eqref{conclusion_1_rescaled} implies
\betagin{equation}
 \lambdabel{eq7.32}
\alpha_0^{\beta} J \apprle  \Yint\tiltlongdash_{\tilde{\alphapha}_2Q}\left[   \left| \frac{u}{(\tilde{\alphapha}_2-\tilde{\alphapha}_1)\rhoo} \right|^2  + \alpha_0^{2-p+\beta} \left| \frac{u}{(\tilde{\alphapha}_2-\tilde{\alphapha}_1)\rhoo} \right|^{p-\beta} + |\theta|^{p-\beta} \right]dz.
\end{equation}

We now apply Young's inequality to the middle term on the right hand side of \eqref{eq7.32} (possible since $p-\beta < 2$) and make use of \eqref{hypothesis_1}, to get 

\betagin{equation}
 \lambdabel{eq7.34}
 \alpha_0^{\beta} J \apprle  \Yint\tiltlongdash_{\tilde{\alphapha}_2Q}   \left| \frac{u}{(\tilde{\alphapha}_2-\tilde{\alphapha}_1)\rhoo} \right|^2 \ dz+  \alpha_0^2.
\end{equation}

 
 
 Substituting \eqref{eq7.34} into \eqref{eq7.31} followed by 

applying Young's inequality with exponent $\frac{1}{1-\tilde{\gamma}a} $ and $\frac{1}{\tilde{\gamma}a}$, along with  using  \eqref{hypothesis_1} and the restriction $\tilde{\alphapha}_1,\tilde{\alphapha}_2 \in [1,8]$, we get
\betagin{equation}
  \lambdabel{eq7.36}
  \betagin{array}{ll}
   \Yint\tiltlongdash_{ \tilde{\alphapha}_1 Q} \lbr \frac{u}{\rhoo} \rbr^{2} \ dz  

   & \leq \frac{1}{2} \Yint\tiltlongdash_{\tilde{\alphapha}_2 Q}   \left( \frac{u}{\rhoo} \right)^2   dz + C\alpha_0^2 (\tilde{\alphapha}_2-\tilde{\alphapha}_1)^{-\frac{2(1-\tilde{\alphapha})}{\tilde{\alphapha}}} + C \alpha_0^2 (\tilde{\alphapha}_2 - \tilde{\alphapha}_1)^2 \\
   &\leq  \frac{1}{2} \Yint\tiltlongdash_{\tilde{\alphapha}_2 Q}   \left( \frac{u}{\rhoo} \right)^2   dz + C \alpha_0^2 (\tilde{\alphapha}_2-\tilde{\alphapha}_1)^{-\frac{2(1-\tilde{\alphapha})}{\tilde{\alphapha}}} + C \alpha_0^2.
  \end{array}
 \end{equation}

 If we set $\varpi(a) := \Yint\tiltlongdash_{a Q}   \left| u \right|^2   dz$, then we see that for $\tilde{\alphapha}_1=:a < b:=\tilde{\alphapha}_2$,  \eqref{eq7.36} can be rewritten as 
 \[
  \varpi(a) \leq \frac12 \varpi(b) + C\alpha_0^2 \rhoo^2 (b-a)^{-\frac{2(1-\tilde{\alphapha})}{\tilde{\alphapha}}} + C\alpha_0^2 \rhoo^2.
 \]
We can now apply Lemma \ref{iteration_lemma} along with the hypothesis $ \tilde{\alphapha}_1, \tilde{\alphapha}_2 \in [1,8]$ to get 
\[
 \Yint\tiltlongdash_{ Q}   \left| u \right|^2   dz \apprle  \alpha_0^2 \rhoo^2.
\]

This completes the proof of the lemma. 
\end{proof}

\subsection{Reverse H\"older inequality}

\betagin{lemma}
 \lambdabel{reverse_holder_inequality}
 Let $\kappa \geq 1$.  Then there exists  $\beta_0(n,p,q,\Lambda_0,\Lambda_1,b_0,\kappa)$ and $C{(n,p,q,\Lambda_0,\Lambda_1,b_0,\kappa)}$  such that the following holds: Let  $u \in L^2(-T,T;L^2(\Omega))\cap L^{p-\beta}(-T,T;W_0^{1,p-\beta}(\Omega))$ for some $\beta \in (0,\beta_0)$ be a very weak solution of \eqref{main}. Furthermore, let  $Q=Q_{\rhoo,s}$ be a parabolic cylinder with $s = \alpha_0^{2-p} \rhoo^2$ as in \eqref{def_s} be given. Suppose the following bounds hold:
 \betagin{gather}
  \frac{\alpha_0^{p-\beta}}{\kappa} \leq \Yint\tiltlongdash_Q (|\nabla u|+|\theta|)^{p-\beta} \ dz \qquad \text{and} \qquad   \Yint\tiltlongdash_{8^3Q} (|\nabla u|+|\theta|)^{p-\beta} \ dz \leq \kappa \alpha_0^{p-\beta}. \lambdabel{hyp_2}
 \end{gather}
Then the following conclusion holds:
\betagin{equation*}
 \lambdabel{conclusion_2}
 \alpha_0^{p-\beta} \leq C \lbr \Yint\tiltlongdash_{8Q} |\nabla u|^q \ dz \rbr^{\frac{p-\beta}{q}} + \Yint\tiltlongdash_{8Q} |\theta|^{p-\beta} \ dz.
\end{equation*}
\end{lemma}

\betagin{proof}
 Since all the hypothesis of  Lemma \ref{Caccioppoli} (note the scaling factor is $8^2$ instead of $8$ as in Lemma \ref{cac1}) is satisfied, we  have
 \betagin{equation}

\lambdabel{eq7.41}
  \alpha_0^{p-\beta}  \apprle \Yint\tiltlongdash_{8^2Q}\lbr  \alpha_0^{p-2-\beta} \left| \frac{u}{\rhoo} \right|^2 + \left| \frac{u}{\rhoo} \right|^{p-\beta} + |\theta|^{p-\beta}\rbr dz.
 \end{equation}
 
 
 Let $\tilde{g}(z) := \max\left\{ \mathcal{M}( |\nabla \tilde{u}|^q )^{\frac{1}{q}}(z), \mathcal{M}(( |\nabla u|^q |\theta|^q)\lsb{\chii}{8^2 Q})^{\frac{1}{q}}(z),\alpha_0\right\}$, where $\tilde{u}(x,t):= u(x,t) \eta(x) \zetata(t)$ with $\eta(x) \zetata(t) \in C_c^{\infty}(8^2Q)$ satisfying $\eta(x) \zetata(t) \equiv 1$ for $(x,t) \in 8Q$.  Let us set 
 \betagin{equation*}
  \lambdabel{def_J_new}
  J= \sup_{t \in8I} \Xint\tilthlongdash_{8B} \left| \frac{u(x,t)}{\rhoo} \right|^2 \tilde{g}(x,t)^{-\beta} \ dx.
 \end{equation*}

 
Applying a rescaled version of Lemma \ref{Caccioppoli}, we have
\betagin{equation}
 \lambdabel{J_bound}
 \alpha_0^{p-2} J \apprle \Yint\tiltlongdash_{8^2Q}\lbr  \alpha_0^{p-2-\beta} \left| \frac{u}{\rhoo} \right|^2 + \left| \frac{u}{\rhoo} \right|^{p-\beta}+|\theta|^{p-\beta}\rbr dz.
\end{equation}

 
 We now have the following three bounds:
 \betagin{itemize}
  \item In the case $2 \geq p-\beta$, we can apply Lemma \ref{lemma6.2} (this is possible since \eqref{hyp_2} is on $8^3Q$) to get 
  \betagin{equation}
   \lambdabel{eest_1}
   \Yint\tiltlongdash_{8^2Q} \lbr \frac{|u|}{\rhoo} \rbr^2 \ dz \apprle \alpha_0^2.
  \end{equation}
  \item In the case $2 \leq p-\beta$, for  $\sigma \in \{2,p-\beta\}$,  we can apply H\"older's inequality and Theorem \ref{sobolev-poincare}  along with \eqref{hyp_2}, to get
\betagin{equation}
 \lambdabel{eest_2}
 \Yint\tiltlongdash_{8Q} \lbr \frac{|u|}{\rhoo} \rbr^{\sigma} \ dz  \leq \lbr \Yint\tiltlongdash_{8Q} \lbr \frac{|u|}{\rhoo} \rbr^{p-\beta} \ dz\rbr^{\frac{\sigma}{p-\beta}} \apprle  \lbr \Yint\tiltlongdash_{8Q} |\nabla u|^{p-\beta} \ dz\rbr^{\frac{\sigma}{p-\beta}} {\apprle} \alpha_0^{\sigma}.
\end{equation}

 \end{itemize}
 
 Substituting \eqref{eest_1} and \eqref{eest_2}  into \eqref{J_bound} gives
  \betagin{equation}

\lambdabel{eq7.42}

 \alpha^{\beta} J \apprle \alpha_0^2. 
 \end{equation}

 
Substituting \eqref{eq7.42} into a rescaled version of Lemma \ref{lemma7.4} and applying Young's inequality followed by substituting the resulting estimate into \eqref{eq7.41}, we get

\betagin{equation*}
 \betagin{array}{@{}r@{}c@{}l@{}}
 \alpha_0^{p-\beta} & \apprle &\alpha_0^{p-2-\beta}\alpha_0^{2(1-\tilde{\gamma}a)}\lbr  \Yint\tiltlongdash_{ 8^2Q} |\nabla u|^{q} \ dz\rbr^{\frac{2\tilde{\gamma}a}{q}} + \alpha_0^{\frac{\varepsilon_0-\beta}{2}}  \lbr \Yint\tiltlongdash_{8^2 Q} |\nabla u|^q \ dz \rbr^{\frac{p-\varepsilon_0}{q}}  + \Yint\tiltlongdash_{8^2Q} |\theta|^{p-\beta} \ dz \\

 &\apprle & \varepsilon_1 \alpha_0^{p-\beta} + C(\varepsilon_1) \lbr  \Yint\tiltlongdash_{ 8^2Q} |\nabla u|^{q} \ dz\rbr^{\frac{p-\beta}{q}} +  C(\varepsilon_1) \lbr  \Yint\tiltlongdash_{ 8^2Q} |\nabla u|^{q} \ dz \rbr^{\frac{p-\beta}{q}} +  \Yint\tiltlongdash_{8^2Q} |\theta|^{p-\beta} \ dz .
 \end{array}
\end{equation*}
Choosing $\varepsilon_1$ small enough, we get the desired conclusion. This completes the proof of the lemma. 
\end{proof}

\section{Proof of the main result}\lambdabel{main_estimates}

The proof of Theorem \ref{main_theorem} now follows exactly as in proof of \cite[Theorem 2.1]{Bog}. For the sake of completeness, we present the details below. Without loss of generality, we can assume $\rhoo = 1$ and $z_0 = (0,0) \in \partial \Omega \times (-T,T)$. We take $Q_1 := Q_{(0,0)}(1,1) $ and  $Q_2 := Q_{(0,0)}(2,2^2)$

and for any  $z \in Q_2$, define the parabolic distance of $z$ to $\partial Q_2$ by 
\betagin{gather*}
d_p(z):= \inf_{ \tilde{z} \in \mathbb{R}^{n+1} \setminus Q_2} \min \{ |x-\tilde{x}| , \sqrt{|t-\tilde{t}|} \} \lambdabel{par_dist_boundary}. 
\end{gather*}

Furthermore, let $\beta_0$ be the constant from Lemma \ref{reverse_holder_inequality} and $\beta \in (0,\beta_0]$ be such that $p-\beta > \frac{2n}{n+2}$. For $z \in Q_2$, let us define the following function:
\betagin{gather}
  \psi(z) := |\nabla u(z)|+ |\theta(z)| \quad \text{and} \quad f(z) := d_p^{\alpha} (z) \psi(z) \quad \text{with} \ \alpha:= \frac{n+2}{d}, \lambdabel{def_g_f} 
\end{gather}
where $d$ is as defined in \eqref{de_d}. 
Finally  we define $\alpha_0$ to be
\betagin{gather}
 \alpha_0^d := \Yint\tiltlongdash_{Q_2} \psi(z)^{p-\beta} \ dz + 1. \lambdabel{def_al_0}
\end{gather}
Let $\lambda$ be any number such that 
\betagin{gather}
 \lambda \geq b^{\frac{1}{d}} \alpha_0 \qquad \text{where} \ b:= 2^{10(n+2)}.\lambdabel{5.5}
\end{gather}
Now suppose that $\mathfrak{z} \in Q_2$ with $f(\mathfrak{z})>0$, then let us denote the parabolic distance of $\mathfrak{z}$ to $\partial Q_2$ by $r_{\mathfrak{z}} := d_p(\mathfrak{z}) $

and define the intrinsic scaling factor as
\betagin{gather}
 \gamma = \gamma(\mathfrak{z}) := (r_{\mathfrak{z}}^{-\alpha} \lambda )^{2-p} = (d_p(\mathfrak{z})^{-\alpha} \lambda)^{2-p} \lambdabel{intrinsic_scaling}.
\end{gather}

In order to prove higher integrability, we want to apply Lemma \ref{gehring_lemma}. So the rest of the proof is devoted to ensuring that all the hypotheses of Lemma \ref{gehring_lemma} are satisfied.
\betagin{description}[leftmargin=*]
 \item[Case $p \geq 2$:]  Let us note that $r_{\mathfrak{z}}^{\alpha} \leq 2^{\alpha} \leq b^{\frac1{d}} \alpha_0 \leq \lambda$, which implies $\gamma = (r_{\mathfrak{z}}^{-\alpha} \lambda)^{2-p} \leq 1$. Hence we shall consider intrinsic cylinders of the type $Q_{\mathfrak{z}}(R, \gamma R^2)$ with $0 < R \leq r_{\mathfrak{z}}$.  
 
 
 In order to apply Lemma \ref{gehring_lemma}, we need to find an appropriate intrinsic parabolic cylinder around $\mathfrak{z}$ on which all the hypotheses of Lemma \ref{reverse_holder_inequality} are satisfied. 
 In order to do this, let us first take $R$ such that $r_{\mathfrak{z}} \leq 2^9 R < 2^9 r_{\mathfrak{z}}$. In this case, there holds:
 \betagin{equation}
  \lambdabel{5.8}
  \betagin{array}{ll}
   \Yint\tiltlongdash_{Q_{\mathfrak{z}}(R,\gamma R^2)} \psi(z)^{p-\beta} \ dz & \leq \frac{|Q_2|}{|{Q_{\mathfrak{z}}(R,\gamma R^2)}|} \Yint\tiltlongdash_{Q_{2}} \psi(z)^{p-\beta} \ dz \\
   &\overset{\text{\eqref{def_al_0}}}{=} \frac{2^{n+2}}{R^{n+2} \gamma}\alpha_0^d
   \overset{\text{\eqref{5.5}}}{\leq} \frac{2^{10(n+2)}}{r_{\mathfrak{z}}^{n+2}} \frac{\lambda^d}{b} 
    \overset{\text{\eqref{intrinsic_scaling}}}{=} (r_{\mathfrak{z}}^{-\alpha} \lambda)^{p-\beta}.
  \end{array}
 \end{equation}

 
 Furthermore, by Lebesgue differentiation theorem, for every $\mathfrak{z} \in Q_2$ with $f(\mathfrak{z}) > \lambda$, there holds
 \betagin{equation}
  \lambdabel{5.9}
  \lim_{r\searrow 0} \Yint\tiltlongdash_{Q_{\mathfrak{z}}(r,\gamma r^2)} \psi(z)^{p-\beta} \ dz = \psi(\mathfrak{z})^{p-\beta} \overset{\text{\eqref{def_g_f}\eqref{intrinsic_scaling}}}{=} (r_{\mathfrak{z}}^{-\alpha} f(\mathfrak{z}))^{p-\beta} >  (r_{\mathfrak{z}}^{-\alpha} \lambda)^{p-\beta} .
 \end{equation}

 Thus from \eqref{5.8} and \eqref{5.9}, we observe that there should exist  $\rhoo \in \lbr0,\frac{r_{\mathfrak{z}}}{2^9}\rbr$, such that 
 \betagin{equation*}
  \lambdabel{def_rho}
  \betagin{array}{c}
  \Yint\tiltlongdash_{Q_{\mathfrak{z}}(\rhoo, \gamma \rhoo^2)} \psi(z)^{p-\beta} \ dz = (r_{\mathfrak{z}}^{-\alpha} \lambda)^{p-\beta}, \\
  \Yint\tiltlongdash_{Q_{\mathfrak{z}}(R, \gamma R^2)} \psi(z)^{p-\beta} \ dz \leq (r_{\mathfrak{z}}^{-\alpha} \lambda)^{p-\beta},  \quad \ \forall\   R \in [\rhoo, r_{\mathfrak{z}}].
  \end{array}
 \end{equation*}

  We now set $Q:= Q_{\mathfrak{z}}(\rhoo, \gamma \rhoo^2)$, then $2^9Q \subset Q_2$, thus all hypotheses of Lemma \ref{reverse_holder_inequality} are satisfied with $(r_{\mathfrak{z}}^{-\alpha} \lambda, 1)$ instead of $(\alpha_0,\kappappa)$, i.e., the following holds:
  \betagin{equation}
  \lambdabel{5.11}
   (r_{\mathfrak{z}}^{-\alpha} \lambda)^{p-\beta}  = \Yint\tiltlongdash_Q \psi(z)^{p-\beta} \ dz \qquad \text{and} \qquad \Yint\tiltlongdash_{2^8Q} \psi(z)^{p-\beta} \ dz \apprle (r_{\mathfrak{z}}^{-\alpha} \lambda)^{p-\beta}.
  \end{equation}

Thus we can apply  Lemma \ref{reverse_holder_inequality}, which gives
\betagin{equation}
 \lambdabel{5.12}
 (r_{\mathfrak{z}}^{-\alpha} \lambda)^{p-\beta} \apprle \lbr \Yint\tiltlongdash_{2^8Q} |\nabla u|^q\ \ dz \rbr^{\frac{p-\beta}{q}} + \Yint\tiltlongdash_{2^8Q} |\theta|^{p-\beta} \ dz.
\end{equation}
 
 
 
 Since $2^9 \rhoo \leq r_{\mathfrak{z}}$ and $\gamma \leq 1$, we also have for all $z \in 2^8Q$ that 
 \betagin{gather}
  d_p(z) \leq \min\{ r_{\mathfrak{z}} + 2^8 \rhoo, \sqrt{r_{\mathfrak{z}}^2 + \gamma (2^8 \rhoo)^2} \} \leq \frac32 r_{\mathfrak{z}}, \lambdabel{up_bound_dpz} \\
  d_p(z) \geq \min\{ r_{\mathfrak{z}} - 2^8 \rhoo, \sqrt{r_{\mathfrak{z}}^2 - \gamma (2^8 \rhoo)^2} \} \geq \frac12 r_{\mathfrak{z}}. \lambdabel{low_bound_dpz}
 \end{gather}

 Now substituting \eqref{up_bound_dpz} and \eqref{low_bound_dpz} into \eqref{def_g_f}, we find
 \betagin{equation}
  \lambdabel{5.15}
  c^{-1} f(z) \leq r_{\mathfrak{z}}^{\alpha} \psi(z) \leq c f(z), \quad \forall \ z \in 2^8Q \ \text{ with  } \ c = c(n)>1.
 \end{equation}

 We now claim the following estimate holds:
 \betagin{equation}
  \lambdabel{5.16}
  \lambda^{p-\beta} \overset{\text{\redlabel{5.16.a}{a}}}{\apprle} \Yint\tiltlongdash_{2^8Q} f(z)^{p-\beta} \ dz \overset{\text{\redlabel{5.16.b}{b}}}{\apprle} \lbr \Yint\tiltlongdash_{2^8Q} f(z)^q \ dz \rbr^{\frac{p-\beta}{q}} + \Yint\tiltlongdash_{2^8Q} (r_{\mathfrak{z}}^{\alpha} |\theta(z)|)^{p-\beta} \ dz \overset{\text{\redlabel{5.16.c}{c}}}{\apprle} \lambda^{p-\beta}.
 \end{equation}
 \betagin{description}[leftmargin=*]
  \item[Estimate \redref{5.16.a}{a}:] This follows easily by making use of \eqref{5.11} and \eqref{5.15} and subsequently enlarging the parabolic cylinder $Q$. 
  \item[Estimate \redref{5.16.b}{b}:] This is obtained by the following chain of estimates:
  \betagin{equation*}
   \betagin{array}{lll}
    \Yint\tiltlongdash_{2^8Q} f(z)^{p-\beta} \ dz & \overset{\text{\eqref{5.15}}}{\apprle} & r_{\mathfrak{z}}^{\alpha(p-\beta)} \Yint\tiltlongdash_{2^8Q} \psi(z)^{p-\beta} \ dz \overset{\text{\eqref{5.11}}}{\apprle} \lambda^{p-\beta} \\
    & \overset{\text{\eqref{5.12}}}{\apprle} &  r_{\mathfrak{z}}^{\alpha(p-\beta)} \lbr \Yint\tiltlongdash_{2^8Q} |\nabla u|^q \ dz \rbr^{\frac{p-\beta}{q}}  + r_{\mathfrak{z}}^{\alpha(p-\beta)} \Yint\tiltlongdash_{2^8Q} |\theta(z)|^{p-\beta} \ dz \\
    & \overset{\text{\eqref{5.15}}}{\apprle} & \lbr \Yint\tiltlongdash_{2^8Q} f(z)^q \ dz \rbr^{\frac{p-\beta}{q}} + \Yint\tiltlongdash_{2^8Q} (r_{\mathfrak{z}}^{\alpha} |\theta(z)|)^{p-\beta} \ dz.
   \end{array}
  \end{equation*}

  \item[Estimate \redref{5.16.c}{c}:] This follows by applying Jensen's inequality (since $q < p-\beta$) along with the bound \eqref{5.11} to get:
  \betagin{equation*}
   \lbr \Yint\tiltlongdash_{2^8Q} f(z)^q \ dz \rbr^{\frac{p-\beta}{q}} + \Yint\tiltlongdash_{2^8Q} (r_{\mathfrak{z}}^{\alpha} |\theta(z)|)^{p-\beta} \ dz \overset{\text{\eqref{5.15}}}{\apprle}  \Yint\tiltlongdash_{2^8Q} f(z)^{p-\beta} \ dz \overset{\text{\eqref{5.11}}}{\apprle} \lambda^{p-\beta}.
  \end{equation*}
 \end{description}

 
 Thus  \eqref{5.16} holds and  as a consequence, we can apply Lemma \ref{gehring_lemma} over $2^8Q$ to see that  for any $\beta \in (0,\beta_0]$, there exists  $\delta_0 = \delta_0(n,p,b_0,\varepsilon_0,\Lambda_0,\Lambda_1)>0$ such that $f \in L_{loc}^{p-\beta+\delta_1}(Q_2)$ with $\delta_1 = \min\{ \delta_0, \tilde{p}-p+\beta\}$. This is quantified by the estimate:
\betagin{equation*}
 \lambdabel{5.19}
 \iint_{Q_2} f(z)^{p-\beta+\delta} \ dz \apprle \alpha_0^{\delta} \iint_{Q_2} f(z)^{p-\beta} \ dz + \iint_{Q_2} (r_{\mathfrak{z}}^{\alpha} |\theta(z)|)^{p-\beta+\delta} \ dz\qquad \forall \delta \in (0,\delta_1].
\end{equation*}

By iterating the previous arguments, for any $\beta \in (0,\varepsilon_{\text{geh}}]$ where $\varepsilon_{\text{geh}}>0$ is the gain in higher integrability coming from Lemma \ref{gehring_lemma}, we obtain the bound
\betagin{equation}
 \lambdabel{5.20}
 \iint_{Q_2} f(z)^{p} \ dz \apprle \alpha_0^{\beta} \iint_{Q_2} f(z)^{p-\beta} \ dz + \Yint\tiltlongdash_{Q_2} |\theta(z)|^p \ dz.
\end{equation}

For any $z \in Q_1$, we have $d_p(z) \geq \min \{ 1 , \sqrt{3}\} \geq 1$, $\frac{|Q_2|}{|Q_1|} = C(n)$, which implies the following bounds  hold:
\betagin{gather}
 |\nabla u (z)| \leq \psi(z) \leq f(z) \qquad \forall\  z \in Q_1, \lambdabel{5.21}\\
 \Yint\tiltlongdash_{Q_1} |\nabla u|^{p+\beta} \ dz \apprle_n \Yint\tiltlongdash_{Q_2} f(z)^p \ dz ,\lambdabel{5.22}\\
 f(z) \leq 2^{\alpha} \psi(z) \qquad \forall z \in Q_2,  \quad \text{since} \ d_p(z) \leq 2.\lambdabel{5.23}
\end{gather}

Using \eqref{5.21}, \eqref{5.22}, \eqref{5.23} along with \eqref{5.20} and making use of \eqref{def_g_f} and \eqref{def_al_0}, we get
\betagin{equation*}
 \Yint\tiltlongdash_{Q_1} |\nabla u|^p \ dz \apprle \alpha_0^{\beta} \Yint\tiltlongdash_{Q_2} \psi(z)^{p-\beta} \ dz  + \Yint\tiltlongdash_{Q_2} \theta^p \ dz \apprle \lbr \Yint\tiltlongdash_{Q_2} \lbr |\nabla u| + |\theta| \rbr^{p-\beta} \ dz \rbr^{1+\frac{\beta}{d}} + \Yint\tiltlongdash_{Q_2} (1+\theta^p) \ dz.
\end{equation*}
This proves the asserted estimate.

 \item[Case $\frac{2n}{n+2} <p < 2$:] The basic change with respect to the case $p\geq 2$ is that, we now switch to the sub-quadratic scaling, i.e., we consider intrinsic cylinders of the type $Q_{\mathfrak{z}}(\gamma^{-\frac12}R,R^2)$. 
 
 The parameter $\alpha_0$ is still given by \eqref{def_al_0} and $\lambda$ is chosen as in \eqref{5.5} and $\gamma$ is again given as in \eqref{intrinsic_scaling}, where $\mathfrak{z} \in Q_2$ with $f(z) > \lambda$. But in contrast to $p\geq 2$ case, we have $\gamma = (r_{\mathfrak{z}}^{-\alpha} \lambda )^{2-p} \geq 1$. Hence for $R \in (0,r_{\mathfrak{z}})$, we have $Q_{\mathfrak{z}}(\gamma^{-\frac12}R,R^2)\subset Q_2$. Now once again, in order to apply Lemma \ref{gehring_lemma}, we need to find a suitable intrinsic parabolic cylinder around $\mathfrak{z}$, which enables us to apply Lemma \ref{reverse_holder_inequality}.  We observe, from the definition \eqref{def_g_f} that $n+2 = \alpha d$ (recall $d$ is defined as in \eqref{de_d}) and $(2-p)\frac{n}{2} + d = p-\beta$, which gives
 \betagin{equation*}

  \betagin{array}{ll}
   \Yint\tiltlongdash_{Q_{\mathfrak{z}}(\gamma^{-\frac12}R, R^2)} \psi(z)^{p-\beta} \ dz & \apprle \frac{|Q_2|}{|Q_{\mathfrak{z}}(\gamma^{-\frac12}R, R^2)|}\Yint\tiltlongdash_{Q_2} \psi(z)^{p-\beta} \ dz  = \frac{2^{n+2}}{R^{n+2} \gamma^{-\frac{n}{2}}} \alpha_0^d  

   \apprle (r_{\mathfrak{z}}^{-\alpha} \lambda )^{p-\beta}. 
  \end{array}
 \end{equation*}
 Now we can continue as in the $p \geq 2$ case to obtain the desired conclusion. 
\end{description}
This completes the  proof. 
 \qed

\section*{References}

\betagin{thebibliography}{10}

\bibitem{AH}
David~R Adams and Lars~I Hedberg.
\newblock {\em Function Spaces and Potential Theory}, volume 314.
\newblock Grundlehren der Mathematischen Wissenschaften, Springer, Berlin,
  1996.

\bibitem{AP1}
Karthik Adimurthi and Nguyen~Cong Phuc.
\newblock Global lorentz and lorentz-morrey estimates below the natural
  exponent for quasilinear equations.
\newblock {\em Calculus of Variations and Partial Differential Equations},
  54(3):3107--3139, 2015.

\bibitem{Anc}
Alano Ancona.
\newblock On strong barriers and an inequality of hardy for domains in
  {$\mathbb{R}^n$}.
\newblock {\em Journal of the London Mathematical Society}, 2(2):274--290,
  1986.

\bibitem{bogelein2007thesis}
Verena B{\"o}gelein.
\newblock {\em Regularity results for weak and very weak solutions of higher
  order parabolic systems}.
\newblock PhD thesis, 2007.

\bibitem{Bog}
Verena B{\"o}gelein.
\newblock Very weak solutions of higher-order degenerate parabolic systems.
\newblock {\em Advances in Differential Equations}, 14(1/2):121--200, 2009.

\bibitem{bogelein2013regularity}
B{\"o}gelein, Verena and Duzaar, Frank and Mingione, Giuseppe.
\newblock The regularity of general parabolic systems with degenerate diffusion.
\newblock {\em Memoirs of American Mathematical Society}, 221, 2013.

\bibitem{Prato}
Guiseppe Da~Prato.
\newblock {$\mathcal{L}^{(p,\thetaeta)}(\omegaega,\deltalta)$} e loro propriet\'a.
\newblock {\em Annali di Matematica Pura ed Applicata}, 69(1):383--392, 1965.

\bibitem{DiB1}
Emmanuele DiBenedetto.
\newblock {\em Degenerate Parabolic Equations}.
\newblock Universitext, Springer-Verlag, New York, 1993.

\bibitem{Seb}
L~Diening, S~Schwarzacher, B~Stroffolini, and A~Verde.
\newblock Parabolic lipschitz truncation and caloric approximation.
\newblock {\em arXiv preprint
  \href{https://arxiv.org/pdf/1606.01706}{arXiv:1606.01706}}, 2016.

\bibitem{diening2010existence}
Lars Diening, Michael R{\u{u}}{\v{z}}i{\v{c}}ka, and J{\"o}rg Wolf.
\newblock Existence of weak solutions for unsteady motions of generalized
  newtonian fluids.
\newblock 2010.

\bibitem{Gia}
Mariano Giaquinta.
\newblock {\em Introduction to Regularity Theory for Nonlinear Elliptic
  Systems}.
\newblock Birkh\"auser Verlag, ETH Z\"urich, Basel.

\bibitem{Gstr}
Mariano Giaquinta and Michael Struwe.
\newblock On the partial regularity of weak solutions of nonlinear parabolic
  systems.
\newblock {\em Mathematische Zeitschrift}, 179(4):437--451, 1982.

\bibitem{PG}
Przemys{\l}aw G{\'o}rka.
\newblock Campanato theorem on metric measure spaces.
\newblock {\em Ann. Acad. Sci. Fenn. Math}, 34(2):523--528, 2009.
\newblock
  \href{http://www.acadsci.fi/mathematica/Vol34/vol34pp523-528.pdf}{http://www.acadsci.fi/mathematica/Vol34/vol34pp523-528.pdf}.

\bibitem{Grafakos}
Loukas Grafakos.
\newblock {\em Classical fourier analysis}, volume~2.
\newblock Springer, 2008.

\bibitem{HKM}
Juha Heinonen, Tero Kilpel\"ainen, and Olio Martio.
\newblock {\em Nonlinear potential theory of degenerate elliptic equations}.
\newblock Oxford Mathematical Monographs, Oxford University Press, New York,
  1993.

\bibitem{KK94}
T~Kilpel{\"a}inen and P~Koskela.
\newblock Global integrability of the gradients of solutions to partial
  differential equations.
\newblock {\em Nonlinear Analysis: Theory, Methods \& Applications},
  23(7):899--909, 1994.

\bibitem{KL1}
Juha Kinnunen and John~L Lewis.
\newblock Higher integrability for parabolic systems of {$p$}-laplacian type.
\newblock {\em Duke Mathematical Journal}, 102(2):253--272, 2000.

\bibitem{KL}
Juha Kinnunen and John~L Lewis.
\newblock Very weak solutions of parabolic systems of {$p$}-laplacian type.
\newblock {\em Arkiv f{\"o}r Matematik}, 40(1):105--132, 2002.

\bibitem{Le88}
John~L Lewis.
\newblock Uniformly fat sets.
\newblock {\em Transactions of the American Mathematical Society},
  308(1):177--196, 1988.

\bibitem{Gary}
Gary~M Lieberman.
\newblock {\em Second order parabolic differential equations}.
\newblock World scientific, 1996.

\bibitem{Maz}
Vladimir~G Maz'ya.
\newblock {\em Sobolev Spaces}.
\newblock Springer Series in Soviet Mathematics, Springer-Verlag, Berlin, 1985.

\bibitem{Mik}
Pasi Mikkonen.
\newblock On the wolff potential and quasilinear elliptic equations involving
  measures.
\newblock {\em Ann. Acad. Sci. Fenn. Math. Diss}, 104, 1996.

\bibitem{Par}
Mikko Parviainen.
\newblock {\em Global higher integrability for nonlinear parabolic partial
  differential equations in nonsmooth domains}.
\newblock Helsinki University of Technology, 2007.
\newblock
  \href{http://lib.tkk.fi/Diss/2007/isbn9789512289400/}{http://lib.tkk.fi/Diss/2007/isbn9789512289400/}.

\bibitem{Mik1}
Mikko Parviainen.
\newblock Reverse h{\"o}lder inequalities for singular parabolic equations near
  the boundary.
\newblock {\em Journal of Differential Equations}, 246(2):512--540, 2009.

\end{thebibliography}

\end{document}