\documentclass[reqno]{amsart}
\usepackage{hyperref}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2013 (2013), No. 93, pp. 1--24.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2013 Texas State University - San Marcos.}
\vspace{9mm}}

\begin{document}
\title[\hfilneg EJDE-2013/93\hfil A gradient estimate for solutions]
{A gradient estimate for solutions  to parabolic equations
 with discontinuous coefficients}

\author[J. Fan, K. Kim, S. Nagayasu, G. Nakamura \hfil EJDE-2013/93\hfilneg]
{Jishan Fan, Kyoungsun Kim, Sei Nagayasu, Gen Nakamura}  % in alphabetical order

\address{Jishan Fan \newline
  Department of Applied Mathematics,
  Nanjing Forestry University,
  Nanjing 210037,  China}
  \email{fanjishan@njfu.edu.cn}

\address{Kyoungsun Kim \newline
  Department of Mathematics, Inha University,
  Incheon 402-751, Korea}
\email{kskim@inha.ac.kr}

\address{Sei Nagayasu \newline
  Department of Mathematical Science,
  Graduate School of Material Science,
  University of Hyogo,
  2167 Shosha, Himeji, Hyogo 671-2280, Japan}
\email{sei@sci.u-hyogo.ac.jp}

\address{Gen Nakamura \newline
  Department of Mathematics, Inha University,
  Incheon 402-751, Korea}
\email{nakamuragenn@gmail.com}

\thanks{Submitted November 11, 2012. Published April 11, 2013.}
\subjclass[2000]{35K10, 35B65}
\keywords{Parabolic equations; discontinuous coefficients;
gradient estimate}

\begin{abstract}
 Li-Vogelius and Li-Nirenberg gave a gradient estimate for solutions
 of strongly elliptic equations and systems of divergence forms
 with piecewise smooth coefficients, respectively. The discontinuities
 of the coefficients are assumed to be given by manifolds of codimension 1,
 which we called them \emph{manifolds of discontinuities}.
 Their gradient estimate is independent of the distances between manifolds
 of discontinuities. In this paper, we gave a parabolic version of their results.
 That is, we gave a gradient estimate for parabolic equations of divergence forms
 with piecewise smooth coefficients. The coefficients are assumed to be
 independent of time and their discontinuities are likewise the previous
 elliptic equations. As an application of this estimate, we also gave a
 pointwise gradient estimate for the fundamental solution of a parabolic operator
 with piecewise smooth coefficients. Both gradient estimates are independent
 of the distances between manifolds of discontinuities.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{example}[theorem]{Example}
\newtheorem{remark}[theorem]{Remark}
\allowdisplaybreaks

\section{Introduction}\label{section:introduction}

For strongly elliptic, second order scalar equations with real coefficients,
it is well known that their solutions have the H{\"o}lder continuity even
in the case that the coefficients are only bounded measurable functions.
However, the solutions do not have the Lipschitz continuity in general.
For example,
Piccinini-Spagnolo~\cite[p.\ 396, Example 1]{PiccininiSpagnolo}
and
Meyers~\cite[p.\ 204]{Meyers}
gave the following case:

\begin{example}[\cite{Meyers,PiccininiSpagnolo}]  \rm
Let $B_1:=\{x \in \mathbb{R}^{n} : | x | < 1\}$
and each $a_{i j} \in L^\infty (B_1)$ be defined as
\[
a_{1 1} = \frac{M x_1^2 + x_2^2}{| x |^2} , \quad
a_{2 2} = \frac{x_1^2 + M x_2^2}{| x |^2} , \quad
a_{1 2} = a_{2 1} = \frac{(M-1) x_1 x_2}{| x |^2}
\]
with a constant $M>1$.
Then, if we define $u$ as
\begin{equation}\label{eq:example of u}
u (x) = \begin{cases}
  | x |^{1 / \sqrt{M}} \frac{x_1}{| x |}
 & \text{if } x \neq 0 , \\
  0 & \text{if } x = 0 ,
 \end{cases}
\end{equation}
it is easy to see that the H{\"o}lder exponent of $u$ is
at least less than or equal to $1 / \sqrt{M}$
(indeed, for $\overline{x} = ( x_1 , 0 )$ we have
\begin{math}
| u( \overline{x} ) - u(0) |= | \overline{x} |^{1 / \sqrt{M}}
\end{math}. Hence we have
\[
\frac{| u( \overline{x} ) - u(0) |}{
 | \overline{x} |^{( 1 / \sqrt{M} ) + \varepsilon}}
= | \overline{x} |^{- \varepsilon}
\to + \infty \quad\text{as } \overline{x} \to 0
\]
for any $\varepsilon > 0$.)
and $u$ satisfies the strongly elliptic scalar equation with real coefficients
\begin{equation}\label{eq:elliptic equation}
\sum_{i,j=1}^2
\frac{\partial}{\partial x_i} \Big(
 a_{i j} \frac{\partial u}{\partial x_{j}}
\Big) = 0.
\end{equation}

The same thing can be said also to the parabolic equation
\begin{equation}\label{eq:parabolic equation}
\frac{\partial u}{\partial t}
- \sum_{i,j=1}^2
\frac{\partial}{\partial x_i} \Big(
 a_{i j} \frac{\partial u}{\partial x_{j}}
\Big) = 0,
\end{equation}
because $u$ given by \eqref{eq:example of u} satisfies this equation.
\end{example}

This example shows that we cannot expect gradient estimates of solutions to
equations \eqref{eq:elliptic equation} and \eqref{eq:parabolic equation}
in the case $a_{ij}\in L^\infty(B_1)$, but we may have the estimates
in the case of piecewise $C^\mu$ (see \eqref{eq:piecewiseCmu} below) coefficients.

The fact that the gradient estimate of solutions is independent of
the distances between manifolds of discontinuities was first observed by
Babu{\v{s}}ka-Andersson-Smith-Levin~\cite{Babuska} numerically
for certain homogeneous isotropic linear systems of elasticity,
that is $| \nabla u |$
is bounded independently of the distances between manifolds of discontinuities.
They considered that this numerical property of solutions is mathematically true.
This is the so-called Babu{\v{s}}ka's conjecture. Recently, proofs for this conjecture
appeared in \cite{LiVogelius} and \cite{LiNirenberg}.
 In elasticity, a small static deformation of an elastic medium with inclusions can be described by an elliptic system of divergence form with piecewise smooth coefficients. The discontinuities of coefficients form the boundaries of inclusions.
Similar physical interpretation is also possible for heat conductors.
Our main theorem \ref{theorem:main} given below ensures that
this property also holds for parabolic equations of the
form \eqref{eq:parabolic equation}. The details of result given
in \cite{LiVogelius} and \cite{LiNirenberg} for scalar equations
 will be given below as Theorem \ref{theorem:LN}.

To state our main theorem, we begin with introducing several notations
 which will be used throughout this paper.
Let $D \subset \mathbb{R}^{n}$ be a bounded domain
with a $C^{1, \alpha}$ boundary for some $0 < \alpha < 1$, which means
that the domain $D$ contains $L$ disjoint subdomains
$D_1 , \ldots , D_{L}$
with $C^{1, \alpha}$ boundaries, i.e.
$D = ( \bigcup_{m=1}^{L} \overline{D_{m}} ) \setminus \partial D$,
and we also assume that
$\overline{D_{m}} \subset D$ for $1 \leq m \leq L-1$. Physically, $D$ is a material and $D_m\,(1\le m\le L-1)$ are
considered as inclusions in $D$.
We define the $C^{1, \alpha}$ norm
(resp.\ $C^{1, \alpha}$ seminorm)
of $C^{1, \alpha}$ domain $D_{m}$
in the same way as in~\cite{LiNirenberg},
that is, as the largest positive number $a$
such that in the $a$-neighborhood of every point of
$\partial D_{m}$, identified as $0$
after a possible translation and rotation
of the coordinates so that $x_n = 0$ is the tangent to
$\partial D_{m}$ at $0$,
$\partial D_{m}$ is given by the graph of a $C^{1, \alpha}$ function $\psi_{m}$,
defined in $| x' | < 2 a$
($x' = ( x_1 , \ldots , x_{n-1} )$),
the $2 a$-neighborhood of $0$ in the tangent plane, and it satisfies the estimate
\begin{math}
\| \psi_{m} \|_{ C^{1, \alpha} ( | x' | < 2 a )
} \leq 1 / a
\end{math}
(resp.\
\begin{math}
[ \psi_{m} ]_{
 C^{1, \alpha} ( | x' | < 2 a )
} \leq 1 / a
\end{math}),
where
\begin{gather*}
[ \psi ]_{
 C^{1, \alpha} ( | x' | < 2 a )}
 := \sup_{
 | x' | , | \xi' | < 2 a} \frac{|
 \nabla' \psi ( x' )
 - \nabla' \psi ( \xi' )
|}{| x' - \xi' |^{\alpha}} ,
 \\
\| \psi \|_{
 C^{1, \alpha} ( | x' | < 2 a )
}
 := \| \psi \|_{C^{1} ( | x' | < 2 a )}
+ [ \psi ]_{
 C^{1, \alpha} ( | x' | < 2 a )
} .
\end{gather*}

Further, let $( a_{i j} )$ be a symmetric, positive definite matrix-valued
function defined on $D$ satisfying
\begin{equation}\label{eq:aijxiixij}
\lambda | \xi |^2
\leq \sum_{i,j=1}^{n} a_{i j} (x) \xi_i \xi_{j}
\leq \Lambda | \xi |^2 .
\end{equation}
Here each $a_{i j}$ is piecewise $C^{\mu}$ in $D$,
$0 < \mu < 1$; that is,
\begin{equation}\label{eq:piecewiseCmu}
a_{i j} (x) = a_{i j}^{(m)} (x) \quad\text{for } x \in D_{m} , \;
1 \leq m \leq L
\end{equation}
with $a_{i j}^{(m)} \in C^{\mu} ( \overline{D_{m}} )$.

As we have already mentioned above, we will discuss in this paper a
gradient estimate for solutions to parabolic equations with piecewise
smooth coefficients.
Our result is a parabolic version for the results of
Li-Vogelius~\cite{LiVogelius} and the scalar equations version
of Li-Nirenberg~\cite{LiNirenberg}.
They showed that solutions $u\in H^1(D)$ to the elliptic equation
\begin{equation}\label{eq:elliptic}
\sum_{i,j=1}^{n} \frac{\partial}{\partial x_i}
\Big( a_{i j} \frac{\partial u}{\partial x_{j}} \Big)
= h + \sum_{i=1}^{n} \frac{\partial g_i}{\partial x_i},
\end{equation}
where $h\in L^\infty(D)$ and each $g_i$ is defined in $D$ such that
 $g_i|_{D_m}$ $(1\le m\le L)$ have continuous extensions
$\in C^\mu(\overline{D_m}),\,0<\mu<1$ up to $\partial D_m$ have global
 $W^{1, \infty}$
and piecewise $C^{1, \alpha'}$ estimates (see \eqref{eq:ellipticest} below).
These estimates are independent of the distances
between inclusions
when a material has inclusions.


We first give the result of Li-Nirenberg~\cite{LiNirenberg}
for scalar equations.


\begin{theorem}[{\cite[Theorem~1.1]{LiNirenberg}}]\label{theorem:LN}
For any $\varepsilon > 0$,
there exists a constant $C_{\sharp} > 0$ such that
for any $\alpha'$ satisfying
\[
0 < \alpha'
< \min \big\{ \mu , \frac{\alpha}{2 ( \alpha + 1 )} \big\} ,
\]
we have
\begin{equation}\label{eq:ellipticest}
\sum_{m=1}^{L} \| u \|_{
 C^{1, \alpha'} ( \overline{D_{m}} \cap D_{\varepsilon} )
} \leq C_{\sharp} \Big(
 \| u \|_{L^2 (D)}
 + \| h \|_{L^{\infty} (D)}
 + \sum_{m=1}^{L} \sum_{i=1}^{n} \| g_i \|_{
  C^{\alpha'} ( \overline{D_{m}} )
 }
\Big) ,
\end{equation}
where we denote
\[
D_{\varepsilon} := \{
 x \in D : \operatorname{dist}( x, \partial D ) > \varepsilon
\}
\]
and a positive constant $C_{\sharp}$ depends only on
\begin{math}
n, L, \mu, \alpha , \varepsilon , \lambda , \Lambda ,
\| a_{i j} \|_{C^{\alpha'} ( \overline{D_{m}})}
\end{math}
and the $C^{1, \alpha'}$ norms of $D_{m}$.
\end{theorem}

\begin{remark} \rm
The constant $C_{\sharp} > 0$
is independent of the distances between inclusions $D_{m}$.
Therefore, the estimate \eqref{eq:ellipticest} holds
even in the case that some of inclusions touch another inclusions
as in Figure~\ref{figure:1}.
\end{remark}

\begin{figure}[ht]
\begin{center}
\begin{picture}(240,120)
%%% D
\qbezier(48,15.6024)(96,3.6)(176,3.6)
\qbezier(176,3.6)(240,3.6)(240,45.5976)
\qbezier(240,45.5976)(240,75.6)(192,93.6)
\qbezier(192,93.6)(144,111.6)(80,111.6)
\qbezier(80,111.6)(0,111.6)(0,69.6024)
\qbezier(0,69.6024)(0,27.5976)(48,15.6024)
%%% D1
\qbezier(34.6,47.8012)(44.2,41.8)(60.2,41.8)
\qbezier(60.2,41.8)(73,41.8)(73,62.7988)
\qbezier(73,62.7988)(73,77.8)(63.4,86.8)
\qbezier(63.4,86.8)(53.8,95.8)(41,95.8)
\qbezier(41,95.8)(25,95.8)(25,74.8012)
\qbezier(25,74.8012)(25,53.7988)(34.6,47.8012)
\put(49,64){\makebox(0,0){$D_1$}}
%%% D2
\qbezier(102.72,38.24096)(95.04,33.44)(82.24,33.44)
\qbezier(82.24,33.44)(72,33.44)(72,50.23904)
\qbezier(72,50.23904)(72,62.24)(79.68,69.44)
\qbezier(79.68,69.44)(87.36,76.64)(97.6,76.64)
\qbezier(97.6,76.64)(110.4,76.64)(110.4,59.84096)
\qbezier(110.4,59.84096)(110.4,43.03904)(102.72,38.24096)
\put(92,56){\makebox(0,0){$D_2$}}
%%% D3
\qbezier(158.72,73.75904)(151.04,78.56)(138.24,78.56)
\qbezier(138.24,78.56)(128,78.56)(128,61.76096)
\qbezier(128,61.76096)(128,49.76)(135.68,42.56)
\qbezier(135.68,42.56)(143.36,35.36)(153.6,35.36)
\qbezier(153.6,35.36)(166.4,35.36)(166.4,52.15904)
\qbezier(166.4,52.15904)(166.4,68.96096)(158.72,73.75904)
\put(148,44){\makebox(0,0){$D_{3}$}}
%%% D4
\qbezier(154.997793341,67.4460357714)(153.9798989874,71.8587214973)(149.4544155878,76.3842048969)
\qbezier(149.4544155878,76.3842048969)(145.8340288681,80.0045916166)(139.8946713174,74.0652340658)
\qbezier(139.8946713174,74.0652340658)(135.651691219,69.8222539674)(135.8213968465,64.5613795154)
\qbezier(135.8213968465,64.5613795154)(135.9911024739,59.3005050634)(139.6114891936,55.6801183437)
\qbezier(139.6114891936,55.6801183437)(144.1369725932,51.1546349441)(150.0763301439,57.0939924948)
\qbezier(150.0763301439,57.0939924948)(156.0167059284,63.0343682793)(154.997793341,67.4460357714)
\put(144,64){\makebox(0,0){$D_{4}$}}
%%% D5
\qbezier(206.976,74.303616)(200.832,76.224)(190.592,76.224)
\qbezier(190.592,76.224)(182.4,76.224)(182.4,69.504384)
\qbezier(182.4,69.504384)(182.4,64.704)(188.544,61.824)
\qbezier(188.544,61.824)(194.688,58.944)(202.88,58.944)
\qbezier(202.88,58.944)(213.12,58.944)(213.12,65.663616)
\qbezier(213.12,65.663616)(213.12,72.384384)(206.976,74.303616)
\put(200,66.4){\makebox(0,0){$D_{5}$}}
%%% D6
\qbezier(206.976,41.1643392)(200.832,43.4688)(190.592,43.4688)
\qbezier(190.592,43.4688)(182.4,43.4688)(182.4,35.4052608)
\qbezier(182.4,35.4052608)(182.4,29.6448)(188.544,26.1888)
\qbezier(188.544,26.1888)(194.688,22.7328)(202.88,22.7328)
\qbezier(202.88,22.7328)(213.12,22.7328)(213.12,30.7963392)
\qbezier(213.12,30.7963392)(213.12,38.8612608)(206.976,41.1643392)
\put(200,31.68){\makebox(0,0){$D_{6}$}}
\put(148,20){\makebox(0,0){$D_{7}$}}
\end{picture}
\end{center}
\caption{The case that an inclusion touches another inclusion.
($L=7$)} \label{figure:1}
\end{figure}

Now, we consider the parabolic equation
\begin{equation}\label{eq:parabolic}
\frac{\partial u}{\partial t}
- \sum_{i,j=1}^{n} \frac{\partial}{\partial x_i}
\Big( a_{i j} \frac{\partial u}{\partial x_{j}} \Big)
= f - \sum_{i=1}^{n} \frac{\partial f_i}{\partial x_i}
\quad\text{in } Q,
\end{equation}
where
\begin{gather*}
 f\in L^\infty(Q),\,\frac{\partial f}{\partial t}\in L^\kappa(Q),
 \\
 f_i\in L^p(Q),\,\frac{\partial f_i}{\partial t}\in L^p(Q) \quad\text{and}\quad
f_i = f_i^{(m)} \quad\text{on } D_{m} \times (0,T],
\end{gather*}
with $p>n+2$, $\kappa = p(n+2)/(n+2+p)$,
$Q := D \times (0,T]$,
$f_i^{(m)} \in L^{\infty} ( 0, T; C^{\mu} ( \overline{D_{m}} ) ) $.

Now we define a weak solution to the equation \eqref{eq:parabolic}.

\begin{definition} \rm
We call
\begin{math}
u \in V_2^{1,0} (Q)
:= L^2 ( 0,T; H^{1} (D) )
\cap C ( [0,T]; L^2 (D) )
\end{math}
a weak solution to the equation \eqref{eq:parabolic}
when
\begin{align}
& \int_{D} u( x, t' ) \, \varphi ( x, t' ) \, d x
- \int_{0}^{t'} \int_{D}
 u( x, t ) \, \frac{\partial \varphi}{\partial t} ( x, t ) \,
d x \, d t  \\
& \text{}
+ \int_{0}^{t'} \int_{D}
 \sum_{i,j=1}^{n} a_{i j} (x) \,
 \frac{\partial u}{\partial x_{j}} (x,t) \,
 \frac{\partial \varphi}{\partial x_i} (x,t) \,
d x \, d t  \\
& = \int_{0}^{t'} \int_{D}
 f(x,t) \, \varphi (x,t) \,
d x \, d t
+ \int_{0}^{t'} \int_{D}
 \sum_{i=1}^{n} f_i (x,t) \,
 \frac{\partial \varphi}{\partial x_i} (x,t) \,
d x \, d t \label{eq:weaksolution}
\end{align}
for any
\begin{math}
\varphi \in L^2 ( 0, T; {H}_{0}^{1} (D) )
\cap H^{1} ( 0, T; L^2 (D) )
\end{math}
with $\varphi ( \cdot , 0 ) = 0$
and $0 < t' \leq T$,
where $H_{0}^{1} (D)$ is the usual $L^2$-Sobolev space
with supports in $\overline D$.
\end{definition}

Our main result is as follows.

\begin{theorem}[Main theorem]\label{theorem:main}
Any weak solutions $u \in V_2^{1,0} (Q)$
to {\rm \eqref{eq:parabolic}}
have the following up to the inclusion boundary regularity estimate:
For any $\varepsilon > 0$,
there exists a constant $C_{\sharp}' > 0$
such that for any $\alpha'$ satisfying
\begin{equation}\label{eq:alphaprime}
0 < \alpha'
< \min \big\{ \mu , \frac{\alpha}{2 ( \alpha + 1 )} \big\} ,
\end{equation}
we have
\[
\sum_{m=1}^{L}
\sup_{\varepsilon^2 < t \leq T}
\| u ( \cdot , t ) \|_{
 C^{1, \alpha'} ( \overline{D_{m}} \cap D_{\varepsilon} )
} \leq C_{\sharp}'
\left(
 \| u \|_{L^2 (Q)} + F_{\ast} + F_{\ast \ast}
\right) ,
\]
where
\begin{gather*}
F_{\ast}
:= \| f \|_{L^{\kappa} (Q)}
+ \| f \|_{
 L^{\max \{ 2, \kappa \}} (Q)
} + \| f \|_{L^{\infty} (Q)}
+ \|
 \frac{\partial f}{\partial t}
\|_{L^{\kappa} (Q)} ,  \\
F_{\ast \ast}
 := \sum_{i=1}^{n} \Bigl(
 \| f_i \|_{L^{p} (Q)}
 + \|
  \frac{\partial f_i}{\partial t}
 \|_{L^2 (Q)}
 + \|
  \frac{\partial f_i}{\partial t}
 \|_{L^{p} (Q)}
   + \sum_{m=1}^{L} \sup_{0 < t \leq T}
 \| f_i ( \cdot , t ) \|_{
  C^{\alpha'} ( \overline{D_{m}} )
 }
\Bigr)
\end{gather*}
and
$C_{\sharp}'$ depends only on
\begin{math}
n, L, \mu, \alpha , \varepsilon , \lambda , \Lambda , p,
\| a_{i j} \|_{C^{\alpha'} ( \overline{D_{m}})}
\end{math}
and the $C^{1, \alpha'}$ norms of $D_{m}$.
\end{theorem}

\begin{remark}\label{remark:main} \rm
(i)
Again,
the constant $C_{\sharp}' > 0$
is independent of the distances between inclusions $D_{m}$.
Then Theorem~\ref{theorem:main} holds
even in the case that an inclusion touches another inclusion
as Figure~\ref{figure:1}.

(ii)
It is easy to obtain
\begin{gather*}
F_{\ast}
 \leq C^{\ast}
\Big(
 \| f \|_{L^{\infty} (Q)}
 + \|
  \frac{\partial f}{\partial t}
 \|_{L^{\kappa} (Q)}
\Big) ,  \\
F_{\ast \ast}
 \leq C^{\ast} \sum_{i=1}^{n} \Big(
 \sum_{m=1}^{L} \sup_{0 < t \leq T}
 \| f_i ( \cdot , t ) \|_{
  C^{\alpha'} ( \overline{D_{m}} )
 }
 + \|
  \frac{\partial f_i}{\partial t}
 \|_{L^{p} (Q)}
\Big) .
\end{gather*}
However, a constant $C^{\ast} > 0$
depends on $T$ and $D$, unfortunately.
\end{remark}

For heat conductive materials with inclusions,
the solution $u$ of the initial boundary value problem for
\eqref{eq:parabolic} with heat flux given
on $\partial D \times (0,T]$
and zero initial temperature at $t = 0$ in $D$
describes the temperature distribution in $D$.
Injecting the heat flux and measuring the temperature distribution on
$\partial D \times (0,T]$ is the measurement of the so called
active thermography.
This is an non-destructive testing to identify unknown cracks,
cavities and inclusions inside a heat conductor from the measurement.
As a mathematical study of the active thermography,
a method called the dynamical probe method has been given (\cite{IKN}).
It can approximately identify for instance inclusions by one
measurement.
For identifying the inclusions precisely, it needs infinitely
many measurements. Further,
it uses the gradient estimate of the fundamental solution
of the heat equation with discontinuous conductivities.

The dynamical probe method has been developed only for the case that
 the inclusions do not touch another inclusions. So, it is interesting
to consider the case when some of them touch. For the first task to handle
this case, we need to have the gradient estimate of the fundamental solution.
Our main result has given an answer to this. Similar situation can be
considered for active thermography and
non-destructive testing using acoustic waves. For example, \cite{NUW}
and \cite{Yoshida}
effectively used a result of Li-Vogelius~\cite{LiVogelius} to give a
procedure of enclosing the inclusions
by the enclosure method (see \cite{Ikehata}, for example). What is interested about their arguments is that, by adding further arguments, we can
even enclose the inclusions in the case that they can touch another
inclusions \cite{Nag-Nak}.
More precisely for an increasing sets of inclusions, these inclusions can touch
at a point of the boundary of the largest inclusion.
Therefore, we believe that
our gradient estimates will be useful for inverse problems identifying
 unknown inclusions.

This kind of gradient estimates stated in Theorem~\ref{theorem:main}
for solutions of parabolic equations was initiated by Li~\cite{HGLi}
in his doctor thesis written in Chinese and was completed recently
in Li-Li~\cite{LiLi}. In~\cite{LiLi}, they even discussed the interior
gradient estimates of solutions of a second order parabolic system of
divergence form with inclusions which can touch another inclusions.
They also allowed that the coefficients can depend on the time.
However, it should be noted that they could not allow the inclusions
to depend on the time. Hence, showing the interior gradient estimate
for this case is still opened.

We have to emphasize the following two things. The one is that we
independently obtained our results. After we finished our paper and
 posted our paper in the preprint server arXiv,
Li sent the paper~\cite{LiLi} to one of us. But we still did not
know the paper~\cite{HGLi} until very recently by a chance.
The results of~\cite{LiLi} was also posted in arXiv after us.

The another is that our proof of Theorem~\ref{theorem:main}
is totally different from the proofs given in \cite{HGLi} and \cite{LiLi}.
We first reduce the problem of the interior gradient estimates
for solutions of parabolic equations to that of elliptic parts of the equations
using by the idea in~\cite{LRU}, and directly apply the result~\cite{LiNirenberg}
for elliptic equations.
By this method, we estimate not only
the $L^{\infty}$-norm of $\nabla u$ but also
\begin{math}
\sup_{\varepsilon^2 < t \leq T}
\| u( \cdot , t ) \|_{
 C^{1, \alpha'} ( \overline{D_{m}} \cap D_{\varepsilon} )
}
\end{math}
more easily. We also remark that the constant
$C_{\sharp}'$
in Theorem~\ref{theorem:main}
is independent of $T$.
In \cite{HGLi},
this property is not obvious at least
even for the case when the right-hand side
$f$ and $f_i$ on \eqref{eq:parabolic}
are identically equal zero.

The rest of this paper is organized as follows.
In Section~\ref{section:proof},
we prove our main theorem, i.e.\ Theorem~\ref{theorem:main}
by applying Lemma~\ref{lemma:keyestimate}.
We prove Lemma~\ref{lemma:keyestimate}
in Section~\ref{section:someestimates}.
In Section~\ref{section:gradestoffundsol},
we consider a pointwise gradient estimate
for the fundamental solution
of parabolic operators
with piecewise smooth coefficients
by applying Theorem~\ref{theorem:main}.

\section{Proof of main result}\label{section:proof}

In this section, we prove our main theorem.
We first state some estimates in Lemma~\ref{lemma:keyestimate}
which we need to prove our main theorem.
We prove Lemma~\ref{lemma:keyestimate}
in Section~\ref{section:someestimates}.

\begin{lemma}\label{lemma:keyestimate}
Let $( a_{i j} )$ be a matrix-valued function defined on $D$.
Assume that $( a_{i j} )$ is symmetric, positive definite,
and satisfies the condition {\rm \eqref{eq:aijxiixij}}.
Let $Q$ as before and
\begin{math}
\widehat{Q}_{\varepsilon}
:= D_{\varepsilon} \times ( \varepsilon^2 , T ]
\end{math}.
Then for $p>n+2$,
a weak solution $u \in V_2^{1,0} (Q)$
to {\rm \eqref{eq:parabolic}} satisfies
the following estimates:
\begin{align}
\sup_{\varepsilon^2 < t \leq T}
\| u ( \cdot , t ) \|_{L^2 ( D_{\varepsilon} )}
& \leq C \left( \| u \|_{L^2 (Q)} + F_{0} \right) ,
\label{eq:LinftyL2}  \\
\| u \|_{L^{\infty} ( \widehat{Q}_{\varepsilon} )}
& \leq C \left( \| u \|_{L^2 (Q)} + F_{0} \right) ,
\label{eq:Linfty}  \\
\| \frac{\partial u}{\partial t} \|_{
 L^2 ( \widehat{Q}_{\varepsilon})
}
& \leq C \left( \| u \|_{L^2 (Q)} + F_1 \right) ,
\label{eq:dudtL2uL2}
\end{align}
where we set
\begin{align}
F_{0} & :=
\| f \|_{L^{\frac{p(n+2)}{n+2+p}} (Q)}
+ \sum_{i=1}^{n} \| f_i \|_{L^{p} (Q)} ,
\label{eq:defOfF0}  \\
F_1 & :=
\| f \|_{
 L^{\max \{ 2, \frac{p(n+2)}{n+2+p} \}} (Q)
}
+ \sum_{i=1}^{n} \Big(
 \| f_i \|_{L^{p} (Q)}
 + \|
  \frac{\partial f_i}{\partial t}
 \|_{L^2 (Q)}
\Big) , \label{eq:defOfF1}
\end{align}
and $C>0$ depends only on $n, \lambda , \Lambda , p$
and $\varepsilon$.
\end{lemma}

Now we prove our main theorem by applying
Lemma~\ref{lemma:keyestimate}.
This proof is inspired by \cite{LRU}.

\begin{proof}[Proof of Theorem~\ref{theorem:main}]
Before going into the proof, we remark that a general constant $C$ which
we used below in our estimates depends only on $n, \lambda , \Lambda , p$
and $\varepsilon_{j}$ ($j = 1, 2, 3$). To begin with the proof,
let $0 < \varepsilon_1 < \varepsilon_2 < \varepsilon_{3}$.
Then we have
\begin{equation}\label{eq:proof11}
\sup_{\varepsilon_2^2 < t \leq T}
\| u ( \cdot , t ) \|_{L^2 ( D_{\varepsilon_2} )}
\leq C \left(
 \| u \|_{L^2 (Q)} + F_{0}
\right)
\end{equation}
and
\begin{equation}\label{eq:proof12}
\|
 \frac{\partial u}{\partial t}
\|_{L^2 ( \widehat{Q}_{\varepsilon_1} )}
\leq C \left(
 \| u \|_{L^2 (Q)} + F_1
\right)
\end{equation}
by \eqref{eq:LinftyL2} and \eqref{eq:dudtL2uL2}
in Lemma~\ref{lemma:keyestimate},
where $F_{0}$, $F_1$ are defined by
\eqref{eq:defOfF0} and \eqref{eq:defOfF1}.
On the other hand,
$u_{t} = \partial u / \partial t$
satisfies the equation
\[
\frac{\partial u_{t}}{\partial t}
- \sum_{i,j=1}^{n} \frac{\partial}{\partial x_i} \Big(
 a_{i j} (x) \frac{\partial u_{t}}{\partial x_{j}}
\Big) = \frac{\partial f}{\partial t}
- \sum_{i=1}^{n} \frac{\partial}{\partial x_i} \Big(
 \frac{\partial f_i}{\partial t}
\Big)
\]
by applying $\partial / \partial t$ to \eqref{eq:parabolic}
(also see Remark~\ref{remark:steklov}).
Hence we have
\begin{equation}\label{eq:proof2}
\| u_{t} \|_{L^{\infty} ( \widehat{Q}_{\varepsilon_2} )}
\leq C \left(
 \| u_{t} \|_{L^2 ( \widehat{Q}_{\varepsilon_1} )}
 + F_{0}'
\right)
\end{equation}
by Lemma~\ref{lemma:keyestimate} \eqref{eq:Linfty},
where we define
\[
F_{0}' :=
\|
 \frac{\partial f}{\partial t}
\|_{L^{\frac{p(n+2)}{n+2+p}} (Q)}
+ \sum_{i=1}^{n}
\|
 \frac{\partial f_i}{\partial t}
\|_{L^{p} (Q)} .
\]
In particular,
$u_t( \cdot , t ) \in L^{\infty} ( D_{\varepsilon_2} )$
holds for a.e.\ $t \in ( \varepsilon_2^2 , T ]$.
Now we regard the equation \eqref{eq:parabolic} as
the elliptic equation
\begin{equation}\label{eq:aselliptic}
\sum_{i,j=1}^{n} \frac{\partial}{\partial x_i} \Big(
 a_{i j} (x) \frac{\partial u}{\partial x_{j}}
\Big)
= \frac{\partial u}{\partial t}
- f + \sum_{i=1}^{n} \frac{\partial f_i}{\partial x_i}
\end{equation}
by fixing $t \in ( \varepsilon_2^2 , T ]$.
We remark that
$\partial u / \partial t - f \in L^{\infty} ( D_{\varepsilon_2} )$.
Then, for any $\alpha'$
with the condition \eqref{eq:alphaprime},
we have the estimate
\begin{equation}
\begin{aligned}
\sum_{m=1}^{L}
\| u ( \cdot , t ) \|_{
 C^{1, \alpha'}
 ( \overline{D_{m}} \cap D_{\varepsilon_{3}} )
}
& \leq C_{\sharp} \Bigl(
 \| u ( \cdot , t ) \|_{L^2 ( D_{\varepsilon_2} )}
 + \| \frac{\partial u}{\partial t} ( \cdot , t )
 \|_{
  L^{\infty} ( D_{\varepsilon_2} )
 }\\
&\quad
 + \| f ( \cdot , t ) \|_{L^{\infty} ( D_{\varepsilon_2} )}
 + \sum_{m=1}^{L} \sum_{i=1}^{n}
 \| f_i ( \cdot , t ) \|_{
  C^{\alpha'} ( \overline{D_{m}})
 }
\Bigr)
\end{aligned} \label{eq:x}
\end{equation}
by Theorem~\ref{theorem:LN}, where
$C_{\sharp} > 0$ depends only on
\begin{math}
n, L, \mu, \alpha , \varepsilon , \lambda , \Lambda ,
\| a_{i j} \|_{C^{\alpha'} ( \overline{D_{m}})}
\end{math}
and the $C^{1, \alpha'}$ norms of $D_{m}$.
Taking the supremum of the inequality
\eqref{eq:x} over $( \varepsilon_2^2 , T ]$
with respect to $t$,
and using \eqref{eq:proof11}, \eqref{eq:proof12}
and \eqref{eq:proof2},
we have
\begin{align*}
& \sum_{m=1}^{L}
\sup_{\varepsilon_2^2 < t \leq T}
\| u ( \cdot , t ) \|_{
 C^{1, \alpha'}
 ( \overline{D_{m}} \cap D_{\varepsilon_{3}} )
} \\
& \leq C_{\sharp} \Bigl(
 \sup_{\varepsilon_2^2 < t \leq T}
 \| u ( \cdot , t ) \|_{L^2 ( D_{\varepsilon_2} )}
 + \| \frac{\partial u}{\partial t}
 \|_{
  L^{\infty} ( \widehat{Q}_{\varepsilon_2} )
 }
 + \| f \|_{L^{\infty} ( \widehat{Q}_{\varepsilon_2} )} \\
 & \quad
 + \sum_{m=1}^{L} \sum_{i=1}^{n}
 \sup_{\varepsilon_2^2 < t \leq T}
 \| f_i ( \cdot , t ) \|_{
  C^{\alpha'} ( \overline{D_{m}})
 }
\Bigr)  \\
& \leq C_{\sharp} C \Bigl(
 \| u \|_{L^2 (Q)}
 + F_{0} + F_1 + F_{0}'
 + \| f \|_{L^{\infty} ( \widehat{Q}_{\varepsilon_2} )} \\
 & \quad
 + \sum_{m=1}^{L} \sum_{i=1}^{n}
 \sup_{\varepsilon_2^2 < t \leq T}
 \| f_i ( \cdot , t ) \|_{
  C^{\alpha'} ( \overline{D_{m}})
 }
\Bigr) ,
\end{align*}
which is the estimate we want to obtain.
\end{proof}

\begin{remark}\label{remark:steklov} \rm
Since we assume that $u$ belongs only in $V_2^{1,0} (Q)$
with respect to the regularity of a weak solution,
one may think that we cannot apply $\partial / \partial t$ directly.
However, it is enough to consider the Steklov mean function
and to make $h$ tend to $0$,
where we define the Steklov mean function
$v_{h}$ of $v$ by
\[
v_{h} (x,t) = \frac{1}{h} \int_{t}^{t + h} v( x, \tau ) \, d \tau.
\]
Hereafter we omit the detail with respect to this remark
although we often apply this argument. Also see
\cite[III \S2 p.\,141]{LSU}
and (62) in \cite[p.\,152]{LRU}, for example.
\end{remark}

\section{Some estimates}\label{section:someestimates}

In this section, we prove Lemma~\ref{lemma:keyestimate}.
The estimates \eqref{eq:LinftyL2} and \eqref{eq:Linfty} are
well known, but we give these proofs in Appendix
for readers' convenience.
To show the estimate \eqref{eq:dudtL2uL2}, we prepare some necessary lemmas
for its proof.

Throughout this section, $C > 0$ denotes a general constant
depending only on $n, \lambda , \Lambda$. Also, we assume that
the coefficient $( a_{i j} )$ is a matrix-value function
defined on $D$, symmetric, positive definite,
and satisfies the condition \eqref{eq:aijxiixij}.
Moreover, we set
$Q_{r} := B_{r} ( x_{0} ) \times ( t_{0} - r^2 , \, t_{0} ]$,
and assume that
$Q_{2 \rho} \subset D \times (0,T]$
with $0 < \rho \leq 1$.


The following two lemmas are essentially shown in \cite{LRU}.
We give their proofs here for the sake of completeness.

\begin{lemma}[{\cite[Lemma~3]{LRU}}]\label{lemma:DuL2}
Let $1 < r < \infty$ and $1/r + 1 / r' = 1$.
Then a solution $u$ to {\rm \eqref{eq:parabolic}}
satisfies the estimate
\begin{equation}\label{eq:DuL2}
\| \nabla u \|_{L^2 ( Q_{\rho} )}
\leq C \Big[
 ( \rho^{n/2} + \rho^{(n+2) / r'} ) \operatorname{osc}_{Q_{2 \rho}} u
 + \| f \|_{L^{r} ( Q_{2 \rho} )}
 + \sum_{i=1}^{n} \| f_i \|_{L^2 ( Q_{2 \rho} )}
\Big] .
\end{equation}
\end{lemma}

\begin{proof}
Let $\zeta$ be a smooth cut-off function on $Q_{2 \rho}$
satisfying $\zeta \equiv 1$ on $Q_{\rho}$,
$\zeta \equiv 0$ on
$Q_{2 \rho} \setminus Q_{3 \rho / 2}$,
$0 \leq \zeta \leq 1$ on $Q_{2 \rho}$,
and
\begin{math}
| \partial \zeta / \partial t |
+ | \nabla \zeta |^2
\leq C \rho^{-2}
\end{math}
on $Q_{2 \rho}$.
Let $u_{0}$ be the average value of $u$ in $Q_{2 \rho}$:
\[
u_{0} := \frac{1}{| Q_{2 \rho} |}
\iint_{Q_{2 \rho}} u(x,t) \, d x \, d t,
\]
where $| Q_{2 \rho} |$ denotes the measure of $ Q_{2 \rho}$.
Testing \eqref{eq:parabolic}
by $( u - u_{0} ) \zeta^2$
and integrating by parts
(i.e.\ taking $\varphi = ( u - u_{0} ) \zeta^2$
for \eqref{eq:weaksolution}.
Also see Remark~\ref{remark:steklov}), we have
\begin{align*}
& \frac{1}{2} \int_{B_{2 \rho} ( x_{0} )}
 \bigl( ( u - u_{0} )^2 \zeta^2 \bigr) ( x, t_{0} ) \,
d x
- \iint_{Q_{2 \rho}}
 ( u - u_{0} )^2 \zeta \frac{\partial \zeta}{\partial t}
d x \, d t \\
& \text{}
+ \iint_{Q_{2 \rho}}
 \sum_{i,j=1}^{n} a_{i j} \frac{\partial u}{\partial x_{j}}
 \frac{\partial u}{\partial x_i} \zeta^2 \,
d x \, d t
+ 2 \iint_{Q_{2 \rho}}
 \sum_{i,j=1}^{n} a_{i j} \frac{\partial u}{\partial x_{j}}
 ( u - u_{0} ) \zeta \frac{\partial \zeta}{\partial x_i}
d x \, d t \\
& = \iint_{Q_{2 \rho}} f ( u - u_{0} ) \zeta^2 \, d x \, d t
+ \sum_{i=1}^{n} \iint_{Q_{2 \rho}}
 \Big[
  f_i \frac{\partial u}{\partial x_i} \zeta^2
  + 2 f_i ( u - u_{0} ) \zeta \frac{\partial \zeta}{\partial x_i}
 \Big]
d x \, d t.
\end{align*}
Hence we have
\begin{align*}
& \frac{1}{2} \int_{B_{2 \rho} ( x_{0} )}
 \bigl( ( u - u_{0} )^2 \zeta^2 \bigr) ( x, t_{0} ) \,
d x + \lambda \iint_{Q_{2 \rho}}
 | \nabla u |^2 \zeta^2 \,
d x \, d t \\
& \leq \frac{1}{2} \int_{B_{2 \rho} ( x_{0} )}
 \bigl( ( u - u_{0} )^2 \zeta^2 \bigr) ( x, t_{0} ) \,
d x + \iint_{Q_{2 \rho}}
 \sum_{i,j=1}^{n} a_{i j} \frac{\partial u}{\partial x_{j}}
 \frac{\partial u}{\partial x_i} \zeta^2 \,
d x \, d t  \\
& = \iint_{Q_{2 \rho}}
 ( u - u_{0} )^2 \zeta \frac{\partial \zeta}{\partial t}
d x \, d t
- 2 \iint_{Q_{2 \rho}}
 \sum_{i,j=1}^{n} a_{i j} \frac{\partial u}{\partial x_{j}}
 ( u - u_{0} ) \zeta \frac{\partial \zeta}{\partial x_i}
d x \, d t \\
& \hspace*{3ex} \text{}
+ \iint_{Q_{2 \rho}} f ( u - u_{0} ) \zeta^2 \, d x \, d t \\
& \hspace*{3ex} \text{}
+ \sum_{i=1}^{n} \iint_{Q_{2 \rho}}
 \Big[
  f_i \frac{\partial u}{\partial x_i} \zeta^2
  + 2 f_i ( u - u_{0} ) \zeta \frac{\partial \zeta}{\partial x_i}
 \Big]
d x \, d t  \\
& \leq \iint_{Q_{2 \rho}}
 ( u - u_{0} )^2 \zeta
 \big| \frac{\partial \zeta}{\partial t} \big|
d x \, d t
+ \varepsilon_1 \iint_{Q_{2 \rho}}
 | \nabla u |^2 \zeta^2 \,
d x \, d t \\
& \hspace*{3ex} \text{}
+ \frac{C}{\varepsilon_1} \iint_{Q_{2 \rho}}
 | u - u_{0} |^2 | \nabla \zeta |^2 \,
d x \, d t
+ \frac{1}{2} \Big(
 \iint_{Q_{2 \rho}}
  | f \zeta |^{r} \,
 d x \, d t
\Big)^{2/r} \\
& \hspace*{3ex} \text{}
+ \frac{1}{2} \Big(
 \iint_{Q_{2 \rho}}
  | ( u - u_{0} ) \zeta |^{r'} \,
d x \, d t
\Big)^{2 / r'}
+ \varepsilon_1 \iint_{Q_{2 \rho}}
 | \nabla u |^2 \zeta^2 \,
d x \, d t \\
& \hspace*{3ex} \text{}
+ \big( \frac{1}{\varepsilon_1} + 1 \big)
\iint_{Q_{2 \rho}}
 \sum_{i=1}^{n} | f_i |^2 \zeta^2 \,
d x \, d t
+ \iint_{Q_{2 \rho}}
 | u - u_{0} |^2
 | \nabla \zeta |^2 \,
d x \, d t .
\end{align*}
We now take $\varepsilon_1 > 0$ small enough.
Then, we have
\begin{align*}
&\iint_{Q_{\rho}} | \nabla u |^2 \, d x \, d t\\
&\leq \iint_{Q_{2 \rho}}
 | \nabla u |^2 \zeta^2 \,
d x \, d t  \\
& \leq C \iint_{Q_{2 \rho}}
 ( u - u_{0} )^2 \Big[
  \zeta | \frac{\partial \zeta}{\partial t} |
  + | \nabla \zeta |^2
 \Big]
d x \, d t
+ C \Big(
 \iint_{Q_{2 \rho}}
  | ( u - u_{0} ) \zeta |^{r'} \,
 d x \, d t
\Big)^{2 / r'} \\
& \quad
+ C \Big(
 \iint_{Q_{2 \rho}} | f \zeta |^{r} \, d x \, d t
\Big)^{2/r}
+ C \iint_{Q_{2 \rho}}
 \sum_{i=1}^{n} | f_i |^2 \zeta^2 \,
d x \, d t  \\
& \leq C \Big[
 \big( \rho^{n} + \rho^{2(n+2) / r'} \big)
 \big( \operatorname{osc}_{Q_{2 \rho}} u \big)^2
 + \| f \|_{L^{r} ( Q_{2 \rho} )}^2
 + \sum_{i=1}^{n} \| f_i \|_{L^2 ( Q_{2 \rho} )}^2
\Big],
\end{align*}
because
$| u(x,t) - u_{0} | \leq \operatorname{osc}_{Q_{2 \rho}} u$
holds for any $(x,t) \in Q_{2 \rho}$. This completes the proof.
\end{proof}

\begin{lemma}[{\cite[Lemma~5]{LRU}}]\label{lemma:dudtL2}
A solution $u$ to  \eqref{eq:parabolic} satisfies the estimate
\begin{equation}
\begin{aligned}
\|
 \frac{\partial u}{\partial t}
\|_{L^2 ( Q_{\rho} )}
& \leq C \Bigl[
 \rho^{-1} \| \nabla u \|_{L^2 ( Q_{2 \rho} )}
 + \| f \|_{L^2 ( Q_{2 \rho} )}  \\
 & \quad  + \sum_{i=1}^{n}
 \Big(
  \rho^{-1}
  \| f_i \|_{L^2 ( Q_{2 \rho} )}
  + \|
   \frac{\partial f_i}{\partial t}
  \|_{L^2 ( Q_{2 \rho} )}
 \Big)
\Bigr]
\end{aligned} \label{eq:dudtL2}
\end{equation}
\end{lemma}

\begin{proof}
We first take the same smooth cut-off function $\zeta$
as in the proof of Lemma~\ref{lemma:DuL2}.
Testing \eqref{eq:parabolic} by
$( \partial u / \partial t ) \zeta^2$
and integrating by parts (also see Remark~\ref{remark:steklov}),
we have
\begin{align*}
& \frac{1}{2} \int_{B_{2 \rho} ( x_{0} )}
 \sum_{i,j=1}^{n}
 \Big(
  a_{i j} \frac{\partial u}{\partial x_i}
  \frac{\partial u}{\partial x_{j}} \zeta^2
 \Big) ( x, t_{0} ) \,
d x \\
& + \iint_{Q_{2 \rho}}
 \Big[
  | \frac{\partial u}{\partial t}|^2
  \zeta^2
  - \sum_{i,j=1}^{n} a_{i j} \frac{\partial u}{\partial x_i}
  \frac{\partial u}{\partial x_{j}}
  \zeta \frac{\partial \zeta}{\partial t}
  + 2 \sum_{i,j=1}^{n} a_{i j} \frac{\partial u}{\partial x_{j}}
  \frac{\partial u}{\partial t} \zeta
  \frac{\partial \zeta}{\partial x_i}
 \Big]
d x \, d t \\
& = \iint_{Q_{2 \rho}}
 f \frac{\partial u}{\partial t} \zeta^2 \,
d x \, d t
+ \sum_{i=1}^{n} \Bigl[
 \int_{B_{2 \rho} ( x_{0} )}
  \Big( f_i \frac{\partial u}{\partial x_i} \zeta^2 \Big)
  ( x, t_{0} ) \,
 d x \\
 & \quad  + \iint_{Q_{2 \rho}}
  \Big(
   - \frac{\partial f_i}{\partial t}
   \frac{\partial u}{\partial x_i} \zeta^2
   - 2 f_i \frac{\partial u}{\partial x_i}
   \zeta \frac{\partial \zeta}{\partial t}
   + 2 f_i \frac{\partial u}{\partial t}
   \zeta \frac{\partial \zeta}{\partial x_i}
  \Big)
 \Bigr]
d x \, d t
\end{align*}
due to
\[
\sum_{i,j=1}^{n} a_{i j}
\frac{\partial^2 u}{\partial t \partial x_i}
\frac{\partial u}{\partial x_{j}} \zeta^2
= \frac{1}{2} \frac{\partial}{\partial t} \Big(
 \sum_{i,j=1}^{n} a_{i j}
 \frac{\partial u}{\partial x_i}
 \frac{\partial u}{\partial x_{j}} \zeta^2
\Big)
- \sum_{i.j=1}^{n} a_{i j}
\frac{\partial u}{\partial x_i} \frac{\partial u}{\partial x_{j}}
\zeta \frac{\partial \zeta}{\partial t}
\]
and
\[
f_i \frac{\partial^2 u}{\partial t \partial x_i} \zeta^2
= \frac{\partial}{\partial t} \left(
 f_i \frac{\partial u}{\partial x_i} \zeta^2
\right)
- \frac{\partial u}{\partial x_i}
\frac{\partial}{\partial t} ( f_i \zeta^2 ) .
\]
Hence we have
\begin{align*}
& \frac{\lambda}{2} \int_{B_{2 \rho} ( x_{0} )}
 \left( | \nabla u |^2 \zeta^2 \right)
 ( x, t_{0} ) \,
d x + \iint_{Q_{2 \rho}}
 | \frac{\partial u}{\partial t}|^2
 \zeta^2 \,
d x \, d t \\
& \leq \frac{1}{2} \int_{B_{2 \rho} ( x_{0} )}
 \Big(
  \sum_{i,j=1}^{n} a_{i j} \frac{\partial u}{\partial x_i}
  \frac{\partial u}{\partial x_{j}} \zeta^2
 \Big) ( x, t_{0} ) \,
d x + \iint_{Q_{2 \rho}}
 | \frac{\partial u}{\partial t}|^2
 \zeta^2 \,
d x \, d t  \\
& = \iint_{Q_{2 \rho}}
 \Big[
  \sum_{i,j=1}^{n} a_{i j} \frac{\partial u}{\partial x_i}
  \frac{\partial u}{\partial x_{j}}
  \zeta \frac{\partial \zeta}{\partial t}
  - 2 \sum_{i,j=1}^{n} a_{i j} \frac{\partial u}{\partial x_{j}}
  \frac{\partial u}{\partial t} \zeta
  \frac{\partial \zeta}{\partial x_i}
 \Big]
d x \, d t \\
& \quad + \iint_{Q_{2 \rho}}
 f \frac{\partial u}{\partial t} \zeta^2 \,
d x \, d t
+ \sum_{i=1}^{n} \Bigl[
 \int_{B_{2 \rho} ( x_{0} )}
  \left( f_i \frac{\partial u}{\partial x_i} \zeta^2 \right)
  ( x, t_{0} ) \,
 d x \\
 & \quad  + \iint_{Q_{2 \rho}}
  \Big(
   - \frac{\partial f_i}{\partial t}
   \frac{\partial u}{\partial x_i} \zeta^2
   - 2 f_i \frac{\partial u}{\partial x_i}
   \zeta \frac{\partial \zeta}{\partial t}
   + 2 f_i \frac{\partial u}{\partial t}
   \zeta \frac{\partial \zeta}{\partial x_i}
  \Big)
 \Bigr]
d x \, d t  \\
& \leq C \iint_{Q_{2 \rho}}
 | \nabla u |^2 \zeta
 |\frac{\partial \zeta}{\partial t}|
d x \, d t
+ \varepsilon_2 \iint_{Q_{2 \rho}}
 | \frac{\partial u}{\partial t}|^2
 \zeta^2 \,
d x \, d t \\
& \quad
+ \frac{C}{\varepsilon_2} \iint_{Q_{2 \rho}}
 | \nabla u |^2
 | \nabla \zeta |^2 \,
d x \, d t
+ \varepsilon_2 \iint_{Q_{2 \rho}}
 | \frac{\partial u}{\partial t}|^2
 \zeta^2 \,
d x \, d t \\
& \quad
+ \frac{C}{\varepsilon_2} \iint_{Q_{2 \rho}}
 | f |^2 \zeta^2 \,
d x \, d t
+ \varepsilon_2 \int_{B_{2 \rho} ( x_{0} )}
 \big( | \nabla u |^2 \zeta^2 \big) ( x, t_{0} ) \,
d x \\
& \hspace*{3ex} \text{}
+ \frac{C}{\varepsilon_2} \int_{B_{2 \rho} ( x_{0} )}
 \Big( \sum_{i=1}^{n} | f_i |^2 \zeta^2 \Big)
 ( x, t_{0} ) \,
d x \\
& \quad
+ C \iint_{Q_{2 \rho}}
 | \nabla u |^2 \zeta^2 \,
d x \, d t
+ C \iint_{Q_{2 \rho}}
 \sum_{i=1}^{n}
 \big| \frac{\partial f_i}{\partial t} \big|^2
 \zeta^2 \,
d x \, d t \\
& \quad
+ C \iint_{Q_{2 \rho}}
 | \nabla u |^2 \zeta
 |\frac{\partial \zeta}{\partial t}|
d x \, d t
+ C \iint_{Q_{2 \rho}}
 \sum_{i=1}^{n}
 | f_i |^2
 \zeta |\frac{\partial \zeta}{\partial t}|
d x \, d t \\
& \quad
+ \varepsilon_2 \iint_{Q_{2 \rho}}
 | \frac{\partial u}{\partial t}|^2
 \zeta^2 \,
d x \, d t
+ \frac{C}{\varepsilon_2} \iint_{Q_{2 \rho}}
 \sum_{i=1}^{n} | f_i |^2
 | \nabla \zeta |^2 \,
d x \, d t .
\end{align*}
We remark that
\begin{align*}
\int_{B_{2 \rho} ( x_{0} )}
 ( f_i \zeta )^2 ( x, t_{0} ) \,
d x
& = \int_{B_{2 \rho} ( x_{0} )}
 \int_{t_{0} - ( 2 \rho )^2}^{t_{0}}
  \frac{\partial}{\partial t}
  \left( ( f_i \zeta )^2 \right) ( x, t ) \,
 d t \,
d x \\
& \leq C \iint_{Q_{2 \rho}}
 \Big[
  | f_i |^2 \Big(
   \zeta^2 + \zeta | \frac{\partial \zeta}{\partial t}|
  \Big)
  + | \frac{\partial f_i}{\partial t} |^2
  \zeta^2
 \Big]
d x \, d t.
\end{align*}
Therefore, by taking $\varepsilon_2 > 0$ small enough, we have
\begin{align*}
\iint_{Q_{\rho}}
 | \frac{\partial u}{\partial t}|^2
d x \, d t
&\leq \int_{B_{2 \rho} ( x_{0} )}
 \left( | \nabla u |^2 \zeta^2 \right)
 ( x, t_{0} ) \,
d x + \iint_{Q_{2 \rho}}
 | \frac{\partial u}{\partial t}|^2
 \zeta^2 \,
d x \, d t \\
& \leq C \iint_{Q_{2 \rho}}
 | \nabla u |^2 \Big(
  \zeta^2
  + \zeta
  |\frac{\partial \zeta}{\partial t}|
  + | \nabla \zeta |^2
 \Big)
d x \, d t
+ C \iint_{Q_{2 \rho}}
 | f |^2 \zeta^2 \,
d x \, d t \\
& \hspace*{3ex} \text{}
+ C \iint_{Q_{2 \rho}}
 \sum_{i=1}^{n} \Big[
  | f_i |^2 \Big(
   \zeta^2
   + \zeta
   |\frac{\partial \zeta}{\partial t}|
   + | \nabla \zeta |^2
  \Big)
  + | \frac{\partial f_i}{\partial t}|^2
  \zeta^2
 \Big]
d x  \, d t  \\
& \leq C \rho^{-2} \| \nabla u \|_{L^2 ( Q_{2 \rho} )}^2
+ C \| f \|_{L^2 ( Q_{2 \rho} )}^2
+ C \rho^{-2} \sum_{i=1}^{n}
\| f_i \|_{L^2 ( Q_{2 \rho} )}^2 \\
& \quad
+ C \sum_{i=1}^{n} \|
 \frac{\partial f_i}{\partial t}
\|_{L^2 ( Q_{2 \rho} )}^2 .
\end{align*}
\end{proof}

We obtain the estimate \eqref{eq:dudtL2uL2}
from Lemmas~\ref{lemma:Linfty} (given in Appendix),
\ref{lemma:DuL2} and
\ref{lemma:dudtL2}.

\section{A gradient estimate of the fundamental solution}
\label{section:gradestoffundsol}

In this section, we consider a gradient estimate
of the fundametal solution of parabolic operators.
We first state some facts.
It is known that if coefficient $( a_{i j} )$ is
a symmetric and positive definite matrix-value
$L^{\infty}(\mathbb{R}^n)$ function
satisfying \eqref{eq:aijxiixij}, then
there exists a fundamental solution
$\Gamma (x,t;y,s)$ of the parabolic operator
\begin{equation}\label{eq:op}
\frac{\partial}{\partial t}
- \sum_{i,j=1}^{n} \frac{\partial}{\partial x_i} \Big(
 a_{i j} \frac{\partial}{\partial x_{j}}
\Big)
\end{equation}
with the estimate
\begin{equation}\label{eq:estoffundsol}
| \Gamma (x,t;y,s) |
\leq \frac{C_{\ast}}{(t-s)^{n/2}} \exp \Big(
 - \frac{c_{\ast} | x - y |^2}{t-s}
\Big) \chi_{[ s, \infty )} (t)
\end{equation}
for all
$t,s \in \mathbb{R}$,
and a.e.\ $x,y \in \mathbb{R}^{n}$,
where $C_{\ast} , c_{\ast} > 0$ depend only on
$n, \lambda , \Lambda$
(see \cite{Aronson} or \cite{FS}, for example).
In particular, the constants $C_{\ast}$ and $c_{\ast}$
are independent of the distance between inclusions.
If the coefficients $( a_{i j} )$ is not piecewise smooth
but H{\"o}lder continuous in the whole space $\mathbb{R}^{n}$,
then the pointwise gradient estimate
\[
| \nabla_{x} \Gamma (x,t;y,s) |
\leq \frac{C_{\ast}}{(t-s)^{(n+1)/2}}
\exp \Big(
 - \frac{c_{\ast} | x - y |^2}{t-s}
\Big) \chi_{[ s, \infty )} (t)
\]
holds for
$t,s \in \mathbb{R}$, a.e.\ $x,y \in \mathbb{R}^{n}$
(see \cite[Chapter IV \S11--13]{LSU}, for example).

Now, the aim of this section is to show
the gradient estimate \eqref{eq:gradientestimate}
in Theorem~\ref{theorem:gradientestimate}
even if the coefficients are piecewise $C^{\mu}$ in $D$.
We assume that $( a_{i j} )$ defined in $D$
satisfies the conditions
\eqref{eq:aijxiixij} and \eqref{eq:piecewiseCmu},
and extend it to the whole $\mathbb{R}^{n}$
by defining $( a_{i j} ) \equiv \Lambda I$
in $\mathbb{R}^{n} \setminus D$,
where $I$ is the identity matrix.
We remark that this extension does not destroy the conditions
 \eqref{eq:aijxiixij} and \eqref{eq:piecewiseCmu}.
Then there exists a fundamental solution
$\Gamma (x,t;y,s)$ of the parabolic operator \eqref{eq:op}
with the estimate \eqref{eq:estoffundsol}
as we stated above.


To prove our gradient estimate of the fundamental solution,
we apply the following corollary from
Theorem~\ref{theorem:main}.

\begin{corollary}\label{corollary:1}
Let $0 < \rho \leq 1$. Then a solution $u$ to the parabolic equation
\begin{equation}\label{eq:parabolicrho}
\frac{\partial u}{\partial t}
- \sum_{i,j=1}^{n} \frac{\partial}{\partial x_i} \Big(
 a_{i j} \frac{\partial u}{\partial x_{j}}
\Big) = 0 \quad\text{in }
B_{\rho} ( x_{0} ) \times ( t_{0} - \rho^2 , \, t_{0} ]
\end{equation}
has the estimate
\begin{equation}\label{eq:estimaterho}
\| \nabla u \|_{L^{\infty} (
  B_{\rho / 2} ( x_{0} ) \times ( t_{0} - ( \rho / 2 )^2 , \, t_{0}] )
 )
} \leq \frac{C_{\sharp}'}{\rho^{n/2+2}}
\| u \|_{L^2 (
  B_{\rho} ( x_{0} ) \times ( t_{0} - \rho^2 , \, t_{0} ]
 )
} ,
\end{equation}
where $C_{\sharp}' > 0$ depends only on
$n, L, \mu, \alpha , \lambda , \Lambda$,
and
\begin{math}
\| a_{i j} \|_{C^{\alpha'} ( \overline{D_{m}})}
\end{math}
and the $C^{1, \alpha'}$ norms of $D_{m}$
for some $\alpha'$ with {\rm \eqref{eq:alphaprime}}.
\end{corollary}

\begin{proof}
It is sufficient to apply the scaling argument. Let
$\rho y = x - x_{0}$,
$\rho^2 ( s - 1 ) = t - t_{0}$
and
\begin{equation} \label{eq:tilde}
\begin{gathered}
\widetilde{u} (y,s) := u(x,t)
= u \bigl( \rho y + x_{0} , \, \rho^2 (s-1) + t_{0} \bigr) ,
  \\
 \widetilde{a}_{i j} (y)
:= a_{i j} (x) = a_{i j} ( \rho y + x_{0} ),
 \\
 \widetilde{D}_{m}
:= \big\{
 \frac{1}{\rho} ( x - x_{0} ) : x \in D_{m}
\big\} .
\end{gathered}
\end{equation}
Then we have
\begin{equation}\label{eq:tildeu}
\frac{\partial \widetilde{u}}{\partial s}
- \sum_{i,j=1}^{n} \frac{\partial}{\partial y_i} \Big(
 \widetilde{a}_{i j} \frac{\partial \widetilde{u}}{\partial y_{j}}
\Big) = 0 \quad\text{in } B_1 (0) \times (0,1].
\end{equation}
Therefore, by noting Remark~\ref{remark:indepofrho}, we have
\[
\| \nabla \widetilde{u} \|_{
 L^{\infty} ( B_{1/2} (0) \times ( 3/4, 1 ] )
} \leq C_{\sharp}'
\| \widetilde{u} \|_{L^2 ( B_1 (0) \times (0,1) )}
\]
by Theorem~\ref{theorem:main},
where $C_{\sharp}'$ depends only on
\begin{math}
n, L, \mu, \alpha , \lambda , \Lambda ,
\| a_{i j} \|_{C^{\alpha'} ( \overline{D_{m}})}
\end{math},
and the $C^{1, \alpha'}$ seminorms of $D_{m}$.
By this estimate and the definition \eqref{eq:tilde},
we obtain the estimate \eqref{eq:estimaterho}.
\end{proof}

\begin{remark}\label{remark:indepofrho} \rm
One may think that  $C_{\sharp}'$
depends also on $\rho$
since
\begin{math}
\| \widetilde{a}_{i j} \|_{
 C^{\alpha'} ( \overline{\widetilde{D}_{m}} )
}
\end{math}
and the $C^{1, \alpha'}$ norms of
$\widetilde{D}_{m}$
depend on $\rho$.
However, we can take $C_{\sharp}'$
independent of $\rho$
by taking the following into consideration.

First we consider
\begin{align*}
\| \widetilde{a}_{i j} \|_{
 C^{\alpha'} ( \overline{\widetilde{D}_{m}} )
}
& = \| \widetilde{a}_{i j} \|_{
 C^{0} ( \overline{\widetilde{D}_{m}} )
} + [ \widetilde{a}_{i j} ]_{
 C^{\alpha'} ( \overline{\widetilde{D}_{m}} )
} \\
& := \sup_{y \in \overline{\widetilde{D}_{m}}}
| \widetilde{a}_{i j} (y) |
+ \sup_{y, \eta \in \overline{\widetilde{D}_{m}}}
\frac{
 | \widetilde{a}_{i j} (y) - \widetilde{a}_{i j} ( \eta ) |
}{| y - \eta |^{\alpha'}} .
\end{align*}
It is easy to show
\[
\| \widetilde{a}_{i j} \|_{
 C^{0} ( \overline{\widetilde{D}_{m}} )
} = \| a_{i j} \|_{C^{0} ( \overline{D_{m}} )}
\]
and
\[
[ \widetilde{a}_{i j} ]_{
 C^{\alpha'} ( \overline{\widetilde{D}_{m}} )
} = \rho^{\alpha'} [ a_{i j} ]_{
 C^{\alpha'} ( \overline{D_{m}} )
} \leq [ a_{i j} ]_{C^{\alpha'} ( \overline{D_{m}} )} .
\]
Then we have
\[
\| \widetilde{a}_{i j} \|_{
 C^{\alpha'} ( \overline{\widetilde{D}_{m}} )
}
\leq \| a_{i j} \|_{C^{\alpha'} ( \overline{D_{m}} )} .
\]

Next we consider the $C^{1, \alpha'}$ norms of
$\widetilde{D}_{m}$.
We need to recall the proofs of the results of
\cite{LiNirenberg} and \cite{LiVogelius}
more carefully.
In the case when we consider
the $L^{\infty}$-norm of $\nabla \widetilde{u}$
for a solution $\widetilde{u}$ to the equation \eqref{eq:tildeu},
the influence of the $C^{1, \alpha'}$ norms of
subdomains $\widetilde{D}_{m}$ appears only
in the following constant $C$ in \eqref{eq:LiVogeliusp118C}:
We estimate
$O \bigl( | x' |^{1 + \alpha} \bigr)$
in the equation (49) in \cite[p.\ 118]{LiVogelius}, i.e.\
\begin{equation}\tag{49}%\tag*{(49) in \cite{LiVogelius}}
f_{m} ( x' )
= f_{m} ( 0' )
+ \nabla f_{m} ( 0' ) x'
+ O \bigl( | x' |^{1 + \alpha} \bigr)
\end{equation}
as
\begin{equation}\label{eq:LiVogeliusp118C}
\bigl|
 O \bigl( | x' |^{1 + \alpha} \bigr)
\bigr|
\leq C | x |^{1 + \alpha}
\end{equation}
(See also \cite[Lemma~4.3]{LiNirenberg}).
Here $C^{1, \alpha}$ functions $f_{m}$
are defined in the cube $(-1,1)^{n}$,
and the graphs of $f_{m}$ describe $\partial D_{m}$.
Now we remark that the constant $C$ in \eqref{eq:LiVogeliusp118C}
depends only on the $C^{1, \alpha}$ seminorms of $f_{m}$.
We consider the variable change $\rho y = x$.
Then the graph $x_n = f_{m} ( x' )$
is changed to $y_n = \widetilde{f}_{m} ( y' )$,
where
\begin{math}
\widetilde{f}_{m} ( y' )
:= \rho^{-1} f_{m} ( \rho y' )
\end{math}, and we have
\begin{align*}
[ \widetilde{f}_{m} ]_{C^{1, \alpha} ( (-1,1)^{n} )}
& \leq [ \widetilde{f}_{m} ]_{
 C^{1, \alpha} ( ( - 1 / \rho , 1 / \rho )^{n} )
} \\
& = \rho^{\alpha} [ f_{m} ]_{C^{1, \alpha} ( (-1,1)^{n} )}
\leq [ f_{m} ]_{C^{1, \alpha} ( (-1,1)^{n} )} .
\end{align*}
Hence, even when we consider the variable change $\rho y = x$,
we can take the constant $C$ in \eqref{eq:LiVogeliusp118C}
independent of $\rho$.

Considering the circumstances mentioned above,
we can take $C_{\sharp}' > 0$
independent of $\rho$.
\end{remark}

Now we state the estimate of
$\nabla_{x} \Gamma (x,t;y,s)$.

\begin{theorem}\label{theorem:gradientestimate}
We have
\begin{equation}\label{eq:gradientestimate}
| \nabla_{x} \Gamma (x,t;y,s) |
\leq \frac{C}{( t - s )^{(n+1)/2}}
\exp \Big(
 - \frac{c | x - y |^2}{t-s}
\Big)
\end{equation}
for a.e.\ $x, y \in \mathbb{R}^{n}$ and $t > s$
with $| x - y |^2 + t - s \leq 16$,
where
$C, c>0$ depend only on
\begin{math}
n, L, \mu, \alpha , \lambda , \Lambda ,
\| a_{i j} \|_{C^{\alpha'} ( \overline{D_{m}})}
\end{math}
and the $C^{1, \alpha'}$ seminorms of $D_{m}$
for some $\alpha'$ with {\rm \eqref{eq:alphaprime}}.
\end{theorem}

We prove Theorem~\ref{theorem:gradientestimate}
in the same way as the proof of
\cite[Proposition 3.6]{DiCristoVessella}.
We first show the following lemmas.

\begin{lemma}\label{lemma:GammaL2a}
Let
$\rho := ( | x_{0} - \xi |^2 + t_{0} - \tau )^{1/2} / 4$.
Then
\[
\int_{t_{0} - \rho^2}^{t_{0}} \int_{B_{\rho} ( x_{0} )}
 | \Gamma( x, t; \xi , \tau ) |^2 \,
d x \, d t
\leq \frac{( C_{\ast}' )^2 \rho^{n}}{( t_{0} - \tau )^{n-1}}
\exp \Big(
 - \frac{2 c_{\ast}' | x_{0} - \xi |^2}
        {t_{0} - \tau}
\Big)
\]
for $t_{0} > \tau$, where
$C_{\ast}' , c_{\ast}' > 0$
depend only on $n, \lambda , \Lambda$.
\end{lemma}

\begin{proof}
By \eqref{eq:estoffundsol}, it is sufficient to obtain the estimate
\begin{equation}
\begin{aligned}
I_{0}
& :=  \int_{t_{0} - \rho^2}^{t_{0}} \int_{B_{\rho} ( x_{0} )}
 \frac{1}{( t - \tau )^{n}}
 \exp \Big(
  - \frac{2 c_{\ast} | x - \xi |^2}{t - \tau}
 \Big) \chi_{[ \tau , \infty )} (t) \,
d x \, d t  \\
& \leq
\frac{( C_{\ast}' )^2 \rho^{n}}{( t_{0} - \tau )^{n-1}}
\exp \Big(
 - \frac{2 c_{\ast}' | x_{0} - \xi |^2}
        {t_{0} - \tau}
\Big).
\end{aligned} \label{eq:I0}
\end{equation}
We consider the following three cases:
\begin{itemize}
\item[(i)]  $t_{0} - \rho^2 \leq \tau < t_{0}$,
\item[(i)]  $t_{0} - 2 \rho^2 \leq \tau \leq t_{0} - \rho^2$,
\item[(i)] $\tau \leq t_{0} - 2 \rho^2$.
\end{itemize}
Now we consider the case (i).
Then we have
$( \sqrt{15} - 1 ) \rho \leq | x - \xi |$
for any $x \in B_{\rho} ( x_{0} )$, because
$| x_{0} - \xi | \geq \sqrt{15} \, \rho$.
Hence we have
\[
I_{0} \leq
\int_{\tau}^{t_{0}} \int_{B_{\rho} ( x_{0} )}
 \frac{1}{( t - \tau )^{n}}
 \exp \Big(
  - \frac{c_1 \rho^2}{t - \tau}
 \Big)
d x \, d t
= | B_1 (0) | \rho^{n}
\int_{0}^{t_{0} - \tau} \varphi_1 (s) \, d s,
\]
where $\varphi_1 (s) := s^{-n} \exp ( - c_1 \rho^2 / s )$ and
$c_1 := 2 ( \sqrt{15} - 1 )^2 c_{\ast}$.
If $0 < t_{0} - \tau \leq c_1 \rho^2 / n$,
then we have
\[
\int_{0}^{t_{0} - \tau} \varphi_1 (s) \, d s
\leq \int_{0}^{t_{0} - \tau} \varphi_1 ( t_{0} - \tau ) \, d s
= ( t_{0} - \tau )^{-n+1} \exp \Big(
 - \frac{c_1 \rho^2}{t_{0} - \tau}
\Big) ,
\]
because $\varphi_1 (s) \leq \varphi_1 ( t_{0} - \tau )$ holds for any
$s \in [ 0, \, t_{0} - \tau ]$.
On the other hand, if
$c_1 \rho^2 / n \leq t_{0} - \tau \leq \rho^2$,
then we have
\begin{align*}
\int_{0}^{t_{0} - \tau} \varphi_1 (s) \, d s
& \leq \int_{0}^{t_{0} - \tau}
 \varphi_1 \big( \frac{c_1 \rho^2}{n} \big)
d s
= \big( \frac{n}{c_1} \big)^{n} ( t_{0} - \tau )
\rho^{- 2 n} \exp ( - n ) \\
& \leq \big( \frac{n}{c_1} \big)^{n} ( t_{0} - \tau )^{1-n}
\exp \Big(
 - \frac{c_1 \rho^2}{t_{0} - \tau}
\Big),
\end{align*}
where we used the properties that
\begin{gather*}
 \varphi_1 (s)
\leq \varphi_1 \Big( \frac{c_1 \rho^2}{n} \Big)
\quad\text{for any } 0 < s \leq t_{0} - \tau ;
 \\
 n \geq \frac{c_1 \rho^2}{t_{0} - \tau} \quad \text{and}\quad
\rho^2 \geq t_{0} - \tau .
\end{gather*}
Summing up, we have
\[
I_{0} \leq \max \big\{
 1, \, \big( \frac{n}{c_1} \big)^{n}
\big\} | B_1 (0) |
\rho^{n} ( t_{0} - \tau )^{1-n} \exp \Big(
 - \frac{c_1 \rho^2}{t_{0} - \tau}
\Big) .
\]

Let us consider the case (ii).
Then we have $( \sqrt{14} - 1 ) \rho \leq | x - \xi |$
for all $x \in B_{\rho} ( x_{0} )$,
because $| x_{0} - \xi | \geq \sqrt{14} \, \rho$.
Hence we have
\[
I_{0}
 \leq \int_{t_{0} - \rho^2}^{t_{0}} \int_{B_{\rho} ( x_{0} )}
 \frac{1}{( t - \tau )^{n}}
 \exp \Big(
  - \frac{c_2 \rho^2}{t - \tau}
 \Big)
d x \, d t \\
 = | B_1 (0) | \rho^{n}
\int_{t_{0} - \rho^2 - \tau}^{t_{0} - \tau}
 \varphi_2 (s) \,
d s,
\]
where $\varphi_2 (s) := s^{-n} \exp ( - c_2 \rho^2 / s )$ and
$c_2 := 2 ( \sqrt{14} - 1 )^2 c_{\ast}$.
In a similarly way as the case (i),
if $\rho^2 \leq t_{0} - \tau \leq c_2 \rho^2 / n$,
then we have
\begin{align*}
\int_{t_{0} - \rho^2 - \tau}^{t_{0} - \tau}
 \varphi_2 (s) \,
d s
& \leq \int_{t_{0} - \rho^2 - \tau}^{t_{0} - \tau}
 \varphi_2 ( t_{0} - \tau ) \,
d s
= \rho^2 ( t_{0} - \tau )^{-n}
\exp \Big( - \frac{c_2 \rho^2}{t_{0} - \tau} \Big) \\
& \leq ( t_{0} - \tau )^{-n+1}
\exp \Big( - \frac{c_2 \rho^2}{t_{0} - \tau} \Big),
\end{align*}
because $\varphi_2 (s) \leq \varphi ( t_{0} - \tau )$
for any $s \in [ t_{0} - \rho^2 - \tau , \, t_{0} - \tau ]$,
and we have $\rho^2 \leq t_{0} - \tau$.
On the other hand,
if $c_2 \rho^2 / n \leq t_{0} - \tau \leq 2 \rho^2$,
then we have
\begin{align*}
\int_{t_{0} - \rho^2 - \tau}^{t_{0} - \tau}
 \varphi_2 (s) \,
d s
& \leq \int_{t_{0} - \rho^2 - \tau}^{t_{0} - \tau}
 \varphi_2 \left( \frac{c_2 \rho^2}{n} \right)
d s
= \left( \frac{n}{c_2} \right)^{n}
\rho^{-2n+2} \exp ( -n ) \\
& \leq 2^{n-1} \left( \frac{n}{c_2} \right)^{n}
( t_{0} - \tau )^{1-n}
\exp \Big( - \frac{c_2 \rho^2}{t_{0} - \tau} \Big),
\end{align*}
where we used the properties that
\begin{gather*}
 \varphi_2 (s)
\leq \varphi_2 \Big( \frac{c_2 \rho^2}{n} \Big)
\quad\text{for any } t_{0} - \rho^2 - \tau \leq s \leq t_{0} - \tau ;
 \\
 n \geq \frac{c_2 \rho^2}{t_{0} - \tau} \text{, and }
\rho^2 \geq \frac{t_{0} - \tau}{2} .
\end{gather*}
Summing up, we have
\[
I_{0} \leq | B_1 (0) |
\max \big\{ 1, \, 2^{n-1} \big( \frac{n}{c_2} \big)^{n} \big\}
\rho^{n} ( t_{0} - \tau )^{1-n}
\exp \Big( - \frac{c_2 \rho^2}{t_{0} - \tau} \Big) .
\]
Finally we consider the case (iii).
We first remark that
\[
\int_{t_{0} - \rho^2}^{t_{0}} ( t - \tau )^{-n} \, d t
\leq
 \begin{cases}
  \frac{1}{n-1} ( t_{0} - \rho^2 - \tau )^{-n+1}
 &  \text{if } n \geq 2, \\
  \log 2 & \text{if } n = 1,
 \end{cases}
\]
because $t_{0} - \tau \leq 2 ( t_{0} - \rho^2 - \tau )$.
In particular, we have
\[
\int_{t_{0} - \rho^2}^{t_{0}} ( t - \tau )^{-n} \, d t
\leq ( t_{0} - \rho^2 - \tau )^{-n+1}
\leq 2^{n-1} ( t_{0} - \tau )^{-n+1} .
\]
Hence we have
\begin{align*}
I_{0}
& \leq | B_1 (0) | \rho^{n}
\int_{t_{0} - \rho^2}^{t_{0}} ( t - \tau )^{-n} \, d t
\leq 2^{n-1} | B_1 (0) | \rho^{n}
( t_{0} -  \tau )^{-n+1} \\
& \leq 2^{n-1} \exp (8) | B_1 (0) | \rho^{n}
( t_{0} - \tau )^{-n+1}
\exp \Big(
 - \frac{| x_{0} - \xi |^2}{t_{0} - \tau}
\Big),
\end{align*}
because
\begin{math}
| x_{0} - \xi |^2 / ( t_{0} - \tau)
\leq ( 4 \rho )^2 / 2 \rho^2 = 8
\end{math}.
Therefore we have the estimate \eqref{eq:I0} in every case.
\end{proof}

\begin{proof}[Proof of Theorem~\ref{theorem:gradientestimate}]
Let $x_{0} , \xi \in \mathbb{R}^{n}$ and $t_{0} > \tau$.
Let
\begin{math}
\rho := ( | x_{0} - \xi |^2 + t_{0} - \tau )^{1/2} / 4
\leq 1
\end{math}.
Then, by Corollary~\ref{corollary:1},
we have
\begin{align*}
& \| \nabla_{x} \Gamma ( \cdot , \cdot ; \xi , \tau ) \|_{
 L^{\infty} (
  B_{\rho / 2} ( x_{0} ) \times ( t_{0} - ( \rho / 2 )^2 , \, t_{0} )
 )
} \\
& \leq \frac{C_{\sharp}'}{\rho^{n/2 + 2}}
\|
 \Gamma ( \cdot , \cdot ; \xi , \tau )
\|_{L^2
 ( B_{\rho} ( x_{0} ) \times ( t_{0} - \rho^2 , \, t_{0} ) )
},
\end{align*}
because we have
\[
\frac{\partial \Gamma}{\partial t} ( x, t; \xi , \tau )
- \sum_{i,j=1}^{n} \frac{\partial}{\partial x_i} \Big(
 a_{i j} (x) \frac{\partial \Gamma}{\partial x_{j}}
 ( x, t; \xi , \tau )
\Big) = 0 \quad\text{in }
B_{\rho} ( x_{0} ) \times ( t_{0} - \rho^2 , t_{0} ] .
\]
By this estimate and Lemma~\ref{lemma:GammaL2a}, we have
\begin{align*}
& \| \nabla_{x} \Gamma ( \cdot , \cdot ; \xi , \tau ) \|_{
 L^{\infty} (
  B_{\rho / 2} ( x_{0} ) \times ( t_{0} - ( \rho / 2 )^2 , \, t_{0} ]
 )
} \\
& \leq \frac{C_{\sharp}'}{\rho^{n/2 + 2}}
\|
 \Gamma ( \cdot , \cdot ; \xi , \tau )
\|_{L^2
 ( B_{\rho} ( x_{0} ) \times ( t_{0} - \rho^2 , \, t_{0} ] )
}  \\
& \leq \frac{C_{\sharp}' C_{\ast}'}{\rho^2}
\frac{1}{( t_{0} - \tau )^{(n-1)/2}}
\exp \Big(
 - \frac{c_{\ast}' | x_{0} - \xi |^2}
        {t_{0} - \tau}
\Big)  \\
& \leq
\frac{16 C_{\sharp}' C_{\ast}'}
     {( t_{0} - \tau )^{(n+1)/2}}
\exp \Big(
 - \frac{c_{\ast}' | x_{0} - \xi |^2}
        {t_{0} - \tau}
\Big),
\end{align*}
because we have
$\rho^{-2} \leq 16 ( t_{0} - \tau )^{-1}$.
Hence the proof is complete.
\end{proof}

\section{Appendix}

In Appendix, we show the estimates
\eqref{eq:LinftyL2} and \eqref{eq:Linfty}
in Lemma~\ref{lemma:keyestimate}
for the sake of completeness.
To begin with,
we give some embedding lemma
which is necessary to show the estimates
\eqref{eq:LinftyL2} and \eqref{eq:Linfty}.
First, the following Gagliardo-Nirenberg's inequality is well known
(see \cite[p.\ 24, Theorem~9.3]{Friedman}, for example).
\begin{lemma}[Gagliardo-Nirenberg's inequality]
\label{lemma:GagliardoNirenberg}
Let
$r, s$
be any numbers satisfying
$1 \leq r, s \leq \infty$,
and let
$j, k$
be any integers satisfying
$0 \leq j < k$.
If $u$ is any function in
$W_{s}^{k} ( \mathbb{R}^{n} ) \cap L^{r} ( \mathbb{R}^{n} )$,
then
\begin{equation}\label{eq:GagliardoNirenberg}
\| D^{j} u \|_{L^{q}( \mathbb{R}^{n} )}
\leq C_1 \| D^{k} u \|_{L^{s}( \mathbb{R}^{n} )}^{\gamma}
\| u \|_{L^{r}( \mathbb{R}^{n} )}^{1 - \gamma} ,
\end{equation}
where
\begin{equation}\label{eq:gamma}
\frac{1}{q} = \frac{j}{n} +
\gamma \Big( \frac{1}{s} - \frac{k}{n} \Big)
+ \frac{1 - \gamma}{r}
\end{equation}
for all $\gamma$ in the interval
\[
\frac{j}{k} \leq \gamma \leq 1,
\]
where a positive constant $C_1$
depends only on $n, k, j, r, s, \gamma$,
with the following exception:
If $k - j - n/s$ is a nonnegative integer,
then {\rm \eqref{eq:GagliardoNirenberg}} holds
only for $j/k \leq \gamma < 1$.
\end{lemma}


Then, as an application of Lemma~\ref{lemma:GagliardoNirenberg},
we have the following embedding lemma.

\begin{lemma}[Embedding lemma]\label{lemma:SobolevV2}
Let
\begin{math}
v \in L^{\infty} \bigl( 0, T; L^2 (D) \bigr)
\cap L^2 \bigl( 0, T; H_{0}^{1} (D) \bigr)
\end{math}.
Then $v \in L^{2(n+2)/n} (Q)$ holds.
Moreover, we have the estimate
\begin{equation}
\begin{aligned}
\| v \|_{L^{2(n+2)/n} (Q)}
& \leq C_1
\| v \|_{L^{\infty} ( 0, T; L^2 (D) )}^{2/(n+2)}
\| \nabla v \|_{L^2 (Q)}^{n/(n+2)}  \\
& \leq C_1 \left(
 \| v \|_{L^{\infty} ( 0, T; L^2 (D) )}
 + \| \nabla v \|_{L^2 (Q)}
\right) , 
\end{aligned} \label{eq:V2}
\end{equation}
where a positive constant $C_1$ depends only on $n$,
and we denote $Q := D \times (0,T]$.
\end{lemma}

\begin{proof}
We apply Lemma~\ref{lemma:GagliardoNirenberg}
with $q = 2(n+2)/n$, $r=2$, $s=2$, $k=1$ and $j=0$.
Then the equation \eqref{eq:gamma} yields $\gamma = n/(n+2)$.
Hence we have
\[
\| v ( \cdot , t ) \|_{L^{2(n+2)/n} (D)}
\leq C_1 \| \nabla v ( \cdot , t ) \|_{L^2 (D)}^{n/(n+2)}
\| v ( \cdot , t ) \|_{L^2 (D)}^{2/(n+2)} .
\]
Therefore, we have
\begin{align*}
\| v \|_{L^{2(n+2)/n} (Q)}^{2(n+2)/n}
& = \int_{0}^{T}
 \| v ( \cdot , t ) \|_{L^{2(n+2)/n} (D)}^{2(n+2)/n} \,
d t \\
& \leq \int_{0}^{T}
 \left(
  C_1 \| \nabla v ( \cdot , t ) \|_{L^2 (D)}^{n/(n+2)}
  \| v ( \cdot , t ) \|_{L^2 (D)}^{2/(n+2)}
 \right)^{2(n+2)/n} \,
d t  \\
& \leq C_1^{2(n+2)/n}
\| v \|_{L^{\infty} ( 0, T; L^2 (D) )}^{4/n}
\| \nabla v \|_{L^2 (Q)}^2 .
\end{align*}
By this inequality and Young's inequality,
we have the estimate \eqref{eq:V2}.
\end{proof}
Based on Di~Giorgi's famous argument, we start to estimate
solutions to the parabolic equation
\eqref{eq:parabolic}.
By testing $( u - k )_{+} \zeta^2$ to \eqref{eq:parabolic}
we have the following lemma,
where $v_{+} (x) := \max \{ v(x), 0 \}$.

\begin{lemma}\label{lemma:DG}
Let $p > 2$.
Let
\begin{math}
Q_{\rho} := B_{\rho} ( x_{0} ) \times ( t_{0} - \rho^2 , t_{0} ]
\subset Q
\end{math}
and
\begin{math}
\zeta \in C^{\infty} \bigl(
 [ t_{0} - \rho^2 , t_{0} ] ; C_{0}^{\infty} ( B_{\rho} ( x_{0} ) )
\bigr)
\end{math}
satisfy
$0 \leq \zeta \leq 1$
and
$\zeta ( \cdot , t_{0} - \rho^2 ) = 0$.
Then a solution $u$ to the parabolic equation
 \eqref{eq:parabolic} satisfies
\begin{equation}
\begin{aligned}
& \| ( u - k )_{+} \zeta \|_{
 L^{\infty} ( t_{0} - \rho^2 , t_{0} ; L^2 ( B_{\rho} ( x_{0} ) ) )
}^2
+ \bigl\|
 \nabla \bigl( ( u - k )_{+} \zeta \bigr)
\bigr\|_{L^2 ( Q_{\rho} )}^2  \\
& \leq C_2 \Bigl[
 \Big(
  \|
   \frac{\partial \zeta}{\partial t}
  \|_{L^{\infty} ( Q_{\rho} )}
  + \| \nabla \zeta \|_{L^{\infty} ( Q_{\rho} )}^2
 \Big) \| ( u - k )_{+} \|_{L^2 ( Q_{\rho} )}^2
  \\
 & \quad  + F_{0, \rho}^2
 \bigl|
  Q_{\rho}  \cap \{ u(x,t) > k \}
 \bigr|^{1 - 2 / p}
\Bigr]
\end{aligned} \label{eq:DG}
\end{equation}
for any $k \in \mathbb{R}$,
where
\begin{equation}\label{eq:F0rho}
F_{0, r} := \| f \|_{L^{\frac{p(n+2)}{n+2+p}} ( Q_{r} )}
+ \sum_{i=1}^{n} \| f_i \|_{L^{p} ( Q_{r} )}
\quad\text{for } r>0
\end{equation}
and $C_2 > 0$ depends only on $n, \Lambda$ and $\lambda$.
\end{lemma}

\begin{proof}
Multiplying \eqref{eq:parabolic} by
$( u - k )_{+} \zeta^2$
and integrating it over
\begin{math}
Q_{\rho}'
:= B_{\rho} ( x_{0} ) \times ( t_{0} - \rho^2 , t' )
\end{math}
(also see Remark~\ref{remark:steklov}),
we have
\begin{align*}
\text{(LHS)}
& = \iint_{Q_{\rho}'}
 \Big( \frac{\partial}{\partial t} ( u - k )_{+} \Big)
 ( u - k )_{+} \zeta^2 \,
d x \, d t \\
& \quad
- \sum_{i,j=1}^{n}
\iint_{Q_{\rho}'}
 \frac{\partial}{\partial x_i} \Big(
  a_{i j} \frac{\partial}{\partial x_{j}} ( u - k )_{+}
 \Big)
 ( u - k )_{+} \zeta^2 \,
d x \, d t  \\
& = \frac{1}{2}
\iint_{Q_{\rho}'}
 \Big( \frac{\partial}{\partial t} ( u - k )_{+}^2 \Big)
 \zeta^2 \,
d x \, d t \\
& \quad
+ \sum_{i,j=1}^{n}
\iint_{Q_{\rho}'}
 a_{i j} \frac{\partial}{\partial x_{j}} ( u - k )_{+}
 \frac{\partial}{\partial x_i} \bigl( ( u - k )_{+} \zeta^2 \bigr)
 \,
d x \, d t  \\
& = \frac{1}{2}
\iint_{Q_{\rho}'}
 \Big[
  \frac{\partial}{\partial t} \bigl(
   ( u - k )_{+}^2 \zeta^2
  \bigr)
  - 2 ( u - k )_{+}^2 \zeta \frac{\partial \zeta}{\partial t}
 \Big]
d x \, d t \\
& \quad + \sum_{i,j=1}^{n}
\iint_{Q_{\rho}'}
 a_{i j}
 \frac{\partial}{\partial x_{j}} \bigl( ( u - k )_{+} \zeta \bigr)
 \frac{\partial}{\partial x_i} \bigl( ( u - k )_{+} \zeta \bigr) \,
d x \, d t \\
& \quad - \sum_{i,j=1}^{n}
\iint_{Q_{\rho}'}
 a_{i j} ( u - k )_{+}^2 \frac{\partial \zeta}{\partial x_{j}}
 \frac{\partial \zeta}{\partial x_i} \,
d x \, d t  \\
& =
\frac{1}{2}
\int_{B_{\rho} ( x_{0} )}
 ( u - k )_{+}^2 \zeta^2 \,
d x \biggr|_{t = t'}
- \iint_{Q_{\rho}'}
 ( u - k )_{+}^2 \zeta \frac{\partial \zeta}{\partial t} \,
d x \, d t \\
& \quad + \sum_{i,j=1}^{n}
\iint_{Q_{\rho}'}
 a_{i j}
 \frac{\partial}{\partial x_{j}} \bigl( ( u - k )_{+} \zeta \bigr)
 \frac{\partial}{\partial x_i} \bigl( ( u - k )_{+} \zeta \bigr) \,
d x \, d t \\
& \quad - \sum_{i,j=1}^{n}
\iint_{Q_{\rho}'}
 a_{i j} ( u - k )_{+}^2 \frac{\partial \zeta}{\partial x_{j}}
 \frac{\partial \zeta}{\partial x_i} \,
d x \, d t .
\end{align*}
Hence we have
\begin{equation}
\begin{aligned}
& \frac{1}{2} \int_{B_{\rho} ( x_{0} )}
 ( u - k )_{+}^2 \zeta^2 \,
d x \biggr|_{t = t'}  \\
& + \sum_{i,j=1}^{n}
\iint_{Q_{\rho}'}
 a_{i j}
 \frac{\partial}{\partial x_{j}} \bigl( ( u - k )_{+} \zeta \bigr)
 \frac{\partial}{\partial x_i} \bigl( ( u - k )_{+} \zeta \bigr) \,
d x \, d t  \\
& = \iint_{Q_{\rho}'}
 ( u - k )_{+}^2 \zeta \frac{\partial \zeta}{\partial t} \,
d x \, d t
+ \sum_{i,j=1}^{n}
\iint_{Q_{\rho}'}
 a_{i j} ( u - k )_{+}^2 \frac{\partial \zeta}{\partial x_{j}}
 \frac{\partial \zeta}{\partial x_i} \,
d x \, d t  \\
& \quad + \iint_{Q_{\rho}'}
 f ( u - k )_{+} \zeta^2 \,
d x \, d t + \sum_{i=1}^{n}
\iint_{Q_{\rho}'}
 f_i
 \frac{\partial}{\partial x_i} \bigl( ( u - k )_{+} \zeta^2 \bigr)
 \,
d x \, d t .
\end{aligned} \label{eq:identity1}
\end{equation}
We remark that
\begin{align*}
& \Big|
 \iint_{Q_{\rho}'}
  f_i
  \frac{\partial}{\partial x_i} \bigl( ( u - k )_{+} \zeta^2 \bigr)
  \,
 d x \, d t
\Big| \\
& = \Bigl|
 \iint_{Q_{\rho}' \cap \{ u(x,t) > k \}}
  f_i \zeta
  \frac{\partial}{\partial x_i} \bigl( ( u - k )_{+} \zeta \bigr) \,
 d x \, d t \\
 & \quad  + \iint_{Q_{\rho}' \cap \{ u(x,t) > k \}}
  f_i ( u - k )_{+} \zeta
  \frac{\partial \zeta}{\partial x_i} \,
 d x \, d t
\Bigr|  \\
& \leq \varepsilon_1
\iint_{Q_{\rho}' \cap \{ u(x,t) > k \}}
 \big|
  \frac{\partial}{\partial x_i} \bigl( ( u - k )_{+} \zeta \bigr)
 \big|^2 \,
d x \, d t
+ \frac{1}{\varepsilon_1}
\iint_{Q_{\rho}' \cap \{ u(x,t) > k \}}
 | f_i \zeta |^2 \,
d x \, d t \\
& \quad + \iint_{Q_{\rho}' \cap \{ u(x,t) > k \}}
 | f_i \zeta |^2 \,
d x \, d t
+ \iint_{Q_{\rho}' \cap \{ u(x,t) > k \}}
 ( u - k )_{+}^2
 \big| \frac{\partial \zeta}{\partial x_i} \big|^2 \,
d x \, d t .
\end{align*}
Hence, by \eqref{eq:aijxiixij} and \eqref{eq:identity1}, we have
\begin{align*}
& \frac{1}{2} \int_{B_{\rho} ( x_{0} )}
 ( u - k )_{+}^2 \zeta^2 \,
d x \biggr|_{t = t'}
+ \lambda
\iint_{Q_{\rho}'}
 \bigl| \nabla \bigl( ( u - k )_{+} \zeta \bigr) \bigr|^2 \,
d x \, d t \\
& \leq \frac{1}{2} \int_{B_{\rho} ( x_{0} )}
 ( u - k )_{+}^2 \zeta^2 \,
d x \biggr|_{t = t'} \\
& \quad + \sum_{i,j=1}^{n}
\iint_{Q_{\rho}'}
 a_{i j}
 \frac{\partial}{\partial x_{j}} \bigl( ( u - k )_{+} \zeta \bigr)
 \frac{\partial}{\partial x_i} \bigl( ( u - k )_{+} \zeta \bigr) \,
d x \, d t  \\
& = \iint_{Q_{\rho}'}
 ( u - k )_{+}^2 \zeta \frac{\partial \zeta}{\partial t} \,
d x \, d t
+ \sum_{i,j=1}^{n}
\iint_{Q_{\rho}'}
 a_{i j} ( u - k )_{+}^2 \frac{\partial \zeta}{\partial x_{j}}
 \frac{\partial \zeta}{\partial x_i} \,
d x \, d t \\
& \quad + \iint_{Q_{\rho}'}
 f ( u - k )_{+} \zeta^2 \,
d x \, d t + \sum_{i=1}^{n}
\iint_{Q_{\rho}'}
 f_i
 \frac{\partial}{\partial x_i} \bigl( ( u - k )_{+} \zeta^2 \bigr)
 \,
d x \, d t  \\
& \leq \|
 \frac{\partial \zeta}{\partial t}
\|_{L^{\infty} ( Q_{\rho} )}
\iint_{Q_{\rho}'}
 ( u - k )_{+}^2 \,
d x \, d t
+ \Lambda \| \nabla \zeta \|_{L^{\infty} ( Q_{\rho} )}^2
\iint_{Q_{\rho}'}
 ( u - k )_{+}^2 \,
d x \, d t \\
& \quad + \iint_{Q_{\rho}'}
 f ( u - k )_{+} \zeta^2 \,
d x \, d t
+ \varepsilon_1
\iint_{Q_{\rho}'}
 \left|
  \nabla \bigl( ( u - k )_{+} \zeta \bigr)
 \right|^2 \,
d x \, d t \\
& \hspace*{5ex} \text{}
+ \left( \frac{1}{\varepsilon_1} + 1 \right)
\iint_{Q_{\rho}' \cap \{ u(x,t) > k \}}
 \sum_{i=1}^{n} | f_i |^2 \,
d x \, d t \\
& \quad + n \| \nabla \zeta \|_{L^{\infty} ( Q_{\rho} )}^2
\iint_{Q_{\rho}'}
 ( u - k )_{+}^2 \,
d x \, d t ;
\end{align*}
that is,
\begin{align*}
& \frac{1}{2} \int_{B_{\rho} ( x_{0} )}
 ( u - k )_{+}^2 \zeta^2 \,
d x \biggr|_{t = t'}
+ ( \lambda - \varepsilon_1 )
\iint_{Q_{\rho}'}
 \bigl| \nabla \bigl( ( u - k )_{+} \zeta \bigr) \bigr|^2 \,
d x \, d t \\
& \leq ( \Lambda + n ) \Big(
 \|
  \frac{\partial \zeta}{\partial t}
 \|_{L^{\infty} ( Q_{\rho} )}
 + \| \nabla \zeta \|_{L^{\infty} ( Q_{\rho} )}^2
\Big)
\iint_{Q_{\rho}'}
 ( u - k )_{+}^2 \,
d x \, d t \\
& \quad + \big( \frac{1}{\varepsilon_1} + 1 \big)
\iint_{Q_{\rho}' \cap \{ u(x,t) > k \}}
 \sum_{i=1}^{n} | f_i |^2 \,
d x \, d t
+ \iint_{Q_{\rho}'}
 f ( u - k )_{+} \zeta^2 \,
d x \, d t.
\end{align*}
Taking the supremum of the inequality over
$( t_{0} - \rho^2 , t_{0} ]$
with respect to $t'$, we have
\begin{equation}
\begin{aligned}
& \max \big\{
 \frac{1}{2}
 \| ( u - k )_{+} \zeta \|_{
  L^{\infty} ( t_{0} - \rho^2 , t_{0} ; L^2 ( B_{\rho} ( x_{0} ) ) )
 }^2 , \
 ( \lambda - \varepsilon_1 )
 \bigl\|
  \nabla \bigl( ( u - k )_{+} \zeta \bigr)
 \bigr\|_{L^2 ( Q_{\rho} )}^2
\big\}  \\
& \leq
( \Lambda + n ) \big(
 \|
  \frac{\partial \zeta}{\partial t}
 \|_{L^{\infty} ( Q_{\rho} )}
 + \| \nabla \zeta \|_{L^{\infty} ( Q_{\rho} )}^2
\big)
\iint_{Q_{\rho}}
 ( u - k )_{+}^2 \,
d x \, d t  \\
& \quad + \big( \frac{1}{\varepsilon_1} + 1 \big)
\iint_{Q_{\rho} \cap \{ u(x,t) > k \}}
 \sum_{i=1}^{n} | f_i |^2 \,
d x \, d t
+ \iint_{Q_{\rho}}
 f ( u - k )_{+} \zeta^2 \,
d x \, d t. 
\end{aligned} \label{eq:proofDG1}
\end{equation}
Now we estimate
the last two terms in the right-hand side of \eqref{eq:proofDG1}.
First we obtain
\begin{equation}\label{eq:proofDG2}
\iint_{Q_{\rho} \cap \{ u(x,t) > k \}}
 | f_i |^2 \,
d x \, d t
\leq \bigl|
 Q_{\rho} \cap \{ u(x,t) > k \}
\bigr|^{1 - 2 / p}
\| f_i \|_{L^{p} ( Q_{\rho} )}^2
\end{equation}
by H{\"o}lder's inequality.
Now we estimate
$\iint_{Q_{\rho}} f ( u - k )_{+} \zeta^2 \, d x \, d t$.
We first recall
\begin{align*}
& \| ( u - k )_{+} \zeta \|_{L^{2(n+2)/n} ( Q_{\rho} )} \\
& \leq C_1 \left(
 \| ( u - k )_{+} \zeta \|_{
  L^{\infty} ( t_{0} - \rho^2 , t_{0} ; L^2 ( B_{\rho} ( x_{0} ) ) )
 }
 + \bigl\|
  \nabla \bigl( ( u - k )_{+} \zeta \bigr)
 \bigr\|_{L^2 ( Q_{\rho} )}
\right)
\end{align*}
by Lemma~\ref{lemma:SobolevV2},
where $C_1 > 0$ depends only on $n$.
Then, by this inequality,
H{\"o}lder's inequality and Young's inequality,
we have
\begin{equation}
\begin{aligned}
& \iint_{Q_{\rho}} f ( u - k )_{+} \zeta^2 \, d x \, d t
 \\
& \leq
\| f \zeta \|_{
 L^{2(n+2)/(n+4)} ( Q_{\rho} \cap \{ u(x,t) > k \} )
}
\| ( u - k )_{+} \zeta \|_{L^{2(n+2)/n} ( Q_{\rho} )}
  \\
& \leq \varepsilon_2
\| ( u - k )_{+} \zeta \|_{L^{2(n+2)/n} ( Q_{\rho} )}^2
+ \frac{1}{\varepsilon_2}
\| f \zeta \|_{
 L^{2(n+2)/(n+4)} ( Q_{\rho} \cap \{ u(x,t) > k \} )
}^2   \\
& \leq 2 \varepsilon_2 C_1^2
\max \left\{
 \| ( u - k )_{+} \zeta \|_{
  L^{\infty} ( t_{0} - \rho^2 , t_{0} ; L^2 ( B_{\rho} ( x_{0} ) ) )
 }^2 , \
 \bigl\|
  \nabla \bigl( ( u - k )_{+} \zeta \bigr)
 \bigr\|_{L^2 ( Q_{\rho} )}^2
\right\}  \\
& \quad + \frac{1}{\varepsilon_2}
\| f \|_{
 L^{\frac{p(n+2)}{n+2+p}} ( Q_{\rho} )
}^2
\bigl|
 Q_{\rho} \cap \{ u(x,t) > k \}
\bigr|^{1 - 2 / p}
\end{aligned} \label{eq:proofDG3}
\end{equation}
because
$2(n+2)/(n+4) < p(n+2)/(n+2+p)$.
By \eqref{eq:proofDG1}, \eqref{eq:proofDG2} and \eqref{eq:proofDG3},
we obtain the estimate \eqref{eq:DG}.
\end{proof}

By the same argument, we obtain the following lemma for the function
$v_{-} (x) := \max \{ - v(x), 0 \}$.

\begin{lemma}
Under the same assumption as in Lemma~\ref{lemma:DG},
a solution $u$ to the parabolic equation
 \eqref{eq:parabolic} satisfies
\begin{equation}
\begin{aligned}
& \| ( u - k )_{-} \zeta \|_{
 L^{\infty} ( t_{0} - \rho^2 , t_{0} ; L^2 ( B_{\rho} ( x_{0} ) ) )
}^2
+ \bigl\|
 \nabla \bigl( ( u - k )_{-} \zeta \bigr)
\bigr\|_{L^2 ( Q_{\rho} )}^2  \\
& \leq C_2 \Bigl[
 \Big(
  \|
   \frac{\partial \zeta}{\partial t}
  \|_{L^{\infty} ( Q_{\rho} )}
  + \| \nabla \zeta \|_{L^{\infty} ( Q_{\rho} )}^2
 \Big) \| ( u - k )_{-} \|_{L^2 ( Q_{\rho} )}^2  \\
 & \quad  + F_{0, \rho}^2
 \bigl|
  Q_{\rho}  \cap \{ u(x,t) < k \}
 \bigr|^{1 - 2 / p}
\Bigr]
\end{aligned} \label{eq:DGprime}
\end{equation}
for any $k \in \mathbb{R}$,
where we define $F_{0, \rho}$ as  \eqref{eq:F0rho},
and $C_2 > 0$ depends only on $n, \Lambda$ and $\lambda$.
\end{lemma}

The estimate \eqref{eq:LinftyL2} easily follows from Lemmas~\ref{lemma:DG}
and \ref{lemma:DG}$'$.
Our next task is to prove the estimate \eqref{eq:Linfty}.
We start by giving a technical lemma which will be used later.

\begin{lemma}\label{lemma:ym}
Let $\widetilde{C} > 0$, $b > 1$ and $\varepsilon > 0$.
If a sequence $\{ y_{m} \}_{m=0}^{\infty}$ satisfies
\begin{equation}\label{eq:ymcond}
y_{0} \leq \theta_{0}
:= \widetilde{C}^{- 1 / \varepsilon} b^{- 1 / \varepsilon^2}
\quad\text{and }
0 \leq y_{m+1} \leq \widetilde{C} b^{m} y_{m}^{1 + \varepsilon} ,
\end{equation}
then
\[
\lim_{m \to \infty} y_{m} = 0
\]
holds.
\end{lemma}

\begin{proof}
We show
\begin{equation}\label{eq:ymr}
y_{m} \leq \frac{\theta_{0}}{r^{m}} , \quad
m = 0, 1, 2, \ldots
\end{equation}
by inductive method,
where we will determine $r > 1$ later.
By assumption, \eqref{eq:ymr} with $m=0$ holds.
Hence we now assume \eqref{eq:ymr} holds,
and show \eqref{eq:ymr} for $m+1$.
By the assumption \eqref{eq:ymcond}
and the induction hypothesis, we have
\[
y_{m+1} \leq \widetilde{C} b^{m} y_{m}^{1 + \varepsilon}
\leq \widetilde{C} b^{m}
\big( \frac{\theta_{0}}{r^{m}} \big)^{1 + \varepsilon}
= \frac{\theta_{0}}{r^{m+1}} \widetilde{C}
b^{m} \frac{\theta_{0}^{\varepsilon}}{r^{m \varepsilon - 1}} .
\]
Now we take $r = b^{1 / \varepsilon}$. Then we have
\[
y_{m+1} \leq
\frac{\theta_{0}}{r^{m+1}} \widetilde{C}
b^{m} \frac{\theta_{0}^{\varepsilon}}{r^{m \varepsilon - 1}}
= \frac{\theta_{0}}{r^{m+1}} \widetilde{C}
r \theta_{0}^{\varepsilon}
= \frac{\theta_{0}}{r^{m+1}} ,
\]
which is \eqref{eq:ymr} for $m+1$.
\end{proof}

Now we are now ready to show the estimate \eqref{eq:Linfty}.
The estimate easily follows if we have the following lemma.

\begin{lemma}\label{lemma:Linfty}
Let $p > n + 2$.
Then a solution $u$ to {\rm \eqref{eq:parabolic}}
satisfies the estimate
\[
\| u \|_{L^{\infty} ( Q_{\rho} )}
\leq C_{\rho}
\left(
 \| u \|_{L^2 ( Q_{2 \rho} )}
 + F_{0, 2 \rho}
\right) ,
\]
where we define $F_{0, 2 \rho}$ by {\rm \eqref{eq:F0rho}},
and $C_{\rho} > 0$ depends only on $n, \lambda , \Lambda , p$
and $\rho$.
\end{lemma}

\begin{proof}
First of all a letter $C$ denotes a general constant
depending only on $n, \Lambda , \lambda$ and $p$.
Now, let
$\rho_{m} := ( 1 + 2^{-m} ) \rho$ and
$k_{m} = k ( 2 - 2^{-m} )$
for $m = 0, 1, 2, \ldots$,
where we will determine $k > 0$ later.
For $m = 0, 1, 2, \ldots$,
we take cut-off functions
\begin{math}
\zeta_{m} \in C^{\infty} ( Q_{\rho_{m}} )
\end{math}
which satisfy
\begin{gather*}
 0 \leq \zeta_{m} \leq 1 \quad\text{in } Q_{\rho_{m}} ,
 \\
 \zeta_{m} =  \begin{cases}
  1 &\text{in } Q_{\rho_{m+1}} , \\
  0 &\text{in }  Q_{\rho_{m}} \setminus Q_{( \rho_{m} + \rho_{m+1} ) / 2} ,
 \end{cases}  \\
 \|
 \frac{\partial \zeta_{m}}{\partial t}
\|_{L^{\infty} ( Q_{\rho_{m}} )}
+ \| \nabla \zeta_{m} \|_{L^{\infty} ( Q_{\rho_{m}} )}^2
\leq \frac{C}{( \rho_{m} - \rho_{m+1} )^2} .
\end{gather*}
We remark that
$\zeta_{m} = 0$
near
\begin{math}
B_{\rho_{m}} ( x_{0} ) \times \{ t_{0}- \rho_{m}^2 \}
\cup \partial B_{\rho_{m}} ( x_{0} )
\times ( t_{0} - \rho_{m}^2 , t_{0} )
\end{math}
in particular.
By Lemmas~\ref{lemma:SobolevV2} and \ref{lemma:DG}, we have
\begin{equation}
\begin{aligned}
& \|
 ( u - k_{m+1} )_{+} \zeta_{m}
\|_{L^{2(n+2)/n} ( Q_{\rho_{m}} )}^2  \\
& \leq C \Bigl(
 \| ( u - k_{m+1} )_{+} \zeta_{m} \|_{
  L^{\infty} (
   t_{0} - \rho_{m}^2 , t_{0} ; L^2 ( B_{\rho_{m}} ( x_{0} ) )
  )
 }^2  \\
 & \hspace*{10ex} \text{}
 + \bigl\|
  \nabla \bigl( ( u - k_{m+1} )_{+} \zeta_{m} \bigr)
 \bigr\|_{L^2 ( Q_{\rho_{m}} )}^2
\Bigr)   \\
& \leq C \Bigl[
 \Big(
  \|
   \frac{\partial \zeta_{m}}{\partial t}
  \|_{L^{\infty} ( Q_{\rho_{m}} )}
  + \| \nabla \zeta_{m} \|_{L^{\infty} ( Q_{\rho_{m}} )}^2
 \Big)
 \| ( u - k_{m+1} )_{+} \|_{L^2 ( Q_{\rho_{m}} )}^2
  \\
 & \quad  + F_{0, \rho_{m}}^2
 \bigl|
  Q_{\rho_{m}}  \cap \{ u(x,t) > k_{m+1} \}
 \bigr|^{1 - 2 / p}
\Bigr]   \\
& \leq C \Bigl[
 \frac{2^{2 m}}{\rho^2}
 \| ( u - k_{m+1} )_{+} \|_{L^2 ( Q_{\rho_{m}} )}^2
 + F_{0, 2 \rho}^2
 \bigl| A_{m} ( k_{m+1} ) \bigr|^{1 - 2 / p}
\Bigr] ,
\end{aligned} \label{eq:a001}
\end{equation}
where $A_{m} (l) := Q_{\rho_{m}} \cap \{ u(x,t) > l \}$
for $l \in \mathbb{R}$.
Now we take $k > 0$ as
\begin{equation}\label{eq:k1}
k \geq \rho^{1 - (n+2) / p}  F_{0, 2 \rho} .
\end{equation}
Then we have
\begin{align*}
& \|
 ( u - k_{m+1} )_{+} \zeta_{m}
\|_{L^{2(n+2)/n} ( Q_{\rho_{m}} )}^2 \\
& \leq C \Bigl[
 \frac{2^{2 m}}{\rho^2}
 \| ( u - k_{m+1} )_{+} \|_{L^2 ( Q_{\rho_{m}} )}^2
 + \frac{k^2}{\rho^{2 ( 1 - (n+2) / p )}}
 \bigl| A_{m} ( k_{m+1} ) \bigr|^{1 - 2 / p}
\Bigr]
\end{align*}
by the estimate \eqref{eq:a001}.
By defining
\begin{math}
\varphi_{m} := \|
 ( u - k_{m} )_{+}
\|_{L^2 ( Q_{\rho_{m}} )}^2,
\end{math}
we have
\begin{equation}
\begin{aligned}
&\varphi_{m+1}\\
& = \|
 ( u - k_{m+1} )_{+} \zeta_{m}
\|_{L^2 ( Q_{\rho_{m+1}} )}^2
\leq \|
 ( u - k_{m+1} )_{+} \zeta_{m}
\|_{L^2 ( Q_{\rho_{m}} )}^2  \\
& \leq | A_{m} ( k_{m+1} ) |^{2/(n+2)}
\|
 ( u - k_{m+1} )_{+} \zeta_{m}
\|_{L^{2(n+2)/n} ( Q_{\rho_{m}} )}^2
 \\
& \leq C | A_{m} ( k_{m+1} ) |^{2/(n+2)}  \\
& \quad\times
\Bigl[
 \frac{2^{2 m}}{\rho^2}
 \| ( u - k_{m+1} )_{+} \|_{L^2 ( Q_{\rho_{m}} )}^2
 + \frac{k^2}{\rho^{2 ( 1 - (n+2) / p )}}
 \bigl| A_{m} ( k_{m+1} ) \bigr|^{1 - 2 / p}
\Bigr]   \\
& \leq C | A_{m} ( k_{m+1} ) |^{2/(n+2)}
\Bigl[
 \frac{2^{2 m}}{\rho^2} \varphi_{m}
 + \frac{k^2}{\rho^{2 ( 1 - (n+2) / p )}}
 \bigl| A_{m} ( k_{m+1} ) \bigr|^{1 - 2 / p}
\Bigr] ,
\end{aligned} \label{eq:a002}
\end{equation}
where we used H{\"o}lder's inequality
and the estimate
\[
\| ( u - k_{m+1} )_{+} \|_{L^2 ( Q_{\rho_{m}} )}^2
\leq \| ( u - k_{m} )_{+} \|_{L^2 ( Q_{\rho_{m}} )}^2
= \varphi_{m} .
\]
On the other hand, we have
\begin{align*}
\varphi_{m}
& = \| ( u - k_{m} )_{+} \|_{L^2 ( Q_{\rho_{m}} )}^2
\geq \iint_{A_{m} ( k_{m+1} )}
 ( u - k_{m} )_{+}^2 \,
d x \, d t \\
& \geq \iint_{A_{m} ( k_{m+1} )}
 ( k_{m+1} - k_{m} )_{+}^2 \,
d x \, d t
= \frac{k^2}{2^{2m+2}}
| A_{m} ( k_{m+1} ) | ;
\end{align*}
that is,
\begin{equation}\label{eq:a003}
| A_{m} ( k_{m+1} ) |
\leq \frac{2^{2m+2}}{k^2} \varphi_{m} .
\end{equation}
By \eqref{eq:a002} and \eqref{eq:a003}, we have
\begin{equation}
\begin{aligned}
\varphi_{m+1}
& \leq C 2^{2 m \left( 1 + \frac{2}{n+2} \right)}  \\
& \quad \times \Big[
 \rho^{-2} k^{- \frac{4}{n+2}} \varphi_{m}^{1 + \frac{2}{n+2}}
 + \rho^{- 2 \left( 1 - \frac{n+2}{p} \right)}
 k^{- \left( \frac{4}{n+2} - \frac{4}{p} \right)}
 \varphi_{m}^{1 + \frac{2}{n+2} - \frac{2}{p}}
\Big] .
\end{aligned} \label{eq:a004}
\end{equation}
We now take $k$ as
\begin{equation}\label{eq:k2}
k \geq \Big(
 \frac{1}{| Q_{2 \rho} |}
 \iint_{Q_{2 \rho}} u^2 \, d x \, d t
\Big)^{1/2} .
\end{equation}
Then we have
\[
\varphi_{m}
\leq \iint_{Q_{\rho_{m}}} u^2 \, d x \, d t
\leq \iint_{Q_{2 \rho}} u^2 \, d x \, d t
\leq | Q_{2 \rho} | k^2 ;
\]
that is,
\[
\varphi_{m}^{2/p} \leq
| Q_{2 \rho} |^{2/p} k^{4/p} .
\]
By this inequality and \eqref{eq:a004}, we have
\begin{equation}
\begin{aligned}
\varphi_{m+1}
& \leq C 2^{2 m \left( 1 + \frac{2}{n+2} \right)}
\varphi_{m}^{1 + \frac{2}{n+2} - \frac{2}{p}}
\left[
 \rho^{-2} k^{- \frac{4}{n+2}}
 \varphi_{m}^{2/p}
 + \rho^{- 2 \left( 1 - \frac{n+2}{p} \right)}
 k^{- \left( \frac{4}{n+2} - \frac{4}{p} \right)}
\right]  \\
& \leq C 2^{2 m \left( 1 + \frac{2}{n+2} \right)}
\varphi_{m}^{1 + \frac{2}{n+2} - \frac{2}{p}}  \\
& \hspace*{3ex} \text{} \times
\left[
 \rho^{-2} k^{- \frac{4}{n+2}}
 | Q_{2 \rho} |^{2/p} k^{4/p}
 + \rho^{- 2 \left( 1 - \frac{n+2}{p} \right)}
 k^{- \left( \frac{4}{n+2} - \frac{4}{p} \right)}
\right]   \\
& \leq C 2^{2 m \left( 1 + \frac{2}{n+2} \right)}
\rho^{- 2 \left( 1 - \frac{n+2}{p} \right)}
k^{- \frac{4}{n+2} \left( 1 - \frac{n+2}{p} \right)}
\varphi_{m}^{1 + \frac{2}{n+2} - \frac{2}{p}} .
\end{aligned} \label{eq:a005}
\end{equation}
Now we denote
$y_{m} := k^{-2} | Q_{2 \rho} |^{-1} \varphi_{m}$.
Then by \eqref{eq:a005}, we have
\begin{equation}\label{eq:a006}
y_{m+1}
\leq C 2^{2 m \left( 1 + \frac{2}{n+2} \right)}
y_{m}^{1 + \left( \frac{2}{n+2} - \frac{2}{p} \right)} ,
\end{equation}
which is the second condition of \eqref{eq:ymcond}
with
\begin{equation}\label{eq:a008}
\widetilde{C} = C, \quad
b = 2^{2 ( 1 + \frac{2}{n+2} )},
\quad
\varepsilon = \frac{2}{n+2} - \frac{2}{p} .
\end{equation}
Then
$\lim_{m \to \infty} y_{m} = 0$
if
\begin{equation}\label{eq:a007}
y_{0} \leq C^{- 1 / \varepsilon} b^{- 1 / \varepsilon^2}
=: \theta_{0}
\end{equation}
by Lemma~\ref{lemma:ym},
where $b$ and $\varepsilon$ are defined by \eqref{eq:a008}
and $C$ is the constant $C$ in \eqref{eq:a006}.
We remark that the condition \eqref{eq:a007} is equivalent to
\begin{equation}\label{eq:a009}
\| ( u - k )_{+} \|_{L^2 ( Q_{2 \rho} )}^2
\leq \theta_{0} k^2 | Q_{2\rho} | .
\end{equation}
Now we take $k$ as
\begin{equation}\label{eq:k3}
k^2 \geq \frac{1}{\theta_{0} | Q_{2\rho} |}
\| u \|_{L^2 ( Q_{2 \rho} )}^2 .
\end{equation}
Then the condition \eqref{eq:a009}; i.e.,
the condition \eqref{eq:a007} is satisfied.

Summing up, if we take $k$ such that
the conditions \eqref{eq:k1}, \eqref{eq:k2} and \eqref{eq:k3} are satisfied,
then we have
\begin{math}
\lim_{m \to \infty} y_{m} = 0
\end{math}.
On the other hand, since
\begin{align*}
y_{m}
& = \frac{1}{k^2 | Q_{2 \rho} |} \varphi_{m}
= \frac{1}{k^2 | Q_{2 \rho} |}
\|
 ( u - k_{m} )_{+}
\|_{L^2 ( Q_{\rho_{m}} )}^2 \\
& \to \frac{1}{k^2 | Q_{2 \rho} |}
\|
 ( u - 2 k )_{+}
\|_{L^2 ( Q_{\rho} )}^2
\quad\text{as } m \to \infty .
\end{align*}
Then we have
$\| ( u - 2 k )_{+} \|_{L^2 ( Q_{\rho} )}^2 = 0$;
that is,
\begin{equation}\label{eq:Linftyresult}
u \leq 2 k \quad\text{a.e. in } Q_{\rho} .
\end{equation}
Now we take $k$ as
\[
k = \frac{1}{\sqrt{\theta_{0} | Q_{2 \rho} |}}
\| u \|_{L^2 ( Q_{2 \rho})}
+ \rho^{1 - (n+2) / p} F_{0, 2 \rho} ,
\]
which satisfies
the conditions \eqref{eq:k1}, \eqref{eq:k2} and \eqref{eq:k3}.
Hence we have \eqref{eq:Linftyresult}, which is
\[
\sup_{Q_{\rho}} u \leq C_{\rho}
\left(
 \| u \|_{L^2 ( Q_{2 \rho} )} + F_{0, 2 \rho}
\right) .
\]
Replacing Lemma~\ref{lemma:DG}
by Lemma~\ref{lemma:DG}$'$
and doing the same argument, we can obtain
\[
- u \leq C_{\rho}
\left(
 \| u \|_{L^2 ( Q_{2 \rho} )} + F_{0, 2 \rho}
\right) \quad\text{in } Q_{\rho}
\]
and thus the proof has been completed.
\end{proof}

\subsection*{Acknowledgments}
We would like to acknowledge the referee's careful reading of the paper
and comments which improved the paper.
The first author is supported
by Postdoctoral Fellowship for Foreign Researchers
from the Japan Society for the Promotion of Science.
The third and fourth authors are supported
by Research Fellowships of the Japan Society
for the Promotion of Science for Young Scientists and
Grant-in-Aid for Scientific Research (B)
(No.\ 22340023) of Japan Society for Promotion of Science, respectively.

\begin{thebibliography}{99}

\bibitem{Aronson} D.\ G.\ Aronson;
\emph{Non-negarive solutions of linear parabolic equations},
Ann.\ Scuola Norm.\ Sup.\ Pisa (3) {\bf 22} (1968), 607--694.
\bibitem{Babuska}

I.\ Babu{\v{s}}ka, B.\ Andersson, P.\ Smith and K.\ Levin;
\emph{Damage analysis of fiber composites,
Part {\rm I}: Statisical analysis on fiber scale},
Comput.\ Methods Appl.\ Mech.\ Engrg.\ {\bf 172} (1999),
no.\ 1-4, 27--77.

\bibitem{DiCristoVessella} M.\ Di\ Cristo and S.\ Vessella;
\emph{Stable determination of the discontinuous conductivity
coefficient of a parabolic equation},
SIAM J.\ Math.\ Anal.\ {\bf 42} (2010), no.\ 1, 183--217.

\bibitem{FS} E.\ B.\ Fabes and D.\ W.\ Stroock;
\emph{A new proof of Moser's parabolic Harnack inequality
using the old ideas of Nash},
Arch.\ Rational Mech.\ Anal.\ {\bf 96} (1986), no.\ 4, 327--338.

\bibitem{Friedman} A.\ Friedman;
\emph{Partial differential equations}, 2008,
Dover Publications.

\bibitem{Ikehata} M.\ Ikehata;
\emph{The enclosure method and its applications},
Analytic extension formulas and their applications
(Fukuoka, 1999/Kyoto, 2000) 87--103,
Int.\ Soc.\ Anal.\ Appl.\ Comput.\ 9,
Kluwer Acad.\ Publ., Dordrecht, 2001.

\bibitem{IKN} V.\ Isakov, K.\ Kim and G.\ Nakamura;
\emph{Reconstruction of an unknown inclusion by thermography}, to appear in
Ann.\ Scuola Norm.\ Sup.\ Pisa.

\bibitem{LRU} O.\ A.\ Lady{\v{z}}enskaja, V.\ Ya.\ Rivkind
and N.\ N.\ Ural'ceva;
\emph{The classical solvability of diffraction problems},
Trudy Mat.\ Inst.\ Steklov {\bf 92} (1966), 116--146.

\bibitem{LSU} O.\ A.\ Lady{\v{z}}enskaja, V.\ A.\ Solonnikov,
N.\ N.\ Ural'ceva;
\emph{Linear and quasi-linear equations of parabolic type,
Translated from the Russian by S.\ Smith}, 1967,
American Mathematical Society.

\bibitem{HGLi} H.\ Li;
\emph{Several problems for partial differential equationes
from astrophysics and physics: existence of compressible
rotating Newtonian stars and linear elasticity of
composite material},
Thesis, Beijing Normal University, 2009.

\bibitem{LiLi} H.\ Li and Y.\ Li;
\emph{Gradient Estimates for Parabolic Systems from Composite Material},
arXiv:1105.1437v1.

\bibitem{LiNirenberg} Y.\ Li and L.\ Nirenberg;
\emph{Estimates for elliptic systems from composite material},
Comm.\ Pure Appl.\ Math.\ {\bf 56} (2003), no.\ 7, 892--925.

\bibitem{LiVogelius} Y.\ Li and M.\ Vogelius;
\emph{Gradient estimates for solutions
to divergence form elliptic equations
with discontinuous coefficients},
Arch.\ Rational Mech.\ Anal. {\bf 153} (2000), 91--151.

\bibitem{Meyers} N.\ G.\ Meyers;
\emph{An $L^{p}$-estimate for the gradient of solutions
of second order elliptic divergence equations},
Ann.\ Scuola Norm.\ Sup.\ Pisa (3) {\bf 17} (1963), 189--206.

\bibitem{Nag-Nak} S.\ Nagayasu and G.\ Nakamura;
\emph{Enclosure method for identifying unknown inclusions
with possible contact}, in preparation.

\bibitem{NUW} S.\ Nagayasu, G.\ Uhlmann and J.-N.\ Wang;
\emph{Reconstruction of penetrable obstacles in acoustics},
SIAM J.\ Math.\ Anal.\ {\bf 43} (2011), 189--211.

\bibitem{PiccininiSpagnolo} L.\ C.\ Piccinini and S.\ Spagnolo;
\emph{On the H{\"o}lder continuity of solutions
of second order elliptic equations in two variables},
Ann.\ Scuola Norm.\ Sup.\ Pisa (3) {\bf 26} (1972), 391--402.

\bibitem{Yoshida} K.\ Yoshida;
\emph{Reconstruction of a penetrable obstacle
by complex spherical waves},
J.\ Math.\ Anal.\ Appl.\ {\bf 369} (2010), 645--657.

\end{thebibliography}

\end{document}